matlab fft2 opencv,基于OPENCV的FFT2变换

// opencv.cpp : 定义控制台应用程序的入口点。

//

#include "stdafx.h"

#include

using namespace std;

IplImage * DoPyrdown(IplImage* old);

void GuassMohu(IplImage *img);

IplImage * doCanny(IplImage* in,double lowThresh,double

highThresh,double aperture);

void fft2(IplImage *src, IplImage *dst);

void fft2shift(IplImage *src, IplImage *dst);

int _tmain(int argc, _TCHAR* argv[])

{

IplImage

*procimg=cvLoadImage("D:\\ImgProc\\lena.bmp",0);

IplImage * Fourier =

cvCreateImage(cvGetSize(procimg),IPL_DEPTH_64F,2);

//CvRect roi=cvRect(10,10,50,50);

//cvAddS(sub_img,cvScalar(100),sub_img);

//cvReleaseImageHeader(&sub_img);

//cvLine(procimg,cvPoint(0,0),cvPoint(100,100),cvScalar(80,80,80),10);

double m,M;

double scale;

double shift;

cvNamedWindow("old");

cvShowImage("old",procimg);

fft2(procimg,Fourier);

IplImage * ImageRe;IplImage * ImageIm;IplImage *

Image;IplImage * ImageDst;

ImageRe =

cvCreateImage(cvGetSize(procimg),IPL_DEPTH_64F,1);

ImageIm =

cvCreateImage(cvGetSize(procimg),IPL_DEPTH_64F,1);

Image =

cvCreateImage(cvGetSize(procimg),procimg->depth,procimg->nChannels);

ImageDst =

cvCreateImage(cvGetSize(procimg),procimg->depth,procimg->nChannels);

fft2shift(Fourier, Image);

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);

cvWaitKey(0);

while(1);

return 0;

}

void GuassMohu(IplImage *img)

{

cvNamedWindow("In");

cvNamedWindow("Out");

cvShowImage("In",img);

IplImage * outimg=cvCreateImage(cvGetSize(img),8,3);

cvSmooth(img,outimg,CV_GAUSSIAN,5,5);

cvShowImage("Out",outimg);

cvReleaseImage(&outimg);

cvWaitKey(0);

cvDestroyWindow("In");

cvDestroyWindow("Out");

}

IplImage * DoPyrdown(IplImage* old)

{

assert(old->width%2==0&&old->height%2==0);

IplImage *

happynew=cvCreateImage(cvSize(old->width/2,old->height/2),old->depth,old->nChannels);

cvPyrDown(old,happynew);

cvNamedWindow("In");

cvNamedWindow("Out");

cvShowImage("In",old);

cvShowImage("Out",happynew);

cvWaitKey(0);

return(happynew);

}

IplImage * doCanny(IplImage* in,double lowThresh,double

highThresh,double aperture)

{

if (in->nChannels!=3)

return(0);

IplImage

*out=cvCreateImage(cvGetSize(in),IPL_DEPTH_8U,1);

cvCanny(in,out,lowThresh,highThresh,aperture);

cvNamedWindow("In");

cvNamedWindow("Out");

cvShowImage("In",in);

cvShowImage("Out",out);

cvWaitKey(0);

return(out);

}

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, 1, 0);

// 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);

}