Форум программистов [Powered by Invision Power Board]

Цитата

Если я правильно понял вас, то ваша задача - это поиск шаблона в изображении

Цитата

Тогда, если мы сделаем сравнение двух полученных векторов, то это получится ответ, как если бы шаблон приложили в точку (0;0). Если мы сдвинем вектор шаблона на единицу, то это будет соответствовать положению шаблона (0;1), и т.д.

Код


/* obtain the FFT of img1 */
    fftw_execute( fft_img1 );
    
    /* obtain the FFT of img2 */
    fftw_execute( fft_img2 );
    
    /* obtain the cross power spectrum */
    for( i = 0; i < fft_size ; i++ ) {
        res[i][0] = ( img2[i][0] * img1[i][0] ) - ( img2[i][1] * ( -img1[i][1] ) );
        res[i][1] = ( img2[i][0] * ( -img1[i][1] ) ) + ( img2[i][1] * img1[i][0] );

        tmp = sqrt( pow( res[i][0], 2.0 ) + pow( res[i][1], 2.0 ) );

        res[i][0] /= tmp;
        res[i][1] /= tmp;
    }

    /* obtain the phase correlation array */
    fftw_execute(ifft_res);

    /* normalize and copy to result image */
    for( i = 0 ; i < fft_size ; i++ ) {
        poc_data[i] = res[i][0] / ( double )fft_size;
    }

Цитата

так все таки надо смотреть все возможные комбинации положения шаблона?

Цитата

еще смущает что 1D массивы не проссаживается ли производительность против 2D массива? (в плане считать фурье туда сюда)

Код

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;
    int height   = src->height;
    int step     = src->widthStep;
    int fft_size = width * height;

    /* setup pointers to images */
    uchar   *src_data = ( uchar* ) src->imageData;
    uchar   *tpl_data = ( uchar* ) tpl->imageData;
    double  *poc_data = ( double* )poc->imageData;

    /* allocate FFTW input and output arrays */
    fftw_complex *img1 = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * width * height );
    fftw_complex *img2 = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * width * height );
    fftw_complex *res  = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * width * height );   
    
    /* setup FFTW plans */
    fftw_plan fft_img1 = fftw_plan_dft_2d( height ,width, img1, img1, FFTW_FORWARD,  FFTW_ESTIMATE );
    fftw_plan fft_img2 = fftw_plan_dft_2d( height ,width, img2, img2, FFTW_FORWARD,  FFTW_ESTIMATE );
    fftw_plan ifft_res = fftw_plan_dft_2d( height ,width, res,  res,  FFTW_BACKWARD, FFTW_ESTIMATE );
    
    /* load images' data to FFTW input */
    for( i = 0, k = 0 ; i < height ; i++ ) {
        for( j = 0 ; j < width ; j++, k++ ) {
            img1[k][0] = ( double )src_data[i * step + j];
            img1[k][1] = 0.0;
            
            img2[k][0] = ( double )tpl_data[i * step + j];
            img2[k][1] = 0.0;
        }
    }
    
    ///* Hamming window */
    //double omega = 2.0*M_PI/(fft_size-1);
    //double A= 0.54;
    //double B= 0.46;
    //for(i=0,k=0;i<height;i++)
    //{
    //    for(j=0;j<width;j++,k++)
    //    {
    //        img1[k][0]= (img1[k][0])*(A-B*cos(omega*k));
    //        img2[k][0]= (img2[k][0])*(A-B*cos(omega*k));
    //    }
    //}

    /* obtain the FFT of img1 */
    fftw_execute( fft_img1 );
    
    /* obtain the FFT of img2 */
    fftw_execute( fft_img2 );
    
    /* obtain the cross power spectrum */
    for( i = 0; i < fft_size ; i++ ) {
        res[i][0] = ( img2[i][0] * img1[i][0] ) - ( img2[i][1] * ( -img1[i][1] ) );
        res[i][1] = ( img2[i][0] * ( -img1[i][1] ) ) + ( img2[i][1] * img1[i][0] );

        tmp = sqrt( pow( res[i][0], 2.0 ) + pow( res[i][1], 2.0 ) );

        res[i][0] /= tmp;
        res[i][1] /= tmp;
    }

    /* obtain the phase correlation array */
    fftw_execute(ifft_res);

    //normalize and copy to result image
    for( i = 0 ; i < fft_size ; i++ ) {
        poc_data[i] = res[i][0] / ( double )fft_size;
    }

    /* deallocate FFTW arrays and plans */
    fftw_destroy_plan( fft_img1 );
    fftw_destroy_plan( fft_img2 );
    fftw_destroy_plan( ifft_res );
    fftw_free( img1 );
    fftw_free( img2 );
    fftw_free( res );
}

Peak FFTW_test(IplImage* src,IplImage* temp)
{
    clock_t start=clock();

    int t_w=temp->width;
    int t_h=temp->height;

    /* create a new image, to store phase correlation result */
    IplImage* poc = cvCreateImage( cvSize(temp->width,temp->height ), IPL_DEPTH_64F, 1 );

    /* get phase correlation of input images */
    phase_correlation2D( src, temp, poc );

    /* find the maximum value and its location */
    CvPoint minloc, maxloc;
    double  minval, maxval;
    cvMinMaxLoc( poc, &minval, &maxval, &minloc, &maxloc, 0 );

    /* IplImage* poc_8= cvCreateImage( cvSize(temp->width, temp->height ), 8, 1 );
    cvConvertScale(poc,poc_8,(double)255/(maxval-minval),(double)(-minval)*255/(maxval-minval));
    cvSaveImage("poc.png",poc_8); */

    cvReleaseImage( &poc );

    clock_t end=clock();

    int time= end-start;

    //fprintf( stdout, "Time= %d using clock() \n" ,time/*dt*/ );
    //fprintf( stdout, "Maxval at (%d, %d) = %2.4f\n", maxloc.x, maxloc.y, maxval );

    CvPoint pt;
    pt.x= maxloc.x;
    pt.y= maxloc.y;
    //4 случая локации точки максимальной кореляции
    //if(maxloc.x>=0&&maxloc.x<=t_w/2&&maxloc.y>=0&&maxloc.y<=t_h/2)
    //{
    //    pt.x= src->width-maxloc.x;
    //    pt.y= -maxloc.y;
    //}
    //if(maxloc.x>=t_w/2&&maxloc.x<=t_w&&maxloc.y>=0&&maxloc.y<=t_h/2)
    //{
    //    pt.x= src->width-maxloc.x;
    //    pt.y= src->height-maxloc.y;
    //}
    //if(maxloc.x>=0&&maxloc.x<=t_w/2&&maxloc.y>=t_h/2&&maxloc.y<=t_h)
    //{
    //    /*pt.x= -maxloc.x;
    //    pt.y= -maxloc.y;*/
    //    pt.x= src->width-maxloc.x; //тут непонятно
    //    pt.y= src->height-maxloc.y;
    //}    
    //if(maxloc.x>=t_w/2&&maxloc.x<=t_w&&maxloc.y>=t_h/2&&maxloc.y<=t_h)
    //{
    //    pt.x= -maxloc.x;
    //    pt.y= src->height-maxloc.y;
    //}

    Peak pk;
    pk.maxval= maxval;
    pk.pt=pt;
    return pk;
}

Код

// adjust shift relative to image center...
    Point2d center((double)padded1.cols / 2.0, (double)padded1.rows / 2.0);

Код

#include "stdafx.h"

#include <opencv.hpp>

using namespace cv;
using namespace std;

int _tmain(int argc, _TCHAR* argv[])
{
    Mat image1= imread("1.png",0);
    Mat image2= imread("2.png",0);

    Mat r1;
    image1.convertTo(r1,CV_64F);
    Mat r2;
    image2.convertTo(r2,CV_64F);

    Point2d phaseShift1;
    Point2d phaseShift2;

    //check size
    if(r1.cols!=r2.cols||r1.rows!=r2.rows)
    {
        //image1 size 883x646
        //image2 size 882x599

        //gives peak -423 -1 (423 1)
        int n_cols= 1024;
        int n_rows= 1024;

        //gives peak  224 -1 (-224 1)
        int n_cols= max(r1.cols,r2.cols);
        int n_rows= max(r1.rows,r2.rows);

        Mat r1_pad;
        copyMakeBorder(r1,r1_pad,0,n_rows-r1.rows,0,n_cols-r1.cols, BORDER_CONSTANT, Scalar::all(0));
        Mat r2_pad;
        copyMakeBorder(r2,r2_pad,0,n_rows-r2.rows,0,n_cols-r2.cols, BORDER_CONSTANT, Scalar::all(0));

        //imwrite("b1.png",r1_pad);
        //imwrite("b2.png",r2_pad);

        //Mat hann;
        //createHanningWindow(hann, r1_pad.size(), CV_64F);
        phaseShift1 = phaseCorrelate(r1_pad, r2_pad/*, hann*/);
        phaseShift2 = phaseCorrelate(r2_pad, r1_pad/*, hann*/);
        int f=0;
    }
    else
    {
        //Mat hann;
        //createHanningWindow(hann, r1.size(), CV_64F);
        phaseShift1 = phaseCorrelate(r1, r2/*, hann*/);
        phaseShift2 = phaseCorrelate(r2, r1/*, hann*/);
    }

    return 0;
}

Код

cv::Point2d cv::phaseCorrelateRes(InputArray _src1, InputArray _src2, InputArray _window, double* response)
{
    Mat src1 = _src1.getMat();
    Mat src2 = _src2.getMat();
    Mat window = _window.getMat();

    CV_Assert( src1.type() == src2.type());
    CV_Assert( src1.type() == CV_32FC1 || src1.type() == CV_64FC1 );
    CV_Assert( src1.size == src2.size);

    if(!window.empty())
    {
        CV_Assert( src1.type() == window.type());
        CV_Assert( src1.size == window.size);
    }

    int M = getOptimalDFTSize(src1.rows);
    int N = getOptimalDFTSize(src1.cols);

    Mat padded1, padded2, paddedWin;

    if(M != src1.rows || N != src1.cols)
    {
        copyMakeBorder(src1, padded1, 0, M - src1.rows, 0, N - src1.cols, BORDER_CONSTANT, Scalar::all(0));
        copyMakeBorder(src2, padded2, 0, M - src2.rows, 0, N - src2.cols, BORDER_CONSTANT, Scalar::all(0));

        if(!window.empty())
        {
            copyMakeBorder(window, paddedWin, 0, M - window.rows, 0, N - window.cols, BORDER_CONSTANT, Scalar::all(0));
        }
    }
    else
    {
        padded1 = src1;
        padded2 = src2;
        paddedWin = window;
    }

    Mat FFT1, FFT2, P, Pm, C;

    // perform window multiplication if available
    if(!paddedWin.empty())
    {
        // apply window to both images before proceeding...
        multiply(paddedWin, padded1, padded1);
        multiply(paddedWin, padded2, padded2);
    }

    // execute phase correlation equation
    // Reference: http://en.wikipedia.org/wiki/Phase_correlation
    dft(padded1, FFT1, DFT_REAL_OUTPUT);
    dft(padded2, FFT2, DFT_REAL_OUTPUT);

    mulSpectrums(FFT1, FFT2, P, 0, true);

    magSpectrums(P, Pm);
    divSpectrums(P, Pm, C, 0, false); // FF* / |FF*| (phase correlation equation completed here...)

    idft(C, C); // gives us the nice peak shift location...

    fftShift(C); // shift the energy to the center of the frame.

    // locate the highest peak
    Point peakLoc;
    minMaxLoc(C, NULL, NULL, NULL, &peakLoc);

    // get the phase shift with sub-pixel accuracy, 5x5 window seems about right here...
    Point2d t;
    t = weightedCentroid(C, peakLoc, Size(5, 5), response);

    // max response is M*N (not exactly, might be slightly larger due to rounding errors)
    if(response)
        *response /= M*N;

    // adjust shift relative to image center...
    Point2d center((double)padded1.cols / 2.0, (double)padded1.rows / 2.0);

    return (center - t);
}