001 #include <opencv2/opencv.hpp>
002 #include <opencv2/opencv_lib.hpp>
003 
004 void cvShiftDFT(CvArr *src_arr, CvArr *dst_arr);
005 
006 void myDftImage(IplImage *imageSrc, CvMat *matComp) {
007     CvMat *matReal, *matImag;
008     CvSize size;
009     int sizeX, sizeY;
010     int sx, sy;
011     double r, g, b;
012 
013     size = cvGetSize(imageSrc);
014     sizeX = size.width; sizeY = size.height;
015 
016     matReal = cvCreateMat(sizeY, sizeX, CV_64FC1);
017     matImag = cvCreateMat(sizeY, sizeX, CV_64FC1);
018 
019     for(sy = 0; sy < sizeY; sy++) {
020         for(sx = 0; sx < sizeX; sx++) {
021             b = ((uchar*)(imageSrc->imageData + imageSrc->widthStep * sy))[imageSrc->nChannels * sx + 0];
022             g = ((uchar*)(imageSrc->imageData + imageSrc->widthStep * sy))[imageSrc->nChannels * sx + 1];
023             r = ((uchar*)(imageSrc->imageData + imageSrc->widthStep * sy))[imageSrc->nChannels * sx + 2];
024             (matReal->data.db + matReal->height * sy)[sx] = (0.299 * r + 0.587 * g + 0.114 * b) / 255.0;
025         }
026     }
027 
028     cvSetZero(matImag);
029     cvMerge(matReal, matImag, NULL, NULL, matComp);
030 
031     cvDFT(matComp, matComp, CV_DXT_FORWARD);
032 
033     cvReleaseMat(&matReal);
034     cvReleaseMat(&matImag);
035 
036     return;
037 }
038 
039 void myNormalizeMat(CvMat *matCompSrc, CvMat *matCompDst) {
040     CvMat *matReal, *matImag;
041     CvMat *matAmp;
042 
043     matReal = cvCreateMat(matCompSrc->height, matCompSrc->width, CV_64FC1);
044     matImag = cvCreateMat(matCompSrc->height, matCompSrc->width, CV_64FC1);
045     matAmp = cvCreateMat(matCompSrc->height, matCompSrc->width, CV_64FC1);
046 
047     cvSplit(matCompSrc, matReal, matImag, NULL, NULL);
048 
049     cvPow(matReal, matReal, 2.0); cvPow(matImag, matImag, 2.0);
050     cvAdd(matReal, matImag, matAmp, NULL);
051     cvPow(matAmp, matAmp, 0.5);
052 
053     cvSplit(matCompSrc, matReal, matImag, NULL, NULL);
054 
055     cvDiv(matReal, matAmp, matReal, 1.0);
056     cvDiv(matImag, matAmp, matImag, 1.0);
057 
058     cvMerge(matReal, matImag, NULL, NULL, matCompDst);
059 
060     cvReleaseMat(&matReal);
061     cvReleaseMat(&matImag);
062     cvReleaseMat(&matAmp);
063 
064     return;
065 }
066 
067 int main(int argc, char *argv[]) {
068     IplImage *imageOrg0, *imageOrg1;
069     IplImage *imagePha;
070     CvMat *matComp0, *matComp1;
071     CvMat *matReal, *matImag;
072     CvSize size;
073     int sizeX, sizeY;
074 
075     if(argc != 3) {
076         printf("Usage: %s <ImageFile0> <ImageFile1>\n", argv[0]);
077         exit(-1);
078     }
079 
080     imageOrg0 = cvLoadImage(argv[1], CV_LOAD_IMAGE_GRAYSCALE);
081     imageOrg1 = cvLoadImage(argv[2], CV_LOAD_IMAGE_GRAYSCALE);
082     if(imageOrg0 == NULL || imageOrg1 == NULL) {
083         fprintf(stderr, "Error: Can not open ImageFile(s)\n");
084         exit(-1);
085     }
086 
087     size = cvGetSize(imageOrg0);
088     sizeX = size.width; sizeY = size.height;
089     imagePha = cvCreateImage(size, IPL_DEPTH_64F, 1);
090 
091     matComp0 = cvCreateMat(sizeY, sizeX, CV_64FC2);
092     matComp1 = cvCreateMat(sizeY, sizeX, CV_64FC2);
093     matReal = cvCreateMat(sizeY, sizeX, CV_64FC1);
094     matImag = cvCreateMat(sizeY, sizeX, CV_64FC1);
095 
096     // (1) 2つの原画像をフーリエ変換し、空間周波数スペクトル（複素平面）を得る
097     myDftImage(imageOrg0, matComp0);
098     myDftImage(imageOrg1, matComp1);
099 
100     // (2) 2つのスペクトルの合成積を取る（2つ目の複素平面は複素共役とする）
101     cvMulSpectrums(matComp0, matComp1, matComp0, CV_DXT_MUL_CONJ);
102 
103     // (3) 合成積を振幅で正規化する
104     myNormalizeMat(matComp0, matComp0);
105 
106     // (4) 合成積を逆離散フーリエ変換し、結果を実数配列と虚数配列に分解する
107     cvDFT(matComp0, matComp0, CV_DXT_INVERSE);
108     cvSplit(matComp0, matReal, matImag, NULL, NULL);
109 
110     // (5) 第1象限と第3象限、第2象限と第4象限を入れ替える
111     cvShiftDFT(matReal, matReal);
112 
113     // (6) 実数配列を正規化し、位相限定相関画像を得る
114     cvNormalize(matReal, imagePha, 0.0, 1.0, CV_MINMAX, NULL);
115 
116     cvNamedWindow("Original Image0", CV_WINDOW_AUTOSIZE);
117     cvNamedWindow("Original Image1", CV_WINDOW_AUTOSIZE);
118     cvNamedWindow("Phase Image", CV_WINDOW_AUTOSIZE);
119 
120     cvShowImage("Original Image0", imageOrg0);
121     cvShowImage("Original Image1", imageOrg1);
122     cvShowImage("Phase Image", imagePha);
123 
124     cvWaitKey(0);
125 
126     cvDestroyWindow("Original Image0");
127     cvDestroyWindow("Original Image1");
128     cvDestroyWindow("Phase Image");
129     cvReleaseImage(&imageOrg0);
130     cvReleaseImage(&imageOrg1);
131     cvReleaseImage(&imagePha);
132     cvReleaseMat(&matReal);
133     cvReleaseMat(&matImag);
134 
135     exit(0);
136 }
137 
138 void cvShiftDFT(CvArr *src_arr, CvArr *dst_arr) {
139     CvMat *tmp = 0;
140     CvMat q1stub, q2stub;
141     CvMat q3stub, q4stub;
142     CvMat d1stub, d2stub;
143     CvMat d3stub, d4stub;
144     CvMat *q1, *q2, *q3, *q4;
145     CvMat *d1, *d2, *d3, *d4;
146 
147     CvSize size = cvGetSize(src_arr);
148     CvSize dst_size = cvGetSize(dst_arr);
149     int cx, cy;
150 
151     if(dst_size.width != size.width || dst_size.height != size.height) {
152         cvError(CV_StsUnmatchedSizes, "cvShiftDFT", "Source and Destination arrays must have equal sizes", __FILE__, __LINE__);
153     }
154 
155     if(src_arr == dst_arr) {
156         tmp = cvCreateMat(size.height / 2, size.width / 2, cvGetElemType(src_arr));
157     }
158     cx = size.width / 2;
159     cy = size.height / 2;
160 
161     q1 = cvGetSubRect(src_arr, &q1stub, cvRect(0, 0, cx, cy));
162     q2 = cvGetSubRect(src_arr, &q2stub, cvRect(cx, 0, cx, cy));
163     q3 = cvGetSubRect(src_arr, &q3stub, cvRect(cx, cy, cx, cy));
164     q4 = cvGetSubRect(src_arr, &q4stub, cvRect(0, cy, cx, cy));
165     d1 = cvGetSubRect(src_arr, &d1stub, cvRect(0, 0, cx, cy));
166     d2 = cvGetSubRect(src_arr, &d2stub, cvRect(cx, 0, cx, cy));
167     d3 = cvGetSubRect(src_arr, &d3stub, cvRect(cx, cy, cx, cy));
168     d4 = cvGetSubRect(src_arr, &d4stub, cvRect(0, cy, cx, cy));
169 
170     if(src_arr != dst_arr) {
171         if(!CV_ARE_TYPES_EQ(q1, d1)) {
172             cvError(CV_StsUnmatchedFormats, "cvShiftDFT", "Source and Destination arrays must have the same format", __FILE__, __LINE__);
173         }
174         cvCopy(q3, d1, 0);
175         cvCopy(q4, d2, 0);
176         cvCopy(q1, d3, 0);
177         cvCopy(q2, d4, 0);
178     } else {
179         cvCopy(q3, tmp, 0);
180         cvCopy(q1, q3, 0);
181         cvCopy(tmp, q1, 0);
182         cvCopy(q4, tmp, 0);
183         cvCopy(q2, q4, 0);
184         cvCopy(tmp, q2, 0);
185     }
186 
187     return;
188 }