Téma cvičení je měření a rekonstrukce dat pozitronové emisní tomografii (PET). Jádrem zobrazování PET je detekce dvojice fotonů, které vzniknou následkem anihilace pozitronu vzniklého při rozpadu radiofarmaka v pacientově těle. Fotony letí od sebe po přímce a v ideálním případě je jejich průlet současně detekován dvojicí detektorů, které jsou uspořádány prstencovitě kolem snímaného subjektu. V závislosti na koncentraci radiofarmaka budeme simulovat rozpady, náhodné dráhy fotonů a jejich detekci na prstenci s $N$ detektory. V druhé části pak z naměřených dat rekonstruovat 2D snímek – matici detekcí z první části si nejprve převedeme na sinogram a výsledný snímek získáme filtrovanou zpětnou projekcí.
Podrobnější popis principu PET najdete např. na webu LF UPOL.
Domácí cvičení má tři části - nejprve si spočteme počet částic, které vzniknou v průběhu měření, tento výpočet (numericky) použijeme v druhé, simulační části. Nakonec z naměřených dat zrekonstruujeme výsledný snímek.
Naprogramujte funkci, která vypočítá výstup měřícího prstence PET. Vstupem funkce nechť je čas podání radiofarmaka $T_i$ [s], čas začátku měření $T_b$ [s], čas konce měření $T_e$ [s], počet detektorů $D$, poloměr detektorové kružnice $r$[px], rozptyl pozitronů $\sigma^2$ [px], množství $n$[mol], poločas rozpadu radiofarmaka $\tau_r$ [s] a fantom aktivit (relativní distribuce radiofarmaka v těle pacienta, pro absolutní roložení celý fantom znormujeme, aby součet byl 1 a vynásobíme celkovám počtem čozpadlých částic $\Delta N$). Výstupem pak bude matice současně aktivovaných detektorů R která má na souřadnici [i,j] počet současných detekcí elementy číslo i a j a celkový počet měřených rozpadů $\Delta N$. Získaná data si můžete rekonstruovat pomocí funkce activity=PETreconstruction(detPair,r,Nd) ze staženého balíčku.
Předpokládejme že v čase $t=0$ bylo vyrobeno radiofarmakum, které pak bylo v čase $T_i$ podáno pacientovi. Počet vyrobených molekul je $N_0 =n \cdot N_A$,kde $N_A$ je Avogadrova konstanta $N_A=6.0221415\cdot 10^{23} \textrm{mol}^{-1}$, ale v čase $T_i$ se už část radioizotopu rozpadla. Uvažujme, že za čas $dt$ se v látce rozpadne $dN$ částic $$-dN =\lambda_r N d t \label{difeq}$$ kde $\lambda_r$ je rozpadová konstanta $\lambda_r=\frac{\ln{2}}{\tau_r}$. Řešením rovnice je známá exponenciální závislost počtu izotopů na čase $$N_1=N_0\textrm{e}^{-\lambda_r t}$$.
Než se dostane radioaktivní látka do místa určení (často je podávána systémově, například intravenózně), uplyne nějaký čas, proto zahájíme měření až v $T_b$.
Naše úloha bude simulovat měření pomocí fludeoxyglukózy (FDG), která je analogem glukózy s inkorporovaným atomem fluóru 18, který se rozpadá za vzniku pozitronu. Protože je však aktivní látka od času $T_i$ do konce experimentu v lidském těle, začne se zákonitě vylučovat.
Přibližně 75% FDG je zachyceno ve tkáni a rozpadá se s $\tau_r=110$ min, 25% látky se vyloučí ledvinami s $\tau_v=16$ min. To je samozřejmě možné modelovat, ale vzhledem ke krátkému $\tau_v$ a známému poměru zachyceného a vyloučeného radiofarmaka si můžeme práci zjednodušit (při zachování dostatečné přesnosti) a uvažovat jen 75% molekul radiofarmaka.
Každý z rozpadajících se atomů vyzáří pozitron, který se pohybuje náhodným směrem tkání a poté co ztratí svou kinetickou energii, anihiluje s elektronem hmoty. Při tom vzniknou dva fotony s energií přibližně 511 keV, vzdalující se od místa srážky na opačné strany po náhodně orientované přímce. Pokud dvojice fotonů poletí podél osy pacientova těla, nemá pro náš zobrazovací systém význam. Z geometrického uspořádání je jasné, že našeho myšleného prstencového detektoru dosáhne jen malá část záření.
Uvažujme pacienta vysokého 180 cm, z něhož pořizujeme přibližně 5 mm řez. Pokud by bylo radiofarmakum zachycováno rovnoměrně podél osy pacientova těla, pak v našem řezu zůstane přibližně 5/1800 částic. Toto zjednodušeni nám výrazně usnadní práci, při zachování řádové přesnosti.
Náš detektor je tenký prstenec kolmý k ose pacienta a z toho je jasně patrné, že nemůže zachytit fotony, jejichž dráha není téměř kolmá k pacientovu tělu. Takové dvojice fotonů minou detektor a jsou z našeho pohledu nepodstatné, protože nepřispívají k tvorbě obrazu. Definujme že rozdílový úhel (mezi rovinou detektoru a dráhou fotonů) může být maximálně $1^\circ$ , detektor tedy zasáhnou jen 2/180 fotonů vzniklých v řezu.
Vstupní balíček lab13_pet_data.zip obsahuje fantom aktivit a strukturální fantom (slouží pro vizualizaci) ve formátu .mat
a dále funkci
[d1, d2] = DetectorPairFromLOR(i, j, alpha, sz, R, N_det, act_phantom)
. Souřadný systém nastavte tak, aby bod (0, 0) ležel ve středu fantomu. Dále pro zjednodušení uvažujte počet detektorů $D$ dělitelný 4 beze zbytku.
Postup:
DetectorPairFromLOR
, která při zadaných souřadnicích vzniku rozpadu a směrovém úhlu rozpadové přímky spočte průsečíky a vrátí indexy detektorů),
- Naprogramovanou funkci aplikujte na fantom aktivit s parametry ze [#Zadání]], poloměr detektorové kružnice volte tak, aby přesně opsala čtverec fantomu ($r=\frac{1}{2}\cdot\sqrt{ M^2 + N^2 }$ kde M, N jsou velikosti obrázku/fantomu).
Na vstupu máme matici $P$ s počty aktivací párů detektorů a dále poloměr detektorového prstence $r$. Rekonstrukci provedeme ve dvou krocích
Transformace na sinogram Převod změřené matice na sinogram lze provést více způsoby, zde si představíme jeden založený na projekcích. Pro každý projekční úhel $\vartheta$ budeme uvažovat ty projekční paprsky, které procházejí dvojící protilehlých detektorů. Pro tyto body máme změřený počet aktivit v matici $P$. Uvědomte si, že body nejsou na úsečce $[-r, r]$ rovnoměrně rozmístěné (blíže k sobě na krajích, dál od sebe okolo počátku). Do sinogramu, který vzorkuje osu $p$ rovnoměrně, přeneseme tato diskrétní data interpolací. Když pootočíme projekci o úhel, který odpovídá polovině úhlu mezi středy detektorů, tj. pro $\vartheta_2 = \vartheta_1 + D / 2\pi $ bude párování detektorů na projekčních přímkách $[11\:10], [12\:9], [1\:8], [2\:7], [3\:6], [4\:5]$.
Rekonstrukce výsledný rekonstruovaný snimek změřených aktivit získáme ze sinogramu pomocí filtrované zpětné projekce (iradon
)
Vizualizujte výslednou aktivaci jako barevný overlay nad strukturním fantomem, vykreslete obrázek aktivit obarvený pomocí colormap hot. Pro nastavení průhlednosti použijte normalizovanou mapu aktivit, jak je to v ukázce níže:
activity = PETReconstruction(detPair, R); [...] figure; imagesc(struct_phantom); colormap gray; hold on; imshow(activity_colored); ch=get(gca,'Children'); set(ch(1),'AlphaData',activity_normalized); set(ch(1),'AlphaDataMapping','scaled');