1. 程式人生 > >opencv 中 快速傅立葉變換 FFT

opencv 中 快速傅立葉變換 FFT

opencv 中 傅立葉變換 FFT,程式碼如下:

void fft2(IplImage *src, IplImage *dst)
{   //實部、虛部
	IplImage *image_Re = 0, *image_Im = 0, *Fourier = 0;
	//   int i, j;
	image_Re = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);  //實部
	//Imaginary part
	image_Im = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);  //虛部
	//2 channels (image_Re, image_Im)
	Fourier = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 2);
	// Real part conversion from u8 to 64f (double)
	cvConvertScale(src, image_Re);
	// Imaginary part (zeros)
	cvZero(image_Im);
	// Join real and imaginary parts and stock them in Fourier image
	cvMerge(image_Re, image_Im, 0, 0, Fourier);

	// Application of the forward Fourier transform
	cvDFT(Fourier, dst, CV_DXT_FORWARD);
	cvReleaseImage(&image_Re);
	cvReleaseImage(&image_Im);
	cvReleaseImage(&Fourier);
}

void fft2shift(IplImage *src, IplImage *dst)
{
	IplImage *image_Re = 0, *image_Im = 0;
	int nRow, nCol, i, j, cy, cx;
	double scale, shift, tmp13, tmp24;
	image_Re = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);
	//Imaginary part
	image_Im = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);
	cvSplit( src, image_Re, image_Im, 0, 0 );
	//具體原理見岡薩雷斯數字影象處理p123
	// Compute the magnitude of the spectrum Mag = sqrt(Re^2 + Im^2)
	//計算傅立葉譜
	cvPow( image_Re, image_Re, 2.0);
	cvPow( image_Im, image_Im, 2.0);
	cvAdd( image_Re, image_Im, image_Re);
	cvPow( image_Re, image_Re, 0.5 );
	//對數變換以增強灰度級細節(這種變換使以窄帶低灰度輸入影象值對映
	//一寬頻輸出值,具體可見岡薩雷斯數字影象處理p62)
	// Compute log(1 + Mag);
	cvAddS( image_Re, cvScalar(1.0), image_Re ); // 1 + Mag
	cvLog( image_Re, image_Re ); // log(1 + Mag)

	//Rearrange the quadrants of Fourier image so that the origin is at the image center
	nRow = src->height;
	nCol = src->width;
	cy = nRow/2; // image center
	cx = nCol/2;
	//CV_IMAGE_ELEM為OpenCV定義的巨集,用來讀取影象的畫素值,這一部分就是進行中心變換
	for( j = 0; j < cy; j++ ){
		for( i = 0; i < cx; i++ ){
			//中心化,將整體份成四塊進行對角交換
			tmp13 = CV_IMAGE_ELEM( image_Re, double, j, i);
			CV_IMAGE_ELEM( image_Re, double, j, i) = CV_IMAGE_ELEM(
				image_Re, double, j+cy, i+cx);
			CV_IMAGE_ELEM( image_Re, double, j+cy, i+cx) = tmp13;

			tmp24 = CV_IMAGE_ELEM( image_Re, double, j, i+cx);
			CV_IMAGE_ELEM( image_Re, double, j, i+cx) =
				CV_IMAGE_ELEM( image_Re, double, j+cy, i);
			CV_IMAGE_ELEM( image_Re, double, j+cy, i) = tmp24;
		}
	}
	//歸一化處理將矩陣的元素值歸一為[0,255]
	//[(f(x,y)-minVal)/(maxVal-minVal)]*255
	double minVal = 0, maxVal = 0;
	// Localize minimum and maximum values
	cvMinMaxLoc( image_Re, &minVal, &maxVal );
	// Normalize image (0 - 255) to be observed as an u8 image
	scale = 255/(maxVal - minVal);
	shift = -minVal * scale;
	cvConvertScale(image_Re, dst, scale, shift);
	cvReleaseImage(&image_Re);
	cvReleaseImage(&image_Im);
}

void CCVMFCView::OnFuliyeTransform()
{
	IplImage *src;
	IplImage *Fourier;   //傅立葉係數
	IplImage *dst ;
	IplImage *ImageRe;
	IplImage *ImageIm;
	IplImage *Image;
	IplImage *ImageDst;
	double m,M;
	double scale;
	double shift;
	//src = workImg;
	if(workImg->nChannels==3)
		OnColorToGray();
	src=cvCreateImage(cvGetSize(workImg),IPL_DEPTH_64F,workImg->nChannels);  //源影象
	imageClone(workImg,&src);
	cvFlip(src);

	Fourier = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,2);
	dst = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,2);
	ImageRe = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);
	ImageIm = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);
	Image = cvCreateImage(cvGetSize(src),src->depth,src->nChannels);
	ImageDst = cvCreateImage(cvGetSize(src),src->depth,src->nChannels);
	fft2(src,Fourier);                  //傅立葉變換
	fft2shift(Fourier, Image);          //中心化
	cvDFT(Fourier,dst,CV_DXT_INV_SCALE);//實現傅立葉逆變換,並對結果進行縮放
	cvSplit(dst,ImageRe,ImageIm,0,0);

	cvNamedWindow("源影象",0);
	cvShowImage("源影象",src);             
	//對陣列每個元素平方並存儲在第二個引數中
	cvPow(ImageRe,ImageRe,2);               
	cvPow(ImageIm,ImageIm,2);
	cvAdd(ImageRe,ImageIm,ImageRe,NULL);
	cvPow(ImageRe,ImageRe,0.5);
	cvMinMaxLoc(ImageRe,&m,&M,NULL,NULL);
	scale = 255/(M - m);
	shift = -m * scale;
	//將shift加在ImageRe各元素按比例縮放的結果上,儲存為ImageDst
	cvConvertScale(ImageRe,ImageDst,scale,shift);

	cvNamedWindow("傅立葉譜",0);
	cvShowImage("傅立葉譜",Image);
	cvNamedWindow("傅立葉逆變換",0);
	cvShowImage("傅立葉逆變換",ImageDst);
	//釋放影象
	cvWaitKey(10000);
	cvReleaseImage(&src);
	cvReleaseImage(&Image);
	cvReleaseImage(&ImageIm);
	cvReleaseImage(&ImageRe);
	cvReleaseImage(&Fourier);
	cvReleaseImage(&dst);
	cvReleaseImage(&ImageDst);
	Invalidate();
}


from: http://blog.csdn.net/abcjennifer/article/details/7359952