001 #include <opencv2/opencv.hpp>
002 #include <opencv2/opencv_lib.hpp>
003 
004 int main(int argc, char* argv[]) {
005     IplImage *imageOrg, *imagePha;
006     CvMat *matReal, *matImag, *matComp;
007     CvMat *matRealNorm, *matImagNorm;
008     CvMat *matAmp;
009     CvSize size;
010     int sizeX, sizeY;
011 
012     if(argc != 2) {
013         printf("Usage: %s <ImageFile>\n", argv[0]);
014         exit(-1);
015     }
016 
017     imageOrg = cvLoadImage(argv[1], CV_LOAD_IMAGE_GRAYSCALE);
018     if(imageOrg == NULL) {
019         fprintf(stderr, "Error: Can not open ImageFile [%s]\n", argv[1]);
020         exit(-1);
021     }
022 
023     size = cvGetSize(imageOrg);
024     sizeX = size.width; sizeY = size.height;
025     imagePha = cvCreateImage(size, IPL_DEPTH_64F, 1);
026 
027     matReal = cvCreateMat(sizeY, sizeX, CV_64FC1);
028     matImag = cvCreateMat(sizeY, sizeX, CV_64FC1);
029     matComp = cvCreateMat(sizeY, sizeX, CV_64FC2);
030 
031     matRealNorm = cvCreateMat(sizeY, sizeX, CV_64FC1);
032     matImagNorm = cvCreateMat(sizeY, sizeX, CV_64FC1);
033 
034     matAmp = cvCreateMat(sizeY, sizeX, CV_64FC1);
035 
036     // (1) 入力画像を実数配列にコピーし虚数配列とマージして、複素平面を作成
037     cvScale(imageOrg, matReal);
038     cvSetZero(matImag);
039     cvMerge(matReal, matImag, NULL, NULL, matComp);
040 
041     // (2) 離散フーリエ変換を行い、結果を実数配列と虚数配列に分解
042     cvDFT(matComp, matComp, CV_DXT_FORWARD);
043     cvSplit(matComp, matReal, matImag, NULL, NULL);
044     cvSplit(matComp, matRealNorm, matImagNorm, NULL, NULL); // 正規化用
045 
046     // (3) スペクトルの振幅を計算 Amp = sqrt(Real^2 + Imag^2)
047     cvPow(matReal, matReal, 2.0); cvPow(matImag, matImag, 2.0);
048     cvAdd(matReal, matImag, matAmp, NULL);
049     cvPow(matAmp, matAmp, 0.5);
050 
051     // (4) 実数配列と虚数配列を、振幅で正規化
052     cvDiv(matRealNorm, matAmp, matRealNorm, 1.0);
053     cvDiv(matImagNorm, matAmp, matImagNorm, 1.0);
054 
055     // (5) 正規化した実数配列と虚数配列をマージして、複素平面を作成
056     cvMerge(matRealNorm, matImagNorm, NULL, NULL, matComp);
057 
058     // (6) 逆離散フーリエ変換を行い、結果を実数配列と虚数配列に分解
059     cvDFT(matComp, matComp, CV_DXT_INVERSE);
060     cvSplit(matComp, matReal, matImag, NULL, NULL);
061 
062     // (7) 実数配列を正規化し、位相画像を作成
063     cvNormalize(matReal, imagePha, 0.0, 1.0, CV_MINMAX, NULL);
064 
065 
066     cvNamedWindow("Original Image", CV_WINDOW_AUTOSIZE);
067     cvNamedWindow("Phase Image", CV_WINDOW_AUTOSIZE);
068 
069     cvShowImage("Original Image", imageOrg);
070     cvShowImage("Phase Image", imagePha);
071 
072     cvWaitKey(0);
073 
074     cvDestroyWindow("Original Image");
075     cvDestroyWindow("Phase Image");
076     cvReleaseImage(&imageOrg);
077     cvReleaseImage(&imagePha);
078     cvReleaseMat(&matReal);
079     cvReleaseMat(&matImag);
080     cvReleaseMat(&matComp);
081     cvReleaseMat(&matRealNorm);
082     cvReleaseMat(&matImagNorm);
083     cvReleaseMat(&matAmp);
084 
085     exit(0);
086 }