使用C中的FFTW高通滤波器 [英] High Pass Filter using FFTW in C

查看：485
本文介绍了使用C中的FFTW高通滤波器的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！

问题描述

我有关于FFT的问题。我已经设法做FFT向前和向后的C.使用FFTW现在，我想申请高通滤波器的边缘检测，我的一些来源说，刚调零幅度的中心。

这是我的输入图像
http://i62.tinypic.com/2wnxvfl.jpg

基本上我做什么是：


  1. 转发FFT

  2. 输出转换为二维数组

  3. 请向前移位FFT

  4. 使real和imag值为0时，从中心的距离是高度
    5. 25％
  5. 生成幅度

  6. 请向后移FFT

  7. 转换成一维数组

  8. 请向后FFT。

这是原始幅值，处理后的幅度，其结果

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

有人可以帮助我，告诉我哪一部分是错误的，如何做高通使用C中的FFTW过滤。

感谢你。

来源$ C ​​$ C：

  unsigned char型** FFT2（INT宽度，高度INT，unsigned char型**像素，焦炭一号线[100]，2号线的char [100]，CHAR 3号线[100]，文件名字符[100]）
{
  fftw_complex *中，* DFT，IDFT *，* DFT2;  // fftw_complex TMP1，TMP2;
  fftw_plan plan_f，plan_i;
  INT I，J，K，W，H，N，W2，H2;  W =宽度;
  H =高度;
  N = W * H;  unsigned char型** pixel_out;
  pixel_out =的malloc（H *的sizeof（无符号字符*））;
  对于（i = 0; I＆LT; H，我++）
    pixel_out [I] =的malloc（W * sizeof的（无符号字符））;  在=（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）;  / *运行正向FFT * /  plan_f = fftw_plan_dft_2d（W，H，中，DFT，FFTW_FORWARD，FFTW_ESTIMATE）;
  对于（I = 0，K = 0; I＆下; H;我+ +）
  {
    为（J = 0; J＆LT;瓦; J ++，K ++）
    {
      在[K] [0] =像素[I] [J]。
      在[k]的[1] = 0.0;
    }
  }
  fftw_execute（plan_f）;
  双maxReal = 0.0;
  对于（i = 0; I＆LT; N;我++）
    maxReal = DFT [I] [0]≥ maxReal？ DFT [I] [0]：maxReal;  的printf（MAX REAL：％F \\ N，maxReal）;
  / * * fftshift /
  //转换为2D
  双皇冠temp1中;
  temp1中的malloc =（H *的sizeof（双**））;
  对于（i = 0; I＆LT; H，我++）{
      temp1中[I] =的malloc（W * sizeof的（双*））;
      为（J = 0; J＆LT;瓦; J ++）{
          temp1中[I] [J] =的malloc（2 * sizeof的（双））;
      }
  }  双皇冠TEMP2;
  TEMP2 =的malloc（H *的sizeof（双**））;
  对于（i = 0; I＆LT; H，我++）{
      TEMP2 [I] =的malloc（W * sizeof的（双*））;
      为（J = 0; J＆LT;瓦; J ++）{
          TEMP2 [I] [J] =的malloc（2 * sizeof的（双））;
      }
  }
  对于（i = 0; I＆LT; H，我++）{
      为（J = 0; J＆LT;瓦; J ++）{
          temp1中[I] [j]的[0] = DFT [我* W + J] [0];
          temp1中[I] [j]的[1] = DFT [我* W + J] [1];
      }
  }  INT 2 = H / 2;
  INT N2 = W / 2;  //前进换档
  对于（i = 0; I＆LT; M2;我++）
  {
     对于（K = 0; K＆LT; N2; k ++）
     {
          双tmp13 [2] = {的temp1 [I] [K] [0]，temp1中由[i] [k]的[1]};
          的temp1 [I] [K] [0] = temp1中第[i + 2] [K + N 2] [0];
          temp1中由[i] [k]的[1] = temp1中第[i + 2] [K + N 2] [1];
          temp1中[1 + 2] [K + N 2] [0] = tmp13 [0];
          temp1中[1 + 2] [K + N 2] [1] = tmp13 [1];          双tmp24 [2] = {的temp1 [1 + 2] [k]的[0]，temp1中[1 + 2] [k]的[1]};
          temp1中[1 + 2] [K] [0] = temp1目录[I] [K + N2] [0];
          temp1中[1 + 2] [k]的[1] = temp1目录[I] [K + N2] [1];
          的temp1 [I] [K + N2] [0] = tmp24 [0];
          的temp1 [I] [K + N2] [1] = tmp24 [1];
     }
  }  //处理  对于（i = 0; I＆LT; H，我++）{
      为（J = 0; J＆LT;瓦; J ++）{
          如果（distance_to_center（I，J，2，N2）下; 0.25 * H）
          {
            temp1中[I] [J] [0] =（双）0.0;
            temp1中[I] [j]的[1] =（双）0.0;
          }
      }
  }  / *复制的幅度* /
  对于（i = 0; I＆LT; H，我++）{
      为（J = 0; J＆LT;瓦; J ++）{
          TEMP2 [I] [J] [0] = temp1中[I] [J] [0];
          TEMP2 [I] [j]的[1] = temp1中[I] [J] [1];
      }
  }
  //转移落后
  对于（i = 0; I＆LT; M2;我++）
  {
     对于（K = 0; K＆LT; N2; k ++）
     {
          双tmp13 [2] = {的temp1 [I] [K] [0]，temp1中由[i] [k]的[1]};
          的temp1 [I] [K] [0] = temp1中第[i + 2] [K + N 2] [0];
          temp1中由[i] [k]的[1] = temp1中第[i + 2] [K + N 2] [1];
          temp1中[1 + 2] [K + N 2] [0] = tmp13 [0];
          temp1中[1 + 2] [K + N 2] [1] = tmp13 [1];          双tmp24 [2] = {的temp1 [1 + 2] [k]的[0]，temp1中[1 + 2] [k]的[1]};
          temp1中[1 + 2] [K] [0] = temp1目录[I] [K + N2] [0];
          temp1中[1 + 2] [k]的[1] = temp1目录[I] [K + N2] [1];
          的temp1 [I] [K + N2] [0] = tmp24 [0];
          的temp1 [I] [K + N2] [1] = tmp24 [1];
     }
  }  //转换回1D
  对于（i = 0; I＆LT; H，我++）{
      为（J = 0; J＆LT;瓦; J ++）{
          DFT [我* W + J] [0] = temp1中[I] [J] [0];          DFT [我* W + J] [1] = temp1中[I] [J] [1];
          DFT2 [我* W + J] [0] = TEMP2 [I] [J] [0];          DFT2 [我* W + J] [1] = TEMP2 [I] [J] [1];      }
  }
  / * *级/  双最大= 0;
  双分钟= 0;
  双MAG = 0;
  对于（I = 0，K = 1; I＆下; H;我++）{
      为（J = 0; J＆LT;瓦; J ++，K ++）{
          MAG =开方（POW（DFT2 [我* W + J] [0]，2）+ POW（DFT2 [我* W + J] [1]，2））;
          如果（MAX＆LT; MAG）
              MAX = MAG;
      }
  }
  双** magTemp;
  magTemp =的malloc（H *的sizeof（双*））;
  对于（i = 0; I＆LT; H，我++）{
      magTemp [I] =的malloc（W *的sizeof（双））;
  }  对于（I = 0，K = 0; I＆下; H;我+ +）
  {
    为（J = 0; J＆LT;瓦; J ++，K ++）
    {      双MAG =开方（POW（DFT2 [我* W + J] [0]，2）+ POW（DFT2 [我* W + J] [1]，2））;
      MAG = 255 *（MAG /最大）;
      // magTemp [I] [J] = 255-MAG; // Putih
      magTemp [I] [J] = MAG; //项目
    }
  }  / *增亮幅度* /  对于（I = 0，K = 0; I＆下; H;我+ +）
  {
    为（J = 0; J＆LT;瓦; J ++，K ++）
    {
      //双温= magTemp [I] [J]。
      双TEMP =（双）（255 /（日志（1 + 255）））*日志（1 + magTemp [I] [J]）;
      pixel_out [I] [J] =（unsigned char型）温度;    }
  }  generateImage（宽度，高度，pixel_out，行1，行2，3号线，文件名，大小）;
  / *落后FFT * /
  plan_i = fftw_plan_dft_2d（W，H，DFT，IDFT，FFTW_BACKWARD，FFTW_ESTIMATE）;
  fftw_execute（plan_i）;
  对于（I = 0，K = 0; I＆下; H;我+ +）
  {
    为（J = 0; J＆LT;瓦; J ++，K ++）
    {
      双TEMP = IDFT [I * W + J] [0] / N;
      pixel_out [I] [J] =（unsigned char型）温度; // +像素[I] [J]。    }
  }
  generateImage（宽度，高度，pixel_out，一号线，2号线，3号线，文件名，落后）;  返回pixel_out;
}

修改新源$ C ​​$ C

我向前移位之前添加这部分，其结果是如预期也

  //散文
  //创建过滤器
  unsigned char型** pixel_filter;
  pixel_filter =的malloc（H *的sizeof（无符号字符*））;
  对于（i = 0; I＆LT; H，我++）
    pixel_filter [I] =的malloc（W * sizeof的（无符号字符））;
  对于（i = 0; I＆LT; H，我++）{
      为（J = 0; J＆LT;瓦; J ++）{
          如果（distance_to_center（I，J，2，N2）小于20）
          {
            pixel_filter [I] [J] = 0;
          }
          其他
          {
            pixel_filter [I] [J] = 255;
          }
      }
  }
  generateImage（宽度，高度，pixel_filter，行1，行2，3号线，文件名，过滤器1）;
  对于（i = 0; I＆LT; M2;我++）
  {
     对于（K = 0; K＆LT; N2; k ++）
     {
          unsigned char型tmp13 = pixel_filter [I] [K];
          pixel_filter [I] [K] = pixel_filter [1 + 2] [K + N2];
          pixel_filter [I + M2] [K + N2 = tmp13;          unsigned char型tmp24 = pixel_filter [I + M2] [K];
          pixel_filter [1 + 2] [K] = pixel_filter [I] [K + N2];
          pixel_filter [I] [K + N2 = tmp24;
     }
  }
  generateImage（宽度，高度，pixel_filter，行1，行2，3号线，文件名，过滤器2）;
  对于（i = 0; I＆LT; H，我++）{
      为（J = 0; J＆LT;瓦; J ++）{
          temp1中[I] [j]的[0] * = pixel_filter [I] [J]。
          temp1中[I] [j]的[1] * = pixel_filter [I] [J]。
      }
  }


解决方案

您总的想法是确定。从输出中，这是很难说是否有直接将程序会计问题，或者这是否可能是预期的结果。尝试用更空的空间填充所述源图像，并在频域中筛选出较小的区域。

作为一个方面说明，在C这样做似乎令人难以置信的痛苦。这里是在Matlab等效实现。不包括绘图，这是大约10 $ C $行C。您也可以尝试数值的Python（numpy的）。

 ％演示在Matlab频域图像滤波％定义网格
X = linspace（-1，1，1001）;
Y = X;
[X，Y] = meshgrid（X，Y）;％做一个正方形（源图像）
矩形=（ABS（X）下0.1）及（ABS（Y）＆LT; 0.1）;％计算变换
rect_hat = FFT2（RECT）;％使高通滤波器
R =开方（十^ 2 + Y. ^ 2）;
FILT =（R GT; 0.05）;％应用过滤器
rect_hat_filtered = rect_hat * ifftshift（FILT）。％计算逆变换
rect_filtered = ifft2（rect_hat_filtered）;%%绘制一切图1）
于imagesc（RECT）;
标题（源）;
轴方
另存为（GCF，'fig1.png'）;图（之二）
于imagesc（ABS（fftshift（rect_hat）））;
标题（FFT（源）'）;
轴方
另存为（GCF，'fig2.png'）;图（之三）
于imagesc（FILT）;
标题（过滤器（频域）'）;
轴方
另存为（GCF，'fig3.png'）;图（4）
于imagesc（fftshift（ABS（rect_hat_filtered）））;
标题（FFT（源）*过滤器。'）;
轴方
另存为（GCF，'fig4.png'）;图（5）
于imagesc（ABS（rect_filtered））
标题（结果）;
轴方
另存为（GCF，'fig5.png'）;

源图像：

傅里叶变换源图像的：

过滤器：

施加（相乘）的滤波器的傅立叶的结果变换源图像的：

以逆变换给出的最终结果：

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

解决方案

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中的FFTW高通滤波器的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！

查看全文
相关文章
C/C++最新文章
热门教程
热门工具
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆