图像处理 之 同态滤波
没什么好讲的 原图 --> 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 \n", 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("\n%s\n", 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("\n"); } }
void image_print(double *pImgData, int w, int h) { return; printf("\n"); for (int i=0; i<h; i++) { for (int j=0; j<w; j++) { printf("[ %7.2f]",pImgData[i*w + j]); } printf("\n"); } }
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("\n==========================================\n");
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);
}