Obrazowanie:Obrazowanie Medyczne/Podstawy Rekonstrukcji Obrazów Tomograficznych

Wstęp

Badany obiekt składa się z dwóch struktur (szare koła). Dokonajmy rzutowania obiektu na płaszczyznę. Prześwietlając obiekt w taki sposób, iż promienie Rentgenowskie padają na niego z lewej strony, otrzymamy rozkład natężenia wskazujący na obecność we wnętrzu tylko jednej struktury. Podobne rzutowanie, tylko dokonane z innego kierunku, wskazuje na obecność we wnętrzu obiektu dwóch struktur.

Poszczególne tkanki i narządy organizmu ludzkiego w różny sposób osłabiają promieniowanie X. Fakt ten został zaobserwowany jeszcze przez Wilhelma Roentgena, zaraz po odkryciu tego promieniowania i zastosowany do wykonywania zdjęć ręki. Trzy miesiące po odkryciu promieniowanie Rentgenowskiego wykorzystano je do prześwietlenia rannej głowy pacjenta. Po pierwszym roku stosowania promieniowania Rentgenowskiego w diagnostyce medycznej wiadomo było, że obraz wnętrze ciała ludzkiego powstaje na zasadzie rzucania przez poszczególne narządy, naświetlane promieniowaniem X, cieni na kliszę fotograficzną. Niestety, narządy organizmu ludzkiego wzajemnie się przesłaniają, co prowadzi od nakładania się na siebie obrazów poszczególnych struktur wewnętrznych człowieka. Ponadto, tworzenie zdjęć w klasycznej radiografii polega na rzutowaniu trójwymiarowego obiektu jakim jest człowiek na płaszczyznę, co dodatkowo zniekształca obraz. Wspominanie czynniki bardzo niekorzystnie wpłynęły na jakość zdjęć radiologicznych. Już w roku 1895 (rok po odkryciu promieniowania X), Elihu Thompson zauważył, iż wykonanie dwóch zdjęć radiologicznych z różnych pozycji lampy i detektora względem badanego obiektu, a następnie obejrzenie zdjęć w stroboskopie prowadzi do poprawienia jakości obrazu. Od tego momentu było wiadomo, iż wprowadzenie ruchu lampy Rentgenowskiej i detektora względem badanego obiektu pozwala otrzymać informację, umożliwiającą poprawę jakość zdjęć ciała ludzkiego. Zagadnienie to zilustrowano na rys. 1. Od tego momentu rozpoczęły się prowadzone na szeroką skalę badania nad problemem ruchu lampy rentgenowskiej i detektora względem pacjenta. Ich efektem było ciekawych rozwiązań, funkcjonujących jeszcze po zakończeniu drugiej wojny światowej. W międzyczasie, w roku 1917 austriacki matematyk Johann Radon, udowodnił pewne twierdzenie, które dzisiaj leży u podstaw tzw. metod rekonstrukcji obrazów. Praktyczna realizacja osiągnięcia Radona jest w zasadzie możliwa tylko przy pomocy komputerów, dlatego na początku dwudziestego wieku twierdzenie to przeszło niezauważone. Twierdzenie tym zajmiemy się na początku kolejnego rozdziału. W latach 50-tych i 60-tych ubiegłego wieku, problem wykonywania przestrzennego obrazu ciała ludzkiego nadal nie był rozwiązany. Funkcjonowało wiele technik wykonywania tego rodzaju zdjęć, jednakże jakość uzyskiwanych obrazów nie była zadowalająca. Każdy z odkrywców nadawał opracowanej przez siebie metodzie własną nazwę, i tak do roku 1962 roku można było spotkać się z takimi terminami jak stratigraphy, planigraphy, body-section radiography, tomography, verigraphy, radiothomy. Wszystkie one dotyczyły rozwiązań, w których lampa rentgenowska i detektory wykonywały ruch dookoła pacjenta. Ostatecznie, w roku 1962 International Commision on Radioation Units and Measurements, zarekomendowała dla tych metod termin ‘’’Tomografia’’’ (gr. thomos – warstwa, graphein – rysować), czyli obrazowanie warstwowe. Od roku 1956 Allan McLeod Cormack, amerykański fizyk pochodzenia południowoafrykańskiego w ramach hobby interesuje się promieniowaniem X i dokonuje wyliczeń będących podwalinami współczesnej tomografii rentgenowskiej. Wyniki publikuje w Journal of Applied Physics w 1963 i 1964 roku. Prace te jednak nie zostają zauważone. W tych samych latach Godfrey Newbold Hounsfield, brytyjski inżynier elektronik pracuje nad urządzeniami tomograficznymi, jednakże bez znaczącego sukcesu. W końcu napotyka na obliczenia Cormack’a i w oparciu o nie buduje pierwszy Rentgenowski Tomograf Komputerowy (X – Ray Computed Tomography, CT, w żargonie lekarskim zwany po prostu tomografem), zaprezentowany w roku 1971 i nazwany EMI. CT znajduje się w użytku od roku 1972 (pierwszy zbadany pacjent). W roku 1979 Hounsfield i Cormack otrzymują nagrodę Nobla, z dziedziny fizjologii i medycyny za opracowanie Komputerowej Tomografii Rentgenowskiej (Hounsfield ponadto otrzymuje od Królowej Brytyjskiej tytuł szlachecki Sir).

Twierdzenie Radona i Transformata Radona

W roku 1905 W. Radon udowodnił następujące twierdzenie: „Obraz obiektu dwuwymiarowego można zrekonstruować na podstawie nieskończone ilości rzutów jednowymiarowych”. Jak okaże się w dalszych częściach poniższego rozdziału, rzutowanie to odpowiada wykonywaniu na obiekcie pewnej transformacji, nazywanej Transformacją Radona. Dokonanie na wynikach rzutowania Odwrotnej Transformacji Radona umożliwia zrekonstruowanie obrazu obiektu. W następnym podrozdziale podamy matematyczne sformułowanie tego twierdzenia, natomiast w tym miejscu pragniemy zwrócić uwagę na pewien fakt. W twierdzeniu pojawia się termin rekonstrukcja. Otóż, w przeciwieństwie do klasycznej radiografii, w której zdjęcie obiektu powstaje natychmiast na filmie, potem zaś obraz jest tylko wywoływany i utrwalany, w przypadku tomografii, nie można na podstawie serii pomiarów od razu uzyskać obrazu. Niezbędny jest dodatkowy etap – wyznaczenie na podstawie zebranych danych obrazu badanego obiektu. Proces ten nazywamy rekonstrukcją.

Podstawy Matematyczne — Delta Diraca

Poprawne zrozumienie Twierdzenia Radona wymaga opanowania pewnych zagadnień matematycznych, które przypomnimy w poniższym podrozdziale. właściwości:
A.

B.

C.

D.

E.

F.

Inne ciekawe właściwości wynikające z właściwości A-F:

G.

Dowód:
H.

Dowód:
[math] \delta(ax+b)=\delta(a(x+\frac{b}{a}))={1}{|a|}\delta(x+\frac{b}{a}) [/math]

I.

Dowód:

[math] \begin{array}{l} \int_{-\infty}^{\infty} g(x)\delta(ax+b)dx=\int_{-\infty}^{\infty} g(x)\delta(a(x+\frac{b}{a}))dx \\ \\ \int_{-\infty}^{\infty} g(x)\delta(a(x+\frac{b}{a}))dx=\frac{1}{|a|}\int_{-\infty}^{\infty} g(x)\delta(x+\frac{b}{a}))dx \\ \\ x+\frac{b}{a}=x'\\ \\ x=x'-\frac{b}{a} \\ \\ dx = dx' \\ \\ \int_{-\infty}^{\infty} g(x)\delta(x+\frac{b}{a})dx = \int_{-\infty}^{\infty} g(x'-\frac{b}{a})\delta(x')dx' \\ \\ \int_{-\infty}^{\infty} g(x'-\frac{b}{a})\delta(x')dx' =\frac{1}{|a|}g(-\frac{b}{a}) \\ \end{array} [/math]

Transformata Radona w kartezjańskim układzie współrzędnych

Ilustracja prezentująca wyznaczanie transformaty Radona. Obiekt opisany jest za pomocą dwuwymiarowej funkcji g(x,y), która reprezentuje rozkład jego pewnej cechy fizycznej (np. gęstości, rozkładu radiofarmaceutyka, liniowego współczynnika osłabienia promieniowania X). Dokonujemy sumowania wartości funkcji g(x,y) wzdłuż pewnej prostej, parametryzowanej za pomocą współczynników a i b. Wynik sumowania — G(a,b) zależy oczywiście oczywiście od przebiegu prostej i jest prezentowany we współrzędnych A-B parametrów prostej.

Transformata Radona to transformata całkowa o następującej postaci: Transformatę Radona zilustrowano na rys. 2. Własności:
A. Transformata Radona jest liniowa.
Dowód:
Niech będzie dana funkcja [math]h(x,y)[/math], będąca kombinacją liniową funkcji [math]g_i(x,y)[/math]: Wtedy jej transformata Radona [math]h(x,y)[/math] wynosi: i jest kombinacją liniową transformat Radona funkcji [math]g_i(x,y)[/math].

B. Transformata Radona funkcji przesuniętej do współrzędnych [math](x_0,y_0)[/math]
Dokonując zamiany zmiennych:

dostajemy:

suma wyrażeń: [math]kx_0+b-y_0[/math] jest pewną stałą, będącą nowym współczynnikiem przesunięcia prostej względem początku układu współrzędnych, wzdłuż której następuje całkowanie funkcji Ostatecznie dostajemy następujący wynik. Jeżeli transformata Radona funkcji [math]g(x,y)[/math] wynosi [math]G(a,b)[/math], to transformata Radona funkcji [math]g(x-x_0,y-y_0)[/math] jest równa:

Przykład 1. Transformata Radona punktu znajdującego się w środku układu współrzędnych.

Ilustracja do Przykładu I — wyznaczenie Transformaty Radona Punktu znajdującego się w środku układu współrzędnych.

Zakładamy, że punkt ma nieskończenie wielką amplitudę. Taki punkt możemy zamodelować za pomocą Delty Diraca: co oznacza, że funkcja [math]g(x,y)[/math] ma następującą postać: Obliczmy transformatę Radona powyższej funkcji: Przekształcimy ten wzór, korzystając z właściwości Delty Diraca.

Ostatecznie otrzymaliśmy zatem następujący wynik:

Jak widzimy, transformata Radona punktu, znajdującego się w środku układu współrzędnych da wynik niezerowy, jeśli współczynnik przesunięcia prostej będzie równy zero. Jest to wynik zgodny z naszą intuicją, bowiem, tylko prosta o takim współczynniku przesunięcia przechodzi przez środek układu współrzędnych, a zatem pokrywa się z rozważanym punktem. Innymi słowy, transformata Radona przekształca punkt leżący w początku układu współrzędnych, na prostą w przestrzeni parametrów Radona biegnącą równolegle do osi pionowej i przechodzącej przez punkt 0. Wynik ten zilustrowano na rys. 3.

Przykład 2. Transfromata Radona punktu o dowolnych współrzędnych [math](x_0,y_0)[/math]

Ilustracja do Przykładu II — wyznaczenie Transformaty Radona punktu o dowolnych współrzędnych.

Jest to bardzo ważny przykład, ponieważ każdy obiekt możemy przedstawić jako sumę punktów o różnych współrzędnych, których amplituda reprezentuje pewną cechę fizyczną obiektu. Wykorzystując dalej liniowość Transformaty Radona, możemy obliczyć Transformatę Radona danego obiektu jako sumę Transformat Radona punktów składających się na ten obiekt. Podobnie jak w poprzednim przykładzie, punkt możemy wyrazić za pomocą Delty Diraca. Postać funkcji [math]g(x,y)[/math] jest zatem następująca: Transformatę Radona powyższej funkcji obliczymy wykorzystując właściwości B transformaty Radona. Transformatę Radona punktu leżącego w środku układu współrzędnych wynosi: [math] g(x,y) = \delta(x)(y)\xrightarrow{\textrm{Transfromata Radona}}G(a,b)=\delta(b) [/math] a zatem transformata Radona punktu przesuniętego jest równa: Transformata Radona przekształca zatem punkt o dowolnych współrzędnych w prostą w przestrzeni parametrów Radona. Wynik ten zilustrowano na rys. 4

Przykład 3. Transformata Radona prostej.

Ilustracja do Przykładu III — wyznaczenie Transformaty Radona linii prostej.

Przykład ten podajemy ze względu na jego znaczenie w grafice komputerowej. Uprzedzając wyniki obliczeń podamy, iż Transformata Radona punktu w kartezjańskim układzie współrzędnych prowadzi do punktu we układzie współrzędnych Transformaty Radona. Punkt te ma współrzędne odpowiadające parametrom prostej. W ten sposób można wykorzystać Transformatę Radona do wynajdywanie w obiekcie, czy obrazach linii prostych.
Prostą [math]y=kx+l[/math] możemy opisać przy pomocy Delty Diraca, zakładając że punkty leżące na prostej mają nieskończoną amplitudę. wtedy: jej transformata Radona jest równa: całkując prawą stronę powyższego równania najpierw po zmiennej dy oraz wykorzystując własność Delty Diraca dostajemy: W zależności od wartości parametrów [math]a,b,k,l[/math] możemy rozważyć następujące przypadki rozwiązań: a. [math]k\neq a[/math] oraz dowolne [math]l[/math] i [math]b[/math] Korzystając z własności B Delty Diraca ostatecznie otrzymujemy: b. [math]k=a[/math] i [math]l\neq b[/math], wtedy: c. [math]k=a[/math] i [math]l=b[/math]

Podsumowując, transformata Radona funkcji prostej przekształca ją w przestrzeni (A-B) w funkcję płaską o amplitudzie [math]\frac{1}{|a-k|}[/math], która osiąga wartość nieskończoną w punkcie o współrzędnych odpowiadających parametrom prostej. Wynik ten zilustrowano na rys. 5

Normalna postać Transformaty Radona

Obiekty spotykane w rzeczywistości cechuje zwykle pewna symetria, w związku z czym dogodniej jest je opisywać w innych niż kartezjańskim układzie współrzędnych. Np. człowiek w przybliżeniu posiada symetrię walcową. Ponadto Twierdzenie Radona mówi o tym, iż do wiernego zrekonstruowania obiektu potrzebna jest nieskończona liczba rzutów tego obiektu. Warunek ten w praktyce jest oczywiście niemożliwy do zrealizowania, w związku z czym zawsze będziemy mieli do czynienia z obrazem będącym pewnym przybliżeniem rzeczywistego obiektu. Niemniej powstaje pytanie, z jakich kierunków dokonać rzutowania badanego obiektu. Okazuje się, że dowolne ruchy występujące w przyrodzie, można uzyskać za pomocą obrotu i przesunięcia. Uwzględniając symetrię ciała człowieka, wspomniane ruchy wygodniej jest opisywać w biegunowym układzie współrzędnych. Przejdziemy teraz zatem do wyrażenia Transformaty Radona w biegunowym układzie współrzędnych, którą będziemy nazywać Normalną Transformatą Radona. Pomocne przy tym będą dwa dodatkowe zagadnienia — parametryzacja proste jw układzie biegunowym oraz obrót układu współrzędnych, które opiszemy w następnych rozdziałach.

Parametryzacja prostej w układzie biegunowym .

Ilustracja do zagadnienia parametryzacji prostej.

Na początku dokonajmy parametryzacji prostej — zamiast opisywać prostą za pomocą współczynnika a (nachylenia prostej względem układu OX) oraz współczynnika b — przecięcia prostej z osią OY wprowadzimy parametry [math]t_p[/math], określający odległość prostej od początku układu współrzędnych oraz kąt φ jaki tworzy normalna do prostej względem osi OX. Znaczenie parametrów [math]t_p[/math] i φ zilustrowano na rys. 3. Przy pomocy parametrów [math]t_p[/math] i φ równanie prostej może być zapisane w następujący sposób:

Obrót układu współrzędnych:

Ilustracja do zagadnienia obrotu układu współrzędnych.

Przypomnijmy znaną z kursu matematyki transformację współrzędnych punktu [math]x_0,y_0[/math] w układzie X-Y do układu T-S, obróconego względem niego o kąt [math]\phi[/math]. W układzie T-S punkt ma współrzędne [math]t_0,s_0[/math]:

Zagadnienie to zilustrowano na rys. 7. Jakobian tego przekształcenia wynosi:

Normalna Transformata Radona

Parametryzacja prostej za pomocą przesunięcia [math]t_p[/math] oraz kąta φ umożliwia następujący zapis Transformaty Radona:

Dokonamy teraz pewnych przekształceń, polegających na:

• wykorzystaniu transformacji układu współrzędnych pomiędzy układem X-Y i układem obróconym T-S,
• wykorzystaniu pewnych właściwości Delty Dirac'a.

Otrzymany wzór to alternatywna postać Normalnej Transformaty Radona. Prosze zauważyć, iż wyrażenie

Przykład IV — Normalna transformata Radona punktu

Rozważmy ponownie ważne zagadnienie, jakim jest Transformata Radona Punktu. Tak jak w przykładzie II punkt o współrzędnych [math]x_0,y_0[/math] opiszemy za pomocą Delty Dirac'a.

Jak wygląda wykres otrzymanej funkcji [math]G(t_0,\phi)[/math]. Wygodniej będzie zapisać współrzędne punktu [math]x_0, y_0[/math] we współrzędnych biegunowych, co zilustrowano na rys. 8.

Punkt w biegunowym układzie współrzędnych.

Wynikiem przekształcenia punktu opisanego w biegunowym układzie współrzędnych za pomocą Transformaty Radona jest krzywa sinusoidalna.

Otrzymaliśmy zatem następujący wynik:

Wyrażenie [math]\delta(t_0\sin(\phi + \phi_0) - t_p)[/math] jest niezerowe, tylko wtedy, gdy: czyli gdy punkty [math]t_p[/math] i [math]\phi[/math] leżą na sinusoidzie. Wynik ten zilustrowano na zilustrowano na rys. 9.

Sinogram

Po lewej stronie rysunku badany obiekt, składający się z dwóch kwadratów. Po prawej stronie rysunku jego sinogram, czyli Transformata Radona.

Jak już to było wspomniane w poprzednich rozdziałach, dowolny obiekt możemy zapisać jako sumę punktów o różnych współrzędnych. W związku z tym, iż Transformatą Radona punktu w biegunowym układzie współrzędnych jest krzywa sinusoidalna, Transformata Radona dowolnego obiektu będzie składała się z bardzo wielu sinusoid. Jest to powód, dla którego wykres Transformaty Radona nazywamy sinogramem. Przykład obiektu i jego sinogramu zaprezentowano na rys. 10.

Odwrotna Transformata Radona.

Poprzednie rozdziały dotyczyły tylko i wyłącznie Transformaty Radona oraz niektórych jej właściwości. Jaki związek ma Transformata Radona z rekonstrukcją obrazu badanego obiektu? Przypominamy, iż Transformata Radona jest pewnym rodzajem jedno-wymiarowego rzutowania. Obiekt badamy wzdłuż pewnej prostej, sparametryzowanej np. w biegunowym układzie współrzędnych. Wynikiem pojedynczego pomiary wzdłuż określonego przez prostą kierunku,jest suma punktów tworzących obiekt i leżących na tej prostej. Punkty te reprezentują pewną wybraną cechę obiektu. Okazuje się, że istnieje transformata odwrotna do Transformaty Radona, która na podstawie serii opisanych pomiarów umożliwia odtworzenie obrazu obiektu. Jeśli transformata obiektu [math]g(x,y)[/math] ma postać: to jej transformata odwrotna wynosi

Podana postać odwrotnej Transformaty Radona jest niewygodna do implementacji numerycznej (wymaga czasochłonnego liczenia wielu pochodnych), niemniej jest dowodem na to, iż możliwa jest na podstawie jednowymiarowych rzutów rekonstrukcja obrazu dwu lub trój wymiarowego obiektu. Pozostaje jeszcze kwestia, jak aparatura medyczna jest w stanie dokonać wyznaczenia Transformaty Radona obiektu. Tematem tym zajmiemy się w kolejnych rozdziałach.