Opencv + FFTW3 图象高斯高低通滤波

Myth丶恋晨 2022-05-15 01:50 508阅读 0赞

本文转自 [https://blog.csdn.net/guoruijiushiwo/article/details/77411491][https_blog.csdn.net_guoruijiushiwo_article_details_77411491]

1：LPF

![Center][]

公式：out(i,j) = exp(-((i-M/2)^2+(j-N/2)^2)/2/sigma^2);

2：HPF:

![Center 1][]

公式：out(i,j) = 1 - exp(-((i-M/2)^2+(j-N/2)^2)/2/sigma^2);

代码：

void HPF_Fliter(vector<double>&data, int cols, int rows, double gamma)
    {
    #if 1
    	float gamma22 = 2 * gamma*gamma;
    	float temp = 0.0;
    	int halfRows = cvRound(rows / 2);
    	int halfCols = cvRound(cols / 2);
    	char debug[12] = { 0 };
    	for (int j = 1; j <= rows; j++)
    	{
    		for (int i = 1; i <= cols; i++)
    		{
    			temp = (pow(j - halfRows, 2) + pow(i - halfCols, 2));
    			data.push_back(1.0 - exp(-temp / gamma22));
    		}
    	}
    #else
    	
    	fstream readFile("C:\\Users\\Tony\\Desktop\\debug.txt", ios::in);
    	string temp;
    	int lineCount = 0;
    	while (getline(readFile, temp))
    	{
    		vector<string>dest;
    		split(temp, "\t", dest);
    		lineCount++;
    		if (dest.size() > 1)
    		{
    			for (int i = 0; i < dest.size(); i++)
    			{
    				//printf_s("%s\tn", dest.at(i));
    				data.push_back(atof(dest.at(i).c_str()));
     
    			}
    		}
     
    		//printf_s("\n");
    	}
    #endif
    }

//生成频谱模板
    void test()
    {
    	int  i;
    	fftw_complex*din, *out;
    	fftw_plan p,backp;
     
    	Mat image = imread("C:\\Users\\Tony\\Desktop\\Image\\lenaCopy.bmp", CV_LOAD_IMAGE_GRAYSCALE);
    	din = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*image.cols*image.rows);
    	out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*image.cols*image.rows);
     
     
    	Mat floatImage(image.size(), CV_32FC1);
    	image.convertTo(floatImage, CV_32FC1);
    	
     
    	for (size_t j = 0; j < image.rows; j++)
    	{
    		for (size_t i = 0; i < image.cols; i++)
    		{
    			din[i + j*image.cols][0] = *floatImage.ptr<float>(j, i);
    			din[i + j*image.cols][1] = 0;
    		}
    	}
    	//forward fft
    	p = fftw_plan_dft_2d(image.rows, image.cols, din, out, FFTW_FORWARD, FFTW_ESTIMATE);
    	fftw_execute(p);
     
    	Mat Resource(image.size(), CV_32FC1);
    	Mat Imsource(image.size(), CV_32FC1);
    	for (size_t j = 0; j < image.rows; j++)
    	{
    		for (size_t i = 0; i < image.cols; i++)
    		{
    			*Resource.ptr<float>(j, i) = out[i + j*image.cols][0];//实部
    			*Imsource.ptr<float>(j, i) = out[i + j*image.cols][1];//虚部
    		}
    	}
    	
    	Mat Redest(image.size(), CV_32FC1);
    	Mat Imdest(image.size(), CV_32FC1);
    	fftshift(Resource, Redest);
    	fftshift(Imsource, Imdest);
    	
    	Mat absImage(image.size(), CV_32FC1);
    	for (size_t j = 0; j < image.rows; j++)
    	{
    		for (size_t i = 0; i < image.cols; i++)
    		{
    			*absImage.ptr<float>(j, i) = sqrt(pow(*Redest.ptr<float>(j, i), 2) + pow(*Imdest.ptr<float>(j, i), 2));
    		}
    	}
    	
    	
    	absImage += Scalar::all(1);
    	Mat logImage(image.size(), CV_32FC1);
    	log(absImage, logImage);
    	//SaveData(logImage);
    	double min = 0.0, max = 0.0;
    	minMaxLoc(logImage, &min, &max);
    	double scale = 255/(max-min);
    	double shift = -min*scale;
    	Mat myImage(image.size(), CV_8UC1);
    	convertScaleAbs(logImage, myImage, scale, shift);
    	imwrite("C:\\Users\\Tony\\Desktop\\lenaCopy.bmp", myImage);
    }

高低通滤波：

void testOne()
    {
    	int  i;
    	fftw_complex*din, *out;
    	fftw_plan p, backp;
     
    	Mat image = imread("C:\\Users\\Tony\\Desktop\\blur.bmp", CV_LOAD_IMAGE_GRAYSCALE);
    	din = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*image.cols*image.rows);
    	out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*image.cols*image.rows);
     
     
    	Mat floatImage(image.size(), CV_32FC1);
    	image.convertTo(floatImage, CV_32FC1);
    	for (size_t j = 0; j < image.rows; j++)
    	{
    		for (size_t i = 0; i < image.cols; i++)
    		{
    			din[i + j*image.cols][0] = *floatImage.ptr<float>(j, i);
    			din[i + j*image.cols][1] = 0.0;
    		}
    	}
    	//forward fft
    	p = fftw_plan_dft_2d(image.rows, image.cols, din, out, FFTW_FORWARD, FFTW_ESTIMATE);
    	fftw_execute(p);
     
    	
    	//图像频谱中心化
    	Mat Resource(image.size(), CV_32FC1);
    	Mat Imsource(image.size(), CV_32FC1);//虚部
    	Mat Redest(image.size(), CV_32FC1);
    	Mat Imdest(image.size(), CV_32FC1);//虚部
    	for (size_t j = 0; j < image.rows; j++)
    	{
    		for (size_t i = 0; i < image.cols; i++)
    		{
    			*Resource.ptr<float>(j, i) = out[i + j*image.cols][0];//实部
    			*Imsource.ptr<float>(j, i) = out[i + j*image.cols][1];//虚部
    		}
    	}
     
     
    	fftshift(Resource, Redest);
    	fftshift(Imsource, Imdest);//虚部
     
     
     
    	//图像频谱乘以不通的数据,构建滤波
    	vector<double>data;
    	HPF_Fliter(data, image.cols, image.rows, 10);
    	//LPF_Fliter(data, image.cols, image.rows, 100);
    	Mat filterRe(image.size(), CV_32FC1);
    	Mat filterIm(image.size(), CV_32FC1);//虚部
    	for (size_t j = 0; j < image.rows; j++)
    	{
    		for (size_t i = 0; i < image.cols; i++)
    		{
    			*filterRe.ptr<float>(j, i) = *Redest.ptr<float>(j, i) * data.at(j*image.cols + i);//实部
    			*filterIm.ptr<float>(j, i) = *Imdest.ptr<float>(j, i) * data.at(j*image.cols + i);//虚部
    		 }
    	}
     
     
     
     
    	//图像频谱逆中心化
    	Mat NotCenterRedest(image.size(), CV_32FC1);
    	Mat NotCenterImdest(image.size(), CV_32FC1);//虚部
    	ifftShift(filterRe, NotCenterRedest);
    	ifftShift(filterIm, NotCenterImdest);//虚部
     
     
    	//重新构造复数，傅里叶反变换
    	fftw_complex *spectrum = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*image.cols*image.rows);
    	for (size_t j = 0; j < image.rows; j++)
    	{
    		for (size_t i = 0; i < image.cols; i++)
    		{
    			spectrum[j*image.cols + i][0] = *NotCenterRedest.ptr<float>(j, i);//实部
    			spectrum[j*image.cols + i][1] = *NotCenterImdest.ptr<float>(j, i);//虚部
    		}
    	}
     
    	
     
    	//傅里叶反变换
    	fftw_complex *result = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*image.cols*image.rows);
    	backp = fftw_plan_dft_2d(image.rows, image.cols, spectrum, result, FFTW_BACKWARD, FFTW_ESTIMATE);
    	fftw_execute(backp);
    	int size = image.cols*image.rows;
     
    	Mat IvsResource(image.size(), CV_32FC1, Scalar::all(0));
    	//Mat IvsImsource(image.size(), CV_32FC1, Scalar::all(0));//虚部
    	//Mat ifftImage(image.size(), CV_32FC1, Scalar::all(0));
    	Mat myImage(image.size(), CV_8UC1);
    	for (size_t j = 0; j < image.rows; j++)
    	{
    		for (size_t i = 0; i < image.cols; i++)
    		{
    			*IvsResource.ptr<float>(j, i) = result[i + j*image.cols][0] / size;//实部
    			//*IvsImsource.ptr<float>(j, i) = result[i + j*image.cols][1] / size;//虚部
    			//*ifftImage.ptr<float>(j, i) = sqrt(pow(result[i + j*image.cols][0], 2) + pow(result[i + j*image.cols][1], 2));
    			*myImage.ptr<byte>(j, i) = cvRound(abs(*IvsResource.ptr<float>(j, i)));//实部
    		}
    	}
    	//SaveData(IvsResource);
    	imwrite("C:\\Users\\Tony\\Desktop\\Copy.bmp", myImage);
     
    	Mat gammaImage(image.size(), CV_8UC1);
    	GammaCorrection(myImage, gammaImage, 0.3);
     
    	//释放资源
    	fftw_destroy_plan(p);
    	fftw_destroy_plan(backp);
    	fftw_free(din);
    	fftw_free(out);
    	fftw_free(result);
    }

[https_blog.csdn.net_guoruijiushiwo_article_details_77411491]: https://blog.csdn.net/guoruijiushiwo/article/details/77411491
[Center]: /images/20220515/c59313944e664319b2a451c3a2b5287c.png
[Center 1]: /images/20220515/125f0876dc7543838194e6b48e9ce051.png