在C中使用FFTW的高通滤波器

我有一个关于FFT的问题。 我已经设法在C中使用FFTW向前和向后进行FFT。现在，我想应用高通滤波器进行边缘检测，我的一些消息来源说只是将幅度的中心归零。

这是我的输入图片http://sofzh.miximages.com/c/2wnxvfl.jpg

基本上我所做的是：

1. 正向FFT
2. 将输出转换为2D数组
3. 进行FFT移位
4. 当距离中心的距离为高度的25％时，将real和imag值设为0
5. 产生幅度
6. 进行反向FFT移位
7. 转换为1D数组
8. 做后向FFT。

这是原始幅度，处理幅度和结果

http://sofzh.miximages.com/c/aysx9s.png

有人可以帮助我，告诉我哪个部分是错的，以及如何在C中使用FFTW进行高通滤波

谢谢。

源代码：

unsigned char **FFT2(int width,int height, unsigned char **pixel, char line1[100],char line2[100], char line3[100],char filename[100]) { fftw_complex* in, * dft, * idft, * dft2; //fftw_complex tmp1,tmp2; fftw_plan plan_f,plan_i; int i,j,k,w,h,N,w2,h2; w = width; h = height; N = w*h; unsigned char **pixel_out; pixel_out = malloc(h*sizeof(unsigned char*)); for(i = 0 ; i<h;i++) pixel_out[i]=malloc(w*sizeof(unsigned char)); in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *N); dft = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *N); dft2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *N); idft = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *N); /*run forward FFT*/ plan_f = fftw_plan_dft_2d(w,h,in,dft,FFTW_FORWARD,FFTW_ESTIMATE); for(i = 0,k = 0 ; i < h ; i++) { for(j = 0 ; j < w ; j++,k++) { in[k][0] = pixel[i][j]; in[k][1] = 0.0; } } fftw_execute(plan_f); double maxReal = 0.0; for(i = 0 ; i  maxReal ? dft[i][0] : maxReal; printf("MAX REAL : %f\n",maxReal); /*fftshift*/ //convert to 2d double ***temp1; temp1 = malloc(h * sizeof (double**)); for (i = 0;i < h; i++){ temp1[i] = malloc(w*sizeof (double*)); for (j = 0; j < w; j++){ temp1[i][j] = malloc(2*sizeof(double)); } } double ***temp2; temp2 = malloc(h * sizeof (double**)); for (i = 0;i < h; i++){ temp2[i] = malloc(w*sizeof (double*)); for (j = 0; j < w; j++){ temp2[i][j] = malloc(2*sizeof(double)); } } for (i = 0;i < h; i++){ for (j = 0; j < w; j++){ temp1[i][j][0] = dft[i*w+j][0]; temp1[i][j][1] = dft[i*w+j][1]; } } int m2 = h/2; int n2 = w/2; //forward shifting for (i = 0; i < m2; i++) { for (k = 0; k < n2; k++) { double tmp13[2] = {temp1[i][k][0],temp1[i][k][1]}; temp1[i][k][0] = temp1[i+m2][k+n2][0]; temp1[i][k][1] = temp1[i+m2][k+n2][1]; temp1[i+m2][k+n2][0] = tmp13[0]; temp1[i+m2][k+n2][1] = tmp13[1]; double tmp24[2] = {temp1[i+m2][k][0],temp1[i+m2][k][1]}; temp1[i+m2][k][0] = temp1[i][k+n2][0]; temp1[i+m2][k][1] = temp1[i][k+n2][1]; temp1[i][k+n2][0] = tmp24[0]; temp1[i][k+n2][1] = tmp24[1]; } } //process for (i = 0;i < h; i++){ for (j = 0; j < w; j++){ if(distance_to_center(i,j,m2,n2) < 0.25*h) { temp1[i][j][0] = (double)0.0; temp1[i][j][1] = (double)0.0; } } } /* copy for magnitude */ for (i = 0;i < h; i++){ for (j = 0; j < w; j++){ temp2[i][j][0] = temp1[i][j][0]; temp2[i][j][1] = temp1[i][j][1]; } } //backward shifting for (i = 0; i < m2; i++) { for (k = 0; k < n2; k++) { double tmp13[2] = {temp1[i][k][0],temp1[i][k][1]}; temp1[i][k][0] = temp1[i+m2][k+n2][0]; temp1[i][k][1] = temp1[i+m2][k+n2][1]; temp1[i+m2][k+n2][0] = tmp13[0]; temp1[i+m2][k+n2][1] = tmp13[1]; double tmp24[2] = {temp1[i+m2][k][0],temp1[i+m2][k][1]}; temp1[i+m2][k][0] = temp1[i][k+n2][0]; temp1[i+m2][k][1] = temp1[i][k+n2][1]; temp1[i][k+n2][0] = tmp24[0]; temp1[i][k+n2][1] = tmp24[1]; } } //convert back to 1d for (i = 0;i < h; i++){ for (j = 0; j < w; j++){ dft[i*w+j][0] = temp1[i][j][0]; dft[i*w+j][1] = temp1[i][j][1]; dft2[i*w+j][0] = temp2[i][j][0]; dft2[i*w+j][1] = temp2[i][j][1]; } } /* magnitude */ double max = 0; double min = 0; double mag=0; for (i = 0, k = 1; i < h; i++){ for (j = 0; j < w; j++, k++){ mag = sqrt(pow(dft2[i*w+j][0],2) + pow(dft2[i*w+j][1],2)); if (max < mag) max = mag; } } double **magTemp; magTemp = malloc(h * sizeof (double*)); for (i = 0;i < h; i++){ magTemp[i] = malloc(w*sizeof (double)); } for(i = 0,k = 0 ; i < h ; i++) { for(j = 0 ; j < w ; j++,k++) { double mag = sqrt(pow(dft2[i*w+j][0],2) + pow(dft2[i*w+j][1],2)); mag = 255*(mag/max); //magTemp[i][j] = 255-mag; //Putih magTemp[i][j] = mag; //Item } } /* brightening magnitude*/ for(i = 0,k = 0 ; i < h ; i++) { for(j = 0 ; j < w ; j++,k++) { //double temp = magTemp[i][j]; double temp = (double)(255/(log(1+255)))*log(1+magTemp[i][j]); pixel_out[i][j] = (unsigned char)temp; } } generateImage(width,height,pixel_out,line1,line2,line3,filename,"magnitude"); /* backward fft */ plan_i = fftw_plan_dft_2d(w,h,dft,idft,FFTW_BACKWARD,FFTW_ESTIMATE); fftw_execute(plan_i); for(i = 0,k = 0 ; i < h ; i++) { for(j = 0 ; j < w ; j++,k++) { double temp = idft[i*w+j][0]/N; pixel_out[i][j] = (unsigned char)temp; //+ pixel[i][j]; } } generateImage(width,height,pixel_out,line1,line2,line3,filename,"backward"); return pixel_out; } 

编辑新的源代码

我在前移之前添加了这部分，结果也如预期的那样。

 //proses //create filter unsigned char **pixel_filter; pixel_filter = malloc(h*sizeof(unsigned char*)); for(i = 0 ; i<h;i++) pixel_filter[i]=malloc(w*sizeof(unsigned char)); for (i = 0;i < h; i++){ for (j = 0; j < w; j++){ if(distance_to_center(i,j,m2,n2) < 20) { pixel_filter[i][j] = 0; } else { pixel_filter[i][j] = 255; } } } generateImage(width,height,pixel_filter,line1,line2,line3,filename,"filter1"); for (i = 0; i < m2; i++) { for (k = 0; k < n2; k++) { unsigned char tmp13 = pixel_filter[i][k]; pixel_filter[i][k] = pixel_filter[i+m2][k+n2]; pixel_filter[i+m2][k+n2] = tmp13; unsigned char tmp24 = pixel_filter[i+m2][k]; pixel_filter[i+m2][k] = pixel_filter[i][k+n2]; pixel_filter[i][k+n2] = tmp24; } } generateImage(width,height,pixel_filter,line1,line2,line3,filename,"filter2"); for (i = 0;i < h; i++){ for (j = 0; j < w; j++){ temp1[i][j][0] *= pixel_filter[i][j]; temp1[i][j][1] *= pixel_filter[i][j]; } } 

Answer 1

你的总体想法是好的。 从输出中，很难判断程序中是否存在简单的会计问题，或者这是否是预期的结果。 尝试使用更多空白空间填充源图像，并过滤出频域中较小的区域。

作为旁注，在C中执行此操作显得非常痛苦。 这是Matlab中的等效实现。 不包括绘图，它是大约10行代码。 您也可以尝试使用Numerical Python（NumPy）。

 % Demonstrate frequency-domain image filtering in Matlab % Define the grid x = linspace(-1, 1, 1001); y = x; [X, Y] = meshgrid(x, y); % Make a square (source image) rect = (abs(X) < 0.1) & (abs(Y) < 0.1); % Compute the transform rect_hat = fft2(rect); % Make the high-pass filter R = sqrt(X.^2 + Y.^2); filt = (R > 0.05); % Apply the filter rect_hat_filtered = rect_hat .* ifftshift(filt); % Compute the inverse transform rect_filtered = ifft2(rect_hat_filtered); %% Plot everything figure(1) imagesc(rect); title('source'); axis square saveas(gcf, 'fig1.png'); figure(2) imagesc(abs(fftshift(rect_hat))); title('fft(source)'); axis square saveas(gcf, 'fig2.png'); figure(3) imagesc(filt); title('filter (frequency domain)'); axis square saveas(gcf, 'fig3.png'); figure(4) imagesc(fftshift(abs(rect_hat_filtered))); title('fft(source) .* filter'); axis square saveas(gcf, 'fig4.png'); figure(5) imagesc(abs(rect_filtered)) title('result'); axis square saveas(gcf, 'fig5.png'); 

源图片：

源图像的傅里叶变换：

filter：

将滤波器应用（乘以）源图像的傅里叶变换的结果：

采用逆变换得出最终结果：