Template Matching Issue

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

Template Matching Issue

alinmallen
Hi:
I am useing Fourier-Mellin method to match 2 same size images for invariance to rotation and following bellow steps for the matching:

1.FFT transformation to SrcImg and TmpImg and calculate their Log-Magnitude and save as SrcImg(2) and TmpImg(2).

2.Log-Polar Transformation to SrcImg(2) and TmpImg(2) and correspondingly saved as SrcImg(3) and TmpImg(3).

3.FFT transformation to SrcImg(3) and TmpImg(3),and calculate their phase correlation peak value location to determine rotation.

The problem is I always get peak value at cooridinate (0,0) on phase correlation plane.
I have my sample image and result show in bellow album.
http://groups.yahoo.com/group/OpenCV/photos/album/759048714/pic/list?mode=tn&order=ordinal&start=1&dir=asc

Can someone help me to tell anything I did wrong on that job?

Appreciate!!

Reply | Threaded
Open this post in threaded view
|

Image matching through rotation, scale using log polar frequency domain

Heather Larson-2
I was just able to successfully correlate rotated images using the an fft'ed log-polar image.  I will include my code below.  I only care about rotation at this point, which shows up in my y axis of the inverse fft'ed phase correlation, but finding the scale isn't hard, its on the x axis, and finding translation would mean going back out of log-polar and getting the phase correlation of the fft'ed images, which wouldn't be hard given this code.

I found many papers where they took the fft of the image and then moved it to log-polar space, but my implementation does the opposite and works really well.  It makes sense to me the way I did it.  Good luck filtering through all my code.

I have this code in two files, a header and c++ file, so hopefully I am including everything that is needed from those two files (and not anything more because I don't want to cause any more confusion)

From the Header file
// Correlation of scene and model maps using DFT in OpenCV libraries
        const double CORRELATION_THRESHOLD = 0.1;
        CvMat *rotationMatrix; // A commonly used rotation matrix
        CvSize imgSize;
        CvSize imgSmallSize;
        CvSize dftSize;
        IplImage *sceneSmallImg;
        IplImage *modelSmallImg;
        IplImage *sceneImg;
        IplImage *modelImg;
        IplImage *scenePolar;
        IplImage *modelPolar;
        IplImage *rotationImg;
        int nDFTHeight;
        int nDFTWidth;
        IplImage *realInput;
        IplImage *imaginaryInput;
        IplImage *complexInput;
        CvMat *sceneDFT;
        CvMat *modelDFT;
        IplImage *imageRe;
        IplImage *imageIm;
        IplImage *imageImMag;
        IplImage *imageMag;


// Image correlation data structures
        rotationMatrix = cvCreateMat(2,3,CV_32F); imgSmallSize = cvSize(CELL_X, CELL_Y);
        imgSize = cvSize(CELL_X, CELL_Y);
        sceneSmallImg = cvCreateImage( imgSmallSize, nBinSize, 1 );
        modelSmallImg = cvCreateImage( imgSmallSize, nBinSize, 1 );
        sceneImg = cvCreateImage( imgSize, nBinSize, 1 );
        modelImg = cvCreateImage( imgSize, nBinSize, 1 );
        scenePolar = cvCreateImage( imgSize, nBinSize, 1 );
        modelPolar = cvCreateImage( imgSize, nBinSize, 1 );
        rotationImg = cvCreateImage( imgSize, nBinSize, 1 );
        // Optimal image size for fastest DFT computation
        nDFTHeight = cvGetOptimalDFTSize( imgSize.height );
        nDFTWidth = cvGetOptimalDFTSize( imgSize.width );
        // Allocate floating point frames used for DFT (real, imaginary, and complex)
        realInput = cvCreateImage( imgSize, IPL_DEPTH_64F, 1 );
        imaginaryInput = cvCreateImage( imgSize, IPL_DEPTH_64F, 1 );
        complexInput = cvCreateImage( imgSize, IPL_DEPTH_64F, 2 );
        // Allocate DFT output images
        dftSize = cvSize(nDFTWidth, nDFTHeight);
        sceneDFT = cvCreateMat( nDFTHeight, nDFTWidth, CV_64FC2 );
        modelDFT = cvCreateMat( nDFTHeight, nDFTWidth, CV_64FC2 );
        imageRe = cvCreateImage( dftSize, IPL_DEPTH_64F, 1 );
        imageIm = cvCreateImage( dftSize, IPL_DEPTH_64F, 1 );
        imageImMag = cvCreateImage( dftSize, IPL_DEPTH_64F, 1 );
        imageMag = cvCreateImage( dftSize, IPL_DEPTH_64F, 1 );
 

// Find correlation between scene and model images (obstacle bitmaps) and return a correcting angle (in radians)
// This function is invariant to translation, rotation, or scaling.  This algorithm is based on the Fourier shift theorem
// which states that translation is invariant in the frequency domain when comparing the phase correlation.  If the images
// are in log-polar space, then rotation and scale are represented as translation and are invariant in the frequency domain.
//
// To understand the theory behind this function, please see the following papers:
// An FFT-Based Technique for Translation, Rotation, and Scale-Invariant Image Registration by B. Srinivasa Reddy and B.N.Chatterji, 1996
// Robust Image Registration Using Log-Polar Transform by G. Wolberg, S. Zokai
// RST-Invariant Digital Image Watermarking Based on Log-Polar Mapping and Phase Correlation by Dong Zheng, Jiying Zhao, and Abdulmotaleb El Saddik
// Search for Fourier-Mellin Transform
//
// If one would rather not go into the frequency domain, then you can also do cross correlation in OpenCV using the
// cvMatchTemplate function.
double WorldModel::imageCorrelation(obstacle_bit_map_incoming *map, double dCurrentHeading)
{
        double dAngleCorrection = 0.0;

        // Before we update the world model with the current sensor view object scene, we should make sure they match up correctly
        // Up to this point, we have the input of some gps or compass, telling us our current heading.  This, however, cannot always
        // be trusted, so we will try to correlate the current sensor view object scene with the world model to find the best fit.
        // We do this by converting the world model and sensor view to log-polar images and correlating the phase

        // Rotate the world model image to the current heading
        double dAngle = dCurrentHeading * 180.0 / M_PI;
        dAngle = 90.0 - dAngle;
        cv2DRotationMatrix(cvPoint2D32f((double)nRobotPixelX, (double)nRobotPixelY), dAngle, 1, rotationMatrix);
        cvWarpAffine(cvWorldModel, tempImg, rotationMatrix, CV_INTER_NN | CV_WARP_FILL_OUTLIERS, cvScalarAll(0));

        // put the data in sensorMap into the sceneImg
        cvZero(sceneImg);
        cvZero(sceneSmallImg);
        BwImage img(sceneSmallImg);
        for (int i=0; i<CELL_Y; i++)
                for (int j=0; j<CELL_X; j++)
                        img[i][j] = ((map->obstacles[i][j/8] & (0x01 << (j%8))) >> (j%8)) * 255; // map->obstacles is a bitmap, with 8 pixel values in each byte
        cvResize(sceneSmallImg, sceneImg); // If we want to take a look at the data in high res, increase the size of this img
       
        cvZero(modelImg);
        cvZero(modelSmallImg);
        // Resize the world model image to the same size as the scene
        CvRect roi = cvRect(nRobotPixelX-CELL_X/2, nRobotPixelY-CELL_Y, CELL_X, CELL_Y);
        cvSetImageROI(tempImg, roi);
        cvResize(tempImg, modelSmallImg);//, CV_INTER_NN);
        cvResetImageROI(tempImg);
        cvResize(modelSmallImg, modelImg);
       
        // Change the values inside the model to be either 255 or 0
        img = BwImage(sceneImg);
        BwImage img2 = BwImage(modelImg);
        for (int i=0; i<modelImg->height; i++)
        {
                for (int j=0; j<modelImg->width; j++)
                {
                        img[i][j] = img[i][j]? 255: 0;
                        img2[i][j] = img2[i][j]? 255: 0;
                }
        }
       
        // The log polar transform in OpenCV
        CvPoint2D32f pntCenter = cvPoint2D32f(imgSize.width/2, imgSize.height);
        double M = 20.0; // Not really sure what the best value of M should be, but this works well
        cvZero( scenePolar );
        cvZero( modelPolar );
        cvLogPolar( sceneImg, scenePolar, pntCenter, M, CV_INTER_LINEAR+CV_WARP_FILL_OUTLIERS);
        cvLogPolar( modelImg, modelPolar, pntCenter, M, CV_INTER_LINEAR+CV_WARP_FILL_OUTLIERS);

        CvMat tmp;
        // Processing of scene
        cvScale( scenePolar, realInput, 1.0, 0 );
        cvZero( imaginaryInput );
        cvMerge( realInput, imaginaryInput, NULL, NULL, complexInput );
        cvGetSubRect( sceneDFT, &tmp, cvRect(0, 0, scenePolar->width, scenePolar->height) );
        cvCopy( complexInput, &tmp, NULL );
        if (sceneDFT->cols > scenePolar->width)
        {
                cvGetSubRect( sceneDFT, &tmp, cvRect(scenePolar->width, 0, sceneDFT->cols - scenePolar->width, scenePolar->height) );
                cvZero( &tmp );
        }
        cvDFT( sceneDFT, sceneDFT, CV_DXT_FORWARD, complexInput->height );
        cvSplit( sceneDFT, imageRe, imageIm, 0, 0 );

        // Processing of model
        cvScale( modelPolar, realInput, 1.0, 0 );
        cvMerge( realInput, imaginaryInput, NULL, NULL, complexInput );
        cvGetSubRect( modelDFT, &tmp, cvRect(0, 0, modelPolar->width, modelPolar->height) );
        cvCopy( complexInput, &tmp, NULL );
        if (sceneDFT->cols > modelPolar->width)
        {
                cvGetSubRect( modelDFT, &tmp, cvRect(modelPolar->width, 0, sceneDFT->cols - modelPolar->width, modelPolar->height) );
                cvZero( &tmp );
        }
        cvDFT( modelDFT, modelDFT, CV_DXT_FORWARD, complexInput->height );

        // At this point, we want to get the phase correlation of these two images
        // If f and g are two images, then F and G are their fourier transforms
        // The cross-power spectrum (phase correlation) is FG*/|FG*| where * is the conjugate of the image
        // There are a few ways to calculate this:
        //
        // 1.
        // We can do this by comparing the difference in phase of each complex image using the following equation
        // pdm (phase difference matrix) = exp(i*(atan2(imag(F), real(F)) - atan2(imag(G), real(G))));
        // pcf (phase correlation function) = real(ifft(pdm));
        //
        // 2.
        // Multiply one image by the conjugate of the other
        // fftr = F .* conj(G); // Matlab code
        // Normalize by dividing by the magnitude of the product
        // pcf = ifft2(fftr./abs(fftr));
        //
        // 3.
        // Multiply one image by the conjugate of the other
        // fftr = F .* conj(G); // Matlab code
        // mag = sqrt(real(fftr).^2 + imag(fftr).^2);
        // imgReal = real(fftr) ./ mag;
        // imgImag = imag(fftr) ./ mag;
        // merge = complex(imgReal, imgImag);
        // pcf = real(ifft2(merge));
        //
        // This algorithm uses method number 3, because it can be done using OpenCV function calls and it is still pretty fast

        // Multiply spectrums of the scene and the model (use CV_DXT_MUL_CONJ to get correlation instead of convolution)
        cvMulSpectrums( sceneDFT, modelDFT, sceneDFT, CV_DXT_MUL_CONJ );

        // Split Fourier in real and imaginary parts
        cvSplit( sceneDFT, imageRe, imageIm, 0, 0 );

        // Compute the magnitude of the spectrum components: Mag = sqrt(Re^2 + Im^2)
        cvPow( imageRe, imageMag, 2.0 );
        cvPow( imageIm, imageImMag, 2.0 );
        cvAdd( imageMag, imageImMag, imageMag, NULL );
        cvPow( imageMag, imageMag, 0.5 );

        // Normalize correlation (Divide real and imaginary components by magnitude)
        cvDiv( imageRe, imageMag, imageRe, 1.0 );
        cvDiv( imageIm, imageMag, imageIm, 1.0 );
        cvMerge( imageRe, imageIm, NULL, NULL, sceneDFT );

        // inverse dft
        cvDFT( sceneDFT, sceneDFT, CV_DXT_INVERSE_SCALE, complexInput->height );
        cvSplit( sceneDFT, imageRe, imageIm, 0, 0 );

        double dMin = 0.0;
        double dMax = 0.0;
        CvPoint pntMin;
        CvPoint pntMax;
        cvMinMaxLoc( imageRe, &dMin, &dMax, &pntMin, &pntMax, NULL );

        CvFont font1;
        cvInitFont( &font1, CV_FONT_VECTOR0, 0.4f, 0.4f, 0.0f,2 );
        char strCorrelation[255];
        sprintf(strCorrelation, "%.2f", dMax);
        cvPutText( modelPolar, strCorrelation, cvPoint(10,10), &font1, CV_RGB(0,0,255));

        if ( dMax > CORRELATION_THRESHOLD )
        {
                // check pntMax location and increase or decrease dAngleCorrection from it
                int x = pntMax.x; // log range
                if (x > (imageRe->width/2))
                        x = x - imageRe->width; // positive or negative values
                int y = pntMax.y; // angle
                if (y > (imageRe->height/2))
                        y = y - imageRe->height; // positive or negative values
               
                // The correct correlation angle for log-polar depends on the ratio of the length of the y axis
                // with 360 degrees in 1 complete revolution
                dAngleCorrection = (double)y / ((double)imageRe->height/360.0);
               
                char strAngleCorrection[255];
                sprintf(strAngleCorrection, "%.2f", (double)dAngleCorrection);
                cvPutText( modelPolar, strAngleCorrection, cvPoint(10,25), &font1, CV_RGB(0,0,255));

                //char strScaleCorrection[255];
                //sprintf(strScaleCorrection, "%.2f", (double)x);
                //cvPutText( modelPolar, strScaleCorrection, cvPoint(10,40), &font1, CV_RGB(0,0,255));
        }

        cvNamedWindow("SceneOriginal", 1);
        cvNamedWindow("ModelOriginal", 1);
        cvNamedWindow("ScenePolar", 1);
        cvNamedWindow("ModelPolar", 1);

        cvShowImage("SceneOriginal", sceneImg);
        cvShowImage("ModelOriginal", modelImg);
        cvShowImage("ScenePolar", scenePolar);
        cvShowImage("ModelPolar", modelPolar);

        // Rotate the image by the amount determined above and display it
        //cv2DRotationMatrix(pntCenter, dAngleCorrection, 1, rotationMatrix);
        //cvWarpAffine(sceneImg, rotationImg, rotationMatrix, /*CV_INTER_NN |*/ CV_WARP_FILL_OUTLIERS, cvScalarAll(0));
        //cvNamedWindow("Rotation",1);
        //cvShowImage("Rotation", rotationImg);

        // I tried to reverse the log polar for this image, but it never worked.  Maybe the image was too small...
        //cvLogPolar( scenePolar, scenePolar, pntCenter, nLogPolarScale, CV_INTER_LINEAR+CV_WARP_FILL_OUTLIERS+CV_WARP_INVERSE_MAP );
        //cvNamedWindow("SceneRPolar",1);
        //cvShowImage("SceneRPolar", scenePolar);

        dAngleCorrection = (double)dAngleCorrection * M_PI / 180.0;

        return dAngleCorrection;
}

// The Matlab code to do the same thing as the imageCorrelation function
//clear all;
//close all;
//
//rot_angle = 68;
//
//%img1 = imread('reflexive174.jpg');
//img2 = img1;
//
//[height, width] = size(img1);
//center = [width/2 height];
//
//figure;
//%img1 = [img1; zeros(height, width)];
//imshow(img1);
//img1polar = logpolar(img1, width/2, height);
//figure;
//imshow(img1polar);
//img1fft = fft2(img1polar);
//figure;
//imshow(img1fft);
//
//figure;
//img2 = [img2; zeros(height, width)];
//img2 = imrotate(img2, rot_angle, 'nearest', 'crop');
//
//img2 = img2(1:height, :);
//imshow(img2);
//img2polar = logpolar(img2, width/2, height);
//figure;
//imshow(img2polar);
//img2fft = fft2(img2polar);
//figure;
//imshow(img2fft);
//
//% Three ways to calculate the phase correlation
//%1.
//% Create phase difference matrix
//%pdm = exp(i*(angle(img1fft)-angle(img2fft)));
//%pdm = exp(i*(atan2(imag(img1fft), real(img1fft)) - atan2(imag(img2fft), real(img2fft))));
//% Solve for phase correlation function
//%pcf = real(ifft2(pdm));
//
//%2.
//%Multiply one image by the conjugate of the other
//fftr = img1fft .* conj(img2fft);
//%Normalize by dividing by the magnitude of the product
//pcf = ifft2(fftr./abs(fftr));
//
//%3.
//%Multiply one image by the conjugate of the other
//%fftr = img1fft .* conj(img2fft);
//%mag = sqrt(real(fftr).^2 + imag(fftr).^2);
//%imgReal = real(fftr) ./ mag;
//%imgImag = imag(fftr) ./ mag;
//%merge = complex(imgReal, imgImag);
//%pcf = real(ifft2(merge));
//
//figure;
//plot(pcf);
//
//[max_pcf, imax] = max(abs(pcf(:)));
//[ypeak, xpeak] = ind2sub(size(pcf),imax(1))
//['Angle ' num2str(rot_angle) ' caused max value of ' num2str(max_pcf) ' located at x=' num2str(xpeak) ' y=' num2str(ypeak)]
//correlation_angle = xpeak / (width/360);
//correlation_scale = ypeak;
//
//% Use this angle and rotate the image - this should match
//img11 = [img1; zeros(height, width)];
//img12 = imrotate(img11, correlation_angle, 'nearest', 'crop');
//img12 = img12(1:height, :);
//figure;
//imshow(img12);


--- In [hidden email], "alinmallen" <alinmallen@...> wrote:

>
> Hi:
> I am useing Fourier-Mellin method to match 2 same size images for invariance to rotation and following bellow steps for the matching:
>
> 1.FFT transformation to SrcImg and TmpImg and calculate their Log-Magnitude and save as SrcImg(2) and TmpImg(2).
>
> 2.Log-Polar Transformation to SrcImg(2) and TmpImg(2) and correspondingly saved as SrcImg(3) and TmpImg(3).
>
> 3.FFT transformation to SrcImg(3) and TmpImg(3),and calculate their phase correlation peak value location to determine rotation.
>
> The problem is I always get peak value at cooridinate (0,0) on phase correlation plane.
> I have my sample image and result show in bellow album.
> http://groups.yahoo.com/group/OpenCV/photos/album/759048714/pic/list?mode=tn&order=ordinal&start=1&dir=asc
>
> Can someone help me to tell anything I did wrong on that job?
>
> Appreciate!!
>


Reply | Threaded
Open this post in threaded view
|

Re: Image matching through rotation, scale using log polar frequency domain

julirobin81
This post has NOT been accepted by the mailing list yet.
@Heather Larson-2: Can you please attach your header and cpp files because the code is missing variables definitionsm, and its unclear.

If you have some example images, that will be great.

Thank you.