Opencv note 21 frequency domain filtering

One 、 summary

        Fourier transform of image and its two important measures ： Amplitude spectrum and phase spectrum . Understand two important concepts ： Low frequency and high frequency . Low frequency refers to graph

The Fourier transform of “ Center position ” Nearby areas . Be careful , Unless otherwise specified , The Fourier transform of the image mentioned later is centralized . High frequency follows “ Center position ” The distance increases , That is, the peripheral region of the center of Fourier transform , there “ Center position ” It refers to the position of the maximum value of the amplitude spectrum corresponding to the Fourier transform .

The presentation of frequency domain filter in program or mathematical operation can be understood as a matrix , The width of the matrix 、 Height and width of Fourier transform of image 、 The height is the same , The following commonly used low-pass 、 qualcomm 、 Bandpass 、 Band stop Wait for the key steps of filtering , It is through certain criteria that the matrix is constructed .

Two 、 step

Steps are as follows ：

The following is a simple example to explain the whole steps of frequency domain filtering in detail .

First step ： Input image matrix I. Assuming that ：

The second step ： Multiply each pixel value of the image matrix by （-1）r+c Get the matrix I′,I′ =I.*（-1）r+c , Its

in r and c Represents the position index of the current pixel value in the matrix .

The third step ： Because the width and height of the image matrix are 7, In order to use the fast algorithm of Fourier transform , Yes I′ repair 0, Use command getOptimalDFTSize（7） Get a value no less than 7 And can be decomposed into 2p ×3q ×5r Minimum integer of , The calculation result is 8. So in the matrix I′ Fill in a line on the right and bottom of 0, Write it down as f：

  Step four ： The complex matrix is obtained by using the fast algorithm of Fourier transform F.

 

        OpenCV The complex matrix is stored in two channels , That is, the first channel stores the real part of the complex matrix , The second channel stores the imaginary part of the complex matrix .

Step five ： Build a frequency domain filter Filter. The frequency domain filter is essentially a fast Fourier transform matrix obtained in the fourth step F Have the same number of rows 、 Complex matrix of column number , In general, it is a real matrix , Let's say it's One is all 1 Matrix ：

Frequency domain filtering mentioned in this chapter , Such as low-pass filtering 、 High pass filtering 、 Custom filtering , The key step is The matrix is constructed by certain standards to complete the filtering of the image in the frequency domain .

Step six ： The fast Fourier transform matrix obtained in the fourth step F And the frequency domain filter obtained in step 5 Filter Multiply by the corresponding position of （ Dot product of matrix ）. Of course , If the filter is a real matrix , Then in the code implementation in , The real and imaginary parts of the Fourier transform are dot multiplied by the frequency domain filter respectively , namely Ffilter =F.*Filter, Because the filter constructed here is all 1 Matrix , therefore F filter =F.

Step seven ： For the dot multiplication matrix obtained in step 6 Ffilter Inverse Fourier transform , Get the complex matrix F′.

Step eight ： Take the complex matrix F′ The real part of .

Step nine ： Similar to step 2 , Multiply the matrix obtained in step 8 by （-1）r+c.

Step 10 ： Because in the step of fast Fourier transform 0 operation , So the real part moment obtained in step 9 The size of the array may be larger than the original , So we need to cut , Take the upper left corner of the real part matrix , The size is the same as the original drawing . The result of clipping , That is, the result of frequency domain filtering . In this example , Because the filter is all 1 Of

matrix , It is equivalent to not doing any processing to the original drawing , That is, the final filtering result is the same as the original image .

The frequency domain filtering algorithm is completed according to the above ten steps , Next, the construction of common filters is introduced in detail Method 、 Code implementation and its effect .

3、 ... and 、 Low pass filtering

  Fourier transform for image , Low frequency information indicates the area in the image where the gray value changes slowly ; High frequency information is just the opposite , Represents the part where the gray value changes rapidly , Like the edge . Low pass filtering , seeing the name of a thing one thinks of its function , Keep the low-frequency information of Fourier transform ; Or weaken the high-frequency information of Fourier transform ; High pass filtering is just the opposite , Retain the high-frequency information of Fourier transform , Remove or weaken the low-frequency information of Fourier transform .

Three commonly used low-pass filters ：

               H、W Respectively represent the high resolution of image fast Fourier transform 、 wide , The maximum of Fourier spectrum （ Center point ） The position of is maxR,maxC,radius Represents the truncation frequency ,D（r,c） Represents the distance to the center

1、 Ideal low-pass filtering ：ilpFilter=[ilpFilter(r,c)]H*w,

 2、 Butterworth low pass filter

The second is Butterworth low-pass filter , remember blpFilter=[blpFilter（r,c）]H×W

 3、 Gaussian low pass filter

remember glpFilter=[glpFilter（r,c）]H×W;

  effect ：

The closer the filter is to the center point, the closer the value is to 1, The farther away from the center, the smaller the value is 1, Multiplied by Fourier transform , It is equivalent to retaining low-frequency information , Weaken or remove high-frequency information

Four 、 Low pass filter code implementation

//  The fast Fourier transform 
void  fft2Image(InputArray  image, OutputArray F)
{
	// To the Mat  The type of 
	Mat  Image = image.getMat();
	int rows = Image.rows;
	int cols = Image.cols;
	//  Get the best two columns 
	int  r = getOptimalDFTSize(rows);
	int  c = getOptimalDFTSize(cols);

	//  Enter on the left and below 10
	Mat  f;
	copyMakeBorder(Image, f, 0, r - rows, 0, c - cols, BORDER_CONSTANT,Scalar::all(0));
	//  Two channels , The first channel stores the real part    The second channel is for storing imaginary parts 
	dft(f, F, DFT_COMPLEX_OUTPUT);

	// Inverse Fourier transform  
	//dft(f, F, DFT_INVERSE+DFT_REAL_OUTPUT+DFT_SCALE);// Return only the real part 
}


//  Amplitude spectrum 
void amplitudeSpectrum(InputArray  _srcfft, OutputArray  _dstSpectrum) 
{
	if (_srcfft.channels() != 2) return ;
	//  If there are two channels, start to separate channels   

	vector<Mat> FFT2Channels;
	split(_srcfft, FFT2Channels);

   //  Calculate the Fourier transform amplitude value     sqrt(pow(r,2)+pow(L,2))
	magnitude(FFT2Channels[0], FFT2Channels[1], _dstSpectrum);

}
//  Gray level display of Fourier spectrum 
Mat graySpectrum(Mat  s) 
{
	Mat  dst;
	log(s + 1, dst);
	//  normalization 
	normalize(dst,dst,0,1,NORM_MINMAX);
	//  To show   Gray scale   
	dst.convertTo(dst,CV_8UC1,255,0);
	return dst;
}
//   Calculate the phase angle of Fourier transform ==phase  operator 
Mat  phaseS(Mat sfft)
{
	//  Phase spectrum 
	Mat  phase;
	phase.create(sfft.size(), CV_64FC1);
	//  The separation channel 
	vector<Mat> fft2Channels;
	split(sfft,fft2Channels);
	//  Start calculating the phase spectrum 
	for (int r = 0; r < phase.rows;r++) {
		for (int c = 0; c < phase.cols; c++) {
			//  Real component    Imaginary part 
			double  real = fft2Channels[0].at<double>(r,c);
			double  img = fft2Channels[1].at<double>(r,c);
			// atan2  The return value of 【0,180】,【-180,0】

			phase.at<double>(r, c) = atan2(img,real);
		}
	}
	return phase;
}
// -----------------------------------------------------------

#if  1 // Low pass filtering 

enum typeFilter 
{
	 ILPFILTER,
	 BLPFILTER,
	 GLPFILTER,
};
Mat   createFilter(Size   size, Point  center,float radius, int  type, int n = 2)
{

	Mat  ilpFilter = Mat::zeros(size,CV_32FC1);  //  Define a size Size filter 
	int  rows = size.height;
	int  cols = size.width;

	if (radius <= 0)return ilpFilter;


	//  Build an ideal low-pass filter 
	if (type== ILPFILTER)
	{
		for (int  r=0;r< rows;r++)
		{
			for (int c = 0; c < cols;c++) 
			{
				//  Calculated distance 
				float  norm2 = pow(abs(float(r - center.y)), 2)+ pow(abs(float(r - center.y)), 2);
				if (norm2 < radius)
				{
					ilpFilter.at<float>(r, c) = 1;
				}
				else 
				{
					ilpFilter.at<float>(r, c) = 0;
				}

			}
		}
	}
	// 
	else if(type == BLPFILTER)
	{
		for (int r = 0; r < rows; r++)
		{
			for (int c = 0; c < cols; c++)
			{
	/*			float m = sqrt(pow(r - center.y, 2.0) + pow(c - center.x, 2.0));
				float  n = m / radius;
				float  o = pow(n,2*n);
				float  s = 1.0 + o;
				float  i = 1.0 / s;*/
				ilpFilter.at<float>(r, c) = float(1.0 / (1.0 + (pow(sqrt(pow(r - center.y, 2.0) + pow(c - center.x, 2.0)) / radius, 2 * n))));

			}
		}
	}
	//  Gauss filtering 
	else if (type == GLPFILTER)
	{
		for (int r = 0; r < rows; r++)
		{
			for (int c = 0; c < cols; c++)
			{
				ilpFilter.at<float>(r, c) = float(exp(-(pow(c-center.x,2.0)+pow(r-center.y,2.0))/(2*pow(radius,2.0))));
			}
		}
	}
	return  ilpFilter;
}

Mat    src;//  Input image 
Mat    F;//  Fast Fourier transform of image 
Point    maxLocation;// The coordinates of the maximum of the Fourier spectrum 
int   radius = 20; //  Truncation frequency 
const   int  max_Radius = 100;//  Maximum cutoff frequency 



Mat  ilpFilter;// low pass filter 
int  filterType=0;//
int  max_FilterType = 2;


Mat   F_ilpFilter;//  Low pass Fourier transform 
Mat   FlpSpetrum;//  Gray level of Fourier spectrum of low-pass Fourier transform 

Mat  result;//  The filtered image 

string  lpFilterSpectrum = " Low pass Fourier spectrum ";// 




void   callback_lpFilter(int ,void *);


int main()
{
	 src = imread("C:\\Users\\19473\\Desktop\\opencv_images\\611.png",CV_LOAD_IMAGE_GRAYSCALE);
	if (!src.data)
	{
		cout << " Failed to read the original graph " << endl;
		return -1;
	}
	imshow(" Original picture ",src);
	//  Data conversion      Convert data to double 
	Mat  f1; 
	src.convertTo(f1,CV_32FC1,1.0,0.0);
	// 2    Number of each   x pow(-1,c+r);
	for (int r = 0; r < f1.rows; r++)
	{
		for (int c = 0; c < f1.cols; c++)
		{
			//  Judge parity 
			if (r+c%2) 
			{
				f1.at<float>(r, c) *= -1;
			}
		}
	}

	//  Start Fourier transform 
	fft2Image(f1,F);// F  After Fourier transform 
	
	Mat as;
	// Get the Fourier spectrum 
	amplitudeSpectrum(F,as);

	// Gray level display of Fourier spectrum 
	Mat  s = graySpectrum(as);
	imshow(" Gray level display of original Fourier spectrum ",s);

	//  Find the coordinates of the maximum value of the Fourier spectrum 
	minMaxLoc(s,NULL,NULL,NULL,&maxLocation);

	//------------------------------------  Low pass filtering ----------------------------------------------
	namedWindow(lpFilterSpectrum,WINDOW_AUTOSIZE);
	createTrackbar(" Low pass filtering ：", lpFilterSpectrum,&filterType, max_FilterType, callback_lpFilter);
	createTrackbar(" radius ：", lpFilterSpectrum, &radius, max_Radius, callback_lpFilter);
	callback_lpFilter(0,0);
	waitKey(0);

	return 0;
}
//  Callback function 
void   callback_lpFilter(int, void*) 
{
	// -----------------  Build a low-pass filter ----------------------------------
	ilpFilter = createFilter(F.size(),maxLocation,radius, filterType,2);
	// --------------------  Low pass filter and Fourier transform of image start point multiplication =====================================
	F_ilpFilter.create(F.size(),F.type());


	for (int r = 0; r < F_ilpFilter.rows; r++)
	{
		for (int c = 0; c < F_ilpFilter.cols; c++)
		{
			//  Take out the values of the fast Fourier transform and the ideal low-pass filter of the current position respectively 
			Vec2f f_rc = F.at<Vec2f>(r,c);
			float  lpFilter_rc = ilpFilter.at<float>(r,c);
			//  The low-pass filter is multiplied by the corresponding position of the fast Fourier transform of the image 
			F_ilpFilter.at<Vec2f>(r, c) = f_rc * lpFilter_rc;
		}
	}

	//  Fourier spectrum of low-pass Fourier transform 
	amplitudeSpectrum(F_ilpFilter, FlpSpetrum);

	FlpSpetrum = graySpectrum(FlpSpetrum);
	imshow(lpFilterSpectrum, FlpSpetrum);



	//  Inverse Fourier transform    And as long as the real part 
	dft(F_ilpFilter,result,DFT_SCALE+DFT_INVERSE+DFT_REAL_OUTPUT);


	//  multiply   （-1） Of （r+c） Power 
	for (int r = 0; r < result.rows; r++)
	{
		for (int c = 0; c < result.cols; c++)
		{
			if ((c + r) % 2) result.at<float>(r, c) *= -1;
		}
	}

	//  An important step    Convert the result to  CV_8u 
	result.convertTo(result, CV_8UC1, 1.0, 0);
	result = result(Rect(0,0,src.cols,src.rows)).clone();
	imshow(" The image after low-pass filtering ",result);

}
#endif

  5、 ... and 、 High pass filtering

maxR,maxC Represents the position of the maximum value of the Fourier spectrum of the image ,D(r,c):

1、 Ideal high pass filtering

2. Butterworth high pass filter

3. Gaussian high pass filter

  6、 ... and 、 Bandpass Band stop filter

Band pass filtering refers to retaining only a certain range of frequency bands . Here are three commonly used bandpass filters . false set up BW Represents bandwidth ,D0 Represents the radial center of the bandwidth , Other symbols and low pass 、 High pass filtering is the same .

1. Ideal bandpass filter and ideal bandstop filter

 2. Butterworth bandpass filter and Butterworth bandstop filter

 3. Gaussian bandpass filter and Gaussian bandstop filter

  belt Resistance The filter is just the opposite of the above

7、 ... and 、 Custom filter

#if  1 // Custom filter 

Mat image;//  Input image 
Mat    ImageFFT;//  Fast Fourier transform of image 
Point    maxLocation;// The coordinates of the maximum of the Fourier spectrum 



Mat  ilpFilter;// low pass filter 


Mat   F_ImageSpetrum;//  Fourier spectrum of Fourier transform 

Mat  result;//  The filtered image 

string  windowName = " Gray level of Fourier amplitude spectrum ";// 

bool  drawing_box = false;//  Mouse events 

Point drawPoint;
Rect rectFilter;
bool  gotRectFilter = false;

void  mouseRectHandler(int event,int x,int  y,int,void* ) 
{
	switch (event)
	{
	case CV_EVENT_LBUTTONDOWN://  The mouse click 
		drawing_box = true;
		//  Record the starting point 
		drawPoint = Point(x, y);
		break;
	case CV_EVENT_MOUSEMOVE:
		if (drawing_box)
		{
			//  Move the mouse to the lower right corner 
			if (x> drawPoint.x&&y>= drawPoint.y)
			{
				rectFilter.x = drawPoint.x;
				rectFilter.y = drawPoint.y;
				rectFilter.width = x - drawPoint.x;
				rectFilter.height= y - drawPoint.y;
			}

			//  Move the mouse to the upper right corner 
			if (x >= drawPoint.x && y <= drawPoint.y)
			{
				rectFilter.x = drawPoint.x;
				rectFilter.y = y;
				rectFilter.width = x - drawPoint.x;
				rectFilter.height = drawPoint.y- y ;
			}

			//  Move the mouse to the upper left corner 
			if (x >= drawPoint.x && y <= drawPoint.y)
			{
				rectFilter.x = x;
				rectFilter.y = y;
				rectFilter.width = drawPoint.x - x;
				rectFilter.height = drawPoint.y - y;
			}

			//  Move the mouse to the lower left corner 
			if (x >= drawPoint.x && y <= drawPoint.y)
			{
				rectFilter.x = x;
				rectFilter.y = drawPoint.y;
				rectFilter.width =  drawPoint.x-x;
				rectFilter.height = y- drawPoint.y ;
			}
		}
		break;
	case CV_EVENT_LBUTTONUP:
		drawing_box = false;
		gotRectFilter = true;
		break;
	default:
		break;
	}
}


int main()
{
	image = imread("C:\\Users\\19473\\Desktop\\opencv_images\\612.png",CV_LOAD_IMAGE_GRAYSCALE);
	if (!image.data)
	{
		cout << " Failed to read the original graph " << endl;
		return -1;
	}
	imshow(" Original picture ", image);
	//  Data conversion      Convert data to double 
	Mat  f1; 
	image.convertTo(f1,CV_32FC1,1.0,0.0);
	// 2    Number of each   x pow(-1,c+r);
	for (int r = 0; r < f1.rows; r++)
	{
		for (int c = 0; c < f1.cols; c++)
		{
			//  Judge parity 
			if (r+c%2) 
			{
				f1.at<float>(r, c) *= -1;
			}
		}
	}

	//  Start Fourier transform 
	fft2Image(f1,ImageFFT);// F  After Fourier transform 
	
	Mat amplSpec;
	// Get the Fourier spectrum 
	amplitudeSpectrum(ImageFFT, amplSpec);

	// Gray level display of Fourier spectrum 
	Mat  spectrum = graySpectrum(amplSpec);
	

	//  Find the coordinates of the maximum value of the Fourier spectrum 
	minMaxLoc(amplSpec,NULL,NULL,NULL,&maxLocation);

	//------------------------------------  Customize    wave filtering ----------------------------------------------
	namedWindow(windowName,WINDOW_AUTOSIZE);

	setMouseCallback(windowName, mouseRectHandler,NULL);
	for (;;) 
	{

		spectrum(rectFilter).setTo(0);
		//  Customize   Filter and Fourier transform point multiplication 
		ImageFFT(rectFilter).setTo(Scalar::all(0));
		imshow(windowName, spectrum);
		//  Press esc  See exit editing 
		if (waitKey(10)==27)
		{
			break;
		}

	}

	Mat  result;

	//  Inverse Fourier transform    And as long as the real part 
	dft(ImageFFT, result, DFT_SCALE + DFT_INVERSE + DFT_REAL_OUTPUT);


	//  multiply   （-1） Of （r+c） Power 
	for (int r = 0; r < result.rows; r++)
	{
		for (int c = 0; c < result.cols; c++)
		{
			if ((c + r) % 2) result.at<float>(r, c) *= -1;
		}
	}

	//  An important step    Convert the result to  CV_8u 
	result.convertTo(result, CV_8UC1, 1.0, 0);
	result = result(Rect(0, 0, image.cols, image.rows)).clone();
	imshow(" Image after custom filtering ", result);
	waitKey(0);

	return 0;
}
//  Callback function 

#endif