====== 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 [[http://old.lf.upol.cz/menu/struktura-lf/kliniky/klinika-nuklearni-mediciny/pedagogicka-cinnost/fyzikalni-zaklady-zobrazovani-v-nuklearni-medicine-a-radiacni-ochrana/podrobny-obsah-webovych-stranek-projektu-frvs-20522007/|LF UPOL]].

===== Domácí cvičení =====

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 {{:courses:a6m33zsl:petspect.pdf|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í <del>podaného</del> vyrobeného radiofarmaka je //5e-12 molu//,
    * poločas rozpadu dané látky je //110 minut//

  * **[2 b]** Vytvořte funkci //{{courses:a6m33zsl:petmeasure.m|[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í <del>podaného</del> 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]]. 


==== Tvorba signálu a jeho měření ====

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$. /* 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.


==== Implementace - Měření ====

Vstupní balíček {{ :courses:zsl:lab13_pet_data.zip |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:**
  - 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]$ 
  - 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$).
  - 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).

/*
Uvědomte si, že matice aktivovaných detektorů <latex>R</latex> by měla být symetrická (//[i,j]// odpovídá //[j,i]//). Pokud bychom však započetli každou detekci pro obě dvojice, bude matice //R// obsahovat <latex>2\Delta N</latex> rozpadů. Každý rozpad je tedy nutné započítat jen jednou, a matici pak upravit tak, aby byla symetrická (pomocí funkcí //tril// a //triu// a transpozice). Zkontrolujte si že <latex>\sum_{i,j}{R}=\Delta N</latex> . Přibližná podoba horní poloviny matice //R// je na obrázku dolní polovina je vynulovaná.



{{:courses:a6m33zsl:pet_fantom.png?375 |Původní fantom hustoty ropadu.}}
{{ :courses:a6m33zsl:pet_projspace.png?460|Výsledná matice spoludetekovaných fotonů (horní trojúhelníková).}}

 

Jelikož se při tvorbě PET signálu uplatňují náhodné jevy, bude i matice //R// pro každý experiment mírně odlišná. Příklad 300 vybraných paprsků pro obrázek se čtyřmi hot-spoty můžete vidět na obrázku, kde je fantom se 4 body vykazujícími <latex>\beta^+</latex> rozpad (po jednom v levém a pravém rohu a dva uprostřed). Vlevo---zelenou čarou je vyznačena kružnice detektorů, zobrazeno je pouze 300 paprsků. Vpravo---průhlednost jednotlivých pásů určených dvojicí detektorů odpovídá počtu zachycených rozpadů v tomto pásu. Tmavší odstíny odpovídají více detekcím. 


{{courses:a6m33zsl:pet_paprsky.png?420 |Visualizace vybraných 300 parsků rozpadů.}}
{{ :courses:a6m33zsl:pet_bprojekce.png?380|Obrázek přibližně odpovídá výsledku rekonstrukce měřených dat se čtyřmi aktivními oblastmi pomocí zpětné projekce}}
*/

==== Implementace - Rekonstrukce ====

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í. {{ :courses:zsl:petrecon.png?600 |}}
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: <code matlab>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');
</code>


{{tag>medical pet matlab}}