wzpan
4/2/2014 - 1:21 PM

Switch two images' phase.

Switch two images' phase.

#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>

using namespace cv;

// Rearrange the quadrants of a Fourier image so that the origin is at
// the image center
void shiftDFT(Mat &fImage )
{
    Mat tmp, q0, q1, q2, q3;

    // first crop the image, if it has an odd number of rows or columns

    fImage = fImage(Rect(0, 0, fImage.cols & -2, fImage.rows & -2));

    int cx = fImage.cols / 2;
    int cy = fImage.rows / 2;

    // rearrange the quadrants of Fourier image
    // so that the origin is at the image center

    q0 = fImage(Rect(0, 0, cx, cy));
    q1 = fImage(Rect(cx, 0, cx, cy));
    q2 = fImage(Rect(0, cy, cx, cy));
    q3 = fImage(Rect(cx, cy, cx, cy));

    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);

    q1.copyTo(tmp);
    q2.copyTo(q1);
    tmp.copyTo(q2);
}

int main()
{
    // Load an image
    Mat I1 = imread("lena.jpg", CV_LOAD_IMAGE_GRAYSCALE);
    Mat I2 = imread("pepper.jpg", CV_LOAD_IMAGE_GRAYSCALE);

    Mat fI1;
    Mat fI2;
    I1.convertTo(fI1, CV_32F);
    I2.convertTo(fI2, CV_32F);

    //expand input image to optimal size
    int m = getOptimalDFTSize( I1.rows );
    int n = getOptimalDFTSize( I1.cols ); 

    Mat padded1, padded2;
    
    // on the border add zero values
    copyMakeBorder(fI1, padded1, 0, m - I1.rows, 0, n - I1.cols, BORDER_CONSTANT, Scalar::all(0));
    copyMakeBorder(fI2, padded2, 0, m - I2.rows, 0, n - I2.cols, BORDER_CONSTANT, Scalar::all(0));
    
    //Perform DFT
    Mat fourierTransform1;
    Mat fourierTransform2;
    Mat planes1[2], planes2[2];

    dft(fI1, fourierTransform1, DFT_SCALE|DFT_COMPLEX_OUTPUT);
    dft(fI2, fourierTransform2, DFT_SCALE|DFT_COMPLEX_OUTPUT);

    shiftDFT(fourierTransform1);
    shiftDFT(fourierTransform2);

    // compute the magnitude and phase
    split(fourierTransform1, planes1);// planes[0] = Re(DFT(I)), planes[1] = Im(DFT(I))
    split(fourierTransform2, planes2);// planes[0] = Re(DFT(I)), planes[1] = Im(DFT(I))
    
    Mat ph1, mag1;
    mag1.zeros(planes1[0].rows, planes1[0].cols, CV_32F);
    ph1.zeros(planes1[0].rows, planes1[0].cols, CV_32F);
    cartToPolar(planes1[0], planes1[1], mag1, ph1);

    Mat ph2, mag2;
    mag2.zeros(planes2[0].rows, planes2[0].cols, CV_32F);
    ph2.zeros(planes2[0].rows, planes2[0].cols, CV_32F);
    cartToPolar(planes2[0], planes2[1], mag2, ph2);

    polarToCart(mag1, ph2, planes1[0], planes1[1]);
    polarToCart(mag2, ph1, planes2[0], planes2[1]);

    merge(planes1, 2, fourierTransform1);
    merge(planes2, 2, fourierTransform2);

    shiftDFT(fourierTransform1);
    shiftDFT(fourierTransform2);
    
    //Perform  IDFT
    Mat inverseTransform1, inverseTransform2;
    
    dft(fourierTransform1, inverseTransform1, DFT_INVERSE|DFT_REAL_OUTPUT);
    dft(fourierTransform2, inverseTransform2, DFT_INVERSE|DFT_REAL_OUTPUT);
        
    namedWindow("original image 1");
    imshow("original image 1", I1);
    namedWindow("original image 2");
    imshow("original image 2", I2);
    waitKey(0);

    cv::Mat out1, out2;
    inverseTransform1.convertTo(out1, CV_8U);
    inverseTransform2.convertTo(out2, CV_8U);
    
    namedWindow("result image 1");
    imshow("result image 1", out1);
    namedWindow("result image 2");
    imshow("result image 2", out2);
    waitKey(0); 
}