High Pass Filter using FFTW in C

2019-06-01 02:06发布

站内文章 / C
36 0 0

问题:

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:

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:

Fourier transform of the source image:

The filter:

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

Taking the inverse transform gives the final result:



标签： c image-processing fft fftw magnitude
举报

收藏的人(0)

Ta的文章 更多文章
登录 后发表评论
0条评论
还没有人评论过~