fft.cs

2维fft用于数字图像处理 C# dll装配件

using System;
using System.Drawing;


namespace MyApp.FFT
{
	/// <summary>
	/// Summary description for Class1.
	/// </summary>
	/// 	
	public class complexd
	{
		public double real;
		public double im;
		public complexd(double x,double y)
		{
			real=x;
			im=y;
		}
	}
	public class MyFFT
	{
		private int myHeight;
		private int myWidth;
		private double[][] BitIll;

		public MyFFT()
		{
			InitialBitIll();
		}
		private void InitialBitIll()
		{
			BitIll=new double[5][];
			BitIll[0]=new double[5];
			BitIll[0][0]=0;
			BitIll[0][1]=1;
			BitIll[0][2]=2;
			BitIll[0][3]=1;
			BitIll[0][4]=0;
			BitIll[1]=new double[5];
			BitIll[1][0]=1;
			BitIll[1][1]=4;
			BitIll[1][2]=8;
			BitIll[1][3]=4;
			BitIll[1][4]=1;
			BitIll[2]=new double[5];
			BitIll[2][0]=2;
			BitIll[2][1]=8;
			BitIll[2][2]=16;
			BitIll[2][3]=8;
			BitIll[2][4]=2;
			BitIll[3]=new double[5];
			BitIll[3][0]=1;
			BitIll[3][1]=4;
			BitIll[3][2]=8;
			BitIll[3][3]=4;
			BitIll[3][4]=1;
			BitIll[4]=new double[5];
			BitIll[4][0]=0;
			BitIll[4][1]=1;
			BitIll[4][2]=2;
			BitIll[4][3]=1;
			BitIll[4][4]=0;
			for(int i=0;i<5;i++)
				for(int j=0;j<5;j++)
				BitIll[i][j]=BitIll[i][j]/80;		
		}
		private complexd[][] GetBits(Bitmap myBitmap)
		{
			int height=myBitmap.Height;
			int width=myBitmap.Width;
			int nHeight=(int) Math.Pow(2,Math.Ceiling(Math.Log(height,2)));
			int nWidth=(int) Math.Pow(2,Math.Ceiling(Math.Log(width,2)));
			myHeight=height;
			myWidth=width;
			complexd[][] array=new complexd[nHeight][];
			int i,j;
			for(i=0;i<=nHeight-1;i++)
			{
				array[i]=new complexd[nWidth];
				for(j=0;j<=nWidth-1;j++)
				{
					if(i<height&&j<width)
						array[i][j]=new complexd(myBitmap.GetPixel(j,i).R,0);
					else array[i][j]=new complexd(0,0);
				}
			}
			return array;
		}
		private double[][] GetArray(Bitmap myBitmap)
		{
			int height=myBitmap.Height;
			int width=myBitmap.Width;

			double[][] array=new double[height][];
			int i,j;
			for(i=0;i<=height-1;i++)
			{
				array[i]=new double[width];
				for(j=0;j<=width-1;j++)
				{
					array[i][j]=myBitmap.GetPixel(j,i).R;
				}
			}
			return array;
		}
		private complexd Mux(complexd x1,complexd x2)
		{		
			double real=x1.real*x2.real-x1.im*x2.im;
			double im=x1.real*x2.im+x1.im*x2.real;
			complexd x=new complexd(real,im);
			return x;
		}
		private complexd Add(complexd x1,complexd x2)
		{
			complexd x=new complexd(0,0);
			x.real=x1.real+x2.real;
			x.im=x1.im+x2.im;
			return x;
		}
		private complexd Del(complexd x1,complexd x2)
		{
			complexd x=new complexd(0,0);
			x.real=x1.real-x2.real;
			x.im=x1.im-x2.im;
			return x;
		}
		private void BitReverse(int[] br, int m,int N)
		{
			for(int i=0;i<N;i++)
			{
				br[i]=i;
			}
			int lh=N/2;
			int j=lh;
			int n1=N-2;

			for(int i=1;i<=n1;i++)
			{
				if(i<j)
				{
					int t=br[i];
					br[i]=br[j];
					br[j]=t;
				}
				int k=lh;
				while(j>=k)
				{
					j=j-k;
					k/=2;
				}
				j=j+k;
			}
		}
		private complexd WNR(int N,int r)
		{
			double x=Math.Cos(((-2)*r*Math.PI)/N);
			double y=Math.Sin(((-2)*r*Math.PI)/N);
			complexd c=new complexd(x,y);
			return c;
		}
		private complexd IWNR(int N,int r)
		{
			double x=Math.Cos(((-2)*r*Math.PI)/N);
			double y=Math.Sin(((-2)*r*Math.PI)/N);
			complexd c=new complexd(x,y);
			return c;
		}
		private complexd[][] Rotate(complexd[][] array)
		{
			int n=array[0].Length;
			complexd[][] ar=new complexd[n][];
			for(int i=0;i<n;i++)
			{
				ar[i]=new complexd[array.Length];
				for(int j=0;j<array.Length;j++)
					ar[i][j]=array[j][i];
			}
			return ar;
		}
		private void Operate(complexd[] array,int i,int j,complexd W)
		{
			complexd ii=array[i];
			complexd jj=array[j];
			array[i]=Add(ii,Mux(W,jj));
			array[j]=Del(ii,Mux(W,jj));
		}
		private void FFT1(complexd[] array,int sig)
		{
			int Width=array.Length;
			int bit=(int)Math.Ceiling(Math.Log(Width,2))+1;
			int[] br=new int[Width];
			BitReverse(br,bit,Width);
			complexd[] ar=(complexd[])array.Clone();
			for(int i=0;i<Width;i++)
				array[i]=ar[br[i]];
			int m=(int) Math.Log(Width,2);
			
			for(int l=1;l<=m;l++)
			{
				int b=(int)Math.Pow(2,l-1);
				for(int j=0;j<=b-1;j++)
				{
					int r=(int)Math.Pow(2,m-l)*j;
					for (int k=j;k<=Width-1;k+=(int)Math.Pow(2,l))
					{
						complexd w;
						if(sig==0)
							w=WNR(Width,r);
						else w=IWNR(Width,r);
						Operate(array,k,k+b,w);
					}
				}
			}
			if(sig==1)
			{
				for(int i=0;i<Width;i++)
				{
					array[i].real=array[i].real/Width;
					array[i].im=array[i].im/Width;
				}
			}
		}
		private complexd[][] FFT2(complexd[][] PicData,int sig)
		{
			int height=PicData.Length;
			int width=PicData[0].Length;
			for(int i=0;i<height;i++)
				FFT1(PicData[i],sig);	
			PicData=Rotate(PicData);
			for(int i=0;i<width;i++)
				FFT1(PicData[i],sig);	
			PicData=Rotate(PicData);
			return PicData;
		}
		private complexd[][] FFT2(Bitmap myBitmap)
		{
			complexd[][] PicData=GetBits(myBitmap);
			PicData=FFT2(PicData,0);
			PicData=Array2Array(PicData);
			return PicData;
		}
		public Bitmap Butterworth(Bitmap myBitmap,double pass)
		{
			complexd[][] PicData=FFT2(myBitmap);
			Bitmap m=new Bitmap(PicData.Length,PicData[0].Length);
			//Bitmap m=new Bitmap(myHeight,myWidth);
			double tp=TotalPower(PicData);
			int D0=0;
			int M=m.Height;
			int N=m.Width;
			while(true)
			{
				if(PartPower(PicData,D0)/tp>=pass)
					break;
				else D0++;
			}
			complexd[][] f=new complexd[M][];

			for(int i=0;i<M;i++)
			{
				f[i]=new complexd[N];
				for(int j=0;j<N;j++)
				{
					double bt=Math.Sqrt(Math.Pow(Math.Abs(i-M/2),2)+Math.Pow(Math.Abs(j-N/2),2))/D0;
					bt=1/(1+Math.Pow(bt,2));
					f[i][j]=new complexd(PicData[i][j].real*bt,PicData[i][j].im*bt);
				}
			}
			f=Array2Array(f);
			f=FFT2(f,1);
			for(int i=0;i<M;i++)
				for(int j=0;j<N;j++)
				{
					int c=(int)r2(f[i][j]);
					if(c>255) c=255;
					if(c<0) c=0;
					m.SetPixel(j,i,Color.FromArgb(c,c,c));
				}
			m.RotateFlip(RotateFlipType.Rotate180FlipNone);
			Rectangle rect=new Rectangle(0,0,myWidth,myHeight);
			m=m.Clone(rect,m.PixelFormat);
			return m;
		}

		public Bitmap IdealLowPass(Bitmap myBitmap,double pass)
		{
			complexd[][] PicData=FFT2(myBitmap);
			Bitmap m=new Bitmap(PicData.Length,PicData[0].Length);
			double tp=TotalPower(PicData);
			int D0=0;
			int M=PicData.Length;
			int N=PicData[0].Length;
			while(true)
			{
				if(PartPower(PicData,D0)/tp>=pass)
					break;
				else D0++;
			}
			complexd[][] f=new complexd[M][];

			for(int i=0;i<M;i++)
			{
				f[i]=new complexd[N];
				for(int j=0;j<N;j++)
				{
					f[i][j]=new complexd(0,0);
					if((Math.Pow(Math.Abs(i-M/2),2)+Math.Pow(Math.Abs(j-N/2),2))<=D0*D0)
					{
						f[i][j]=PicData[i][j];
					}
				}
			}
			f=Array2Array(f);
			f=FFT2(f,1);
			for(int i=0;i<m.Height;i++)
				for(int j=0;j<m.Width;j++)
				{
					int c=(int)r2(f[i][j]);
					if(c>255) c=255;
					if(c<0) c=0;
					m.SetPixel(j,i,Color.FromArgb(c,c,c));
				}
			m.RotateFlip(RotateFlipType.Rotate180FlipNone);
			Rectangle rect=new Rectangle(0,0,myWidth,myHeight);
			m=m.Clone(rect,m.PixelFormat);
			return m;
		}
		private double TotalPower(complexd[][] array)
		{
			double re=0;
			for(int i=0;i<array.Length;i++)
				for(int j=0;j<array[0].Length;j++)
				{
					double f=r2(array[i][j]);
					re+=f*f;
				}
			return re;
		}
		private double PartPower(complexd[][] array,int D0)
		{
			double re=0;
			int M=array.Length;
			int N=array[0].Length;
			for(int i=0;i<array.Length;i++)
				for(int j=0;j<array[0].Length;j++)
				{
					if((Math.Pow(Math.Abs(i-M/2),2)+Math.Pow(Math.Abs(j-N/2),2))<=D0*D0)
					{
						double f=r2(array[i][j]);
						re+=f*f;
					}
				}
			return re;
		}
		public complexd[][] Array2Array(complexd[][] array)
		{
			
			int M=array[0].Length;
			int N=array.Length;
			complexd[][] re=new complexd[N][];
			for(int i=0;i<N/2;i++)
			{
				re[i]=new complexd[M];
				for(int j=0;j<M/2;j++)//1
				{
					re[i][j]=array[i+N/2][j+M/2];
				}
				for(int j=M/2;j<M;j++)//2
				{
					re[i][j]=array[i+N/2][j-M/2];
				}
			}
			for(int i=N/2;i<N;i++)//4
			{
				re[i]=new complexd[M];
				for(int j=0;j<M/2;j++)
				{
					re[i][j]=array[i-N/2][j+M/2];
				}
				for(int j=M/2;j<M;j++)
				{
					re[i][j]=array[i-N/2][j-M/2];
				}
			}
			return re;
		}
		
		private double r2(complexd c)
		{
			return Math.Sqrt(c.real*c.real+c.im*c.im);
		}
		public Bitmap GetFFTPic(Bitmap myBitmap)
		{
			complexd[][] PicData=GetBits(myBitmap);
			PicData=FFT2(PicData,0);
			complexd[][] PicData1=Array2Array(PicData);
			double[][] data=GetPowerData(PicData1);
			double[] maxmin=GetMaxMinData(data);
			maxmin[0]=(int)Math.Log(1+maxmin[0],Math.E);
			maxmin[1]=(int)Math.Log(1+maxmin[1],Math.E);
			Bitmap m=new Bitmap(PicData1.Length,PicData1[0].Length);
			for(int i=0;i<m.Height;i++)
				for(int j=0;j<m.Width;j++)
				{
					int c=(255/(int)(maxmin[0]-maxmin[1]))*((int)Math.Log(1+data[i][j],Math.E));
					if(c>255) c=255;
					if(c<0) c=0;
					m.SetPixel(j,i,Color.FromArgb(c,c,c));
				}
			return m;
		}
		private double[][] GetPowerData(complexd[][] array)
		{
			double[][] re=new double[array.Length][];
			for(int i=0;i<array.Length;i++)
			{
				re[i]=new double[array[0].Length];
				for(int j=0;j<array[0].Length;j++)
				{
					re[i][j]=r2(array[i][j]);
				}
			}
			return re;
		}
		private double[] GetMaxMinData(double[][] array)
		{
			double[] re=new double[2];
			re[0]=array[0][0];
			re[1]=array[0][0];
			for(int i=0;i<array.Length;i++)
				for(int j=0;j<array[0].Length;j++)
				{
					if(array[i][j]>re[0]) re[0]=array[i][j];
					else re[1]=array[i][j];
				}
			return re;		
		}
		public double GetError(Bitmap myBitmap1,Bitmap myBitmap2)
		{
			double[][] array1=GetArray(myBitmap1);
			double[][] array2=GetArray(myBitmap2);
			double re=0;
			for(int i=0;i<=array1.Length-1;i++)
				for(int j=0;j<=array1[0].Length-1;j++)
				{
					re+=Math.Pow((array1[i][j]-array2[i][j]),2);
				}
			re=re/(array1.Length*array1[0].Length);
			re=Math.Sqrt(re);
			return re;
		}
		
		private double GetError(double[][] array1,double[][] array2)
		{
			double re=0;
			for(int i=0;i<=array1.Length-1;i++)
				for(int j=0;j<=array1[0].Length-1;j++)
				{
					re+=Math.Pow((array1[i][j]-array2[i][j]),2);
				}
			re=re/(array1.Length*array1[0].Length);
			re=Math.Sqrt(re);
			return re;
		}
		


	}
	
}