Search
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 článek. Ze školních počítačů si můžete článek stáhnout kliknutím na ikonku IEEE full text pdf.
Ú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 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.
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.
Existují tři řešení:
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í.
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).
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.
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.
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.
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.
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
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.