10 Podstawowa matematyka rekonstrukcji tomograficznych.pdf

(301 KB) Pobierz
Microsoft Word - 10_mat_rek.doc
X. PODSTAWOWA MATEMATYKA REKONSTRUKCJI
TOMOGRAFICZNYCH
10.1 Definicje; metoda wstecznej projekcji w tomografii transmisyjnej
Zakładamy, że mamy wiązkę równoległą o rozmiarach w x h , gdzie w – szerokość, a h -
wysokość wiązki. Rzecz będzie więc dotyczyć tomografów pierwszej generacji. Przyjmijmy,
że dokonaliśmy pomiaru wzdłuż przerywanej linii, równoległej do osi y r obracającego się
układu (x r ,y r ), związanego z układem źródło-detektor, podczas gdy z nieruchomym pacjentem
związany jest układ (X,Y). W danej odległości x r od osi y r ustawionej pod kątem Φ
w stosunku do osi Y pacjenta mierzone natężenie wynosi:
+∞
I
(
Φ
,
x
)
=
I
exp
μ
(
x
,
y
)
dy
,
r
0
r
r
r
Y
y r
x r
(10.1)
Φ
(x,y)
gdzie relacja pomiędzy współrzędnymi
punktu (x,y) oraz (x r ,y r ) jest następująca:
Θ
X
t
x
r
=
x
cos
Φ
+
y
sin
Φ
(10.2)
y
r
=
x
sin
Φ
+
y
cos
Φ
a μ (x r ,y r ) oznacza liniowy współczynnik
Rys. 10.1 Przyjęty do opisu układ
współrzędnych
pochłaniania związany z punktem (x r , y r )
= (x,y). Wzór (10.1) łatwo jest
przekształcić na:
I
+∞
p
(
Φ
,
x
)
=
ln
0
=
μ
(
x
,
y
)
dy
(10.3)
r
I
(
Φ
,
x
)
r
r
r
r
1
73222147.001.png
Wielkość p(Φ,t) nosi nazwę transformaty Radona wielkości μ. Wielkość tę nazywamy dla
prostoty nazywali projekcją . Łatwo stwierdzić, że wystarczy ją zmierzyć tylko na półokręgu,
gdyż zamiana detektora i źródła miejscami nie może zmieniać wartości projekcji, tj.
p
Φ
(
x
r
)
p
Φ
,
x
r
)
=
p
π
+
Φ
,
x
r
)
(10.4)
Zadaniem tomografii jest, jak już mówiliśmy we wcześniejszym wykładzie, zrekonstruowanie
funkcji μ (x r ,y r ), a więc i μ (x,y). Rekonstrukcja ta wcale nie jest prosta, nie tylko ze względów
czysto matematycznych. Przede wszystkim trzeba mieć świadomość, że obrazując rozkład
współczynnika pochłaniania w danej płaszczyźnie zakładamy, że dane pomiarowe
rzeczywiście dotyczą nieskończenie cienkich przekrojów, tak że zamiast trójwymiarowych
voxeli możemy mówić o dwuwymiarowych pixelach. Po drugie, zakładamy, że wszystkie
rejestrowane fotony poruszały się po liniach prostych pomiędzy źródłem a detektorem.
W rzeczywistości wiązka promieniowania X ma skończone rozmiary i rozbieżność kątową,
a w trakcie przechodzenia przez obiekt wiązka ulega „stwardnieniu”, gdyż promieniowanie
o niższych energiach jest silniej pochłaniane i do odleglejszych warstw przechodzi relatywnie
więcej promieniowania o wyższej energii. Dodatkowo jeszcze zakładamy, że rozkład
współczynnika absorpcji w ramach rozmiaru wiązki i jej rozbieżności kątowej jest
jednorodny.
Charakterystyczną odmiennością metody SPECT od transmisyjnej tomografii komuterowej
(CT) jest badanie nie tyle współczynnika pochłaniania w danym obszarze, ile aktywności
wychodzącej z danego miejsca, tak więc mierzymy wielkość
+∞
p
A
(
Φ
,
x
r
)
=
A
(
x
r
,
y
r
)
dy
r
,
(10.5)
gdzie A(x r ,y r ) – aktywność wychodząca z punktu (x r ,y r ), którą wyznaczamy w oparciu
o pomiar wielkości p A ( Φ ,x r ). Zarówno w metodzie CT, jak i SPECT dążymy do wyznaczenia
rozkładu dwuwymiarowego z serii mierzonych danych jednowymiarowych.
Niech pomiary będą wykonywane w serii kroków, w których kąt Φ zmienia się o δΦ, a
odległość x r zmienia o δx r = δt. Efektywnie wyznaczamy zatem wielkości p ij , gdzie
2
p ij
=
p
(
δ
,
j
δ
t
)
(10.6)
Nas interesuje odpowiednia funkcja podcałkowa we wzorze (10.3) lub (10.5). Funkcję tę
także rekonstruujemy na dyskretnej siatce pixeli o rozmiarach np. w x w , a więc zmierzamy
do znalezienia wartości
μ =
ij
μ
( jw
iw
,
)
(10.7)
lub
A ij =
A
( jw
iw
,
)
(10.8)
W praktyce numerycznej wygodnie jest posługiwać się raczej macierzą jednowymiarową niż
dwuwymiarową. Jeśli rozmiar interesującej nas macierzy wynosi N = n x n , to można zapisać
pierwsze n współczynników pierwszego wiersza, następnie przypisać pierwszemu
współczynnikowi drugiego wiersza element (n+1)- szy, pierwszemu współczynnikowi
trzeciego wiersza element (2n+1)- szy itd. W podobny sposób możemy opisać zarówno
wielkości mierzone, jak i rekonstruowane. W takiej zdyskretyzowanej formie nasze równania
mają postać:
=
N
p
j
=
r
jk
μ
k
(10.9)
k
1
Jeśli dysponujemy J pomiarami projekcji p j , to wielkości {μ k } teoretycznie można otrzymać
przez proste odwrócenie macierzy r jk . W praktyce nie jest to wcale takie proste. Po pierwsze,
liczba elementów tej macierzy jest znaczna. Jeśli n = 256, to liczba elementów wektora μ j
wynosi 65536, a więc – przy identycznej długości wektora p j macierz r jk zawiera 65536 x
65536 elementów, tj. ponad 4 miliardy elementów. Odwracanie tak wielkiej macierzy jest
trudnością samą w sobie – problemem numerycznym (zapewnienie odpowiedniej dokładności
odwracania), ale też i czasowym. Po drugie, zanim zabierzemy się do rekonstrukcji musimy
3
wykonać wszystkie pomiary, co wydłuża czas uzyskiwania obrazu. Wreszcie niebagatelną
sprawą są szumy pomiarowe, które mogą bardzo zdeformować wynik.
W dużym przybliżeniu można uzyskiwać obrazy metodą wstecznej projekcji , która polega na
przypisaniu każdemu pixelowi znajdującemu się na danej linii tego samego ułamka
zmierzonej wartości natężenia, tj. jeśli zmierzone natężenie wynosi I , a na tej linii znajduje się
m pixeli, to każdemu pixelowi przypisujemy wartość I/m (alternatywnie możemy każdemu
pixelowi przypisać wartość I , gdyż w końcu sprowadzi się to tylko do normalizacji natężeń
wobrazie). Suma natężeń w każdym pixelu, uzyskanych z każdego pomiaru daje
wyobrażenie o interesującym nas obrazie. Na rys. 10.2 pokazano zastosowanie tej metody do
rekonstrukcji świecenia w wypadku istnienia dwóch świecących punktów i tylko dwóch,
prostopadłych projekcji. Mierząc natężenia wzdłuż kierunków prostopadłych otrzymamy np.
identyczne wartości, powiedzmy jedynkowe, jak na rysunku z lewej strony.
0
0 1/6 0 1/6 0
1
1/6 1/6 1/3 1/6 1/3 1/6
0
0 1/6 0 1/6 0
1
0
0 1/6 0 1/6 0
1/6 1/6 1/3 1/6 1/3 1/6
1 1
0
0 1/6 0 1/6
Rys. 10.2 Odtwarzanie obrazu punktów świecących (czerwone) z dwóch prostopadłych
projekcji. Po lewej stronie rysunku pokazano miejsca świecenia (czerwone punkty)
i natężenia zmierzone wzdłuż projekcji. Po prawej stronie pokazujemy wynik rekonstrukcji
metodą wstecznej projekcji.
Postępując zgodnie z algorytmem wstecznej projekcji, pixelom na liniach drugiej i piątej
(licząc od góry) przypiszemy wartości 1/6 i podobnie w kolumnach 3 i 5 (od lewej). Łatwo
sprawdzić, że po dodaniu obu wartości miejscom świecenia (punkty) przypisze się w ten
sposób wartości 1/3. Taka sama wartość pojawi się w pixelach (2,3) i (5,5). Pixelom nie
leżącym wzdłuż mierzonej projekcji przypiszemy oczywiście natężenia zerowe.
4
73222147.002.png
Dysponowanie tak ograniczoną informacją i prostym algorytmem doprowadziło nas zatem do
znalezienia miejsc świecenia, ale także i artefaktów: smug w wierszach 2 i 5 oraz kolumnach
3 i 5, a także nie istniejących miejsc świecenia o natężeniach identycznych z rzeczywistymi
miejscami świecenia. Łatwo sprawdzić, że wykonanie dodatkowych projekcji pod kątami 45
stopni także pozostawi te artefakty.
W ogólnym przypadku położenia miejsc „gorących” będą silnie rozmyte, a przy okazji
pojawią się inne artefakty, choć o słabszych natężeniach. Wraz ze wzrostem liczby projekcji,
rozmycie miejsca świecenia spada i otrzymywany rozkład natężenia zmienia się jak 1/r (rys.
10.3), co wykażemy nieco później, niemniej jednak jest on na ogół trudny do zaakceptowania
w rzeczywistej diagnostyce przypadków.
Rys. 10.3 Wynik rekonstrukcji punktu świecącego metodą wstecznej projekcji, gdy
wykona się dużą liczbę projekcji
W tym szczególnym wypadku jednego punktu świecącego (pochłaniającego) matematyka
rekonstrukcji wygląda prościej.
=
n
μ
x
,
y
)
=
p
(
j
δφ
,
r
)
δφ
,
(10.11)
j
1
5
(
73222147.003.png
Zgłoś jeśli naruszono regulamin