Laboratorium „Elektronika Medyczna”

dr inz. Jerzy Ihnatowicz

Instrukcja do cwiczenia „Tomograf komputerowy”

1. Wstep

Medyczny tomograf komputerowy (CT) jest urzadzeniem do nieinwazyjnej wizualizacji obrazów przekrojów (2D) ciala czlowieka przy wykorzystaniu promieniowania przenikliwego (zasadniczo promieniowania rentgenowskiego). „Nieinwazyjnosc” w pojeciu przyjetym przez srodowisko lekarskie okresla te procedury i techniki diagnostyczne, podczas których nie zostaje naruszona ciaglosc powlok skórnych ciala ludzkiego; diagnostyka nieinwazyjna jest wiec dzialaniem bezkrwawym. Tak rozumiana diagnostyka nieinwazyjna moze jednak stanowic obciazenie organizmu pacjenta niepozadanymi skutkami promieniowania lub/i dzialaniem srodków chemicznych podawanych wziewnie, doustnie lub dorektalnie.

W przypadku tomografii komputerowej CT (CT jest akronimem okreslenia angielskiego Computed Tompgraphy) glównym obciazeniem organizmu diagnozowanego pacjenta jest stosunkowo duza wartosc pochlonietej dawki promieniowania rentgenowskiego, zalezna od energii fotonów promieniowania X oraz zadanej dokladnosci wizualizacji obrazu przekroju.

Nadto, jesli scenariusz badania przewiduje uzyskanie obrazu wielu róznych przekrojów to obciazenie organizmu pacjenta moze byc duze.

Przez „obraz przekroju” rozumie sie rozklad liniowego wspólczynnika pochlaniania promieniowania rentgenowskiego µ( x, y, z ) [cm2/g] dla badanych struktur ciala ludzkiego 0

lezacych wewnatrz dwóch równoleglych plaszczyzn oddalonych miedzy soba o stala odleglosc s okreslana zwykle jako „grubosc warstwy” (lub z angielskiego „slice”). Parametry ( ,

x y, z ) okreslaja polozenie badanej struktury w kartezjanskim ukladzie wspólrzednych 0

0XYZ. Celem uproszczenia procesu tworzenia obrazu przyjmuje sie zazwyczaj stala wartosc z oraz zaklada sie, ze s = 0. Jesli wartosci µ( x, y) zostana przedstawione w postaci pikseli 0

obrazu cyfrowego o jasnosci zaleznej od wartosci µ( x, y) , to tak uzyskany obraz parametryczny (który bedziemy okreslac jako obraz tomograficzny CT) stanowi najczestsza postac informacji diagnostycznej z urzadzen CT.

Obraz tomograficzny CT stanowi dwuwymiarowe odwzorowanie stosunków topograficznych struktur o róznych wartosciach wspólczynnika µ. Dokladnosc odwzorowania przekroju anatomicznego jest wiec ograniczona zmiennoscia wartosci µ dla róznych rodzajów tkanek znajdujacych sie w strefie przekroju.

2. Równanie rzutu

Podstawowym zadaniem tomografii CT jest odtworzenie obrazu przekroju, to znaczy wyznaczenie estymat skonczonego zbioru wartosci

µ( x, y) . Danymi wejsciowymi sa

wartosci rzutów, to znaczy wartosci calek z pochlaniania promieniowania rentgenowskiego wzdluz prostych lezacych w plaszczyznie przekroju. Rzuty zawieraja scalkowane wartosci µ( x, y) . Korzystajac z typowego modelu procesu tworzenia rzutów, przedstawionego na rysunku 1, mozemy zapisac tak zwane równanie rzutu (lub „projekcji”) p( l, θ ) ,

wyznaczonego wzdluz prostej L lezacej w odleglosci l od poczatku ukladu wspólrzednych X0Y, pod katem ? wzgledem osi Y:

y

z

Θ

µ ( r,φ)

z

z = 0

.

r

l

Φ

Θ

0

x

L (prosta rzutowania )

Rys. 1. Model tworzenia rzutu wg Radona.

+

2

2

z

p ( l , ? ) =

µ( l( + z ), θ + ar ctg dz

)

∫ ∞

dla l ≠ 0

∞

−

l

[ 1 ]

+

π

p ( l , ? ) =

µ( z, θ +

dz

)

∫ ∞

dla l = 0

∞

−

2

Zauwazmy, ze para ( l,? ) w sposób jednoznaczny wyznacza polozenie prostej rzutowania L.

Wartosciami rzutów sa calki z pochlaniania promieniowania rentgenowskiego wzdluz prostej L. Przyjety model, zaproponowany przez Radona, nie stoi w sprzecznosci z prawem Lamberta-Beer’a gdyz rozpatrujac pochlanianie wnoszone przez kolejne warstwy absorbenta mozemy – zgodnie z oznaczeniami jak na rysunku 2 zapisac: I0

x1

µ 1

x

µ

2

2

X

x

µ

3

3

..

.

.

..

x

µ

n

n

In

L

Rys. 2. Pochlanianie promieniowania przez kolejne, homogeniczne warstwy absorbenta.

− x ⋅1 µ 1

− x 2 µ

⋅ 2

− x ⋅3 µ 3

− x µ

⋅

n

n

I = I

.

0 ⋅ e

⋅ e

⋅ e

⋅L ⋅ e

n

[ 2 ]

Wyrazenie [2] mozna przeksztalcic do postaci:

− x ⋅ µ + x ⋅ µ + +

L x ⋅ µ

I

(

)

1

1

2

2

e

n

n

n

=

,

I 0

a nastepnie

I

x ⋅ µ + x ⋅ µ +L + x 0

⋅ µ = ln

.

1

1

2

2

n

n

I n

[ 3 ]

Jezeli zalozymy, ze bedziemy rozpatrywali warstwy o grubosciach x

, to równanie [3]

i → 0

mozemy zapisac w postaci calki

x

I

( x) dx = ln 0 = (

p L)

∫ µ

.

I

0

n

[ 4 ]

I

Porównujac uzyskany rezultat z równaniami rzutu [1], latwo zauwazyc, ze wartosc 0

ln

jest

I n

wartoscia rzutu p(L) wzdluz prostej rzutowania L. Oczywiscie, wartosc ta moze byc latwo zmierzona w rzeczywistym ukladzie tomografu komputerowego. Róznica postaci równania rzutu [4] oraz równania rzutu [1] wynika z tego, ze w modelu Radona uwzgledniono dodatkowo zapis polozenia prostej rzutowania L w ukladzie wspólrzednych X0Y.

3. Zadanie rekonstrukcji

Zadanie rekonstrukcji obrazu tomograficznego sprowadza sie do znalezienia zbioru wartosci µ( r,Φ) lub µ( x, y) na podstawie zbioru wartosci rzutów p( l, θ ) . Istnieje kilka podstawowych metod rozwiazanie tego problemu: metoda wykorzystujaca transformacje Fouriera, metoda wykorzystujaca transformacje Radona oraz grupa metod iteracyjnych, z których najbardziej popularna jest tak zwana metoda ART (Algebraic Reconstruction Technique).

Metoda ART polega na rozwiazaniu ukladu równan liniowych, w których wyrazy wolne sa wartosciami rzutów p( l, θ ) , szukane wartosci zmiennych to µ( x, y) zas wspólczynniki przy niewiadomych

µ( x, y) tworza tak zwana siatke modelu rekonstrukcji. Siatka modelu rekonstrukcji jest zwykle siatka n × n jednakowych, kwadratowych oczek o boku d, której postac narzuca konstrukcje ukladu równan i wybór zarejestrowanych przez detektory

promieniowania rentgenowskiego rzutów (wartosci logarytmów naturalnych stosunku I

natezenia promieniowania dochodzacego do absorbenta I i detektora I : 0

ln

).

0

n

I n

Rozwazmy przykladowa siatke n × n oczek o boku d. Liczba nieznanych wartosci µ( x, y) wynoszaca 2

n wymaga do ich okreslenia 2

n równan oraz 2

n wartosci rzutów. W zaleznosci

od polozenia prostych rzutowania L, wspólczynniki a przy zmiennych µ przyjmuja ij

ij

okreslone wartosci z przedzialu

,

0

2 ⋅ d

. Ilustruje to ponizszy rysunek 3.

d

µ

. . . µ

11

1 n

. . .

. . .

d

n

. . .

... ... ... ...

...

µ

. . .

n 1

µnn

n

Lj

Lk

Rys. 3. Przechodzenie prostej rzutowania L przez siatke.

Zapiszmy kilka równan dla zadania rekonstrukcji obrazu o rozdzielczosci 3 × 3 piksele.

Oznaczajac nieznane (poszukiwane) wartosci µ odpowiednio µ , µ ,L, µ , jak na rysunku ij

1

2

9

4, latwo jest ulozyc równania rzutów L oraz L : 2

5

µ 1 µ 2 µ 3

µ 4 µ 5 µ 6

µ 7 µ 8 µ 9

L

L

2

5

Rys. 4. Przykladowe rzuty w modelu rekonstrukcji 3x3.

L : µ ⋅ d + µ ⋅ d + µ ⋅ d = p ; 2

2

5

8

2

[ 5 ]

L : µ ⋅ 2 ⋅ d + µ ⋅ 2 ⋅ d = p .

5

4

8

5

Zauwazmy, ze jesli wartosci rzutów podzielimy przez stala siatki d, to wartosci wspólczynników przy zmiennych µ mozna latwo okreslic jako znormalizowana dlugosc i

odcinka prostej rzutowania L przechodzacego przez j-te oczko siatki, w której przyjeto stala i

wartosc wspólczynnika pochlaniania promieniowania µ . Przykladowo, dla rzutu L mamy k

2

p

L

2

: 0 ⋅ µ + 1⋅ µ + 0 ⋅ µ + 0 ⋅ µ + 1⋅ µ + 0 ⋅ µ + 0 ⋅ µ + 1⋅ µ + 0 ⋅ µ =

.

2

1

2

3

4

5

6

7

8

9

d

Wektor a =

bedziemy okreslali jako wektor wspólczynników siatki dla 2

[ ,1,

0 0,

,

1

,

0 0,

0

,

1

,

0

]

rzutu L .

2

p

1

I

Wyraz wolny

2

0

P =

= ⋅ln

, okreslali bedziemy jako znormalizowana wartosc rzutu L .

2

2

d

d

I 2

W przypadku tomografii komputerowej CT zachodzi, ze liczba dostepnych wartosci rzutów jest zwykle wieksza lub znacznie wieksza od liczby poszukiwanych wartosci niewiadomych µ. Dla siatki n × n oczek (czyli dla rekonstruowanego obrazu o rozdzielczosci n × n pikseli) bedzie zatem

2

k 〉 n ,

[ 6 ]

gdzie k jest liczba rzutów.

Ta nadmiarowosc, wynikajaca zwykle z geometrii urzadzenia tomograficznego (liczby polozen ukladu zródlo-detektory, liczby detektorów, odleglosci miedzy nimi, itp.) jest

niezbedna, gdyz w rzeczywistym ukladzie urzadzenia do tomografii komputerowej CT kazda ze zmierzonych wartosci rzutów jest obarczona pewnym nieznanym bledem. Poszukiwany zbiór rozwiazan { µ , µ , µ ,L, µ

winien jednakowo dobrze spelniac kazde z k równan 2

1

2

3

n }

rzutów. Mamy wiec do czynienia z zagadnieniem poszukiwania wspólrzednych globalnego minimum w przestrzeni

2

n wymiarowej. Jedna z bardziej skutecznych i szybko zbieznych metod iteracyjnego rozwiazania tego zadania jest metoda ART.

4. Metoda ART rekonstrukcji obrazu n×n pikseli

Niech 2

n elementowy wektor µ zawiera wartosci poszukiwanych wspólczynników µ

i

µ = [ µ , µ , µ ,L, µ , µ ; 2

2

1

2

3

n

1

−

n ]

[ 7 ]

oraz niech k elementowy wektor p okresla wartosci znormalizowanych rzutów p = [ p , p , p ,L, p L

,

2 ,

, p

1

2

3

k

n

]

[ 8 ]

gdzie

2

k > n .

Z kazdym z k rzutów skojarzony jest 2

n elementowy wektor a . Zbiór k wektorów i

{ a , a , a , ,

L a

1

2

3

k }

zawiera lacznie

2

k ⋅ n elementów, stanowiacych pelnie informacji o geometrii rzutów w tomografie komputerowym:

a = a , a , ,

L a

1

[

2

1

,

1

,

1 2

,

1 n ]

a = a , a , ,

L a

2

[

2

1

,

2

,

2 2

,

2 n ]

M

a = a a L a

i

[ , , , 2

i 1

,

i,2

i, n ]

M

a = a

a

L a

k

[ , , , 2

k 1

,

k ,2

k, n ]

[ 9 ]

Oznaczajac dodatkowo, ze µ jest wektorem wartosci wspólczynników µ = µ µ L µ

i

[ , , , 2

1

2

n ]

i

wyznaczonych na podstawie i-tego równania rzutu ( i =

,

2

,

1

L, k ) oraz, ze indeks górny (l) jest

numerem kolejnej iteracji, algorytm rekonstrukcji obrazu tomograficznego z rzutów mozna zapisac nastepujaco:

)

1

(

µ

= 0 ;

0

( l

T

)

( l)

l

( )

p − a

⋅ µ

i

i

i

µ

= µ

;

1

+

−1

−

⋅ a

i

i

T

i

a

⋅ a

i

i

[ 10 ]

( l 1

+ )

( l )

µ

= µ

.

0

k

Oczywiscie, i =

,

2

,

1

L, k .

5. Cel cwiczenia

Celem cwiczenia jest zapoznanie sie cwiczacych z metoda ART na przykladzie rekonstrukcji obrazu przekroju z rozdzielczoscia 3 × 3 piksele, na podstawie wartosci 16 rzutów uzyskanych w tomografie komputerowym o konfiguracji skanera przyrostowego z wiazka szpilkowa, co schematycznie ilustruje ponizszy rysunek 5: S

S

D

β

O

O

D

d 1

β = 45o

β = 90o

Rys. 5. Skaner CT przyrostowy z wiazka szpilkowa.

Skaner przyrostowy z wiazka szpilkowa charakteryzuje sie pojedynczym ukladem zródlo (S)

– detektor (D). Oslabienie skolimowanej wiazki promieniowania emitowanej przez zródlo (S) jest mierzone przy wykorzystaniu detektora (D). Odczyt z detektora sluzy do wyznaczenia wartosci rzutu p (→ równanie [4]).

Srodek obszaru rekonstruowanego obrazu przekroju O pokrywa sie ze srodkiem obrotu ukladu zródlo-detektor. Dla kazdego z polozen katowych β, uklad zródlo-detektor przesuwa sie z krokiem d miedzy dwoma skrajnymi polozeniami, tak ustalonymi, aby mozna bylo 1

zmierzyc wartosci rzutów dla calego obszaru przekroju.

Latwo zauwazyc, ze dla M polozen katowych oraz (2 ⋅ N + ) 1 zarejestrowanych rzutów w

kazdym z tych polozen uzyskamy M ⋅ (2 ⋅ N + ) 1 róznych wartosci rzutów, które pozwola na

odtworzenie obrazu przekroju o rozdzielczosci n × n [pikseli] przy czym musi zachodzic, ze 2

n < M ⋅ (2 ⋅ N + )

1 .

W rozpatrywanym w cwiczeniu przykladzie przyjeto, ze kat β przyjmowal wartosci o

o

o

o

β = 0 45

,

90

,

135

,

oraz ze krok przesuwu d ukladu zródlo-detektor wyniósl 1

d

d =

,

1

2

gdzie d jest bokiem oczek siatki.

Poszczególne rzuty (lacznie 16 rzutów p , p ,

,

L p ) zarejestrowano w konfiguracji jak na

1

2

16

rysunku 6.

d

2

d

2

p 8

p 7

p

p

p

p

p

p

1

2

3

4

5

6

β = 0o

β = 45o

p

p

p

16

15

14

p 13

p 11

p

p

10

12

p 9

β = 90o

β = 135o

d

µ

Rekonstruowany

1

µ 2

µ 3

obraz przekroju

d

µ 4 µ µ 6

3x3 piksele

5

µ

µ

µ

7

8

9

Rys. 6. Model tomografu komputerowego rozpatrywany w cwiczeniu.

Obraz zrekonstruowanego przekroju ( 3 × 3 piksele) na podstawie 16 rzutów ( p , p ,

,

L p )

1

2

16

cwiczacy maja wyznaczyc na podstawie tablicy C liczb rzeczywistych zawierajacej wartosci 16 rzutów.

Wyznaczanie wartosci rzutów dla dowolnie zdefiniowanego obrazu przekroju oraz proces rekonstrukcji tego obrazu na podstawie rzutów, metoda ART. odbywa sie w srodowisku programowym Vidas 2.5, przy wykorzystaniu gotowych procedur. Uruchamianie tych procedur sprowadza sie do napisania ich nazwy oraz nacisniecia klawisza ‘Enter’. Konieczne

modyfikacje procesu rekonstrukcji sa realizowane przez wpisanie prostych konstrukcji programowania w uproszczonym jezyku C.

6. Program cwiczenia

• Zapoznac sie z obsluga stanowiska laboratoryjnego i uruchamianiem procedur wyznaczania wartosci rzutów i obrazów rekonstruowanych;

• Korzystajac z procedury proj3x3 wygenerowac rzuty obrazu przekroju;

• Napisac uklad 16 równan rzutów (dla przypadku omawianego w instrukcji do cwiczenia) postaci:

µ ⋅aT

=

pi

dla i = 1 , 2 , 3 , . . . , 16;

• Korzystajac z procedury art3x3 przeprowadzic rekonstrukcje obrazu przekroju i zanotowac wartosci bledu sredniokwadratowego (SRK) dla kolejnych iteracji;

• Zamodelowac wystepowanie bledów detekcji promieniowania rentgenowskiego (skutkujacych bledami wartosci rzutów) o wartosciach +1 , 2 , 5 , 10 , 20 , 30 %

wzgledem wartosci poprawnej. Modelowanie wystepowania bledów zrealizowac za pomoca ponizszej sekwencji polecen:

for i = 1 , i <= 9 , i = i + 1 : R[i] = 0.0

# warunek poczatkowy

dc := 0.0 # zmienna dc okresla blad dc = 0.01 # blad 1 %, zmieniac odpowiednio 0.02, itd.

for i = 1 , i <=16 , i = i + 1

C[i] = C[i] * (1.0 + dc)

endfor

• Przeprowadzic rekonstrukcje obrazu przekroju ze znieksztalconych wartosci rzutów; zanotowac wartosci bledu sredniokwadratowego dla kolejnych iteracji;

• Okreslic dopuszczalny (w [%]) poziom bledów detekcji (pomiaru wartosci rzutów) nie powodujacy znieksztalcen obrazu rekonstruowanego;

• Sformulowac i zapisac wnioski z przeprowadzonych doswiadczen.

7. Pytania sprawdzajace

• Co rozumie sie przez nieinwazyjne medyczne techniki diagnostyczne?

• Czy (a jesli tak, to dlaczego?) nieinwazyjne medyczne techniki diagnostyczne moga wplywac niekorzystnie na organizm pacjenta?

• Okreslic podstawowe zadanie tomografii CT.

• Co oznacza pojecie „przekrój anatomiczny”?

• Co i dlaczego stanowi podstawowe ograniczenie dokladnosci odwzorowania obrazu przekroju anatomicznego w tomografii CT?