Co jest złego w tym kodzie do rekonstrukcji tomograficznej metodą Fouriera?

Ostatnio bawiłem się algorytmami rekonstrukcji tomograficznej. Mam już ładne działające implementacje FBP, ART, iteracyjny schemat SIRT / SART, a nawet używam prostej algebry liniowej (powoli!). To pytanie nie dotyczy żadnej z tych technik ; odpowiedzi w formie „dlaczego ktoś miałby to zrobić w ten sposób, zamiast tego podajemy kod FBP”, nie są tym, czego szukam.

Następną rzeczą, którą chciałem zrobić z tym programem, było „ uzupełnienie zestawu ” i wdrożenie tak zwanej „ metody rekonstrukcji Fouriera ”. Rozumiem to w zasadzie, że zastosujesz FFT 1D do sinogramów „ekspozycji”, ustawisz je jako promieniowe „szprychy koła” w przestrzeni Fouriera 2D (że jest to użyteczna rzecz do zrobienia bezpośrednio z twierdzenia o wycinku środkowym) , interpoluj te punkty do regularnej siatki w tej przestrzeni 2D, a następnie powinna istnieć możliwość odwrotnej transformacji Fouriera w celu odzyskania pierwotnego celu skanowania.

Brzmi prosto, ale nie miałem dużo szczęścia, kiedy otrzymałem jakieś rekonstrukcje, które wyglądają jak oryginalne cele.

Poniższy kod Python (numpy / SciPy / Matplotlib) dotyczy najbardziej zwięzłego wyrażenia, jakie mogłem wymyślić z tego, co próbuję zrobić. Po uruchomieniu wyświetla następujące informacje:

Rysunek 1: cel

Rycina 2: sinogram celu

Rycina 3: Wiersze sinogramu edytowane FFT

Rysunek 4: górny rząd to przestrzeń 2D FFT interpolowana z rzędów sinogramów w dziedzinie Fouriera; dolny rząd to (dla celów porównawczych) bezpośrednia 2D FFT celu. W tym momencie zaczynam się podejrzewać; wykresy interpolowane z sinogramu FFT wyglądają podobnie do wykresów wykonanych bezpośrednio 2D-FFTingiem celu ... a jednak różnią się.

Ryc. 5: odwrotna transformata Fouriera z ryc. 4. Mam nadzieję, że będzie to bardziej rozpoznawalne jako cel niż w rzeczywistości.

Jakieś pomysły, co robię źle? Nie jestem pewien, czy moje rozumienie rekonstrukcji metody Fouriera jest zasadniczo wadliwe, czy jest tylko jakiś błąd w moim kodzie.

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, wreszcie to złamałem.

Sztuczka polegała w zasadzie na umieszczeniu niektórych fftshift/ ifftshifts we właściwym miejscu, aby reprezentacja przestrzeni 2D Fouriera nie była gwałtownie oscylacyjna i skazana na niemożność dokładnej interpolacji. Tak przynajmniej myślę to naprawiło. Większość tego, co mam ograniczone rozumienie teorii Fouriera, opiera się na ciągłym całkowaniu i zawsze uważam tę dyskretną domenę i FFT za nieco ... dziwaczną.

Chociaż uważam, że kod Matlab jest raczej tajemniczy, muszę przyznać, że ta implementacja przynajmniej dała mi pewność, że ten algorytm rekonstrukcji może być wyrażony dość zwięźle w tego rodzaju środowisku.

Najpierw pokażę wyniki, a następnie kod:

Rysunek 1: nowy, bardziej złożony cel.

Rysunek 2: sinogram (OK OK, to transformacja Radona) celu.

Rycina 3: Wiersze sinogramu z edycją FFT (wykreślone z DC w środku).

Rycina 4: Sinogram FFT przekształcony w przestrzeń 2D FFT (DC w środku). Kolor jest funkcją wartości bezwzględnej.

Rysunek 4a: Powiększ środek 2D przestrzeni FFT, aby lepiej pokazać promieniowy charakter danych sinogramu.

Rysunek 5: Górny rząd: przestrzeń 2D FFT interpolowana z rozmieszczonych promieniowo rzędów sinogramów FFT. Dolny rząd: oczekiwany wygląd po prostu 2D FFT-ing celu.

Ryc. 5a: Powiększ środkowy obszar wykresów podrzędnych na ryc. 5, aby pokazać, że wyglądają one na całkiem dobrą zgodność jakościową.

Rysunek 6: Test kwasowy: odwrotna FFT 2D interpolowanej przestrzeni FFT odzyskuje cel. Lena nadal wygląda całkiem nieźle pomimo wszystkiego, przez co właśnie ją przeszliśmy (prawdopodobnie dlatego, że jest wystarczająco dużo „szprych” sinogramu, aby dość gęsto pokryć płaszczyznę 2D FFT; rzeczy stają się interesujące, jeśli zmniejszysz liczbę kątów ekspozycji, więc nie jest to już prawdą ).

Oto kod; przedstawia wątki w czasie krótszym niż 15 sekund na 64-bitowym SciPy Debiana / Wheezy na 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()

Aktualizacja 2013-02-17: Jeśli jesteś wystarczająco zainteresowany, aby przejść przez to miejsce, więcej wyników z programu do nauki, którego był częścią, można znaleźć w postaci tego plakatu . Treść kodu w tym repozytorium może również być interesująca (chociaż należy pamiętać, że kod nie jest tak usprawniony jak powyższy). W pewnym momencie mogę spróbować spakować go jako „notebook” IPython.

Nie wiem dokładnie, gdzie jest problem, ale twierdzenie o przekroju oznacza, że ​​te dwa specjalne przypadki powinny być prawdziwe:

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

Postępuj zgodnie z kodem i spróbuj znaleźć punkt, w którym przestają one być równoważne, pracując do przodu od sinogramu i wstecz od wygenerowanego 2D FFT.

To nie wygląda dobrze:

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]])

Generowany przez Ciebie 2D FFT jest obrócony o 90 stopni względem tego, co powinien być?

Sugeruję pracę z wielkością i fazą, a nie z rzeczywistością i wyobrażeniami, aby łatwiej zobaczyć, co się dzieje:

(Białe rogi to nie wszystko log(abs(0)), nie stanowią problemu)

Wierzę, że rzeczywista teoretyczny powód, dla którego pierwsze rozwiązanie nie działa pochodzi z faktu, że obroty są wykonywane w odniesieniu do centrów obrazach, powodując przesunięcie [S/2, S/2], co oznacza, że każdy z wierszy waszych sinogramnie pochodzi od 0do S, ale raczej od -S/2do S/2. W twoim przykładzie przesunięcie jest w rzeczywistości offset = np.floor(S/2.). Zauważ, że to działa na Sparzyste lub nieparzyste, a jest odpowiednikiem tego, co zrobiłeś w kodzie S/2(chociaż jest bardziej wyraźne problemy unika, kiedy Sto float, na przykład).

Domyślam się, że opóźnienia fazowe, które wprowadza ta zmiana w transformacie Fouriera (FT), są źródłem tego, o czym mówisz w drugim komunikacie: fazy są pomieszane i trzeba to zrekompensować, aby móc zastosuj odwrócenie transformaty Radona. Trzeba pogłębić tę teorię, aby mieć pewność, czego dokładnie potrzeba, aby odwrotność działała zgodnie z oczekiwaniami.

Aby zrekompensować to przesunięcie, możesz użyć fftshift tak, jak to robiłeś (co umieszcza środek każdego wiersza na początku, a ponieważ użycie DFT faktycznie odpowiada obliczeniu transformacji Fouriera sygnału okresowego S, otrzymujesz odpowiednie rzeczy ) lub jawnie kompensuj ten efekt w złożonej transformacie Fouriera podczas obliczania sinogramFT. W praktyce zamiast:

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

możesz usunąć ifftshifti pomnożyć każdy wiersz przez wektor korekcyjny:

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)

Wynika to z właściwości transformacji Fouriera, jeśli chodzi o przesunięcie czasowe (sprawdź na stronie wikipedii FT „twierdzenie o przesunięciu” i zastosuj przesunięcie równe - offset- ponieważ umieszczamy obraz z powrotem wokół środka).

Podobnie, możesz zastosować tę samą strategię do rekonstrukcji i zastąpić fftshiftją korektą faz, w obu wymiarach, ale w innym kierunku (kompensacja wstecz):

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))
        )
    )

Cóż, to nie poprawia twojego rozwiązania, ale rzuca kolejne światło na teoretyczne aspekty twojego pytania. Mam nadzieję, że to pomaga!

Ponadto nie lubię go używać, fftshiftponieważ ma on tendencję do zakłócania sposobu fftobliczania. W takim przypadku należy jednak umieścić środek FT na środku obrazu przed interpolacją, aby uzyskać fft2(lub przynajmniej zachować ostrożność podczas ustawiania r- aby można było zrobić to całkowicie za fftshiftdarmo!), I fftshiftnaprawdę się przydaje tam. Wolę jednak nadal używać tej funkcji do celów wizualizacji, a nie w ramach „rdzenia” obliczeniowego. :-)

Z poważaniem,

Jean-Louis

PS: próbowałeś zrekonstruować obraz bez przycinania koła? co daje całkiem fajny efekt rozmycia na rogach, fajnie byłoby mieć taką funkcję w programach takich jak Instagram, prawda?