10 Podstawowa matematyka rekonstrukcji tomograficznych, Nauka - różności, Fizyka medyczna, Medycyna ...
[ Pobierz całość w formacie PDF ]
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
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
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
(
[ Pobierz całość w formacie PDF ]