====== Pozitronová emisní tomografie - rekonstrukce ======

Existuje několik rekonstrukčních algoritmů pro tomografické metody. V dnešním cvičení se podrobněji seznámíte s iterativní rekonstrukcí pomocí EM algoritmu, někdy též označované jako multiplikativní iterativní rekonstrukce. Prostudujte přednášku o PET nebo [[http://www.ncbi.nlm.nih.gov/pubmed/18238264|článek]]. Ze školních počítačů si můžete článek stáhnout kliknutím na ikonku IEEE full text pdf. 

===== Zadání =====

Úkolem dnešního cvičení je demonstrovat rekonstrukci obrazu z dat měřených pomocí pozitronové emisní tomografie. Potřebná data a kódy si stáhněte {{courses:a6m33zsl:pet_reconstruction_datacode.zip|zde}}.

V kódu z minulého cvičení nahraïte námi dodanou funkci //reconstruct// funkcí //reconstructEM// kterou naprogramujete. Vstupem funkce bude matice současně aktivovaných detektorů //R//, poloměr detektorové kružnice //r// a počet detektorů na kružnici //D//. Výstupem rekonstruovaný obraz řezu.

===== Rekonstrukce =====

Označme si rekonstruovaný obrázek //x//, souřadnice pixelů <latex>\mathbf{r}\in V</latex>, //V// je diskrétní podprostor <latex>\mathbb{R}^2</latex>. Symbolem <latex>\mathbf{d}</latex> označme dvojici detektorů <latex>\mathbf{d}=(d_1,d_2)</latex> a množinu všech těchto dvojic označme <latex>\mathbb{D}^2</latex>.

Jádrem rekonstrukčního algoritmu je iterativní vzorec:

<latex>
	x^{new}(\mathbf{r})=x(\mathbf{r})\sum_{\mathbf{d}\in \mathbb{D}^2}{\frac{s(\mathbf{d})\pi(\mathbf{d}|\mathbf{r})}{\sum_{\mathbf{r}'\in V}{x(\mathbf{r}')\pi(\mathbf{d}|\mathbf{r}')}}}\, ,
</latex>

kde //x// je aktuální a <latex>x^{new}</latex> nová rekonstrukce, <latex>s(\mathbf{d})</latex> je počet detekcí dvojicí detektorů <latex>\mathbf{d}</latex> a <latex>\pi(\mathbf{d}|\mathbf{r})</latex> je pravděpodobnost, že byla událost z pixelu <latex>\mathbf{r}</latex> detekovaná párem detektorů <latex>\mathbf{d}</latex>. Matice <latex>s(\mathbf{d})</latex> odpovídá měření, které jste naprogramovali na minulém cvičení. Pravděpodobnost <latex>\pi(\mathbf{d}|\mathbf{r})</latex> je dána geometrií přístroje.

Všimněte si i významu jmenovatele, který vyjadřuje co bychom změřili pokud by měřeným objektem byla naše současná rekonstrukce a čitatele který rozděluje skutečné měření do jednotlivých pixelů na základě <latex>\pi(\mathbf{d}|\mathbf{r})</latex>.

Iterování podle tohoto předpisu není obtížné, problém je efektivně vypočítat a reprezentovat <latex>\pi(\mathbf{d}|\mathbf{r})</latex>. Matice <latex>\pi(\mathbf{d}|\mathbf{r})</latex> má rozměr rekonstruovaného obrázku x počet pásů definovaných dvojicí detektorů, tedy označíme-li počet řádků a sloupců obrazu <latex>m</latex> je rozměr <latex>m\times m\times D(D-1)/2</latex>, což je velká matice (pro <latex>m=256</latex>, <latex>D=100</latex> a data reprezentovaná typem //double// je to 2.42 GB na jednu vrstvu. Takovou matici není 32b program schopen alokovat). To znamená, že pro velké objemy dat nejsme schopni tuto matici uchovávat.

{{ :courses:a6m33zsl:pet_pas.png?400 |Příklady tří detektorových pásů definovaných detektory číslo 2 a 11, 3 a 8, 12 a 13.}}

Existují tři řešení:
  - Počítat jednotlivé části <latex>\pi(\mathbf{d}|\mathbf{r})</latex> vždy znova když jsou potřeba.
  - Omezit se na menší rozlišení v prostoru, bitové hloubce, popřípadě obojím.
  - Využít řídkosti matice <latex>\pi(\mathbf{d}|\mathbf{r})</latex> (většina prvků je nulových).

První řešení prodlouží čas potřebný na jednu iteraci, druhé sníží kvalitu výstupu a třetí je obtížnější na programování (dají se využít sparse knihovny, které tento problém řeší). My použijeme druhou metodu, protože nás nezajímá ani tak vysoké rozlišení, jako spíše princip rekonstrukce, navíc je tato metoda nejsnazší na programování.


==== Výpočet funkci π ====

Parametr <latex>\pi(\mathbf{d}|\mathbf{r}_i)</latex> (přesněji jeho odhad) pro pixel <latex>\mathbf{r}_i</latex> bychom mohli hypoteticky získat vložením bodového zdroje zaření na pozici <latex>\mathbf{r}_i</latex> v měřeném prostoru a normalizací výsledku počtem zaznamenaných rozpadů (samozřejmě že čím více rozpadů, tím přesnější odhad). Jak je zřejmé z tohoto myšleného experimentu, <latex>\pi(\mathbf{d}|\mathbf{r})</latex> je dán jen geometrickým uspořádáním měřidla a dá se z něj vypočítat.

Přesný výpočet <latex>\pi(\mathbf{d}|\mathbf{r}_i)</latex> je komplikovaný a pomalý proto použijeme pro dnešní cvičení dvě aproximace. První z aproximací považuje pravděpodobnost rozpadu v pásu za konstantní (rychlejší a méně přesná, zvláště pokud jsou pásy výrazně širší než 1 px), druhý neuvažuje velikost pixelů (pomalejší a přesnější, zvláště pro široké pásy).

==== Aproximace pomocí konstantních pásů ====


Předpokládejme, že každý pixel vyzáří fotony do všech dvojic detektorů v jejichž pásu se nachází se stejnou pravděpodobností. Pokud leží pixel na hranici dvou pásů, pak je pravděpodobnost pro každý z nich vynásobena plochou pixelu obsaženou v daném pásu. Pro pixely mimo pás je pravděpodobnost nulová. Pro lepší představu se podívejte na obrázek. Matici <latex>\pi_s(\mathbf{d}|\mathbf{r}_i)</latex> a vektor indexů přiřazující každé dvojici detektorů indexy do matice <latex>R</latex> získáte pomocí funkce //makePi_stripes//.

{{ :courses:a6m33zsl:pet_probs1.png?400 |Pravděpodobnost pixelů v jednom detektorovém pásu, aproximace pomocí konstantních pásů.}}

==== Aproximace pomocí úhlového překryvu ====

Reprezentujme každý pixel obrázku jako bod ležící v jeho středu. Takový bod definuje společně se dvěma krajními body detektoru úhel <latex>\alpha</latex>. Uvažujme nyní jeden bod (pixel) <latex>\mathbf{r}</latex> a dvojici detektorů. Úhly definované bodem a každým z detektorů označme <latex>\alpha_1</latex> a <latex>\alpha_2</latex>. Dále definujme úhly <latex>\beta_1</latex> a <latex>\beta_2</latex> jako úhly středově souměrné k <latex>\alpha_1</latex> a <latex>\alpha_2</latex> podle <latex>\mathbf{r}</latex>. Pravděpodobnost vyzáření fotonů z bodu <latex>\mathbf{r}</latex> do dvojice detektorů <latex>\mathbf{d}</latex> je úměrná menšímu z průniků úhlů <latex>\alpha_1\cap\beta_2</latex> či <latex>\alpha_2\cap\beta_1</latex>. Přesněji je to průnik úhlů dělený <latex>2\pi</latex>. Daleko lepší představu o situaci si uděláte z obrázku. Matici <latex>\pi_a(\mathbf{d}|\mathbf{r}_i)</latex> a vektor indexů přiřazující každé dvojici detektorů indexy do matice <latex>R</latex> získáte pomocí funkce //makePi_angles//.

{{ :courses:a6m33zsl:pet_probs2.png?400 |Pravděpodobnost jednoho pixelu pro jeden pár detektorů, aproximace pomocí úhlového překryvu.}}

===== Iterování =====

Naprogramujte výpočet iterativního vzorce. Pro cyklus přes jednotlivé iterace použijte funkci //while//. Pro vlastní výpočet nepoužívejte další cykly, vystačíte si s funkcemi //repmat// a //reshape//. Parametr <latex>s(\mathbf{d})</latex> získáte přerovnáním prvků matice <latex>R</latex> podle mapovací proměnné z našich funkcí opět nepoužívejte cyklus, ale funkci //sub2ind// a následné indexování <latex>R</latex>. Před výpočtem převeïte <latex>R</latex> na horní trojúhelníkovou. Ve výpočtu <latex>x^{new}</latex> se může vyskytnout výraz 0/0, definujme tedy 0/0=0. 

Iterování zastavte když bude <latex>\sum_r{\big( x^{new}(\mathbf{r})-x(\mathbf{r})\big)^2}<\epsilon</latex>. Práh <latex>\epsilon</latex> zvolte tak, aby byl výsledek rekonstrukce uspokojivý, svou volbu zdůvodněte v reportu.  


===== Experimenty =====

Změřte a rekonstruujte virtuální fantom vámi vytvořenou funkcí //reconstructEM// i námi dodanou funkcí //reconstruct// pro čtyři různá nastavení z tabulky~\ref{parametry} a pro obě aproximace <latex>\pi(\mathbf{d}|\mathbf{r}_i)</latex>. Do reportu vložte matici //R// společně s aktivitami rekonstruovanými filtrovanou zpětnou projekcí (//reconstruct//) a obrázky rekonstruované EM algoritmem s využitím obou aproximací <latex>\pi_a</latex>, <latex>\pi_s</latex> (celkem 16 obrázků).

Pozor, obrázky se budou jevit převrácené kolem osy //x//, vlivem rozdílného počátku souřadného systému obrázku a dat. Do reportu vkládejte obrázky orientované stejně jako vstup. Protože bude report plný podobných rekonstrukcí, pečlivě popisujte všechny vyobrazení.

Zhodnoťte vliv parametrů na rekonstrukce a porovnejte kvalitu obou rekonstrukčních funkcí pro obě aproximace. Porovnávejte kontrast, na pozadí rekonstruovannou aktivitu, šum, případné artefakty a váš subjektivní dojem.

^  Proměnná  ^  1 ^  2  ^  3  ^  4  ^
^  <latex>T_i</latex>  |  0:30:00  |  0:30:00  |  0:30:00  |  0:30:00  |
^  <latex>T_b</latex>  |  0:35:00  |  0:35:00  |  0:35:00  |  0:35:00  |
^  <latex>T_e</latex>  |  0:42:00  |  0:42:00  |  0:35:30  |  0:42:00  |
^  //D//  |  90  |  90  |  90  |  45  |
^  <latex>P_{size}</latex>  |  <latex>70\times 70</latex> px  |  <latex>70\times 70</latex> px  |  <latex>70\times 70</latex> px  |  <latex>70\times 70</latex> px  |  
^  <latex>\sigma^2</latex>  |  2  |  0  |  0  |  0  |
^  //n//  |  5e-12 mol  |   5e-12 mol  |   5e-12 mol  |   5e-12 mol  |  
^  <latex>\tau_r</latex>  |  110 min  |  110 min  |  110 min  |  110 min  |  

===== Bonus +1b. =====

Zobrazte rekonstruované aktivity společně (přes sebe) s anatomickým fantomem. Postup bude podobný jako v úloze fMRI, jen s tím rozdílem, že do //AlphaData// přiřadíte místo binární masky matici rekonstruovaných aktivit. Budete také muset nastavit parametr //AlphaDataMapping//. Pro mapování aktivit na jasy použijte lut //hotmy//. Pro kontrolu se podívejte na obrázek

{{courses:a6m33zsl:pet_rspace.png?400|}} {{courses:a6m33zsl:pet_vizual.png?300|}}

Matice //R// a rekonstruované aktivity přeložené přes anatomický fantom pro nastavení: <latex>T_i=0:30:00</latex>, <latex>T_b=0:35:00</latex>, <latex>T_e=0:40:00</latex>, <latex>D=230</latex>, <latex>P_{size}=256\times 256</latex> px, <latex>\sigma^2=2</latex> px, <latex>n=5\cdot 10^{-12}</latex> mol, <latex>\tau_r=110</latex> min.

{{tag>medical pet reconstruction}}