Previous tutorial: Anisotropic image segmentation using gradient structure tensors
| Original author | Karpushin Vladislav |
|---|---|
| Compatibility | OpenCV >= 3.0 |
Goals
In this tutorial you will learn
- How to remove periodic noise in Fourier domain
Theory
Comments
This tutorial is based on the book [98]. The images on this page are real world images.
Periodic noise produces spikes in the Fourier domain and can often be detected through visual analysis.
How to eliminate periodic noise in Fourier domain?
Periodic noise can be significantly reduced through frequency domain filtering. In this page, we use a notch filter with an appropriate radius to completely surround the noise spike in the Fourier domain. Notch filters filter out frequencies within a predefined neighborhood around a center frequency. The number of notch filters is arbitrary. The shape of the notch area can also be arbitrary (such as rectangular or circular). In this page we use three circular notch suppression filters. The power spectral density of the image is used for visual detection of noise spikes.
Source code
You can find the source code in the OpenCV Source code repository at samples/cpp/tutorial_code/ImgProc/periodic_noise_removing_filter/periodic_noise_removing_filter.cpp.
#include <iostream>
#include "opencv2/highgui.hpp"
#include <opencv2/imgcodecs.hpp>
#include <opencv2/imgproc.hpp>
using namespace cv;
using namespace std;
void fftshift(const Mat & amp; inputImg, Mat & amp; outputImg);
void filter2DFreq(const Mat & amp; inputImg, Mat & amp; outputImg, const Mat & amp; H);
void synthesizeFilterH(Mat & amp; inputOutput_H, Point center, int radius);
void calcPSD(const Mat & amp; inputImg, Mat & amp; outputImg, int flag = 0);
const String keys =
"{help h usage ? | | print this message }"
"{@image |period_input.jpg | input image name }"
;
int main(int argc, char* argv[])
{<!-- -->
CommandLineParser parser(argc, argv, keys);
string strInFileName = parser.get<String>("@image");
samples::addSamplesDataSearchSubDirectory("doc/tutorials/imgproc/periodic_noise_removing_filter/images");
Mat imgIn = imread(samples::findFile(strInFileName), IMREAD_GRAYSCALE);
if (imgIn.empty()) // Check whether the image has been loaded
{<!-- -->
cout << "ERROR : Image cannot be loaded...!!" << endl;
return -1;
}
imshow("Image corrupted", imgIn);
imgIn.convertTo(imgIn, CV_32F);
//Only need to process even-numbered images
Rect roi = Rect(0, 0, imgIn.cols & amp; -2, imgIn.rows & amp; -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);
//Calculate H (stop)
// Filter (start)
Mat imgOut;
fftshift(H, H);
filter2DFreq(imgIn, imgOut, H);
// filter (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);
imshow("Debluring", imgOut);
imwrite("filter.jpg", H);
waitKey(0);
return 0;
}
void fftshift(const Mat & amp; inputImg, Mat & amp; 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 & amp; inputImg, Mat & amp; outputImg, const Mat & amp; 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 & amp; 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);
}
// The function calculates PSD (power spectral density) through fft, there are two flags
// Flag = 0 means return PSD
// flag = 1 means return log(PSD)
void calcPSD(const Mat & amp; inputImg, Mat & amp; 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;
// Calculate PSD = sqrt(Re(DFT(I))^2 + Im(DFT(I))^2)^2
Mat imgPSD;
magnitude(planes[0], planes[1], imgPSD); //imgPSD = sqrt(power spectral 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;
}
}
Description
Periodic noise reduction with frequency domain filtering includes power spectral density calculation (for visual detection of noise spikes), notch suppression filter synthesis and frequency filtering:
// Only need to process even-numbered images Rect roi = Rect(0, 0, imgIn.cols & amp; -2, imgIn.rows & amp; -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); //Calculate H (stop) // Filter (start) Mat imgOut; fftshift(H, H); filter2DFreq(imgIn, imgOut, H); // filter (stop)
The function calcPSD() calculates the power spectral density of an image:
void calcPSD(const Mat & amp; inputImg, Mat & amp; 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;
// Calculate PSD = sqrt(Re(DFT(I))^2 + Im(DFT(I))^2)^2
Mat imgPSD;
magnitude(planes[0], planes[1], imgPSD); //imgPSD = sqrt(power spectral 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;
}
}
The function synthesizeFilterH() forms the transfer function of an ideal circular notch suppression filter based on the center frequency and radius:
void synthesizeFilterH(Mat & amp; 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);
}
The function filter2DFreq() performs frequency domain filtering on the image. The functions fftshift() and filter2DFreq() are copied from the Defocus Filter tutorial.
Filter results
The image below shows an image severely corrupted by periodic noise of various frequencies.
Image corrupted by periodic noise
In the power spectral density shown below, noise components can easily appear as bright points (spikes).
Display the power spectral density of periodic noise
The image below shows a notch rejection filter with an appropriate radius to completely surround the noise spike.
Notch Filter
The result of using the notch filter to process the image is shown below.
Filter results
The improvement is very noticeable. There is significantly less periodic noise visible in this image compared to the original image.
You can also find a quick video demonstration of this filtering method on YouTube.