高通濾波器使用FFTW在C

Question

我有一個關於FFT的問題。我已經設法使用C中的FFTW進行FFT向前和向後的FFT。現在，我想應用高通濾波器進行邊緣檢測，我的一些消息來源說，只是將幅度的中心置零。

這是我的輸入圖像 http://i62.tinypic.com/2wnxvfl.jpg

基本上我做什麼都：

1. 正向FFT
2. 轉換輸出到二維數組
3. 做正向FFT移
4. 充分利用當距離中心的距離爲高度的25％時，真實和成像值爲0
5. 生成大小
6. 執行向後移位FFT
7. 轉換成一維數組
8. 執行向後FFT。

這是原來的大小，處理後的幅度，結果

http://i58.tinypic.com/aysx9s.png

有人可以幫助我，告訴我哪一部分是錯誤的，如何做高通使用FFTW濾波在C.

謝謝。

的源代碼：

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 < N ; i++) 
    maxReal = dft[i][0] > 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; 
} 

EDIT新的源代碼

我向前移位前添加這部分，其結果是作爲也在意料之中。

//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行代碼。你也可以嘗試數值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'); 

源圖像：

傅立葉變換的源圖像：

濾波器：

施加（乘以）過濾器與所述傅立葉結果源圖像的變換：

以逆變換給出了最終結果：