Edukacyjna rekonstrukcja, czyli rekonstrukcja obrazu za pomocą frameworku DAPHNE. Przykład DFR (Direct Fourier Reconstruction) i CST (Central Section Theorem)

marca 04, 2023

Edukacyjna rekonstrukcja, czyli rekonstrukcja obrazu za pomocą frameworku DAPHNE. Przykład DFR (Direct Fourier Reconstruction) i CST (Central Section Theorem)

Dzisiejszy post będzie postem edukacyjnym. Co prawda w moim doktoracie używam tylko pakietów typowo matematycznych i do wyświetlania obrazu, ale istnieje wiele ścieżek nauki rekonstrukcji obrazu (ASTRA toolbox, MATLAB etc.). Jednym z nich jest właśnie framework DAPHNE.

DAPHNE (Didactics tomogrAPHic recoNstruction framEwork) to framework służący do przeprowadzenia prostych eksperymentów z najbardziej popularnymi algorytmami rekonstrukcji dla PET i CT.

W czym będziemy pracować?

Będziemy pracować w Jupyter. Jeśli chcesz skorzystać z wersji online zapraszam Cię tu: https://jupyter.org/. Wybieramy opcję Try it on your browser. Wybieramy opcję JupyterLab i Python.

Możesz też zainstalować środowisko Anaconda, gdzie Jupyter znajdzie się już w tej instalacji: https://www.anaconda.com/ lub miniconda: https://repo.anaconda.com/miniconda/. Zaznajom się koniecznie z Jupyter przez realizacją tego ćwiczenia. Szybkie wykonanie kodu w komórce jest możliwe za pomocą skrótu klawiszowego Shift+ENTER.

Serializacja i deserializacja danych

Python dostarcza moduł, który pozwala użytkownikowi na serializację i deserializację obiektów: pickle. W DAPHNE można uzyskać dostęp do serializacji pickle poprzez narzędzia Pickle i Unpickle, które są zawarte w pakiecie Misc. Obiekty serializowane za pomocą pickle nie są czytelne dla człowieka ani nie mają standardowego formatu, który można by otworzyć za pomocą zewnętrznych programów.

Rekonstrukcja za pomocą wiązki równoległej

Algorytmy dostępne w ramach DAPHNE opierają się na twierdzeniu centralnego przekroju , łączącym transformatę Fouriera (FT) obiektu z jego RT. Z tego powodu metody te znane są również jako metody rekonstrukcji obrazu oparte na transformacie Fouriera. Najprostszym algorytmem opartym na CST jest tzw. bezpośrednia rekonstrukcja fourierowska (Direct Fourier Reconstruction - DFR), choć obecnie najczęściej stosowaną w praktyce metodą jest filtrowana projekcja wsteczna (Filtered Backprojection - FBP).

Istotną kwestią będzie przyswojenie teorii na temat rekonstrukcji za pomocą wiązki równoległej. W omawianym przykładzie jako dane wejściowe przyjmiemy wstępnie obliczony zestaw całek liniowych dwuwymiarowego fantomu głowy Shepp-Logana. Wykonamy następujące czynności:

Import biblioteki.
Definicja parametrów geometrycznych
Wczytanie danych fantomu i wstępnie obliczonych projekcji w geometrii wiązki równoległej
Ustawienie algorytmu DFR i wykonanie rekonstrukcji
Wizualizacja danych
Zaimportujmy najpierw główne (w istocie bardzo nieliczne) biblioteki wymagane dla tego prostego projektu.

Import bibliotek

Biblioteki, z których będziemy korzystać to m. in. numpy, matplotlib i scipy. Zaczynamy od ich importu. Jeśli nie posiadasz tych bibliotek zainstalowanych, to zrób to przez Anaconda Navigator lub https://anaconda.org/.

import numpy as np
import matplotlibimport matplotlib.pyplot as pltplt.rcParams.update({'font.size': 10})

import scipy.interpolateimport scipy.miscimport scipy.ndimage.interpolation
%matplotlib notebook

Wczytanie i wyświetlenie danych fantomu i sinogramu

Teraz możemy wczytać dane z dysku dla fantomu. Niektóre parametry geometryczne, takie jak liczba binów radialnych i kątowych, a także liczba pikseli obrazu docelowego, są definiowane i inicjalizowane na początku kolejnej komórki. Dane można pobrać stąd: DFdemo.

N_ang = 180           # N. of angles
N_rad = 256           # N. of radial samples
Imsize_X = 256
Imsize_Y = 256

phantom = np.fromfile("../Data/SL_HC_256x256.raw",dtype=np.float32).reshape((Imsize_Y,Imsize_X))
sino = np.fromfile("../Data/SL_HC_180x256_paralsino_theta_r.raw",dtype=np.float32).reshape((N_rad,N_ang))
sino = np.rot90(sino,2)

# Zero-padding. See function reference at https://numpy.org/doc/stable/reference/generated/numpy.pad.html
padded_sino = np.pad(sino, [(0,N_rad),(0,0)], mode='constant')

print("Shape of the phantom object: " + str(phantom.shape))
print("Shape of the sinogram object: " + str(sino.shape))
print("Shape of the padded sinogram object: " + str(padded_sino.shape))

W efekcie otrzymujemy następujące dane:

Shape of the phantom object: (256, 256)
Shape of the sinogram object: (256, 180)
Shape of the padded sinogram object: (512, 180)

Wyświetlamy fantom i dane projekcyjne.

plt.figure(figsize=(9,4))
#plt.figure()
plt.gray()
plt.subplot(121)
plt.title("Phantom")
plt.xlabel(r"$x$ (pixel)")
plt.ylabel(r"$y$ (pixel)")
plt.imshow(phantom)

plt.subplot(122)
plt.title("Sinogram")
plt.xlabel(r"$\theta_{id}$")
plt.ylabel(r"$x^{\prime}$ (pixel)")
plt.imshow(sino)
plt.tight_layout()

Po uruchomieniu tego kodu otrzymujemy:



Twierdzenie o przekroju centralnym (Central Section Theorem)

Teraz wykorzystamy CST, aby wykonać rekonstrukcję DFR dla pojedynczego wycinka. W tym celu wykorzystamy FFT (Fast Fourier Transform), które znajduje się w bibliotece NumPy. Pierwszym krokiem jest oczywiście obliczenie jednowymiarowej FFT sinogramu linia po linii wzdłuż kierunku radialnego. Użycie funkcji fftshift (https://numpy.org/doc/stable/reference/generated/numpy.fft.fftshift.html) jest wymagane, aby umieścić ujemne częstotliwości w odpowiednim miejscu podczas obliczania FFT.

sino_fft1d=np.fft.fftshift(
                np.fft.fft(
                    np.fft.ifftshift(
                        sino,
                        axes=0
                    ),
                axis=0
                ),
            axes=0
            )

Zgodnie z CST zawartość informacyjna stosu 1D FFT linii radialnych sinogramu jest taka sama jak FFT 2D obiektu wzdłuż linii radialnych. Wizualnie zweryfikujemy to później w kolejnych komórkach. Pozostało nam wyświetlić część rzeczywistą i urojoną sinogramu FFT.

plt.figure(figsize=(9,5))
plt.subplot(121)
plt.title("Sinogram FFT (real)")
plt.imshow(np.real(sino_fft1d))
plt.xlabel(r"$\theta_{id}$")
plt.ylabel(r"$\nu$ (pixel$^{-1}$)")

plt.subplot(122)
plt.title("Sinogram FFT (imag)")
plt.imshow(np.imag(sino_fft1d))
plt.xlabel(r"$\theta_{id}$")
plt.ylabel(r"$\nu$ (pixel$^{-1}$)")

Co nam daje:

CST jest faktycznie ekspolatowany, aby dokonać rekonstrukcji. Utworzymy pomocnicze tablice zawierające odpowiednie współrzędne w dziedzinie częstotliwości. Dokładniej, th i nu będą przechowywały współrzędne kątowe i radialne 𝜈,𝜃
sinogramu poddanego transformacji Fouriera. Dwie tablice będą natomiast przechowywać współrzędne kartezjańskie dwuwymiarowej płaszczyzny Fouriera, oznaczane przez 𝜉,𝜐.

th=np.array([np.pi*i/N_ang for i in range(N_ang)])
nu=np.arange(N_rad)-N_rad/2
th,nu=np.meshgrid(th,nu)
nu=nu.flatten()
th=th.flatten()

# Coordinates of 2D frequency space
xi_pol=(N_rad/2)+nu*np.cos(th)
upsilon_pol=(N_rad/2)+nu*np.sin(th)

Porównamy teraz próbki transformaty 1D FFT sinogramu, ułożonego w siatce radialnej, z FFT 2D samego obiektu. CST mówi, że te dwa byty powinny być równe w przestrzeni ciągłej, w okręgu dziedziny częstotliwości takiej, że 𝜉2+𝜐2≤𝜈2𝑚𝑎𝑥
  (gdzie 𝜈𝑚𝑎𝑥
  jest częstotliwością Nyquista).

# Display the 1D radial FFT of the sinogram 
plt.gray()
plt.figure(figsize=(9,3))
plt.subplot(131)
plt.title(r"Sinogram FFT along $x^{\prime}$")
plt.xlabel(r"$\theta$ (rad)", labelpad=0)
plt.ylabel(r"$\nu$", labelpad=-10)
plt.imshow(np.log(np.abs(sino_fft1d)), extent=[0,np.pi,-0.5,0.5], aspect="auto")

# Display the remapping of the 1D FFT of the sinogram on 2D frequency space
plt.subplot(132)
plt.title("Remapped sino FFT in 2D freq. domain")
plt.scatter(
    (xi_pol-Imsize_X/2)/Imsize_X,
    (upsilon_pol-Imsize_Y/2)/Imsize_Y,
        c=np.log(np.abs(sino_fft1d.flatten())),
    marker='.',
    edgecolor='none'
    )
plt.xlabel(r"$\xi$")
plt.ylabel(r"$\upsilon$", labelpad=-8)

# Display the 2D FFT of phantom for comparison
plt.subplot(133)
plt.title("2D FFT of the phantom")
plt.imshow(np.log(np.abs(phantom_fft2d)), extent=[-0.5,0.5,-0.5,0.5], aspect="auto")
plt.xlabel(r"$\xi$")
plt.ylabel(r"$\upsilon$", labelpad=-8)

plt.tight_layout()

Otrzymujemy następujące wyniki:


Siatkowanie danych w dziedzinie częstotliwości i rekonstrukcja obrazu
Figura w prawej ramce jest FFT 2D rekonstruowanego obiektu, to również figura w środkowej ramce powinna być taka sama przynajmniej według (przestrzeni nieciągłej) CST. Rzeczywiście, wyświetlony obraz powyżej w środkowej ramce jest tylko wykresem rozproszenia, więc nadal musimy przemapować go na siatkę kartezjańską przed wysłaniem go do odwrotnej 2D FFT i próbą rekonstrukcji. W tym celu użyjemy funkcji griddata biblioteki scipy.interpolate. Więcej na temat samego działania interpolate zapraszam do przeczytania dokumentacji Pythona https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html.

# Create a rectangular grid of points in 2D frequency domain
xi_cart,upsilon_cart=np.meshgrid(np.arange(N_rad),np.arange(N_rad))
xi_cart=xi_cart.flatten()
upsilon_cart=upsilon_cart.flatten()

# Interpolate the 2D Fourier space grid from the transformed sinogram
remapped_fft2d=scipy.interpolate.griddata(
        (xi_pol,upsilon_pol),
        sino_fft1d.flatten(),
        (xi_cart,upsilon_cart),
        method='cubic',
        fill_value=0.0
    ).reshape((N_rad,N_rad))

Teraz próbki sinogramu poddanego transformacji 1D-FFT zostały przemapowane do siatki kartezjańskiej przy użyciu interpolacji sześciennej. Przekształcona tablica jest przechowywana w remapped_fft2d. Ponieważ tablica ta przechowuje próbki FFT w siatce kartezjańskiej, jesteśmy teraz gotowi do ich rekonstrukcji przy użyciu odwrotnej FFT 2D.

# The final step - reconstruct the image
# Inverse transform from 2D FFT to the image space
recon=np.real(
    np.fft.fftshift(
        np.fft.ifft2(
            np.fft.ifftshift(remapped_fft2d)
            )
        )
    )

Wyświetlamy oryginalny obraz i obraz zrekonstruowany oraz błąd, który osiągnęliśmy w tej rekonstrukcji.

# Compute the error image (original - recon)
difference=np.subtract(phantom,recon)

# Display the reconstructed image and difference with original
plt.figure(figsize=(9,5))
plt.subplot(131)
plt.title("Reconstruction (DFR) \n min\max: " + str(np.around(np.min(recon),3)) + "\\" + str(np.around(np.max(recon),3)))
plt.imshow(recon, interpolation='bicubic')
plt.xlabel(r"$x$ (pixel)")
plt.ylabel(r"$y$ (pixel)", labelpad = 0)
plt.subplot(132)
plt.title("Original \n min\max: " + str(np.around(np.min(phantom),3)) + "\\" + str(np.around(np.max(phantom),3)))
plt.imshow(phantom, interpolation='bicubic')
plt.xlabel(r"$x$ (pixel)")
plt.ylabel(r"$y$ (pixel)", labelpad = 0)
plt.subplot(133)
plt.title("Difference  \n min\max: " + str(np.around(np.min(difference),3)) + "\\" + str(np.around(np.max(difference),3)))
plt.imshow(difference, interpolation='bicubic')
plt.xlabel(r"$x$ (pixel)")
plt.ylabel(r"$y$ (pixel)", labelpad = 0)
plt.tight_layout()


 
I to tyle, jeśli chodzi o rekonstrukcję tą metodą. W następnym wpisie zrobimy rekonstrukcję FBP z wykorzystaniem bibliotek. Pełen kod znajduje się na GitHubie: https://github.com/jolapodolszanska/DFR-reconstruction.

Bibliografia
Panetta D, Camarlinghi N. 3D Image Reconstruction for CT and PET : A Practical Guide with Python. CRC Press; 2020. Available from: https://www.taylorfrancis.com/books/9780429270239

Szukaj na tym blogu

Medycznie w IT

Edukacyjna rekonstrukcja, czyli rekonstrukcja obrazu za pomocą frameworku DAPHNE. Przykład DFR (Direct Fourier Reconstruction) i CST (Central Section Theorem)

W czym będziemy pracować?

Serializacja i deserializacja danych

Rekonstrukcja za pomocą wiązki równoległej

Import bibliotek

Wczytanie i wyświetlenie danych fantomu i sinogramu

Twierdzenie o przekroju centralnym (Central Section Theorem)

Siatkowanie danych w dziedzinie częstotliwości i rekonstrukcja obrazu

Bibliografia

Komentarze

Prześlij komentarz

Popularne posty

Wearables – jak urządzenia noszone zmieniają nasze zdrowie

Jak odczytywać tagi w pliku DICOM?