====== Pozitronová emisní tomografie - měření dat ======

Pozitronová emisní tomografie je zobrazovací metoda nukleární medicíny. Tomografický obraz je rekonstruován díky současné detekci dvou fotonů, které byly vyzářeny při anihilaci pozitronu vzniklého při rozpadu radiofarmaka v pacientově těle.
/* (<latex>\beta^+</latex> rozpad). Tyto dva fotony se od místa anihilace pohybují ve vzájemném úhlu téměř <latex> 180^\circ</latex> a proto lze snadno zjistit přímku na které k jevu došlo.

{{ :courses:a6m33zsl:pet_anihilace.png?300 |Schéma anihilace páru elektron-pozitron}}

Pokud se podaří nasbírat dostatek rozpadů, můžeme z takových dat rekonstruovat rozložení jevů (anihilací) v prostoru. Je dobré si uvědomit, že měříme místa anihilace, nikoli skutečné pozice rozpadu radiofarmaka, což v principu omezuje rozlišovací schopnost PET. Pozitron v praxi ulétne většinou 2-3 mm, tato systematická chyba tedy není příliš velká.

Měření signálu je v PET prováděno pomocí scintilačních krystalů, které přemění vysokoenergetický foton na světelné, detekovatelné běžnými fotodetektory. Nastudujte si {{:courses:a6m33zsl:petspect.pdf|přednášku}} o PET.
*/
Názorný popis principu PET najdete na webu [[http://www.lf.upol.cz/menu/struktura-lf/kliniky/klinika-nuklearni-mediciny/pedagogicka-cinnost/fyzikalni-zaklady-zobrazovani-v-nuklearni-medicine-a-radiacni-ochrana/pozitronova-emisni-tomografie/principy-pet/|Principy PET]]

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

Úkolem následujícího cvičení je demonstrovat tvorbu signálu v pozitronové emisní tomografii. Potřebná data a kódy si stáhněte {{courses:a6m33zsl:pet_datacode.zip?lïnkonly|zde}}.

  - Určete počet (analyticky i numericky) rozpadlých částic po dobu měření (viz {{:courses:a6m33zsl:petspect.pdf|přednáška}} nebo část [[#Tvorba signálu a jeho měření]]). [0.5b.]
    * 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 radiofarmaka je //5e-12 molu//,
    * poločas rozpadu dané látky je //110 minut//
  - Vytvořte funkci //{{courses:a6m33zsl:petmeasure.m|[deltaN,detPair]=petMeasure(fantom,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í. Fantom si vytvořte pomocí funkce //phantom(p_def);// a matice definující fantomy, které jsou v souboru //ph_def.mat// (který jste si stáhli již dříve). Ke každé veličině jsou v závorkách uvedeny hodnoty, které použijete (více viz popis v [[#Tvorba signálu a jeho měření]]). [3b.]
    * //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 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)//
  - Do reportu vložte obrázek hustoty rozpadů (//detPair//) vytvořený z fantomu aktivit a symetrické matice //detPair// (pro číselnou reprezentaci jednotlivých barevných odstínů použijte funkci //[[http://www.mathworks.com/help/matlab/ref/colorbar.html|colorbar]]//). [0.5b.]
  - Rekonstruujte naměřenou matici //detPair// funkcí //PETreconstruction(detPair,r,Nd);// kterou jste si stáhli na začátku cvičení. Spočtete a visualizujte relativní chybu mezi rekonstruovaným fantomem z vašeho měření a originálním fantomem. [1b.]

Funkce které jste si stáhli na začátku cvičení neodevzdávejte.

===== 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 <latex>T_i</latex> [s], čas začátku měření <latex>T_b</latex> [s], čas konce měření <latex>T_e</latex> [s], počet detektorů //D//, poloměr detektorové kružnice //r//[px], rozptyl pozitronů <latex>\sigma^2</latex> [px], množství //n//[mol], poločas rozpadu radiofarmaka <latex>\tau_r</latex> [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 <latex>\Delta N?</latex>). 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ů <latex>\Delta N</latex> . 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 <latex>t=0</latex>  bylo vyrobeno radiofarmakum, které pak bylo v čase <latex>T_i</latex>  podáno pacientovi.

Počet vyrobených molekul je <latex>N_0 =n \cdot N_A</latex>,kde <latex>N_A</latex> je Avogadrova konstanta <latex>N_A=6.0221415\cdot 10^{23} \textrm{mol}^{-1}</latex>, ale v čase <latex>T_i</latex> se už část radioizotopu rozpadla.

Uvažujme, že za čas <latex>dt</latex> se v látce rozpadne <latex>dN</latex> částic <latex>-\ud N=\lambda_r N \ud t \label{difeq}</latex>, kde <latex>\lambda_r</latex> je rozpadová konstanta <latex>\lambda_r=\frac{\ln{2}}{\tau_r}</latex>. Řešením rovnice je známá exponenciální závislost počtu izotopů na čase <latex>N_1=N_0\textrm{e}^{-\lambda_r t}</latex>.

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 <latex>T_b</latex>. 
/* Protože je však aktivní látka od času <latex>T_i</latex> do konce experimentu v lidském těle, začne se zákonitě vylučovat. Analogicky k rovnici se vylučuje rychlostí <latex>-d N=\lambda_v N d t</latex> tentokrát však s vylučovací konstantou <latex>\lambda_v</latex>. Oba tyto jevy probíhají současně, úbytek aktivních částic tedy bude <latex>-d N=(\lambda_r + \lambda_v) N d t</latex>.

Stejným postupem jako v případě předešlé rovnice snadno získáme <latex>N_2=N_1e^{-(\lambda_r+\lambda_v) t}</latex>. Naše úloha však 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. Tato bioaktivní látka se váže v tkáni a proto je její farmakokinetika jiná než popisuje dříve zmíněný zjednodušený model.
*/

Protože je však aktivní látka od času <latex>T_i</latex> do konce experimentu v lidském těle, začne se zákonitě vylučovat. Naše úloha však 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.

Přibližně 75% FDG je zachyceno ve tkáni a rozpadá se s <latex>\tau_r=110</latex> min, 25% látky se vyloučí ledvinami s <latex>\tau_v=16</latex> min. To je samozřejmě možné modelovat, ale vzhledem ke krátkému <latex>\tau_v</latex> 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ě <latex>1^\circ</latex> , detektor tedy zasáhnou jen 2/180 fotonů vzniklých v řezu.


===== Implementace =====

  - Rozdělte kružnici opsanou fantomu (s poloměrem //r// a středem uprostřed fantomu/obrázku) na <latex>D</latex> částí (detektorů). Detektory číslujte od jedničky (standard Matlabu) proti směru hodinových ručiček od úhlu <latex>\alpha=0</latex> viz obrázek. {{ :courses:a6m33zsl:pet_deleni.png?300 |Číslování jednotlivých detektorů PET pro D=6.}}
  - Vypočítejte kolik atomů radioizotopu se rozpadne za měřený interval v jedné vrstvě pacientova těla. Korigujte toto číslo pro dráhy fotonů mířící mimo detektor (výsledkem je <latex>\Delta N</latex> ). Rozdělte spočtené rozpady podle fantomu (vynásobte fantom normalizovaný na jednotkový součet celkovým počtem rozpadů).
  - 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 <latex>\sigma^2</latex> px (tento náhodný jev se dá snadno modelovat pomocí filtrace fantomu Gaussovským filtrem s danou kovarianční maticí. Filtr vytvoříte snadno pomocí funkce //fspecial// a aplikujete ho na data pomocí //imfilter//),
    * směr letu fotonů (každému rozpadu přiřaïte náhodné číslo od 0 do <latex>\pi</latex> ),
    * průsečík trajektorie fotonů s detektorovou kružnicí (určete které detektory zachytily letící fotony),
    * vytvořte výstupní matici <latex>R</latex> funkce o rozměrech <latex>[D\times D]</latex>, která bude obsahovat počet excitací pro každou dvojici synchronně aktivovaných detektorů.
  - Naprogramovanou funkci vyzkoušejte na fantomu aktivit s parametry (viz [[#Zadání]]), poloměr detektorové kružnice volte tak, aby přesně opsala čtverec fantomu (<latex>r=\frac{1}{2}\cdot\sqrt{ m^2 + n^2 }</latex> 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}}





{{tag>medical pet matlab}}