使用fftw3对图片进行fft2操作

81 阅读 0 评论 54 点赞

我是靠谱客的博主结实衬衫，最近开发中收集的这篇文章主要介绍使用fftw3对图片进行fft2操作，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

在matlab中有个非常好用的二维傅里叶变换函数fft2,c++没有只能自己重写但是时间效率远远不如matlab,可以通过调用fftw3库来进行fft2函数的重写。

64位编译好的fftw3库：链接：https://pan.baidu.com/s/1-0LYnkh0pNsZCDewk3RS9w
提取码：zvgz

记得使用fftw3库之前要把环境配好，include、lib都要配好。

fft2方法：

输入：自定义的复数矩阵src,图片长宽nr和nc以及标识符inverse，当inverse=1时函数为二维傅里叶变换，当inverse=0时函数为二维傅里叶逆变换即为ifft2

输出：傅里叶变换后的复数矩阵dst

void fft2(complex1 **src, complex1 **dst, int nr, int nc, bool inverse)
{
	fftw_complex *signal_img = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*nr*nc);
	fftw_complex *out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*nr*nc);


	for (int i = 0; i < nr; i++)
	{
		for (int j = 0; j < nc; j++)
		{
			signal_img[i*nc + j][0] = src[i][j].r;
			signal_img[i*nc + j][1] = src[i][j].i;

		}
	}
	fftw_plan signal_forward_plan;
	if (inverse)
	{
		signal_forward_plan = fftw_plan_dft_2d(nr, nc, signal_img, out, FFTW_FORWARD, FFTW_ESTIMATE);
		fftw_execute(signal_forward_plan);
		
		for (int i = 0; i < nr; i++)
		{
			for (int j = 0; j < nc; j++)
			{
				dst[i][j] = complex1(out[i*nc + j][0], out[i*nc + j][1]);
			}
		}
	}
	else
	{
		signal_forward_plan = fftw_plan_dft_2d(nr, nc, signal_img, out, FFTW_BACKWARD, FFTW_ESTIMATE);
		fftw_execute(signal_forward_plan);
		
		for (int i = 0; i < nr; i++)
		{
			for (int j = 0; j < nc; j++)
			{
				dst[i][j] = complex1(out[i*nc + j][0] / (double)(nr*nc), out[i*nc + j][1] / (double)(nr*nc));
			}
		}
		
	}
	fftw_destroy_plan(signal_forward_plan);
	fftw_free(signal_img);
	fftw_free(out);

}

程序主函数

#include <iostream>
#include "opencv2/core/core.hpp"  
#include "opencv2/highgui/highgui.hpp"  
#include "opencv2/imgproc/imgproc.hpp"  
#include<math.h>
#include<fftw3.h>
struct complex1 {
	double r;
	double i;

	complex1(double real = 0, double imag = 0)
	{
		r = real;
		i = imag;
	}
	double abs() { return sqrt(r*r + i * i); }
	complex1 operator-(const complex1 &o) { return complex1(r - o.r, i - o.i); }
	complex1 operator*=(const complex1 &o) {
		return complex1(r * o.r - i * o.i, r*o.i + i * o.r);
	}
	complex1 &operator/=(const double o) {
		r /= o;
		i /= o;
		return *this;
	}
	complex1 &operator+=(const complex1 &o) {
		r += o.r;
		i += o.i;
		return *this;
	}
	complex1 &conj(const complex1 &o)
	{
		r = o.r;
		i = -o.i;
		return *this;
	}
	

};//定义一个复数
int main()
{
     cv::Mat Image_refernce = cv::imread("C:\Users\25353\Desktop\123\GF4_PMS_E119.8_N33.7_20160625_L1A0000117115_band_1.tif", cv::IMREAD_UNCHANGED);//输入一张图片，我这是16位深的图片。
    //图像的长宽
	int nr = Image_refernce.rows;
	int nc = Image_refernce.cols;
    complex1 **buf1ft = new complex1*[nr], **x = new complex1*[nr];
	for (int i = 0; i < nr; i++) {
		buf1ft[i] = new complex1[nc];
		x[i] = new complex1[nc];
	}
    for (int i = 0; i < nr; i++) {
		for (int j = 0; j < nc; j++) {
		
			x[i][j] = complex1(Image_refernce.at<ushort>(i, j), 0);//将像素值转换成复数complex1格式
		}
	}
    fft2(x, buf1ft, nr, nc, 1);//二维傅里叶变换
    for (int i = 0; i < nr; i++)
	{
		delete[] buf1ft[i];
		delete[] x[i];
	}
	delete[] buf1ft;
	delete[] x;
	
    return 0；
}