【OpenCV学习】DFT变换

作者：gnuhpc
出处：http://www.cnblogs.com/gnuhpc/

#include "cv.h"
#include "highgui.h"
#include "cxcore.h"
void cvShiftDFT(CvArr *src_arr,CvArr *dst_arr)
{
    CvMat * tmp;
    CvMat q1stub,q2stub;
    CvMat q3stub,q4stub;
    CvMat d1stub,d2stub;
    CvMat d3stub,d4stub;
    CvMat *q1,*q2,*q3,*q4;
    CvMat *d1,*d2,*d3,*d4;
    CvSize size = cvGetSize(src_arr);
    CvSize dst_size = cvGetSize(dst_arr);
    int cx,cy;
    if ((dst_size.width!= size.width)||(dst_size.height!=size.height))
    {
        cvError(CV_StsUnmatchedSizes,"cvShiftDFT",
            "Source and Destination arrays must have the same sizes",
            __FILE__,__LINE__);
    }
    if (src_arr == dst_arr)
    {
        tmp = cvCreateMat(size.height/2,size.width/2,cvGetElemType(src_arr));
    }
    cx=size.width/2;//取出图像的原点
    cy=size.height/2;
    q1=cvGetSubRect(src_arr,&q1stub,cvRect(0,0,cx,cy));
    //取出图像的第一象限，由q1指针指向它
    q2=cvGetSubRect(src_arr,&q2stub,cvRect(cx,0,cx,cy));
    //取出图像的第二象限，由q2指针指向它
    q3=cvGetSubRect(src_arr,&q3stub,cvRect(cx,cy,cx,cy));
    //取出图像的第三象限，由q3指针指向它
    q4=cvGetSubRect(src_arr,&q4stub,cvRect(0,cy,cx,cy));
    //取出图像的第四象限，由q4指针指向它
   
    d1=cvGetSubRect(src_arr,&d1stub,cvRect(0,0,cx,cy));
    d2=cvGetSubRect(src_arr,&d2stub,cvRect(cx,0,cx,cy));
    d3=cvGetSubRect(src_arr,&d3stub,cvRect(cy,cy,cx,cy));
    d4=cvGetSubRect(src_arr,&d4stub,cvRect(0,cy,cx,cy));
    if (src_arr!=dst_arr)
    {
        if (!CV_ARE_TYPES_EQ(q1,d1))
        {
            cvError(CV_StsUnmatchedSizes,"cvShiftDFT",
                "Source and Destination arrays must have the same sizes",
                __FILE__,__LINE__);       
        }
        //以图像中心为原点，调整傅里叶变换图像的四个象限区，
        //即第一与第三象限交换，第二与第四象限交换
        cvCopy(q3,d1,0);
        cvCopy(q4,d2,0);
        cvCopy(q1,d3,0);
        cvCopy(q2,d4,0);
    }
    else
    {//若源矩阵和目的矩阵相同则直接在源矩阵中进行操作
        cvCopy(q3,tmp,0);
        cvCopy(q1,q3,0);
        cvCopy(tmp,q1,0);
        cvCopy(q4,tmp,0);
        cvCopy(q2,q4,0);
        cvCopy(tmp,q2,0);
    }
}
int main(int argc,char ** argv)
{
    const char* filename =(argc>=2?argv[1]:"lena.jpg");
    IplImage *im;
    IplImage *realInput,*imaginaryInput,*complexInput;
    IplImage *image_Re,*image_Im;
    int dft_M,dft_N;
    CvMat *dft_A;
    CvMat tmp;
    double m,M;
    im = cvLoadImage(filename, CV_LOAD_IMAGE_GRAYSCALE  );//加载图像
    if (!im)
    {
        return -1;
    }
    //分配空间
    realInput = cvCreateImage(cvGetSize(im),IPL_DEPTH_64F,1);//单通道
    imaginaryInput =cvCreateImage(cvGetSize(im),IPL_DEPTH_64F,1);//单通道
    complexInput = cvCreateImage(cvGetSize(im),IPL_DEPTH_64F,2);//双通道
    cvScale(im,realInput,1.0,0.0);
    //#define cvScale cvConvertScale=>readInput=im
    cvZero(imaginaryInput);
    //清空这个图像的内容
    cvMerge(realInput,imaginaryInput,NULL,NULL,complexInput);
    //混合这两个图像作为complexInput的两个通道
    /*得到最优DFT尺寸 */
    dft_M = cvGetOptimalDFTSize(im->height-1);
    dft_N = cvGetOptimalDFTSize(im->width-1);
   
    dft_A = cvCreateMat(dft_M,dft_N,CV_64FC2);
    image_Re = cvCreateImage(cvSize(dft_N,dft_M),IPL_DEPTH_64F,1);//实部
    image_Im = cvCreateImage(cvSize(dft_N,dft_M),IPL_DEPTH_64F,1);//虚部
    cvGetSubRect(dft_A,&tmp,cvRect(0,0,im->width,im->height));
    cvCopy(complexInput,&tmp,NULL);
    if( dft_A->cols > im->width )//若得到的最优DFT尺寸在宽度上大于原图，则重新取
    {
        cvGetSubRect(dft_A,&tmp,cvRect(im->width,0,dft_A->cols-im->width,im->height));
        cvZero(&tmp);
    }
    cvDFT(dft_A,dft_A,CV_DXT_FORWARD,complexInput->height);
   
    cvNamedWindow("win",0);
    cvNamedWindow("magnitude",0);
    cvShowImage("win",im);
   
    //分割出实部和虚部
    cvSplit(dft_A,image_Re,image_Im,0,0);
   
    //计算功率谱 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,NULL);//image_Re<=image_Re+image_Im
    cvPow(image_Re,image_Re,0.5);
    //计算log(1+Mag)
    cvAddS(image_Re,cvScalarAll(1.0),image_Re,NULL);
    cvLog(image_Re,image_Re);
    cvShiftDFT(image_Re,image_Re);
    cvMinMaxLoc(image_Re,&m,&M,NULL,NULL,NULL);
    cvScale(image_Re,image_Re,1.0/(M-m),1.0*(-m)/(M-m));
    cvShowImage("magnitude",image_Re);
    cvWaitKey(-1);
    return 0;
} 

 

作者：gnuhpc
出处：http://www.cnblogs.com/gnuhpc/