图像处理 之 同态滤波

 没什么好讲的  原图 --> YUV  --> fft2 -->FFTSHIFT  --> 指数滤波 -->IFFTSHIFT --> ifft2 --> YUV--> RGB24

上面看懂就看明白了，如果看不明白，就度娘一点点查....

#define printf #define PI (3.1415926535897932384626433)

#include <math.h>

typedef struct _plural {  double real;  double img; } plural;

plural EE( plural b1, plural b2) {  plural b3;  b3.real= b1.real*b2.real - b1.img*b2.img;  b3.img = b1.real*b2.img  + b1.img*b2.real;  return(b3); }

void iOrder(plural *pData,int levelN) {  plural tmp;  int k,j=1;

 //变变址运算，实现数据运算前的位倒序  for(int i=1;i< levelN;i++)  {

  if(i>j)   {    tmp=pData[j-1];    pData[j-1]=pData[i-1];    pData[i-1]=tmp;   }

  k=levelN/2;   while(k<j)   {    j=j-k;    k=k/2;   }   j=j+k;  }

}

void FFT_Init(int N, plural *&pFftW) {  int f = N;  int m = 0;  for(m=1;(f=f/2)!=1;m++) ;

 int i;  pFftW=new plural[m];  int len = 1;  for(i=0;i<m;i++)  {   len = len<<1;   int tmp = len>>1;

  //此处是把360度 划分成 2^m 份 -->此处只取一份   pFftW[i].real = cos(PI/tmp);   pFftW[i].img  =-sin(PI/tmp);

  printf("%4d %11.4f %11.4f ", i, pFftW[i].real, pFftW[i].img);

 }

}

void FFT_DFT(plural *pData, int N, plural *pFftW) {  int m = 0;

 int f = N;  for(m=1;(f=f/2)!=1;m++) ;

 // 倒序  iOrder(pData,N);

 plural v,w;  int len = 1;  for(int k = 0; k<m; k++) // 一共得计算多少次  {   len = len<<1;

  v.real=1.0;   v.img=0.0;

  int tmp = len>>1;

  //此处是把360度 划分成 2^m 份 -->此处只取一份   //w.real=cos(PI/tmp);   //w.img=-sin(PI/tmp);

  w.real =pFftW[k].real;   w.img  =pFftW[k].img;

  for (int j = 0; j<len/2; j++)   {    for (int i=0; i< N; i+=len)    {     //取出二对复数     plural tmp1 = pData[i+j];     plural tmp2 = pData[i+j + len/2];

    // *W(0,N)     // 复数相乘 (a+bj)(c+dj) = (ac-bd) + (bc+ad)j     plural t=EE(pData[i+j + len/2],v);

    pData[i+j + len/2].real = pData[i+j].real - t.real;     pData[i+j + len/2].img = pData[i+j].img - t.img;     pData[i+j].real = pData[i+j].real + t.real;     pData[i+j].img = pData[i+j].img + t.img;

   }    v=EE(v,w); //此处实际上是把Z坐标角度相加    v.real = v.real;    v.img  = v.img;   }

 } }

void IFFT_DFT(plural *pData, int N, plural *pFftW) {  int m = 0;

 int f = N;  for(m=1;(f=f>>1)!=1;m++) ;

 for (int i=0;i<N; i++)  {   pData[i].img = -pData[i].img;  }

 FFT_DFT(pData, N, pFftW);

 for(int k = 0; k<m; k++) // 一共得计算多少次  {

  for (int i=0; i< N; i++)   {    pData[i].real /= 2;    pData[i].img /= 2;   }  }

 for (int i=0;i<N; i++)  {   pData[i].img = -pData[i].img;  }

}

void fft2shifft(plural *pImgDataPlural, int PluralWidth, int PluralHeight) {

}

void ifft2shifft(plural *pImgDataPlural, int PluralWidth, int PluralHeight) {

}

void plural2D_transpose(plural *pImgDataPlural, int PluralWidth, int PluralHeight, plural *pDstDataPlural) {  // 转置  for (int i=0; i<PluralHeight; i++)  {   for (int j=0; j<PluralWidth; j++)   {    pDstDataPlural[j*PluralHeight+i]=pImgDataPlural[i*PluralWidth+j];   }  } }

void plural_print(plural *pData2D, int w, int h, char *pbuff = "") {  return;  printf(" %s ", pbuff);  for (int i=0; i<h; i++)  {   for (int j=0; j<w; j++)   {    printf("[ %7.3f  %7.3f]",pData2D[i*w + j].real, pData2D[i*w + j].img);   }   printf(" ");  } }

void image_print(double *pImgData, int w, int h) {  return;  printf(" ");  for (int i=0; i<h; i++)  {   for (int j=0; j<w; j++)   {    printf("[ %7.2f]",pImgData[i*w + j]);   }   printf(" ");  } }

void fft2(plural *pImgDataPlural, int PluralWidth, int PluralHeight, plural *pFftWidth, plural *pFftHeight, plural *pTmpDataPlural) {  // FFT2D  for (int i=0; i<PluralHeight; i++)  {   FFT_DFT(&pImgDataPlural[i*PluralWidth],PluralWidth, pFftWidth);  }

 plural_print(pImgDataPlural, PluralWidth,PluralHeight , "行FFT");

 // 转置  plural2D_transpose(pImgDataPlural, PluralWidth,PluralHeight, pTmpDataPlural);;

 plural_print(pTmpDataPlural, PluralHeight ,PluralWidth, "转置之后 ");

 for (int i=0; i<PluralWidth; i++)  {   FFT_DFT(&pTmpDataPlural[i*PluralHeight],PluralHeight, pFftHeight);  }  plural_print(pTmpDataPlural,PluralHeight, PluralWidth , "再做行FFT");

 plural2D_transpose(pTmpDataPlural,PluralHeight, PluralWidth, pImgDataPlural);

 fft2shifft(pImgDataPlural,PluralWidth, PluralHeight);

}

void ifft2(plural *pImgDataPlural, int PluralWidth, int PluralHeight, plural *pFftWidth, plural *pFftHeight, plural *pTmpDataPlural) {  ifft2shifft(pImgDataPlural,PluralWidth, PluralHeight);

 plural2D_transpose(pImgDataPlural,PluralWidth, PluralHeight, pTmpDataPlural);

 // IFFT2D  for (int i=0; i<PluralWidth; i++)  {   IFFT_DFT(&(pTmpDataPlural[i*PluralHeight]),PluralHeight, pFftHeight);  }

 plural_print(pTmpDataPlural,PluralHeight, PluralWidth , "再做行IFFT");

 // 转置  for (int i=0; i<PluralWidth; i++)  {   for (int j=0; j<PluralHeight; j++)   {    pImgDataPlural[j*PluralWidth+i]=pTmpDataPlural[i*PluralHeight+j];   }  }

 plural_print(pImgDataPlural,PluralWidth, PluralHeight , "再次转置");

 for (int i=0; i<PluralHeight; i++)  {   IFFT_DFT(&(pImgDataPlural[i*PluralWidth]),PluralWidth, pFftWidth);  }

 plural_print(pImgDataPlural, PluralWidth,PluralHeight , "行IFFT");

}

// double d0=9; // double Rh=3.5; // double Rl=0.6; // double c=2.8; // double n=1;

double d0=3; double Rh=4; double Rl=0.5; double c=3.8; double n=1;

void IMG_filter(plural *pImgDataPlural, int width, int height) {

 int cenW = width/2 + 0;  int cenH = height/2 + 0;

 for (int i=0; i<height; i++)  {   for (int j=0; j<width; j++)   {    //double d = sqrt( (i-height>>1)*(i-height>>1) + (j-width>>1)*(j-width>>1) *1.0);    double d = (i-cenH)*(i-cenH) + (j-cenW)*(j-cenW);    //double d = (j-cenW)*(j-cenW);    d = sqrt(d);    //    if (d<0.000001)    //     d = 0.000001;    //pCoff[i*width + j] = (Rh-Rl)* pow(exp(-c*(d0/d)),n)+Rl;    //pCoff[i*width + j] = (Rh-Rl)* exp(-c*(d0/d))+Rl;

   double tmp = (Rh-Rl)* exp(-c*(d0/d))+Rl;

   pImgDataPlural[i*width + j].real = pImgDataPlural[i*width + j].real * tmp;    pImgDataPlural[i*width + j].img  = pImgDataPlural[i*width + j].img  * tmp;

  }  }

}

void img_filter(double *pImgData, int w, int h) {  static plural *pImgDataPlural = NULL;  static plural *pTmpDataPlural = NULL;

 static plural *pFftWidth = NULL;  static plural *pFftHeight = NULL;

 static int PluralWidth  = 0;  static int PluralHeight = 0;

 // 1. 申请复数空间  // 计算  int m,n;  int k = 0;  int f = w - 1;

 for(k=1;(f=f>>1)!=0;k++) ;  m = 1<<k;

 f = h - 1;  for(k=1;(f=f>>1)!=0;k++) ;  n = 1<<k;

 if ( (m != PluralWidth) || (n != PluralHeight) )  {   if (pImgDataPlural != NULL) delete pImgDataPlural;   if (pTmpDataPlural != NULL) delete pTmpDataPlural;   if (pFftWidth != NULL) delete pFftWidth;   if (pFftHeight != NULL) delete pFftHeight;

  pImgDataPlural = new plural[m*n];   pTmpDataPlural = new plural[m*n];

  PluralWidth = m;   PluralHeight = n;

  FFT_Init(m, pFftWidth);   FFT_Init(n, pFftHeight);

 }

 image_print(pImgData, w, h);

 // 2. 填充数据

 plural tmp;

 for (int i=0; i<PluralHeight; i++)  {   for (int j=0; j<PluralWidth; j++)   {    tmp.img = 0.0;    if ( (i<h) && (j< w) )    {     tmp.real = pImgData[i*w + j];    }    else    {     tmp.real = 0;    }

   //频域坐标位移    if ((j / PluralWidth + j % PluralWidth) % 2 == 0)    {     //tmp.real = 0;    }    else    {     tmp.real = -tmp.real;    }

   if ((i / PluralHeight + i % PluralHeight) % 2 == 0)    {     //tmp.real = 0;    }    else    {     tmp.real = -tmp.real;    }

   pImgDataPlural[i*PluralWidth + j] = tmp;

  }  }

 plural_print(pImgDataPlural, PluralWidth,PluralHeight , "填充之后");  //printf(" ========================================== ");

 fft2(pImgDataPlural, PluralWidth,PluralHeight, pFftWidth, pFftHeight, pTmpDataPlural);

 IMG_filter(pImgDataPlural, PluralWidth,PluralHeight);

 ifft2(pImgDataPlural, PluralWidth,PluralHeight, pFftWidth, pFftHeight, pTmpDataPlural);

 // 反向填充数据

 for (int i=0; i<PluralHeight; i++)  {   for (int j=0; j<PluralWidth; j++)   {    tmp = pImgDataPlural[i*PluralWidth + j];

   //频域坐标位移    if ((j / PluralWidth + j % PluralWidth) % 2 == 0)    {     //tmp.real = 0;    }    else    {     tmp.real = -tmp.real;    }

   if ((i / PluralHeight + i % PluralHeight) % 2 == 0)    {     //tmp.real = 0;    }    else    {     tmp.real = -tmp.real;    }

   if ( (i<h) && (j< w) )    {     pImgData[i*w + j] = tmp.real;    }

  }  }

 plural_print(pImgDataPlural, PluralWidth,PluralHeight );

}

void RGB2YUV(BYTE *pImgDataRGB, int width, int height, double *pImgDataY, double *pImgDataU, double *pImgDataV) {  int pos;  double Y;  BYTE R,G,B;  for (int i=0; i<height; i++)  {   for (int j=0; j<width; j++)   {    pos = (i*width+j);    B = pImgDataRGB[pos * 3 + 0];    G = pImgDataRGB[pos * 3 + 1];    R = pImgDataRGB[pos * 3 + 2];

   //    Y =  0.299*R + 0.587*G + 0.114*B;    //    pImgDataY[pos] = Y;    //    pImgDataU[pos] = 0.492*(B-Y);//-0.147*R - 0.289*G + 0.436*B =    //    pImgDataV[pos] = 0.877*(R-Y);//0.615*R - 0.515*G - 0.100*B    //

   pImgDataY[pos] = R;

  }  }

}

void YUV2RGB(double *pImgDataY, double *pImgDataU, double *pImgDataV, int width, int height, BYTE *pImgDataRGB) {

 //  memset(pImgDataY,0,width*sizeof(double) * 30 * 3);

 int pos;  double Y,U,V;  double R,G,B;  for (int i=0; i<height; i++)  {   for (int j=0; j<width; j++)   {    pos = (i*width+j);

   Y = pImgDataY[pos];    //    U = pImgDataU[pos];    //    V = pImgDataV[pos];

   //    R = Y + 1.140*V;    //    G = Y - 0.394*U - 0.581*V;    //    B = Y + 2.032*U;    //    //    if (R > 255 ) R = 255;    //    if (G > 255 ) G = 255;    //    if (B > 255 ) B = 255;    //    //    if (R < 0 ) R = 0;    //    if (G < 0 ) G = 0;    //    if (B < 0 ) B = 0;

   if (Y > 255) Y = 255;    if (Y < 0) Y = 0;

   pImgDataRGB[pos * 3 + 0] = (BYTE)Y;    pImgDataRGB[pos * 3 + 1] = (BYTE)Y;    pImgDataRGB[pos * 3 + 2] = (BYTE)Y;

   //    pImgDataRGB[pos * 3 + 0] = (BYTE)B;    //    pImgDataRGB[pos * 3 + 1] = (BYTE)G;    //    pImgDataRGB[pos * 3 + 2] = (BYTE)R;

  }  }

}

IMGPRO_API void video_imgpro(BYTE *pImgData, int width, int height) {

 //  memset(pImgData + width*3*10 ,255, width*3*1);

 // 0 创建参数

 // 1. CSC  static double *pY = new double[width*height];  static double *pU = new double[width*height];  static double *pV = new double[width*height];

 static int oldWidth = 0;  static int oldHeight = 0;  if ( (oldWidth != width) || (oldHeight != height) )  {   if (pY) delete pY;   if (pU) delete pU;   if (pV) delete pV;

  pY = new double[width*height];   pU = new double[width*height];   pV = new double[width*height];

 }

 RGB2YUV(pImgData, width, height, pY, pU, pV);

 img_filter(pY,width,height);

 YUV2RGB( pY, pU, pV, width, height, pImgData);

}