Lab 13 Pozitronová emisní tomografie (PET) - měření a rekonstrukce dát

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.

• [1 b] Určete analyticky i numericky počet rozpadlých částic po dobu měření (viz přednáška nebo část Tvorba signálu a jeho měření) pokud
  • ve 30. minutě byla podána látka pacientovi,
  • ve 35. minutě začal být pacient snímán,
  • ve 39. minutě skončilo snímání pacienta,
  • množství podaného vyrobeného radiofarmaka je 5e-12 molu,
  • poločas rozpadu dané látky je 110 minut

• [2 b] Vytvořte funkci [deltaN,detPair]=petMeasure(act_map,ti,tb,te,Nd,r,n,sigma2,taur), kde vstupem vaší funkce bude virtuální fantom reprezentující aktivitu v každém pixelu a další parametry měření. Ke každé veličině jsou v závorkách uvedeny hodnoty pro naši konkrétní úlohu, podrobněji jsou parametry vysvětleny v části Tvorba signálu a jeho měření, postup pak v části Implementace - Měření). Do reportu vložte obrázek hustoty rozpadů (detPair) vytvořený z fantomu aktivit a dále obrázek symetrické matice detPair (rozdílné hustoty zakódujte barevně pomocí vhodně zvolené colormap). Parametry funkce jsou:
  • deltaN - počet rozpadlých částic po celou dobu snímán pacienta
  • detPair - mřížka s kumulovanými počty rozpadů pro každou kombinaci dvojic detektorů (stejná velikost jako je počet detektorů [Nd Nd])
  • fantom - fantom koncentrace radiofarmaka (obrázek o velikosti [128 128])
  • ti - čas ve kterém byla látka pacientovy podána (30. minuta)
  • tb - čas kdy se začal pacient snímat (35. minuta)
  • te - čas kdy skončilo snímání pacienta (39. minuta)
  • Nd - počet detektoru na celém obvodu kružnice (100)
  • r - poloměr kružnice na které jsou rozmístěny detektory (pro celý fantom)
  • n - množství podaného vyrobeného radiofarmaka (5e-12 molu)
  • sigma2 - rozptyl, jak daleko doletí částice před svým rozpadem (3px)
  • taur - poločas rozpadu dané látky (110 minut)

• [2 b] Napiště funkci Img = PETReconstruction(detPair, R) a zrekonstruujte naměřená data z matice detPair, rekonstruovanou matici I vizualizujte, idálně jako overlay nad strukturálním fantomem. Dále spočtěte a visualizujte relativní chybu mezi fantomem aktivit z předcházejícího kroku a zrekonstruovaným fantomem získaným z PETReconstruction. Detailní popis k rekonstrukci je v sekci Implementace - Rekonstrukce.

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 rozpadlý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$.

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}$, 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:

1. Rozdělte kružnici opsanou fantomu (s poloměrem $r$ a středem uprostřed fantomu/obrázku) na $D$ částí (detektorů). Detektory číslujte od jedničky (standard Matlabu) po směru hodinových ručiček od pozice $[0, r]$
2. Vypočítejte kolik atomů radioizotopu se rozpadne za měřený interval v jedné vrstvě pacientova těla, pak korigujte toto číslo pro dráhy fotonů mířící mimo detektor a obdržíte $\Delta N$). Rozdělte celkový počet rozpadů na jednotlivé pozice podle hodnot fantomu (fantom aktivit znormalizujte a přenásobte $\Delta N$).
3. Pro každý rozpad určete:
   • směr letu pozitronu a vzdálenost anihilace, předpokládejte že místo anihilace je určeno Gaussovou funkcí s diagonální kovarianční maticí. Oba prvky diagonály nechť jsou $\sigma^2$ px (tento náhodný jev se dá snadno modelovat pomocí filtrace fantomu Gaussovským filtrem s danou kovarianční maticí. Filtr vytvoříte pomocí funkce fspecial a aplikujete přes metodu imfilter),
   • směr letu fotonů (každému rozpadu přiřaďte náhodný úhel od 0 do $\pi$ ),
   • určete průsečík trajektorie fotonů s detektorovou kružnicí (můžete použít funkci 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ů),
   • vytvořte výstupní matici $P$ o rozměrech $D \times D$, která bude obsahovat počet excitací pro každou dvojici synchronně aktivovaných detektorů.

- Naprogramovanou funkci aplikujte na fantom aktivit s parametry ze Domácí cvičení, 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');