2D Korrelation - Hilfe bei der Implementierung mit FFT (FFTW und OpenCV)

[BoondockDuck]

Hallo!
Ich verwende FFTW und OpenCV und versuche eine zweidimensionale Korrelation im Frequenzbereich zu implementieren. Es gibt ein Bild (image) und ein kleineren Bildausschnitt (templ). Die Position von templ im image soll gefunden werden. (Position des Maximalwertes des Lösungsraumes der Korrelation).

Ich weiß nicht was ich falsch mache. Soweit ich die Theorie verstanden habe ist eine Korrelation eine Faltung zweier Signale. Eines der beiden Signale muss zeitlich invertiert ("rumgedreht") sein. Das erreicht man entweder im Zeitbereich durch tatsächliches umdrehen (ich weiß gar nicht wie man das auf ein 2D-Signal anwenden soll?) oder im Frequenzbereich durch elementweise Multiplikation der beiden Fouriertransformierten der Signale. Eins der Signale ist dabei komplex konjugiert.

Lösungsraum R = IF ( F(image) x F*(templ) )
(IF(): inverse FFT, F() : FFT, x : elementweise multiplikation, * : komplex konjugiert)

Dabei ist zu beachten dass der Lösungsraum größer ist. Er hat die Dimensionen von Image, der relevante bereich ist aber (image.höhe-templ.höhe / image.breite-Templ.breite)

Da image und templ prinzipbedingt nicht gleich groß sind wird templ in x und y richtung mit nullen erweiterd (zero padding) um anschließend korrekt multiplizieren zu können. Die Ergebnisse des restlichen Lösungsraum kommen nach meinem Verständnis durch das "wrapping" bzw. überlappen der Faltung zustande.

Wenn ich meinen Code auf ein Bild anwende bekomme ich das richtige Ergebnis
also
Image = IF ( F ( Image ) )
Multipliziere ich die beiden FFTs wie oben beschrieben kommt Murks raus.

Das ist der Code (sorrry, lang). Bei zwei Bildern bekomme ich ein Ergebnis was sich graphisch darstellen lässt (siehe Anhang). Man erkennt dass "Wrapping" der Faltung. Vermutlich lässt sich das mit richtigem Padding beheben. Man erkennt auch ein Maximum oben links welches fast (auf 2 Pixel) genau der Position meiner Verschiebung entspricht. D.h. das Ergebnis ist fast richtig. Was ich nicht verstehe ist woher dieses Streifenmuster kommt. Es scheint als wäre das Bild falsch skaliert. Das ist es aber nicht der Fall.

Ich weiß nicht was ich falsch mache.
Und jetzt auf die plätze, fertich, los! ... hoffentlich

Das ist ech ein fasriges thema ... denn irgendwie steht nicht an jeder Ecke jemand rum den man da um Hilfe bitten kann.

Gruß
Christoph

    fftw_complex    *image_data, *image_fft_result, *templ_data, *templ_fft_result, *mul_result, *ifft_result;
    double            *mat_result;
    fftw_plan        image_plan_forward_2d, templ_plan_forward_2d, plan_backward_2d;
    int                i,j,k;
    int                size_w = image.size().width;
    int                size_h = image.size().height;
    int                size =    size_w * size_h;
    cv::Size        orig_image_size = image.size();
    cv::Size        orig_templ_size = templ.size();
    
    ofstream myfile;
    // Mittelwert von beiden Signalen abziehen
    image -= cv::mean(image);
    templ -= cv::mean(templ);
    
    //Padding of template to fit image size
    cv::Mat templ_padded(cv::Size(size_w,size_h), image.type(), cv::Scalar(0));
    cv::Mat templ_mask(templ_padded,cv::Rect(cv::Point(0,0),templ.size()));
    templ.copyTo(templ_mask);
    templ = templ_padded;
    
    // data alloc
    image_data        = ( fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * size );
    image_fft_result  = ( fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * size );
    templ_data        = ( fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * size );
    templ_fft_result  = ( fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * size );
    mul_result          = ( fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * size );
    ifft_result       = ( fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * size );
    mat_result        = ( double* )       fftw_malloc( sizeof( double )       * size );

    image_plan_forward_2d  = fftw_plan_dft_2d( size_w, size_h, image_data, image_fft_result, FFTW_FORWARD,  FFTW_ESTIMATE );
    templ_plan_forward_2d  = fftw_plan_dft_2d( size_w, size_h, templ_data, templ_fft_result, FFTW_FORWARD,  FFTW_ESTIMATE );
    plan_backward_2d       = fftw_plan_dft_2d( size_w, size_h, mul_result, ifft_result,      FFTW_BACKWARD, FFTW_ESTIMATE );
    
    // alles null setzen (weglassen ?)
    for (i = 0; i < size; i++) {
        image_data[i][0] = 0.0;
        image_data[i][1] = 0.0;
        
        templ_data[i][0] = 0.0;
        templ_data[i][1] = 0.0;
    }
    
    /* load images' data to FFTW input */
    for( i = 0, k = 0 ; i < image.size().height ; i++ ) {
        for( j = 0 ; j < image.size().width ; j++, k++ ) {
            image_data[k][0] = image.at<double>(i, j);
            image_data[k][1] = 0.0;
        }
    }
   
    /* load templ data to FFTW input */
    for( i = 0, k = 0 ; i < templ.size().height ; i++ ) {
        for( j = 0 ; j < templ.size().width ; j++, k++ ) {
            templ_data[k][0] = templ.at<double>(i, j);
            templ_data[k][1] = 0.0;
        }
    }
    
    fftw_execute( image_plan_forward_2d );
    fftw_execute( templ_plan_forward_2d );
    
    
    // komplex konjugiert von a mit b multopilizerein
    // (a - bi)*(c + di) = (ac + bd)*(ad - bc)i
    double a,b,c,d;
    for( i = 0 ; i < size ; i++ ) {
        a = templ_fft_result[i][0];
        b = templ_fft_result[i][1];
        c = image_fft_result[i][0];
        d = image_fft_result[i][1];
        
        mul_result[i][0] = a*c+b*d;
        mul_result[i][1] = a*d-b*c;
    }
    
    fftw_execute( plan_backward_2d );
    
    for( i = 0 ; i < size ; i++ ) {
        mat_result[i] = (double)(ifft_result[i][0] / size);
    }
    
    /* free memory */
    fftw_destroy_plan( image_plan_forward_2d );
    fftw_destroy_plan( templ_plan_forward_2d );
    fftw_destroy_plan( plan_backward_2d );
    
    fftw_free( image_data );
    fftw_free( image_fft_result );
    fftw_free( templ_data );
    fftw_free( templ_fft_result );
    fftw_free( ifft_result );