High Pass Filter using FFTW in C

2019-06-01 01:30发布

站内问答 / C
2015 1 6

I have a question regarding FFT. I already manage to do FFT forward and backward using FFTW in C. Now, I want to apply high pass filter for edge detection, some of my source said that just zeroing the centre of the magnitude.

This is my input image http://i62.tinypic.com/2wnxvfl.jpg

Basically what I do are :

  1. Forward FFT
  2. Convert the output to 2D array
  3. Do forward FFT shifting
  4. Make the real and imag value to 0 when the distance from the centre is 25% of the height
  5. Generate the magnitude
  6. Do backward FFT shifting
  7. Convert into 1D array
  8. Do Backward FFT.

This is the original magnitude, the processed magnitude, and the result

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

can someone help me, to tell me which part is wrong and how to do the high pass filtering using FFTW in C.

Thank You.

The Source Code:

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 new source code

I add this part before the forward shifting, the result is as expected also.

//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];
      }
  }

1条回答
男人必须洒脱
2楼-- · 2019-06-01 02:19

Your general idea is OK. From the output, it's hard to tell whether there's simply an accounting problem in your program, or whether this is perhaps the expected result. Try padding the source image with much more empty space, and filter out a smaller area in the frequency domain.

As a side note, doing this in C appears incredibly painful. Here is an equivalent implementation in Matlab. Not including plotting, it's around 10 lines of code. You might also try 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');

The source image: enter image description here

Fourier transform of the source image: enter image description here

The filter: enter image description here

Result of applying (multiplying) the filter with the fourier transform of the source image: enter image description here

Taking the inverse transform gives the final result: enter image description here

查看更多
0人赞 添加讨论(0) 举报
登录 后发表回答
相关问题
查看全部
相关文章
查看全部
收藏的人(6)