OpenCV  4.5.1
Open Source Computer Vision
Periodic Noise Removing Filter

Prev Tutorial: Anisotropic image segmentation by a gradient structure tensor

Original author Karpushin Vladislav
Compatibility OpenCV >= 3.0

Goal

In this tutorial you will learn:

  • how to remove periodic noise in the Fourier domain

Theory

Note
The explanation is based on the book [gonzalez]. The image on this page is a real world image.

Periodic noise produces spikes in the Fourier domain that can often be detected by visual analysis.

How to remove periodic noise in the Fourier domain?

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.

Source code

You can find source code in the samples/cpp/tutorial_code/ImgProc/periodic_noise_removing_filter/periodic_noise_removing_filter.cpp of the OpenCV source code library.

#include <iostream>
using namespace cv;
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()
{
Mat imgIn = imread("input.jpg", IMREAD_GRAYSCALE);
if (imgIn.empty()) //check whether the image is loaded or not
{
cout << "ERROR : Image cannot be loaded..!!" << endl;
return -1;
}
imgIn.convertTo(imgIn, CV_32F);
// it needs to process even image only
Rect roi = Rect(0, 0, imgIn.cols & -2, imgIn.rows & -2);
imgIn = imgIn(roi);
// PSD calculation (start)
Mat imgPSD;
calcPSD(imgIn, imgPSD);
fftshift(imgPSD, imgPSD);
normalize(imgPSD, imgPSD, 0, 255, NORM_MINMAX);
// PSD calculation (stop)
//H calculation (start)
Mat H = Mat(roi.size(), CV_32F, Scalar(1));
const int r = 21;
synthesizeFilterH(H, Point(705, 458), r);
synthesizeFilterH(H, Point(850, 391), r);
synthesizeFilterH(H, Point(993, 325), r);
//H calculation (stop)
// filtering (start)
Mat imgOut;
fftshift(H, H);
filter2DFreq(imgIn, imgOut, H);
// filtering (stop)
imgOut.convertTo(imgOut, CV_8U);
normalize(imgOut, imgOut, 0, 255, NORM_MINMAX);
imwrite("result.jpg", imgOut);
imwrite("PSD.jpg", imgPSD);
fftshift(H, H);
normalize(H, H, 0, 255, NORM_MINMAX);
imwrite("filter.jpg", 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);
dft(complexI, complexI, DFT_SCALE);
Mat planesH[2] = { Mat_<float>(H.clone()), Mat::zeros(H.size(), CV_32F) };
Mat complexH;
merge(planesH, 2, complexH);
Mat complexIH;
mulSpectrums(complexI, complexH, complexIH, 0);
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;
c4 = Point(c3.x,c2.y);
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);
}
// Function calculates PSD(Power spectrum density) by fft with two flags
// flag = 0 means to return PSD
// flag = 1 means to return log(PSD)
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);
dft(complexI, complexI);
split(complexI, planes); // planes[0] = Re(DFT(I)), planes[1] = Im(DFT(I))
planes[0].at<float>(0) = 0;
planes[1].at<float>(0) = 0;
// compute the PSD = sqrt(Re(DFT(I))^2 + Im(DFT(I))^2)^2
Mat imgPSD;
magnitude(planes[0], planes[1], imgPSD); //imgPSD = sqrt(Power spectrum density)
pow(imgPSD, 2, imgPSD); //it needs ^2 in order to get PSD
outputImg = imgPSD;
// logPSD = log(1 + PSD)
if (flag)
{
Mat imglogPSD;
imglogPSD = imgPSD + Scalar::all(1);
log(imglogPSD, imglogPSD);
outputImg = imglogPSD;
}
}

Explanation

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:

// it needs to process even image only
Rect roi = Rect(0, 0, imgIn.cols & -2, imgIn.rows & -2);
imgIn = imgIn(roi);
// PSD calculation (start)
Mat imgPSD;
calcPSD(imgIn, imgPSD);
fftshift(imgPSD, imgPSD);
normalize(imgPSD, imgPSD, 0, 255, NORM_MINMAX);
// PSD calculation (stop)
//H calculation (start)
Mat H = Mat(roi.size(), CV_32F, Scalar(1));
const int r = 21;
synthesizeFilterH(H, Point(705, 458), r);
synthesizeFilterH(H, Point(850, 391), r);
synthesizeFilterH(H, Point(993, 325), r);
//H calculation (stop)
// filtering (start)
Mat imgOut;
fftshift(H, H);
filter2DFreq(imgIn, imgOut, H);
// filtering (stop)

A function calcPSD() calculates power spectrum density of an image:

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);
dft(complexI, complexI);
split(complexI, planes); // planes[0] = Re(DFT(I)), planes[1] = Im(DFT(I))
planes[0].at<float>(0) = 0;
planes[1].at<float>(0) = 0;
// compute the PSD = sqrt(Re(DFT(I))^2 + Im(DFT(I))^2)^2
Mat imgPSD;
magnitude(planes[0], planes[1], imgPSD); //imgPSD = sqrt(Power spectrum density)
pow(imgPSD, 2, imgPSD); //it needs ^2 in order to get PSD
outputImg = imgPSD;
// logPSD = log(1 + PSD)
if (flag)
{
Mat imglogPSD;
imglogPSD = imgPSD + Scalar::all(1);
log(imglogPSD, imglogPSD);
outputImg = imglogPSD;
}
}

A function synthesizeFilterH() forms a transfer function of an ideal circular shape notch reject filter according to a center frequency and a radius:

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;
c4 = Point(c3.x,c2.y);
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);
}

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.

Result

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.

You can also find a quick video demonstration of this filtering idea on YouTube.

cv::Mat::rows
int rows
the number of rows and columns or (-1, -1) when the matrix has more than 2 dimensions
Definition: mat.hpp:2092
cv::Point_< int >
cv::imread
CV_EXPORTS_W Mat imread(const String &filename, int flags=IMREAD_COLOR)
Loads an image from a file.
cv::Mat::clone
Mat clone() const CV_NODISCARD
Creates a full copy of the array and the underlying data.
cv::Scalar_< double >::all
static Scalar_< double > all(double v0)
returns a scalar with all elements set to v0
cv::Rect
Rect2i Rect
Definition: types.hpp:462
cv::NORM_MINMAX
@ NORM_MINMAX
flag
Definition: base.hpp:207
cv::idft
void idft(InputArray src, OutputArray dst, int flags=0, int nonzeroRows=0)
Calculates the inverse Discrete Fourier Transform of a 1D or 2D array.
cv::Mat::zeros
static MatExpr zeros(int rows, int cols, int type)
Returns a zero array of the specified size and type.
cv::split
void split(const Mat &src, Mat *mvbegin)
Divides a multi-channel array into several single-channel arrays.
cv::Mat::at
_Tp & at(int i0=0)
Returns a reference to the specified array element.
CV_8U
#define CV_8U
Definition: interface.h:73
cv::Point_::y
_Tp y
y coordinate of the point
Definition: types.hpp:187
cv::Point_::x
_Tp x
x coordinate of the point
Definition: types.hpp:186
cv::log
Quat< T > log(const Quat< T > &q, QuatAssumeType assumeUnit=QUAT_ASSUME_NOT_UNIT)
cv::Scalar_< double >
cv::dft
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.
cv::Rect_::size
Size_< _Tp > size() const
size (width, height) of the rectangle
cv::Mat::convertTo
void convertTo(OutputArray m, int rtype, double alpha=1, double beta=0) const
Converts an array to another data type with optional scaling.
CV_32F
#define CV_32F
Definition: interface.h:78
cv::Mat::empty
bool empty() const
Returns true if the array has no elements.
cv::Mat::cols
int cols
Definition: mat.hpp:2092
cv::Mat::size
MatSize size
Definition: mat.hpp:2114
imgcodecs.hpp
cv::IMREAD_GRAYSCALE
@ IMREAD_GRAYSCALE
If set, always convert image to the single channel grayscale image (codec internal conversion).
Definition: imgcodecs.hpp:71
cv::magnitude
void magnitude(InputArray x, InputArray y, OutputArray magnitude)
Calculates the magnitude of 2D vectors.
cv::mulSpectrums
void mulSpectrums(InputArray a, InputArray b, OutputArray c, int flags, bool conjB=false)
Performs the per-element multiplication of two Fourier spectrums.
cv::Rect_
Template class for 2D rectangles.
Definition: types.hpp:421
cv::merge
void merge(const Mat *mv, size_t count, OutputArray dst)
Creates one multi-channel array out of several single-channel ones.
cv::Scalar
Scalar_< double > Scalar
Definition: types.hpp:669
cv::DFT_SCALE
@ DFT_SCALE
Definition: base.hpp:231
cv::Point
Point2i Point
Definition: types.hpp:194
cv::Mat
n-dimensional dense array class
Definition: mat.hpp:798
cv::Mat::copyTo
void copyTo(OutputArray m) const
Copies the matrix to another one.
cv
"black box" representation of the file storage associated with a file on disk.
Definition: affine.hpp:52
imgproc.hpp
cv::pow
softfloat pow(const softfloat &a, const softfloat &b)
Raising to the power.
cv::Mat_< float >
cv::imwrite
CV_EXPORTS_W bool imwrite(const String &filename, InputArray img, const std::vector< int > &params=std::vector< int >())
Saves an image to a specified file.
cv::normalize
static Vec< _Tp, cn > normalize(const Vec< _Tp, cn > &v)
cv::datasets::circle
@ circle
Definition: gr_skig.hpp:62
cv::circle
void circle(InputOutputArray img, Point center, int radius, const Scalar &color, int thickness=1, int lineType=LINE_8, int shift=0)
Draws a circle.