mikrocontroller.net

Forum: Digitale Signalverarbeitung / DSP Radontransformation über DFT berechnen


Autor: Alexander Liebhold (lippi2000)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hey Leute, brauch mal ne gedankliche Hilfestellung.

Mir ist klar wie die Radontransformation funktioniert und kann diese 
auch interpretieren. In DFT/FFT bin ich auch fit. Nur fehlt mir die 
Vorstellung, warum eine Radontransformation über eine 2D FFT 
anschließenden Koordinatentransformation (kartesisch -> polar + 
Interpolation) und anschließende 1D IFFT berechenbar ist. Also von den 
mathematischen Formelumstellung ist das klar.

Nur wie man sich den Zusammenhang vorstellt ist mir noch nicht 100%'ig 
klar. Liegt aber alles am Fourier Zentralscheiben Theorem.

Hat jemand ne gute Erklärung?

Gruß Alexander.

Autor: Nils (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo Alexander,

ich versuche mal Deine Fragestellung zu raten, indem ich das für mich 
anhand meiner alten Literatur noch mal nachvollziehe.

Zunächst: die Grundlage der Radon-Transformation ist das 
Projektionstheorem.
Die Radon-Transormierte ist ja die Gesamtheit aller Projektionen, z. B. 
bei einer CT.

* Nur kurz zusammengefasst, damit die Begriffe klar sind:
1) Grundlage ist zunächst die Tatsache, dass die 
Absorptionskoeffizienten eines Körpers bestimmt werden können, wenn man 
den Körper von allen Seiten durchleuchtet, das Intensitätsprofil 
aufzeichnet und das sich ergebene Gleichungssystem löst.
Der diskrete Fall ...
I_0 -> |µ11|µ12| -> I_1
I_0 -> |µ21|µ22| -> I_2
führt nach Drehung um jeweils 90° auf ein Gleichungssystem:
I_1 = I_0 exp(-(µ11-µ12)dx)
I_2 = I_0 exp(-(µ21+µ22)dx)
... (Kombinationen für alle Winkel)
Die Nachteile dieser diskreten Betrachtungsweise (großes Gl.-System, 
num. Instabilität) führen zu 2).
2) Im kontinuierlichen Fall werden diese Gleichungen zu Scharen aus 
Integralen längs eines winkelabhängigen Weges (Projektionen):
Intensitätsfunktion(Abstand, Winkel) = Integral(Dämpfungsfunktion(x,y)) 
längs des winkelabhängigen Weges
Also durch die Messung gegeben:
Eine Intensitätsfunktion, die von (Abstand, Winkel) abhängt
Zu bestimmen:
Die Dämpfungsfunktion der durchleuchteten Fläche in Abhängigkeit der 
kartes. Koordinaten (x,y)
Das Linienintegral (besser die Schar) lässt sich nun nach 
Koordinatentransformation in ein Produkt von Fouriertransformationen 
schreiben.
Soweit das Bekannte.

* Interpretationen:

Projektionstheorem
Es geht um Transformationen. Transformationen vermitteln zwischen 
mathematischen Räumen. Hier dem Projektionsraum und dem 
Transformationsraum.
Die Fouriertransformierten FÜR einen gegeben Winkel im Projektionsraum 
entsprechen der Fouriertransformierten im Transformationsraum ENTLANG 
des gegebenen Winkels - das ist das Projektionstheorem. Im Klartext: 
Führe eine FT für einen Messwert FÜR bestimmten Winkel durch und Du 
erhältst die FT des Intensitätsprofils ENTLANG dieses Winkels. 
(Natürlich immer in der Gesamtheit der Messwerte gesehen).

Koordinatensysteme
Sie haben mit dem eigentlichen Geschehen zunächst nicht viel zu tun. Die 
Messung des Intensitätsprofils erfolgt in Polarkoordinaten; für die 
Berechnung der Rückprojektion sind kartes. Koordinaten erforderlich, 
weil die DFT im Kartesischen 'spielt', d. h. das DFT-Raster muss auf 
Werte 'treffen', die aus einer Messung per Rotation stammen. Daher 
werden die Intensitätsprofile aus den Messungen interpoliert.

Was wird daraus?
Wir reden nun von inversen Transformationen, d. h. aus Projektionen 
(Intensitätsmessungen) werden Rekonstruktionen eines 2-Dim. Gebiets 
(Rückprojektionen). Das war ja das Ziel - ein Objekt von allen Seiten 
durchleuchten und dann auf sein Innenleben zu schließen.
Das das Linienintegral - also die kontinuierliche Schreibweise des 
diskreten Gleichungssystems als Produkt zweier FTs schreibbar war, wurde 
schnell klar. Durch die Berücksichtigung unterschiedlicher 
Koordinatensysteme kamen in den FTs radiusabhängige Koeffizienten ins 
Spiel und wir so zu unserer Radon-Transformierten. Dabei mussten wir 
Federn lassen: Die ursprüngliche Integration über -Unendlich bis 
+Unendlich im kartes. Raum, zwang uns schrittweise zu 
Integrationsbereichen 0..2Pi, -Pi..Pi, bis zu Koeffizienten K in 
Abhängigkeit des Durchleuchtungsradius.

Was soll das olle K in den FTs?
Wir reden ja nun eigentlich über ein Intergral mal Koeffizienten K. 
Diese K sind nichts anderes als Filter. Per Rückprojektionen projizieren 
wir Messwertprofile in dem xy-Gebiet, indem wir die 
Fouriertransformierte mit Filterkoeffizienten multiplizieren und die 
Gesamtheit überlagern.

Etwas elektrotechnischer oder warum die Koeffizienten K Filter sind:
Die Messung der Intensitäten des rotierenden Gebiets ergibt nichts 
anderes als ein Spektrum. Dieses Spektrum soll uns etwas über das 
Innenleben (die infinitesimalen Absorbtionskoeffizeinten) des Gebietes 
sagen. Das ist Projektionstheorem. Wir finden ein Spektrum im 
Ortsbereich, denn unsere Messung ist eine Messung Amplitude/Ort. Die 
(dikretisierten) Bereiche sind Amplituden im Spektrum multipliziert mit 
Filterkoeffizienten K. Das ist eine Filterung.

Was ist dise Filterung im Frequenzraum?
Eine Faltung. Warum? Weil die inverse FT genau der Faltungskern der 
Filterf'on |K| ist.

Was bringt das?
Spiele an der Filterfunktion und Du modifizierst, beispielsweise eine 
CT:
- TP: Grundstrukturen
- HP: Feinheiten
- BP: Gezielte Untersuchungen

Mit Sicherheit keine gute Erklärung - trotzdem - ich hoffe mein 
Geschreibsel bringt Dir trotzdem was.

Gruß,
Nils

Autor: Alexander Liebhold (lippi2000)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hey Danke Nils,
sehr ausführlich. Was ich brauche ist die Hintransformation. Die inverse 
Radon (wie bei CT notwendig) benötige ich im Moment nicht.

Also laß mich kurz zusammenfassen:

Ich mache eine 2D FFT/DFT und erhalte ein 2D Spektrum. Der Grundgedanke 
ist, dass eine Projektion mit Winkel 0° genau der Fouriertransformierten 
entspricht. (Zentralscheiben-Theorem)

Mal angenommen ich habe genau eine Linie parallel zur x-Achse. Dann 
sollte folgendes Spektrum entstehen:

Spektrum x-Richtung: Nur Zeile 0 enthällt Werte (DC Anteil)
Spektrum y-Richtung: Spaltfunktion da Linie einer Deltafunktion 
entspricht.

Beide Richtungen überlagert ergibt dann mein Ausgangsspektrum.

Nun führe ich die Projektion bei 0° durch und erhalte für Zeile 0 meine 
Funktionswerte (Zeile 0 enthält im Spektrum DC-Anteil).

Anschließend drehe ich meine Projektionseben. Neue Projektion 
durchführen und unter genau dem Winkel in mein Spektrum eintragen. Bei 
45° dürft ich als Projektion etwas annäherndes wie ein rechtwinkliges 
Dreieck erhalten.

Was ich nicht verstehe ist: Bei genau 90° habe ich als Projektion nur 
noch einen Peak. Im Spektrum müsste ich doch nun aber die Spaltfunktion 
sehen???


Bild im Ortsraum: ooo = Linie

########################
#          |           #
#          |           #
#          |  ooooooo  #
#          |           #
#          |           #
#----------------------#
#          |           #
#          |           #
#          |           #
#          |           #
#          |           #
########################

Gruß Alex

Autor: Nils (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo Alex,

> Spektrum y-Richtung: Spaltfunktion da Linie einer Deltafunktion entspricht.
Wenn mich meine Anschauung nicht veräppelt, würdest Du die Spaltfunktion 
erhalten, wenn Du für die Linie eine Linie mit endlicher Breite nimmst 
(was ja numerisch immer der Fall ist). Wenn Du das mit der 'echten' 
Delta-F'on rechnest, sollte eine Konstante rauskommen.
Liege ich da falsch?

Gruß,
Nils

Autor: Alexander Liebhold (lippi2000)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Absolut richtig!
Da die Deltafunktion unendlich kurz ist, besitzt sie Ausblendeigenschaft 
und hat ein unendliches Spektrum, also eine Linie da alle Frequenzen 
enthalten sind. Erst wenn ich ein Rechteck draus mache wird es zu einer 
Spaltfunktion. Kleiner Denkfehler. Danke.

Nach dem ich also die 2D FFT durchgeführt habe, habe ich alle 
Projektionen bestimmt. Diese findet man wieder, wenn man eine Ebene 
durch den Ursprung legt und diese um einen Winkel Phi dreht, dann kann 
man die Projektion aus dem Spektrum ablesen. Lieg ich richtig?

Wenn man die Koordinatenumwandlung weg lässt, müsste man eine inverse 1D 
FFT auf genau dieser gedrehten Ebene für alle Winkel Phi durchführen. 
Für diesen Schritt fehlt mir noch die Vorstellung. Denn wenn jede 
Projektion bereits auf einer um Phi gedrehten Ebene liegt, bräuchte ich 
diese doch nur in ein Diagramm r,Phi zu schreiben und hätte die 
Radontransformierte.

Oder sollten etwa die komplexen Werte in Polarkoordinaten aufgetragen 
werden? Und dann ist immer noch die inverse 1D FFT offen.

Autor: Alexander Liebhold (lippi2000)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Hier mal ein Beispiel.

Autor: Nils (Gast)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
> Wenn man die Koordinatenumwandlung weg lässt, müsste man eine inverse 1D
> FFT auf genau dieser gedrehten Ebene für alle Winkel Phi durchführen.
Da kann ich Dir nicht ganz folgen. Es ist doch so: Deine Messwerte, mit 
denen Du in die FFT rein gehst sind in Polarkoordinaten. Die FFT selbst 
setzt kartesische Koordinaten voraus. Die kleine (schlampige) Skizze 
deutet das an:
Die Fouriertransformierte der Projektion sind durch die roten kleinen 
Kreise gekennzeichnet. Die diskrete FT benötigt aber Werte im kartes. 
Raster (blau). In der Praxis muss daher interpoliert werden, bevor man 
in die FFT reingeht. Mit anderen Worten: die Interpolation findet im 
Koordinatenraum statt und nicht im Transformationsraum. Erst dann geht 
es in die FFT.
Kann es sein, das da das Problem liegt?

Die Sache mit der inversen 1D FFT folgt ja aus dem Projektionstheorem:
1-Dim FT der Projektion = 2-Dim FT
für einen gegebenen Projektionswinkel Phi.
Die Gesamtheit der Projektionen für alle Winkel Phi ist dann die 
Radon-Transformierte.
Ich überlege, ob man sich das vorstellen kann:
Wenn ich auf Fourierentwicklungen für bestimmte Winkel aus bin, kann ich 
natürlich das fragliche Gebiet f(x,y) für jeden Winkel entwickeln. Dabei 
integriert man das gesamte Gebiet für jede Koordinate und benötigt also 
eine 2-Dim FT. Das ist aber überflüssig, weil es reicht, die Projektion 
für den fraglichen Winkel zu betrachten. Ich benötige also nur eine FT 
entlang der Projektion. Das Wissen über das gesamte Gebiet entsteht aus 
allen Winkeln - der Radon-Transformierten.

Ich hoffe, ich habe jetzt nicht an Deinem Problem vorbeigeredet. 
Ansonsten schildere doch mal, was Du genau vor hast.

Autor: Alexander Liebhold (lippi2000)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Also ich meine das schon etwas anders.
Mein Ausgangsdatum ist ein Bild 2D und kartesisch.

Jetzt führe ich eine 2D-FFT durch. Im Transformationsraum bin ich immer 
noch kartesisch. Da ja eine 1D-FFT von einer Projektion mit Winkel phi 
gleich der 2D-FFT entlang einer Linie durch den Koordinatenursprung um 
den Winkel phi gedreht ist, folgt nun:

Die 2D-FFT habe ich bereits, nur entlang einer Linie ebend noch nicht. 
Dies geht einfach durch Wandlung in Polarkoordinaten. OK. Jetzt muss das 
noch mit der inversen 1D-FFT zurücktransformiert werden. Übrig bleibt 
die Projektion. Die 1D-FFT wird auf jeden Polarwinkel meiner 2D-FFT 
angewendet. Ich glaub das ist des Rätsels Lösung.

Interpolieren muss ich natürlich bei den Polarkoordinaten, da mit 
größerwerdendem Abstand natürlich immer mehr Werte für benachbarte 
Winkel entfallen. Bzw.gibt es Winkel, die nicht genau auf meinen 
gewünschten Projektionswinkelstücken liegen.

Fraglich ist jetzt nur wie man das am besten umsetzt???

Also 2DFFT auf ein Bild ist kein Problem. Einfach 2-mal ne 1D-FFT laufen 
lassen. Ein wenig schmalz muss man natürlich in die Polarkoordinaten und 
der Interpolation stecken. Ich würde dann einfach spaltenweise meine 
Polarkoordinaten anreihen.(Also jede Spalte für ein bestimmten Winkel 
0...360°) Da man nur eine Auswahl von Winkeln verwendet, müssen 
Ergebnisse mit nicht 100% zutreffenden Winkelangaben interpoliert 
werden. OK. Letzter Schritt ist dann eine spaltenweise 1D-IFFT um wieder 
aus dem Transformationsraum zurück zu kehren.

Ich glaub jetzt hab ich es verstanden.

Danke für die Diskussion.

Gruß Alexander

Autor: Detlef _a (detlef_a)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
>>Ich glaub jetzt hab ich es verstanden.

Ich noch nicht. Kannst/willst Du Deinen Code posten für die Bilder vom 
15. !?

THX
Cheers
Detlef

: Wiederhergestellt durch Admin
Autor: Alexander Liebhold (lippi2000)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Oh, es gibt auch Mitleser, wie schön.

Anbei das M-File. Das Eingangsbild muss ein Grauwertbild der Größe 
50x50x3 Pixel sein, da ich das Intensitätsbild aus Farbkanal 1 sofort 
übernehme. Ansonsten musst du den Code anpassen.

Wie gesagt, das Prinzip der Radontransformation ist sehr einfach. Kann 
man sich auch gut erklären. Ich habe für die Erklärungsversuche nicht 
die Projektionsebene gedreht, sondern das Bild an sich. Anschließend 
immer alle Zeilen von links nach rechts integriert (aufsummiert). Diese 
Spalte habe ich dann in ein Koordinatensystem gelegt f(r,phi). Und nun 
das Ursprungsbild um den Winkel phi drehen und immer von links nach 
rechts aufsummieren und diese Spalte an der entsprechenden Stelle im 
f(r,phi) Koordinatensystem ablegen. Anbei dazu noch ein kleines Bild.

Nur die Berechnung über FFT ist eher schwierig. Der eigentliche geistige 
Knacks ist das Fourier Slice Theorem.

Gruß Alexander

Autor: Harald M. (mare_crisium)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Alex,

Du wirst's nicht glauben, aber es gibt noch einen Mitleser ;-) !

Irgendwie erinnert mich die Radon-Transformierte an die 
Hough-Transformierte einer um 45 Grad geneigten Geraden. Ist das Zufall?

Gruss, Harald

Autor: Nils (Gast)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Hallo Alex,

ok, ich denke ich habe verstanden was Du vor hast:
Die 2Dim-FT bildet den Koordinatenraum für Deine Radon-Transformierte.
Du führst nun eine inverse FT entlang r für jeden Winkel Phi durch, um 
Dein Gebiet zu rekonstruieren.
Korrekt?

Dann müsstest Du auf einem komplexwertigem Gebiet interpolieren. Der 
einfachste Ansatz wäre eine lineare Interpolation. Bei den oben 
erwähnten CT-Anwendungen arbeitet man mit Polynomen oder Splines - aber 
da ist die Problemstellung ja anders (man hat Messungen aus einem 
rotierenden System und muss in kartes. Koord. um die FFT durchzuführen).

Vielleicht ist es ja eine Schnapsidee, aber man könnte so eine Art 
Mittelwertbildung machen. In der angefügten Skizze ist das 
2Dim-FT-Raster blau.
Es geht darum, einen Wert in dem Raster aus C1, C2, C3, C4 für P(r,Phi) 
zu interpolieren. Die C1 bis C4 sind komplexwertig aus der 2D-FT.
Dann könnte man die Werte mit Hilfe einer Abstandsfunktion D gewichten:
P(r,Phi) = (D(C1,P)*C1 + D(C2,P)*C2 + D(C3,P)*C3 + D(C4,P)*C4)/4
Dabei ist das D auf 1 normiert, d. h. D(C,P) = 1, falls Abstand minimal.
Natürlicher wäre es vielleicht mit einem quadratischen Abstands-Begriff 
zu spielen.
Ich bin mir nicht sicher, inwieweit das pratikabel für Dich ist.

Gruß,
Nils

Autor: Alexander Liebhold (lippi2000)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
@Nils: Du liegst absolut richtig. So meinte ich es.

Wenn ich deine Interpolation verstehe, machst du folgendes:

Berechnung des Schwerpunktes an der Stelle P. Für die entsprechende 
Wichtung verwendest du den Abstand zwischen dem gesuchten Punkt P und 
den jeweileigen C1-C4. Anschließend den komplexen Wert in 
Polarkoordinatensystem eintragen. Somit kann man sich selbst die 
Winkelaufteilung für die Berechnung wählen. Dürfte eigentlich so gehen.

Echt super. Wenn ich mal bissel mehr Zeit hätte würde ich versuchen das 
Ganze mal zu proggen. Hab zwar schon einige FFT-Algorithmen in C/Basic 
und ASM geschrieben aber en M-File solltes für die Interpolation sollte 
es erst mal tun. Möchte das ganze am Ende mal auf nem µC-Integrieren.

@Harald: Hey Mitleser :-) Also du hast das Erfasst. Bin über Hough auf 
die Radontransformation gekommen. Und musste dann feststellen, daß der 
Herr Hough sich vielleicht doch ganz schön mit fremden Blumen schmückt. 
Die Hough-Transformation ist lediglich eine simple Abwandlung der 
Radontr. Bzw.hat sich Hough einen Spezialfall zu nutzen gemacht.

Gruß Alexander

Autor: Nils (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
@Alex:
> Berechnung des Schwerpunktes an der Stelle P.
Das trifft es.
> Möchte das ganze am Ende mal auf nem µC-Integrieren.
Das meinte ich mit 'praktikabel'. Könnte auf einem µC zeitkritisch 
werden.
> Wenn ich mal bissel mehr Zeit hätte...
Geht mir genauso. Berichte doch weiter, wenn Du weiter gekommen bist. 
Hört sich sehr interessant an, was Du da machst.

Gruß,
Nils

Autor: Alexander Liebhold (lippi2000)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Moin, Moin,

also hab mir das wie folgt gedacht:

Winkelschritt festlegen: delta phi = 5°

Wertebereich der Y-Achse in Radonebene:

Jetzt müssen alle r für jeden Winkelschritt von phi durchlaufen werden.
delta_phi=5;
   r_max = sqrt(power(x-1,2)+power(y-1,2));
 for(phi=0;phi<360;phi+= delta_phi)
    {for(r=1;r<= r_max;r++)
       {x1=r*cos(phi);
        y1=r*sin(phi);
        //Interpolation
       }
    }

Auf die Punkte C1...C4 kommt man eigentlich auch ganz leicht.

C1(x,y) = [RoundDown(x1);RoundUp(y1)]
C2(x,y) = [RoundUp(x1);RoundUp(y1)]
C3(x,y) = [RoundDown(x1);RoundDown(y1)]
C4(x,y) = [RoundUp(x1);RoundDown(Y1)]

Nun die Abstände:

Nun deine Formel für die Gewichtsfunktion im Punkt P:

Alle Punkte C1-C4 und der Punkt P(x,y) sind natürlich komplex. Der 
Ergebniswert wird dann in eine Matrix POL[r,phi] eingetragen.

Jetzt noch die IFFT spaltenweise (also für alle phi).

Na mal schauen was Matlab zu der Überlegung sagt...

Autor: Alexander Liebhold (lippi2000)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Ach ja, in meinem M-File ist ein Fehler. Ich habe dort nicht kartesisch 
zu Polar transformiert sondern die komplexen Frequenzen nach Betrag und 
Phase ausgegeben. Dies ist natürlich falsch.

Autor: Alexander Liebhold (lippi2000)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Der Ursprung für die Polarkoordinaten muss ja nun eigentlich in der 
Mitte des Bildes(Frequenzraum) liegen???

FT=FFT2(image)
[xmax ymax] = size(FT)

Für die Polarkoordinatenberechnung des Punktes P(r,phi) als Beispiel sei 
r=2 und phi = 45° --> x=1.41 y=1.41 also liegt der Punkt zwischen 
C1(1,2);C2(2,2);C3(1,1);C4(2,1).

Die Wichtung kann auch berechnet werden. Der Ursprung liegt doch nun in 
der Bildmitte. Also ist P(2,45°)=(dC1*C1(xmax/2+xC1, 
ymax/2+yC1)..../(dC1+dC2+dC3+dC4)

Liege ich da richtig? Muss ich das komplette Spektrum verwenden? Ich 
denke mal ja. Somit wäre im Quadrant I das Spektrum für x,y. Im Quadrant 
IV dann der gespiegelte Spektralteil (gespiegelt an der x-Achse).

Autor: Harald M. (mare_crisium)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Alex,

Deine FFT sollte eigentlich die Frequenzen von -f_max bis +f_max 
umfassen, d.h. der Nullpunkt sollte in der Mitte der Realwert-Achse 
(x-Achse) liegen. Genauso müssten die Phasen, also der Imaginärteil, 
alle Werte zwischen -180° und +180° umfassen, d.h. der Nullpunkt sollte 
in der Mitte der Imaginärwert-Achse (y-Achse) liegen. In der angehängten 
pdf-Datei findest Du auch den Vorschlag für eine alternative 
Interpolationsmethode. Bin ein bisschen in Eile, deshalb alles nur in 
Kurzform; ich hoffe Du kannst es trotzdem lesen ;-).

Ciao

Harald.

Autor: Nils (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
> Der Ursprung für die Polarkoordinaten muss ja nun eigentlich in der
> Mitte des Bildes(Frequenzraum) liegen???
Ja, hätte ich so vorausgesetzt. Ansonsten kriegst Du einen 
Verschiebungsvektor bei der Polarumrechnung rein, der Dir Deine IFFT 
ziemlich erschweren würde. Das hat Harald ja beschrieben.

> Muss ich das komplette Spektrum verwenden?
Wenn Du es prinzipiell meinst: Für die Rekonstruktion des Gebiets auf 
jeden Fall. Die Information über das Gebiets ist ja durch die 2D-FT 
'integral formuliert'.
Wenn Du das in Richtung Optimierung meinst:
Da grüble ich auch drüber nach. Eigentlich reicht es ja die benachbarten 
Gitterpunkte entlang des Radiusvektors für jeden Winkel zu betrachten.
Dafür fällt mir aber kein schlauer Algorithmus ein (außer ganz brachial 
die interpolierten Werte in ein Array zu schreiben und dann der IFFT 
mundgerecht zu servieren. Aber das ist nicht wirklich toll). Aber in 
einem Schritt?

Dein P(x,y) ist richtig normiert. Wie gesagt, wenn das nicht gut 
funktioniert, könnte man auch ein quadratisches Mittel versuchen. 
Vielleicht liefert Haralds lineare Interpolation ja schon gute 
Ergebnisse - das müsste man probieren. Mein 'Gefühl' bei der Sache ist 
aber, dass ein gewichtetes quadratisches Mittel optimal wäre, weil dies 
(salopp gesagt) zum Abstandsbegriff der Fouriertransformation passt.

> Somit wäre im Quadrant I das Spektrum für x,y. Im Quadrant
> IV dann der gespiegelte Spektralteil (gespiegelt an der x-Achse).
Wieso? Gibt es Symmetrien, die ich übersehe?

Autor: Harald M. (mare_crisium)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Alex,

wie Nils schon sagt: Die gesamte verfügbare Information ist in den x*y 
Zahlen der 2D FFT vollständig enthalten - mehr gibt's nicht. Auch durch 
Interpolation kann man nicht mehr herausholen: Das ist genauso, als wenn 
Du ein Foto vergrösserst und dabei zwischen den Bildpunkten 
interpolierst. Es kommen keine zusätzlichen Details zum Vorschein, 
sondern das Bild wird immer flauer, weil sich dieselbe Informationsmenge 
auf eine grössere Bildfläche verteilt.

@Nils, ich hab' mal versucht, Deinen Vorschlag der quadratischen 
Gewichtung umzusetzen, weil er tatsächlich mit Ausdrücken hantiert, die 
der FFT angemessener sind. Die Schwierigkeit entsteht beim Normieren: 
Der grösste mögliche Abstand ist die Diagonale, darauf muss die Distanz 
normiert werden. Das hat aber zur Folge, dass z.B. im Punkt (x1,y1) mit 
dem Wert c1, die Gewichte der Punkte c2 bei (x1+1,y1) und c3 bei 
(x1,y+1) nicht gleich Null sind. Betrachte ich jetzt das Quadrat links 
daneben und berechne den interpolierten Wert im Punkt (x1-e,y1), dann 
bekomme ich Beiträge von c1 und c3 aber auch vom Punkt c1' bei 
(x1-1,y1). Deshalb gibt's auf dem Rand zwischen den benachbarten 
"Punktquadraten" bei starken Gradienten einen Sprung. Oder mach' ich 
einen Fehler?

Ciao,

Harald.

Autor: Nils (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
@Harald
Nur kurz von mir, da es ja schon spät/früh ist: In der Diskussion habe 
ich die Interpolation nur kurz angerissen und etwas nebulös auf die 
Umsetzung hingewiesen. Mein erster Gedanke dabei war, dass man die 
Ränder jeweils vorbesetzen muss, bzw. gesondert betrachtet. Im Grunde 
genommen werden die Gewichte auf dem Rand maximal und fallen zur Mitte 
eines Rasters hin ab - dass muss stetig passieren.
Da hast Du den Finger genau in die Wunde gelegt - ich denke Deine 
Überlegungen sind richtig.
Ich hoffe ich finde etwas Zeit das genauer zu überlegen - wahrscheinlich 
sind Alex und Du bei der Lösung wieder schneller.

Gruß,
Nils

Antwort schreiben

Die Angabe einer E-Mail-Adresse ist freiwillig. Wenn Sie automatisch per E-Mail über Antworten auf Ihren Beitrag informiert werden möchten, melden Sie sich bitte an.

Wichtige Regeln - erst lesen, dann posten!

  • Groß- und Kleinschreibung verwenden
  • Längeren Sourcecode nicht im Text einfügen, sondern als Dateianhang

Formatierung (mehr Informationen...)

  • [c]C-Code[/c]
  • [avrasm]AVR-Assembler-Code[/avrasm]
  • [code]Code in anderen Sprachen, ASCII-Zeichnungen[/code]
  • [math]Formel in LaTeX-Syntax[/math]
  • [[Titel]] - Link zu Artikel
  • Verweis auf anderen Beitrag einfügen: Rechtsklick auf Beitragstitel,
    "Adresse kopieren", und in den Text einfügen




Bild automatisch verkleinern, falls nötig
Bitte das JPG-Format nur für Fotos und Scans verwenden!
Zeichnungen und Screenshots im PNG- oder
GIF-Format hochladen. Siehe Bildformate.
Hinweis: der ursprüngliche Beitrag ist mehr als 6 Monate alt.
Bitte hier nur auf die ursprüngliche Frage antworten,
für neue Fragen einen neuen Beitrag erstellen.

Mit dem Abschicken bestätigst du, die Nutzungsbedingungen anzuerkennen.