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

[BoondockDuck]

Hallo!
Ich verwende FFTW (http://www.fftw.org/) und OpenCV (http://opencv.willowgarage.com/wiki/) 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),te mpl.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 );

[ePyx]

Also zur Theorie :
Die KKF kann man im spektralen Bereich durch eine Multiplikation zweier Signale, wobei das eine konjungiert komplex ist, realisiert werden. Das stimmt schon einmal. Problem dabei ist meiner Meinung nach nur, dass die Verschiebung nicht größer pi/2 sein darf (siehe konjungiert komplexe zahl grafisch dargestellt).

Das gilt für den 1D- und 2D-Bereich gleichermaßen.

Hilfreich wären beide Bilder, also der Ausschnitt und das zu untersuchende Bild. Stellt die Grafik von dir den Betrag, die Phase oder einfach nur das komplexe Ergebnisfeld dar?

Hab bis jetzt immer nur im Zeitbereich die KF angewendet, daher würde ich aus dem Bauch heraus sagen, dass die beiden weißen Flecken die Punkte sind, wo der Korrelationskoeffizient sein Maximum erreich.

[BoondockDuck]

Hi ,
das oben zu sehende Bild zeigt die Inverse FFT der Korrelation dar (Lösungsraum R = IF ( F(image) x F*(templ) )) im Zeit-bzw Bildbereich dar. Die Darstellung ist ja eh nicht aussagekräftig.
Zum Vergleich habe ich die 2D Korrelation mit OpenCV-hauseigenen Methoden durchgeführt. (Anhang result_opencv.jpg). Das Bild ist im Prinzip mit meinem identisch, nur habe ich das Bild nicht abgeschnitten und bei mir sind die komischen Streifen drin.

Wie soll ich mir denn die Verschiebung um Pi/2 im 2D - Fall vorstellebn? Aber selbst eine Verschiebung um wenige Pixel ändert nichts am Streifenmuster.

Anbei auch mal das Testbild und den Ausschnitt.

Gruß
Christoph

[ePyx]

Wie soll ich mir denn die Verschiebung um Pi/2 im 2D - Fall vorstellebn?

die DFT ist ja nur eine Differenzengleichung (Summe) und wird im 2D-Fall nicht wirklich zweidimensional angewendet. Viel eher wird sie auf die Spalten und anschließend auf die Zeilen deines Bildes angewendet.

Zeitbereich/Diskreter Bereich:
Stellt man sich den Farbverlauf einer Zeile als Funktion in Abhängigkeit der Zeit oder des Indexes vor, so ist auch dort eine Phasenverschiebung möglich.

Dein Ergebnisbild lässt darauf schließen, dass du die zyklische Fortsetzung der FFT nicht unterdrückst. Meinst du das mit Abschneiden ?

[BoondockDuck]

Dein Ergebnisbild lässt darauf schließen, dass du die zyklische Fortsetzung der FFT nicht unterdrückst. Meinst du das mit Abschneiden ?Genau das meine ich. Aber das sollte doch eigentlich kein Problem ergeben wie bei mir zu sehen?

[ePyx]

Naja es sieht zu mindestens periodisch aus.



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;
}
}


Macht doch eigentlich so mehr Sinn oder ?



for( i = 0 ; i < image.size().height ; i++ ) {
for( j = 0, k = 0; j < image.size().width ; j++, k++ ) {
image_data[k][0] = image.at<double>(i, j);
image_data[k][1] = 0.0;
}
}


Schließlich durch läufst du zeilenweise und pro Zeile schaufelst du Daten. wenn k nicht zurückgesetzt wird, könnte das die Ursache sein.

[BoondockDuck]

Ich wäre ja froh wenn es das wäre.

Die Daten liegen sequentiell im Speicher. K ist der Zähler für das kontinuierliche Array, i und j die für Zeilen und Spalten.
Dabei gilt : k = i * image.size().width + j


Das hier:



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;
}
}


entspricht dem hier:



/* load images' data to FFTW input*/
for( k = 0 ; k < size ; k++ ) {
image_data[k][0] = ((double *) image.data)[k];
image_data[k][1] = 0.0;
}


Der fancy image.at operator erspart einem nur den Pointer-Typecast und man kann direkt Pixel adressieren ohne sich über die Repräsentation im Speicher im klaren zu sein. (Deswegen muss hier auch explizit dach (double*) gecasted werden weil ich den Datentyp kenne.) Zweiteres ist halt schneller aber hässlicher.

[ePyx]

Naja, um ehrlich zu sein, ist die zweite weniger hässlich. Da sieht man wenigstens gleich was Sache ist.

Jedenfalls stimmt irgend etwas in deinen Zeilen nicht. Entweder beim Zusammensetzen oder schon beim Ergebnis der FFT benutzt du deine Hilfsvariablen a bis d garnicht ? Wunderte mich gerade darüber, dass du dort im Kommentar schreibst, dass du die konjungiert-komplexe Form berechnen willst, es aber nicht tust. Oder bin ich blind?

[BoondockDuck]

Oder bin ich blind?
Nein vollkommen richtig... Copy und Paste Fehler, Ich habe es oben korrigiert.

Und ... Ich hab den Fehler gefunden :p. Der Code funktioniert, es hing (hängt) am theretischen Verständnis. Das Bild (image) muss anscheinend eine Zweierpotenz als Seitenlänge aufweisen. Also 2,4,8,16 ... px hoch und breit.
Ich hatte das bisher überlesen oder ignoriert.

Soweit so schön. Und jetzt geh ich ins Bett.

[ePyx]

Ok, das hatte ich irgendwie auch vorausgesetzt. Da die FFT ja bekanntlich nur für 2er Potenzen richtig funktioniert. Das würde dann die Artefakte erklären. Aber schön, dass es funktioniert.