opencv傅里叶变换实例

OpenCV快速傅里叶变换实例

int main(){    Mat I = imread("ted_cruz.jpg", CV_LOAD_IMAGE_GRAYSCALE);    if (I.empty())        return -1;    cout << I.size() << endl;    Mat padded;                            //expand input image to optimal size    int m = getOptimalDFTSize(I.rows);    int n = getOptimalDFTSize(I.cols); // on the border add zero values    copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0));    cout << padded.size() << endl;    Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };    Mat complexI;    merge(planes, 2, complexI);         // Add to the expanded another plane with zeros    dft(complexI, complexI);            // this way the result may fit in the source matrix    // compute the magnitude and switch to logarithmic scale    // => log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))    split(complexI, planes);                   // planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))    magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude    Mat magI = planes[0];    magI += Scalar::all(1);                    // switch to logarithmic scale    log(magI, magI);    // crop the spectrum, if it has an odd number of rows or columns    magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));    // rearrange the quadrants of Fourier image  so that the origin is at the image center    int cx = magI.cols / 2;    int cy = magI.rows / 2;    Mat q0(magI, Rect(0, 0, cx, cy));   // Top-Left - Create a ROI per quadrant    Mat q1(magI, Rect(cx, 0, cx, cy));  // Top-Right    Mat q2(magI, Rect(0, cy, cx, cy));  // Bottom-Left    Mat q3(magI, Rect(cx, cy, cx, cy)); // Bottom-Right    Mat tmp;                           // swap quadrants (Top-Left with Bottom-Right)    q0.copyTo(tmp);    q3.copyTo(q0);    tmp.copyTo(q3);    q1.copyTo(tmp);                    // swap quadrant (Top-Right with Bottom-Left)    q2.copyTo(q1);    tmp.copyTo(q2);    normalize(magI, magI, 0, 1, CV_MINMAX); // Transform the matrix with float values into a    // viewable image form (float between values 0 and 1).    imshow("Input Image", I);    // Show the result    imshow("spectrum magnitude", magI);    waitKey();    return 0;}OpenCV 中 傅里叶变换 FFT，代码如下：
[cpp] view plain copy
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
离散傅里叶变换就是将图像从空间域转换到频域，这一转换基本原理为：
任一函数都可以表示成无数个正弦和余弦函数的和的形式，二维图像的傅里叶变换可用公式表示为：
其中，f是空间域，F是频域，转换之后的频域值是复数，因此显示傅里叶变换之后的结果需要使用实物图像加虚数图像或者幅度图像加相位图像的形式。
示例;
[cpp] view plain copy
#include"stdafx.h"  
#include <opencv2/core/utility.hpp>  
#include "opencv2/imgproc.hpp"  
#include "opencv2/imgcodecs.hpp"  
#include "opencv2/highgui.hpp"  
#include"opencv2/core/core.hpp"  
#include <iostream>  
  
using namespace cv;  
using namespace std;  
  
void DFT_Change(Mat &, string);  
// 定义全局变量  
  
  
int main() {  
    system("color 5E");  
    //读取图片  
    Mat image1 = imread("E:\\pictures\\For_Project\\Opencv_Test\\哈尔施特塔\\05.jpg", 0);  
    namedWindow("【原始图像1】", 1);  
    imshow("【原始图像1】", image1);  
  
    Mat image2 = imread("E:\\pictures\\For_Project\\Opencv_Test\\哈尔施特塔\\b4740bd7f5ee276803fdaa95bfe128d6.jpg", 0);  
    namedWindow("【原始图像2】", 1);  
    imshow("【原始图像2】", image2);  
  
    Mat image3 = imread("E:\\pictures\\For_Project\\Opencv_Test\\哈尔施特塔\\cbf4d74061b713bb5da0bdafedbf1178.jpg", 0);  
    namedWindow("【原始图像3】", 1);  
    imshow("【原始图像3】", image3);  
      
    double time1 = static_cast<double>(getTickCount());  
    //进行DFT傅里叶变换操作   
    string windows_name1 = "频谱幅值1";  
    DFT_Change(image1,windows_name1);  
    //计算运行时间并输出    
    time1 = ((double)getTickCount() - time1) / getTickFrequency();  
    cout << "1.该方法运行时间为：" << time1 << "秒" << endl;  
  
    double time2 = static_cast<double>(getTickCount());  
    //进行DFT傅里叶变换操作   
    string windows_name2 = "频谱幅值2";  
    DFT_Change(image2,windows_name2);  
    //计算运行时间并输出    
    time2 = ((double)getTickCount() - time2) / getTickFrequency();  
    cout << "2.该方法运行时间为：" << time2 << "秒" << endl;  
  
    double time3 = static_cast<double>(getTickCount());  
    //进行DFT傅里叶变换操作   
    string windows_name3 = "频谱幅值3";  
    DFT_Change(image3,windows_name3);  
    //计算运行时间并输出    
    time3 = ((double)getTickCount() - time3) / getTickFrequency();  
    cout << "3.该方法运行时间为：" << time3 << "秒" << endl;  
  
    //等待键盘按键‘q’退出  
    while (char(waitKey(1)) != 'q') {}  
    return 0;  
}  
  
void DFT_Change(Mat &image,string windows_name) {  
      
    Mat srcImage = image.clone();  
    // ShowHelpText();  
  
    // 将输入图像延扩到最佳的尺寸，边界用0补充q  
    int m = getOptimalDFTSize(srcImage.rows);  
    int n = getOptimalDFTSize(srcImage.cols);  
  
    //将添加的像素值初始化为0  
    Mat padded;  
    copyMakeBorder(srcImage, padded, 0, m - srcImage.rows, 0, n - srcImage.cols, BORDER_CONSTANT, Scalar::all(0));  
  
    //为傅里叶变换的结果（实部和虚部）分配存储空间  
    //讲planes数组组合成一个多通道的数组complexI  
    Mat planes[] = { Mat_<float>(padded),Mat::zeros(padded.size(),CV_32F) };  
    Mat complexI;  
    merge(planes, 2, complexI);  
  
    //就地进行离散傅里叶变换  
    dft(complexI, complexI);  
  
    //将复数转换为幅值  
    //讲多通道的数组complexI分离成几个单通道数组  
    split(complexI, planes);  
    magnitude(planes[0], planes[1], planes[0]);  
    Mat magnitudeImage = planes[0];  
  
    //进行对数尺度缩放  
    magnitudeImage += Scalar::all(1);  
    log(magnitudeImage, magnitudeImage);  
  
    //剪切和重分布幅度图像象限  
    //若有奇数或奇数行，进行频谱裁剪  
    magnitudeImage = magnitudeImage(Rect(0, 0, magnitudeImage.cols&-2, magnitudeImage.rows&-2));  
  
    //重新排列傅里叶图像中的象限，使得原点在图像中心  
    int cx = magnitudeImage.cols / 2;  
    int cy = magnitudeImage.rows / 2;  
    //ROI区域的左上  
    Mat q0(magnitudeImage, Rect(0, 0, cx, cy));  
    //ROI区域的右上  
    Mat q1(magnitudeImage, Rect(cx, 0, cx, cy));  
    //ROI区域的左下  
    Mat q2(magnitudeImage, Rect(0, cy, cx, cy));  
    //ROI区域的右下  
    Mat q3(magnitudeImage, Rect(cx, cy, cx, cy));  
  
    //交换象限，左上与右下交换，右上和左下交换  
    Mat tmp1, tmp2;  
    q0.copyTo(tmp1);  
    q3.copyTo(q0);  
    tmp1.copyTo(q3);  
  
    q1.copyTo(tmp2);  
    q2.copyTo(q1);  
    tmp2.copyTo(q2);  
  
    //归一化处理，用0-1之间的浮点值将矩阵变换为可视的图像格式  
    normalize(magnitudeImage, magnitudeImage, 0, 1, NORM_MINMAX);  
  
    imshow(windows_name, magnitudeImage);  
}  


http://www.cnblogs.com/tornadomeet/archive/2012/07/26/2610414.html