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

最新推荐文章于 2023-06-28 12:51:01 发布
原创 最新推荐文章于 2023-06-28 12:51:01 发布 · 2.3k 阅读 · 5 ·
CC 4.0 BY-SA版权
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。

1：LPF
公式：out(i,j) = exp(-((i-M/2)^2+(j-N/2)^2)/2/sigma^2);
代码：
         
void LPF_Fliter(vector<double>&data, int cols, int rows, double gamma)
{
	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(exp(-temp / gamma22));
		}
	}
}

2：HPF:


公式：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);
}


FFTW图像处理之频域滤波和频域分析
深之JohnChen的专栏
05-24 458
重用FFTW计划(plan)避免重复创建开销。对于大图像，考虑使用多线程FFTW。：应用汉宁窗等可以减少频谱泄漏。：在显示频谱前，可以使用。对于实数输入图像，使用。
DFT的matlab源代码-filter:matlab中filter函数的C语言实现，利用了fft替代了卷积
05-26
DFT的matlab源代码filter.c 这是matlab过滤功能的C实现 fft.c 这将提取fftwf_plan_dft_1d接口 依存关系 fft.c需要库。 测试用例 filter.c的性能优于matlab。 确切的测试用例敬请期待。
fftw库 vs2019_FFTW vs. OpenCV cvDFT
weixin_35823441的博客
02-15 470
I have used both (OpenCV and FFTW) and you can expect FFTW to be faster than the simpler implementation in OpenCV (how much depends a lot on your processor and image sizes of course). However, if you ...
matlab z变换离散化_傅里叶变换在MATLAB中的应用(频域滤波）
weixin_39897449的博客
11-19 435
图像滤波分为空间域滤波和频域滤波，空间滤波的内容见本人的另一篇文章：清逸：MATLAB中的图像变换之线性空间滤波​zhuanlan.zhihu.com本文主要讲述如何在MATLAB中实现频域滤波，那么，怎么实现呢，我们这里讲的所有的滤波都是通过傅里叶变换在频域中实现的，所有这部分和傅里叶变换渊源很深，至于傅里叶变换本身，我自己也不能解释的很清楚，我们只讲他如何在matlab中实现和应用。深入了解傅...
fftw3图片傅里叶变换_开发者都在使用的OpenCV FFT傅里叶检测上传到服务器上的模糊图片...
weixin_39828847的博客
11-20 453
FFT 快速傅里叶变换傅里叶变换快速傅里叶变换(Fast Fourier Transform，FFT)是一种可在 O(nlogn)时间内完成的离散傅里叶变换(Discrete Fourier transform，DFT)算法。具体什么是傅里叶变换，这里不再讲述，有兴趣的小伙伴们可以自行百度。在计算机视觉方面，我们经常将FFT视为代表两个域中的图像的图像处理工具：傅立叶(即频域)域空间域因...
fftw3图片傅里叶变换_使用傅里叶变换进行图像边缘检测
weixin_39820185的博客
11-28 458
今天我们介绍通过傅里叶变换求得图像的边缘什么是傅立叶变换？简单来说，傅里叶变换是将输入的信号分解成指定样式的构造块。例如，首先通过叠加具有不同频率的两个或更多个正弦函数而生成信号f(x)，之后，仅查看f(x)的图像缺无法了解使用哪种或多少原始函数来生成f(x)。这就是傅立叶变换最神奇的地方。将f(x)函数通过一个傅立叶变换器，我们就可以得到一个新的函数F(x)。F(x)的是最初生成f(x)函数的频...
基于OpenCV和fftw3库的MOSSE目标跟踪代码
使用Visual Studio C++开发MOSSE代码时，需要正确配置OpenCV库和fftw3库。这些库必须在编译前安装到系统中，并在项目中正确地包含它们的头文件和库文件。当MOSSE代码被封装为DLL（动态链接库）时，意味着可以被其他...
【图像处理高效实现】：FFTW在图像分析中的案例研究
随后深入探讨了FFTW的核心特性，如多线程优化技术和计划（Plans）对FFT加速的作用，以及在图像处理中的应用方法，包括频域分析、图像去噪、边缘检测和图像压缩等。最后，分析了FFTW在大规模图像处理
C++图像处理：小波变换与复原滤波技术
在C++中实现图像的小波变换、滤波和复原，需要涉及到诸如OpenCV、FFTW（快速傅里叶变换库）等第三方库。以下是一个简单的代码示例框架，用于说明如何利用C++和相关库进行小波变换的操作： ```cpp #include <opencv2...
#include <fftw3.h> #include <opencv2/opencv.hpp> #include <cmath> // 1. 读取灰度图像 cv::Mat readGrayscale(const char* path) { cv::Mat img = cv::imread(path, cv::IMREAD_GRAYSCALE); if(img.empty()) throw std::runtime_error("Failed to load image"); return img; } // 2. FFT变换 (返回复数矩阵) fftw_complex* imageToFFT(const cv::Mat& img) { int w = img.cols, h = img.rows; // 分配内存 fftw_complex* in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * w * h); fftw_complex* out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * w * h); // 初始化实数转复数 for(int i=0; i<h; i++) { uchar* p = img.ptr<uchar>(i); for(int j=0; j<w; j++) { in[i*w+j][0] = p[j] / 255.0; // 实部 [0,1] in[i*w+j][1] = 0; // 虚部 } } // 创建FFT计划并执行 fftw_plan plan = fftw_plan_dft_2d(h, w, in, out, FFTW_FORWARD, FFTW_ESTIMATE); fftw_execute(plan); fftw_destroy_plan(plan); fftw_free(in); return out; } // 3. 频域可视化（幅度谱） cv::Mat visualizeSpectrum(fftw_complex* freqData, int w, int h) { cv::Mat mag(h, w, CV_32F); for(int i=0; i<h; i++) { float* p = mag.ptr<float>(i); for(int j=0; j<w; j++) { double re = freqData[i*w+j][0]; double im = freqData[i*w+j][1]; p[j] = log(1 + sqrt(re*re + im*im)); // 对数增强 } } // 归一化并转换类型 cv::normalize(mag, mag, 0, 1, cv::NORM_MINMAX); cv::Mat display; mag.convertTo(display, CV_8U, 255); // 中心化 (FFTShift) int cx = w/2, cy = h/2; cv::Mat q0(display, cv::Rect(0, 0, cx, cy)); // 左上 cv::Mat q1(display, cv::Rect(cx, 0, cx, cy)); // 右上 cv::Mat q2(display, cv::Rect(0, cy, cx, cy)); // 左下 cv::Mat q3(display, cv::Rect(cx, cy, cx, cy)); // 右下 // 交换象限 cv::Mat tmp; q0.copyTo(tmp); q3.copyTo(q0); tmp.copyTo(q3); q1.copyTo(tmp); q2.copyTo(q1); tmp.copyTo(q2); return display; } // 4. 低通滤波（高斯滤波器） void gaussianLowPass(fftw_complex* freqData, int w, int h, double cutoff) { int cx = w/2, cy = h/2; for(int y=0; y<h; y++) { for(int x=0; x<w; x++) { // 计算距离中心点距离 double dx = (x - cx) * 1.0 / cx; double dy = (y - cy) * 1.0 / cy; double d = sqrt(dx*dx + dy*dy); // 应用高斯函数 double factor = exp(-d*d / (2*cutoff*cutoff)); freqData[y*w+x][0] *= factor; freqData[y*w+x][1] *= factor; } } } // 5. IFFT逆变换 cv::Mat ifftToImage(fftw_complex* freqData, int w, int h) { fftw_complex* out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*w*h); fftw_plan plan = fftw_plan_dft_2d(h, w, freqData, out, FFTW_BACKWARD, FFTW_ESTIMATE); fftw_execute(plan); // 转换为灰度图像 cv::Mat result(h, w, CV_8U); for(int i=0; i<h; i++) { uchar* p = result.ptr<uchar>(i); for(int j=0; j<w; j++) { // 取实数部分并缩放到[0,255] double val = out[i*w+j][0]/(w*h); // FFTW未归一化 p[j] = cv::saturate_cast<uchar>(val*255); } } fftw_destroy_plan(plan); fftw_free(out); return result; } // 6. 计算MSE和PSNR void evaluate(const cv::Mat& orig, const cv::Mat& recon, double& mse, double& psnr) { mse = 0; const int N = orig.total(); for(int i=0; i<orig.rows; i++) { const uchar* p1 = orig.ptr<uchar>(i); const uchar* p2 = recon.ptr<uchar>(i); for(int j=0; j<orig.cols; j++) { double diff = p1[j] - p2[j]; mse += diff*diff; } } mse /= N; if(mse < 1e-10) psnr = 100; // 防止除0 else psnr = 10 * log10(255*255 / mse); } int main() { // 1. 读取灰度图像 cv::Mat img = readGrayscale("input.jpg"); // 2. 执行FFT fftw_complex* freq = imageToFFT(img); // 3. 可视化频域 cv::Mat spectrum = visualizeSpectrum(freq, img.cols, img.rows); cv::imwrite("spectrum.jpg", spectrum); // 4. 低通滤波 (0.2为截止频率参数) gaussianLowPass(freq, img.cols, img.rows, 0.2); // 5. 逆变换重建图像 cv::Mat reconstructed = ifftToImage(freq, img.cols, img.rows); cv::imwrite("reconstructed.jpg", reconstructed); // 6. 计算质量指标 double mse, psnr; evaluate(img, reconstructed, mse, psnr); printf("MSE: %.2f\tPSNR: %.2f dB\n", mse, psnr); fftw_free(freq); return 0; }
最新发布
06-05
但FFTW要求行优先，而OpenCV也是行优先，所以可以直接按行复制）步骤5：创建FFT计划并执行步骤6：将频域数据移位（将零频移到中心）以便滤波步骤7：构建低通滤波器（例如高斯低通滤波器）步骤8：将滤波器应用于频域...
1. 使用 C/C 编写代码，读取灰度图像。 2. 对图像进行 FFT 运算，生成频域复数矩阵并可视化。 3. 对频域数据直接进行 IFFT，重建原始图像，计算重建图像与原图的均 方误差（MSE）或峰值信噪比（PSNR），验证 FFT/IFFT 的可逆性。 4. 对频域数据进行低通滤波，再做 IFFT，比较滤波后图像与原图像的均 方误差（MSE）或峰值信噪比（PSNR）。
06-05
计算原始图像和滤波后图像的MSE和PSNR由于用户要求使用C/C++，我们将采用以下方案：-使用FFTW进行FFT和IFFT-使用OpenCV进行图像读取、显示和保存（方便且常用）-如果不允许使用OpenCV，则需要自己读取图像文件（例如...
fftw
trentczp的博客
06-19 298
fftw使用例子：https://blog.youkuaiyun.com/grafx/article/details/38750107
fftw3图片傅里叶变换_数字图像处理（二）—— 离散非周期傅里叶变换
weixin_39960710的博客
11-28 421
在数字图像处理中碰到的信号基本上全部都是离散非周期信号，所以我单独把这部分拿出来，介绍一下傅里叶变换是怎么应用在图像领域中的。1. 一维离散非周期信号的傅里叶变换设 为离散非周期有限长的信号，长度为 。那么，需构造周期信号 ，其周期长度为 ，且 ，原信号不足的之处用 补齐。则在某周期内 。信号 的傅里叶变换为 （1）。逆傅里叶变换为 （2）。1.1 离散傅里叶变换本质将傅里叶...
fftw3图片傅里叶变换_DFT-去除图片周期性噪声
weixin_39729272的博客
11-27 399
DFT：离散傅里叶变换要处理的图片如下周期性噪声图片处理结果如下# -*- coding: utf-8 -*- 参考链接OpenCV-Python-Tutorial[3]_人工智能_风吴痕的博客-优快云博客​blog.youkuaiyun.com...
OpenCV快速傅里叶变换(FFT)用于图像和视频流的模糊检测
qq_50380073的博客
06-28 1335
OpenCV快速傅里叶变换(FFT)用于图像和视频流的模糊检测
FFTW3
majinguang2004的专栏
11-30 6327
<br />在做ITK的配准时发现了document的一个bug，itkcurvatureRegistrationFilter无法显示,在Nabble里看了一下，原来06年就有人提出这个问题了，不过在10年11-24这个问题才解决。解决的方法是：在加上预编译#ifndef USE_FFTWF #define USE_FFTWF<br />#ifndef USE_FFTWD  #define USE_FFTWD。但是如果你的电脑里没有fftw3（世界上最快的傅立叶变化）那么这个filter还是无法进行，下面的
fftw的使用
lightblueme的专栏
04-20 988
fftw_complex in[1024], out[1024]; fftw_plan p; p=fftw_plan_dft_1d(1024,in,out,FFTW_FORWARD,FFTW_MEASURE); for(int i=0;i <1024;i ++) { in[i][0]=data[i]; in[i][1]=0.0;
FFTW的使用方法 - 实数输入的数字滤波示例
weixin_38713073的博客
11-05 7613
主要记录FFTW的安装过程，FFTW的简单示例代码。并实现一个输入和输出为一维实数的数字滤波器
对信号进行滤波处理
jk_101的博客
05-25 8981
首先，建立三个模拟信号，它们的频率分别为5、20、50Hz，并将这个三个信号相加合成一个信号，即： y1 = sin(2*pi*5*t); y2 = cos(2*pi*20*t); y3 = sin(2*pi*50*t); %%合并信号 y = y1 + y2 + y3; plot(y); 如图所示： 对该信号进行傅里叶变换得到频谱图得： ...
评论 6
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值