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>
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()
{
{
cout << "ERROR : Image cannot be loaded..!!" << endl;
return -1;
}
imgIn = imgIn(roi);
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);
fftshift(H, H);
filter2DFreq(imgIn, imgOut, H);
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));
q3.copyTo(q0);
q1.copyTo(tmp);
q2.copyTo(q1);
}
void filter2DFreq(
const Mat& inputImg,
Mat& outputImg,
const Mat& H)
{
merge(planes, 2, complexI);
merge(planesH, 2, complexH);
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)
{
merge(planes, 2, complexI);
planes[0].
at<
float>(0) = 0;
planes[1].
at<
float>(0) = 0;
outputImg = imgPSD;
if (flag)
{
imglogPSD = imgPSD + Scalar::all(1);
log(imglogPSD, imglogPSD);
outputImg = imglogPSD;
}
}
n-dimensional dense array class
Definition: mat.hpp:811
CV_NODISCARD_STD Mat clone() const
Creates a full copy of the array and the underlying data.
MatSize size
Definition: mat.hpp:2138
void copyTo(OutputArray m) const
Copies the matrix to another one.
_Tp & at(int i0=0)
Returns a reference to the specified array element.
int cols
Definition: mat.hpp:2116
bool empty() const
Returns true if the array has no elements.
int rows
the number of rows and columns or (-1, -1) when the matrix has more than 2 dimensions
Definition: mat.hpp:2116
void convertTo(OutputArray m, int rtype, double alpha=1, double beta=0) const
Converts an array to another data type with optional scaling.
_Tp y
y coordinate of the point
Definition: types.hpp:202
_Tp x
x coordinate of the point
Definition: types.hpp:201
Template class for 2D rectangles.
Definition: types.hpp:444
Size_< _Tp > size() const
size (width, height) of the rectangle
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 normalize(InputArray src, InputOutputArray dst, double alpha=1, double beta=0, int norm_type=NORM_L2, int dtype=-1, InputArray mask=noArray())
Normalizes the norm or value range of an array.
void log(InputArray src, OutputArray dst)
Calculates the natural logarithm of every array element.
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.
void pow(InputArray src, double power, OutputArray dst)
Raises every array element to a power.
@ NORM_MINMAX
flag
Definition: base.hpp:207
@ DFT_SCALE
Definition: base.hpp:231
Rect2i Rect
Definition: types.hpp:489
Point2i Point
Definition: types.hpp:209
#define CV_8U
Definition: interface.h:73
#define CV_32F
Definition: interface.h:78
@ circle
Definition: gr_skig.hpp:62
@ 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.
"black box" representation of the file storage associated with a file on disk.
Definition: core.hpp:106
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.