SZAKDOLGOZAT
GPU ALAPÚ PET REKONSTRUKCIÓS KÓD TOVÁBBFEJLESZTÉSE IDŐFÜGGŐ MÉRÉSI REKONSTRUKCIÓKRA
KOLESZÁR ÁDÁM
Témavezető:
Dr. Légrády Dávid egyetemi docens BME Nukleáris Technikai Intézet Nukleáris Technika Tanszék
Budapesti Műszaki és Gazdaságtudományi Egyetem 2012
FIZIKA BSC :: SZAKDOLGOZAT :: KIÍRÁS Témavezető: Neve: Tanszéke: E-mail címe: Telefonszáma:
Légrády Dávid NTI
[email protected] 1254
Azonosító: Sz-2012-105 Szakdolgozat
GPU alapú PET rekonstrukciós kód továbbfejlesztése időfüggő mérési
címe:
rekonstrukciókra
Melyik
''Alkalmazott fizika'' ''Fizika''
szakiránynak ajánlott? A
jelentkezővel c++ és CUDA programozási tapasztalat, Monte Carlo módszerek
szemben támasztott elvárások: Leírása:
A hallgatónak a Teratomo projekt keretében kifejlesztett Monte Carlo alapú GPU-n működő PET rekonstrukciós kód fejlesztésébe kell bekapcsolódnia különös tekintettel a C++ környezetben a kód további fejlesztésére a következő új funkciók algorimtusának kifejlesztésére, implementációjára és tesztelésére: - random koincidenciák okozta instabilitások kiküszöbölése - időfüggő felvételek rekonstrukciója magasabb szintű szűrési eljárások
2
ÖNÁLLÓSÁGI NYILATKOZAT Alulírott , a Budapesti Műszaki és Gazdaságtudományi Egyetem Fizikus alapszak (BSc) fizika szakirányának hallgatója kijelentem, hogy ezt a szakdolgozatot meg nem engedett segítség igénybevétele nélkül, saját magam készítettem. Minden olyan szövegrészt, adatot, diagramot, ábrát, vagy bármely más elemet, amelyet vagy szó szerint, vagy azonos értelemben, de átfogalmazva másoktól vettem át, a forrás megadásával egyértelműen megjelöltem.
Budapest, 2012. június 1. ………………………… Koleszár Ádám hallgató
TARTALOM I. Bevezető ......................................................................................................................................................... 1 I.1 PET ............................................................................................................................................................ 1 I.2 Tomográfiás eljárások ....................................................................................................................... 2 I.3 Monte Carlo módszerek .................................................................................................................... 3 I.4 Nagy teljesítményű számítási rendszerek ................................................................................. 4 II. Rekonstrukció ............................................................................................................................................ 5 II.1 A programról ....................................................................................................................................... 5 II.2 Az algoritmusról ................................................................................................................................. 5 III. Fejlesztői dokumentáció ....................................................................................................................... 7 III.1 Szoftverkörnyezet ............................................................................................................................ 7 III.2 A Program felépítése ....................................................................................................................... 7 III.2.1 HPC SDK ....................................................................................................................................... 8 III.2.2 MCP SDK....................................................................................................................................... 8 III.2.3 PETprocessing SDK.................................................................................................................. 8 III.2.4 Teratomo SDK ............................................................................................................................ 8 III.3 A program működése ..................................................................................................................... 8 III.4 Főbb osztályok................................................................................................................................... 9 III.4.1 AlgorithmSettings osztály ..................................................................................................... 9 III.4.2 DetectorEfficiency .................................................................................................................... 9 III.4.3 DetectorGeometry osztály .................................................................................................... 9 III.4.4 Array2D osztály ...................................................................................................................... 10 III.4.5 Volume3D osztály................................................................................................................... 10 III.4.6 KleinNishinaDistribution osztály ..................................................................................... 10 III.4.7 Timer osztály............................................................................................................................ 11 III.4.8 DataMemoryObject osztályok ........................................................................................... 11 4
III.4.9 LorSet osztály........................................................................................................................... 11 III.4.10 SourceDistribution osztály............................................................................................... 11 III.4.11 MaterialDistribution osztály ........................................................................................... 11 III.4.12 VoxelSourceList osztály .................................................................................................... 12 III.4.13 XcomDatabase ....................................................................................................................... 12 III.4.14 Reconstruction osztály ...................................................................................................... 12 III.5 Osztály kapcsolatok ....................................................................................................................... 12 III.6 Rekonstrukció folyamata ............................................................................................................ 13 III.7 GPU kód .............................................................................................................................................. 16 III.7.1 Vetítés ......................................................................................................................................... 17 III.7.2 Szűrések ..................................................................................................................................... 18 IV. Időfüggő rekonstrukció ....................................................................................................................... 19 IV.1 Kivitelezési mód .............................................................................................................................. 19 IV.2 Megvalósítás ..................................................................................................................................... 19 IV.3 Mérési adatok................................................................................................................................... 20 IV.4 Eredmények ...................................................................................................................................... 20 IV.4.1 Teljes adatsor ........................................................................................................................... 20 IV.4.2 Időszeletek rekonstrukciója ............................................................................................... 21 IV.4.3 Funkcionális elemzés ............................................................................................................ 22 V. Összefoglalás ............................................................................................................................................. 26 Idézett forrásmunkák ................................................................................................................................. 27
I. BEVEZETŐ A sugárzások detektálásán alapuló berendezések elengedhetetlenné váltak napjaink orvostudományában, elsősorban a diagnosztikában. Röntgen sugárzáson alapuló készülékek, mellyel csonttörést, tüdőgyulladást lehet diagnosztizálni, már a XX. század elejétől működtek és napjainkra rengeteget fejlődtek. A tomográfia fejlődésével pedig egészen új berendezések jelentek meg (CT, PET, MRI), melyek már sokkal pontosabb, akár 3 dimenziós képet is adhatnak a betegről, így mondhatni forradalmasította a tumoros megbetegedések és egyéb a berendezéssel kimutatható elváltozással járó betegségek gyors és fájdalommentes diagnosztikáját. Ezenkívül a PET és a SPECT lehetővé tette a molekuláris és a sejtszintű anyagcsere vizsgálatát. A képalkotás azonban a később ismertetett okoknál fogva számítás igényes, ezért új, nagyobb számítási kapacitással rendelkező eszközök bevonására van szükség. Ezen közül az egyik, az elsősorban a szórakoztató elektronikában használt grafikus processzorok, melyek nagy teljesítménye gyorsíthatja a képalkotási eljárásokat. Szakdolgozatom célja egy már régóta fejlesztés alatt álló videokártyás gyorsítást használó PET rekonstrukciós eljárás megismerése és továbbfejlesztése. I.1 PET A rövidítés a pozitron emissziós tomográfiát takarja, mely pozitronok annihilációjából származó gamma fotonok detektálásán alapszik. A β+ sugárzásból származó pozitron lelassulva, majd egy elektronnal egyesülve kettő, az elektron nyugalmi tömegének (511keV) megfelelő, egymással ellentétes irányban kilépő gamma fotonná alakul. A foton párt koincidenciába kapcsolt szcintillációs detektorokkal érzékeljük, körben a beteg körül akárcsak más tomográfiás eljárások során. Az együttes beütésből érzékelt egyenesekből (Line of response, LOR) a forrás helye megfelelő képalkotó eljárásokkal visszaállítható. Az eljáráshoz pozitron-bomló anyagot kell a szervezetbe juttatni. Ezek az anyagok mesterségesen felaktivált atomok, melyek ismert, élettanilag fontos molekulába beágyazhatóak (pl. cukrok, fehérjék). A radioaktív anyag, már a beadástól kezdve nyomon követhető, ezáltal alkalmas a PET, anyagcsere folyamatok tanulmányozására. 1
Megfelelő molekulába ágyazva, elérhető, hogy a kontrasztanyag egy bizonyos szervben, vagy egy rákos daganatban dúsuljon fel, ezáltal rákdiagnosztikában is fontos szerepe van. Ilyen célra használják például a víz protonokkal történő bombázása során keletkező 18F atomot, melyet cukormolekulába ágyaznak. A fotonokat szcintillációs kristálytűkkel detektáljuk, melyből több (pl. 16x16) csatlakozik 4 foto-elektron sokszorozóhoz, melyek együttes jelének arányából meghatározható, hogy melyik kristályban történt a beütés. Elméletben ez mind nagyon jól hangzik, azonban gyakorlatban koránt sem ilyen egyszerű. A betegnek nem lehet beadni akármekkora aktivitású anyagot (ált. 200400MBq), és napokig sem tarthat egy vizsgálat, ezért a mérés során nagyon kevés beütésből kell gazdálkodni, ami rontja a pontosságot és/vagy a felbontást. A dolgot súlyosbítja, hogy a gamma fotonok elnyelődhetnek, vagy szóródhatnak az emberi szervezetben. Mindkét eset a detektált párok számát csökkenti, ezért a megfelelő felbontású képalkotáshoz fejlett és nagy számítási kapacitást igénylő módszerek szükségesek. I.2 TOMOGRÁFIÁS ELJÁRÁSOK A PET, ahogy a neve is mutatja, tomográfiás képalkotó eljárás, ezért alapvetően a Radon-transzformációra és inverzére épül. Az előrevetítés során a céltárgy körül forgatva a detektort, vagy több detektort alkalmazva a körön vizsgáljuk az átmenő (CT), vagy a kijutó (PET) sugárzás intenzitását. A keletkező képet szinogramnak hívják (egy pont képe a Radon transzformáció után egy sinus-hullám lesz). A visszavetítés során ebből a szinogramból, visszatranszformálható, az eredeti forráseloszlás, például a szűrt visszavetítés eljárással. A PET esetében, azonban a már említett kis beütésszámok jelentősen rontják az eljárás során keletkező kép minőségét, továbbá az anyagban (a páciensben) történő gyengülés és egyéb nemlineáris hatások miatt korrekciók szükségesek. Egy a fizikai jelenségeket a Radon-transzformációnál
részletesebben
figyelembe
venni
képes
matematikai
keretben, a pácienst, vagy mintát és a detektorokat egy rendszerválasz-mátrixszal ( Ajk ) írjuk le, melynek celláiban található egy adott térfogatelemből (voxel, xk ) egy adott
2
detektorpárba (válaszegyenesek, LOR-ok, y j ) beérkező beütéseinek valószínűsége, melyet az alábbi egyenlet foglal össze [1].
y j Ajk xk
(1.1)
k
Még ha csak egyetlen detektorkört alkalmazunk és a páciens térfogatát sem osztjuk nagyon kis darabokra, a rendszermátrix mérete akkor is nagy méretű, ezáltal jelentősen megnöveli a számítás memóriaigényét. A rendszermátrixon belül továbbá figyelembe kell és lehet venni azt, hogy a pozitronnak lassulnia kell mielőtt egy elektronnal egyesülne, mely során a keletkezés helyétől elvándorol (pozitronvándorlás), a keletkező fotonok szóródását, az elnyelődését és a sugárzás gyengülését is. Szintén ide tartozik a detektorok érzékenysége és az is, hogy a foton a detektorkristályokon belül is szóródhat, nehezítve a beérkezés pontos helyének detektálását. A
rendszermátrix
ismeretében
alkalmazhatunk
iteratív
eljárásokat
a
képrekonstrukcióhoz, melyek azon alapulnak, hogy egy előre megbecsült képből iterációval, vagyis a kapott eredményeinket újra feldolgozva kaphatjuk meg a tényleges képet. A rendszermátrix előállítása történhet ismert forrás vizsgálatának segítségével is, de egyre gyakrabban alkalmazzák a Monte Carlo módszereket. I.3 MONTE CARLO MÓDSZEREK Ezen módszerek lényege, hogy valószínűségi eloszlásokkal leírható fizikai rendszereket modellezhetünk velük. Megfelelő eloszlású véletlen számokkal tetszőleges statisztikus jellegű fizikai jelenség modellezhető. Alapvető eljárás az analóg lejátszás, mely során a valóságos történéseket úgy követjük nyomon, hogy a valószínűségi eloszlással jellemezhető, elemi történéseket lejátsszuk. Esetünkben a teljes gamma foton transzportot, mely tartalmazza a fent említett szórási, elnyelési és gyengülési jelenségeket. A rendszer geometriai jellemzőit és a vizsgált térfogat (pl. páciens) sűrűségeloszlását pedig CT berendezés segítségével kaphatjuk meg, mely már számos berendezés esetén egybe van építve a PET-tel. A Monte Carlo módszereknek azonban nagy hátránya, hogy nagyon számítás igényesek, sok részecske „életét” kell lejátszanunk ahhoz, hogy megfelelő statisztikát és pontos 3
eredményt kapjunk. Ez normál egy-két processzoros rendszereken, a vizsgálandó térfogattól függően akár napokig is eltarthatna. I.4 NAGY TELJESÍTMÉNYŰ SZÁMÍTÁSI RENDSZEREK A napjainkban használt általános processzor mellett, mely a PC-inkben van, a tudományos világban számos nagyobb teljesítményű megoldás is kínálkozik. Ilyenek a különböző szuperszámítógépek, vagy pedig a hálózatban kötött kisebb teljesítményű gépekből épített klaszterek. Ezeket régóta használják különböző számításigényes feladatok elvégzésére, modellezésre stb. Méretük és fogyasztásuk miatt azonban nem alkalmazhatóak hatékonyan PET/SPECT berendezésekkel. A technikai fejlődés azonban számos jóval alkalmasabb eszköz használatát teszi lehetővé, melyeket ma még elsősorban a szórakoztató ipar használ, de hála a fejlesztéseknek már általános célokra is elérhetőek. Ezek az eszközök a számítógépbe szerelhető videokártyák, melyek grafikus processzora kiválóan alkalmas egyszerű számítások párhuzamos elvégzésére. Azonban az architektúra "egy utasítás sok adaton" jellege miatt korlátozottan használható, ugyanis minden olyan művelet mely nem illik be az architektúra által támasztott követelményekben felesleges órajel ciklusokat emészt fel. Ezért kell többek között kerülni az elágazásokat és a nem meghatározott méretű ciklusokat. Ha a megfelelő szabályokat betartjuk, akkor egyes algoritmusok akár nagyságrendekkel felgyorsulhatnak, mint ha hagyományos processzorokon futtattuk volna őket.
4
II. REKONSTRUKCIÓ II.1 A PROGRAMRÓL A PET rekonstrukciós szoftver egy koprodukciós program, melyet közösen fejleszt a BME NTI, a BME IIT, az SE Izotópdiagnosztikai és Onkoterápiás Klinikája és a Mediso Kft. A programot évek óta fejlesztik a két tanszéken különböző megközelítésű Monte Carlo algoritmusokkal, de mindkét kód az NVIDIA CUDA programfejlesztő környezetet használja. Ez egy C programnyelven alapuló fejlesztési környezet, mely az NVIDIA grafikus processzorok programozására alkalmas. Mivel elterjedt programozási nyelven alapszik, ezért viszonylag könnyen elsajátítható a használata, ha az ember megtanulja az alapvető architekturális jellemzőket és korlátokat. Az évek során a kód rengeteget fejlődött. A kezdeti proof of concept 1 megvalósításból mára egy, a mai divatos programozási metodikákat is követő, könnyen fejleszthető szoftvercsomaggá vált. A folyamat azonban nem volt egyszerű. Különböző okok miatt sem az eredeti fejlesztők, sem a nagyobb átírásokat végző fejlesztők nem dolgoznak már a projekten, ezért programozói tapasztalataimra építve az is a feladataim közé tartozott, hogy némiképp megértsem a jelenlegi kód főbb részeit és ezeket dokumentáljam is, későbbi felhasználás céljából. II.2 AZ ALGORITMUSRÓL A BME NTI által fejlesztett kód Monte Carlo alapú számításokat végez és a rekonstrukció során minden fizikai jelenséget igyekszik figyelembe venni, a pozitron vándorlástól, a Compton-szórásig. A vizsgált részecskék életét az aktivitást tartalmazó térfogattól a detektorokig követi, ellentétben a BME IIT által fejlesztett kóddal, mely a jelenséget visszafelé vizsgálja. Mivel a mérési adataink a detektorokban érkezett beütések és a koincidenciák alapján válaszegyesek, ezért az IIT-s innen indul ki és visszafelé vizsgálja az eseményeket a forrásig. Mindkét módszernek megvannak a maga előnyei és hátránya. Az NTI-s kód teljes mértékben 1
leírja
a
jelenségek
fizikáját,
azonban
több
helyen
adat
és
POC: a kutatás-fejlesztés legelső fázisa, melyben kiderül, hogy egy új ötlet megvalósítható-e vagy sem
5
számításigényesebb, mint másik kód mely viszont pont a fizikai jelenségek, mint például Compton-szórás miatt kénytelen mesterséges elhanyagolásokkal és korrekciókkal operálni a helyes eredmény érdekében. Az algoritmusnak lényegében a rendszermátrixos egyenlet (1.1) inverzét kell megoldania, de a mátrix a térfogatelemek és a válaszegyenesek számából következő dimenziója igen nagy ezért nem tároljuk magát mátrixot, hanem elemeit minden iterációban újraszámoljuk, az egyenlet inverzének megoldása helyett pedig az ML-EM2 iteratív eljárást alkalmazzuk.
Maximum Likelihood – Expectation Maximization: Legnagyobb valószínűségű becslés maximálás, részletesebben az Orvos-fizikus tankönyv (2.4.5) fejezetében [2] 2
6
III. FEJLESZTŐI DOKUMENTÁCIÓ III.1 SZOFTVERKÖRNYEZET A program futtatásához szükség van NVIDIA típusú videokártyára, de ha megfelelő beállításokkal fordítjuk le akkor anélkül is működőképes. Ekkor a működése természetesen sokkal lassabb, annak ellenére, hogy a felhasznált OpenMP könyvtár segítségével többprocesszoros rendszert is hatékonyan ki tud használni. Amennyiben videokártyára fordítjuk le a programot, szükség van a CUDA fordítóra, mely a videokártyás kódokat maga fordítja, míg a processzora írt kódok fordítását átadja a GNU C/C++ fordítónak. A program csomagot Linux operációs rendszerre fejlesztették, mely elsősorban a fordító rendszert érinti. A kezdetekben még volt törekvés arra, hogy platform független legyen, de a kód jelenlegi állapotában ez már csak nagyon nehezen kivitelezhető. A fordító rendszeren kívül a program pár kód ellenőrző eszközt is használ, melyet szintén Linux rendszerekre fejlesztettek. A kód igen nagy mérete miatt a teljes fordítás automatizált szkriptek segítségével zajlik, így új rendszerekre különösebb gond nélkül átvihető. A program elsődleges dokumentációja Doxygen-nel készül, mely a kódba írt, megfelelően
formázott
megjegyzésekből
képes
dokumentációt
generálni.
A
dokumentáció kevésbé kód specifikus része, pedig korábbi szakdolgozatokban és projektbeszámolókban található. A kód karbantartását és a változások követését Subversion verziókezelő rendszerrel végzik. III.2 A PROGRAM FELÉPÍTÉSE A program fejlesztése során a kód mennyiség nagy mértékben megnőtt. Mivel már nem elsősorban egyetemi fejlesztésről van szó, hanem egy konkrét cég végzi a kód karbantartását, ezért több dolgot az ő igényeiknek megfelelően kellett megalkotni, hogy a kód később is feldolgozható és fejleszthető maradjon. Ehhez szükségessé váltak olyan segédcsomagok is, melyeknek nincs algoritmikus szerepe a rekonstrukcióban.
7
III.2.1 HPC SDK3 Az általánosítás és a szállíthatóság érdekében volt törekvés arra, hogy kód számos részét absztrahálják és egy olyan fejlesztési eszközkészletet hozzanak létre, melyre aztán számos egyéb nagy számítási kapacitású rendszerre készített rekonstrukciós eljárás is építhetett volna. Azonban ez a csomag túlságosan általános és a hagyományos programozástól elrugaszkodott szemléletű ezért nem lett kihasználva a benne rejlő potenciál. Most pedig törekvés van a csomag teljes eltávolítására. III.2.2 MCP SDK Általános célú programelemek gyűjteménye. Ebből a csomagból a rekonstrukció csak nagyon kevés elemet használ. Ilyenek például a naplózás és a standard kimenet formázása. III.2.3 PETPROCESSING SDK A HPC SDK-hoz hasonlóan elsősorban PET berendezések geometriáját és tulajdonságait leíró absztrakt objektumok, valamint adatfeldolgozást segítő rutinok gyűjteménye. III.2.4 TERATOMO SDK Ez a csomag a tényleges rekonstrukciós kódokat. Mind a BME IIT, mind a BME NTI által fejlesztett kódot. A továbbiakban ezért ezzel a fogok részletesebben foglalkozni, ezen belül is az NTI által fejlesztett kóddal. III.3 A PROGRAM MŰKÖDÉSE A program rekonstrukciós motorjához lényegében külső egységként kapcsolódik a felhasználó számára elérhető alkalmazás. A program elindítható a saját parancssoros alkalmazással is, de a konfigurálás könnyebb elvégzéshez rendelkezésre állnak külső szkriptek, mely a beállítások nagy részét automatikusan elvégzik. Erre elsősorban szoftver fejlesztés, és a könnyebb felhasználhatóság érdekében van szükség. Később ugyanis olyan szintre kell fejlődnie ennek, hogy a rendszert használó orvosok is kényelmesen használhassák a programot. 3
SDK: Software Development Kit - szoftverfejlesztő eszközkészlet/csomag
8
Jelen állapotában azonban a programnak számos paraméterét kell konfigurálni egy-egy rekonstrukció előtt. Ilyenek például a geometriai jellemzők, a detektor tulajdonság, a mérési adatok tulajdonságai, az anyag jellemzők, vagy például a lehetséges kimeneti adatok. A program többféle kimenettel rendelkezik. Egyrészt a standard kimeneten folyamatosan tájékoztat az aktuális rekonstrukciós folyamat állapotáról, másrészt többféle szöveges és bináris alapú módon kaphatunk belőle adatokat. A kiértékeléshez legtöbbször három dimenziós adattömböket választunk kimenetnek, melyet megfelelő képszerkesztő szoftverrel vizsgálhatunk. III.4 FŐBB OSZTÁLYOK A program belépési pontja a main() függvény, mely elvégzi a parancssori paraméterek és a felhasználó által megadott konfigurációs fájlok kiértékelést és a rekonstrukció előkészítését. Itt kerül ellenőrzésre maga rendszer is, hogy alkalmas-e a kód futtatására. A konfiguráció után következik a megadott adatok beolvasása. Ilyenek lehetnek a mérési adatok, kristály hatásfok együtthatók, hatáskeresztmetszetek. Az összes beállítás és adatbeolvasás után indulhat a rekonstrukció. III.4.1 ALGORITHMSETTINGS OSZTÁLY Ez az osztály tartalmazza és kezeli a program összes megadható paraméterét. Minden paraméterhez tartoznak beállító, kiolvasó és ellenőrző metódusok. A paramétereket vagy a parancssori paraméterekből, vagy pedig a konfigurációs fájlban megadott értékekből állítja elő. A programban számos más osztály is támaszkodik rá. III.4.2 DETECTOREFFICIENCY Ez az osztály tartalmazza a detektor válaszfüggvényhez szükséges adatokat. III.4.3 DETECTORGEOMETRY OSZTÁLY Ez az osztály tartalmazza a detektorok geometriai tulajdonságait:
detektor modulok száma
detektorkristályok száma modulonként
detektorok fizikai mérete 9
Ezekből az adatok aztán kiszámolható a LOR-ok száma, az aktív térfogat és egyéb geometriai adatok. Az osztály erre is tartalmaz metódusokat. Tervbe van véve, hogy ezzel az osztállyal tetszőleges detektor megadható legyen, de egyelőre három előre megadott detektor adatait tartalmazza a program, melyek közül a NanoPET/CT [3] rendszerrel foglalkoztam behatóbban, melyet pre-klinikai mérésekhez használnak. III.4.4 ARRAY2D OSZTÁLY Minta osztály két dimenziós tömbök kezelésére, mely a tényleges, tetszőleges típusú adatot egy dimenziós tömbben tárolja. Tartalmaz metódusokat a memóriafoglalásra, alapvető aritmetikai műveletek elvégzésre és a 2D-1D indexelés számolására. III.4.5 VOLUME3D OSZTÁLY Ez az osztály szerkezetében és metódusaiban szinte teljesen azonos az Array2D osztállyal, csak három dimenzióban. A térfogati forrás- és anyageloszlást leíró osztályok őse. III.4.6 KLEINNISHINADISTRIBUTION OSZTÁLY Ez az osztály gondoskodik a Compton-szórás modellezéséről, azáltal, hogy a KleinNishina formula alapján kiszámolja a részecske energiáját és eltérülési szögét a szóródás után. A számolás a Kahn algoritmussal történik, de mivel az több elágazást is tartalmaz, nem teszi alkalmassá arra, hogy párhuzamos számítást végző algoritmusban használjuk, ezért nem tehető bele a videokártyás kódba. Ezért egy táblázatba a kód előre kiszámolja a szórási szögeket és energiákat. Az értékeket a program a gyors elérés érdekében textúra memóriába tárolja. Ahhoz, hogy az adatok felbontása minél inkább közelítse magát a Kahn algoritmussal számolható értékeket, két szintű táblázat készül. Az első viszonylag durvább felbontással tartalmazza a teljes szögtartományt, míg a második finomítva tartalmazza a kisebb szögtartományokat.
10
III.4.7 TIMER OSZTÁLY A program futásidejének mérésére szolgáló osztály, mely az időt az eltelt órajelek alapján számolja. III.4.8 DATAMEMORYOBJECT OSZTÁLYOK Egy, kettő és három dimenziós adat objektumokat leíró osztályok. Segítségükkel ugyanaz az adatmennyiség lefoglalható a rendszer és a videokártya memóriájában is, valamint a kettő között másolható. Ezek is mintaosztályok, így tetszőleges típusú adattal megtölthetőek. Ezek az osztályok kétféleképpen működhetnek: példányosításkor vagy mindkét memóriában lefoglalódik a hely az adat számára, vagy csak a videokártya memóriában. A Tex végződésű osztályok az adatokat a videokártyán nem a globális memóriában tárolják, hanem a textúráknak fenntartott memóriában. Az itt tárolt adatok bizonyos műveletek során gyorsabban elérhetőek. A Const végződésű osztály a futás során konstansnak szánt adatok tárolója. A kódban való könnyebb azonosítás érdekében ezek az adatok névvel is elláthatók. III.4.9 LORSET OSZTÁLY A LorSet osztály írja le a bemeneti adatot, vagyis a detektorok koincidenciájából kapott egyeneseket (LOR). Az adatsort a rendszermemóriában foglaljuk le, de ahhoz, hogy grafikus processzorral feldolgozzuk a videokártya dedikált memóriájába kell másolni, ezért az osztály mindkét memória kezelésére tartalmaz metódusokat. III.4.10 SOURCEDISTRIBUTION OSZTÁLY Ez osztály reprezentálja a térfogati forráseloszlást, mellyel az eredményt és számos köztes térfogatot is leírunk. III.4.11 MATERIALDISTRIBUTION OSZTÁLY Ez az osztály írja le az anyag eloszlását a térfogatban, melynek a helyes szórásszámításban van szerepe. 11
III.4.12 VOXELSOURCELIST OSZTÁLY Ez az osztály írja le az egyes voxelek forráseloszlását, vagyis, hogy mennyi pozitron található egy térfogategységben. A VoxelSourceExList osztály lényegében ugyanezt csinálja, de ez az osztály a térfogatot kisebb részekre osztva tartalmazza. Ez biztosítja, hogy a grafikus processzoron futó algoritmus kisebb, gyorsabban feldolgozható adatcsomagokra osztva kapja meg a teljes térfogatot. III.4.13 XCOMDATABASE Ez az osztály végzi az egyes anyagok hatáskeresztmetszetének beolvasását és tárolását. A hatáskeresztmetszet adatokat a NIST XCOM programmal állítjuk elő. Egy-egy anyag adatait az XcomDatabaseMaterial osztály példányaiban tároljuk és ezeket fogja össze az XcomDatabase osztály. III.4.14 RECONSTRUCTION OSZTÁLY A Reconstruction osztály tartalmazza a rekonstrukcióhoz szükséges adatokat és metódusokat. Legfőbb metódusa a doReconstrucion(), a rekonstrukció belépési pontja, mely az megadott bemeneti adatokból (LOR) kiszámolja a térfogat aktivitás eloszlását. Ha referencia bemenet is meg van adva egy korábbi kiszámolt térfogat, vagy matematikai fantom formájában, akkor minden iteráció után kiszámolja a különböző normákat, mellyel az új eredmény összehasonlítható. A rekonstrukció folyamatáról részletesebben a későbbiekben lesz szó. III.5 OSZTÁLY KAPCSOLATOK A lentebb látható osztálykapcsolati diagramon azt próbáltam szemléltetni, hogy az egyes osztályok milyen adatkapcsolatban állnak egymással. A diagram az UML (Unified Modelling Language - Egyesített modell nyelv) szerint készült, mely az objektum orientált programozás általános ábrázoló és tervező eszközévé vált. Az ábrán üres nyíllal az osztályok közötti leszármazást jelöljük. Egyszerű vonallal az általános kapcsolat jelöljük, melyhez szoktak számosságot is jelölni. Ezzel lehet például azt ábrázolni, hogy ha egy osztály több, vagy csak egy példányt tartalmaz egy másik
12
osztályból. Az UML ezen kívül több más kapcsolatot is definiál az osztálydiagramon, de ezekre most nincs szükség.
T : typename Volume3D
<
> srclst
T : typename Array2D
1..*
1
T : typename DataMemoryObject3DTex
VoxelSourceExList
1
T : typename LorSet
1 MaterialDistribution
T : typename DataMemoryObject1D
1
SourceDistribution
VoxelSourceList
1
1
4
1
1..* <> VoxelSource
T : typename DataMemoryObject3D
ReconstructionNTI
AlgorithmSettingsNTI 1
7
T : typename DataMemoryObjectConst
1 DetectorGeometryNTI
1
1 DetectorEfficiencyNTI
1
1
KleinNishinaDistribution
2 T : typename DataMemoryObject2DTex 1
XcomDatabaseMaterial
1 XcomDatabaseNTI
<>
1. ÁBRA OSZTÁLYKAPCSOLATI DIAGRAM
III.6 REKONSTRUKCIÓ FOLYAMATA A rekonstrukció adatmentéssel kezdődik. Amennyiben a konfigurációban megadtunk szinogram vagy detektorkép kimenetet, akkor ezeket a képeket a program itt számolja ki a LOR-okból. A szinogram a Radon-transzformáció eredménye, míg a detektorkép pedig a detektorokban tapasztalt beütéseket ábrázolja kiterítve.
13
2. ÁBRA KITERÍTETT DETEKTORKÉP
3. ÁBRA SZINOGRAM
Ezután következik a forrás- és anyageloszlás előkészítése és az esetleg kezdő térfogati eloszlás vagy referencia eloszlás beolvasása és ez alapján az L2 norma kiszámolása, mely jól jellemzi a referencia és a rekonstruált adat közti kvantitatív különbségeket. Az L2 normát a következőképpen számoljuk, ahol Ai és Bi az egyes voxelek aktivitása, n pedig a voxelek száma [1]: n
normL 2
A B i 1
i
2
i
(4.1)
n
Ezután következik a LOR-ok feldolgozása. Amennyiben egy válaszegyenesnek sok olyan szomszédja van ahol nem volt beütés, akkor feltételezhető, hogy zaj, ezért ezeket a LORokat kiszűrjük. Ha nem volt megadva kezdeti térfogat eloszlás akkor a rekonstrukció a visszavetítéssel kezdődik, melynek során elkészül az első térfogat rekonstrukció a LOR-ok alapján. Ha volt megadva kezdeti érték akkor kezdődhet az iteráció. Az iteráció során először végrehajtjuk a térfogati eloszlásra az előrevetítést. Az előrevetítés során a térfogati eloszlásból kiszámoljuk a detektorválaszt. Ez az alábbi képlettel írható le legegyszerűbben, melyben y j a detektorpárokban detektált fotonok száma, xk a voxelek aktivitása és Ajk a rendszermátrix [1]. 14
y j Ajk xk
(4.2)
k
Ezek után ismét elkészülnek a detektorkép, szinogram és MLOR kimenetek, melyek kiírása után következik a visszavetítés, mely lényegében az előrevetítés műveletét használja, de segítségével korrekciós tényezőket számolunk ki a következő módon, ahol ck jelöli a tényezőket, pedig a detektorpárhoz tartozó aránytényező [1].
xk n
ck
n1
xk ck
1 Ajk
n
A j
n
y jn mért jk
y jn számolt
(4.3) (4.4)
j
A kapott adatot, vagyis a térfogati aktivitás eloszlást különböző módokon szűrjük (medián szűrés, tüskék eltávolítása). Végül a térfogat eloszlás kimenetei következnek. Amennyiben a konfiguráció során megadtuk, a program képes az adatokat szöveges MVOL fájlban, vagy bináris RAW fájlban kiírni, vagy pedig a középső szeletet PNG képfájlban. A kimeneti adatok után újabb iterációs lépés következik.
15
[ kimenet ] LOR kimenetek
forrás- és anyageloszlás elokészítése szomszédvizsgálat
elorevetítés [ kezdeti térfogat ] [ iteráció vége ]
elorevetítés
LOR kimenetek [ kimenet ] visszavetítés
szurések térfogat kimenetek [ kimenet ]
4. ÁBRA REKONSTRUKCIÓ FOLYAMATA
III.7 GPU KÓD Ahogy fentebb már említettem a grafikus processzorra nem lehet hagyományos módon programozni. Ezek a processzorok arra vannak kitalálva, hogy több adaton végezzék el ugyanazt a műveletet párhuzamosan. A mai videokártyák már számos olyan általános programozásban használt elemet támogatnak, melyeket a 10 évvel ezelőtti modellek még nem is ismertek, így ma már sokkal egyszerűbb ezekre az eszközökre programozni, csak be kell tartani néhány szabályt és ismerni kell a hardver korlátait. A legfontosabb, hogy a programban kerülni kell az elágazásokat és a nem meghatározott méretű ciklusokat. A párhuzamosan futó szálak ugyanis akkor működnek hatékonyan, ha 16
mindegyik futásideje azonos, szinkronban dolgoznak. Egy elágazás elronthatja a szinkront és a lassabb szálak feleslegesen lassítják a rendszert, ugyanis bizonyos megadott és automatikusan generált pontokon az egyes szálak bevárják egymást. Ezen kívül figyelnünk kell arra, hogy a videokártyák jóval korlátozottabb mennyiségű mint a rendszermemória, de gyorsabb is. A mai csúcskártyákon 3-4GB dedikált memória található, míg rendszermemóriából már egy egyszerűbb alaplapba is tehetünk akár 16GB-ot. Ha az algoritmus jó és hatékonyan párhuzamosítható, akkor még mindig figyelnünk kell arra, hogy nehogy memórialimitált legyen az alkalmazásunk. Ha sokkal több adatot kell feldolgoznunk mint amennyi elfér a videomemóriába akkor annak a rendszer és a videokártya közti másolgatás lesz az eredménye, mely sebesség szempontjából szintén szűk keresztmetszett. III.7.1 VETÍTÉS Az előrevetítés és a visszavetítés műveletének hasonlósága miatt a GPU-n egyetlen függvényt definiálunk erre célra. Ez a függvény szálanként végzi az egyes foton párok életútjának követését. Hogy a fentebb említett szinkronizálást fenntartsuk olyan Monte Carlo eljárásokat alkalmazunk, melyek során a fotonok nem halhatnak meg. Ez ugyanis egy szál idő előtti végét jelentené, ezen túl pedig egy fölösleges számítást, mely 2∙1010 indított részecske esetén igen nagy mennyiségű elpazarolt órajel ciklus. Az első lépésben sorsolunk egy haladási irányt a foton pár egyik fotonjának. A párja ennek az iránynak az ellentettjét kapja. Ezután a Woodcock módszerrel kisorsoljuk a foton szabad úthosszát. A szabad úthossz végén mintavételezünk egy új irányt és energiát a Klein Nishina formula alapján készített táblákból. Ezt egy ciklusban végezzük addig míg a foton el nem hagyta a vizsgált térfogatot, vagy el nem érte a detektort. Amennyiben a eléri a detektort a fotonpár másik felére is végrehajtjuk a műveletsorozatot. A vetítés egy GPU-n futó függvény, melyhez a processzor számára több befoglaló függvény is van. Egyrészt az mely eldönti, hogy előre vagy visszavetítésről van szó, másrészt az amelyik a különböző bemeneti paramétereket alapján felparaméterezi a GPU-s függvényt is. Mivel a GPU-s függvény egy mintafüggvény ezért fordításkor
17
minden a futás során használt példánya külön fordul bele a kódba, ezzel is elősegítve a fordító program számára a hatékonyabb optimalizálást. III.7.2 SZŰRÉSEK A vetítésen kívül az LOR-ok szűrése az amit lehet GPU-val gyorsítani, hiszen jól behatárolható adatsoron végez azonos műveleteket. A LOR-ok szűrésére az egyik alapvető fontosságú lépés a már fentebb is említett szomszédvizsgálat. Ezen kívül alkalmazhatunk még medián szűrést is a LOR-okra, mind az előrevetítés előtt mind utána.
18
IV. IDŐFÜGGŐ REKONSTRUKCIÓ A program megismerése és dokumentálása után feladatom az időfüggő rekonstrukció kipróbálása volt és a programba építési lehetőségének előkészítése. Az időfüggő rekonstrukció fontos, mert a PET eljárást alapvetően funkcionális tomográfiás eljárásnak alkották meg, vagyis, hogy időbeli folyamatok is követhetőek legyenek rajta. Így például a kontrasztanyag felszívódásának, vagy egyes szervekben történő felhalmozódásának vizsgálata. IV.1 KIVITELEZÉSI MÓD Erre a legkézenfekvőbb megoldás az algoritmus gyorsítás, hogy már vizsgálat után minél előbb látható képet kapjunk és a testben lezajló folyamatokat is követhessük, de a jelenlegi technológián a futásidő nem csökkenthető minden határon túl. Ezen túl a kis aktivitás is jelentős szerepet játszik. Mikor a kontraszt anyagot beadják, a felszívódásig jó eséllyel nagyon kevés beütésből kell jó rekonstrukciót alkotni, mely idő és a számítás igényes. Ezért első körben egy megoldással foglalkoztunk, mely javíthatja az utólagos rekonstrukciót. A módszer lényege, hogy a rekonstrukciót először elvégezzük a teljes adatsoron, vagyis a teljes vizsgálat eredményét feldolgozzuk. Ezek után az ebből a rekonstrukcióból származó eredmény bemeneti adatként használjuk az egyes időszeletek rekonstrukciójához. Ily módon elkerülhetjük azt, hogy főleg a kezdeti alacsony beütésszámok esetén ne tudjunk látható képet rekonstruálni,
ugyanis az végeredmény már megvan az
algoritmusnak csak az a dolga, hogy az iterációk során eltávolítsa azokat a részleteket, melyek nem abban az időszeletben keletkeztek. Ezzel az eljárással már a korai időszeletekben is értelmezhető képet kaphatunk. IV.2 MEGVALÓSÍTÁS Köszönhetően a program objektum orientált fejlesztési szemléletének és annak, hogy a mérési eredmények azonos szerkezetű, tömörítetlen szövegfájlokban tárolhatók, nagyon könnyű az amúgy időszeletekre osztott mérési eredményeket ismét összegezni. 19
Ezt egyelőre a programon kívül végezzük, egy nagyon egyszerű alkalmazással. Az összegzett adatsorra ezután egy viszonylag hosszabb és részletekre törekvő konfigurációval futtatjuk az algoritmust. A kapott eredmény ezután bemeneti adatként megadható az egyes időszeletek rekonstrukciójához. Az egyes időszeletek rekonstrukciója adatmennyiség hiányában nem igényel annyi iterációt mint a teljes adatsor, ahhoz hogy értékelhető eredményt kapjunk. IV.3 MÉRÉSI ADATOK Az ötlet teszteléséhez a Mediso Kft-től kaptunk tényleges mérési adatokat. A mérés egy egér több mint egy órás vizsgálatából származik. A tényleges 3780 másodperces mérési időt 24 különböző hosszúságú időszeletre osztották. Az egérnek FDG kontrasztanyagot adtak be, mely
18F
pozitron bomló izotópot tartalmazó cukorszerű vegyület és rögtön
indították a mérést. IV.4 EREDMÉNYEK IV.4.1 TELJES ADATSOR Először vizsgáljuk meg a teljes adatsorra készült rekonstrukciót. A vizsgálandó térfogatot 256×256×256 felbontással rekonstruáltam. A térfogat valós mérete 8,8×8,8×9,4 cm3 volt, mely lényegében a NanoPET belső mérete [3]. Ahogy az alábbi képeken látható a kontrasztanyag leginkább a vesékében és a húgyhólyagban dúsult fel, nagy mennyiségben megtalálható a szívben és az agyban is.
20
gerincoszlop
vese
agy
szív 5. ÁBRA OLDALNÉZETŰ KERESZTMETSZET
vese
gerincoszlop
agy 6. ÁBRA FELÜLNÉZETI KERESZTMETSZET
IV.4.2 IDŐSZELETEK REKONSTRUKCIÓJA Összehasonlításként kiválasztottunk az azonos időszeletekből olyan rekonstrukciót, mely rendelkezett bemeneti térfogateloszlással és olyat is mely nem, hogy valóban látható-e különbség a két eljárás között. A szemléltetéshez a térfogat egy olyan szeletét 21
választottuk, melyen jól látható az egér agya. Ehhez a térfogatból a 210. szeletet választottam. A képen látható polipszerű karok, a rekonstrukció során keletkeznek a random koincidenciák szűrésének hiánya miatt. A random koincidenciákat olyan beütések, melyek nem egy fotonpárból érkeztek, hanem egyszerűen csak egy időben érkeztek, de különböző párkeltésből. Ilyen gyakran előfordul ha a detektorgyűrűn kívül is van aktivitás, mely az árnyékolások ellenére is bejuthat a detektorokba. Ez a jelenség nagy valószínűséggel el fog tűnni, ha a random szűrést magába az előrevetítési lépésbe építjük bele.
7. ÁBRA EGÉRAGY - TELJES MÉRÉS 210. SZELET
A rekonstrukciót tehát több időszeletre lefuttattam úgy is, hogy nem adtam meg az iterációnak kezdő térfogatot. Ekkor azonban több próbálkozás és paraméter állítás ellenére sem sikerült értékelhető eredményt kapni, így nyugodtan mondhatjuk azt, hogy az ötlet sikeres volt és megéri az extra számolást választani, a jó minőségű rekonstrukcióhoz. IV.4.3 FUNKCIONÁLIS ELEMZÉS Mivel a PET funkcionális vizsgálati módszer, ezért az eredményeket ilyen szempontból is megvizsgáltam, hogy végigkövethető-e egy folyamat az egyes időszeletek rekonstrukciója alapján. Az adatsorokkal együtt azt az információt kaptam, hogy általában az egér agyában az aktivitás először lassabban halmozódik fel, majd lassan cseng le, ellentétben mondjuk a kiválasztó szervrendszerével (vese, húgyhólyag), 22
melyben már a kezdetektől fogva nagy aktivitás egyenletesen bomlik el. Ezért vizsgáltam az agya és a vesék körüli térfogatot számszerűen is. IV.4.3.1
AGY
Az alábbi képeken jól látható az aktivitás csökkenése. A rekonstruált térfogatból ismét a fentebb is választott 210. szeletet választottam. A 13., 15., 18. és 24. időszelet rendre a 1080s, 1440s, 2340s, 3780s eltelt időt jelentik.
8. ÁBRA EGÉR AGY - 13,15,18,24 IDŐSZELET
Hogy az adatokat számszerűen is feldolgozzam, a kapott térfogatból kiválasztottam a 200-232 szelet közti térfogatot és összegeztem benne a pozitron bomlások számát, majd ezt leosztva az egyes időszeletek méretével megkaptam az aktivitás időbeli változását a mérés során, melyet az alábbi grafikonon tüntettem fel. A térfogati eloszlásból a képeket az ImageJ programmal szerkesztettem és jelenítettem meg. A grafikon Origin szoftverrel készült.
9. ÁBRA EGÉR AGY - ÖSSZEGZÉSI TÉRFOGAT
23
agyi aktivitás
aktivitás [Bq]
1,50E+015
1,00E+015
5,00E+014
0,00E+000 0
500
1000
1500
2000
2500
3000
3500
4000
eltelt idő [s]
10. ÁBRA EGÉR AGY - AKTIVITÁS
IV.4.3.2
EGÉR KIVÁLASZTÓ SZERVEK
A kiválasztó szervrendszerhez tartozik a vese és a húgyhólyag, melyből elsősorban a mérésből a vesék látszanak jól. A kiértékeléshez az előzőekhez hasonlóan jártam el. A térfogatból a 0-28 szeleteket összegeztem.
11. ÁBRA EGÉR VESE - ÖSSZEGZÉSI TÉRFOGAT
24
vesék aktivitása 2,00E+016
aktivitás [Bq]
1,50E+016
1,00E+016
5,00E+015
0
500
1000
1500
2000
2500
3000
3500
4000
eltelt idő [s]
12. ÁBRA EGÉR VESÉK – AKTIVITÁS
A grafikonon jól látható, hogy a vesékben jól hamarabb felhalmozódik az aktivitás. A rekonstrukció során is ez a terület volt mindig a legfényesebb.
25
V. ÖSSZEFOGLALÁS A szakdolgozatom elkészítése során volt alkalmam bekapcsolódni egy PET rekonstrukciós program fejlesztésébe. Megismerhettem a keretrendszereket melyek rekonstruáló motort körbeveszik és alkalmam volt komolyabban foglalkozni a kód elemzésével. A kód megismerése után a funkcionális elemzéssel foglalkoztam. Ennek során egy egérről készült PET felvételt értékeltem ki. A mérés 3780 másodpercig tartott, de az adatokat eltérő hosszúságú időszeletekre osztották. A mérés elején sűrűbben, a mérés vége felé már csak 6 percenként rögzítették az adatokat. A kiértékelés lényege a funkcionális elemzés volt, vagyis, hogy időben hogy változott az aktivitás az egér szervezetében. Mivel a mérési adatok, elsősorban mérés elején még nagyon kevés beütés tartalmaznak, ezért a rekonstrukció nagyon nehéz. Mint kiderült a program jelenlegi állapotában nem is alkalmas ezen adatsorok kiértékelésére. Ezért választottam azt a megoldást, hogy először a teljes mérést értékeltem ki és ennek eredményét az iteráció kezdeti értékeként használva rekonstruáltam a kisebb adatsorokat. Ennek segítségével sikerült időbeli rekonstrukciót készíteni, mely már alkalmas funkcionális elemzése.
26
IDÉZETT FORRÁSMUNKÁK [1] Cserkaszky, Á. (2011). Monte-Carlo based PET image reconstruction on GPU. Master's Thesis . Budapest, Hungary. [2] Kári, B., Kinga, K., Légrády, D., Bérczi, V., & Czifrus, S. (2012). Tankönyv Fizikusoknak. Forrás: Elektronikus oktatási anyag kialakítása az élő szervezet strukturális összetevőinek
és
biokémiai
folyamatainak
képalkotó
elemzésére:
http://oftankonyv.reak.bme.hu [3] Mediso. (2012). NanoPET/CT product page. Forrás: Mediso Medical Imaging Systems: http://www.mediso.hu/products.php?type=hw&fid=2&pid=55
27
