2D Cross Correlation mit FFTW

Snotgun

Mitglied
Hallo!

Ich bastel gerade an einem Projekt, in dem ich in Bildern bestimmte Features erkennen muss. Ich würde gerne über eine Kreuzkorrelation im Frequenzraum ermitteln, wo bestimmte Features liegen. Geschrieben ist mein Programm bis dato in C#, als Wrapper verwende ich http://www.sdss.jhu.edu/~tamas/bytes/fftwcsharp.html

Die FFT selbst funktioniert schon halbwegs zuverlässig. Ich kriege Bilder in den Frequenzraum und auch wieder zurücktransformiert - und zwar sowohl das Eingangsbild als auch das kleinere Bild, das die zu findenden Features darstellt. Problematisch wird es bei der Convolution / Cross Correlation. Wenn ich das richtig verstanden habe (und sicher bin ich mir inzwischen nirgendwo mehr), müssen die Bilder im Frequenzraum miteinander multipliziert werden, Zahl für Zahl. Entsprechende Tutorials sagen das gleiche und zwei komplexe Zahlen zu multiplizieren ist eigentlich auch nicht so schwer. Merkwürdig ist dabei mein Ergebnis, das scheinbar nur ein Rauschen ist. Die Werte explodieren dabei teilweise in den siebenstelligen Wertebereich, deswegen gehen beim Eintragen in das byte-Array wahrscheinlich alle Muster vollends verloren.

Falls jemand schon mal etwas ähnliches gemacht hat und/oder Ahnung von der Materie hat, wäre ich für jede Idee dankbar, wo mein Fehler liegen könnte.

Meine Umrechnungsmethode sieht wie folgt aus. Ich hoffe das ist nicht zu viel Code und die Kommentare helfen ein wenig beim Verständnis:

Code:
public override void Apply(byte[] toFilter, BitmapInfo info)
        {

            for (int i = 0; i < findWidth; i++) //Kopiere Bild, das gefunden werden soll, in größeres Array
            {
                for (int j = 0; j < findHeight; j++)
                {
                    ffind[(i * info.Height + j) * 2] = imageToFind[i * findHeight + j];
                }
            }

            int n = toFilter.Length / bytesPerPixel; //Anzahl der Pixel

            toRowMajorOrder(toFilter, fin, info.Width, info.Height, info.Stride, bytesPerPixel); //Kopiere Bild in "Row Major" Anordnung, also x und y vertauscht

            Marshal.Copy(fin, 0, pin, n * 2); // Kopiere Eingangsbild in unmanaged array
            Marshal.Copy(ffind, 0, pfind, n * 2); // Kopiere Bild, das gefunden werden soll, in unmanaged array

            fplan1 = fftwf.dft_2d(info.Width, info.Height, pin, pout, fftw_direction.Forward, fftw_flags.Estimate); //Plan zum Übersetzen des großen Bildes in Frequenzraum
            fplanfind = fftwf.dft_2d(info.Width, info.Height, pfind, pfindout, fftw_direction.Forward, fftw_flags.Estimate); // Plan zum Übersetzen des kleinen Bildes in Frequenzraum
            fplan2 = fftwf.dft_2d(info.Width, info.Height, pout, pin, fftw_direction.Backward, fftw_flags.Estimate); // Plan zum Rückübersetzen des großen Bildes
            fplanfind2 = fftwf.dft_2d(info.Width, info.Height, pfindout, pfind, fftw_direction.Backward, fftw_flags.Estimate); // Plan zum Rückübersetzen des kleinen Bildes, wird nicht genutzt, nur zum Debuggen bei Bedarf

            fftwf.execute(fplan1);
            fftwf.execute(fplanfind);

            Marshal.Copy(pout, fout, 0, n * 2); // Kopiere großes Bild in gemanageten Bereich
            Marshal.Copy(pfindout, ffindout, 0, n * 2); // Analog das kleine Bild

            // Multiplikation der beiden frisch zurückgeholten Bilder im Frequenzraum
            // Achtung: Komplexe Zahlen am Werk! Jedes zweite Element ist imaginärer Anteil!
            float[] tempResult = new float[ffindout.Length];
            float scale = 1f / (float)(info.Width * info.Height);
            for (int i = 0; i < ffindout.Length; i += 2)
            {
                tempResult[i] = (fout[i] * ffindout[i]
                                    - fout[i + 1] * ffindout[i + 1]) * scale;
                tempResult[i + 1] = (fout[i] * ffindout[i + 1]
                                    + fout[i + 1] * ffindout[i]) * scale;
            }


            Marshal.Copy(tempResult, 0, pout, n * 2); // Zurückkopieren des Multiplikationsergebnisses in ungemanageten Bereich

            fftwf.execute(fplan2); // Rückübersetzen aus dem Frequenzraum

            // Aufräumen
            fftwf.destroy_plan(fplan1);
            fftwf.destroy_plan(fplan2);
            fftwf.destroy_plan(fplanfind);
            fftwf.destroy_plan(fplanfind2);

            Marshal.Copy(pin, fout, 0, n * 2); // Zurückholen des Ergebnisses in gemanageten Bereich

            toColumnMajorOrder(fout, OutputData, info.Width, info.Height, info.Stride, bytesPerPixel); // Alles noch Row Major, muss umstrukturiert werden
        }
 
A

Amegon

Hi, habe ein ähnliches Programm und die Crosscorrelation klappt auhc soweit (nutze auhc einen Wrapper, aber dann ansteuerung per Vb, und genauso einen wrapper für OpenCV).
Bei mir sieht die CrossCorrelation folgendermaßen aus:

Code:
    Private Function CrossCorrelation(ByVal F As Double(,), ByVal G As Double(,)) As Double(,)
        Dim _lastElement As Integer = F.GetUpperBound(0)

        Dim _result(_lastElement, 1) As Double 'bei VB gibt man hier oberstes Element an, nciht count
'bei dir wäre das _result(_lastElement+1, 2)

        For i As Integer = 0 To _lastElement
            Dim a As Double = G(i, 0) * F(i, 0) + G(i, 1) * F(i, 1)
            Dim b As Double = -G(i, 0) * F(i, 1) + G(i, 1) * F(i, 0)
 'Laufzeit ware mehr als 3 mal so hoch, als ich a^2+b^2 geschrieben habe
            Dim c As Double = Math.Sqrt(a * a + b * b)
            _result(i, 0) = a / c
            _result(i, 1) = b / c
        Next

        Return _result

    End Function

und wenn du das Fertige Bild erhälst, dann ist noch folgendes notwendig:
Code:
 Private Function Normalize(ByVal vector As Double(,), Optional ByVal removeImgPart As Boolean = False) As Double(,)
        Dim _result(vector.GetUpperBound(0), vector.GetUpperBound(1)) As Double
        Dim _size As Integer = vector.GetUpperBound(0) + 1
        If removeImgPart Then
            For i As Integer = 0 To vector.GetUpperBound(0)
                _result(i, 0) = vector(i, 0) / _size
                _result(i, 1) = 0
            Next
        Else
            For i As Integer = 0 To vector.GetUpperBound(0)
                _result(i, 0) = vector(i, 0) / _size
                _result(i, 1) = vector(i, 1) / _size
            Next
        End If
        Return _result
    End Function
wobei die Fallunterscheidung nicht wirklcih relevant ist. Wichtig ist nur, das jeder Wert standardmäßig ein x-faches ist. und dieser faktor ist ganz genau von höhexbreite vom bild abhängig. wenn du 1000 pixel insgesamt hast, so hast du werte von 0-255 -> 0 bis 255000 oder bei der eingabe in den bereichen 0-1 -> 0 bis 1000

viel erfolg beim anwenden, und wenn du weiterkommst, wie man die kreuzkorrelation noch genau genug anwendet, wenn die bilder unterschiedlcih groß sind, dann bitte bscheid geben :)

zu den rrechenfehlern bei dir, ich gehe davon aus, das es folgendermaßen der richtige weg ist:
bei kreuzkorrellation jedes pixel durch Wurzel(real^2+imag^2) teilen
dann inverse fouriertranformation, und dann erst die normalisierung, indem du jedes pixel durch höhe*breite vom bild teilst. Da der vektor des Bildes genau höhe*breite elemente hat bzw. wenn du einen 1 dimensionalen vektor benutzt, dann hat er ja höhe*breite*2 drin, daher musst du dann nur durch anzahl elemente oder anzahlelemente/2 teilen
 
Zuletzt bearbeitet von einem Moderator: