OPENCV计算机视觉图像处理频域傅里叶 DFT 变换低通滤波逆变换IDFT

OpenCV计算机视觉图像频域傅里叶 DFT 变换低通滤波逆变换IDFT

实验室做图像的，经常用到这部分，为了检测屏幕，看过好多博客，试用过许多代码，这个算是我找到的比较好用的，也容易改。

傅里叶变换（DFT）变换

    Mat padded;                            //expand input image to optimal size
	int a = getOptimalDFTSize(I.rows);
	int b = getOptimalDFTSize(I.cols); // on the border add zero values
	copyMakeBorder(I, padded, 0, a - I.rows, 0, b - I.cols, BORDER_CONSTANT, Scalar::all(0));


	for (int i = 0; i < padded.rows; i++)//中心化操作
	{
		float *ptr = padded.ptr<float>(i);
		for (int j = 0; j < padded.cols; j++)	ptr[j] *= pow(-1, i + j);
	}

	Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };
	Mat complexI;
	merge(planes, 2, complexI);         

	dft(complexI, complexI);            // this way the result may fit in the source matrix
	complexI.convertTo(complexI, CV_32F);//转为浮点型数据*

处理DFT 变换的频谱图

    Mat spectrogram1;
	spectrogram1 = complexI.clone();
	split(spectrogram1, planes);//分离通道 
	planes[0] = Re(DFT(I),planes[1]=Im(DFT(I)));

	magnitude(planes[0], planes[1], planes[0]);
	Mat magnitudeImage1 = planes[0];
	//进行对数尺度变换
	magnitudeImage1 += Scalar::all(1);
	log(magnitudeImage1, magnitudeImage1);

	normalize(magnitudeImage1, magnitudeImage1, 0, 1, NORM_MINMAX);
	magnitudeImage1 = magnitudeImage1 * 255;

想显示频谱图的话需要归一化到 0-255

高斯低通滤波

Mat gaussianBlur(padded.size(), CV_32FC2);//GuussianBlur low filter

	float D0 = 2 * 100 * 100;//该参数关乎滤波器的截止频率
	for (int i = 0; i < padded.rows; i++)
	{
		float*p = gaussianBlur.ptr<float>(i);

		for (int j = 0; j < padded.cols; j++)
		{
			float d = pow(i - padded.rows / 2, 2) + pow(j - padded.cols / 2, 2);
			p[2 * j] = expf(-d / D0);
			p[2 * j + 1] = expf(-d / D0);
		}
	}
	multiply(complexI, gaussianBlur, gaussianBlur);//滤波器和图像矩阵相乘

处理低通滤波后的频谱图

如果你想看一下滤波前后的频谱图有什么区别的话可以这样做

    Mat spectrogram;
	spectrogram = gaussianBlur.clone();
	split(spectrogram, planes);//分离通道 planes[0] = Re(DFT(I),planes[1]=Im(DFT(I)));

	magnitude(planes[0], planes[1], planes[0]);
	Mat magnitudeImage = planes[0];
	//进行对数尺度变换
	magnitudeImage += Scalar::all(1);
	log(magnitudeImage, magnitudeImage);

	normalize(magnitudeImage, magnitudeImage, 0, 255, NORM_MINMAX);
	//magnitudeImage = magnitudeImage * 255;
	imwrite("\\\\Mac\\Home\\Downloads\\测试\\2低通滤波后的频谱图.bmp", magnitudeImage);

	//imwrite((*sendImg1).sNowTimeFileDir + to_string(cameraFlag) + to_string(pictureFlag) + "gaussian_spectrogram.png", magnitudeImage);

逆变换 IDFT

    Mat iDft[] = { Mat::zeros(padded.size(),CV_32F),Mat::zeros(padded.size(),CV_32F) };
	idft(gaussianBlur, gaussianBlur, DFT_SCALE | DFT_REAL_OUTPUT); // IDFT进行DFT逆变换
	split(gaussianBlur, iDft);
	magnitude(iDft[0], iDft[1], iDft[0]);
	// normalize(iDft[0], iDft[0], 0, 1, CV_MINMAX);
	//iDft[0] = 255 * iDft[0];//去掉归一化



	Mat dst(iDft[0], Rect(0, 0, I.cols, I.rows));//之前扩过尺寸，现在截回来
	// Mat I(I, Rect(8, 5, I.cols - 10, I.rows - 10));
	blur(dst, dst, Size(7, 7));

完整版代码

void dftprocess(Mat &imgIn, Mat &imgOut, int cameraFlag, int pictureFlag)
{



	long t0 = GetTickCount();
	Mat I;
	cvtColor(imgIn, I, CV_BGR2GRAY);


	//Mat I = imread(filename, CV_LOAD_IMAGE_GRAYSCALE);
	I.convertTo(I, CV_32F);
	blur(I, I, Size(5, 5));
	Mat padded;                            //expand input image to optimal size
	int a = getOptimalDFTSize(I.rows);
	int b = getOptimalDFTSize(I.cols); // on the border add zero values
	copyMakeBorder(I, padded, 0, a - I.rows, 0, b - I.cols, BORDER_CONSTANT, Scalar::all(0));


	for (int i = 0; i < padded.rows; i++)//中心化操作
	{
		float *ptr = padded.ptr<float>(i);
		for (int j = 0; j < padded.cols; j++)	ptr[j] *= pow(-1, i + j);
	}

	Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };
	Mat complexI;
	merge(planes, 2, complexI);         // Add to the expanded another plane with zeros

	dft(complexI, complexI);            // this way the result may fit in the source matrix
	complexI.convertTo(complexI, CV_32F);//转为浮点型数据*


  //****************************写DFT变换后的频谱图*********************
	Mat spectrogram1;
	spectrogram1 = complexI.clone();
	split(spectrogram1, planes);//分离通道 planes[0] = Re(DFT(I),planes[1]=Im(DFT(I)));

	magnitude(planes[0], planes[1], planes[0]);
	Mat magnitudeImage1 = planes[0];
	//进行对数尺度变换
	magnitudeImage1 += Scalar::all(1);
	log(magnitudeImage1, magnitudeImage1);

	normalize(magnitudeImage1, magnitudeImage1, 0, 1, NORM_MINMAX);
	magnitudeImage1 = magnitudeImage1 * 255;

	//imwrite((*sendImg1).sNowTimeFileDir + to_string(cameraFlag) + to_string(pictureFlag) + "DFT_spectrogram.png", magnitudeImage1);
	imwrite("\\\\Mac\\Home\\Downloads\\测试\\1DFTimg.bmp", magnitudeImage1);


	/************************gaussian****************************/
	Mat gaussianBlur(padded.size(), CV_32FC2);//GuussianBlur low filter

	float D0 = 2 * 100000 * 100000;//低通截止频率设为100最合适
	for (int i = 0; i < padded.rows; i++)
	{
		float*p = gaussianBlur.ptr<float>(i);

		for (int j = 0; j < padded.cols; j++)
		{
			float d = pow(i - padded.rows / 2, 2) + pow(j - padded.cols / 2, 2);
			p[2 * j] = expf(-d / D0);
			p[2 * j + 1] = expf(-d / D0);
		}
	}
	multiply(complexI, gaussianBlur, gaussianBlur);




	//****************************写高斯滤波后的频谱图*********************
	Mat spectrogram;
	spectrogram = gaussianBlur.clone();
	split(spectrogram, planes);//分离通道 planes[0] = Re(DFT(I),planes[1]=Im(DFT(I)));

	magnitude(planes[0], planes[1], planes[0]);
	Mat magnitudeImage = planes[0];
	//进行对数尺度变换
	magnitudeImage += Scalar::all(1);
	log(magnitudeImage, magnitudeImage);

	normalize(magnitudeImage, magnitudeImage, 0, 255, NORM_MINMAX);
	//magnitudeImage = magnitudeImage * 255;
	imwrite("\\\\Mac\\Home\\Downloads\\测试\\2低通滤波后的频谱图.bmp", magnitudeImage);

	//imwrite((*sendImg1).sNowTimeFileDir + to_string(cameraFlag) + to_string(pictureFlag) + "gaussian_spectrogram.png", magnitudeImage);

	
	/**********************逆变换****************************/

	Mat iDft[] = { Mat::zeros(padded.size(),CV_32F),Mat::zeros(padded.size(),CV_32F) };
	idft(gaussianBlur, gaussianBlur, DFT_SCALE | DFT_REAL_OUTPUT); // IDFT进行DFT逆变换
	split(gaussianBlur, iDft);
	magnitude(iDft[0], iDft[1], iDft[0]);
	// normalize(iDft[0], iDft[0], 0, 1, CV_MINMAX);
	//iDft[0] = 255 * iDft[0];//去掉归一化



	Mat dst(iDft[0], Rect(0, 0, I.cols, I.rows));
	// Mat I(I, Rect(8, 5, I.cols - 10, I.rows - 10));
	blur(dst, dst, Size(7, 7));
	//imwrite("C:\\dsm0730高斯测试\\对逆变换截取图.png", dst);

	imwrite("\\\\Mac\\Home\\Downloads\\测试\\3IDFTimg.bmp", dst);

行吧，第一篇博客诞生了。冈萨雷斯的数字图像处理太深奥，我就是个菜菜。实验室项目中多次用到频域变换，干脆我就写个博客吧。欢迎留言