使用fftw3对图片进行fft2操作

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
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);

}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
#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；
}