Was stimmt mit diesem Code für die tomographische Rekonstruktion nach der Fourier-Methode nicht?

Ich habe in letzter Zeit mit tomographischen Rekonstruktionsalgorithmen herumgespielt. Ich habe bereits gute Implementierungen von FBP, ART, einem SIRT / SART-ähnlichen iterativen Schema und sogar der Verwendung von linearer Algebra (langsam!). Bei dieser Frage geht es nicht um eine dieser Techniken . Antworten in der Form "Warum sollte es jemand so machen, hier ist ein FBP-Code" sind nicht das, wonach ich suche.

Das nächste, was ich mit diesem Programm machen wollte, war " Komplettieren " und Implementieren der sogenannten " Fourier-Rekonstruktionsmethode ". Mein Verständnis ist, dass Sie eine 1D-FFT auf das Sinogramm "Exposures" anwenden, diese als radiale "Speichen eines Rades" im 2D-Fourier-Raum anordnen (dies ist eine nützliche Sache, die sich direkt aus dem zentralen Slice-Theorem ergibt). interpolieren Sie von diesen Punkten zu einem regelmäßigen Gitter in diesem 2D-Raum, und dann sollte es möglich sein, die Fourier-Transformation umzukehren, um das ursprüngliche Scan-Ziel wiederherzustellen.

Klingt einfach, aber ich habe nicht viel Glück gehabt, Rekonstruktionen zu bekommen, die in etwa wie das ursprüngliche Ziel aussehen.

Der folgende Python-Code (numpy / SciPy / Matplotlib) handelt von dem prägnantesten Ausdruck, den ich mir ausdenken konnte. Beim Ausführen wird Folgendes angezeigt:

Abbildung 1: das Ziel

Abbildung 2: Ein Sinogramm des Targets

Abbildung 3: Die FFT-ed-Sinogrammzeilen

Fig. 4: Die oberste Zeile ist der 2D-FFT-Raum, der aus den Fourier-Domain-Sinogrammzeilen interpoliert wird. Die untere Reihe ist (zu Vergleichszwecken) die direkte 2D-FFT des Ziels. Dies ist der Punkt, an dem ich anfange, misstrauisch zu werden; Die aus den Sinogramm-FFTs interpolierten Diagramme ähneln den durch direktes 2D-FFTing des Ziels erstellten Diagrammen ... und sind dennoch unterschiedlich.

Abbildung 5: Die inverse Fourier-Transformation von Abbildung 4. Ich hätte gehofft, dass dies als Ziel etwas besser erkennbar ist als es tatsächlich ist.

Irgendwelche Ideen, was ich falsch mache? Ich bin mir nicht sicher, ob mein Verständnis der Fourier-Methodenrekonstruktion grundlegend fehlerhaft ist oder ob mein Code nur einen Fehler enthält.

import math
import matplotlib
import matplotlib.pyplot as plt
import numpy as np

import scipy.interpolate
import scipy.fftpack
import scipy.ndimage.interpolation

S=256  # Size of target, and resolution of Fourier space
A=359  # Number of sinogram exposures

# Construct a simple test target
target=np.zeros((S,S))
target[S/3:2*S/3,S/3:2*S/3]=0.5
target[120:136,100:116]=1.0

plt.figure()
plt.title("Target")
plt.imshow(target)

# Project the sinogram
sinogram=np.array([
        np.sum(
            scipy.ndimage.interpolation.rotate(
                target,a,order=1,reshape=False,mode='constant',cval=0.0
                )
            ,axis=1
            ) for a in xrange(A)
        ])

plt.figure()
plt.title("Sinogram")
plt.imshow(sinogram)

# Fourier transform the rows of the sinogram
sinogram_fft_rows=scipy.fftpack.fftshift(
    scipy.fftpack.fft(sinogram),
    axes=1
    )

plt.figure()
plt.subplot(121)
plt.title("Sinogram rows FFT (real)")
plt.imshow(np.real(np.real(sinogram_fft_rows)),vmin=-50,vmax=50)
plt.subplot(122)
plt.title("Sinogram rows FFT (imag)")
plt.imshow(np.real(np.imag(sinogram_fft_rows)),vmin=-50,vmax=50)

# Coordinates of sinogram FFT-ed rows' samples in 2D FFT space
a=(2.0*math.pi/A)*np.arange(A)
r=np.arange(S)-S/2
r,a=np.meshgrid(r,a)
r=r.flatten()
a=a.flatten()
srcx=(S/2)+r*np.cos(a)
srcy=(S/2)+r*np.sin(a)

# Coordinates of regular grid in 2D FFT space
dstx,dsty=np.meshgrid(np.arange(S),np.arange(S))
dstx=dstx.flatten()
dsty=dsty.flatten()

# Let the central slice theorem work its magic!
# Interpolate the 2D Fourier space grid from the transformed sinogram rows
fft2_real=scipy.interpolate.griddata(
    (srcy,srcx),
    np.real(sinogram_fft_rows).flatten(),
    (dsty,dstx),
    method='cubic',
    fill_value=0.0
    ).reshape((S,S))
fft2_imag=scipy.interpolate.griddata(
    (srcy,srcx),
    np.imag(sinogram_fft_rows).flatten(),
    (dsty,dstx),
    method='cubic',
    fill_value=0.0
    ).reshape((S,S))

plt.figure()
plt.suptitle("FFT2 space")
plt.subplot(221)
plt.title("Recon (real)")
plt.imshow(fft2_real,vmin=-10,vmax=10)
plt.subplot(222)
plt.title("Recon (imag)")
plt.imshow(fft2_imag,vmin=-10,vmax=10)

# Show 2D FFT of target, just for comparison
expected_fft2=scipy.fftpack.fftshift(scipy.fftpack.fft2(target))

plt.subplot(223)
plt.title("Expected (real)")
plt.imshow(np.real(expected_fft2),vmin=-10,vmax=10)
plt.subplot(224)
plt.title("Expected (imag)")
plt.imshow(np.imag(expected_fft2),vmin=-10,vmax=10)

# Transform from 2D Fourier space back to a reconstruction of the target
fft2=scipy.fftpack.ifftshift(fft2_real+1.0j*fft2_imag)
recon=np.real(scipy.fftpack.ifft2(fft2))

plt.figure()
plt.title("Reconstruction")
plt.imshow(recon,vmin=0.0,vmax=1.0)

plt.show()

OK, ich habe das endlich geknackt.

Der Trick bestand im Grunde darin, einige fftshift/ ifftshifts an der richtigen Stelle zu platzieren, damit die 2D-Fourier-Raumdarstellung nicht wild schwankt und es unmöglich ist, genau zu interpolieren. Zumindest glaube ich, dass dies behoben wurde. Der größte Teil meines begrenzten Verständnisses für die Fouriertheorie basiert auf der kontinuierlichen Integralformulierung, und ich finde die diskrete Domäne und die FFTs immer etwas ... skurril.

Obwohl ich Matlab-Code als eher kryptisch empfinde, muss ich dieser Implementierung das Vertrauen zuschreiben, dass dieser Rekonstruktionsalgorithmus in einer solchen Umgebung relativ kompakt ausgedrückt werden kann.

Zuerst zeige ich Ergebnisse, dann Code:

Abbildung 1: Ein neues, komplexeres Ziel.

Abbildung 2: Das Sinogramm (OK OK, es ist die Radon-Transformation) des Ziels.

Abbildung 3: Die FFT-ed-Zeilen des Sinogramms (mit DC in der Mitte aufgetragen).

Abbildung 4: Das in den 2D-FFT-Raum transformierte FFT-Ed-Sinogramm (DC in der Mitte). Farbe ist eine Funktion des absoluten Wertes.

Abbildung 4a: Vergrößern Sie die Mitte des 2D-FFT-Raums, um die Radialität der Sinogrammdaten besser zu veranschaulichen.

Abbildung 5: Obere Zeile: Der aus den radial angeordneten FFT-ed-Sinogrammzeilen interpolierte 2D-FFT-Raum. Untere Reihe: Das erwartete Erscheinungsbild, wenn einfach eine 2D-FFT auf das Ziel angewendet wird.

Abbildung 5a: Zoomen Sie in den mittleren Bereich der Unterzeichnungen in Abbildung 5, um zu zeigen, dass diese qualitativ ziemlich gut übereinstimmen.

Abbildung 6: Der Säuretest: Die inverse 2D-FFT des interpolierten FFT-Raums stellt das Ziel wieder her. Trotz allem, was wir gerade durchgemacht haben, sieht Lena immer noch ziemlich gut aus (wahrscheinlich, weil es genügend Sinogramm- "Speichen" gibt, um die 2D-FFT-Ebene ziemlich dicht abzudecken. Interessant wird es, wenn Sie die Anzahl der Belichtungswinkel verringern, sodass dies nicht mehr zutrifft ).

Hier ist der Code; zeigt die Handlungen in weniger als 15 Sekunden auf Debian / Wheezys 64-Bit-SciPy auf einem i7.

import math
import matplotlib
import matplotlib.pyplot as plt
import numpy as np

import scipy.interpolate
import scipy.fftpack
import scipy.misc
import scipy.ndimage.interpolation

S=256 # Size of target, and resolution of Fourier space
N=259 # Number of sinogram exposures (odd number avoids redundant direct opposites)

V=100 # Range on fft plots

# Convenience function
def sqr(x): return x*x

# Return the angle of the i-th (of 0-to-N-1) sinogram exposure in radians.
def angle(i): return (math.pi*i)/N

# Prepare a target image
x,y=np.meshgrid(np.arange(S)-S/2,np.arange(S)-S/2)
mask=(sqr(x)+sqr(y)<=sqr(S/2-10))
target=np.where(
    mask,
    scipy.misc.imresize(
        scipy.misc.lena(),
        (S,S),
        interp='cubic'
        ),
    np.zeros((S,S))
    )/255.0

plt.figure()
plt.title("Target")
plt.imshow(target)
plt.gray()

# Project the sinogram (ie calculate Radon transform)
sinogram=np.array([
        np.sum(
            scipy.ndimage.interpolation.rotate(
                target,
                np.rad2deg(angle(i)), # NB rotate takes degrees argument
                order=3,
                reshape=False,
                mode='constant',
                cval=0.0
                )
            ,axis=0
            ) for i in xrange(N)
        ])

plt.figure()
plt.title("Sinogram")
plt.imshow(sinogram)
plt.jet()

# Fourier transform the rows of the sinogram, move the DC component to the row's centre
sinogram_fft_rows=scipy.fftpack.fftshift(
    scipy.fftpack.fft(
        scipy.fftpack.ifftshift(
            sinogram,
            axes=1
            )
        ),
    axes=1
    )

plt.figure()
plt.subplot(121)
plt.title("Sinogram rows FFT (real)")
plt.imshow(np.real(sinogram_fft_rows),vmin=-V,vmax=V)
plt.subplot(122)
plt.title("Sinogram rows FFT (imag)")
plt.imshow(np.imag(sinogram_fft_rows),vmin=-V,vmax=V)

# Coordinates of sinogram FFT-ed rows' samples in 2D FFT space
a=np.array([angle(i) for i in xrange(N)])
r=np.arange(S)-S/2
r,a=np.meshgrid(r,a)
r=r.flatten()
a=a.flatten()
srcx=(S/2)+r*np.cos(a)
srcy=(S/2)+r*np.sin(a)

# Coordinates of regular grid in 2D FFT space
dstx,dsty=np.meshgrid(np.arange(S),np.arange(S))
dstx=dstx.flatten()
dsty=dsty.flatten()

plt.figure()
plt.title("Sinogram samples in 2D FFT (abs)")
plt.scatter(
    srcx,
    srcy,
    c=np.absolute(sinogram_fft_rows.flatten()),
    marker='.',
    edgecolor='none',
    vmin=-V,
    vmax=V
    )

# Let the central slice theorem work its magic!
# Interpolate the 2D Fourier space grid from the transformed sinogram rows
fft2=scipy.interpolate.griddata(
    (srcy,srcx),
    sinogram_fft_rows.flatten(),
    (dsty,dstx),
    method='cubic',
    fill_value=0.0
    ).reshape((S,S))

plt.figure()
plt.suptitle("FFT2 space")
plt.subplot(221)
plt.title("Recon (real)")
plt.imshow(np.real(fft2),vmin=-V,vmax=V)
plt.subplot(222)
plt.title("Recon (imag)")
plt.imshow(np.imag(fft2),vmin=-V,vmax=V)

# Show 2D FFT of target, just for comparison
expected_fft2=scipy.fftpack.fftshift(
    scipy.fftpack.fft2(
        scipy.fftpack.ifftshift(
            target
            )
        )
    )

plt.subplot(223)
plt.title("Expected (real)")
plt.imshow(np.real(expected_fft2),vmin=-V,vmax=V)
plt.subplot(224)
plt.title("Expected (imag)")
plt.imshow(np.imag(expected_fft2),vmin=-V,vmax=V)

# Transform from 2D Fourier space back to a reconstruction of the target
recon=np.real(
    scipy.fftpack.fftshift(
        scipy.fftpack.ifft2(
            scipy.fftpack.ifftshift(fft2)
            )
        )
    )

plt.figure()
plt.title("Reconstruction")
plt.imshow(recon,vmin=0.0,vmax=1.0)
plt.gray()

plt.show()

Update 2013-02-17: Wenn Sie interessiert genug waren, um durch dieses Los zu waten, finden Sie in Form dieses Posters weitere Ergebnisse aus dem Selbststudienprogramm, zu dem es gehört . Der Hauptteil des Codes in diesem Repository kann ebenfalls von Interesse sein (obwohl der Code nicht annähernd so rationalisiert ist wie der oben genannte). Ich kann versuchen, es irgendwann als IPython "Notebook" neu zu verpacken.

Ich weiß nicht genau, wo das Problem liegt, aber der Slice-Satz besagt, dass diese beiden Sonderfälle zutreffen sollten:

fft2(target)[0] = fft(sinogram[270])
fft2(target)[:,0] = fft(sinogram[0])

Befolgen Sie also Ihren Code und versuchen Sie, den Punkt zu finden, an dem diese nicht mehr äquivalent sind, indem Sie vom Sinogramm aus vorwärts und von der generierten 2D-FFT aus rückwärts arbeiten.

Das sieht nicht richtig aus:

In [47]: angle(expected_fft2[127:130,127:130])
Out[47]: 
array([[-0.07101021,  3.11754929,  0.02299738],
       [ 3.09818784,  0.        , -3.09818784],
       [-0.02299738, -3.11754929,  0.07101021]])

In [48]: fft2_ = fft2_real+1.0j*fft2_imag

In [49]: angle(fft2_[127:130,127:130])
Out[49]: 
array([[ 3.13164353, -3.11056554,  3.11906449],
       [ 3.11754929,  0.        , -3.11754929],
       [ 3.11519503,  3.11056604, -2.61816765]])

Die 2D-FFT, die Sie generieren, ist um 90 Grad von der Soll-Position gedreht.

Ich würde vorschlagen, eher mit Betrag und Phase als mit real und imaginär zu arbeiten, damit Sie leichter sehen können, was passiert:

(Die weißen Ecken sind -inf davon log(abs(0)), sie sind kein Problem)

Ich glaube, dass der eigentliche theoretische Grund, warum die erste Lösung nicht funktioniert hat, in der Tatsache liegt, dass die Rotationen in Bezug auf die Bildmitten ausgeführt werden und einen Versatz von induzieren [S/2, S/2], was bedeutet, dass jede Ihrer Reihen sinogramnicht von 0bis ist S, sondern von -S/2bis S/2. In Ihrem Beispiel ist der Versatz tatsächlich offset = np.floor(S/2.). Beachten Sie, dass dies bei Sgeraden oder ungeraden Werten funktioniert und dem entspricht, was Sie in Ihrem Code getan haben S/2(wenngleich dies expliziter ist, werden Probleme vermieden, wenn Ses sich floatbeispielsweise um ein a handelt).

Ich vermute, dass die Phasenverzögerungen, die diese Verschiebung in der Fouriertransformation (FT) verursacht, die Ursache für das sind, worüber Sie in Ihrer zweiten Nachricht sprechen: Die Phasen sind durcheinander, und man muss diese Verschiebung kompensieren, um dies zu können Wenden Sie die Inversion der Radon-Transformation an. Man muss sich eingehender mit dieser Theorie befassen, um sicherzugehen, was genau erforderlich ist, damit das Inverse wie erwartet funktioniert.

Um diesen Versatz zu kompensieren, können Sie entweder fftshift wie zuvor verwenden (wodurch die Mitte jeder Zeile am Anfang steht. Da die Verwendung von DFT tatsächlich der Berechnung der Fourier-Transformation eines S-periodischen Signals entspricht, erhalten Sie das richtige Material ), oder kompensieren Sie diesen Effekt explizit in der komplexen Fourier-Transformation, wenn Sie die sinogramFT berechnen . In der Praxis anstelle von:

sinogram_fft_rows=scipy.fftpack.fftshift(
    scipy.fftpack.fft(
        scipy.fftpack.ifftshift(
            sinogram,
            axes=1
            )
        ),
    axes=1
    )

Sie können ifftshiftjede Zeile entfernen und mit einem Korrekturvektor multiplizieren:

offset = np.floor(S/2.)
sinogram_fft_rows = scipy.fftpack.fftshift(
    scipy.fftpack.fft(sinogram, axis=1)
    * (np.exp(1j * 2.* np.pi * np.arange(S) * offset / S)),
    axes=1)

Dies ergibt sich aus den Fouriertransformationseigenschaften, wenn eine Zeitverschiebung in Betracht gezogen wird (überprüfen Sie die FT-Wikipedia-Seite für das "Verschiebungssatz" und wenden Sie die Verschiebung gleich an, - offsetda wir das Bild wieder in die Mitte setzen).

Ebenso können Sie die gleiche Strategie auf die Rekonstruktion anwenden und die fftshiftdurch die Korrektur der Phasen in beiden Dimensionen, aber in der anderen Richtung (Rückkompensation) ersetzen :

recon=np.real(
    scipy.fftpack.ifft2(
        scipy.fftpack.ifftshift(fft2)
        *  np.outer(np.exp(- 1j * 2.* np.pi * np.arange(S) * offset / S),
                    np.exp(- 1j * 2.* np.pi * np.arange(S) * offset / S))
        )
    )

Nun, dies verbessert Ihre Lösung nicht, sondern wirft ein anderes Licht auf die theoretischen Aspekte Ihrer Frage. Hoffentlich hilft das!

Darüber hinaus bin ich nicht so gern mit, fftshiftweil es dazu neigt, mit der Art und Weise, wie die fftberechnet wird, herumzuspielen. In diesem Fall müssen Sie jedoch die Mitte der FT vor der Interpolation in die Bildmitte setzen, um zu erhalten fft2(oder zumindest vorsichtig bei der Einstellung zu sein r- damit Sie sie vollständig frei machen können fftshift!), Und dies ist fftshiftwirklich praktisch Dort. Ich bevorzuge es jedoch, diese Funktion für Visualisierungszwecke zu verwenden und nicht innerhalb des Berechnungs- "Kerns". :-)

Freundliche Grüße,

Jean-Louis

PS: Haben Sie versucht, das Bild zu rekonstruieren, ohne den Kreis zu beschneiden? Das ergibt einen ziemlich coolen Unschärfeeffekt in den Ecken. Wäre es nicht schön, eine solche Funktion in Programmen wie Instagram zu haben?