Ask Your Question
1

phase correlation for image registration(image stitching)

asked 2012-08-20 04:27:44 -0600

mrgloom gravatar image

updated 2019-12-09 07:56:23 -0600

Akhil Patel gravatar image

I'm trying to stitch 2 images using cross correlation (phase correlation).Images are the same size.Only shift is present. I tried opencv with cvDFT and opencv+FFTW it seems to work, but for some reason coorelation peak coordinate is not the shift coordinate, maybe it depends on quadrant where correlation point is. So the question is how to obtain shift point from correlation point.

Note: I don't need "interest points" SIFT-like approaches.

My suggestion commented in code. Here is the code that I use:

class Peak
{
public:
    CvPoint pt;
    double  maxval;
};
Peak old_opencv_FFT(IplImage* src,IplImage* temp)
{
    CvSize imgSize = cvSize(src->width, src->height);
    // Allocate floating point frames used for DFT (real, imaginary, and complex) 
    IplImage* realInput = cvCreateImage( imgSize, IPL_DEPTH_64F, 1 ); 
    IplImage* imaginaryInput = cvCreateImage( imgSize, IPL_DEPTH_64F, 1 ); 
    IplImage* complexInput = cvCreateImage( imgSize, IPL_DEPTH_64F, 2 );
    int nDFTHeight= cvGetOptimalDFTSize( imgSize.height ); 
    int nDFTWidth= cvGetOptimalDFTSize( imgSize.width ); 
    CvMat* src_DFT = cvCreateMat( nDFTHeight, nDFTWidth, CV_64FC2 ); 
    CvMat* temp_DFT = cvCreateMat( nDFTHeight, nDFTWidth, CV_64FC2 );
    CvSize dftSize = cvSize(nDFTWidth, nDFTHeight);
    IplImage* imageRe = cvCreateImage( dftSize, IPL_DEPTH_64F, 1 ); 
    IplImage* imageIm = cvCreateImage( dftSize, IPL_DEPTH_64F, 1 );
    IplImage* imageImMag = cvCreateImage( dftSize, IPL_DEPTH_64F, 1 ); 
    IplImage* imageMag = cvCreateImage( dftSize, IPL_DEPTH_64F, 1 ); 

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

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

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

    // Split Fourier in real and imaginary parts 
    cvSplit(src_DFT,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,src_DFT); 

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

    double minval = 0.0; 
    double maxval = 0.0; 
    CvPoint minloc; 
    CvPoint maxloc; 
    cvMinMaxLoc(imageRe,&minval,&maxval,&minloc,&maxloc,NULL); 

    int x=maxloc.x; // log range 
    //if (x>(imageRe->width/2)) 
    //        x = x-imageRe->width; // positive or negative values 
    int y=maxloc.y; // angle 
    //if (y>(imageRe->height/2)) 
    //        y = y-imageRe->height; // positive or negative values 

    Peak pk;
    pk.maxval= maxval;
    pk.pt=cvPoint(x,y);
    return pk;
}
void phase_correlation2D( IplImage* src, IplImage *tpl, IplImage *poc )
{
    int     i, j, k;
    double  tmp;

    /* get image properties */
    int width    = src->width ...
(more)
edit retag flag offensive close merge delete

2 answers

Sort by ยป oldest newest most voted
3

answered 2012-08-22 08:57:55 -0600

Michael Burdinov gravatar image

updated 2012-08-30 01:30:12 -0600

Here is how I did it (and it worked properly):

1) Create 2 images (fft1 and fft2) according to getOptimalDFTSize().

int width = getOptimalDFTSize(max(image1.cols,image2.cols));
int height = getOptimalDFTSize(max(image1.rows,image2.rows));
Mat fft1(Size(width,height),CV_32F,Scalar(0));
Mat fft2(Size(width,height),CV_32F,Scalar(0));

2) Filled them with values from original images (image1 and image2).

3) Do the following calculations:

dft(fft1,fft1,0,image1.rows);
dft(fft2,fft2,0,image2.rows);
mulSpectrums(fft1,fft2,fft1,0,true);
idft(fft1,fft1);
double maxVal;
Point maxLoc;
minMaxLoc(fft1,NULL,&maxVal,NULL,&maxLoc);
resX = (maxLoc.x<width/2) ? (maxLoc.x) : (maxLoc.x-width);
resY = (maxLoc.y<height/2) ? (maxLoc.y) : (maxLoc.y-height);

Also your question is not about 'image stitching' but about a step done before that - 'image registration'. Title of your question is somewhat misleading.

Edit:

mrgloom, it is already almost full code. The only part I didn't posted yet is how I copied the data between the images. This is how I did it:

for(j=0; j<image1.rows; j++)
for(i=0; i<image1.cols; i++)
fft1.at<float>(j,i) = image1.at<unsigned char>(j,i);

for(j=0; j<image2.rows; j++)
for(i=0; i<image2.cols; i++)
fft2.at<float>(j,i) = image2.at<unsigned char>(j,i);

I didn't knew back than about functions like convertTo.

edit flag offensive delete link more

Comments

how you copy image? I tried you method but it doesn't work, maybe I.m doing something wrong, look my update.

mrgloom gravatar imagemrgloom ( 2012-08-28 01:51:13 -0600 )edit

Your third option seems to be correct. Don't know why it is not working.

Michael Burdinov gravatar imageMichael Burdinov ( 2012-08-29 01:19:24 -0600 )edit

maybe I need to use imaginary part?

mrgloom gravatar imagemrgloom ( 2012-08-29 08:24:58 -0600 )edit

It worked properly as is in my program.

Michael Burdinov gravatar imageMichael Burdinov ( 2012-08-29 12:24:18 -0600 )edit

Maybe it doesn't work because my images have only 20% overlap? Can you post full code?

mrgloom gravatar imagemrgloom ( 2012-08-30 01:05:22 -0600 )edit

See my edit. Couldn't post code in comment.

Michael Burdinov gravatar imageMichael Burdinov ( 2012-08-30 01:31:09 -0600 )edit

here your code that I use http://codepaste.ru/12808/ and it give me (-60,-22) or (60,22) , but ImageJ stittching 2D plugin gives me (-20,-23) or (20,23). here are images http://dl.dropbox.com/u/8841028/FFT%20template%20matching/test%20images/im1.pnghttp://dl.dropbox.com/u/8841028/FFT%20template%20matching/test%20images/im2.png

mrgloom gravatar imagemrgloom ( 2012-12-13 05:27:25 -0600 )edit

Hi mrgloom,

Did you keep the full code? I try run this code, but doesn't work.. I wan't to implement this code on the case parallax object. Thanks.

Rgds, Mohammad Badrullami

Mohammad Badrullami gravatar imageMohammad Badrullami ( 2015-06-02 20:54:04 -0600 )edit

Hello, I am trying to get your code, but the page is not able to open. I have also tried to copy whole code, but which files have to import in code, it is missing. So, I got error as:

64  IntelliSense: identifier "fftw_complex" is undefined    c:\Users\Administrator\Documents\Visual Studio 2013\Projects\OM\OM\ImageStiching_Using_PhaseCorrelation.cpp 124 2   OM

ans so on.

Kindly provide me link from where I will get whole code.

Thank you.

darshana.mistry gravatar imagedarshana.mistry ( 2016-02-24 01:00:24 -0600 )edit
0

answered 2016-03-31 20:34:37 -0600

g9164313 gravatar image

https://github.com/g9164313/sample-re...

try this, it will translate and rotate reference image to fit source image. I use only functions provided by OpenCV.

edit flag offensive delete link more

Question Tools

Stats

Asked: 2012-08-20 04:27:44 -0600

Seen: 14,367 times

Last updated: Feb 24 '16