C#图像快速傅立叶变换及二维傅立叶变化

来源:互联网 发布:井下电钳网络考试平台 编辑:程序博客网 时间:2024/05/20 16:00
有关傅立叶变化资料:

http://www.academictutorials.com/graphics/graphics-2d-fourier-transform.asp

http://blog.csdn.net/lnstree/article/details/7412837

了解完傅立叶变化,如果想了解具体的算法实现的话,自行查找相关资料。

傅立叶变化对于图像处理的意义可以简单的描述成,一道阳光通过三棱镜分解成不同波段的光波。而傅立叶变换的功能就是三棱镜。

许多在空间域无法较好解决的问题,转换成频率域就可以很简单的解决。

 

以下是网上FFT(快速傅立叶变换的VC代码)(基于2的幂数)

 

/************************************************************************* * * 函数名称: *   FFT() * * 参数: *   complex<double> * TD- 指向时域数组的指针 *   complex<double> * FD- 指向频域数组的指针 *   r-2的幂数,即迭代次数 * * 返回值: *   无。 * * 说明: *   该函数用来实现快速付立叶变换。 * ************************************************************************/VOID WINAPI FFT(complex<double> * TD, complex<double> * FD, int r){// 付立叶变换点数LONGcount;// 循环变量inti,j,k;// 中间变量intbfsize,p;// 角度doubleangle;complex<double> *W,*X1,*X2,*X;// 计算付立叶变换点数count = 1 << r;// 分配运算所需存储器W  = new complex<double>[count / 2];X1 = new complex<double>[count];X2 = new complex<double>[count];// 计算加权系数for(i = 0; i < count / 2; i++){angle = -i * PI * 2 / count;W[i] = complex<double> (cos(angle), sin(angle));}// 将时域点写入X1memcpy(X1, TD, sizeof(complex<double>) * count);// 采用蝶形算法进行快速付立叶变换for(k = 0; k < r; k++){for(j = 0; j < 1 << k; j++){bfsize = 1 << (r-k);for(i = 0; i < bfsize / 2; i++){p = j * bfsize;X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2]) * W[i * (1<<k)];}}X  = X1;X1 = X2;X2 = X;}// 重新排序for(j = 0; j < count; j++){p = 0;for(i = 0; i < r; i++){if (j&(1<<i)){p+=1<<(r-i-1);}}FD[j]=X1[p];     }// 释放内存delete W;delete X1;delete X2;}

 

 

 

以下是相应的二维傅立叶变换源码

/************************************************************************* * * 函数名称: *   Fourier() * * 参数: *   LPSTR lpDIBBits    - 指向源DIB图像指针 *   LONG  lWidth       - 源图像宽度(象素数) *   LONG  lHeight      - 源图像高度(象素数) * * 返回值: *   BOOL               - 成功返回TRUE,否则返回FALSE。 * * 说明: *   该函数用来对图像进行付立叶变换。 * ************************************************************************/BOOL WINAPI Fourier(LPSTR lpDIBBits, LONG lWidth, LONG lHeight){// 指向源图像的指针unsigned char*lpSrc;// 中间变量doubledTemp;// 循环变量LONGi;LONGj;// 进行付立叶变换的宽度和高度(2的整数次方)LONGw;LONGh;intwp;inthp;// 图像每行的字节数LONGlLineBytes;// 计算图像每行的字节数lLineBytes = WIDTHBYTES(lWidth * 8);// 赋初值w = 1;h = 1;wp = 0;hp = 0;// 计算进行付立叶变换的宽度和高度(2的整数次方)while(w * 2 <= lWidth){w *= 2;wp++;}while(h * 2 <= lHeight){h *= 2;hp++;}// 分配内存complex<double> *TD = new complex<double>[w * h];complex<double> *FD = new complex<double>[w * h];// 行for(i = 0; i < h; i++){// 列for(j = 0; j < w; j++){// 指向DIB第i行,第j个象素的指针lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;// 给时域赋值TD[j + w * i] = complex<double>(*(lpSrc), 0);}}for(i = 0; i < h; i++){// 对y方向进行快速付立叶变换FFT(&TD[w * i], &FD[w * i], wp);}// 保存变换结果for(i = 0; i < h; i++){for(j = 0; j < w; j++){TD[i + h * j] = FD[j + w * i];}}/*for(i = 0; i < w; i++){// 对x方向进行快速付立叶变换FFT(&TD[i * h], &FD[i * h], hp);}*/// 行for(i = 0; i < h; i++){// 列for(j = 0; j < w; j++){// 计算频谱dTemp = sqrt(FD[j * h + i].real() * FD[j * h + i].real() +          FD[j * h + i].imag() * FD[j * h + i].imag()) / 100;// 判断是否超过255if (dTemp > 255){// 对于超过的,直接设置为255dTemp = 255;}// 指向DIB第(i<h/2 ? i+h/2 : i-h/2)行,第(j<w/2 ? j+w/2 : j-w/2)个象素的指针// 此处不直接取i和j,是为了将变换后的原点移到中心//lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - (i<h/2 ? i+h/2 : i-h/2)) + (j<w/2 ? j+w/2 : j-w/2);// 更新源图像* (lpSrc) = (BYTE)(dTemp);}}// 删除临时变量delete TD;delete FD;// 返回return TRUE;}


转成相应C#源码:

       //  快速傅立叶变换       public static void FFT(COMPLEX[] t,  COMPLEX[] f, int r)   // t为时域,f为频域 r为2的幂数       {           long count;           int i, j, k, p, bsize;           COMPLEX[] W;           COMPLEX[] X1;           COMPLEX[] X2;           COMPLEX[] X;           COMPLEX comp;           double angle;  // 计算加权时所需角度           count = 1 << r;           W = new COMPLEX[count / 2];           X1 = new COMPLEX[count];           X2 = new COMPLEX[count];           X = new COMPLEX[count];           for (i = 0; i < count / 2; i++)           {               angle = i * Math.PI * 2 / count;               W[i].re = (double)Math.Cos(angle);               W[i].im = -(double)Math.Sin(angle);           }           t.CopyTo(X1, 0);           for (k = 0; k < r; k++)           {               for (j = 0; j < 1 << k; j++)               {                   bsize = 1 << (r - k);                   for (i = 0; i < bsize/2; i++)                   {                       p = j * bsize;                       X2[i + p].im = X1[i + p].im + X1[i + p + bsize / 2].im;                       X2[i + p].re = X1[i + p].re + X1[i + p + bsize / 2].re;                       comp.im = X1[i + p].im - X1[i + p + bsize / 2].im;                       comp.re = X1[i + p].re - X1[i + p + bsize / 2].re;                       X2[i + p + bsize / 2].re = comp.re * W[i * (1 << k)].re -                           comp.im * W[i * (1 << k)].im;                       X2[i + p + bsize / 2].im = comp.re * W[i * (1 << k)].im +                        comp.im * W[i * (1 << k)].re; ;                   }               }               X = X1;               X1 = X2;               X2 = X;           }           for (j = 0; j < count; j++)           {               p = 0;               for (i = 0; i < r; i++)               {                   if ((j & (1 << i)) != 0)                   {                       p += 1 << (r - i - 1);                   }               }               f[j].re = X1[p].re;               f[j].im = X1[p].im;           }          }
 // 二维傅立叶变换       public static Bitmap Fourier(Bitmap tp)       {           // 原图像的宽与高           int w = tp.Width;           int h = tp.Height;           // 傅立叶变换的实际宽高           long lw = 1;           long lh = 1;           // 迭代次数           int wp = 0; int hp = 0;           long i, j;           long n,m;           double temp;           byte [] ky = new byte[w*h];           ky = ChangeByte(tp);           COMPLEX[] t; COMPLEX[] f;           while (lw * 2 <= w)           {               lw *= 2;               wp++;           }           while (lh * 2 <= h)           {               lh *= 2;               hp++;           }           t = new COMPLEX[lw * lh];           f = new COMPLEX[lw * lh];           COMPLEX[] tw = new COMPLEX[lw];           COMPLEX[] th = new COMPLEX[lw];           for (i = 0; i < lh; i++)           {               for (j = 0; j < lw; j++)               {                   t[i*lw + j].re = ky[ i*w+ j];                   t[i *lw + j].im = 0;               }           }           for (i = 0; i < lh; i++) // 垂直方向傅立叶变换           {                     Array.Copy(t, i * lw, tw, 0, lw);             Array.Copy(f, i * lw, th, 0, lw);              FFT(tw, th, wp);           // Array.Copy(tw, 0, t, i * lw, lw);              Array.Copy(th, 0, f, i * lw, lw);           }           for (i = 0; i < lh; i++)           {               for (j = 0; j < lw; j++)               {                   t[j * lh + i].re = f[i * lw + j].re;                   t[j * lh + i].im = f[i * lw + j].im;               }           }           COMPLEX[] ow = new COMPLEX[lh];           COMPLEX[] oh = new COMPLEX[lh];           for (i = 0; i < lw; i++)           {               Array.Copy(t, i * lh, ow, 0, lh);               Array.Copy(f, i * lh, oh, 0, lh);               FFT(ow, oh, hp);               //Array.Copy(ow, 0, t, i * lh, lh);               oh.CopyTo(f, i * lh);           }           for (i = 0; i < lh; i++)           {               for (j = 0; j < lw; j++)               {                   temp = Math.Sqrt(f[j * lh + i].re * f[j * lh + i].re + f[j * lh + i].im * f[j * lh + i].im) / 100;                   if (temp > 255)                   {                       temp = 255;                   }                   n = i < lh / 2 ? i + lh / 2 : i - lh / 2;                   m = j < lw / 2 ? j + lw / 2 : j - lw / 2;                   ky[n* w + m] = (byte)(temp);               }           }           tp = ChangeBitmap(ky, tp);           return tp;       }


 

测试图片:       

原图傅立叶变换图ok了。

                                                               
	
				
		
原创粉丝点击