Periodic noise produces spikes in the Fourier domain that can often be detected by visual analysis.
Periodic noise can be reduced significantly via frequency domain filtering. On this page we use a notch reject filter with an appropriate radius to completely enclose the noise spikes in the Fourier domain. The notch filter rejects frequencies in predefined neighborhoods around a center frequency. The number of notch filters is arbitrary. The shape of the notch areas can also be arbitrary (e.g. rectangular or circular). On this page we use three circular shape notch reject filters. Power spectrum densify of an image is used for the noise spike’s visual detection.
#include <iostream>
using namespace std;
void fftshift(const Mat& inputImg, Mat& outputImg);
void filter2DFreq(const Mat& inputImg, Mat& outputImg, const Mat& H);
void synthesizeFilterH(Mat& inputOutput_H,
Point center,
int radius);
void calcPSD(const Mat& inputImg, Mat& outputImg, int flag = 0);
int main()
{
if (imgIn.empty())
{
cout << "ERROR : Image cannot be loaded..!!" << endl;
return -1;
}
imgIn.convertTo(imgIn,
CV_32F);
Rect roi =
Rect(0, 0, imgIn.cols & -2, imgIn.rows & -2);
imgIn = imgIn(roi);
Mat imgPSD;
calcPSD(imgIn, imgPSD);
fftshift(imgPSD, imgPSD);
const int r = 21;
synthesizeFilterH(H,
Point(705, 458), r);
synthesizeFilterH(H,
Point(850, 391), r);
synthesizeFilterH(H,
Point(993, 325), r);
Mat imgOut;
fftshift(H, H);
filter2DFreq(imgIn, imgOut, H);
imgOut.convertTo(imgOut,
CV_8U);
fftshift(H, H);
return 0;
}
void fftshift(const Mat& inputImg, Mat& outputImg)
{
outputImg = inputImg.clone();
int cx = outputImg.cols / 2;
int cy = outputImg.rows / 2;
Mat q0(outputImg,
Rect(0, 0, cx, cy));
Mat q1(outputImg,
Rect(cx, 0, cx, cy));
Mat q2(outputImg,
Rect(0, cy, cx, cy));
Mat q3(outputImg,
Rect(cx, cy, cx, cy));
Mat tmp;
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);
}
void filter2DFreq(const Mat& inputImg, Mat& outputImg, const Mat& H)
{
Mat planes[2] = { Mat_<float>(inputImg.clone()),
Mat::zeros(inputImg.size(),
CV_32F) };
Mat complexI;
merge(planes, 2, complexI);
Mat complexH;
merge(planesH, 2, complexH);
Mat complexIH;
idft(complexIH, complexIH);
split(complexIH, planes);
outputImg = planes[0];
}
void synthesizeFilterH(Mat& inputOutput_H,
Point center,
int radius)
{
Point c2 = center, c3 = center, c4 = center;
c2.
y = inputOutput_H.rows - center.y;
c3.x = inputOutput_H.cols - center.x;
circle(inputOutput_H, center, radius, 0, -1, 8);
circle(inputOutput_H, c2, radius, 0, -1, 8);
circle(inputOutput_H, c3, radius, 0, -1, 8);
circle(inputOutput_H, c4, radius, 0, -1, 8);
}
void calcPSD(const Mat& inputImg, Mat& outputImg, int flag)
{
Mat planes[2] = { Mat_<float>(inputImg.clone()),
Mat::zeros(inputImg.size(),
CV_32F) };
Mat complexI;
merge(planes, 2, complexI);
planes[0].at<float>(0) = 0;
planes[1].at<float>(0) = 0;
Mat imgPSD;
outputImg = imgPSD;
if (flag)
{
Mat imglogPSD;
log(imglogPSD, imglogPSD);
outputImg = imglogPSD;
}
}
static MatExpr zeros(int rows, int cols, int type)
Returns a zero array of the specified size and type.
_Tp y
y coordinate of the point
Definition: types.hpp:187
static Scalar_< double > all(double v0)
returns a scalar with all elements set to v0
void split(const Mat &src, Mat *mvbegin)
Divides a multi-channel array into several single-channel arrays.
void mulSpectrums(InputArray a, InputArray b, OutputArray c, int flags, bool conjB=false)
Performs the per-element multiplication of two Fourier spectrums.
void magnitude(InputArray x, InputArray y, OutputArray magnitude)
Calculates the magnitude of 2D vectors.
void merge(const Mat *mv, size_t count, OutputArray dst)
Creates one multi-channel array out of several single-channel ones.
void idft(InputArray src, OutputArray dst, int flags=0, int nonzeroRows=0)
Calculates the inverse Discrete Fourier Transform of a 1D or 2D array.
void dft(InputArray src, OutputArray dst, int flags=0, int nonzeroRows=0)
Performs a forward or inverse Discrete Fourier transform of a 1D or 2D floating-point array.
@ NORM_MINMAX
flag
Definition: base.hpp:207
@ DFT_SCALE
Definition: base.hpp:231
Rect2i Rect
Definition: types.hpp:462
static Vec< _Tp, cn > normalize(const Vec< _Tp, cn > &v)
Point2i Point
Definition: types.hpp:194
Scalar_< double > Scalar
Definition: types.hpp:669
#define CV_8U
Definition: interface.h:73
#define CV_32F
Definition: interface.h:78
softfloat pow(const softfloat &a, const softfloat &b)
Raising to the power.
Quat< T > log(const Quat< T > &q, QuatAssumeType assumeUnit=QUAT_ASSUME_NOT_UNIT)
@ IMREAD_GRAYSCALE
If set, always convert image to the single channel grayscale image (codec internal conversion).
Definition: imgcodecs.hpp:71
CV_EXPORTS_W bool imwrite(const String &filename, InputArray img, const std::vector< int > ¶ms=std::vector< int >())
Saves an image to a specified file.
CV_EXPORTS_W Mat imread(const String &filename, int flags=IMREAD_COLOR)
Loads an image from a file.
void circle(InputOutputArray img, Point center, int radius, const Scalar &color, int thickness=1, int lineType=LINE_8, int shift=0)
Draws a circle.
"black box" representation of the file storage associated with a file on disk.
Definition: affine.hpp:52
Periodic noise reduction by frequency domain filtering consists of power spectrum density calculation (for the noise spikes visual detection), notch reject filter synthesis and frequency filtering:
A function synthesizeFilterH() forms a transfer function of an ideal circular shape notch reject filter according to a center frequency and a radius:
A function filter2DFreq() filters an image in the frequency domain. The functions fftshift() and filter2DFreq() are copied from the tutorial Out-of-focus Deblur Filter.
The figure below shows an image heavily corrupted by periodical noise of various frequencies. 
The noise components are easily seen as bright dots (spikes) in the Power spectrum density shown in the figure below. 
The figure below shows a notch reject filter with an appropriate radius to completely enclose the noise spikes. 
The result of processing the image with the notch reject filter is shown below. 
The improvement is quite evident. This image contains significantly less visible periodic noise than the original image.