MTA DOKTORA PÁLYÁZAT DOKTORI ÉRTEKEZÉS KÉPALKOTÓ ÉS KÉPFELDOLGOZÓ ELJÁRÁSOK HATÉKONYSÁGÁNAK NÖVELÉSE AZ ORVOSINFORMATIKÁBAN

dc_592_12

MTA DOKTORA PÁLYÁZAT DOKTORI ÉRTEKEZÉS

KÉPALKOTÓ ÉS KÉPFELDOLGOZÓ ELJÁRÁSOK HATÉKONYSÁGÁNAK NÖVELÉSE AZ ORVOSINFORMATIKÁBAN

Benyó Balázs a műszaki tudomány kandidátusa

Budapest 2013.

dc_592_12

dc_592_12

Tartalomjegyzék

Bevezetés

1

Célkitűzések

3

1. Párhuzamos SPECT rekonstrukció

5

1.1. Nukleáris képalkotó módszerek alkalmazása az orvosi gyakorlatban . . . . .

6

1.1.1. Egy-fotonos számítógépes tomográfia (SPECT) . . . . . . . . . . .

7

1.2. Párhuzamos végrehajtást lehetővé tevő módszerek és platformok . . . . . .

9

1.2.1. Graphical Processing Unit (GPU) . . . . . . . . . . . . . . . . . . .

10

1.3. Probléma ismertetése . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

1.4. Megoldási módszerek . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

1.4.1. Iteratív képrekonstrukció . . . . . . . . . . . . . . . . . . . . . . . .

17

1.4.2. Detektor távolságfüggő felbontásának kompenzációja . . . . . . . .

24

1.4.3. GPU-ra szabott rekonstrukciós algoritmus . . . . . . . . . . . . . .

26

1.5. Az elért eredmények és azok értékelése . . . . . . . . . . . . . . . . . . . .

32

1.6. Összefoglalás . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

1.7. A fejezetben ismertetett eredmények tézisszerű megfogalmazása . . . . . .

52

1.7.1. A fejezet eredményeit bemutató publikációk . . . . . . . . . . . . .

52

2. Foggyökér csatorna meghatározása CT képeken

53


54

2.2. Megoldási módszerek . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

56

2.2.1. 3D objektumok térbeli vonalvázának meghatározása . . . . . . . . .

56

i

dc_592_12 ii 2.2.2. Micro-CT felvételek feldolgozása . . . . . . . . . . . . . . . . . . . .

57

2.2.3. CBCT felvételek feldolgozása . . . . . . . . . . . . . . . . . . . . .

65

2.3. Eredmények . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

2.3.1. Micro-CT felvételek feldolgozása . . . . . . . . . . . . . . . . . . . .

69

2.3.2. Cone beam CT felvételek feldolgozása . . . . . . . . . . . . . . . . .

71

2.4. Az elért eredmények értékelése . . . . . . . . . . . . . . . . . . . . . . . . .

72

2.5. Összefoglalás . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

77


78


78

3. Agyi régiók elkülönítése MR felvételeken

79


80

3.2. Alkalmazott módszerek . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

82

3.2.1. Az FCM algoritmus . . . . . . . . . . . . . . . . . . . . . . . . . . .

82

3.2.2. Inhomogén zaj kompenzálására szolgáló modellek . . . . . . . . . .

83

3.2.3. Az FCM algoritmus alkalmazása inhomogén zaj becslésére . . . . .

84

3.2.4. Gyors osztályozás hisztogram alapján . . . . . . . . . . . . . . . . .

84

3.2.5. A javasolt gyors osztályozó algoritmus . . . . . . . . . . . . . . . .

85

3.3. Az elért eredmények és azok értékelése . . . . . . . . . . . . . . . . . . . .

88

3.4. Összefoglalás . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

94


95


95

4. Eredmények hasznosítása

97

4.1. A kidolgozott képrekonstrukciós algoritmus hasznosítása . . . . . . . . . .

98

4.2. Foggyökér csatorna felismerésére kidolgozott módszerek hasznosítása . . . .

99

4.3. Agyi területek szegmentálására kidolgozott algoritmus hasznosítása . . . . 100 Irodalomjegyzék

i

dc_592_12

Bevezetés Az orvosbiológia és mérnöki tudományok határterülete egyike a legdinamikusabban fejlődő tudományterületeknek. Ennek az interdiszciplináris tudományterületnek a gyors fejlődését több tényező is motiválja. Egyfelől a tudományterülettel határos mérnöki területek (informatika, elektronika, automatizálás, nanotechnológia stb.) fejlődése olyan új technológiák és eszközök megjelenéséhez vezetett, amelyek lehetőséget teremtettek új orvosi diagnosztikai és terápiás eljárások kidolgozására, illetve a meglevő eljárások hatékonyabb változatainak kifejlesztésére. Másfelől a népesség mai társadalmakban tapasztalható elöregedése egyre növekvő követelményeket támaszt az egészségügyi ellátórendszerrel szemben, ami gyorsabb, olcsóbb, kevesebb emberi munkát igénylő, vagyis hatékonyabb orvosi diagnosztikai és terápiás módszerek kidolgozását követeli meg. Ezen egymást erősítő hatások eredményezték az interdiszciplináris tudományterület napjainkban tapasztalható gyors fejlődését mind tudományos, mind gazdasági szempontból. A hatékonyabb diagnosztikához, illetve terápiához szükséges új módszerek kidolgozása különös kihívást jelent a területen dolgozó kutatók számára több okból is. Az első kihívás, amellyel leggyakrabban kell szembenézni a kutatóknak a biológiai rendszerek modellezése során jelentkezik. A megfigyelt élettani, biológiai rendszerek az átlagos mérnöki rendszerekkel összehasonlítva lényegesen nagyobb bonyolultságú, összetett rendszerek, így azok leírása, modellezése általában nem lehetséges a mérnöki gyakorlatban szokásos módszerek közvetlen alkalmazásával. Az élettani, biológiai rendszerek kezeléséhez a fiziológiai folyamatok pontos megértése szükséges annak érdekében, hogy azoknak a megfigyelés szempontjából lényeges jellemzőit ki tudjuk emelni. A kutatóknak tehát nemcsak együtt kell dolgozni az interdiszciplináris terület más határterületén dolgozókkal, de azok tudásának java részét is el kell sajátítaniuk a hatékony diagnosztikai módszerek kidolgozásához szükséges modellek megalkotásához. A másik gyakran előforduló kihívás az orvosbiológiai mérnöki kutatások során, hogy a kidolgozott módszerek és algoritmusok rendkívül nagy számításigényűek. Ez természetesen más tudományterületeken is számtalanszor előforduló probléma. A különbség az, hogy az orvosi gyakorlat igen szigorú követelményeket támaszt egy-egy kezelés vagy diagnosz1

dc_592_12 2 tikai művelet végrehajtási idejére vonatkozóan, vagyis a kidolgozott algoritmusok végrehajtási idejének csökkentése ezen a területen igen fontos gyakorlati jelentőséggel bír. Azaz, amennyiben gyakorlatban felhasználható módszereket kívánunk kidolgozni, nem kerülhető el az algoritmusok számítási igényének olyan mértékben történő csökkentése, hogy azok futási ideje összhangba kerüljön az orvosi gyakorlat által támasztott követelményekkel. Kutatómunkám célja – ezen kihívásokkal összhangban – olyan új mérnöki módszerek és algoritmusok elméleti és gyakorlati kidolgozása, amelyek közvetlenül, a korábbiaknál hatékonyabban használhatók az orvosi gyakorlat diagnosztikai feladatainak megoldására. A doktori értekezésemben ismertetett munkámat a következő területeken végeztem: • Képminőség javítását eredményező nukleáris képalkotó berendezésekben alkalmazott képrekonstrukciós eljárások továbbfejlesztése. • Számítógépes tomográfiai (CT) felvételek alapján a foggyökér csatorna térbeli elhelyezkedésének meghatározására szolgáló, automatikusan végrehajtható, új képfeldolgozó eljárások kidolgozása. • Mágneses rezonancián alapuló (MR) képalkotó berendezéssel készített koponyafelvételek szegmentálására szolgáló képfeldolgozó módszerek létrehozása.

dc_592_12

Kutatási célkitűzések Kutatómunkám során a következő problémák megoldását tűztem ki célul: • Elnyelési korrekcióval és detektorválasz kompenzációval kiegészített SPECT képrekonstrukcióra szolgáló MLEM algoritmus párhuzamosított, GPU-ra szabott változatának kidolgozása. • Foggyökér csatorna térbeli elhelyezkedésének meghatározására szolgáló algoritmusok kidolgozása micro-CT, illetve cone beam CT (CBCT) felvételek alapján. • Térben lassan változó inhomogén zajjal terhelt mágneses rezonanciás (MR) képalkotóval készített felvételeket feldolgozó, a különböző agyi régiók elkülönítésére alkalmas, a korábbiaknál gyorsabb szegmentáló algoritmus kidolgozása.

3

dc_592_12 4

dc_592_12

1. fejezet Elnyelési korrekcióval és detektorválasz kompenzációval kiegészített párhuzamos SPECT képrekonstrukciós MLEM algoritmus GPU-ra adaptált változatának kidolgozása A SPECT (Single Photon Emission Computed Tomography – egy-fotonos számítógépes tomográfia) eljárás a mindennapi orvosi diagnosztika egyik alapvető képalkotó módszere [109]. A SPECT vizsgálatoknál számos, a képalkotás során fellépő fizikai jelenség befolyásolja negatívan az eredményül kapott kép minőségét, és ezáltal csökkenti azok klinikai diagnosztikai értékét [58]. Az orvosi diagnosztikai munka minőségének javítása érdekében szükséges ezen jelenségek negatív hatásának ellensúlyozása. A képminőség javítására szolgáló módszerek azonban olyan mértékben megnövelik a képrekonstrukciós algoritmusok számításigényét, és ezáltal azok végrehajtási idejét, hogy azok gyakorlati alkalmazása igen nehézkessé, számos esetben lehetetlenné válik [75]. Az elmúlt néhány évben a GPU (Graphical Processing Unit) kártyák rohamos teljesítménynövekedésével olyan párhuzamos számítást lehetővé tevő hardver eszközök jelentek meg az egészségügyi intézmények által megfizethető árkategóriában, amelyek megteremtik az elvi lehetőségét a képrekonstrukciós algoritmusok párhuzamos végrehajtásának. Ezzel a képrekonstrukció futási ideje olyan mértékben csökkenthető, amely lehetővé teszi azok

5

dc_592_12 6 gyakorlati alkalmazását. Ezen elvi lehetőség kiaknázásához azonban ki kell dolgozni a képrekonstrukciós algoritmusok olyan párhuzamos változatát, amely képes a hagyományos számítógépek működésétől eltérő, SIMD (Single Instruction Multiple Data) elven működő GPU-k számítási erőforrásainak hatékony kihasználására [57]. A kutatás keretében a következő két, a képalkotás során fellépő fizikai jelenség által okozott képminőség-romlás kompenzálásával foglalkoztam: • Gamma foton elnyelődés inhomogén elnyelő közegben; • SPECT detektor távolságfüggő felbontásából eredő torzítás. A fejezet ennek a kutatómunkának az eredményét mutatja be. A fejezet a következő részekből épül fel: a) a kutatási terület és a kutatás hátterének rövid ismertetése; b) probléma megfogalmazása; c) a kidolgozott algoritmus, illetve módszer leírása; d) eredmények, illetve azok elemzésének és validációjának összefoglalása.

1.1.

Nukleáris képalkotó módszerek alkalmazása az orvosi gyakorlatban

Az orvosi diagnosztikában alkalmazott képalkotó módszerek és berendezések máig tartó rohamos fejlődése töretlen volt az elmúlt húsz évben [75]. Ma már számos modalitás szolgálja a klinikai diagnosztikát, illetve az élettani kutatásokat. Az anatómiai struktúrákat leképező úgynevezett strukturális képalkotó modalitások (CT, MRI stb.) használata általános az orvosi gyakorlatban. Ezek mellett alkalmazzák az úgynevezett funkcionális képalkotó eljárásokat (SPECT, PET, fMRI stb.), amelyek egy adott szövet vagy szerv működéséről adnak diagnosztikai információt. A közelmúltban az orvosi gyakorlatban megjelentek a multimodális képalkotó berendezések, amelyek egyidőben képesek több, különböző eljárással készült képet készíteni az élő szervezetről [58]. Ezek alkalmazása nem csak a betegek kényelme miatt előnyös, hanem diagnosztikai hasznuk is van, hiszen egymással pontosan összeregisztrált, egyidőben készült képeket szolgáltatnak a betegről. A multimodális képalkotó berendezések általában funkcionális és strukturális képalkotókat ötvöznek (pl. PET/CT), hiszen diagnosztikai szempontból a két eltérő típusú modalitás jól kiegészíti egymást. A multimodális képalkotás megjelenése a képminőség javítása szempontjából is nagy jelentőségű, mert lehetőséget ad olyan képminőség javító eljárások kidolgozására, amelyek különböző modalitású felvételeket használnak.

dc_592_12 Párhuzamos SPECT rekonstrukció

7

1.1. ábra A SPECT vizsgálat során gamma sugárzó nyomjelző anyagot juttatunk a szervezetbe, majd a beteg teste körül adott pozíciókban vetületi képeken detektáljuk az aktivitás elhelyezkedését. A nyomjelző radioaktív anyag térbeli elhelyezkedését a vetületi képeket feldolgozó képrekonstrukciós algoritmusok határozzák meg.

A továbbiakban ismertetem a kutatáshoz kapcsolódó orvosi képalkotó eljárást.

1.1.1.

Egy-fotonos számítógépes tomográfia (SPECT)

A Single Photon Emission Computed Tomography (SPECT, egy-fotonos számítógépes tomográfia) képalkotó eljárás a modern orvosi gyakorlatban igen elterjedt nukleáris képalkotó módszer, amelyet különböző szakorvosi területeken (onkológia, neurológia, kardiológia stb.) [109] alkalmaznak. Az eljárás lényege, hogy gamma fotonokat kibocsátó radiopharmakont juttatnak a szervezetbe. A nyomjelző anyag a vizsgálat szempontjából fontos szervbe vagy szövetben fog feldúsulni. A SPECT vizsgálat célja, hogy meghatározzuk a radiopharmakon térbeli eloszlását a szervezetben, és ezáltal következtessünk a vizsgált szerv vagy szövet állapotára, funkcionális zavarára. Orvosi szempontból a vizsgálat legnagyobb előnye annak nem-invaziv volta, vagyis a beteg szervezetét annak feltárása nélkül lehet vizsgálni. A képalkotás során a páciens teste körül adott pozíciókban egy gamma fotonok érzékelésére kifejlesztett, úgynevezett gamma kamera fordul körbe (lásd 1.1. ábra). A gamma kamera projekciós képeket, a rekonstruálandó aktivitáseloszlás vetületi képeit rögzíti, mivel a kamera csak egy adott irányból – pl. csak a felületére merőlegesen – érkező fotonokat érzékeli. A képrekonstrukció célja, hogy a vetületi képekből meghatározzuk az eredeti aktivitáseloszlást a megfigyelt térrészben.

dc_592_12 8 A rekonstruált kép minősége, annak felbontása orvosdiagnosztikai szempontból igen fontos tényező. Jobb minőségű felvételek lehetőséget adhatnak pl. daganatos betegségek korai felismerésére és azok kezelésére, vagy kardiológiai vizsgálatoknál a szív állapotának pontosabb meghatározására. A képalkotás minősége sajnos számos tényező miatt nem tökéletes: • a beteg szervezetében a gamma fotonok elnyelődhetnek, szóródhatnak, ezeket a gamma kamera egyáltalán nem, vagy nem az eredeti becsapódás helyén fogja érzékelni; • a detektor véges felbontásából adódóan az aktivitáseloszlást csak egy adott pontossággal lehet rekonstruálni; • nem csak a kívánt irányból érkező gamma fotonokat detektálja a kamera, így nem létező képi elemek (más néven műtermékek) jelennek meg a képen; • a beteg vagy annak szervei elmozdulnak a vizsgálat során stb. A SPECT képalkotás számítógép tomográfiás képalkotás, vagyis a megfigyelt térrész vetületi képei alapján rekonstruáljuk a térrészben elhelyezkedő anyagok struktúráját. A vizsgált háromdimenziós térrész rekonstruálására szolgáló algoritmusok alapvetően két típusba tartoznak: analitikus és iteratív eljárások. A klasszikus vagy analitikus módszer a háromdimenziós térrész kétdimenziós vetületeit előállító Radon transzformáció inverzét megvalósító úgynevezett szűrt visszavetítést (Filtered Back Projection, FBP) alkalmazza [62]. Az ide tartozó algoritmusok számításigénye viszonylag alacsony, azonban csak részben képesek a képalkotás során fellépő képminőség rontó hatásokat kompenzálni [12, 61]. Az iteratív módszerek a megfigyelt térrészben elhelyezkedő anyagok struktúráját lépésről lépésre közelítve állítják elő. Egy kezdeti feltételezett anyageloszlásból kiindulva a képalkotás során fellépő fizikai jelenségeket modellezve létrehozzák a vetületi képeket, ezeket összevetik a valóságban detektált adatokkal, majd az eltérés alapján valamilyen algoritmussal módosítják a feltételezett anyageloszlást. Ezen lépéseket addig ismétlik, míg a leképezés eredményeként kapott és a megfigyelt vetületi képek közötti eltérés egy adott hibakorlátnál kisebb nem lesz. Az emissziós tomográfiában a leggyakrabban használt iteratív algoritmus az MLEM (Maximum Likelihood Expectation Maximization) algoritmus [51, 76]. Az iteratív algoritmusok alkalmazásának legnagyobb nehézsége azok nagy számításigénye [41]. Ennek ellenére – a számítástechnika rohamos fejlődése következtében – a használatuk egyre jobban terjed, mert lehetőséget adnak a különböző, képalkotás során fellépő képminőséget rontó jelenségek hatásának kompenzálására.


1.2.

9

Párhuzamos végrehajtást lehetővé tevő módszerek és platformok

A bevezetőben említett futási idő csökkentésének egyik módja a számítások párhuzamos elvégzése. Egy feladat elvégzéséhez szükséges számítások párhuzamosítása a számítástechnika, illetve informatika születése óta alkalmazott módszer. A párhuzamos végrehajtásnak a futási idő csökkentése mellett számos további előnye is lehet a megbízható, adott esetben hibatűrő rendszer létrehozásától az egyenletes erőforrás-kihasználásig. A módszert mégis akkor alkalmazzák leggyakrabban, mikor nagy számításigényű, így nagy futási idejű algoritmusok gyors végrehajtása a cél. Egy adott számításhoz tartozó műveletsor párhuzamosításakor végezhetünk úgynevezett adatpárhuzamosítást és feladatpárhuzamosítást. Adatpárhuzamosítás esetén a számítás során feldolgozandó adatokat osztjuk fel több részletre, és ezeken végezzük el egy időben a szükséges műveleteket. Feladatpárhuzamosítás esetén az eredeti számítási feladatot bontjuk kisebb részekre, majd ezeket hajtjuk végre párhuzamosan, különböző számítási erőforrásokon [50]. Egy feladat vagy algoritmus párhuzamosítása – legyen szó akár adat,- akár feladatpárhuzamosításról – közel sem triviális feladat. Egy új algoritmus megalkotásával jár, amely az eredeti, szekvenciális algoritmussal azonos eredményt szolgáltat, azonban adott lépéseit egymással egy időben, párhuzamosan lehet végrehajtani. Kiindulásként ehhez meg kell találni az eredeti algoritmus azon lépéseit, amelyek nem függenek egymástól, nincs köztük adatfüggés, így azok azonos időben végrehajthatóak. Ez általában csak a végrehajtandó algoritmus alapos elemzése után lehetséges. A párhuzamosítás második lépése a feladatot párhuzamosan végrehajtó algoritmust megvalósító program elkészítése. Párhuzamos programozást számos programnyelv, illetve programozási környezet támogat. Egy adott feladat implementációjához szükséges programozási környezet kiválasztása nagymértékben függ a megvalósítás platformjától, vagyis milyen szoftver-hardver környezetben szeretnénk, illetve tudjuk a párhuzamos algoritmusunkat megvalósító programunkat futtatni. A párhuzamos végrehajtást lehetővé tevő platformok gyors fejlődésen mentek át a közelmúltban, és ez a fejlődés jelenleg is tart. A kialakult rendszerek nagymértékben különböznek egymástól a fizikai méretükben, illetve a műveletvégző csomópontok csatolásának módjában. A párhuzamos rendszerek széles skálája jött létre, léteznek lazán csatolt, nagy méretű, általában térben is elosztott rendszerek pl. GRID, és koncentrált párhuzamos rendszerek pl. FPGA (Field-programmable Gate Array), vagy az eredetileg a számítógép

dc_592_12 10 grafikus funkcióinak gyorsítására szánt GPU (Graphical Processing Unit). A klasszikus szorosan, illetve lazán csatolt rendszertípusok határa mára elmosódott. Számos, e kettő típus között elhelyezkedő megoldás került alkalmazásra pl. különböző típusú cluster-ek, cloud megoldások stb. [50]. A következőkben röviden ismertetem az általam vizsgált probléma megoldásához használt párhuzamos alkalmazások futtatására alkalmas GPU számítási platformot.

1.2.1.

Graphical Processing Unit (GPU)

A Graphical Processing Unit (GPU) a számítógép azon komponense, amely a két (2D) és háromdimenziós (3D) képek előállításáért és megjelenítéséért felelős. Manapság egyre több számítógép rendelkezik önálló GPU egységgel, ami legtöbbször önálló kivehető kártyaként jelenik meg, amely a számítógép processzorától függetlenül képes ellátni feladatát, rendelkezik saját memóriával és művelet végrehajtó egységgel. A GPU-k hardver architektúrája a tipikus grafikus funkciók (pl. különböző térbeli alakzatok adott nézetben való megjelenítése, térbeli alakzatok forgatása, transzformációja, textúrák megjelenítése stb.) gyors ellátására lett optimalizálva (1.2. ábra). A GPU-kban nagyszámú, azonos felépítésű úgynevezett streaming processzor van integrálva, amelyek párhuzamosan képesek működni. A GPU-k memória felépítése hierarchikus, tartalmaz számos memória elérést gyorsító megoldást [66], mérete gyakran vetekszik a számítógép központi memóriájának kapacitásával. A 3D játékok, a számítógépes tervező rendszerek gyors terjedése motiválta a GPU gyártókat egyre gyorsabb és nagyobb teljesítményű egységek kifejlesztésére. A GPU-k architektúrája a SIMD (Single Instruction Multiple Data) feldolgozási logikát támogatja, amely esetén a párhuzamosan dolgozó processzorok azonos műveletsort hajtanak végre különböző adatelemeken. Ennek a feldolgozásnak a sémáját mutatja az 1.3. ábra. A grafikai számításokat általában nagyszámú lineáris algebrai művelet végrehajtása jellemzi. A SIMD architektúra ezen műveletek végrehajtására optimális platform. Mátrixműveleteket – amelyek során a mátrixelemeken külön-külön lehet az aritmetikai műveleteket elvégezni – igen gyorsan lehet GPU-k segítségével végrehajtani. A grafikai számítások igényeinek megfelelően a streaming processzorok aritmetikai egységei elsősorban egész számokon végzett műveleteket, illetve a lebegőpontos műveleteket támogatják. A bittérképes képek különböző felületekre illesztéséhez szükséges interpoláció műveletének elvégzését a GPU-k ugyancsak hardver támogatással gyorsítják.


11

1.2. ábra Az nVidia GPU belső felépítésének sémája. A GPU-ban elhelyezkedő műveletvégző elemek, a Streaming Processzorok (SP) úgynevezett Streaming Multiprocesszorokba (SM) szerveződnek egy dedikált cache memóriával, egy osztott elérésű (shared) memóriával és néhány úgynevezett Special Function Unit-tal (SFU, speciális utasításokat végrehajtó egységgel) együtt. Három SM blokk közösen használ egy második szintű úgynevezett texture (textúra) cache-t. Ezen egységek együtt alkotják a texture processing (textúra feldolgozó) cluster-t (TPC-t). A különböző típusú GPU kártyákon különböző számú TPC helyezkedik el, meghatározva az adott GPU teljesítményét.

1.3. ábra A GPU-kban alkalmazott SIMD feldolgozás logikai vázlata. A feldolgozó processzorokhoz azonos utasításkód kerül az osztottan használt utasítás buszon, míg a processzorok dedikált adatbuszon érkező különböző adatelemeken hajtják végre az adott utasítást.

dc_592_12 12 A GPU-k gyártói a közelmúltban felismerték a termékben rejlő lehetőséget, nevezetesen azt, hogy a GPU-kat nemcsak grafikai műveletek gyorsítására lehet használni, hanem általános célú számítási környezetként is. 2007-ben elsőként az nVidia cég jelent meg Compute Unified Device Architecture (CUDA) környezetével [111], amelyhez olyan segédprogramokat, és a használatukhoz szükséges dokumentációt is adott, amelyek lehetővé tették alkalmazások fejlesztését GPU-kra [66]. Időközben a teraflop kapacitásúra nőtt számítási teljesítmény igen népszerűvé tette a nagy számításigényű, de párhuzamosítható algoritmusok fejlesztői számára a GPU-kat. A GPU-k programozására lehetőség van C nyelvű alkalmazásokból a nyelv GPU specifikus kiterjesztése segítségével, de hasonló módon Fortran nyelven is lehetséges GPU programokat készíteni. Az OpenCL (Open Computing Language) 2009-ben való megjelenésével [63], illetve az időközben történt továbbfejlesztésével (legutóbbi változat OpenCL 1.2 2011. októberében) a GPU, mint általános célú számítási környezet (GPGPU-General Purpose GPU) egyértelműen a számítás-intenzív feladatok elsődleges megvalósítási platformja lett nem utolsó sorban a verhetetlen ár/teljesítmény arány miatt. Alkalmazások fejlesztése GPU környezetben Mint említettem a GPU-kra fejlesztett alkalmazások köre rohamosan bővül. Találunk a klasszikusnak számító képfeldolgozási alkalmazásoktól kezdve áramlásokat szimuláló alkalmazásokat, különböző célú Monte-Carlo szimulációkat, jelfeldolgozó alkalmazásokat, gépi látást megvalósító alkalmazásokat, de a matematikai szimbolikus számításokat támogató környezetek (Maple, Mathematica stb.) is megjelentek GPU platformon [39]. A jelenleg elérhető GPU-k (1.2. ábrán az nVidia GPU architektúráját láthatjuk) rendelkeznek – körülbelül 2008 óta – dupla pontosságú IEEE lebegőpontos aritmetikával, memóriájuk akár 1,5-3 Gbyte is lehet, és a PCI-Express bus-nak köszönhetően az adatcsere a központi memória és a GPU között igen gyors, ráadásul annak végrehajtása a számításokkal párhuzamosan is történhet. Különböző programozási felületek léteznek GPU-k programozására. Legelőször a Shader Model jelent meg, amely elsősorban grafikai alkalmazások végrehajtására volt specializálva. Ezt mindkét nagy GPU gyártó, az nVidia és az AMD (ATI) is támogatta. Mint említettem elsőként az nVidia jelent meg a CUDA (Compute Unified Device Architecture) általános célú számításokat is támogató modelljével, amelyre válaszul az AMD megjelentette a CMT (Close-to-Metal) megoldását, lehetővé téve a GPU-k magas szintű (pl. C vagy R) nyelven történő programozását. Ezeket a gyártóspecifikus környezeteket vélhetően a már említett OpenCL fel fogja váltani, amennyiben fejlesztése a mostani ütemben folytatódik, és ezzel a GPGPU-k programozására egy egységes kvázi szabvány környezet fog elterjedni.

dc_592_12 13

Párhuzamos SPECT rekonstrukció

1.3.

Probléma ismertetése

Gamma foton elnyelődés inhomogén elnyelő közegben A SPECT vizsgálatok során gamma sugárzó radioaktív izotóppal megjelölt anyagot juttatunk a szervezetbe. A SPECT berendezésben levő gamma kamerákkal a bomlás során keletkező gamma fotonokat érzékeljük. A bomlás helye és a detektorba történő becsapódás között a gamma fotonok áthaladnak ezen két pont között elhelyezkedő térrészen, amely térrészben levő anyag részecskéivel kölcsönhatásba kerülhetnek, ütközhetnek velük. Ezen kölcsönhatások eredményeként egy adott detektorterület által detektált események száma csökken. A csökkenés mértéke függ a gamma foton útja során érintett anyag jellemzőitől. Egyes anyagok gamma foton elnyelő hatását az anyag elnyelési együtthatójával jellemezzük. Az elnyelődés mértéke a gamma foton adott anyagban megtett út hosszától, és az anyag elnyelési együtthatójától függ. Az elnyelődést is figyelembe véve egy detektor elemben érzékelhető fotonok számát a Beer-Lambert törvény alapján számíthatjuk: I = I0 · e−

R l

µ(x)dx

(1.1)

Az összefüggésben I a mért események száma a detektor felbontása által meghatározott egységnyi felületén (bin-jén), I0 a tér egy adott voxeléből a detektor bin irányába induló gamma fotonok száma, µ(x) az adott voxelből a bin-be vezető egyenesen levő voxelekben elhelyezkedő anyagot jellemző, adott irányban érvényes gyengítési együtthatókat tartalmazza. Inhomogén elnyelő közeg esetén az elnyelődés erősen irányfüggő lesz. Az inhomogén elnyelő közeg ezért a SPECT rekonstrukció során előállított képen igen súlyos torzítást eredményezhet (lásd 1.4. ábra) [99]. A SPECT vizsgálatoknál általánosan használt radioaktív izotópok esetén az általuk kibocsátott gamma fotonok energiája általában a 70 keV és 400 keV közötti energiatartományba esik. Ebben az energiatartományban a gamma fotonok és a gamma fotonok útja során érintett anyagok között fellépő kölcsönhatások közül a foton abszorpció és a Compton szórás lép fel, amelyek közül a Compton szórás dominál. Tehát a Beer-Lambert törvény segítségével kiszámítható gamma foton szám csökkenés oka lényegében a gamma fotonok kiszóródása az adott irányból, illetve térrészből. A Beer-Lambert közelítés azt feltételezi, hogy a Compton szórás után továbbhaladó gamma foton detektorba érkezésekor keletkező események közül egyet sem érzékelünk. Ezzel ellentétben a valóságban a detektor energiafelbontása, vagyis az általa érzékelt fotonok energiája nem egyetlen energiaszint, hanem egy energiatartomány. Így ha a szórás után egy foton energiája még a detektor érzékenységi tartományában marad, azt a detektor

dc_592_12 14

(a) Elnyelési térkép és az eredeti aktivitáseloszlás

(b) Rekonstruált szelet

1.4. ábra A képen megfigyelhetjük a különböző elnyelési állandójú anyaggal kitöltött tér okozta torzító hatást a SPECT rekonstrukció során. Egy henger falában viszonylag magas aktivitású sugárzó anyagot helyeztünk el (piros, illetve sárga területek). A hengert egy kis elnyelési állandójú (ábrán középkék színű) hengerbe helyeztük, amelyben aszimmetrikusan egy nagy elnyelési állandójú (világoskék színnel jelölt) henger is volt. Az eredeti aktivitáseloszlás metszeti képét a bal oldali ábra, míg a rekonstruált kép metszeti képét a jobb oldali ábra mutatja. A jobb oldali ábrán a henger aktivitást tartalmazó falának a szórás okozta torzulása, kiszélesedése látható, illetve az eltűnése a nagy elnyelési állandójú térrészben.

érzékelni fogja. Az így kapott eseményeket nevezzük szórt eseménynek. A szórt események okozta torzító hatás kiküszöbölése érdekében a Beer-Lambert törvényre alapuló elnyelési korrekció alkalmazása esetén meg kell valósítani a szórt események hatását kompenzáló szóráskorrekciót. A szórt események torzító hatásának kiküszöbölése elsősorban humán SPECT vizsgálatok esetén kap nagy jelentőséget, mert ott gyakran akár 50% körüli is lehet a szórt események aránya. A szórt események okozta hatás csökkentésének egyik lehetséges módszere a SPECT berendezés módosításával a több energiaablakos adatgyűjtés [46], [48] használata. Ekkor a különböző energiatartományokban mért események száma alapján bizonyos közelítéssel megbecsülhető a szórt események száma. Az adatgyűjtés során ugyan lesznek szórt eseményekből eredő becsapódások az elsődleges bomlás okozta detektált események energiatartományában is, de a szórt eseményekre vonatkozó becslés alapján ezek számát csökkentjük. A módszer alkalmazásához szükség van a speciális adatgyűjtésre alkalmas berendezésre, amely viszonylag drága, ráadásul a szórt események becslése is csak adott körülmények között pontos, így az ilyen berendezések nem terjedtek el széleskörűen a gyakorlatban. Amennyiben tehát a foton elnyelődés és szórás torzító hatását ki akarjuk küszöbölni, azt a rekonstrukciós algoritmus végrehajtása során kell figyelembe venni.


15

1.5. ábra A gamma kamera felépítése. A megfigyelt térrész felől egy adott irányból (általában merőlegesen) érkező gamma fotonokat a kollimátor hivatott kiszűrni. A szcintillációs kristályba becsapódó gamma fotonok látható tartományban levő optikai fotonokat váltanak ki, míg a PMT-k az így keletkezett fotonok számát sokszorozzák meg, hogy azokat az elektronika érzékelni tudja.

SPECT detektor távolságfüggő felbontásából eredő torzítás SPECT készülékben a gamma kamera (1.5. ábra) feladata a radioaktív nyomjelző anyag bomlása során kibocsátott gamma fotonok érzékelése. A gamma kamera három fő részből áll: • kollimátor, • detektorban elhelyezett szcintillációs kristály, • fotoelektron-sokszorozó csövek (PMT-k, photomultiplier tube-ok). A kollimátor feladata, hogy csak egy adott irányból (általában a detektor felületére merőlegesen) beeső gamma fotonok jussanak el a kollimátor mögött elhelyezkedő szcintillációs kristályhoz. A kollimátor leggyakrabban egy viszonylag vastag ólomlemez, amelybe a kiszűrni kívánt iránynak megfelelően lyukak vannak fúrva. A leggyakrabban használt úgynevezett „parallel hole” (párhuzamos vetítésű) kollimátor lyukai – amelyek a felületre merőleges fotonokat engedik át – az ólomtömb vízszintes felületére merőlegesen, egymással párhuzamosan helyezkednek el.

dc_592_12 16 Ha a gamma foton túljut a kollimátoron, akkor a mögötte elhelyezkedő szcintillációs kristály tömbbe csapódik. A kristályban az elnyelődés során a gamma fotonból optikai fotonok keletkeznek. Az így keletkező optikai fotonok száma viszonylag kicsi, ezért van szükség a PMT-kre. A PMT-k a fotonok számát akár több nagyságrenddel is megsokszorozzák, amely fotonok okozta fényfelvillanást már képes a PMT-k mögött elhelyezett fénydetektáló elektronika érzékelni. Ha egy gamma kamerában elhelyezünk egy pontszerű sugárforrást, akkor annak detektált képe, vagyis a pontforrásból a detektorban levő elektronika által érzékelt becsapódások helye egy kissé elmosódott korong lesz, amely korongon a becsapódások kétdimenziós Gauss eloszlást fognak követni (1.6. ábra). A gamma kamerákban tapasztalható úgynevezett Gauss-os elkenődésnek három oka is van: 1. A detektorban levő szcintillációs kristályba becsapódó gamma foton a kristályon belül véletlenszerűen valamilyen kis távolságba eljut, tehát ha egy foton azonos pontban, azonos szögben is érkezik a kristály felületére, akkor sem mindig ugyan abban a pontban, hanem a kristály egy régiójában történik az elnyelődés, ami önmagában a felvillanások helyének – Gauss-os eloszlással jellemezhető – bizonytalanságot ad. 2. A parallel kollimátor lyukai bár nagyon kicsik, de véges méretűek. Ez lehetővé teszi, hogy ne csak egy egyenes mentén, hanem egy igen kis nyílásszögű kúpszerű térrészből legyen elérhető a kristály felszíne a gamma fotonok számára. 3. A kollimátor anyaga ólom, ami ugyan nagyon hatékonyan, azonban nem tökéletesen árnyékol, vagyis nem nyeli el az összes gamma fotont. A kollimátorban tehát létrejöhet a kollimátor lyukai között a gamma fotonok penetrációja. A gamma kamerán a pontforrás képe tehát egy, az eredeti pont méreténél nagyobb folt lesz, amely folt nagysága függ a pontforrás detektorfelülettől mért távolságától. Ezt a jelenséget láthatjuk az 1.7. ábrán. A jelenség okozza a kamera távolságfüggő felbontását.

Probléma megfogalmazása A fent leírt SPECT képalkotás során fellépő fizikai jelenségek lényegesen rontják a rekonstruált képek minőségét. Ez indokolja, hogy – főként inhomogén elnyelési közeg esetén – mind az elnyelés okozta torzítást [99], mind a detektor távolságfüggő felbontásának hatását [98] korrigálni kell a képalkotás során a jó minőségű rekonstruált kép előállítása érdekében [121].


17

1.6. ábra A gamma kamera távolságfüggő felbontását a foton becsapódások Gauss-os eloszlásával lehet jellemezni a detektor felületén, amely eloszlás (Psf1, Psf2) paraméterei a sugárforrás detektorfelülettől mért távolságának függvényében (P1, P2) változnak.

A fenti két fizikai jelenség torzító hatását korrigáló módszerek beépítése a SPECT képrekonstrukcióra általánosan alkalmazott Maximum Likelihood Expectation Maximization (MLEM) algoritmusba [51, 76] drasztikusan megnöveli az algoritmus számítási bonyolultságát. Ez a fokozott számításigény a jelenlegi klinikai alkalmazás gátja, mert a rekonstrukció futási ideje már nem tolerálható valós alkalmazási környezetben. A kutatómunkám célja a SPECT vizsgálatoknál alkalmazott MLEM rekonstrukciós algoritmus elnyelési korrekciót és a távolságfüggő detektorfelbontást is korrigáló olyan változatának kidolgozása, amely képes a klinikai gyakorlat támasztotta futási idő korlátokat teljesíteni.

1.4.

Megoldási módszerek

1.4.1.

Iteratív képrekonstrukció

A háromdimenziós emisszió tomográfiás képalkotás során a képrekonstrukciós algoritmus alkalmazásának célja, hogy meghatározzuk az emberi szervezetbe juttatott nyomjelző anyag térbeli elhelyezkedését a testben, a testről készült projekciós képek alapján. Mint korábban azt már említettem, az emissziós tomográfiában erre a célra leggyakrabban használt robusztus, iteratív algoritmus az MLEM [51, 76] algoritmus. A gyakorlatban az MLEM algoritmus módosított változatát használjuk, amely a projekciós képeket úgynevezett subset-ekbe (csoportokba) szervezi, és ezeken külön hajtja végre az iterációkat [41]. Az így végrehajtott MLEM algoritmust Ordered Subset Expectation Maximization (OSEM) algoritmusnak hívják, előnye az eredeti változatnál jobb konvergencia. A képre-

dc_592_12 18

(a) Referencia

(b) OSEM algoritmussal rekonstruált kép

(c) Profilgörbék

1.7. ábra Matematikai henger fantom szimulált SPECT felvételének keresztmetszeti képe. A referencia képen (1.7(a)) látható az eredeti aktivitáseloszlás, míg a rekonstruált képen (1.7(b)) a kamera távolságfüggő felbontása okozta elkenődés. Az ábra alsó felén (1.7(c)) a két felső ábrán látható gyűrűk aktivitáseloszlása van ábrázolva egy, a középponton áthaladó egyenes mentén.

konstrukció során végrehajtott számítások menete a két algoritmus esetén azonos, csak a projekciós képek feldolgozásának sorrendjében térnek el egymástól. Emiatt a továbbiakban az általam javasolt módosításokat az MLEM algoritmusra, illetve annak egyes lépéseire hivatkozva mutatom be, attól függetlenül, hogy a tényleges implementáció során az MLEM OSEM változatát használtam mind az új, mind a referenciaként használt eredeti algoritmus esetén. Az MLEM a j-dik voxelben levő fj aktivitás becslésére a (k + 1)-dik iterációban a következő összefüggést adja ([76] alapján):

dc_592_12 19


fjk+1

=

fjk

· PI

1

i=1

aij

·

I X i=1

bi · aij k n=1 fn · ani

PJ

(1.2)

ahol J a rekonstruálandó tér voxeleinek halmazát jelöli, I a detektor pixeleinek halmaza, aij annak a valószínűsége, hogy egy adott j voxelben keletkezett gamma fotont az i detektor pozícióban érzékelünk. Parallel hole (párhuzamos vetítésű) kollimátor alkalmazása esetén ez a valószínűség arányos az i-vel jelölt detektor pozícióból kiinduló merőleges egyenesnek az adott, j-vel jelölt voxelbe eső (vagyis a voxel által kimetszett) szakaszának hosszával. Az összefüggés lényegében minden iterációban kiszámolja, hogy mi lenne a detektált projekciós kép a korábbi iterációban kapott aktivitáseloszlás alapján – ezt a lépést nevezzük előrevetítésnek –, majd a számolt értékek és a detektor egyes pixeleiben mért értékek különbsége alapján korrigálja az egyes voxelek aktivitás értékeit. Ez utóbbi lépést nevezzük visszavetítésnek. Az (1.2) összefüggés előrevetítő operátora az összefüggés első tényezője (a nevezőben szereplő szumma), míg a visszavetítő operátor az összefüggés második része. A visszavetítő operátorban bi az i-dik detektorpixelben mért fotonszám. Az előrevetítő és a visszavetítő operátor használatát egy egyszerű rekonstrukciós példán az 1.8. ábra szemlélteti. A példa az ábra bal oldalán (1.8(a)) egy adott detektor állásban, egy 7x7 voxelt tartalmazó térrészben, adott aktivitáseloszlás esetén a detektor pixeleiben számolt (becsült) beütésszám számítását demonstrálja. A példán a detektált beütésszám számításakor az egyes voxelekben elhelyezkedő aktivitás értékek hatását az egyes voxelek középpontjából induló, detektorra merőleges egyenes mentén „vetítjük” a detektorra. Jobb oldalon (1.8(b)) az aktivitáseloszlás módosításának logikáját jelöli a becsült és a valóságban detektált (ábrán külön nem jelölt) beütésszám arányában. A módosítás mértéke függ az egyes voxelek, illetve az azokban elhelyezkedő aktivitás detektorpixeleken érvényesülő hatásától (ezt az ábra külön nem jelöli). Az MLEM algoritmus egyik előnyös tulajdonsága, hogy kiegészíthető az előre-, illetve visszavetítő operátor a korábban említett fizikai jelenségek (elnyelés, távolságfüggő felbontás) hatásának modelljével, így azok torzító hatása jórészt kiküszöbölhető. Analitikus rekonstrukció során igen nehéz az elnyelő közeg hatását figyelembe venni. A Radon transzformáció inverzét homogén elnyelőközeg esetén használják [61], elsősorban kísérleti környezetben. Az orvosi gyakorlatban azonban nem beszélhetünk homogén elnyelő közegről. Inhomogén közeg esetén is születtek közelítő megoldások, azonban ezek sem bizonyultak alkalmazhatónak a gyakorlatban [65]. Az MLEM alkalmazásakor azonban az egyes térrészekben elhelyezkedő anyagokra jellemző elnyelési faktort figyelembe tudjuk venni, ha az (1.2) összefüggésben szereplő aij valószínűség értékeket az (1.1) összefüggés alapján határozzuk meg. Ennek megfelelően az

dc_592_12 20

(a) előrevetítő operátor

(b) visszavetítő operátor

1.8. ábra Az előre és a visszavetítő operátor használatának szemléltetése egy egyszerű rekonstrukciós példán. Bal oldalon (1.8(a)) egy adott detektor állásban egy 7x7 voxelt tartalmazó térrészben adott, szürke árnyalatokkal jelölt aktivitáseloszlás esetén a detektorpixelekben számolt (becsült), az ábrán hisztogrammal jelölt beütésszám számítását demonstrálja. Jobb oldalon (1.8(b)) az aktivitáseloszlás módosításának logikáját jelöli a becsült és a valóságban detektált beütésszám függvényében.

egy detektorpixel által érzékelt események számát (diszkrét esetben) a következő összefüggés adja:

bi =

J X

−

xj · aij · e

P n∈Nij

cin ·µn

(1.3)

.

j=1

Ezt az MLEM algoritmus eredeti (1.2) formulájába behelyettesítve a következő összefüggést kapjuk:

fjk+1

=

fjk

·P I

1 −

i=1 aij e

P n∈Nij

cin ·µn

·

−

I X i=1

bi e PJ

n=1

P n∈Nij

−

fnk · ani e

cin ·µn

P h∈Nin

cih ·µh

· aij .

(1.4)

Az algoritmus számításakor élhetünk néhány gyorsítással, mert az összefüggés néhány részszámolásának eredményét újra felhasználhatjuk. A rendszer geometriai viszonyaiból – vagyis abból, hogy melyik detektorpixel mely voxelekből „látható” – adódó összefüggéseket, valamint az egyes voxelekben a foton elnyelés valószínűségét előre kiszámíthatjuk, és egy – [voxelek száma] × [detektorpixelek száma] értéket tartalmazó – mátrixban eltárolhatjuk, amely mátrixot a szakirodalom leggyakrabban rendszermátrixnak nevez. Ekkor az (1.4) összefüggés kiszámítása a rendszermátrix egyes sor, illetve oszlopvektorainak az

dc_592_12 21


aktivitás értékekből, illetve az egyes detektorpixelek beütésszámából képzett sor, illetve oszlopvektorok szorzatára egyszerűsödik. A rendszermátrixnak annyi sora van, ahány voxelben szeretnénk rekonstruálni, és annyi oszlopa, ahány detektált pixel található az összes vetületi képben. Ha a rekonstrukciót egy minden oldalán l voxelt tartalmazó voxeltömbön végezzük, és n darab k × k pixelt tartalmazó projekciós képünk van, akkor az S rendszermátrixunk felépítése a következő lesz



S=

        



s1,1

s1,2

s2,1

s2,2

... s1,k×n×n   ... s2,k×n×n  

...

...

...

...

sl×l×l,1 sl×l×l,2 ... sl×l×l,k×n×n

    

(1.5)

ahol si,j az i-dik voxel j-dik pixelre való hatását mutatja. A napjainkban gyakorlatban alkalmazott átlagos képalkotó rendszerben jellemző, hogy 64 db 128x128-as projekciós képsorozat alapján szeretnénk egy minden oldalán 128 voxelt tartalmazó térrészt rekonstruálni. Egy ilyen rendszerben a rendszermátrix egyik dimenzióban 64x128x128, másik dimenzióban 128x128x128 méretű lesz, ami összességében 1012 nagyságrendű érték tárolását jelenti. Mint említettem, minden mátrixelem azt adja meg, hogy az adott vetületi kép adott pixelébe milyen valószínűséggel érkezhet sugárzás egy adott voxelből. Mivel minden voxel csak a detektorok kis részének adott pixeleire van/lehet hatással a képalkotó rendszer geometriai felépítése miatt, így a rendszermátrix egy ritka mátrix lesz. Ez teszi lehetővé azt, hogy ezt a viszonylag tekintélyes méretű mátrixot el tudjuk tárolni, és a rekonstrukció során újból és újból felhasználni a ma hozzáférhető számítógépeken. A korábban használt jelölésekkel összhangban legyen bk+1 a (k + 1)-dik iterációban i az i-vel jelölt detektorvoxelben számított (becsült) beütésszám, S a rendszermátrix, ej a j-dik pozícióban egy értéket tartalmazó egységvektor, ekkor a bk+1 értékét a következő i összefüggéssel számíthatjuk ki:

bk+1 = (S × ej ) × f k i

(1.6)

ahol f k a k-dik iterációban számolt, l × l × l darab fik értéket (i = 1, 2, . . . , (l × l × l)), vagyis a k-dik iterációban becsült összes aktivitás értéket tartalmazó sorvektor.

dc_592_12 22 Fizikai torzító hatásokat is figyelembe vevő korábbi rekonstrukciós módszerek A korábbi fejezetben leírt, az elnyelési együttható iteratív rekonstrukció előre, illetve visszavetítő operátorának számításakor figyelembe vett hatásának rekonstrukcióba történő beépítésétől eltérő, alternatív lehetőség a Bellini által publikált, analitikus rekonstrukció esetén alkalmazható módszer [15]. Ez azonban csak homogén elnyelő közeg esetére ad kielégítő megoldást, így esetünkben nem alkalmazható. Chang és társai által javasolt [26] módszer előnye, hogy inhomogén elnyelő közeg esetén is alkalmazható, azonban ez a módszer is számos közelítéssel él, így gyakran pontatlan eredményre vezet. A számításigényes iteratív SPECT rekonstrukciós módszerek erőforrásigényének az algoritmusok párhuzamosításával történő kompenzálására korábban csak néhány kísérlet történt. Ezek közül megemlítendő McCarthy és Miller [59], akik speciális vektor processzoros architektúrára készítettek elnyelési korrekciót és távolságfüggő detektorfelbontást is figyelembe vevő párhuzamos rekonstrukciós algoritmust. A módszer – elsősorban speciális hardver igénye miatt – kereskedelmi rendszerekbe nem került alkalmazásra. Ez a publikáció említi először a „slice-by-slice blurring” (szeletenkénti elkenés) eljárást, amely a távolságfüggő detektorfelbontás hatását a detektorral párhuzamos síkok mentén levő szeleteken lépésről lépésre számolja, így véve figyelembe a hatás távolságfüggését. Fizikai torzító hatásokat is figyelembe vevő iteratív SPECT rekonstrukció GPU-n történő megvalósítására egyetlen példa a Vetter és társai által [104] publikált munka, amely az úgynevezett slab-by-slab [11] módszert használta, amely nem szeletenként, mint a „slice-byslice blurring” módszer, hanem azok csoportjain iterálva veszi figyelembe a távolságfüggő viselkedést. A javasolt módszer validációjáról nem áll rendelkezésre irodalmi adat. A nukleáris medicinában már korábban is alkalmaztak párhuzamosított módszereket. Néhány az általunk is alkalmazott MLEM rekonstrukció OSEM változatát is használta, azonban ezek közül egyik sem SPECT felvételek rekonstrukciójára alkalmazta a párhuzamos algoritmust. A területen született fontosabb publikációk a következők: • Fang Xu és Mueller [115] inhomogén elnyelő közeg hatását is figyelembe vevő párhuzamos rekonstrukciós algoritmust dolgoztak ki. Nem vették figyelembe a detektor távolságfüggő felbontását. Elsősorban CT rekonstrukcióval foglalkoztak [116], néhány algoritmusuk alkalmazható PET (Pozitron Emissziós Tomográfiás) képek rekonstrukciójához [114]. • Bergner és társai [20] egy Monte Carlo módszert dolgoztak ki PET rekonstrukcióra. Nem foglalkoztak a detektor távolságfüggő felbontásával.


23

• Xu és társai [117] CT rekonstrukcióhoz dolgoztak ki GPU-n futó párhuzamos algoritmust. Ugyan felismerték a „gyűjtő típusú” algoritmus kidolgozásának szükségességét, de meglepő módon az általuk kidolgozott párhuzamos algoritmus OSEM változata lassabban futott, mint az eredeti MLEM algoritmus.

Elnyelési térkép meghatározása Elnyelési korrekciót tartalmazó rekonstrukciós algoritmus alkalmazásához szükség van egy pontos elnyelési térképre. Az elnyelési térkép előállításának több alternatív módja használatos. Gyakori módszer külső sugárforrással (úgynevezett transzmissziós) felvétel készítése, amelyet rekonstruálva megkaphatjuk az elnyelési térképet [119]. A közelmúltban megjelent multimodalitású SPECT készülékek terjedésével – amelyekben a SPECT készülék mellett egy CT berendezés is rendelkezésre áll – egyszerűbb módja kínálkozik az elnyelési térkép előállításának. Ezekben az eszközökben a CT szolgáltatja a morfológiai információt a vizsgált testről, a CT képből pedig kinyerhető az elnyelési térkép [119]. Az elnyelési térkép kinyerésének több módszere is ismert: • skálázás, • szegmentálás, • hibrid módszer. Skálázás esetén a CT felvételekben levő információ alapján közvetlenül határozzuk meg térrészek elnyelési együtthatóját. A CT felvételeken az egyes voxelekben elhelyezkedő anyag jellemzésére az úgynevezett Hounsfield skálát használjuk, az anyag sűrűségét Hounsfield Unit-ban (HU) mérjük. A HU-ban mért sűrűség egyszerű konverzióval átszámítható az előre-, illetve visszavetítő operátorokban szereplő µ elnyelési együtthatóra. Ez a becslés a víz sűrűségével közel azonos sűrűségű anyagokra viszonylag pontos elnyelési együttható értékeket szolgáltat, azonban porózus, kevés vizet tartalmazó szöveteknél már nem használható. Szegmentálás esetén a CT felvételeken az egyes szerveket, szöveti részeket azonosítjuk, meghatározzuk azok határait, majd az egyes szervekre, szövetekre jellemző elnyelési értéket társítunk az adott szerv vagy szövet által elfoglalt térrészhez. A hibrid módszer a fenti két módszert kombinálja a pontosabb elnyelési együttható meghatározása érdekében. Az elnyelési térkép meghatározására az általunk végzett kísérletek során a hibrid módszert alkalmaztuk.

dc_592_12 24

1.4.2.

Detektor távolságfüggő felbontásának kompenzációja

Detektor távolságfüggő felbontásának kompenzációját megvalósító algoritmusokat három típusba sorolhatjuk: • A projekciós képen jelentkező torzítást a rekonstrukció előtt a projekciós képen alkalmazott transzformációval kompenzálhatjuk. A transzformáció általában a frekvencia térben, valamilyen szűrő (Winer filter, Metz filter) [47], [96] alkalmazásával történik. Ennek a módszernek a hátránya, hogy használatával nem lehet figyelembe venni a fizikai hatás távolságfüggését, így pontatlan megoldáshoz vezet, amely pontatlanság a megfigyelt térrész méretével arányosan növekszik. • Hasonlóan a projekciós képeket eredményező detektált adathalmazon végzünk transzformációt a Frequency Distance Relation (FDR) [33] összefüggésen alapuló úgynevezett FDR szűrős módszerrel, amely módszer előnye a korábban említett szűréssel szemben, hogy képes figyelembe venni a távolságfüggő torzító hatást [113], [38], [53]. Ennél a módszernél a transzformációt a detektált adatok úgynevezett sinogramján (a detektált adatok 2D rekonstrukciókhoz átrendezett változatán) végezzük. A módszer hiányossága pont ebből adódik. A távolságfüggő torzító hatás kompenzációja csak egy-egy rekonstruált szeleten belül fog érvényesülni, hiszen a rekonstrukció minden lépésében csak egy-egy szeletét állítjuk elő a rekonstruált térrésznek. • Detektor távolságfüggő felbontása kompenzációját megvalósíthatjuk a rekonstrukciós algoritmus végrehajtása során, az eredeti detektált adatok módosítása nélkül. Erre iteratív rekonstrukció esetén van lehetőség. Ebben az esetben lehetséges a torzító hatás távolságfüggésének pontos figyelembe vétele, azonban a rekonstrukció számításigénye jelentősen megemelkedik [100, 11, 59, 120, 46]. A detektor távolságfüggő felbontása kompenzációjának megvalósításakor a harmadik, iteratív rekonstrukciós algoritmusok esetén alkalmazható változata mellett döntöttünk annak – a többi módszerhez képest viszonyított – pontossága miatt. A módszer számításigényét annak párhuzamosításával fogjuk kompenzálni. A távolságfüggő detektorfelbontás számítását az előrevetítő operátorba építettük be. A detektorválasz függvény módosulását úgy lehet egyszerűen figyelembe venni, hogy a detektor felületén a Gauss eloszlás szerint elosztjuk, „elkenjük” a detektorhoz érkező fotonok hatását. Amikor az aktuális aktivitáseloszlásból becsüljük a rendszer által érzékelt beütésszámokat (projekciós képeket), akkor az adott előrevetített aktivitás által keltett pontválasz értékét lényegében szétosztjuk egy térben, kétdimenziós Gauss eloszlásnak megfelelően a projekciós kép pixelei között. A távolságfüggő detektorfelbontás kompenzálása nélkül a

dc_592_12 25


1.9. ábra A detektor távolságfüggő felbontásának kompenzációjához szükséges kalibráció menetének szemléltetése. A kép felső részén a különböző távolságban elhelyezett vonalaktivitás detektált képe látható, míg az ábra alsó részén az egyes mérések alapján készített beütésszám hisztogramok. Ezek alapján határozhatók meg a különböző távolságokhoz tartozó eloszlások paraméterei.

LOR-on számolt teljes hatás azt a detektorpixelt érné, amelyiket az adott LOR kijelöl. A kompenzáció alkalmazása esetén is ez a pixel fogja megkapni a LOR mentén számított beütésszám legnagyobb részét, azonban a Gauss eloszlásnak megfelelően a környező pixelekben is módosítjuk a beütésszámot. Az elkenés mértéke – mint azt az 1.6. ábra is jól szemlélteti – függ a távolságtól. Egy adott detektor pontválasz függvénye, vagyis az, hogy milyen távolság esetén pontosan milyen Gauss eloszlásfüggvényt kell alkalmazni, kimérhető egy kalibrációs mérés sorozattal, amely során egy kapillárist helyezünk a detektortól adott távolságokra (1.9. ábra). Az egyes planáris képeken a kapilláris képére Gauss függvényt illesztve (1.10. ábra) meghatározható a kollimátor okozta – Gauss eloszlással jellemezhető – foton becsapódás gyakoriságának várható értéke, illetve szórása. A Gauss eloszlás σ szórása lineárisan függ a detektortól mért távolságtól. A kalibrációs mérésből meghatározhatjuk az egyenes meredekségét jellemző psf b, illetve az x tengellyel való metszéspontját jelző psf a paramétereket, amelyek felhasználásával a detektortól bármilyen távolságra elhelyezkedő pontban elhelyezett pontforrás okozta foton becsapódásokat leíró Gauss eloszlás szórása kiszámolható:

σ = psf a + psf b · d ahol d a detektor felületétől (a kollimátortól) számított távolság.

(1.7)

dc_592_12 26

1.10. ábra A kalibráció során felvett mérési pontokhoz tartozó eloszlások szórása (σ). A pirossal jelölt egyenes a szórás értékekre illesztett egyenes. Jól látható a lineáris függés a detektorfelülettől mért távolság függvényében. A pontokra illesztett egyenes x tengellyel való metszéspontja adja a psf a paramétert, míg a meredeksége a psf b paramétert.

Az ilyen módon, a távolság függvényében meghatározott eloszlásfüggvényt a rekonstrukció során úgy vesszük figyelembe, hogy a ténylegesen kalkulált becsapódás helye alapján meghatározzuk a Gauss eloszlás várható értékét, majd az (1.7) összefüggés alapján számított szórás, illetve az eloszlásfüggvény alapján megadjuk, milyen mértékben kell az adott hatást elosztani a becsapódás környékén elhelyezkedő detektorpixelek között. A következő fejezetben bemutatom a javasolt képrekonstrukciós algoritmust, amely megoldást ad a célkitűzésben ismertetett problémára.

1.4.3.

GPU-ra szabott rekonstrukciós algoritmus

A különböző, korábban ismertetett párhuzamos feldolgozásra alkalmas platformok közül a SPECT rekonstrukciós algoritmus implementálására a GPU platformot, az nVidia cég különböző GPU-it választottam. A GPU platform választásának egyértelmű indoka a teljesítmény/ár arány, amely tekintetben ezek az eszközök lényegesen előnyösebbek más lehetséges megoldásoknál. Esetünkben, mivel az algoritmust futtató eszközök várhatóan számos SPECT berendezésbe beépítésre kerülnek, szempont volt a kis fizikai kiterjedés, amely ugyancsak a GPU-k alkalmazását indokolta. Az nVidia cég GPU-i mellett elsősorban a CUDA nyelven történő fejlesztés lehetősége szólt, amely elterjedt és a gyártó által jól támogatott. Megjegyezzük, hogy a kidolgozott párhuzamosított algoritmust bármilyen


27

más GPU implementációs környezetben meg lehet valósítani úgy, hogy kihasználhatók a párhuzamos algoritmus nyújtotta előnyök. A GPU-ra optimalizált algoritmus kidolgozásához először vizsgáljuk meg a hagyományos, CPU-n futó – azon a platformon párhuzamosított – algoritmust, és nézzük meg, miért nem lehet azt optimálisan implementálni, illetve futtatni GPU-n. A hagyományos számítási környezetre készült algoritmus az előrevetítő operátor által meghatározott becsapódásszám értéket az aktivitást tartalmazó voxelekből kiindulva számítja ki. Ez azt jelenti, hogy a párhuzamosan futó szálak egymás után megvizsgálják az egyes voxeleket és kiszámítják, hogy az egyes voxelek, illetve a bennük elhelyezkedő aktivitás milyen mértékben járul hozzá az egyes detektorpixelekben érzékelhető becsapódási események számához. Így minden számítási szál módosítja mindegyik detektorpixelhez tartozó projekciós képet, illetve számolt becsapódási számot. Ez maga után vonja, hogy minden szál minden detektorpixelhez tartozó mátrixelemet el kell, hogy érjen, illetve az egyes szálaknak a detektorpixelhez tartozó becsapódásszám értéket tartalmazó mátrixelemek elérésekor szinkronizálniuk kell az adatelérést az adatkonzisztencia megőrzése érdekében. Ennek a működésnek a GPU-n történő változatlan megvalósítása két szempontból is hátrányos. Egyrészt a GPU-n, egyes számítási magokon futó szálak csak viszonylag kis memóriaterületet tudnak gyorsan elérni. Ez megnehezíti a teljes rekonstruálandó térrészben elhelyezkedő aktivitást leíró adatmátrix elérését, hiszen az lényegesen nagyobb, mint a szálak által gyorsan elérhető memóriaterület. Másrészt, ha a szálaknak minden alkalommal szinkronizálódniuk kell, amikor módosítanak egy mátrixelemet, amely egy-egy detektorpixel becsapódásszám értékét tartalmazza, akkor azok végrehajtása lényegesen lassul az ütközések várható nagy száma miatt. A GPU-n történő optimális végrehajtás érdekében az algoritmust célszerű átszervezni úgynevezett gyűjtő típusú algoritmussá. Az átszervezés lényege, hogy az előrevetítő operátor számításakor az egyes detektorpixelekben előálló projekciós képpontokat úgy számítjuk ki, hogy az egyes detektorpixelekhez rendeljük a számítási szálakat. Ez azt jelenti, hogy a párhuzamosan futó számítási szálak végigmennek azokon a voxeleken, amelyekben lévő aktivitás foton becsapódást okozhat az adott detektorpixelen, és kiszámítják, hogy az egyes voxelek, illetve a bennük elhelyezkedő aktivitás milyen mértékben járul hozzá egy adott detektorpixelben érzékelhető becsapódási esemény számához. A számítás végén a szálhoz rendelt detektorpixelben érzékelt becsapódásszámot kapjuk meg, de minden detektorpixelhez tartozó mátrixelemet csak egy szál fog elérni, frissíteni, ezzel elkerülve a szálak szinkronizációját. A fenti leírásból látható, hogy a kidolgozott megoldás kiküszöböli a szálak gyakori szinkronizációs kényszerét, azonban még nem oldottuk meg azt a problémát, hogy a szálaknak egy időben csak viszonylag kicsi, a szálat futtató maghoz tartozó osztott elérésű

dc_592_12 28

(a) 1.ütem

(b) 2.ütem

(c) 3.ütem

1.11. ábra Az osztott memória hatékony kihasználása érdekében kialakított előrevetítő operátor számításának szemléltetése az első három számítási ütem bemutatásával. A kis képek felső részén látható hisztogramok a projekciós képek számolásának eredményét, az egyes detektorpixelekben mért beütésszám növekedését szemléltetik. Az ábrák alsó részén látható négyzethálós négyzet az aktivitáseloszlást tartalmazó tömböt reprezentálja, a nyilak a projekció számolását.

(shared) memóriában elférő memóriatartalmat kelljen elérni. Erre azért van szükség, mert a GPU-kban a regiszterek mellett két különböző elérési idejű memória található. A minden számítási mag által elérhető központi (main) memória mellett elhelyeztek úgynevezett osztott elérésű (shared) memóriát is, amelyet a számítási magok egy csoportja használ közösen. Az osztott memória mérete ugyan korlátozott (64 kB volt az általunk használt GPU esetén), azonban a központi memóriánál körülbelül két nagyságrenddel gyorsabb az elérése. A gyorsabb végrehajtás érdekében a számítás során az aktivitásokat tartalmazó voxeltömböt olyan részekre bontottuk, amelyek elférnek az egyes magokhoz tartozó osztott memóriában. Az egyes detektorpixelekhez tartozó szálak természetesen dolgozhatnak azonos voxeltömb darabon, így az egyes osztott memóriaegységekhez tartozó számítási magok párhuzamosan futtathatják azokat. Ezt a működést szemlélteti az 1.11. ábra. Az első három számítási ütem látható az ábrán, ahogyan az algoritmus az egyes voxelekben elhelyezkedő aktivitás értékeket tároló mátrix egyes, detektorral párhuzamos rétegeiből származó beütésszámot számítja. A fent leírt módszer alkalmas sugárzás hatásának párhuzamos számítására az osztott memória hatékony kihasználásával. Nem beszéltünk azonban még az elnyelés és a távolságfüggő detektorfelbontás hatásának figyelembe vételéről. Ezen hatások számítását is hozzá


29

kellett hangolni a kialakított új számítási logikához. Mindkét esetben egy-egy optimalizáló lépést építettünk be az algoritmusba, kihasználva a fent kialakított rétegenkénti számolást. Az elnyelés számítását úgy tettük hatékonnyá, hogy minden detektorpixel esetén nemcsak az aktuális pixelben detektált sugárzást (beütésszámot) tartottuk nyilván a számolás során, hanem az adott pixel által látott LOR-on fekvő voxelek akkumulált elnyelési faktorát az éppen feldolgozott voxel rétegig. Így lényegében minden alkalommal fel tudtuk használni a korábbi számításaink eredményét az elnyelő közeg hatásának számításakor. Az egyes voxelekből származó sugárzás detektorpixelen jelentkező hatásának számítása így nagyon leegyszerűsödött, először az adott feldolgozott rétegben levő voxelben elhelyezkedő anyag okozta elnyelés hatásával frissítettük az adott LOR detektorfelület irányába eső szakaszára vonatkozó akkumulált elnyelési értéket, majd kiszámoltuk ez alapján az adott voxelben levő aktivitás okozta sugárzást a detektorpixelben. Ehhez természetesen az akkumulált elnyelési tényezőt használtuk. A távolságfüggő detektorfelbontás hatásának hatékony számítását az tette lehetővé, hogy kihasználtuk azt, hogy az osztott memóriába minden számítási ütemben a detektor felületével párhuzamosan elhelyezkedő voxel réteg adatait töltjük. Így ki tudjuk azt használni, hogy az egyes detektorpixelektől azonos távolságra levő voxelek – korábban említett – elkenődése azonos mértékű, hiszen az csak a voxel és a detektor távolságától függ. Így tehát a minden egyes voxel réteget feldolgozó számítási ütemben azonos – elkenődést jellemző – Gauss eloszlást kell használnunk. Ezt felismerve a számítási ütem megkezdésekor az aktuálisan használt Gauss eloszlást jellemző adatokat eltároljuk – ugyancsak az osztott memóriában –, majd minden detektor voxel (illetve LOR) számításakor újra és újra felhasználhatjuk azt. Az egyes osztott memóriaegységekhez tartozó számítási magok elvileg dolgozhatnak azonos voxeltömb részleten, de különbözőn is. Ezt a voxeltömb, a detektorpixel szám, illetve a használt GPU architektúra függvényében célszerű meghatározni. Az eddigiekben ismertetett, GPU erőforrásainak hatékony kihasználására kidolgozott párhuzamosított SPECT rekonstrukciós algoritmus előrevetítő operátora számításának a sematikus illusztrációját látjuk az 1.12. ábrán. Az ábrán a voxeltömböt, illetve az abban levő aktivitást tároló mátrixot a középső négyzethálós négyzet szemlélteti. Az osztott memória a kis négyzetek sorozatából álló téglalap. Képletesen ebbe kerül bele az aktivitást tároló mátrix egy sora egy adott számítási ütemben, illetve a voxelekben elhelyezkedő aktivitásból származó sugárzás elkenésének számításához használt Gauss eloszlás, amelyet sematikus Gauss eloszlás függvény jelez az ábrán.

dc_592_12 30

1.12. ábra Az előrevetítő operátort kiszámító GPU platformra optimalizált algoritmus menetének szemléltetése. Az ábra közepén szereplő négyzethálós négyzet szemlélteti a rekonstruálandó térrészt reprezentáló adattömböt.

A rekonstruálandó térrészt minden detektor más és más szögből látja. Ahhoz, hogy az egyes processzor csoportok azonos módon feldolgozható adathalmazokat kapjanak, minden detektorhoz ennek az adattömbnek egymáshoz képest transzformált – ha képletesen akarjuk leírni, akkor elforgatott – nézetét kell elkészíteni. Ezt reprezentálja az 1.12. ábrán az (1) lépés. Ezt a transzformációt a kidolgozott algoritmus egy paraméterezett indexelő függvénnyel végzi, amely időről időre átmásolja a mátrixelemeket az elforgatásnak megfelelő pozícióba. Ezután az egyes detektor pozíciókhoz tartozó projekciós képek pixeleit számoló processzor csoportok számára elforgatott mátrixból az algoritmus leszeleteli azt a darabot, ami elfér a processzor csoportok gyors elérésű memória tárába. Ezt a lépést szemlélteti a (2) lépés. Az algoritmus ezt követően a processzor csoport processzorai számára kiosztja a detektor egyes pixeleit, amelyek párhuzamosan számolják a pixelekben becsült beütésszámot a különböző fizikai hatások figyelembe vételével. Amennyiben minden detektorpixelhez tartozó értéket kiszámoltak a processzorok, az algoritmus lecseréli a gyors elérésű memória tartalmát, és megint iterál a detektorpixeleken, párhuzamosan kiosztva azokat a processzorok között. Ezt a folyamatot szemlélteti az 1.12. ábra (3) lépése. Az újraszervezett algoritmus kidolgozásához még egy további lépés hatékonyságán javítottunk. Elvileg a detektorok a voxeltömb által reprezentált térrész körül forognak. Ezen forgás közben az egyes detektorpixelek mindig más és más voxeltömb „részletekre látnak rá”, vagyis a voxeltömb más és más részei alapján határozható meg az adott pixelben


31

1.1. táblázat Az előrevetítő operátort GPU-n hatékonyan számoló algoritmus 01 Inicializáljuk a központi memóriában az aktuális aktivitáseloszlást tároló mátrixot (T ) egyenletes t0 aktivitást helyezve minden voxelbe: tij = t0 , ∀tij ∈ T . 02 Ismételjük a következő lépéseket (i = {1, 2, . . . , K}), ahol K a projekciós képek száma: 03

Kiválasztjuk a következő előállítandó projekciós képet (Pi ).

04

T -t beforgatjuk a Pi -hez tartozó pozícióba (T → Ti ).

05

Ti -t felosztjuk L szeletre úgy, hogy minden szelet elhelyezhető legyen az osztott memóriában és a szeletben a detektortól csak azonos távolságra elhelyezkedő voxelekhez tartozó adatok szerepeljenek.

06

Ismételjük a következő lépéseket (j = {1, 2, . . . , L}), ahol L a szeletek száma:

07

A Ti osztott memóriába elhelyezhető szeletét, Tij -t átmozgatjuk az osztott memóriába.

08

Kiszámoljuk Tij szelethez tartozó Gauss eloszlást (1.7) összefüggés szerint.

09

Az osztott memóriához tartozó magokon futtatjuk az előrevetítő operátor számítását minden Tij -ben levő voxelre.

10

Addig, amíg az egész voxeltömböt fel nem dolgoztuk, vagyis j = L.

11 Addig, amíg az összes projekciós képet elő nem állítottuk, vagyis i = K.

keletkező beütésszám. Annak érdekében, hogy a voxeltömb darabolása során mindig azonos memória rekeszben elhelyezkedő memória részleteket kelljen az egyes osztott elérésű memóriaegységekbe, illetve az egységekhez tartozó szálakhoz – áttételesen pedig a szálakhoz tartozó detektorpixelekhez – másolni, a központi memóriában tárolt voxeltömböt reprezentáló mátrixot időről időre „elforgattuk”, vagyis az elemeit úgy másoltuk át, hogy a különböző detektor pozíciókban az egyes detektorpixelekből látható voxeltömb elemek mindig azonos memória rekeszben legyenek (1.12. ábra). Ezzel kiküszöböltük a voxeltömböt tartalmazó memóriához történő hozzáférés során használt indexek detektor pozíciókhoz tartozó újraszámolását, tovább gyorsítva ezzel az algoritmus végrehajtását. Az előrevetítő operátort számító algoritmus lépéseit az 1.1. táblázat foglalja össze. Az előrevetítő operátorhoz teljesen hasonló módon optimalizáltuk a visszavetítő operátort számoló algoritmust is. A különbség annyi, hogy a „gyűjtő típusú” algoritmust itt úgy definiáltuk, hogy az egyes szálakat a voxeltömb egyes voxeleihez tartozó mátrixelemekhez rendeltük. A visszavetítő operátor számítása során – az előrevetítő operátornál bemutatottakhoz hasonlóan – a távolságfüggő detektorfelbontás kompenzációjához használt, Gauss kernelt reprezentáló mátrix elemeit helyeztük el a gyors elérésű memóriába. A GPU számítási erőforrásait hatékonyan kihasználni képes algoritmus sémáját az 1.13. ábra szemlélteti.

dc_592_12 32

1.13. ábra A visszavetítő operátort kiszámoló, GPU platformra szabott algoritmus szemléltetése. A visszavetítés során az egyes voxelekhez tartozó aktivitásértékek számítását osztottuk ki a processzorok között. A processzor csoportokhoz tartozó gyors elérésű memória részekbe a távolságfüggő detektorfelbontás kompenzációjához használt, Gauss kernelt leíró adatokat töltöttük fel.

1.5.

Az elért eredmények és azok értékelése

A kidolgozott rekonstrukciós algoritmus verifikálását, illetve validációját több lépésben végeztem el. Mind a verifikációt, illetve mind a validációt a nukleáris medicinában a képrekonstrukciós algoritmusok ellenőrzésére bevett gyakorlat szerint hajtottam végre. Így a hagyományos programfejlesztés során bevett verifikációs tesztelés után a megvalósított algoritmusokat először szimulációval előállított matematikai fantomok segítségével teszteltem. A tesztelés során a rekonstruálandó térrészben különböző, rutinszerűen használt geometriai alakzatokba aktivitást definiáltunk és szimulációs szoftverekkel előállítottuk a vetületi képeket, majd ezek rekonstrukcióját hajtottuk végre a kidolgozott algoritmusokkal. Ezután valós SPECT berendezésekkel ismert geometriájú, valós emberi szöveteket és viszonyokat emuláló fantomokról készített vetületi képeket állítottunk elő, és azokon hajtottuk végre a rekonstrukciót. Mindkét esetben az ismert geometriájú aktivitás elhelyezkedését hasonlítottuk össze a rekonstrukció során nyert aktivitás sűrűséggel. A validáció utolsó lépése valós betegekről készített felvételek feldolgozása volt, illetve az eredmények összehasonlítása más rekonstrukciós algoritmusokkal.


33

1.14. ábra Mediso AnyScan kamera.

A következő verifikációs, illetve validációs lépések eredményeit ismertetem: • szimulált matematikai henger fantom vizsgálata, • valósághű mellkas fantom matematikai modelljének vizsgálata, • valósághű mellkas fantomról [74] készült felvételek vizsgálata, • betegekről készült felvételek vizsgálata. A validáció során használt vetületi képeket előállító szimulációkat minden alkalommal a GATE szimulációs toolkit-el [42] végeztük. A szimulátorban definiált rendszer egy átlagos paraméterekkel rendelkező úgynevezett LEHR kollimátoros (Low Energy High Resolution) detektorral rendelkező SPECT berendezés volt, amelynek a paraméterei megegyeznek a Mediso Anyscan kamerájának paramétereivel. A szimulációk során, illetve a méréseknél használt izotóp mindig Tc99m volt. A szimulációk paraméter beállításait az 1.2. táblázat foglalja össze. A GATE szimulációk nagy számításigénye miatt azokat grid-en futtattuk. Erre a Hungrid [4] és a biomed VO [3] rendszereket használtuk. A valós mérések (ide értve az anatómiai fantom és beteg felvételeket is) a SOTE Radiológiai és Onkoterápiás Klinikáján [2] készültek a Mediso Kft. AnyScan SPECT/CT kamerájával (1.14. ábra). A mérések előtt elvégeztük a rendszer – a távolságfüggő detektorválasz kompenzációhoz szükséges – kalibrációs pontválasz méréseit mind a valós kamera, mind a szimulált SPECT berendezés esetében.

dc_592_12 34

1.2. táblázat A validáció során használt képek paraméterei Henger fantom szimuláció psfa

2.33 mm

psfb

0.033

RoR (Radius of Rotation)

281 mm

adatgyűjtés módja

128 db 128x128-as projekciós kép

szkennelt szögtartomány

360 fok

összbeütés-szám

9.000.000 db, kb. 70.000 db/kép NACT szimuláció

psfa

2.33 mm

psfb

0.033


281 mm




360 fok

összbeütés-szám

11.400.000 db, kb. 80.000 db/kép Biodex fantom (Anyscan LEHR kollimátorral)

psfa

2.4 mm

psfb

0.0313


356 mm




360 fok

összbeütés-szám

17.700.000 db, kb. 140.000 db/kép Beteg felvétel (Anyscan LEHR kollimátorral)

psfa

2.4 mm

psfb

0.0313


356 mm




360 fok

összbeütés-szám

6.200.000 db, kb. 97.000 db/kép


35

1.15. ábra Henger alakú aktivitást tartalmazó fantom, amely inhomogén elnyelési állandójú közeggel van kitöltve. Az ábra a fantom keresztmetszeti képe, amely aktivitás és elnyelési térkép is egyben.

GATE szimulációval előállított geometriai fantom vizsgálata A kidolgozott algoritmus verifikációjának első lépéseként arra voltunk kíváncsiak, hogy mennyire hatékonyan sikerül az új algoritmusnak a célul megjelölt fizikai jelenségek hatását kompenzálni. Ezért szabályos geometriai elemeket tartalmazó fantomot definiáltunk. A fantomban az aktivitást egy viszonylag keskeny falú henger falába definiáltuk. A megfigyelt térrészben két különböző elnyelési állandóval rendelkező anyaggal kitöltött hengert helyeztünk el. A hengerek fala párhuzamos volt, azonban a tengelyük egymástól el volt tolva. A fantom metszeti képén jól látható a gyűrű alakú aktivitás a két különböző elnyelési állandójú elnyelő közegben (1.15. ábra). A fantom definiálása után elvégeztük a SPECT képalkotás szimulációját a GATE szimulációs tool segítségével. A szimuláció során összegyűjtöttük az úgynevezett single eventeket. A szimuláció „Radius of Rotation” paramétere 281mm volt. A single event-ek közül eltávolítottuk a Compton szórt eseményeket. Ez lényegében megfelel egy ideális szóráskorrekciónak. Ezzel elértük, hogy a szimulált mérések eredményét a szórás nem torzítja, így a rekonstruált kép minősége a szórás jelenségétől nem, hanem elsősorban az általunk kidolgozott algoritmus hatékonyságától fog függeni. Végül a single event-eket levetítettük 128 db 128x128-as projekciós képre (pixelméret 3mm), 360 fokos szögtartományra. A projekciós képeket rekonstruáltuk az új GPU alapú OSEM algoritmussal, kompenzálva a gamma foton elnyelést és a detektor távolságfüggő felbontását. Az OSEM rekonstrukció során nyolc subset-et használtunk és húsz iterációt hajtottunk végre. A rekonstruált szeletek az 1.16. ábrán láthatók. A Mediso Cardiac Stress/Rest szoftver [60] lehetőséget ad egy adott középpontból a

dc_592_12 36

1.16. ábra A henger fantom három egymás melletti transzaxiális szelete. Felső sor: referencia, az eredeti fantom képe; középső sor: 2D-OSEM algoritmussal rekonstruált szeletek; alsó sor: új GPU alapú 3D-OSEM algoritmussal rekonstruált szeletek elnyelési korrekcióval és detektorválasz kompenzációval. A rekonstrukció eredménye szemléletesen igazolja, hogy a kidolgozott algoritmus hatékonyan képes csökkenteni a képalkotás során fellépő fizikai hatások okozta – korábbi algoritmusok által nem kiküszöbölt – képminőség romlást amellett, hogy a párhuzamosítás eredményeként futási ideje nem növekedett.

tér egyes irányaiban mért aktivitás ábrázolására, illetve két kiválasztott képen mérhető aktivitás összehasonlítására. Ennek a segítségével összehasonlítottuk a hagyományos kétdimenziós, fizikai hatások kompenzációját nem tartalmazó (2D-OSEM) rekonstrukciót az általunk fejlesztett új (3D-OSEM) algoritmussal a két rekonstrukció eredményeként kapott képek alapján. Az összehasonlítás alapjául választott 2D-OSEM algoritmus körülbelül azonos teljesítményű a jelenlegi orvosi gyakorlatban alkalmazott rekonstrukciós algoritmusokkal. A Cardiac Stress/Rest szoftver által generált, szívvizsgálatoknál gyakran alkalmazott úgynevezett bullseye megjelenítést – amely a szívfalnak megfelelő csonka kúp palástját kivetíti egy kör felszínére – láthatjuk az 1.17. ábrán. Ezeken az ábrákon jól látható a rekonstruált aktivitáseloszlás. Az 1.16. és az 1.17. ábrákon látható, hogy az algoritmus a gamma foton elnyelődésből adódó virtuális aktivitáskiesést megszüntette, a metszeti képen a teljes gyűrűt visszaállította. A detektor távolságfüggő elkenéséből adódó hatást, amely a henger falának vastagodásaként jelentkezett és elsősorban a kis elnyelési együtthatójú térrészben volt jól megfigyelhető, szintén korrigálta az algoritmus: az aktivitást tartalmazó gyűrű falvastagsága közel azonos a referencia képen látható vastagsággal.


(a) Referencia

37

(b) 2D-OSEM rekon- (c) 3D-OSEM rekonstrukció eredménye strukció eredménye elnyelési korrekcióval és detektorválasz kompenzációval

1.17. ábra Henger alakú aktivitást tartalmazó geometriai fantom rekonstruált képe bullseye reprezentációban.

GATE szimulációval előállított NCAT fantom vizsgálata

A verifikáció második lépésében az emberi mellkas MRI képek alapján készített matematikai fantomját, az úgynevezett NCAT fantomot [74] használtuk. A fantom segítségével egy valós vizsgálatot, a myokardiális perfúziós vizsgálatot szimuláltuk. Ennek megfelelően az NCAT fantomban aktivitást definiáltunk a myokardiumba és a májba, mivel ezen vizsgálatoknál az alkalmazott radiopharmakon egy része a májban is szokott dúsulni. A szimuláció során a Radius of Rotation 281mm volt, a Compton szórt eseményeket – az előző kísérletben leírt indokok alapján – kiszűrtük a single event-ek közül, ideális szóráskorrekciót végezve ezzel. Végül a single event-eket levetítettük 128 db 128x128-as projekciós képre (pixelméret 3mm), 360 fokos szögtartományra. A projekciós képek rekonstrukciója hasonló módon történt, mint az előző kísérletnél, a képeket rekonstruáltuk mind 2D-OSEM algoritmussal, mind az új GPU alapú OSEM algoritmussal, kompenzálva a gamma foton elnyelést és a detektor távolságfüggő felbontását. Az OSEM rekonstrukció során nyolc subset-et használtunk és húsz iterációt hajtottunk végre, ami megfelel egy átlagos szívvizsgálat paramétereinek. Az 1.18. ábrán egy transzaxiális referencia szelet és a megfelelő rekonstruált szelet látható. A két rekonstrukció eredményét összehasonlítottuk a Cardiac Stress/Rest szoftverkörnyezet segítségével. A rekonstruált szív különböző irányú metszetei az 1.19. ábrán láthatók, míg a bullseye képeket az 1.20. ábra mutatja. Az új GPU alapú 3D-OSEM algoritmussal rekonstruált képek kontrasztosabbak, a bennük levő aktivitáseloszlás egyenletesebb.

dc_592_12 38

(a) Referencia

(b) Rekonstruált szelet

1.18. ábra NCAT fantom referencia metszeti képe (1.18(a)) és rekonstruált megfelelője 3D-OSEM elnyelési korrekcióval és detektorválasz kompenzációval (1.18(b)).

1.19. ábra Az NCAT fantom rekonstruált képeinek különböző metszetei. Mindhárom sorozatban a felső sor (AC-Co jelöléssel) az általunk kidolgozott 3D-OSEM algoritmussal rekonstruált, míg az alsó sor (NonAC felirattal) 2D-OSEM algoritmussal rekonstruált képek metszetei. A képeket összevetve látható a kontraszt javulása, és a homogénebb aktivitáseloszlás, amellett, hogy az algoritmus futási ideje csökkent a párhuzamosítás eredményeként.


39

1.20. ábra NCAT fantom rekonstruált képeinek bullseye megjelenítése. Balra fent a 3DOSEM algoritmus, míg alatta a 2D-OSEM algoritmus által rekonstruált kép. A két rekonstrukció által előállított aktivitáseloszlás különbsége látható a két jobboldali ábrán.

Mellkas fantom mérés SPECT/CT kamerával A verifikáció következő lépésében a kereskedelmi forgalomban kapható Biodex Lung Spine SPECT fantomot [23] használtuk. A Biodex fantom az emberi mellkas anatómiailag helyes modellje. A lapított henger alakú, mellkast formázó tartályban találhatók a belső szerveket modellező úgynevezett tüdő és szív betétek. A tüdő betét két tüdő formájú rekesz, amely hungarocellgolyókkal van feltöltve. Ennek elnyelési együtthatója közel azonos a tüdő szövettel. A tüdő betét mellett található szívpatkó alakú rekesz a szív bal kamráját modellezi. A szív betét „myokardium” része, illetve a belső rekesz is feltölthető aktivitással. A mérések során a Biodex Lung Spine fantom myokardium rekeszét aktivitással töltöttük fel, így készítettük a fantomról a SPECT felvételeket. A SPECT vizsgálat során 360 fokos szkennelési tartományban 128 db 128x128 pixeles projekciós képet gyűjtöttünk (pixelméret: 4mm). A gyűjtési idő projekciós képenként 20 másodperc volt, amely átlagosan kb. 140.000 eseményt indukált projekciós képenként. A projekciós képeket – a korábbi vizsgálatokhoz hasonlóan – rekonstruáltuk mind 2DOSEM algoritmussal, mind az új GPU alapú OSEM algoritmussal, kompenzálva a gamma foton elnyelést és a detektor távolságfüggő felbontását. A rekonstrukciót a korábbiakhoz hasonló paraméterekkel végeztük: nyolc subset-et használtunk az OSEM algoritmusnál és 20 iterációt hajtottunk végre.

dc_592_12 40

1.21. ábra Biodex Lung-Spine SPECT mellkas fantom rekonstruált szeleteinek különböző irányú metszetei. Mindhárom sorban a felső képek a 3D-OSEM, míg az alsó képek a 2DOSEM rekonstrukció eredményét mutatják.

A két rekonstrukciót összehasonlítottuk a már ismertetett Cardiac Stress/Rest szoftverkörnyezet segítségével. A rekonstruált szív különböző irányú metszetei az 1.21. ábrán láthatók, míg a bullseye képeket az 1.22. ábra mutatja. Az új GPU alapú 3D-OSEM algoritmussal rekonstruált képek kívülálló számára is jól látható módon kontrasztosabb eredményt szolgáltattak, mint a 2D-OSEM algoritmussal előállított képek. Orvosdiagnosztikai szempontból fontos jelenség a 2D-OSEM algoritmussal rekonstruált képek bulsseye megjelenítésén (1.22. ábra) látható úgynevezett hátsó szívfal kiesés, amely a képeken kisebb aktivitású foltként jelenik meg – tévesen – azt az információt adva, hogy a szívfal ezen a területen elvékonyodott, illetve eltűnt. Ez tipikusan az inhomogén közeg okozta gamma foton elnyelés hatása. Az új GPU alapú 3D-OSEM algoritmus által rekonstruált képeken ez a hatás már lényegesen csökkent, részben eltűnt. A szív apikális részén levő hasonló jellegű csökkent aktivitás az új algoritmus által rekonstruált képen is részben megmaradt. Ez a szívcsúcs közelében látható aktivitás csökkenés – az úgynevezett szívcsúcs-kiesés – ismert hiányossága a berendezés által készített felvételeknek.


41

1.22. ábra Biodex Lung-Spine SPECT mellkas fantom rekonstruált képei bullseye megjelenítésben. Balra fent a 3D-OSEM algoritmus, míg alatta a 2D-OSEM algoritmus által rekonstruált kép. Mindkét rekonstrukció esetén az eredeti fantomban definiált és a rekonstrukció által előállított aktivitáseloszlás különbsége látható a két kép jobboldalán.

Ezt a jelenséget nagy valószínűséggel a Partial Volume Effect, valamint a szórt események nagy száma okozza. Ennek elemzéséhez további mérések szükségesek. Betegekről készített SPECT/CT felvételek A validáció következő lépésében egy női beteg szívéről készített SPECT felvételt rekonstruáltuk – a korábbiakhoz hasonlóan – két rekonstrukciós algoritmussal. Az AnyScan SPECT/CT berendezés beállításai a következők voltak: 360 fokos szkennelési tartomány; 64 db 128x128 pixeles projekciós kép; pixelméret: 4mm. A gyűjtési idő projekciós képenként 30 másodperc volt, amely átlagosan kb. 97.000 eseményt indukált projekciós képenként. A projekciós képeket rekonstruáltuk 2D-OSEM algoritmussal, valamint az új GPU alapú OSEM algoritmussal, kompenzálva a gamma foton elnyelést és a detektor távolságfüggő felbontását. Az OSEM rekonstrukció során nyolc subset-et használtunk és 10 iterációt hajtottunk végre. A két rekonstrukciót összehasonlítottuk a Cardiac Stress/Rest szoftverkörnyezettel. A rekonstruált szív különböző irányú metszetei az 1.23. ábrán láthatók, míg a bullseye-okat az 1.24. ábra mutatja. Az új GPU alapú 3D-OSEM algoritmussal rekonstruált képek szabad szemmel láthatóan kontrasztosabb eredményt szolgáltattak.

dc_592_12 42

1.23. ábra Valós szív felvétel rekonstrukciójának különböző orientációjú metszeti képei. Felső sorok a 3D-OSEM, míg az alsó sorok a 2D-OSEM rekonstrukció eredményét mutatják.

1.24. ábra Valós szív felvétel rekonstrukciójának bullseye megjelenítése. Balra fent a 3DOSEM algoritmus, míg alatta a 2D-OSEM algoritmus által rekonstruált kép.


43

1.25. ábra Agyról készült SPECT felvétel 2D-OSEM algoritmussal végzett rekonstrukciójának eredménye InterViewXP alkalmazással [60] megjelenítve.

Agyról készült SPECT/CT felvételek rekonstrukciója Agyi elváltozások vizsgálatára rendszeresen alkalmazzák a multi-modalitású SPECT/CT berendezéseket. Vizsgálatunkhoz az AnyScan készüléket használtuk dual modalitású üzemmódban egy férfi páciens agyának vizsgálatára. A rekonstrukció eredményét a klinikumban szokásos megjelenítési módszerek szerint mutatjuk be. Az 1.25. ábra a beteg agyáról készült felvétel 2D-OSEM algoritmussal rekonstruált változatát, míg az 1.26. ábra az általunk kidolgozott új GPU alapú 3D-OSEM algoritmussal előállított rekonstrukció eredményét mutatja. A 3D-OSEM eljárással kapott képek élesebbek, kontrasztosabbak. A 2D-OSEM rekonstrukcióhoz tartozó képeken látható véletlen háttérzajként jelentkező elmosódott felületek a 3D-OSEM rekonstrukció során szinte eltűntek. Párhuzamos környezetre szabott algoritmus hatékonyságának értékelése A GPU környezetre szabott párhuzamos SPECT rekonstrukciós algoritmus hatékonyságát három benchmark feladat végrehajtásával értékeltük ki, amelyek a következők: A(128) Egy 128 × 128 × 128 voxelt tartalmazó térrész rekonstrukciója 128 projekciós kép alapján. Ez felel meg egy ma, a gyakorlatban használt átlagos képrekonstrukciós feladatnak.

dc_592_12 44

1.26. ábra Agyról készült SPECT felvétel 3D-OSEM algoritmussal végzett rekonstrukciójának eredménye InterViewXP alkalmazással [60] megjelenítve.

B(256) Egy 256×256×256 voxelt tartalmazó térrész rekonstrukciója 256 projekciós kép alapján. Ez a benchmark méret már egy ma elérhető nagyfelbontású képalkotó eszközben szükséges képrekonstrukciónak felel meg. C(384) Egy 384 × 384 × 384 voxelt tartalmazó térrész rekonstrukciója 256 projekciós kép alapján. A képrekonstrukciós eszközök továbbfejlesztésével a belátható jövőben ilyen méretű képrekonstrukciós feladatok lehetségesek. A benchmark rekonstrukciós feladatok paramétereit az 1.3. táblázat foglalja össze. A rekonstrukciós tesztjeinket négy különböző GPU kártyán hajtottuk végre: • 9600-GT, • GTX-550-Ti, • GTX-480, • GTX-580.

dc_592_12 45


1.3. táblázat A benchmark feladatok paraméterei Paraméter

A(128)

B(256)

C(384)

Voxel tér mérete

128 × 128 × 128

256 × 256 × 256

384 × 384 × 384

256

256

1282

2562

3842

Iterációk száma (db)

1

1

1

Számítások száma1

1283 × 128 × 2563 × 256 × 202 = 1,07×1011 402 = 6,87×1012

3843 × 256 × 602 = 5,22×1013

Voxel méret (mm)

4

2

1

psfa (mm)

1,73

1,73

1,73

psfb

0,0358

0,0358

0,0358

RoR (mm)

281

281

281

Kollimátor típusa

AnyScan/LEHR

AnyScan/LEHR

AnyScan/LEHR

Projekciós (db)

képek

Projekciós képek száma (db)

1

száma 128 pixel-

Minden voxelt levetítünk minden vetületre, és ott egy Gauss kernel alapján elosztjuk a beütések

számát – így a számítások száma körülbelül a voxelek száma × vetületek száma × Gauss kernel mérete. A Gauss kernel egy 2D mátrix, amely oszlopainak és sorainak száma a detektorfelülettől mért legnagyobb előforduló voxel távolságnál ( = átló fele + RoR) számolt szórás (σ) 3,03-szorosa.

dc_592_12 46 1.4. táblázat A grafikus kártyák legfontosabb paraméterei Paraméter

9600-GT GTX-550 GTX-480 GTX-580

CUDA magok száma

64

192

480

512

Órajel (MHz)

650

900

700

777

Elméleti számítási teljesítmény (GFlops, single-precision)

208

691

1345

1581

Memóriabusz szélessége (bit)

256

192

384

384

Memória sebessége (GB/s)

58

98

177

192

1.5. táblázat A párhuzamos és a hagyományos algoritmus futási idejének összehasonlítása CPU Magok száma

1

9600GT GTX-580 64

512

Egy előrevetítés Futási idő (s)

2 068,86

9,34

1,90

Gyorsulás1 (%)

100

0,45

0,09

Teljes rekonstrukció (20 iteráció)

1

Futási idő (s)

21 964 43

582,16

88,60

Gyorsulás1 (%)

100

2,65

0,40

A CPU-n futtatott algoritmus idejéhez viszonyítva.

A GPU kártyák jellemző paramétereit az 1.4. táblázat foglalja össze. A kártyák közötti legfontosabb különbség a végrehajtó magok számának növekedése, azonban az egyre modernebb modellekben az egyéb hardver jellemzők is javultak. A benchmark kísérletek során elsősorban azt fogjuk megvizsgálni, hogy a nagyobb kártyák kínálta többleterőforrásokat a párhuzamos algoritmus ki tudja-e használni, illetve azt, hogy az általunk javasolt, hardver erőforrásokat hatékonyan kihasználó módszer mekkora gyorsító hatással bír. Annak érdekében, hogy demonstrálni tudjuk a kidolgozott megoldás gyorsító hatását legelőször lefuttattuk a rekonstrukció előrevetítő operátorának kiszámolását, illetve egy 20 iterációt tartalmazó teljes rekonstrukciót hagyományos CPU-t tartalmazó számítógépen, valamint a GPU-ra szabott algoritmusunkat GPU környezetben. Minden esetben az A(128)-as benchmark feladat elvégzése volt a teszt. Az eredmények az 1.5. táblázatban láthatók.


47

A hagyományos algoritmus C++ környezetben készült. Az implementált kód egy jelenleg általánosan használt MLEM rekonstrukciós algoritmust valósít meg, amelyet kiegészítettünk az általunk megcélzott két képminőség javító lépéssel, az elnyelési korrekcióval és a távolságfüggő detektorfelbontás kompenzációjával. Az egy szálon futtatott rekonstrukció során is implementáltuk az adott környezetben lehetséges számítási lépésszámot csökkentő megoldásokat. Ennek hatását láthatjuk a táblázat CPU-s környezetre vonatkozó oszlopában, ahol az öt iterációt tartalmazó rekonstrukció lényegesen kevesebb idő alatt futott le, mint az előrevetítés ötszöröse, pedig a rekonstrukció tartalmaz öt előre- és öt visszavetítést. Ez azért van így, mert az előrevetítésnél a CPU algoritmus nem számolja ki azon voxelekből származó hatásokat, amelyekben nem helyezkedik el aktivitás. A GPU-ra szabott és a CPU-n futó algoritmus futási időinek összehasonlításánál látható a gyorsulás mértéke. A 64 magot tartalmazó GPU esetén az előrevetítés több, mint 200-szorosára, a teljes rekonstrukció több, mint 38-szorosára gyorsult. Ugyanilyen sorrendben a gyorsulás aránya az 512 magot tartalmazó processzornál 1000-szeres, illetve közel 250-szeres. Ez alapján határozottan kijelenthetjük, hogy a kidolgozott módszer képes dinamikusan kihasználni a GPU-k kínálta egyre nagyobb számítási erőforrásokat. A gyorsulás mértéke a processzorszám növekedésénél nagyobb mértékű, ami részben alátámasztja az előző azon állítást, hogy a kidolgozott algoritmus hatékonyan képes kihasználni a GPU erőforrásait. A futási idő abszolút értéke is figyelemre méltó. Egy teljes, 20 iterációt tartalmazó rekonstrukció valamivel kevesebb, mint másfél perc alatt futott le az A(128)-as benchmark esetén. Mint azt említettem, a B(256) rekonstrukciós benchmark mérete egyezik meg körülbelül a ma gyakorlatban alkalmazott legnagyobb felbontású eszközökben elvégzett rekonstrukciós feladatok méretével. Ezzel a benchmarkkal is végeztünk tesztfuttatásokat. Az algoritmus segítségével lehetséges volt egy előrevetítő iterációt valamivel több, mint 1 perc alatt végrehajtani mind a GTX-480-as (75,287 sec), mind a GTX-580-as (62,908 sec) kártyán. Ez azt jelenti, hogy egy gyakorlatban ma használt legnagyobb felbontású eszközben alkalmazott rekonstrukció 20 iterációja végrehajtható kb. fél óra alatt, amely lehetővé teszi a kidolgozott rekonstrukciós algoritmus klinikai használatát. A B(256) benchmark futási eredményeit a későbbiekben még elemezzük, azokat az 1.6. táblázat tartalmazza. A hatékony hardver erőforrás kihasználás másik fontos tesztje, ha különböző számú számítási magot tartalmazó eszközökön futtatjuk az algoritmusunkat, és megnézzük, hogyan alakul a futási ideje egy adott feladat kiszámításának. Mi ezt az A(128), illetve a B(256) benchmark feladatok különböző GPU-kon történő végrehajtásával tettük meg. A futási eredmények az 1.6. táblázatban láthatók.

dc_592_12 48

1.6. táblázat Elért és várt gyorsulás összehasonlítása 9600GT GTX-550 GTX-480 GTX-580 Magok száma

64

192

480

512

Magok aránya1

1

3

7,5

8

100

33,3

13,3

12,5

Várható gyorsulás2 (%)

A(128) benchmark GPU futási idő3 (ms) Elért gyorsulás4 (%)

9 941

3 974

2 039

1 496

100

39,98

20,51

15,05

B(256) benchmark GPU futási idő (ms) Elért gyorsulás4 (%)

441 641

203 292

75 287

62 908

100

46,03

17,05

14,24

1

A 9600GT magok számához képest.

2

A magok számának növekedésével arányos gyorsulás a futási idők arányával jellemezve.

3

A program betöltés közel állandó idejét – a könnyebb összehasonlítás érdekében – levontuk a teljes futási

időből. 4

A 9600GT kártyán futtatott algoritmus idejéhez viszonyítva.

dc_592_12 49

Párhuzamos SPECT rekonstrukció 1.7. táblázat A három benchmark futási idejének összehasonlítása1

Számítások száma Számítások aránya2

A (128)

B (256)

C (384)

1,07 × 1011

6,87 × 1012

2,61 × 1013

1

64,21

487,85

Osztott memóriát kihasználó gyorsítás nélkül Futási idő (ms) Futási idő aránya A(128)-hoz3

10 328

479 141

4 661 644

1

46,39

451,36

Osztott memóriát kihasználó gyorsítással Futási idő (ms) Futási idő aránya A(128)-hoz3 Gyorsítás a javasolt módszerrel4 (%)

7 751

154 308

1 601 057

1

19,91

206,56

75,05

32,21

34,35

1

GTX-580-as GPU kártyán végzett futtatások.

2

Az A(128) benchmark számítási műveleteinek számához viszonyítva.

3

Futási idő aránya az A(128) benchmark futási idejéhez viszonyítva.

4

Gyorsítás mértéke (futási idő aránya) az osztott memória használat nélküli változathoz viszonyítva.

Az eredményeket megvizsgálva egyértelműen látszik, hogy a processzorszám növekedésével egyenes arányban csökken az algoritmus futási ideje. Ez egyrészt megerősíti, hogy a kidolgozott algoritmus képes volt egyre több szálon futni, kihasználva az egyre több rendelkezésre álló processzort, másrészt mutatja, hogy a hatékony erőforrás kihasználást eredményező módszer processzorszámtól függetlenül alkalmazható. Végezetül bemutatom a három benchmark futási idejének összehasonlítását GTX-580as GPU kártyán (1.7. táblázat) végrehajtva. A táblázat tartalmazza az egyes benchmark feladatok végrehajtásához szükséges számítások becsült számát, azok növekedési arányát a legkisebb benchmarkhoz (A(128)) viszonyítva. A táblázat tartalmazza egy darab előrevetítő operátor kiszámításához szükséges időket a kidolgozott algoritmussal, mindhárom benchmark esetén. Annak érdekében, hogy demonstráljuk, az osztott memória hatékony kihasználását lehetővé tevő, általunk kidolgozott megoldás hatékonyságát, az algoritmusunkat lefuttattuk ezen megoldás kikapcsolásával is. A táblázatban ezen kívül szerepel egyrészt az osztott memóriával, illetve az osztott memória nélkül született, különböző benchmarkok esetén mért futási idők egymással történő összehasonlítása, másrészt a benchmarkok futási idejének egymással történő összevetése.

dc_592_12 50 A futási időket megvizsgálva megállapíthatjuk, hogy mindhárom benchmark esetén egyértelműen látható a gyorsulás az osztott memória hatékony kihasználását lehetővé tevő megoldás hatására, valamint az is látható, hogy a futási idő a nagyobb benchmarkok esetén kevésbé növekedett, mint a feladat mérete: a B(256) esetén kb. 64-szeres feladat méret növekedés ellenére a futási idő csak kb. 46, illetve kb. 20-szorosára emelkedett; a C(384) esetén közel 500-szoros feladat növekedés esetén a futási idők kb. 450-szeres, illetve a javasolt módszer alkalmazásával 206-szoros növekedést mutatnak. Az 1.7. táblázatban azt látjuk, hogy a javasolt módszer gyorsító hatása függ a feladat méretétől. Míg a futási idő az A(128) benchmark esetén a javasolt gyorsító módszer alkalmazásakor a 75%-ára, úgy a nagyobb benchmarkok (B(256) és C(384)) esetén a futási idő közel egyharmadára csökkent. Ennek a viselkedésnek oka a következő. A javasolt módszer tartalmaz egy olyan többletadminisztrációs lépést – a voxelekben elhelyezkedő aktivitás értékeket reprezentáló mátrix forgatását az egyes detektorpozícióknak megfelelően –, ami nem szerepel az eredeti algoritmusban, és növeli annak futási idejét. Vizsgálataink alapján ezen lépés végrehajtási ideje az eredeti algoritmus futási idejének akár 25%-át is kiteheti és nagyjából arányos a projekciós képek számával. A megfigyelt feladatmérettől függő gyorsulási mérték növekedés legfontosabb oka tehát, hogy a nagyobb méretű feladatok esetén a javasolt gyorsítási módszer jobban tudja kompenzálni ennek a többletadminisztrációs lépésnek a hatását. Azt azonban meg kell jegyezni, hogy a processzorszám növekedés miatt várt gyorsulási arányt meghaladó gyorsulás okozta futási idő csökkenést elképzelhető, hogy részben a processzorszámtól eltérő hardver jellemzők (pl. memória elérés idejének csökkenése) okozta.

1.6.

Összefoglalás

A fejezetben ismertettem az általam kidolgozott, SPECT eljárással készített képek rekonstrukciójára alkalmazható párhuzamos képrekonstrukciós algoritmust. Az algoritmus iteratív rekonstrukciót hajt végre, a Maximum Likelihood Expectation Maximization (MLEM) elvet követve. Annak érdekében, hogy az algoritmus hatékonyan tudja kihasználni a SIMD típusú GPU hardver erőforrásait, a projekciós képeket előállító előrevetítő operátor, illetve az aktivitáseloszlást módosító visszavetítő operátor számításának hagyományos menetét megváltoztattam. Nem a voxeltömbben tárolt aktivitás értékeket reprezentáló mátrixelemeken történik az iteráció és annak számítása, hogy azok milyen hatást váltanak ki az egyes detektorpixeleken, hanem a detektorpixeleken iterálva az egyes detektorpixelekhez történik a beütésszám számítása. Ez a számítás több detektorpixelre párhuzamosan kerül végrehajtásra a GPU processzorain.


51

A gyors végrehajtás érdekében a hardver lehetőségeit figyelembe véve az algoritmus lépéseit módosítottam, illetve átszerveztem. Ezek a módosítások elsősorban a processzorokon futó szálak által gyorsan elérhető, osztott elérésű memória hatékony kihasználását teszik lehetővé. Az új algoritmus által megvalósított képrekonstrukciós eljárást validáltam, hogy igazoljam a beépített képminőség javító módszerek hatékonyságát. A képrekonstrukciós algoritmus – a képminőséget negatívan befolyásoló fizikai jelenségek hatásának csökkentése eredményeképpen – a korábban használt képrekonstrukciós algoritmusoknál jobb minőségű képet eredményez. Az eredmények elemzése során tesztfeladatok végrehajtásával megvizsgáltam mennyire képes kihasználni az algoritmus a GPU biztosította hardver erőforrásokat. Igazoltam, hogy az algoritmus képes a GPU adta hardver erőforrások hatékony kihasználására, ezáltal sikerült a rekonstrukció idejét az orvosi gyakorlat által támasztott követelményeknek megfelelő határok alá csökkenteni. Továbbfejlesztési lehetőségek A kutatást több irányban is lehetséges folytatni. A hatékony erőforrás kihasználást lehetővé tevő algoritmussal olyan mértékben sikerült csökkenteni a képrekonstrukció futási idejét, hogy abba további képjavító módszerek is beépíthetők anélkül, hogy az algoritmus végrehajtási ideje meghaladná a klinikai alkalmazás által elfogadható időkorlátot. Itt kézenfekvő opció annak vizsgálata, hogy a szóráskorrekciós módszerek miként építhetők be az algoritmusba. Ezenkívül lehetséges az ismertetett módszert más képi modalitások esetén is alkalmazni, itt a Pozitron Emissziós Tomográfia (PET) lehet a kutatás elsődleges célpontja.

dc_592_12 52

1.7.

A fejezetben ismertetett eredmények tézisszerű megfogalmazása

1. tézis: Iteratív SPECT (Single Photon Emission Computed Tomography – egy-fotonos számítógépes tomográfia) képrekonstrukció párhuzamos megvalósítása során alkalmazható módszereket dolgoztam ki, amelyek biztosítják a rekonstrukciót végző algoritmus SIMD (Single Instruction, Multiple Data) típusú GPU hardver eszközön történő végrehajtásakor a rendelkezésre álló számítási erőforrások hatékony kihasználását. A javasolt módszerek alkalmazásával sikerült a létrehozott rekonstrukciós algoritmus futási idejét az orvosi gyakorlat által támasztott követelményeknek megfelelő határok alá csökkenteni. A párhuzamos, képminőség javító eljárásokat tartalmazó algoritmust implementáltam, majd számos, szimulációval, illetve valós méréssel előállított rekonstrukciós feladaton validáltam, és elemeztem annak futási idejét. Bebizonyosodott, hogy a kidolgozott algoritmus – a képminőséget negatívan befolyásoló fizikai jelenségek hatásának csökkentése eredményeképpen – a korábban használt képrekonstrukciós algoritmusoknál nem csupán gyorsabb, de jobb minőségű képet eredményez, így lehetővé teszi az orvosi diagnosztikai munka minőségének jelentős javítását.

1.7.1.

A fejezet eredményeit bemutató publikációk

A fejezetben ismertetett eredményeket bemutató fontosabb publikációk: [91], [90], [95], [93], [92], [94], [14], [44], [7].

dc_592_12

2. fejezet A foggyökér csatorna geometriájának CT felvételek alapján történő automatikus meghatározására szolgáló eljárások kidolgozása A tartósan sikeres fogászati beavatkozás elvégzéséhez igen fontos információ a kezelendő fog geometriájának, elsősorban a gyökércsatorna elhelyezkedésének, illetve alakjának pontos ismerete. A fog belső szerkezetének megismerése a kezelés megkezdése előtt azért lényeges, mert segítséget ad a kezelés megfelelő megtervezéséhez. A kezelés menete és időtartama igen eltérő lehet egy egyenes, illetve egy többszörös görbülettel rendelkező gyökércsatorna esetén. Sajnos a fogak gyökércsatornája – anatómiai helyzetüktől függően – igen nagy változatosságot mutat [31]. Mára a fogkonzerváló központok által megfizethető árkategóriába kerültek olyan új, 3D képalkotó technológiát alkalmazó berendezések, mint amilyen például a cone beam computer tomográfiás (CBCT) berendezés. Ezek a berendezések megteremtik a lehetőségét a fogak belső szerkezetének meghatározására [69, 81, 73]. A CT felvételeket automatikusan feldolgozó algoritmusok kidolgozása azonban komoly kihívást jelent. A berendezések által alkalmazható sugárdózis szigorúan szabályozott [35], más berendezésekkel összehasonlítva azoknál általában alacsonyabb, így az általuk készített felvételek jel-zaj viszonya is rosszabb. Ebben a fejezetben a fogászati beavatkozások támogatására, azok hatékonyságának növelése érdekében két, egymástól eltérő, automatikusan végrehajtható eljárást mutatok be a gyökércsatorna meghatározására és annak leírására, két különböző berendezéssel készített

53

dc_592_12 54 CT felvételek alapján. Az első berendezés egy micro-CT, a másik egy cone beam CT. A micro-CT előnye, hogy más CT berendezésekkel összehasonlítva viszonylag nagy felbontású képet szolgáltat, emiatt intenzíven használják fogászati témájú kutatások során. A cone beam CT berendezés választását az indokolja, hogy ez az a berendezés, ami várhatóan a leggyorsabban fog elterjedni a fogászati gyakorlatban, mert ennek ára a legalacsonyabb a hasonló képességű berendezések között. A kidolgozott eljárások gyakorlatba ültetése jelentős változást hozhat a fogászati beavatkozások során. A javasolt hatékony automatikus eljárások lényegesen egyszerűsíthetik a kezelés megtervezését, illetve kiszámíthatóvá tehetik, adott esetben csökkenthetik a kezelés időtartamát, amely előnyös mind a beavatkozást végző szakember, mind pedig a páciens számára. További előnye az automatikus végrehajtást lehetővé tevő eljárásoknak, hogy azokat nagy tömegű CT felvételen végrehajtva pontosabb képet kaphatunk az egyes fogak jellemző foggyökér csatorna felépítéséről, illetve a lehetséges extrém csatorna alakokról, amely igen hasznos információ lehet a fogászati eszközök tökéletesítése során.

2.1.


A különféle fogak gyökerei változatosak lehetnek, de két páciens anatómiailag egymásnak megfelelő fogának gyökere is jelentős különbségeket mutathat. Emiatt valahányszor fogászati beavatkozást terveznek, az adott fog gyökerét pontosan fel kell térképezni, amennyiben erre lehetőség van. A korszerű orvosi képalkotó berendezések lehetővé teszik, hogy számos metszeti képet készítsünk a fogról, amelyekből képfeldolgozó eljárással ki lehet nyerni a gyökér tényleges alakját. Erre a problémára a közelmúltban számos megoldás készült, amelyek változatos képalkotó eszközök felvételeiből dolgoztak: • Analoui és társai [8] egy geometriai módszert dolgoztak ki a gyökércsatorna mérésére és modellezésére digitális sztereó radiográf felvételekből. • Hong és társai [40] 2D röntgen és endoszkópos felvételeket használtak a fogak 3D modelljének megalkotására. • Endo és társai [34] ultrahangos felvételekből ismerték fel a gyökércsatornát egy fuzzy logikára épülő eljárással. • Lee és társai [52], valamint Park és társai [68] micro-CT felvételeket és egy 3D rekonstrukciós szoftvert alkalmaztak gyökércsatornák 3D görbületének vizsgálatára, elsősorban őrlőfogakban.

dc_592_12 Foggyökér csatorna meghatározása CT képeken

55

• Willershausen és társai [110] röntgen felvételeket, van Soest és társai [101] pedig optikai átvilágítással készített tomográfiás képeket feldolgozó eljárásokat dolgoztak ki a fogak 3D struktúrájának rekonstruálására.

• Germans és társai [36] interaktív navigációt lehetővé tevő virtuális valóság környezetet hoztak létre a fogak 3D struktúrájának és a gyökércsatornák görbületének megjelenítésére, illetve vizsgálatára.

• A gyökércsatorna morfológiájának vizsgálati módszereit fejlesztette Verma és Love [103], valamint Yamada és társai [118].

• Kaya és társai metszőfogak gyökércsatornájának öregedés során történő geometriai változásait vizsgálták [45].

• További korszerű gyökérmodellezési technikákat javasoltak Li és társai [54], valamint Frisardi és társai [37].

További részletes szakirodalmi áttekintést találunk a [69, 32, 81] munkákban.

Probléma megfogalmazása A fog gyökércsatornájának pontos térbeli elhelyezkedése és alakja fontos információ a fogászati beavatkozás tervezésekor. A CT berendezések fogászati gyakorlatban történő terjedésével lehetőség nyílik ezen információ előállítására. A fent ismertetett ma létező eljárások általános hiányossága azonban az, hogy a képfeldolgozásnak csak egyes lépéseire adnak megoldást. Ahhoz, hogy egy-egy eljárást a gyakorlatban is alkalmazni tudjunk olyan módszereket kell kidolgozni, amelyek megvalósítják a képfeldolgozás minden egymás után illeszkedő lépését, és amelyek egymással kompatibilis eredményt szolgáltatnak. A kutatómunkám célja olyan eljárások kidolgozása volt, amelyek hatékonyan alkalmazhatók micro-CT és CBCT berendezésekkel készített fogászati felvételek feldolgozására, a foggyökér csatorna fogon belüli elhelyezkedésének és alakjának meghatározására. Az eljárásoknak megfelelően pontosnak és automatikusan végrehajthatónak kell lenni annak érdekében, hogy azok a gyakorlatban alkalmazhatók legyenek.

dc_592_12 56

2.2. 2.2.1.

Megoldási módszerek 3D objektumok térbeli vonalvázának meghatározása

Definíció szerint egy 2D tárgy vázát azon pontok összessége alkotja, amelyek a tárgy körvonalának legalább két pontjától azonos, minimális távolságra helyezkednek el. Ezen pontok helyét általában vékonyítással, vagy a legnagyobb körök módszerével szokták megállapítani. Ez a definíció térbeli objektumok esetén olyan vázat eredményez, amely mediális elhelyezkedésű vonalakból és felületekből áll. Mivel a gépi látás legtöbb feladatában a 3D tárgyakat vonalszerű görbék gyűjteményével célszerű reprezentálni, bevezették a térbeli vonalváz fogalmát, amelynek ugyan nem létezik precíz matematikai definíciója, de számos, a tárgyakkal végzett műveletek végrehajtása során hasznos tulajdonsággal rendelkezik [27]. A térbeli vonalvázat meghatározó algoritmusok négy fő kategóriába tartoznak: 1. Vékonyítás során a tárgyak felületéről úgynevezett egyszerű pontokat tüntetünk el, amelyek jelenléte nem befolyásolja a tárgy topológiáját. Ezt legegyszerűbben a „hitor-miss” transzformáció térbeli kiterjesztésével [56] érhetjük el, de léteznek további gyorsított, párhuzamos eljárások is [67]. 2. A távolság-mezőket alkalmazó algoritmusok minden belső pontnak kiszámítják a tárgy felületétől mért minimális távolságát, majd a kiszámított értékekből becsülik a vázat [71]. 3. A geometrikus módszerek általában gráfos reprezentáció segítségével becsülik a felület- vagy vonalvázat [112]. 4. Az általánosított potenciálmezők módszere egy belső potenciálmezőt (pl. elektrosztatikus mezőt) definiál a következőképpen: a tárgy külső pontjaira egységnyi töltéseket helyez el, majd minden diszkrét belső pontban kiszámolja az indukált potenciál értékét [28]. A potenciálmező teljes kioltási pontjai és nyeregpontjai, valamint a tárgy felületén levő nagy görbülettel rendelkező pontok segítségével egy szakaszokból álló, hierarchikus vonalvázat épít fel. Lévén, hogy a térbeli vonalváz a síkbeli váz térbeli kiterjesztése, a 3D váz is olyan pontok gyűjteménye lesz, amelyek legalább két felszíni ponttól azonos, minimális távolságra helyezkednek el. Ez a tulajdonsága teszi alkalmassá a vonalvázat, hogy a vonalváz alapján becslést adjunk a foggyökér középvonalára. A térbeli vonalváz megőrzi a tárgy topológiáját és magába foglalja a tárgyat alkotó elemek hierarchiáját, amely fontos lesz számunkra az elágazási pontok becslésénél.


57

2.1. ábra 3D téglatest vonalváza [84]. Baloldalon a vonalváz 2D-s definíciójának 3D kiterjesztésével kapott vázat látjuk a téglatesten belül, pirossal jelölve. Ebben a vázban találunk 2D alakzatot, a téglatest közepén levő téglalapot. Ennek további erodálásával kapjuk a jobboldali ábrát, amelyen a vonalváz ténylegesen csak vonalakat tartalmaz. A vonalváz ilyen módon történő generálásának pontos leírása [27] irodalomban található.

A térbeli vázkereső módszereknek egy alapos összefoglalója megtalálható a [27] irodalomban.

2.2.2.

Micro-CT felvételek feldolgozása

A fogakról készített micro-CT felvétel egy képsorozat, amelyben minden kép nagy felbontású és egycsatornás intenzitáskép. A sorozat képei valójában párhuzamos metszeteket mutatnak a fogból. Egy sorozatban, a felbontás függvényében, sok száz vagy akár ezer kép is lehet. A szomszédos metszetek közötti távolság és egy metszetben a szomszédos képpontok közötti távolság egyenlő, értéke általában 7-15 µm között van. A micro-CT ezen tulajdonsága számos esetben leegyszerűsíti a képek alapján történő számítási műveletek végrehajtását. A képeken az anatómiai struktúrákat a képpontok intenzitása tükrözi. Normál esetben egy ilyen képen egy világos folt található (dentin), amely világosabb a szélein a fogzománc miatt. A világos folt körülvehet egy vagy több sötét foltot, amelyek a foggyökér és a lágy szövetek helyét jelzik a képen. A fogzománc, ahol látható, valamivel világosabb, mint a dentin. Továbbá nemkívánatos zavarok (úgynevezett műtermékek) is vannak a képeken: körkörös mintázatú textúra, világos területek a dentin körül stb. A kidolgozott képfeldolgozási eljárás fő célja a foggyökér térbeli struktúrájának rekonstrukciója, amelyet a megtalált sötét foltok alapján tudunk megtenni. Ezt követően lokalizáljuk azt a térbeli görbét, amely megfelel a foggyökér középvonalának. A felismert középvonalnak követnie kell a gyökér topológiáját, beleértve a görbületeket és elágazásokat.

dc_592_12 58

2.2. ábra Micro-CT metszetek: metszőfog (bal) és őrlőfog (jobb). A bal oldali ábrán láthatunk egy átlós irányban megnyúlt foggyökér csatornát. Jól látható a fog nagyobb részét kitöltő, a képen világos szürke színű dentin, illetve annak külső szélén a fogzománc. A jobb oldali képen három, viszonylag szabályos alakú gyökércsatorna látható. Mindkét képen láthatunk – elsősorban a fogak körül – a képalkotás tökéletlenségéből adódó úgynevezett műtermékeket.

A 2.3. ábra összefoglalja a képfeldolgozási eljárás fontosabb lépéseit, és megmutatja, hogyan épülnek egymásra az egyes lépések. Az egymás után végrehajtott lépések eredményét sematikusan a 2.4. ábra szemlélteti. MCT1 - Előfeldolgozás Az automatikus képszegmentálást megelőzi néhány képminőség javító lépés, amelyeknek célja a zavaró tényezők eltávolítása vagy elnyomása a képen. A szegmentálás gyorsítása, illetve minőségének javítása érdekében a következő előfeldolgozási lépéseket alkalmaztuk: 1. Egyszerű medián szűrő (12 µm-es képpont méret esetén 7 × 7 pixel méretű), amellyel elnyomjuk a képen jelenlevő magas frekvenciájú zajokat (dentin textúrája) és a körkörös textúrát is. 2. A vizsgált terület (RoI-Region of Interest) kijelölése. Levágjuk a kép képfeldolgozás szempontjából érdektelen széleit, amelyek a külső teret reprezentáló sötét képpontokból állnak. Fontos eltárolni a RoI pontos helyzetét, koordinátáit, ezen adatok alapján lehet az egyes szeletek közötti relatív pozíciókat meghatározni. 3. Egyszerű morfológiai operátorok (pl. dilatáció, erózió) segítségével a képen levő feliratok eltávolítása. Ezen műveletek további előnye, hogy szabályosabbá teszik a foggyökér határvonalát a metszetben. Ezen előfeldolgozási lépések után a kép készen áll a szegmentálásra.

dc_592_12 59

Foggyökér csatorna meghatározása CT képeken

Előfeldolgozás (MCT1) ?

?

FCM osztályozás

Globális küszöb ?

?

Metszet bináris szegmentálása (MCT2-3) ?

Sötét foltok megkeresése régiónövesztéssel (MCT4-5) ?

Korreláció vizsgálat (MCT6)

? -

?

Kézi beavatkozás (MCT7) ?

Gyökér 3D struktúrájának rekonstrukciója (MCT8) ?

Térbeli vonalváz megkeresése (MCT9) ?

Utólagos korrekció a vonalvázon (MCT10) 2.3. ábra A micro-CT felvételeket kiértékelő eljárás lépései és azok logikai sorrendje

2.4. ábra Micro-CT felvételek feldolgozási lépéseinek szemléltetése. Láthatjuk egy metszeti képen a szegmentálást, a gyökércsatorna felismerését, a középpont meghatározását, majd a középpontokból a középvonal meghatározását.

dc_592_12 60 MCT2 - Metszetek szegmentálása A metszet szegmentálásának végső kimenete egy bináris kép lesz. Az elnyomott, de még mindig jelenlévő zajok miatt két particionálást hajtunk végre. Mindkét esetben a kép hisztogramját felhasználó, gyorsított FCM (fuzzy c-means) klaszterező algoritmust [86, 82] alkalmazzuk. Az úgynevezett lokális particionálást c = 4 osztályra való klaszterezéssel hozzuk létre, egy adott metszet képpontjainak intenzitás-hisztogramjából. A globális particionálást egyszerű küszöböléssel hajtjuk végre, ahol a küszöbértéket a teljes 3D adathalmaz képpontjainak intenzitás-hisztogramjából határozzuk meg gyorsított FCM algoritmussal. A valóságban nem mindegyik metszet képpontját számoljuk bele a globális hisztogramba, hanem a metszeteknek csak egy reprezentatív részhalmazát, amely a metszetek kb. 2%-át tartalmazza, egyenletesen elosztva a felvétel hossztengelyén. Ezzel is csökkentjük a művelet végrehajtásának időigényét. A globális küszöbértékkel való küszöbölés egyből bináris képet állít elő. A lokális particionálás eredményeként azonban 4 különböző intenzitás értéket kapunk: ezek az osztályok reprezentatív elemei, jelöljük őket ν1 − ν4 -el. Rendezzük ezeket az értékeket növekvő sorrendbe, azaz νi+1 > νi , ∀i = 1 . . . c − 1. Ezt követően a négy klaszter képpontjait a lokális küszöbértékkel vágjuk, ahol a lokális küszöbértéket a következő összefüggés adja meg: τlocal = (νi+1 − νi )/2, i = arg max{νj+1 − νj , j = 1 . . . c − 1}. j

Az esetek többségében mindkét osztályozás eredményeként kapott küszöbérték jó minőségű szegmentálást végez, de vannak kivételek, amikor valamelyik osztályozás helytelen eredményre vezet. Emiatt szükséges, hogy a két különböző módon előállított bináris kép közül kiválasszuk azt, amelyik várhatóan a helyes partíciót adja. Ezt a választást végezzük el a következő lépésben. MCT3 - Döntési fa Ahhoz, hogy a globális és lokális particionálással előállított képek közül intelligens módon válasszuk ki a megfelelőt, egy döntési fát hoztunk létre, amelyet 250 kritikus metszeti kép segítségével tanítottunk meg döntést hozni. A tanításhoz felhasznált metszetek számának meghatározása kísérletek alapján alakult ki. Ez a szám elegendőnek bizonyult a helyes döntést eredményező döntési fa tanításához. Fontos megjegyezni, hogy ez a szám nagyságrendekkel kisebb a feldolgozott szeletek számánál. A szeletek kiválasztásának menetét a későbbiekben ismertetem.


61

A döntést egy négydimenziós keresési térben kell meghozni a következő négy paraméternek megfelelően: τ1 =

ν3 − ν2 ν4 − ν3 ν2 − ν1 ; τ2 = ; τ3 = ; τglobal . 2 2 2

A döntési fa kimenete megmondja, hogy a globális vagy a lokális partíciót vegyük-e figyelembe az adott metszetnél. A döntési fa tanításánál az egyes tanító lépéseknél a tanító mintáknak minél homogénebb csoportjait igyekszünk kialakítani a minták paraméterei alapján definiált entrópia függvény minimalizálásával [97]. Ezt addig folytatjuk, amíg a fa minden levele homogén nem lesz. A tanítás után az általunk készített döntési fa 7-11 egyszerű összehasonlítás során meg tudja hozni a döntést. A döntési fa által hozott döntés után végül minden metszet esetében egyetlen bináris képet kapunk, amelyen a sötét foltokat kell lokalizálni a képfeldolgozó eljárás következő lépése során. Mint említettem, a döntési fának a létrehozásához 250 darab CT metszeti képet használtunk. A 250 tanító kép kiválasztása a következőképpen történt. Először néhány, különböző CT felvételből származó képsort választottunk ki, amelyek választásakor csak arra ügyeltünk, hogy azok között szerepeljen az összes fogtípus közül legalább egy. A felhasznált fog CT képek számát a képek választása során határoztuk meg, hét képsort kiválasztva értük el, hogy a kiválasztott csoportban az összes fogtípusból legalább egy szerepeljen. Miután kiválasztottuk a hét képsort, elvégeztük a hét képsor minden szeletének szegmentálását a korábban említett módosított FCM algoritmussal [86] mind c = 2, mind c = 4 osztály (globális és lokális osztályozás) esetén. Ezután azt a 250 tanító képet választottuk ki, amelyek esetén a c = 4 osztályra történő osztályozás eredményeként kapott osztályok egyik kombinációban sem korreláltak a c = 2 osztályra végzett osztályozás eredményeként kapott osztályokkal. Ezután az így kiválasztott 250 kép alapján tanított döntési fát leteszteltük 25 db CT felvételből származó képsor, felvételenként kb. 1000 darab metszeti képén végrehajtott osztályozásával kapott osztályokon végrehajtott döntésekkel. A teszt eredményeként kevesebb, mint az esetek 0.1%-ban hozott hibás döntést az így előállított döntési fa.

MCT4 - Sötét foltok megkeresése és azonosítása A bináris kép sötét foltjainak megkeresését és azonosítását egy régiónövesztő algoritmussal végezzük [5]. Mindaddig, amíg sötét képpont van a képen, véletlenszerűen kiválasztunk egy sötét képpontot és növesztünk körülötte egy régiót. A külső teret, amely szintén sötét, figyelmen kívül hagyjuk. A megtalált sötét foltokat egyenként elmentjük. A foggyökér minden olyan ága, amely jelen van az adott metszetben, egy sötét folt formájában jelenik

dc_592_12 62 2.1. táblázat A gyökércsatornát potenciálisan reprezentáló sötét foltok kiválasztásának megvalósított eljárásai Megnevezés

Eljárás leírása

P1 - Csak egy folt

Mindig csak a legnagyobb sötét foltot keresi meg a képen.

P2 - Legfeljebb 2/3/4 folt

Megkeressük a második/harmadik/negyedik legnagyobb sötét foltot is, ha az nagyobb egy bizonyos küszöbértéknél.

P3 - Adaptív

Megkeresünk minden olyan sötét foltot, amelynek mérete megfelel az előre rögzített szabályoknak.

meg a bináris képen. Az esetek nagy többségében egy gyökércsatornához – illetve ha több gyökércsatorna van egy képen, akkor minden egyes gyökércsatornához – csak egy sötét folt tartozik. Néhány esetben azonban a képen megjelenő zajok miatt előfordul, hogy egy gyökércsatorna ág több, egymáshoz közel eső foltként látszik a képen. Ezeket az eseteket a kidolgozott, a szeletek közötti korrelációt vizsgáló módszer – az egymáshoz közel eső szeleteken előforduló foltok számának és helyzetének összehasonlítása alapján – automatikusan jelzi. Ezen esetekben a gyökércsatornát ténylegesen reprezentáló helyes folt kiválasztásához kézi beavatkozásra is szükség van. A gyakorlati tapasztalataink alapján – mint azt az eredmények ismertetésekor részletesen bemutatjuk – az ilyen esetek száma viszonylag alacsony. A sötét foltok középpontjának meghatározására két módszert használhatunk: a) meghatározhatjuk a folt súlypontját, b) meghatározhatjuk a középpontot morfológiai vékonyítással is. Az első módszer viszonylag kis számításidejű, egyszerűen implementálható. Hátránya, hogy nem garantálható, hogy a középpont folton belülre esik. A morfológiai vékonyítás garantáltan az alakzaton belül elhelyezkedő pontot ad eredményül, emiatt a rendszer ezt a módszert alkalmazza annak ellenére, hogy ez egy nagyobb számítási igényű feladat. A sötét foltok automatikus kiválasztására három eljárást dolgoztunk ki és építettünk be a leírt módszereket és algoritmusokat megvalósító rendszerbe. Az eljárásokat a 2.1. táblázat foglalja össze. A rendszerben ki lehet választani a használni kívánt eljárást. A választott eljárás lényegesen befolyásolja a futási időt, az P1, P2, P3 sorrendben növekszik. A leggyorsabb, P1 protokollt csak metszőfogaknál ajánlott használni, mert csak ott nem valószínű az elágazások előfordulása. Minden más esetben, amikor a fogon belül gyökércsatorna elágazásra lehet számítani a P2, illetve P3 eljárások használata vezet helyes eredményre. A futási idő mindegyik eljárás alkalmazása esetében a megengedett foltok számával kvázi lineárisan növekedett.


63

MCT5 - Foltok automatikus feldolgozása Az eredeti micro-CT képeken jelen lévő zaj miatt a felismert sötét foltok bizonyos metszetekben nagyon szabálytalan alakúak lesznek. Ezen esetek nagy részét automatikusan, jól ismert képfeldolgozó módszerekkel le lehet kezelni. Például egy világos „sziget” egy sötét folt belsejében könnyen felismerhető és eltávolítható egyszerű szűréssel. További olyan esetek is vannak, ahol egy morfológiai nyitás/zárás segíteni tud. Vannak viszont olyan esetek is, amikor a gyökércsatorna képe több, egymáshoz közel levő, sötét foltként jelenik meg egy metszeti képen. Ilyenkor kézi beavatkozásra van szükség. Szerencsére az ilyen esetek automatikusan felismerhetők. Az ehhez szükséges vizsgálatot valósítja meg a következő, korreláció vizsgálat lépés. MCT6 - Korreláció vizsgálat A micro-CT felvételek kiértékelése során előfordulhat, hogy kézi beavatkozással javítható a feldolgozás kimenetele. Az esetek túlnyomó többségében (az elvégzett tesztek esetén az esetek több, mint 95%-ában) a kézi beavatkozás szükségessége automatikusan felismerhető az alapján, hogy a szomszédos metszetekben megtalált sötét foltok helye lényegesen eltérő, illetve azok egy adott távolságnál messzebb helyezkednek el egymástól. Emiatt a képfeldolgozó eljárás lépései közé beépítettünk egy olyan lépést, amely megvizsgálja az egymás szomszédságában levő szeleteken a sötét foltok elhelyezkedését. Amennyiben a szomszédos metszetek sötét foltjainak középpontjai egy előre meghatározott küszöbértéknél nagyobb távolságra vannak egymástól, akkor ezt az algoritmus jelzi a feldolgozás során. Ezekben az esetekben valószínűsíthető, hogy kézi beavatkozással javítható a szegmentálás. Egy konkrét példát láthatunk egymás mellett elhelyezkedő metszeteken, egymástól távol elhelyezkedő középpontokra a 2.5. ábrán. A vizsgáló módszer kezeli az elágazásokat is. Amennyiben egy elágazás miatt távolodnak el a foltok középpontjai a szomszédos metszeteken, akkor az elágazástól indulva az egyik irányban szisztematikusan nagyobb lesz a foltok száma. Ezt a korreláció figyelő algoritmus felismeri. MCT7 - Kézi beavatkozás Az automatikus képfeldolgozást megvalósító rendszerben a felhasználónak lehetőséget adunk a feldolgozás egyes lépéseiben kézi beavatkozásra, amelyekkel a rendszer által automatikusan hozott döntéseket bizonyos esetekben felülbírálhatja. A legfontosabb ilyen kézi beavatkozási lehetőség, amikor a felhasználó megváltoztathatja a metszeten belüli szegmentálás kimenetét. A rendszer, a korábban leírt módon, a szegmentálás eredményének elemzése alapján automatikusan jelzi, hogy mely képek esetén valószínű, hogy a szeg-

dc_592_12 64

2.5. ábra Szomszédos metszetek sötét foltjainak korrelációs vizsgálata, két automatikusan felismert jelenséggel: a függőleges tengelyen 70-es számmal jelölt metszetnél elágazás van, míg a 310-es számmal jelölt metszetnél kézi beavatkozással javítható a szegmentálás.

2.2. táblázat Kézi beavatkozási lehetőségek az automatikus feldolgozásba Megnevezés Beavatkozás módja M1

A döntési fa által hozott döntés felülbírálata.

M2

A lokális szegmentálás küszöbértékének megváltoztatása.

M3

Egy vagy több automatikusan megtalált sötét folt eltüntetése.

M4

Több megtalált sötét folt egyesítése egy aktív kontúrmodell (kígyó) segítségével, amelyet a felhasználó kézzel inicializálhat.

mentálás javítható kézi beavatkozással. Ekkor a felhasználó maga dönthet, hogy tényleg kívánatos-e kézi beavatkozás vagy sem. Természetesen kézi beavatkozás nélkül is végrehajthatóak a feldolgozás következő lépései. Ugyanakkor a felhasználó saját kezdeményezésből is beavatkozhat akármelyik metszet esetében. Az összes, rendszerben megvalósított kézi beavatkozási lehetőséget a 2.2. táblázat ismerteti.

MCT8 - Foggyökér objektum rekonstrukciója A metszetekben fellelt sötét foltokat térben egymáshoz illesztve megkapjuk a gyökércsatorna térbeli alakját. Ennek az objektumnak a középvonalát fogjuk térbeli vonalváz formájában megkeresni.


65

MCT9 - Gyökércsatorna vonalvázának meghatározása Voxel-halmazként értelmezett objektumok térbeli vonalvázának becslésére számos algoritmus létezik [27]. Ezek közül egy olyan eljárást kellett választani, amely egy sima, nem túl sok elágazást tartalmazó görbét eredményez. Ilyen eredményt leginkább a potenciál mezőt használó vázkereső algoritmusok szolgáltatnak. Jelen alkalmazásban Cornea és társai módszerét [28] implementáltuk.

MCT10 - Utólagos korrekció A térbeli vonalváz pontosan követi a foggyökér elágazásait és megközelítőleg helyesen kezeli azokat a gyökérszakaszokat is, ahol a gyökér vonulata nem merőleges a metszetekre. Ezeken a helyeken a vonalváz jól becsli a foggyökér középvonalát. Viszont vannak olyan speciális helyek, mint pl. a vonalváz végződései, ahol a váz rövidebb lesz a szükségesnél. Ezeket a végződéseket a sötét foltok középpontjai segítségével pótoljuk. A kezdeti futtatások során nagy görbületű felszíni pontok közelében a vonalváznak parazita ágai jelentek meg. Ezeket az implementált vázkereső algoritmus úgynevezett divergencia küszöbérték paraméterének helyes beállításával tudtuk eltávolítani. A javasolt eljárás lépéseit a 2.3. ábra foglalja össze. Az egyetlen, szaggatott vonallal rajzolt téglalap kézi beavatkozást jelképez, az összes többi lépés automatikus.

2.2.3.

CBCT felvételek feldolgozása

Egy cone beam CT felvétel sorozat általában néhány száz párhuzamos metszetből áll, amelyekben a szomszédos képpontok és szomszédos metszetek közötti távolság általában egyenlő és - a felvételt készítő eszköztől függően - 100-300µm között van. A 2.6. ábrán két CBCT felvételből egy-egy metszeti kép látható. Ilyen alacsony felbontásnál a foggyökér csupán 3-4, néha a vékonyabb szakaszain akár csak egyetlen képpontnyi átmérőjű. Emiatt a gyökércsatorna középvonalának pontos meghatározásához szükség van a képpontnál kisebb struktúráknak a becslésére. A metszeti képeken a képpontokban tárolt intenzitás érték, mint általában a CT felvételek esetén, a Hounsfield egységekben (HU) mért sugárgyengítési együtthatóját adja meg az adott képpont által reprezentált térrészben elhelyezkedő anyagnak.

dc_592_12 66

2.6. ábra Állkapocs metszetek CBCT felvételeken: a képpont mérete 200µm (bal) és 300µm (jobb).

CBCT1 - Előfeldolgozás A felvételek jel/zaj viszonya fordított arányban van a pácienst érő sugárdózissal. Mivel a dózist alacsony szinten igyekeznek tartani, a képek zajszintje gyakran magas. Ilyen körülmények között a zajszint csökkentése érdekében numerikus szűrőkre van szükség. A CBCT felvételek szűrésére egy kontextusfüggő aluláteresztő szűrőt [83] alkalmaztunk, amelynek jellemző tulajdonsága az, hogy a látható éleket jelentősen nem befolyásolja, miközben a „só és bors” jellegű, magas frekvenciájú zajokat hatékonyan elnyomja. CBCT2 - Szegmentálás A régiónövesztő módszerek általában egy előre kiválasztott kiindulási pont körül foglalják el azt a maximális méretű területet a képen, amely teljesíti az előre felállított homogenitási kritériumot. A leggyakrabban használt homogenitási vizsgálat a képpontok intenzitásán alapszik. A vizsgálat során általában a régióhoz tartozó képpontok intenzitásának négyzetes szórását hasonlítják össze egy előre megadott küszöbértékkel. A régiónövesztő módszerek gyakran esnek a zajok és inhomogenitások áldozatául: ilyenkor a megtalált régiók aprók lesznek, és későbbi automatikus egyesítésük is nehezen kivitelezhető. A javasolt módszer nagyon hasonlít a konvencionális régiónövesztéshez, azonban a klasszikus régiónövesztés módszerét úgy módosítottuk, hogy az kezelni tudjon térben lassan változó, illesztve nagyfrekvenciás zajokkal terhelt képeket is. A módszer lényegét a továbbiakban ismertetem: Legyen X az osztályozandó képpontok halmaza: X = {x1 , x2 , . . . , xN },

(2.1)

ahol N a képpontok száma. Az X halmaznak egy fuzzy részhalmazát a következőképpen


67

értelmezzük: F = {(xi , µF (xi ))|i = 1 . . . N },

(2.2)

ahol µF : X → [0, 1] az F fuzzy halmaz karakterisztikus függvénye, µF (xi ), (xi ∈ X), az xi elem F halmazhoz tartozásának tagsági függvénye. Az X halmazban bevezethetünk egy fuzzy relációt, amely valójában az X 2 halmaznak lesz egy fuzzy részhalmaza: Ψ = {((xi , xj ), µΨ (xi , xj ))|i, j = 1 . . . N },

(2.3)

ahol µΨ : X 2 → [0, 1]. Az F fuzzy részhalmaz α-vágása definíció szerint az alábbi halmaz: Xα(F ) = {x ∈ X|µF (x) ≥ α}.

(2.4)

Ψα legyen egy fuzzy reláció, amely fennáll xi és xj között, ha ∃α ∈ (0, 1] : µΨ (xi , xj ) ≥ α.

(2.5)

Ha a Ψα fuzzy reláció érvényes az X 2 = {x1 , x2 , . . . , xN }2 halmazon, akkor azt így is felírhatjuk: xi Ψα xj , ∀xi , xj ∈ X. Az X halmaz xi és xj elemei α-láncot alkotnak, ha létezik ξ1 , ξ2 , . . . ξk ∈ X, amelyre xi Ψα ξ1 Ψα ξ2 Ψα . . . Ψα ξk−1 Ψα ξk Ψα xj .

(2.6)

A CBCT felvételek szegmentálására javasolt módszerünkben két képpontot akkor tekintünk azonos régióhoz (szövettípushoz) tartozónak, ha szomszédsági kapcsolataik által ugyanazon az α-láncon vannak. A szegmentálás elvégzéséhez definiálnunk kell a µF () tagsági függvényt, azt, hogy milyen fuzzy relációt alkalmazzunk, és azt, hogy hogyan definiáljuk az α értékét úgy, hogy hatékonyan meg tudjuk különböztetni a különféle szövettípusokat. A feldolgozás végső célja az, hogy egy fog csontállományát lehetőleg egyetlen összefüggő régióként ismerjük fel, amelynek a belsejében a foggyökér egy másik, sötétebb színű régió lesz. A fenti cél elérése érdekében, a [83] munkánkban megfogalmazottak alapján egy összetett fuzzy tagsági függvényt vezetünk be. Ennek megfelelően két képpontot hasonlónak fogunk tekinteni, ha fizikai távolság és intenzitás szerint egyaránt közel állnak egymáshoz. A javasolt fuzzy relációnak ezek alapján két komponense van: µΨ (xi , xj ) = δ(xi , xj )σ(xi , xj ).

(2.7)

dc_592_12 68 A fenti összefüggésben δ(xi , xj ) a két képpont közötti euklideszi távolság alapján meghatározott hasonlóságot definiálja, míg σ(xi , xj ) a képpontok intenzitásának hasonlóságát jellemzi. Mindkét függvényt úgy kell definiálni, hogy a bennük szereplő elemek nulla különbsége esetén a maximális hasonlóságot jellemző értéket, 1-et adjanak, így lesz µΨ (xi , xj ) a két képpont közötti hasonlóság mutatója. Értelmezzük a két komponenst az alábbiak szerint: δ(xi , xj ) = [1 + κδ dist(xi , xj )]−1/2 ,

"

σ(xi , xj ) =

1 + κσ log

(2.8)

#−1/2

v(xi ) v(xj )

.

(2.9)

A κδ , κσ és α paraméterek segítségével lehetőség nyílik az algoritmus működésének finomhangolására. A κδ magas értéke hasonlóvá teszi az algoritmust a konvencionális régiónövesztéshez, míg az alacsony értékei lehetővé teszik, hogy egymással nem közvetlenül szomszédos képpontokat egy régióba soroljunk. Túl nagy κσ esetén csak teljesen azonos intenzitású képpontokat tudunk azonos régióba sorolni, míg alacsony értékek esetén besorolhatunk egymástól valamelyest eltérő intenzitású képpontokat ugyanabba a régióba. Az összefüggő régiók mérete fordítottan arányos az α választott értékével.

CBCT3 - Térbeli rekonstrukció A szegmentált térbeli tartományokból a sétáló kocka módszer pontosított verziójával [64] előállítjuk a fogak és foggyökerek külső felületét úgynevezett háromszögekből álló háló (mesh) formájában.

CBCT4 - Középvonal meghatározása A CBCT képeken levő foggyökér csatorna középvonalának meghatározásakor is a foggyökér vonulatát (középvonalát) egy térbeli vonalvázként becsüljük. Bár az eredményként kapott görbék minősége indokolná a potenciálmezős algoritmus használatát, az alacsony felbontás miatt ezt nem tehetjük meg a CBCT felvételek esetén. Az objektumok mesh alapú felületeiből kiindulva, az úgynevezett felület kontrakció módszerét [9] alkalmazva határozzuk meg a középvonalat.


69

2.7. ábra Metszeti képek szegmentálásának részeredményei. Mindegyik sorban egy-egy metszeti kép feldolgozásának részletei láthatóak. Az első oszlopban (a,d,g) az eredeti képeket, a második oszlopban (b,e,h) a 4 osztályra klaszterezett képeket, míg a harmadik oszlopban (c,f,i) a bináris szegmentált képeket láthatjuk. A felismert foltok középpontját kis kereszt jelöli az utolsó oszlop képein.

2.3. 2.3.1.

Eredmények Micro-CT felvételek feldolgozása

A micro-CT képfeldolgozó eljárást 25 darab CT képen teszteltük, amelyek között 17 darab metszőfog és 8 darab őrlőfog volt. Minden fog külön CT képen szerepelt. A javasolt eljárás teljesen automatikusan dolgozta fel a rendelkezésünkre álló felvétel sorozatok 92%-át, míg a többi felvétel esetén kézi beavatkozásra volt szükség. Egy i5-ös processzort tartalmazó asztali számítógépen egy metszet szegmentálása 0,2-0,3 másodpercig, míg a gyökércsatorna rekonstrukciója 1 másodpercig tartott. A futtatások további részleteit a 2.4. táblázat tartalmazza. A 2.7. ábrán a metszeti képek szegmentálásának részeredményeit és végső kimenetét láthatjuk. Az ábrán három metszeti kép feldolgozásának részletei szerepelnek, ezek különböző mértékű kihívást jelentettek az automatikus feldolgozó eljárásnak. Az első sorban egy egyszerű esetet látunk, két könnyedén felismerhető sötét folttal, amelyek a gyökércsatorna két különböző ágának keresztmetszetei (tehát volt egy elágazás valahol egy távoli metszetben). A második sorban egy furcsa formájú sötét folt van a képen, amelyet az algoritmus

dc_592_12 70 sikeresen felismert gyökércsatornaként. A harmadik sorban szereplő metszet egy komplikáltabb eset: három sötét foltot talált az algoritmus, de ezek valójában a csatornának csak két ágához tartoznak a valóságban. Ebben az esetben hasznos a szeletek közötti korrelációs vizsgálat, amelyből kiderül, hogy a szomszédos metszetekben is ugyanilyen foltok találhatóak, vagy a két folt egyesítésére kézi beavatkozásra van szükség. A 2.8. ábra néhány jellemző szegmentálási eredményt mutat be. Az ábrán két foggyökér négy-négy különböző nézete látható, az ábrákon fekete vonal jelzi az automatikusan felismert gyökércsatorna középvonalát. Mindkét foggyökér struktúráját több, mint 900 metszetből rekonstruáltuk. A képen szereplő fogak metszeteinek több, mint 99%-át globális küszöbérték alkalmazásával szegmentáltuk.

2.8. ábra Két metszőfog gyökerének különféle nézetei a felismert középvonallal.


71

2.9. ábra Egy rekonstruált metszőfog különböző nézetekből: (a) a metszőfog külső felülete; (b) a fog belsejében elhelyezkedő gyökércsatorna külön megjelenítve; (c)-(e) különböző irányokból nézve a fog belsejében levő gyökércsatorna elhelyezkedése a fog felületének virtuális „felmetszése” után.

2.3.2.

Cone beam CT felvételek feldolgozása

A CBCT felvételeken több fog is szerepelt, így a validáció során a gyökércsatornát felismerő eljárás elindítása előtt megadtuk melyik fog gyökércsatornájának felismerését szeretnénk végrehajtani a szegmentálandó fog (metsző- vagy őrlőfog) egyetlen belső pontjának kijelölésével. Ezt követően a feldolgozás teljesen automatikusan történt. Az algoritmus megkereste a legnagyobb intenzitású képpontokat a megjelölt kezdőpont közelében, és ezeket használta a fuzzy láncok korábban ismertetett módszerénél kiinduló pontokként. Egy fog térfogata kb. 0, 5 − 1, 5 × 106 képpontot foglal magában. Ennyi adat kiértékelése egy i5-ös processzort tartalmazó asztali számítógépen kb. 1 másodpercet vett átlagosan igénybe. A futtatások további részleteit a 2.4. táblázat tartalmazza. A szegmentált fog és foggyökér felületét egy mesh hálózat (háromszög gyűjtemény) formájában kapjuk meg. Ez a struktúra bármilyen irányból és bármilyen metszetű nézetben viszonylag egyszerűen megjeleníthető. A 2.9. ábrán egy metszőfog, a 2.10. ábrán pedig egy őrlőfog különféle nézetei láthatóak. Mindkét ábrán a következő részletek szerepelnek: (a) a fog külső felülete, (b) a foggyökér csatorna a foggal azonos felbontásban, (c)-(e) metszeti képek, amelyeken különböző irányokból nézve a fog belsejében levő gyökércsatorna elhelyezkedését láthatjuk a fog felületének virtuális „felmetszése” után. A 2.11. ábrán egy metszőfog struktúráját, a gyökerét és annak középvonalát láthatjuk. Mindhárom nézet ugyanazt a fogat ábrázolja. A koordinátatengelyeken szereplő számok képpontokban vannak megadva. A felbontás ismeretében ezeket egyszerűen át lehet számolni milliméterbe. Zölddel ábrázoltuk a fog metszeteit, vörössel a foggyökeret, míg a fekete vonal az eljárás által meghatározott középvonalat jelöli. A 2.12. ábrán egy őrlőfog két különböző nézetét láthatjuk, ugyanazokat a jelöléseket használva, mint az előző ábránál.

dc_592_12 72

2.10. ábra Egy rekonstruált őrlőfog különböző nézetekből: (a) az őrlőfog külső felülete; (b) a fog belsejében elhelyezkedő gyökércsatorna külön megjelenítve; (c)-(e) különböző irányokból nézve a fog belsejében levő gyökércsatorna elhelyezkedése a fog felületének virtuális „felmetszése” után.

Megállapíthatjuk, hogy a szegmentálás általában pontosabban sikerül a metszőfogak esetén, mert a őrlőfog gyökerének számos elágazása van, és ezek az ágak gyakran sokkal szűkebbek egy metszőfog gyökerénél. A felbontás javításával az őrlőfogak szegmentálásának minőségén változtatni lehetne. A 2.11. ábrán egy sikeresen szegmentált metszőfogat látunk, míg a 2.12. ábrán látható őrlőfognak valójában hosszabbak a gyökerei az ábrán láthatónál. A hiányzó részeket a felvétel viszonylag gyenge felbontása miatt nem volt lehetséges jobb minőségben szegmentálni a teljesen automatikus eljárással. A 2.13. ábra ugyanezen fog rekonstrukciójának eredményét mutatja, azonban itt az automatikus eljárást kézi beavatkozással javítottuk. Az automatikus eljárás megtalálta a gyökércsatorna több, a 2.12. ábrán nem látható szakaszát, azonban a gyökércsatorna - a CT felvétel felbontásánál kisebb mértékű - szűkülete miatt azokat nem kötötte össze. A manuális beavatkozásnál ezeket összekötve a rekonstrukciós algoritmus láthatóan pontosabb eredményt adott. A bemutatott jellemző példák alapján láthatjuk, hogy a gyökércsatorna középvonalának meghatározását a kidolgozott és megvalósított képfeldolgozó eljárás hatékonyan és pontosan végezte el. Mivel a gyökércsatorna átmérője ritkán több 10 képpontnál, a vonalváz kereső algoritmusnak 5-6 iterációt kell végrehajtani. A kapott vonalváz sima, a gyökér vonulatának közepén helyezkedik el, és a gyökér elágazásainál hűen követi a csatorna tényleges alakját. A szakértői értékelés részletei a következő alfejezetben kerülnek ismertetésre.

2.4.

Az elért eredmények értékelése

Mind a micro-CT, mind a CBCT felvételek feldolgozását végző eljárás képes volt a képek több, mint 90%-át automatikusan és pontosan feldolgozni. A megmaradt képek ese-


73

2.11. ábra Egy rekonstruált metszőfog különféle nézetei gyökérrel és középvonallal.

2.12. ábra Egy rekonstruált őrlőfog térbeli ábrázolása háromágú foggyökérrel és középvonallal.

dc_592_12 74

2.13. ábra Egy őrlőfog gyökércsatornájának felismerése (a) kizárólag automatikus eljárással és (b) kézi beavatkozással kiegészített eljárással.

tén a szegmentálás eredményét lényegesen lehetett javítani a manuális beavatkozással. A micro-CT felvételeket feldolgozó algoritmus 25 darab felvétel feldolgozásával került tesztelésre. Két őrlőfog esetén több, egymással nem korreláló foltokat tartalmazó szomszédos metszeti kép fordult elő, ahol a korrelációt figyelő algoritmus nem tudta eldönteni, hogy elágazásról van-e szó, vagy a felvételeken szereplő zaj, illetve a gyökércsatornában elhelyezkedő csontkinövés okozta manuálisan javítandó hibás szegmentálás történt. A többi 23 fog esetén összesen 227 esetben (az összesen 23.391 feldolgozott metszeti képből) jelezte az algoritmus, hogy potenciálisan szükség lehet kézi beavatkozásra. Ezek közül 119 esetben a szegmentálást ellenőrző szakértő is úgy ítélte meg, hogy szükség van javítani, tehát helyes döntést hozott a szegmentált foltok elhelyezkedését figyelő módszer. Az automatikus kiértékelést végző eljárás szakértői validációjának eredményét a 2.3. táblázat foglalja össze. Látható, hogy a szeletek 1%-ánál kevesebb esetben jöhet csak szóba, hogy a kézi beavatkozással érdemben lehet javítani az automatikus feldolgozás eredményét, és ezen esetek kb. felénél az eljárás erre helyesen fel is hívja a felhasználó figyelmét. A CBCT felvételeket feldolgozó eljárást 36 darab fog gyökércsatornájának felismerésével validáltuk. Ezek közül 29 fog esetében a gyökércsatorna felismerését, illetve annak középvonalának meghatározását a szakértő véleménye szerint sem lehetett tovább javítani. A megmaradt 7 fog közül 4 esetében egyszerű, az implementált eljárást megvalósító környezetben egyetlen művelettel megvalósítható kézi beavatkozással, illetve az automatikus eljárás megváltoztatott adatokon történő újrafuttatásával tökéletes gyökércsatorna meghatározás volt végrehajtható. A kiértékelés részleteit a 2.4. táblázat tartalmazza.

dc_592_12 75

Foggyökér csatorna meghatározása CT képeken

2.3. táblázat Az automatikus végrehajtás validációjának részletei micro-CT képek feldolgozása során Szakértők értékelése Nem szükséges

Kézi beavatkozás

beavatkozás

szükséges

23150

14

108

119

Nem szükséges Automatikus beavatkozás feldolgozás

Kézi beavatkozás szükséges

2.4. táblázat A két automatikus képfeldolgozó eljárás futtatásainak eredményei Paraméter

Micro-CT

CBCT

A feldolgozott teszt képek száma

25

36

Kizárólag automatikus eljárással sikeresen feldolgozott képek

23

33

Sikeres feldolgozás aránya

92.0%

91.7%

Sikeres feldolgozás aránya metszőfogak esetén

100.0%

100.0%

Sikeres feldolgozás aránya őrlőfogak esetén

75.0%

83.3%

Legrövidebb feldolgozási idő (másodperc)

288

0.69

Átlagos feldolgozási idő (másodperc)

341

0.93

Leghosszabb feldolgozási idő (másodperc)

490

1.61

dc_592_12 76 A validáció eredménye alapján megállapíthatjuk, hogy a CBCT felvételeket feldolgozó eljárás – hasonlóan a micro-CT felvételeket feldolgozó eljáráshoz – pontosan és hatékonyan képes a fogak gyökércsatornájának automatikus meghatározására. Mindkét modalitást feldolgozó eljárás eredményéül szolgáltatott adatok alapján elkészíthető a fogak gyökércsatorna középvonalának matematikai leírása a középvonalat közelítő spline görbékkel, ahogy azt a [31, 30] irodalmak ismertetik. Mindkét kidolgozott képfeldolgozó eljárás legnagyobb előnye az, hogy azok automatikusan végrehajthatók, és lefedik a CT felvételek feldolgozásának minden lépését. Ez lehetővé teszi az eljárások közvetlen gyakorlati alkalmazását. Az eljárások megvalósításához szükséges egyes módszerek kidolgozásánál nemcsak az egyes módszerek hatékonyságára, hanem arra is ügyelni kellett, hogy azok egymással kompatibilis módon működjenek, egymással kompatibilis eredményt szolgáltassanak. Micro-CT felvételek feldolgozása esetén az automatizálás egyik legfontosabb eleme az eljárás azon lépése, amely automatikusan felismeri a kézi beavatkozás szükségességét. A korreláció figyelő algoritmus figyelembe veszi, hogy nemcsak szegmentálási hiba esetén lehet alacsony két szomszédos cellán szereplő foltok elhelyezkedésének korreláltsága, hanem a gyökércsatorna elágazási helyének közelében is, így a döntése meghozatalakor mind a szegmentált gyökércsatorna középpontok elhelyezkedését, mind azok számát figyelembe veszi a vizsgált szelet környezetében elhelyezkedő metszeti képeken. A CBCT felvételeket feldolgozó eljárás CBCT2 lépése hivatott a fog dentinje által kitöltött térrészt, régiót meghatározni, amelynek kimenete alapján a CBCT3 lépésben állítjuk elő a fog térbeli modelljét. Ezt használjuk fel majd a gyökércsatorna meghatározásához. A CBCT2 lépés sajnos nem minden esetben hoz létre egymással összefüggő dentin, illetve gyökércsatorna régiókat, így a gyökércsatorna keskeny szakaszán néhol nem lehetséges a gyökércsatorna teljes hosszán az automatikus feldolgozás. Mint azt már a 2.13. ábrán bemutatott példába láthattuk, ilyen esetekben a régiók kézi egyesítése segít. Az automatikus működés megvalósításának fontos eleme volt micro-CT feldolgozás során az MCT3 lépés, mikor döntési fát hoztunk létre, amely a két különböző módon végrehajtott szegmentálás eredménye közül kiválasztja a helyes gyökércsatorna meghatározást lehetővé tevőt. Ennek a döntési fának a létrehozásához 7 különböző CT felvételből nyert képből származó összesen 250 darab CT metszeti képet, illetve azok szegmentálási eredményét használtuk. A 250 tanító kép kiválasztása is automatikusan történt. Miután kiválasztottuk a 7 képsort – amelyek választásakor ügyeltünk arra, hogy azok között szerepeljen az összes fogtípus közül legalább egy –, elvégeztük a 7 képsor minden szeletének szegmentálását a korábban említett módosított FCM algoritmussal [86] mind c = 2, mind c = 4 osztály


77

(globális és lokális osztályozás) esetén. Ezután azt a 250 tanító képet választottuk ki, amelyeknél a c = 4 osztályra történő osztályozás eredményeként kapott osztályok egyik kombinációban sem korreláltak a c = 2 osztályra végzett osztályozás eredményeként kapott osztályokkal. Ezután az így kiválasztott 250 kép alapján tanított döntési fát leteszteltük 25 CT kép, képenként kb. 1000 darab metszeti képén végrehajtott osztályozásával kapott osztályokon végrehajtott döntésekkel. A teszt eredményeként kevesebb, mint az esetek 0.1%-nál hozott hibás döntést az így előállított döntési fa.

2.5.

Összefoglalás

Ebben a fejezetben javasoltam két, összetett automatikus eljárást a foggyökér csatornájának felismerésére és középvonalának meghatározására micro-CT, illetve CBCT felvételek alapján. Az eljárások az automatikus végrehajtás mellett megadják a lehetőséget a kézi beavatkozásra, amennyiben a felhasználó így szeretne javítani a gyökércsatorna felismerésének eredményén. A két eljárás kidolgozása során képfeldolgozó módszereket dolgoztam ki, illetve módosítottam annak érdekében, hogy azokat lehetséges legyen automatikusan végrehajtani a két eljárás során a foggyökér csatorna pontos meghatározása céljából. Mindkét kidolgozott eljárást implementáltam, majd validáltam valós micro-CT és CBCT felvételeken. A felvételek több, mint 90%-a esetén az eljárások pontosan meghatározták a gyökércsatorna középvonalát automatikusan, kézi beavatkozás nélkül. Az eljárások eredménye alapján elkészíthető a fogak gyökércsatorna középvonalának matematikai leírása. A megvalósított eljárások a jövőben hozzájárulhatnak a foggyökér csatornát érintő fogászati beavatkozások hatékonyabb tervezéséhez, ezáltal azok sikerességének növeléséhez. Továbbfejlesztési lehetőségek A bemutatott kutatás kézenfekvő folytatása lehet annak vizsgálata, hogy a kidolgozott algoritmusok milyen módon építhetők be képalkotó berendezések képrekonstrukciós környezetébe. Ezáltal olyan szolgáltatással bővülhet a képalkotó berendezések funkcionalitása, amely igen nagy segítség a fogászati diagnosztikai munkában.

dc_592_12 78

2.6.


2. tézis: Olyan képfeldolgozó eljárásokat dolgoztam ki, amelyek fogakról készített CT felvételek alapján képesek automatikusan a foggyökér csatornát pontosan azonosítani és a foggyökér csatorna középvonalának térbeli alakját meghatározni. Az egyik eljárás micro-CT felvételek, a másik eljárás cone beam CT felvételek feldolgozása során alkalmazható. Az eljárásokat validáltam valós micro-CT és cone beam CT felvételek felhasználásával. A validáció során igazoltam, hogy mindkét eljárás az esetek több, mint 90%-ban automatikusan képes volt a gyökércsatorna kézi feldolgozással azonos pontosságú felismerésére, annál lényegesen rövidebb idő alatt.

2.6.1.


A fejezetben ismertetett eredményeket bemutató legfontosabb publikációk: [16], [19], [17], [87], [18], [13].

dc_592_12

3. fejezet Térben lassan változó intenzitású zajjal terhelt MR koponyafelvételeken agyi régiók elkülönítésére szolgáló algoritmus kidolgozása A mágneses rezonanciás (MR) képalkotás az orvosi gyakorlatban gyakran alkalmazott vizsgálati módszer elsősorban az eredményül kapott képek jó felbontása és a különböző szöveti régiók kontrasztos megjelenítése miatt [108]. Az MR képek automatikus feldolgozása, szegmentálása azonban közel sem triviális feladat. A képeken megjelenő, különböző forrásból származó zaj gyakran komoly nehézséget okoz a feldolgozás során. Az egyik leggyakoribb zajtípus az úgynevezett inhomogén intenzitású zaj [107], amely a tér egy adott irányába lassan változó, alacsony frekvenciájú, a hasznos jelhez viszonyítva viszonylag nagy amplitudójú additív zaj. E jelenség következményeként előfordul, hogy azonos anyagú lágyszöveti területek egymástól lényegesen különböző intenzitású pontokként jelennek meg egy MR képen. Míg a magas frekvenciájú zajok eltávolítására számos gyors és viszonylag pontos eljárás létezik [25], addig az inhomogén intenzitású zaj kompenzálása lényegesen nehezebb és számításigényesebb feladat [107], így annak megvalósítása komoly kihívást jelent az MR képek szegmentálása során. A fejezetben az MR képeken megjelenő inhomogén intenzitású zaj kompenzálására, illetve a képek szegmentálására alkalmas, a korábbi megoldásoknál lényegesen alacsonyabb számításigényű algoritmust ismertetek. 79

dc_592_12 80

3.1.


Forrásuk alapján az MR képalkotás során megjelenő inhomogén intenzitású zajnak két típusát különböztethetjük meg: a képalkotó berendezésből eredő, illetve a vizsgált térrészben elhelyezkedő anyagból eredő zajt. A képalkotó berendezésből eredő inhomogén zaj kompenzálására léteznek hatékony kalibrációs eljárások, amelyek egy homogén anyagú fantom segítségével készített képből nyerik ki az eszköz helyes beállításához szükséges a priori információt [10]. Ezzel szemben, a vizsgált beteg testének alakjából, pozíciójából, illetve a szövetek struktúrájából vagy elhelyezkedéséből [77] eredő inhomogén zajként jelentkező hatásokat sokkal nehezebb kezelni [107]. Ilyen esetben a képek utólagos feldolgozására van szükség. A szakirodalomban számos eljárást találunk az inhomogén intenzitású zajjal terhelt MR képeken megjelenő zavarok kompenzálására. Ilyenek a • homomorf szűrők [24, 43], • polinom vagy B-spline felület illesztése intenzitás [102] vagy gradiens [105] segítségével, • maximum likelihood becslés [72], • Markov valószínűségi mezők [122], • fuzzy c-means (FCM) osztályozás [22, 70, 78], • nem parametrikus becslés [29], • magas frekvenciák maximalizálása hisztogram alapján [79], • információ maximalizálás hisztogram alapján [106], • hisztogram illesztés módszere [80]. A leggyakrabban alkalmazott inhomogén zajt becslő módszerek az FCM (fuzzy cmeans) algoritmusból [21] származtatott eljárások, az FCM algoritmus könnyű áttekinthetősége és egyszerű implementációja miatt. Számos adaptáció létezik: Pham és Prince bevezetett egy módosított költségfüggvényt, amelyben az additív zajbecslés mellett megtalálható egy, a becsült additív zajt lassú térbeli változásra kényszerítő tag is [70]. A javasolt módosításhoz a szerzők mellékeltek egy multirezolúciós optimalizálási eljárást is, amelynek segítségével felgyorsították az amúgy rendkívül lassú becslő algoritmust. Ahmed bevezetett egy átlagoló operátort, amelynek segítségével egy képpont adott osztályba sorolását a szegmentálás során a képpont közvetlen szomszédainak besorolása is befolyásolja [6]. Ez

dc_592_12 Agyi régiók elkülönítése MR felvételeken

81

a megoldás kissé lecsökkentette az algoritmus számítási komplexitását, de a zaj megállapítására használt kritériumok gyakran vezettek hibás osztályozáshoz. Siyal és Yu [78] bevezettek egy új szűrési eljárást a becsült zaj amplitudo változásának simítására az FCM algoritmus minden iterációjában, de az osztályozási algoritmusuk nem determinisztikus a simító szűrő természete miatt. A felsorolt szegmentáló algoritmusok egy része képes viszonylag pontos szegmentálásra, azonban ezek hosszú futási ideje megnehezíti azok gyakorlati alkalmazását. A gyors futási idejű algoritmusok ezzel szemben nem teszik lehetővé a gyakorlati alkalmazások által megkövetelt minőségű szegmentálást [107]. A disszertációban ismertetett kutatást ezen kettős követelmény kielégítése vezérelte, a cél egy gyors és pontos szegmentálást megvalósító algoritmus kidolgozása volt. A korábbi évek során több FCM algoritmusra alapuló megoldást is javasoltunk inhomogén zajjal terhelt képeken történő zajelnyomásra, és ezáltal a képeken végrehajtott szegmentálás eredményének javítására. Több olyan eredeti szűrési eljárás hatékonyságát vizsgáltuk meg, amelyekben figyelembe vettük a becsült zajok természetét [86]. Kidolgoztunk egy, a korábbi megoldásoktól eltérő elemeket tartalmazó FCM algoritmust, amelyben a partíciókat valószínűségi fuzzy tagsági függvények [21] és lehetőségfüggvények [49] lineáris kombinációjával modelleztük [89], és kimutattuk ennek a megoldásnak mind a szegmentálás pontosságában, mind pedig a futási idő csökkenésében mérhető kedvező hatását. Végül kidolgoztunk egy gyorsított módszert, amelyben a költségfüggvény átfogalmazásával a pillanatnyi korrigált kép hisztogramját felhasználva, lényegesen csökkenteni tudjuk a szegmentáláshoz, illetve osztályozáshoz szükséges számítási műveletek számát [85]. A módszer alkalmazásával az algoritmus futási ideje az eredetinek 3-5%-ára csökken anélkül, hogy a szegmentálás pontossága romlana. A fejezet ezt a gyors algoritmust ismerteti, amely egycsatornás intenzitásképek c-means típusú klaszterezésével történő szegmentálására szolgál. A gyorsítás alapötlete a következő. A fuzzy klaszterezés során alkalmazott négyzetes költségfüggvény (hibafüggvény) iteratív minimalizálása során megkeressük azokat a számítási műveleteket, amelyek csak a képpontok intenzitás értékétől – és nem azok helyétől – függenek, és a költségfüggvény minimalizálását ezekre az intenzitás értékekre végezzük el. Mivel egy képen megjelenő intenzitás értékek száma jóval kisebb a képpontok számánál, az eljárás számításigénye ezzel jelentősen csökkenthető. A következő fejezetben bemutatom a fent leírt elven alapuló módszert, illetve a módszer alkalmazását a fuzzy c-means klaszterező algoritmus gyorsítására. Megjegyezzük, hogy a módszer elvileg alkalmazható bármilyen c-means típusú osztályozó eljárás gyorsítására.

dc_592_12 82

3.2. 3.2.1.

Alkalmazott módszerek Az FCM algoritmus

A hagyományos FCM algoritmus a bemeneti vektorhalmaz elemeit optimálisan osztályozza egy előre megválasztott számú, c darab osztályba egy kvadratikus költségfüggvény minimalizálásával. Amikor egy egycsatornás, intenzitás értékeket tartalmazó kép képpontjait osztályozzuk, az FCM algoritmus a kép képpontjainak xk , k = 1 . . . n intenzitás értékeit fogja particionálni. Az FCM algoritmus költségfüggvénye (minimalizálandó hibafüggvénye) JFCM =

c X n X

2 um ik (xk − vi ) ,

(3.1)

i=1 k=1

amelyet a következő, úgynevezett valószínűségi megszorítás c X

uik = 1 ∀k = 1 . . . n ,

(3.2)

i=1

figyelembe vételével optimalizálunk, ahol: • uik ∈ [0, 1] jelöli a fuzzy tagsági függvények értékeit, amelyek azt írják le, hogy a k-dik képpont milyen mértékig tartozik az i-dik osztályba, • vi az i-dik osztály prototípusa, reprezentatív eleme vagy centroidja, • m > 1 az úgynevezett Bezdek-féle fuzzy kitevő [21], amellyel az algoritmus fuzzy jellegét lehet szabályozni. A költségfüggvény optimalizálása a gradiensek zérus átmeneteinek és Lagrange szorzók alkalmazásával történik. A minimum megkereséséhez felváltva optimalizáljuk a JFCM költségfüggvényt {uik } szerint, i = 1 . . . c, k = 1 . . . n, rögzített vi értékek mellett, és {vi } szerint, i = 1 . . . c, rögzített uik partíciós mátrix mellett [21]. Minden iterációban az optimális fuzzy partíciót és a klaszterek prototípusát az alábbi összefüggésekkel számítjuk ki: (xk − vi )−2/(m−1) uik = P ∀ i = 1 . . . c, ∀ k = 1 . . . n , (3.3) c −2/(m−1) (xk − vj ) j=1

és

n P

vi =

um ik xk

k=1 n P

k=1

∀i = 1...c .

(3.4)

um ik

A vi osztályprototípusok megfelelő inicializálása után a (3.3) és a (3.4) összefüggéseket felváltva alkalmazzuk mindaddig, amíg ki nem elégítjük a konvergencia kritériumot, vagyis


83

amíg az osztályok prototípusainak változása egy ciklus alatt nem lesz kisebb, mint egy előre rögzített ε küszöbérték. Ekkor a k-dik képpontot (k = 1 . . . n) besoroljuk abba az osztályba, amelybe az uik tagsági függvény alapján a leginkább beletartozik, vagyis abba a wk -dik osztályba, amelyre igaz a következő összefüggés wk = arg max{uik , i = 1 . . . c} . i

3.2.2.

(3.5)

Inhomogén zaj kompenzálására szolgáló modellek

A hagyományos FCM algoritmus fent ismertetett formájában a k-dik képpontot annak xk tényleges intenzitás értéke szerint osztályozza, amelyet csak ideális, zajmentes körülmények között mérhettünk volna meg. Viszont a valóságban mért yk adat eltér a tényleges xk adattól. Jelen fejezetben azt feltételezzük, hogy a képeinket lassan változó intenzitású zaj terheli, eltekintünk minden egyéb (pl. magas frekvenciájú) zajtól. A szakirodalomban három különböző modell [70] segítségével írják le az ilyen inhomogén intenzitású zaj hatását. Amennyiben az inhomogén intenzitású zajt additív zajként fogjuk fel, a k-dik képpont (k = 1 . . . n) intenzitását a következőképpen modellezhetjük: yk = xk + bk ,

(3.6)

ahol bk jelöli az additív zajt az adott képpontban [6, 70, 78]. Amennyiben erősítési mezőként modellezzük az inhomogén zajt [86], akkor bármely k-dik képponthoz hozzárendelünk egy gk erősítési értéket, azaz y k = gk x k .

(3.7)

Végül lehetséges alkalmazni a következő logaritmikus összefüggést, amely alapján az erősítési mező additív zajként jelenik meg [55]: log yk = log gk + log xk .

(3.8)

Bármelyik zajmodellt is alkalmazzuk, a becsült zaj inhomogén tulajdonságának, vagyis változó amplitudojának térben lassan kell változnia. A korábban említett, gradiensek zérus átmeneteiből becsült inhomogenitás értékek ezt a követelményt nem teljesítik. Ezért minden iteráció végén a becsült additív zajra egy erős simító (átlagoló) szűrőt kell alkalmazni.

dc_592_12 84

3.2.3.

Az FCM algoritmus alkalmazása inhomogén zaj becslésére

A fenti additív zajmodellt bevezetve az FCM algoritmus költségfüggvényébe a következő kvadratikus összefüggést kapjuk: JFCM−b =

c X n X

2 um ik (yk − bk − vi ) .

(3.9)

i=1 k=1

Ennek optimalizálási formulái a következők: • Bármely k-dik képpont (k = 1 . . . n) és i-dik osztály (i = 1 . . . c) esetében a megfeleltetési fuzzy tagsági függvény értéke: (yk − bk − vi )−2/(m−1) uik = P . c (yk − bk − vj )−2/(m−1)

(3.10)

j=1

• Bármely i-dik osztály (i = 1 . . . c) esetében, a prototípus becsült értéke: n P

vi =

k=1

um ik (yk − bk ) .

n P k=1

(3.11)

um ik

• Bármely k-dik képpont (k = 1 . . . n) esetén a becsült additív zaj értéke: c P

um ik vi

bk = yk − i=1 c P

i=1

.

(3.12)

um ik

A hagyományos FCM osztályozó algoritmust a 3.1. táblázat foglalja össze [78].

3.2.4.

Gyors osztályozás hisztogram alapján

Amikor egy klaszterező algoritmussal gyors osztályozást akarunk elérni, a hasonló bemeneti adatok együttes kezelése (aggregálása) egy kézenfekvő megoldásnak tűnik. A hagyományos FCM algoritmus – a képek szegmentálási módszertanában – az úgynevezett globális információ alapján szegmentálást végző módszerek közé tartozik. Ez azt jelenti, hogy az azonos színű (intenzitású) képpontok azonos fuzzy tagsági függvény értékek szerint lesznek osztályokba sorolva, függetlenül a képen elfoglalt helyüktől. Ez alapján könnyen belátható, hogy egy egycsatornás intenzitás értékeket tartalmazó kép szegmentálását végrehajthatjuk a kép hisztogramja alapján, gyakorlatilag a képen előforduló intenzitás értékeket osztályozva a képpontok helyett [82]. Ezzel – mivel a képen szereplő intenzitás értékek és a

dc_592_12 85

Agyi régiók elkülönítése MR felvételeken 3.1. táblázat A hagyományos FCM osztályozó algoritmus 01 t = 0 . (t=0)

02 Inicializáljuk az additív zajt bk

= 0, ∀k = 1 . . . n . (t=0)

03 Adjunk az osztályprototípusoknak kezdeti értéket: vi

(t=0)

6= vj

, ∀i, j = 1 . . . c; i 6= j .

04 Ismételjük a következő lépéseket 05

t←t+1 ;

06

Számoljunk új uik tagsági függvény értékeket, ∀i = 1 . . . c és ∀k = 1 . . . n esetén, a (3.10) összefüggés szerint;

07

Számoljunk új vi osztályprototípusokat, ∀i = 1 . . . c, a (3.11) összefüggés szerint;

08

Számoljunk új bk becsült additív zajt , ∀k = 1 . . . n, a (3.12) összefüggés szerint;

09

Alkalmazzunk egy átlagoló szűrőt a becsült zaj simítására;

(t)

(t)

(t)

10 amíg a konvergencia bekövetkezik, azaz

c P i=1

(t)

(t−1)

|vi − vi

|<ε.

11 Soroljuk be a k-dik képpontot (k = 1 . . . n) abba az osztályba, amelynek indexe arg max{uik , i = 1 . . . c}. i

képpontok száma között általában nagyságrendi különbség van – az osztályozás lényeges gyorsítását érhetjük el. A fenti állítás sajnos csak addig igaz, amíg a képen nincs változó intenzitású zaj, és nem akarjuk a zajt becsülni a kép alapján. A gondot az okozza, hogy a zaj inhomogén, amplitudoja helyfüggő, vagyis minden képpont becsült zaja és besorolása függ annak helyétől, illetve a szomszédos képpontok besorolásától. A 3.1. táblázat szerint a változó intenzitású zaj modelljét tartalmazó FCM algoritmus (nevezzük ezt FCM-b algoritmusnak) minden iterációja során egyenként számoljuk a képpontok osztályokba való besorolását, az osztályok prototípusait és a becsült zajt egyaránt. A bemeneti kép azonos intenzitású képpontjait nem lehet azonos módon kezelni, mert a bennük levő zaj különbözővé teszi őket. Emiatt tehát nem lehetséges a fent leírt módszert közvetlenül alkalmazni a számítási lépések számának csökkentésére. A következő fejezetben megmutatjuk, hogy a problémát át lehet fogalmazni úgy, hogy a hasonló intenzitású képpontokon végzett műveletek jelentős része aggregálható legyen, jelentősen csökkentve azzal az algoritmus számításigényét.

3.2.5.

A javasolt gyors osztályozó algoritmus

Induljunk ki a változó intenzitású zaj modelljét tartalmazó FCM algoritmus (FCM-b) költségfüggvényéből, a (3.9) összefüggésből. A bemeneti kép képpontjainak száma legalább 104 -105 nagyságrendű, míg a színtartományban jellemzően 102 -103 intenzitás szintet különböztetünk meg egymástól. Az FCM-b költségfüggvényének a hagyományos módszerrel

dc_592_12 86 történő optimalizálása során, minden egyes iterációban vonjuk össze, vagyis kezeljük azonos elemként azokat a képpontokat, amelyek a pillanatnyi becsült bk inhomogén intenzitású zaj komponense eltávolítása után azonos intenzitásúaknak tekinthetők. Másképpen fogalmazva, vizsgáljuk meg az yk − bk összetett változó eloszlását az intenzitás tartományban, (t) ami iterációról iterációra változni fog. Jelölje hl azon képpontok számát a t-dik iterációban, amelyeknek kompenzált intenzitása egyenlő l-lel, azaz amelyek esetében igaz az yk − bk = l összefüggés. Belátható, hogy ha Ω(t) -vel jelöljük a t-dik iterációban az yk − bk összetett változó által felvehető értékek halmazát, akkor X

(t)

hl = n

(3.13)

l∈Ω(t) (t)

lesz. Valójában a hl értékek – ahol l ∈ Ω(t) – a kompenzált, zaj nélküli kép hisztogramját írják le a t-dik iterációban. A fenti jelölések szerint a t-dik iterációban aggregálhatjuk azokat a képpontokat, amelyeknek az yk − bk értékük azonos, anélkül, hogy megváltozna a költségfüggvény értéke: JFCM−qb =

c X X

(t)

2 hl um il (l − vi ) ,

(3.14)

i=1 l∈Ω(t)

ahol uil ∈ [0, 1] jelöli a fuzzy tagsági függvények értékeit, amelyek azt írják le, hogy az l kompenzált intenzitású képpont milyen mértékig tartozik az i-dik osztályhoz. Az új költségfüggvény optimalizálási formuláit az alábbi összefüggés parciális deriváltjainak zérus átmeneteiből számítjuk ki Lagrange szorzók alkalmazásával: LFCM−qb = JFCM−qb +

X l∈Ω(t)

λl

c X i=1

1−

c X

!

uil

,

(3.15)

i=1

ahol a λl változók a Lagrange szorzók. Az uil szerinti parciális deriváltakból bármely i = 1 . . . c és l ∈ Ω(t) esetén a következő összefüggést kapjuk: ∂LFCM−qb = 0 ⇒ mum−1 (l − vi )2 − λl = 0 il ∂uil

(3.16)

és következésképpen "

λl uil = m

#−1/(m−1)

(l − vi )−2/(m−1) .

(3.17)

A (3.2) összefüggésben ismertetett valószínűségi megszorítás szerint: c X

"

λl ujl = 1 ⇒ 1 = m j=1

#−1/(m−1) c X

(l − vj )−2/(m−1) .

j=1

(3.18)

dc_592_12 87

Agyi régiók elkülönítése MR felvételeken

A (3.17) és a (3.18) összefüggésekből következik, hogy bármely l ∈ Ω(t) kompenzált intenzitás érték, és bármely i = 1 . . . c indexű osztály esetén, a következő összefüggés alapján számíthatjuk az optimális fuzzy tagsági függvény értékeket: (l − vi )−2/(m−1) uil = P . c (l − vj )−2/(m−1)

(3.19)

j=1 (t)

A fenti összefüggés egyszeri kiszámításával egyszerre határozzuk meg hl darab képpont egy adott osztályhoz tartozási valószínűségét. Megjegyzendő, hogy a kapott tagsági függvény értéke nem függ a hisztogramtól. Az osztályprototípusok szerinti parciális deriváltak zérus átmeneteire a következő összefüggést kapjuk: X (t) ∂LFCM−qb (3.20) = 0 ⇔ −2 hl um il (l − vi ) = 0 , ∂vi (t) l∈Ω amelyet így is felírhatunk X

(t)

hl um il l = vi

X

(t)

hl um il .

(3.21)

l∈Ω(t)

l∈Ω(t)

Ezek alapján bármely i = 1 . . . c indexű osztály esetén az optimális prototípust meghatározó összefüggés: P (t) hl um il l (t) l∈Ω vi = P (t) . (3.22) hl um il l∈Ω(t)

A fenti számítást egy iterációban c alkalommal kell elvégezni, ahogyan azt az FCM-b algoritmus esetében is kellett, de itt a számlálóban és a nevezőben szereplő szummák lényegesen kevesebb összeadandó tagot tartalmaznak. Könnyen belátható, hogy az inhomogén zaj lokális értékét minden képpont esetében külön kell becsülni. Viszont ebben az esetben is felgyorsíthatjuk a számítás végrehajtását azáltal, hogy olyan segédváltozókat vezetünk be, amelyeket a számítások során többször felhasználhatunk. Ezen változókat csak egyszer számoljuk ki, majd elhelyezzük őket egy táblázatban, amelyből a számítások során kikeressük, mikor az adott értékre megint szükségünk van. Ennek megfelelően bármely l ∈ Ω(t) kompenzált, vagyis zajmentes intenzitás érték esetében vezessünk be egy ql segédváltozót a következőképpen c P

ql =

um il vi

i=1 c P

i=1

. um il

(3.23)

dc_592_12 88 3.2. táblázat A gyorsított algoritmus 01 t = 0 . (t=0)

02 Inicializáljuk az additív zajt bk

= 0, ∀k = 1 . . . n . (t=0)

03 Adjunk az osztályprototípusoknak kezdeti értéket: vi

(t=0)

6= vj

, ∀i, j = 1 . . . c; i 6= j .

04 Ismételjük a következő lépéseket 05

t←t+1 ;

06

Számítsuk ki az új hisztogramot hl , ∀l ∈ Ω(t) ;

07

Számoljunk új uil tagsági függvény értékeket, ∀i = 1 . . . c és ∀l ∈ Ω(t) esetén, a (3.19) összefüggés szerint;

08

Számoljunk új vi osztályprototípusokat, ∀i = 1 . . . c, a (3.22) összefüggés szerint;

09

Számoljuk ki a ql , ∀l ∈ Ω(t) segédváltozók új értékeit a (3.23) összefüggés szerint;

10

Számoljunk új bk becsült additív zajt, ∀k = 1 . . . n, a (3.24) összefüggés szerint;

11

Alkalmazzunk átlagoló szűrőt a becsült zaj simítására;

(t)

(t)

(t)

(t)

(t)

12 amíg a konvergencia bekövetkezik, azaz

c P i=1

(t)

(t−1)

|vi − vi

|<ε.

13 Soroljuk be a k-dik képpontot (k = 1 . . . n) abba az osztályba, amelynek indexe arg max{ui,yk −bk , i = 1 . . . c} . i

Ezen segédváltozókkal bármely k (k = 1 . . . n) indexű képpont esetében a zaj becslését a következő összefüggés alapján végezhetjük el: (t)

(t−1)

bk = yk − qlk , ahol lk = yk − bk

.

(3.24)

Az így definiált gyors, zajkompenzációt tartalmazó FCM (FCM-qb) algoritmus végrehajtása során az optimális osztályok meghatározását végző iterációkban felváltva alkalmazzuk a (3.19), (3.22) és (3.24) összefüggéseket mindaddig, amíg ki nem elégítjük a korábban említett konvergencia kritériumot, vagyis az osztályprototípusok intenzitás értékének változása egy ciklus alatt nem lesz egy előre definiált ε értéknél kisebb. A szegmentáláshoz a k-dik képpontot az yk − bk kompenzált, zajmentes intenzitás alapján soroljuk be a hozzá legközelebb eső prototípus által képviselt osztályba. A gyors, zajkompenzációt tartalmazó osztályozó algoritmust a 3.2. táblázat foglalja össze.

3.3.

Az elért eredmények és azok értékelése

Először becsüljük meg a hagyományos és gyorsított eljárás számítási komplexitását. Jelölje ω az Ω halmaz számosságát. A 3.3. táblázatban összefoglaltuk a két FCM algoritmus egy-egy iterációban elvégzendő műveleteinek számítási komplexitását a feldolgozandó adatok függvényében. Ha figyelembe vesszük, hogy a képen jelenlevő egymástól különböző intenzitás értékek száma (ω) legalább két nagyságrenddel kisebb a képpontok számánál (n), egyértelműen állíthatjuk, hogy az első három művelet futási ideje várhatóan jóval kisebb

dc_592_12 89

Agyi régiók elkülönítése MR felvételeken 3.3. táblázat A hagyományos és a javasolt FCM algoritmus számítási komplexitása Művelet

Hagyományos

Gyors

algoritmus (FCM-b)

algoritmus (FCM-qb)

Tagsági függvények újraszámítása

O(nc2 )

O(ωc2 )

Osztályprototípusok újraszámítása

O(nc)

O(ωc)

Additív zaj becslése

O(nc)

O(n + ωc)

Zaj simítása

O(n)

O(n)

–

O(n)

Hisztogram újraszámolása

lesz a hagyományos algoritmus futási idejénél. A javasolt gyorsított algoritmus validációjának következő lépéseként mind a hagyományos, mind pedig a javasolt gyorsított algoritmust leteszteltük mesterségesen előállított fantomok és valódi MR képek szegmentálásával. A fantomokat könnyen osztályozható kétszínű képekhez mesterségesen előállított inhomogén zaj hozzáadásával hoztuk létre. A valódi MR képek az Internet Brain Segmentation Repository [1] adatbázisból származnak. A továbbiakban ismertetett futási időket egy 2GHz-en futó Athlon64 processzoros számítógépen mértük. A 3.1. ábra mindkét algoritmus 24 fantomon lefuttatott tesztjének átlagidejét tartalmazza. Két osztály esetén a javasolt gyors eljárás a kép méretének függvényében 15-25szörös gyorsulást ért el. Szintén megfigyelhető, hogy a gyorsulási arány lassan növekszik a kép méretének növekedésével arányban. A 3.2. ábra egy inhomogén fantom kompenzált hisztogramját mutatja be az iteratív algoritmus különböző fázisaiban. A két osztály tökéletesen szétválaszthatóvá válik 10-12 iteráció alatt. A kompenzálás ugyan nem hozza teljesen helyre a fantomot eredeti állapotára, de az FCM-qb algoritmus hibátlanul szét tudja választani a fantom osztályait. A 3.3. ábra szemlélteti egy inhomogén fantom kompenzálását és szegmentálását. Kompenzálás nélkül az FCM alapú szegmentálás hibás eredményre vezet, azonban az inhomogenitást kompenzálva a javasolt algoritmus a régiókat tökéletesen elkülöníti. Negyven különböző valódi agyi MR képet alkalmaztunk a javasolt algoritmus hatékonyságának és pontosságának felmérésére. A képeket mesterségesen inhomogén zajjal terheltük, majd végrehajtottuk rajtuk a hagyományos FCM-b és a javasolt FCM-qb algoritmust. Az osztályozásnál c = 3 osztály keresését írtuk elő, amelyek a fehérállománynak, a szürkeállománynak és a cerebrospinalis folyadéknak felelnek meg.

dc_592_12 90

3.1. ábra Egy osztályozási ciklus időtartama a hagyományos és gyorsított algoritmussal (felül), és a gyorsulási arány (alul) a képen található képpontok számának (n) függvényében. Ezeket az eredményeket c = 2 osztályra szegmentált fantomképek feldolgozásával kaptuk.

3.2. ábra Egy eredetileg két osztályra szegmentálható inhomogén intenzitású zajt tartalmazó fantomkép hisztogramja az optimalizálás különböző iterációiban.

3.3. ábra Fantom osztályozása két osztályra: (a) szinuszos inhomogén zajjal terhelt eredeti kép, (b) osztályozás kompenzálás nélkül, (c) osztályozás az inhomogén zaj kompenzálásával.


91

3.4. ábra Valós MR kép osztályozása három osztályra: (a) szimulált zajjal terhelt eredeti kép, (b) osztályozás kompenzálás nélkül, (c) osztályozás az inhomogén zaj kompenzálásával.

3.5. ábra Egy optimalizálási ciklus időtartama a hagyományos és gyorsított algoritmussal (felül), és a gyorsulási arány (alul), a képen található képpontok számának (n) függvényében. Ezeket az eredményeket c = 3 osztályra szegmentált valós MR képek feldolgozásával kaptuk.

Az osztályozás eredményére látunk példát a 3.4 ábrán. A futási időket a 3.5. ábra foglalja össze. A javasolt algoritmus alkalmazásával a szegmentálás gyorsulásának aránya magasabb, mint azt a fantom képek esetében tapasztaltuk, az új algoritmus 28-32-szer gyorsabban hajtotta végre a szegmentálást, mint a hagyományos FCM algoritmus. A 3.6. (a)-(e) ábra egy agyi MR kép hisztogramjának változását mutatja be a zaj kompenzálása során. A konvergencia általában 60-100 iteráció alatt következik be. Mivel a fehér- és szürkeállomány intenzitása közel áll egymáshoz, ezek hisztogramjai részben átfedik egymást, és emiatt lehetetlen őket tökéletesen szétválasztani. Ez a hibás osztályozások elsődleges forrása, amely egyformán jelen van a hagyományos és a gyorsított eljárás esetében. A 3.6.(f) ábra bemutatja a három osztályprototípus intenzitás értékének változását az

dc_592_12 92

3.6. ábra (a)-(e) Egy inhomogén intenzitású zajjal terhelt agyi MR kép hisztogramja az optimalizálás különböző iterációiban; (f) Az osztályprototípusok változása az iterációk során. iterációk során. Mint látható, a konvergencia a 60-dik iteráció környékén következett be. Mivel c ω n, a 3.3. táblázatban bemutatott számítási komplexitás értékek alapján azt feltételeztük, hogy a gyorsított algoritmus futásideje kisebb mértékben növekszik az osztályok számának emelkedésével, mint a hagyományos algoritmus futási ideje. Ennek igazolására a hagyományos, illetve a javasolt gyors algoritmust lefuttattuk különböző (2től 8-ig változó) számú osztályra történő szegmentálási feladatokon. A 3.7. ábra összegzi a kapott gyorsulási arányokat. Ezen a grafikonon azt látjuk, hogy a gyorsulási arány valóban növekszik, méghozzá nagyjából egyenes arányban az osztályok számával. A 3.8. ábra megmutatja, hogy a hagyományos és a javasolt gyorsított algoritmusnak

3.7. ábra A gyorsulási arány átlagának és szórásának alakulása az osztályok számának függvényében.


93

3.8. ábra Egy optimalizálási ciklus részfeladatainak relatív időigénye az FCM-b és az FCMqb algoritmusok esetén. egyetlen iterációjában hogyan aránylik az egyes számítási lépések relatív időigénye. Ebből az ábrázolásból is kitűnik a gyorsulás ténye: a becsült inhomogenitást simító átlagoló szűrő mindkét esetben ugyanaz a művelet, vagyis annak futási ideje mindkét esetben azonos, azonban míg a hagyományos algoritmus esetében ez az idő a teljes futási idő 1%-nál kisebb, addig a gyorsított algoritmus esetén ez a teljes futási idő 19%-a. A hagyományos algoritmus esetén a partíció újraszámolása tart a leghosszabb ideig, míg a gyorsított algoritmusban a képpontonként számolt zajbecslés és a hisztogram frissítése a leginkább időigényes feladat. Mivel a javasolt gyors algoritmus a hagyományos algoritmussal azonos osztályozási elvre épül, ugyanazokat a számításokat végzi el, emiatt az osztályozás eredményében és pontosságában elvileg nem lehet eltérés. A gyakorlatban azonban létezik egy másodlagos hibaforrás a gyorsított algoritmus végrehajtása során. A gyors algoritmus esetén minden iterációban a becsült zajt kvantálni kell annak érdekében, hogy meg tudjuk határozni, mely képpontokon végzett műveletek kezelhetők közösen. A hagyományos algoritmus esetében erre nincs szükség. Az algoritmus tesztelése során megvizsgáltuk, hogy az iterációnkénti kvantálás miatt mekkora az eltérés az osztályozás pontosságában. A 3.9. ábra mutatja az összehasonlítás eredményét. Láthatjuk, hogy a gyorsított módszer osztályozási hibáinak száma jóval kevesebb, mint 1%-al haladja csak meg a hagyományos megoldás hibáinak számát, tehát a gyakori kvantálás miatt keletkezett hiba elhanyagolható. A javasolt gyors osztályozási algoritmus további, gyakorlati alkalmazás szempontjából előnyös tulajdonságokkal is rendelkezik: • Kompatibilis a 3.2.2. fejezetben bemutatott összes zajmodellel. • A javasolt osztályozási eljárás társítható a szakirodalomban leírt bármelyik szűrési eljárással (pl. [78, 86]).

dc_592_12 94

3.9. ábra A hibás osztályozások arányának alakulása egy valós agyi MR kép esetében, a képen található képpontok számának függvényében.

3.4.

Összefoglalás

A fejezetben az MR felvételeken alkalmazható, hagyományos FCM algoritmus inhomogén, térben lassan változó intenzitású zajt kompenzáló és képszegmentáló módszerének olyan új változatát, amely az eredeti algoritmusnál lényegesen kisebb számítási komplexitású algoritmust eredményezett. A javasolt algoritmust implementáltam, majd benchmark kísérletekkel megmutattam, hogy a javasolt algoritmus alkalmazásával akár két nagyságrenddel is csökkenhet a nagy számításigényű részműveletek futási ideje úgy, hogy közben a szegmentálás minősége nem romlik. Az agyról készült MR képek zaj kompenzálási és szegmentálási feladatának végrehajtása során a javasolt gyors algoritmus alkalmazása 25-30-szoros gyorsulást eredményezett a hagyományos FCM megoldáshoz képest anélkül, hogy számottevően csökkent volna a szegmentálás pontossága. Továbbfejlesztési lehetőségek Az ismertetett kutatás legkönnyebben megvalósítható folytatása lehet a javasolt módszer alkalmazása különböző típusú FCM algoritmusokkal, pl. annak integrálása lehetőségfüggvényeket alkalmazó vagy hibrid c-means algoritmusokkal.


3.5.

95


3. tézis: Új csoportos osztályozáson alapuló módszert adtam egycsatornás intenzitás képek gyors feldolgozására, amely alkalmazható térben lassan változó intenzitású zajjal terhelt felvételek szegmentálása során. A módszer alapján kidolgoztam a zajos MR koponyafelvételeken agyi régiók elkülönítésére szolgáló fuzzy c-means (FCM) szegmentáló algoritmus gyorsított változatát. A javasolt szegmentáló algoritmust mesterségesen előállított fantom képeken, illetve valós MR felvételeken validáltam, valamint összehasonlítottam annak futási idejét a hagyományos FCM szegmentáló algoritmus futási idejével. Beigazolódott, hogy a javasolt algoritmus futási ideje a hagyományos algoritmus futási idejénél legalább egy nagyságrenddel kisebb amellett, hogy a hagyományos algoritmussal azonos minőségű szegmentálást biztosít.

3.5.1.


A fejezetben ismertetett eredményeket bemutató legfontosabb publikációk: [85], [89], [88], [86].

dc_592_12 96

dc_592_12

4. fejezet Eredmények hasznosítása A gyakorlati hasznosítás szempontjából az ismertetett kutatás legfontosabb eredményei a következők:

• SPECT képrekonstrukciós algoritmus: Képrekonstrukciós eszközökbe közvetlenül beépíthető, a rendelkezésre álló hardver erőforrásokat hatékonyan kihasználó, korábbi megoldásokkal összehasonlítva azoknál lényegesen jobb minőségű képet eredményező képrekonstrukciós eljárás.

• Foggyökér csatorna felismerésére kidolgozott módszerek: Olyan, különböző CT modalitásokon működő két képfeldolgozó eljárás, amelyek automatikus működésüket, illetve gyakorlatban történő alkalmazhatóságukat valós felvételeken igazolták.

• Agyi területek szegmentálásra kidolgozott algoritmus: Olyan módszer, illetve az alapján létrehozott algoritmus, amely a hagyományos FCM szegmentáló algoritmusnál legalább egy nagyságrenddel rövidebb idő alatt képes végrehajtani az osztályozás feladatát a hagyományos algoritmussal megegyező pontossággal.

Ebben a fejezetben ismertetem, hogy ezen eredmények mely gyakorlati alkalmazásokban, illetve mely kutatási, valamint kutatásfejlesztési projektekben kerültek felhasználásra és hasznosításra. 97

dc_592_12 98

4.1.

A kidolgozott képrekonstrukciós algoritmus hasznosítása

A kidolgozott SPECT képrekonstrukciós algoritmus beépítésre került a MEDISO Orvosi Berendezés Fejlesztő és Szerviz Kft. InterViewXP alkalmazásába [60]. Az alkalmazás felhasználói választhatják a detektált projekciós képek rekonstrukcióját az 1. tézispontban bemutatott rekonstrukciós algoritmussal. Az algoritmus beépítésre került a Mediso Kft. szabványos rekonstrukciós szoftver könyvtárába, így a jelenleg fejlesztés alatt álló SPECT berendezések a javasolt algoritmus által megvalósított rekonstrukciós opcióval is rendelkeznek. A kutatás eredményei bevezetésre kerültek az egyetemi oktatásba. Az eredmények bemutatásra kerülnek a BME Villamosmérnöki és Informatikai Kar Műszaki és biológiai rendszerek elmélete, valamint Bevezetés az egészségügyi mérnöki tudományokba című tárgyak keretében. A képrekonstrukciós algoritmusok kidolgozása témában született eredmények a következő alapkutatási és kutatásfejlesztési projektek keretében kerültek alkalmazásra: • 2008 – 2011, TeraTomo, NKTH: TECH_08_A2 (2008) NTP program: Különböző modalitású orvosi diagnosztikai tomográfiás berendezésekbe építhető teraflop kapacitású képrekonstrukciós rendszer kifejlesztése (BME konzorciumi projektvezető: Benyó Balázs). Konzorciumi tagok: – MEDISO Orvosi Berendezés Fejlesztő és Szerviz Kft., – Budapesti Műszaki és Gazdaságtudományi Egyetem ∗ Irányítástechnika és Informatika Tanszék, ∗ Nukleáris Technika Intézet, – Semmelweis Egyetem, Radiológiai és Onkoterápiás Klinika. • 2006 – 2009, PETCT_06, Jedlik Ányos program: Multi-modalitású képalkotórendszer sorozatgyártásra történő kifejlesztése orvos-biológiai kutatás és humán orvos-diagnosztika céljára (BME konzorciumi projektvezető: Benyó Balázs). Konzorciumi tagok: – MEDISO Orvosi Berendezés Fejlesztő és Szerviz Kft., – Debreceni Egyetem, Nukleáris Medicina Intézet,

dc_592_12 Eredmények hasznosítása

99

– MTA Atommag Kutató Intézet, Debrecen, – Budapesti Műszaki és Gazdaságtudományi Egyetem ∗ Irányítástechnika és Informatika Tanszék, ∗ Atomfizika Tanszék.

• 2004 – 2008, OTKA 46726: Biztonságkritikus diagnosztikai célú informatikai rendszerek kutatása (témavezető: Benyó Balázs).

4.2.

Foggyökér csatorna felismerésére kidolgozott módszerek hasznosítása

A kutatás eredményei bevezetésre kerültek az egyetemi oktatásba. A Semmelweis Egyetem Fogorvostudományi Karon oktatott Endodoncia című tárgy keretében a kidolgozott módszer segítségével a hallgatók számára lehetővé vált a gyökércsatorna térbeli szerkezetének megismertetése, amely korábban csak síkbeli vetületi képek alapján volt lehetséges. Jelenleg folyik a cone beam CT-re kidolgozott módszer gyakorlati alkalmazásának bevezetése, amely a meglevő cone beam CT készülékkel készített 3D képek feldolgozásával lényegesen megváltoztathatja a foggyökér-betegségek kezeléséhez kapcsolódó beavatkozásokat. Az eredmények ugyancsak bemutatásra kerülnek a BME Villamosmérnöki és Informatikai Kar Műszaki és biológiai rendszerek elmélete, valamint Bevezetés az egészségügyi mérnöki tudományokba című tárgyak keretében. A foggyökér csatorna felismerése és rekonstrukciója témában született eredmények a következő alapkutatási és kutatásfejlesztési projektek keretében kerültek alkalmazásra:

• 2010 – 2014, OTKA K80266: Új módszerek kidolgozása az orvosi diagnosztika hatékonyságának növelésére (Témavezető: Benyó Balázs),

• 2010 – 2012, Minőségorientált, összehangolt oktatási és K+F+I stratégia, valamint működési modell kidolgozása a Műegyetemen (TÁMOP-4.2.1/B-09/1/KMR-20100002): IKT-P1-T9: Modell alapú mérnöki módszerek kidolgozása orvosi és műszaki alkalmazásokhoz (Altémavezető: Benyó Balázs).

dc_592_12 100

4.3.

Agyi területek szegmentálására kidolgozott algoritmus hasznosítása

A kutatás eredményei bevezetésre kerültek az egyetemi oktatásba. Az eredmények bemutatásra kerülnek a BME Villamosmérnöki és Informatikai Kar Műszaki és biológiai rendszerek elmélete, valamint Bevezetés az egészségügyi mérnöki tudományokba című tárgyak keretében. Az agyi területek szegmentálása témában született eredmények a következő alapkutatási projektek keretében kerültek alkalmazásra: • 2003 – 2007, OTKA F046726: Biztonságkritikus diagnosztikai célú informatikai rendszerek kutatása (Témavezető: Benyó Balázs), • 2006 – 2010, OTKA T69055: Új mérési, szabályozási eljárások kidolgozása, orvosinformatikai alkalmazása betegségek korai diagnosztizálására és az optimális terápia megvalósítására (Témavezető: Benyó Zoltán – Benyó Balázs résztvevő), • 2010 – 2014, OTKA K80266: Új módszerek kidolgozása az orvosi diagnosztika hatékonyságának növelésére (Témavezető: Benyó Balázs).

dc_592_12

Köszönetnyilvánítás Köszönöm Dr. Szirmay-Kalos László, egyetemi tanár, a BME Irányítástechnika és Informatika Tanszék jelenlegi tanszékvezetőjének, és Dr. Arató Péter, egyetemi tanár, korábbi tanszékvezetőjének odafigyelő támogatását. Köszönet illeti a – jelenleg általam, korábban édesapám, Dr. Benyó Zoltán által vezetett – Orvosinformatikai Laboratórium minden jelenlegi és volt tagjának a kutatómunkámhoz nyújtott segítségét. Külön köszönöm Szlávecz Ákosnak, Hesz Gábornak és Szilágyi Lászlónak a közreműködését, akik közvetlenül is hozzájárultak kutatásaim sikeréhez. Hasonlóan hálás vagyok a kutatási programokban résztvevő partnerintézmények kutatóinak, akik közül Dr. Dobó-Nagy Csaba, Dr. Bükki Tamás és Dr. Kári Béla nevét szeretném feltétlenül megemlíteni. Köszönöm Dr. Paláncz Béla hathatós segítségét, aki szakmai tanácsaival segített az értekezés elkészítésében. A családom állhatatos támogatása nélkül nem tudtam volna ezeket az eredményeket elérni, köszönöm elsősorban szüleimnek és feleségemnek, Katinak, valamint lányaimnak Júliának és Dorottyának az ösztönzést és a türelmet, amellyel elviselték távollétemet. Végül, de nem utolsó sorban szeretnék köszönetet mondani a BME Irányítástechnika és Informatika Tanszék valamennyi munkatársának azért a baráti légkörért, amelyben öröm volt dolgozni.

101

dc_592_12

Irodalomjegyzék [1] Internet brain segmentation repository. http://www.cma.mgh.harvard.edu/ibsr. [2] SOTE, ÁOK, Radiológiai és Onkoterápiás Klinika. intezetek/?inst_id=55.

http://www.sote.hu/

[3] Biomed Virtual Organization, 2009. http://wiki.healthgrid.org/LSVRC:Biomed. [4] Hungrid Virtuális Szervezet, 2009. http://www.lcg.kfki.hu/. [5] R. Adams and L. Bischof. Seeded region growing. IEEE T Pattern Analisys and Machine Intelligence, 16(6):641–647, 1994. [6] Mohamed N. Ahmed, Sameh M. Yamany, Nevin Mohamed, Aly A. Farag, and Thomas Moriarty. A modified fuzzy c-means algorithm for bias field estimation and segmentation of mri data. Medical Imaging, IEEE Transactions on, 21(3):193–199, 2002. [7] Ákos Szlávecz and Balázs Benyó. Optimal scanning protocol for 180 degree data acquisition in parallel spect imaging. In Proc. of IFAC BMS 2012 – 8th IFAC Symposium on Biological and Medical Systems, pages 12–17, Budapest, 2012. ISBN: 978-3-902823-10-6. [8] Mostafa Analoui, Satthya Krisnamurthy, and Cecil Brown. Modeling and measurement of root canal using stereo digital radiography. In Proceedings of SPIE, volume 3976, pages 306–314, 2000. [9] Oscar Kin-Chung Au, Chiew-Lan Tai, Hung-Kuo Chu, Daniel Cohen-Or, and TongYee Lee. Skeleton extraction by mesh contraction. In ACM Transactions on Graphics (TOG), volume 27, page 44. ACM, 2008. [10] Leon Axel, Jay Costantini, and John Listerud. Intensity correction in surface-coil mr imaging. American Journal of Roentgenology, 148(2):418–420, 1987. i

dc_592_12 ii [11] Chuanyong Bai, Gengsheng L. Zeng, G.T. Gullberg, F. DiFilippo, and S. Miller. Slab-by-slab blurring model for geometric point response correction and attenuation correction using iterative reconstruction algorithms. Nuclear Science, IEEE Transactions on, 45(4):2168–2173, aug 1998. [12] Guillaume Bal and Philippe Moireau. Fast numerical inversion of the attenuated radon transform with full and partial measurements. Inverse Problems, 20:1137– 1164, 2004. [13] Balázs Benyó, László Szilágyi, Zsolt Németh, Gábor Molnár Csaba, and Csaba DobóNagy. Identification of the root canal and its centreline from dental cone beam ct records. In Proc. of IFAC BMS 2012 – 8th IFAC Symposium on Biological and Medical Systems, pages 1–5, Budapest, 2012. ISBN: 978-3-902823-10-6. [14] Kári Béla, Szlávecz Ákos, Hesz Gábor, Bukki Tamás, Pártos Oszkár, Gyorke Tamás, and Benyó Balázs. 3d gyors iteratív kép rekonstrukciós eljárás parallel leképezésű spect modalitáshoz. Magyar Radiologia, 86(1):26–27, 2012. [15] S. Bellini, M. Piacentini, C. Cafforio, and F. Rocca. Compensation of tissue absorption in emission tomography. Acoustics, Speech and Signal Processing, IEEE Transactions on, 27(3):213–218, jun 1979. [16] B Benyó, L Szilágyi, and C Dobó-Nagy. Root canal and center line identification from cone-beam CT volume data. International Journal of Computer Assisted Radiology and Surgery, 5(Suppl 1):S236–S237, 2010. [17] B Benyó, L Szilágyi, Cs G Molnár, Z Németh, and Cs Dobó-Nagy. Automatic identification of the root canal from CBCT. International Journal of Computer Assisted Radiology and Surgery, 6(suppl 1):S210–S211, 2011. [18] Balázs Benyó. Identification of dental root canals and their medial line from micro-ct and cone-beam ct records. BioMedical Engineering OnLine, 11(81), October 2012. doi: 10.1186/1475-925X-11-81. [19] Balázs Benyó, László Szilágyi, and Csaba Dobó-Nagy. A skeletal approach to root canal centreline detection from dental micro-CT records. In CONTROL 2010 – UKACC International Conference on Control 2010, Coventry, pages 1–6, Coventry, 2010. [20] S. Bergner, E. Dagenais, A. Celler, and T. Moller. Using the physics-based rendering toolkit for medical reconstruction. In Nuclear Science Symposium Conference Record, 2005 IEEE, volume 5, pages 3022–3026, oct. 2005.

dc_592_12 IRODALOMJEGYZÉK

iii

[21] J.C. Bezdek. Pattern recognition with fuzzy objective function algorithms. Kluwer Academic Publishers, 1981. [22] J.C. Bezdek, Hall L.O., and Clarke L.P. Review of mr image segmentation techniques using pattern recognition. Medical Physics, 20(4):1033–1048, 1993. [23] Biodex Medical Systems Inc. Lung-Spine SPECT Phantom. http://www.biodex. com/nuclear-medicine/products/phantoms/lung-spine-spect-phantom. [24] B.H. Brinkmann, A. Manduca, and R.A. Robb. Optimized homomorphic unsharp masking for mr grayscale inhomogeneity correction. IEEE Transactions on Medical Imaging, 17(2):161–171, 1998. [25] W. Cai, S. Chen, and D. Zhang. Fast and robust fuzzy c-means clustering algorithms incorporating local information for image segmentation. Pattern Recognition, 40(3):825–838, 2007. [26] Lee-Tzuu Chang. A method for attenuation correction in radionuclide computed tomography. Nuclear Science, IEEE Transactions on, 25(1):638–643, feb. 1978. [27] Nicu D Cornea, Deborah Silver, and Patrick Min. Curve-skeleton properties, applications, and algorithms. Visualization and Computer Graphics, IEEE Transactions on, 13(3):530–548, 2007. [28] Nicu D Cornea, Deborah Silver, Xiaosong Yuan, and Raman Balasubramanian. Computing hierarchical curve-skeletons of 3d objects. The Visual Computer, 21(11):945– 955, 2005. [29] J. Derganc, B. Likar, and F. Pernuš. Nonparametric segmentation of multispectral mr images incorporating spatial and intensity information. Progress in Biomedical Optics and Imaging, 3(1):391–400, 2002. [30] Cs. Dobó-Nagy, G. Keszthelyi, J. Szabó, P. Sulyok, G. Ledeczky, and J. Szabó. A computerized method for mathematical description of three-dimensional root canal axis. J Endodontics, 26(11):639–643, 2000. [31] Cs. Dobó-Nagy and J. Szabó. A mathematically based classification of root canal curvatures on natural human teeth. J Endodontics, 21(11):557–560, 1995. [32] J Dong, SY Hong, and G Hasselgren. Theories and algorithms for 3-d root canal model construction. Computer-Aided Design, 37(11):1177–1189, 2005.

dc_592_12 iv [33] P. R. Edholm, R. M. Lewitt, and B. Lindholm. Novel properties of the fourier decomposition of the sinogram. In In. Proc Int. Workshop Physics and Engineering of Computerized Multidimensional Imaging and Processing, volume 671, pages 8–18, 1986. [34] Maki Endo, Syoji Kobashi, Katsuya Kondo, and Yutaka Hata. Dentistry support ultrasonic system for root canal treatment aided by fuzzy logic. In Systems, Man and Cybernetics, 2005 IEEE International Conference on, volume 2, pages 1494– 1499. IEEE, 2005. [35] European Commission. Radiation Protection No. 172. Cone Beam CT for Dental and Maxillofacial Radiology: Evdence Based Guidelines, 2012. [36] Desmond M Germans, Hans JW Spoelder, Luc Renambot, Henri E Bal, Sander van Daatselaar, and Paul van der Stelt. Measuring in virtual reality: a case study in dentistry. Instrumentation and Measurement, IEEE Transactions on, 57(6):1177– 1184, 2008. [37] Frisardi Gianni, Chessa Giacomo, Barone Sandro, Paoli Alessandro, Razionale Armando, and Frisardi Flavio. Integration of 3d anatomical data obtained by ct imaging and 3d optical scanning for computer aided implant surgery. BMC Medical Imaging, 11(5):1–7. [38] S.J. Glick, B.C. Penney, M.A. King, and C.L. Byrne. Noniterative compensation for the distance-dependent detector response and photon attenuation in spect imaging. Medical Imaging, IEEE Transactions on, 13(2):363–374, jun 1994. [39] GPGPU. General-purpose computation on graphics processing units, 2011. http: //www.gpgpu.org/. [40] Shane Y Hong and Janet Dong. 3d root canal modeling for advanced endodontic treatment. In NDE For Health Monitoring and Diagnostics, pages 321–330. International Society for Optics and Photonics, 2002. [41] H.M. Hudson and R.S. Larkin. Accelerated image reconstruction using ordered subsets of projection data. Medical Imaging, IEEE Trans. on, 13(4):601–609, December 1994. [42] S Jan et al. Gate: a simulation toolkit for pet and spect. Physics in Medicine and Biology, 49(19):4543, 2004.


v

[43] B. Johnston, Atkins M.S., Mackiewich B., and Anderson M. Segmentation of multiple sclerosis lesions in intensity corrected multispectral mri. IEEE Transactions on Medical Imaging, 15(2):154–169, 1996. [44] Béla Kári, Ákos Szálvecz, Gábor Hesz, Tamás Bűkki, Oszkár Pártos, Tamás Győrke, and Balázs Benyó. Novel high speed 3d iterative reconstruction algorithm for parallel projection based spect imaging. In 10th Asia-Oceania Nuclear Medicine and Biology Congress, page 45, Teheran, Iran, 2012. Paper OP066. [45] Sadullah Kaya, Ozkan Adiguzel, Izzet Yavuz, Emin-Caner Tumen, and Zeki Akkus. Cone-beam dental computerized tomography for evaluating changes of aging in the dimensions central superior incisor root canals. Med Oral Patol Oral Cir Bucal, 16(3):e463–e466, 2011. [46] M.A. King, D.J. DeVries, T.-S. Pan, P.H. Pretorius, and J.A. Case. An investigation of the filtering of tew scatter estimates used to compensate for scatter with ordered subset reconstructions. Nuclear Science, IEEE Transactions on, 44(3):1140–1145, jun 1997. [47] Michael A. King, Paul W. Doherty, Ronald B. Schwinger, and Bill C. Penney. A wiener filter for nuclear medicine images. Medical Physics, 10(6):876–880, 1983. [48] Kenneth F. Koral, Fayez M. Swailem, Steven Buchbinder, Neal H. Clinthorne, W. Leslie Rogers, and Benjamin M.W. Tsui. Spect dual-energy-window compton correction: Scatter multiplier required for quantification. Journal of Nuclear Medicine, 31(1):90–98, 1990. [49] Raghu Krishnapuram and James M. Keller. Possibilistic approach to clustering. IEEE Transactions on Fuzzy Systems, 1(2):98–110, 1993. [50] Vipin Kumar. Introduction to Parallel Computing. Addison-Wesley Longman Publishing Co., Inc., Boston, MA, USA, 2nd edition, 2002. [51] Kennth Lange and Richard Carson. EM reconstruction algorithms for emission and transmission tomography. Journal of Computer Assisted Tomography, 8(2):306–316, April 1984. [52] Jong-Ki Lee, Byung-Hyun Ha, Jeong-Ho Choi, Seok-Mo Heo, and Hiran Perinpanayagam. Quantitative three-dimensional analysis of root canal curvature in maxillary first molars using micro-computed tomography. Journal of Endodontics, 32(10):941–945, 2006.

dc_592_12 vi [53] Robert M. Lewitt, Paul R. Edholm, and Weishi Xia. Fourier method for correction of depth-dependent collimator blurring. pages 232–243, 1989. [54] Hua Li, Anthony Yezzi, and Laurent Cohen. 3d multi-branch tubular surface and centerline extraction with 4d iterative key points. In Medical Image Computing and Computer-Assisted Intervention–MICCAI 2009, pages 1042–1050. Springer, 2009. [55] Alan Wee-Chung Liew and Hong Yan. An adaptive spatial fuzzy clustering algorithm for 3-d mr image segmentation. Medical Imaging, IEEE Transactions on, 22(9):1063– 1075, 2003. [56] C Min Ma and Milan Sonka. A fully parallel 3d thinning algorithm and its applications. Computer vision and image understanding, 64(3):420–433, 1996. [57] M Magdics, B Tóth, L Szécsi, B Csébfalvi, L Szirmay-Kalos, A Szlavecz, G Hesz, B Benyó, A Cserkaszky, D Légrády, et al. Detector modeling techniques for preclinical 3d pet reconstruction on the gpu. Full 3-D Concerence, pages 375–378, 2011. [58] G. Mariani, L. Bruselli, T. Kuwert, E.E. Kim, A. Flotats, O. Israel, M. Dondi, and N. Watanabe. A review on the clinical uses of spect/ct. Eur J Nucl Med Mol Imaging, 37:1959–1985, 2010. [59] A.W. McCarthy and M.I. Miller. Maximum likelihood spect in clinical computation times using mesh-connected parallel computers. Medical Imaging, IEEE Transactions on, 10(3):426–436, sep 1991. [60] Mediso. Interview xp software package, 2010. http://www.mediso.hu/. [61] F Natterer. Inversion of the attenuated radon transform. Inverse Problems, 17(1):113, 2001. [62] Frank Natterer. The mathematics of computerized tomography. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2001. [63] B Neelima and Prakash S Raghavendra. Recent trends in software and hardware for gpgpu computing: a comprehensive survey. In Industrial and Information Systems (ICIIS), 2010 International Conference on, pages 319–324. IEEE, 2010. [64] Gregory M Nielson and Bernd Hamann. The asymptotic decider: resolving the ambiguity in marching cubes. In Proceedings of the 2nd conference on Visualization’91, pages 83–91. IEEE Computer Society Press, 1991.


vii

[65] RomanG. Novikov. An inversion formula for the attenuated x-ray transformation. Arkiv för Matematik, 40:145–167, 2002. [66] nVidia. Cuda programming guide. http://developer.nvidia.com/. [67] Kálmán Palágyi and Attila Kuba. A parallel 3d 12-subiteration thinning algorithm. Graphical Models and Image Processing, 61(4):199–221, 1999. [68] Jong-Wook Park, Jong-Ki Lee, Byung-Hyun Ha, Jeong-Ho Choi, and Hiran Perinpanayagam. Three-dimensional analysis of maxillary first molar mesiobuccal root canal configuration and curvature using micro–computed tomography. Oral Surgery, Oral Medicine, Oral Pathology, Oral Radiology, and Endodontology, 108(3):437–442, 2009. [69] Ove A Peters. Current challenges and concepts in the preparation of root canal systems: a review. Journal of Endodontics, 30(8):559–567, 2004. [70] Dzung L. Pham and Jerry L. Prince. Adaptive fuzzy segmentation of magnetic resonance images. Medical Imaging, IEEE Transactions on, 18(9):737–752, 1999. [71] Chris Pudney. Distance-ordered homotopic thinning: a skeletonization algorithm for 3d digital images. Computer Vision and Image Understanding, 72(3):404–413, 1998. [72] Jagath C Rajapakse and Frithjof Kruggel. Segmentation of mr images with intensity inhomogeneities. Image and Vision Computing, 16(3):165–180, 1998. [73] William C Scarfe, Martin D Levin, David Gane, and Allan G Farman. Use of cone beam computed tomography in endodontics. Int J Dent, 2009(634567), 2009. 1-20. [74] W.P. Segars, D.S. Lalush, and B.M.W. Tsui. A realistic spline-based dynamic heart phantom. Nuclear Science, IEEE Transactions on, 46(3):503–506, jun. 1999. [75] Y. Seo, C Mari, and BH Hasegawa. Technological development and advances in single-photon emission computed tomography/computed tomography. Semin Nucl Med, 38:177–198, 2008. [76] L. A. Shepp and Y. Vardi. Maximum likelihood reconstruction for emission tomography. Medical Imaging, IEEE Transactions on, 1(2):113–122, oct. 1982. [77] A. Simmons, Tofts P.S., Barker G.J., and Arridge S.R. Sources of intensity nonuniformity in spin echo images at 1.5 t. Magnetic Resonance in Medicine, 32(1):121–128, 1994.

dc_592_12 viii [78] MY Siyal and Lin Yu. An intelligent modified fuzzy c-means based algorithm for bias estimation and segmentation of brain mri. Pattern Recognition Letters, 26(13):2052– 2062, 2005. [79] John G Sled, Alex P. Zijdenbos, and Alan C. Evans. A nonparametric method for automatic correction of intensity nonuniformity in mri data. Medical Imaging, IEEE Transactions on, 17(1):87–97, 1998. [80] M. Styner. Parametric estimate of intensity inhomogeneities applied to mri. IEEE Transactions on Medical Imaging, 19(3):153–165, 2000. [81] Michael V Swain and Jing Xue. State of the art of micro‐ ct applications in dental research. International journal of oral science, 1(4):177–188, 2009. [82] L Szilagyi, Z Benyo, SM Szilagyi, and HS Adam. Mr brain image segmentation using an enhanced fuzzy c-means algorithm. In Engineering in Medicine and Biology Society, 2003. Proceedings of the 25th Annual International Conference of the IEEE, volume 1, pages 724–726. IEEE, 2003. [83] L. Szilágyi, S. M. Szilágyi, and Z. Benyó. A modified fuzzy c-means algorithm for mr brain image segmentation. Lect Notes Comp Sci, 4633:866–877, 2007. [84] L. Szilágyi, S. M. Szilágyi, D. Iclănzan, and L. Szabó. Efficient 3d curve skeleton extraction from large objects. Lect Notes Comp Sci, 7042:133–140, 2011. [85] L. Szilágyi, S.M. Szilágyi, and B. Benyó. Efficient inhomogeneity compensation using fuzzy c-means clustering models. Computer Methods and Programs in Biomedicine, 108(1):80–89, 2012. doi: 10.1016/j.cmpb.2012.01.005. [86] László Szilágyi, László Dávid, Sándor M Szilágyi, Balázs Benyó, and Zoltán Benyó. Improved intensity inhomogeneity correction techniques in mr brain image segmentation. In Proceedings of the 17 th World Congress The International Federation of Automatic Control, pages 9625–9630, 2008. [87] László Szilágyi, Csaba Dobó-Nagy, and Balázs Benyó. Identification of the root canal from dental micro-ct records. In César San Martin and Sang-Woon Kim, editors, Progress in Pattern Recognition, Image Analysis, Computer Vision, and Applications, volume 7042 of Lecture Notes in Computer Science, pages 339–346. Springer Berlin Heidelberg, 2011. [88] László Szilágyi, Sándor M Szilágyi, Balázs Benyó, and Zoltán Benyó. Application of hybrid c-means clustering models in inhomogeneity compensation and mr brain


ix

image segmentation. In Modeling and Control in Biomedical Systems, volume 7, pages 204–209, 2009. [89] László Szilágyi, Sándor M. Szilágyi, Balázs Benyó, and Zoltán Benyó. Intensity inhomogeneity compensation and segmentation of MR brain images using hybrid c-means clustering models. Biomedical Signal Processing and Control, 6(1):3–12, January 2011. [90] Á Szlávecz, B Benyó, G Hesz, and T Bükki. An optimized SPECT reconstruction algorithm with attenuation correction and resolution recovery. International Journal of Computer Assisted Radiology and Surgery, 4(Suppl 1):S353–S354, 2009. [91] A. Szlávecz, B. Benyó, P. Várady, and T. Bukki. Reconstruction of myocardial short-scan SPECT images. In Intelligent Engineering Systems, 2005. INES ’05. Proceedings. 2005 IEEE International Conference on, pages 37–41, 2005. [92] Á Szlávecz, T Bükki, C Steinbach, and B Benyó. A novel model-based PET detector block simulation approach. Biomedical Signal Processing and Control, 6(1):27–33, 2011. [93] Á Szlávecz, G Hesz, Z Puskás, B Kári, O Pártos, T Györke, T Bükki, B Domonkos, and B Benyó. A fast iterative GPU-based reconstruction algorithm for SPECT imaging involving collimator and attenuation compensation. European Journal of Nuclear Medicine and Molecular Imaging, 37(Suppl 2):S284, 2010. [94] Ákos Szlávecz, Gábor Hesz, Tamás Bükki, Béla Kári, and Balázs Benyó. GPU-based acceleration of the MLEM algorithm for SPECT parallel imaging with attenuation correction and compensation for detector response. In Proceedings of the 18th IFAC World Congress. Milan, Italy., pages 6195–6200, August 28. - September 2. 2011. [95] Ákos Szlávecz, Zoltán Puskás, and Balázs Benyó. GPGPU based acceleration of the single scatter simulation algorithm for positron emission tomography. Buletinul Stiimtific al Universitattii Politehnica din Timisoara-seria Automatica si Calculatoare, 55(69)(4):185–195, 2010. [96] Yasuyuki Takahashi, Kenya Murase, Teruhito Mochizuki, Yoshifumi Sugawara, Hisato Maeda, and Akiyoshi Kinda. Simultaneous 3-dimensional resolution correction in spect reconstruction with an ordered-subsets expectation maximization algorithm. Journal of Nuclear Medicine Technology, 35(1):34–38, 2007. [97] Sergios Theodoridis. Pattern recognition. Elsevier, London, UK, 2nd edition edition, 2003.

dc_592_12 x [98] B M W Tsui, E C Frey, Xide Zhao, D S Lalush, R E Johnston, and W H McCartney. The importance and implementation of accurate 3d compensation methods for quantitative spect. Physics in Medicine and Biology, 39(3):509, 1994. [99] Benjamin M. W. Tsui, Grant T. Gullberg, Eric R. Edgerton, J. Glen Ballard, J. Randolph Perry, William H. McCartney, and Jan Berg. Correction of Nonuniform Attenuation in Cardiac SPECT Imaging. J Nucl Med, 30(4):497–507, 1989. [100] B.M.W. Tsui, H.-B. Hu, D.R. Gilland, and G.T. Gullberg. Implementation of simultaneous attenuation and detector response correction in spect. Nuclear Science, IEEE Transactions on, 35(1):778–783, feb. 1988. [101] G Van Soest, H Shemesh, M-K Wu, LWM van der Sluis, and PR Wesselink. Optical coherence tomography for endodontic imaging. In Biomedical Optics (BiOS) 2008, pages 68430F–68430F. International Society for Optics and Photonics, 2008. [102] P. Vemuri, Kholmovski E.G., Parker D.L., and Chapman B.E. Coil sensitivity estimation for optimal snr reconstruction and intensity inhomogeneity correction in phased array mr imaging. Information processing in medical imaging : proceedings of the ... conference, 19:603–614, 2005. [103] P Verma and RM Love. A micro ct study of the mesiobuccal root canal morphology of the maxillary first molar tooth. International Endodontic Journal, 44(3):210–217, 2011. [104] C. Vetter and R. Westermann. Spect reconstruction on the gpu. In Society of PhotoOptical Instrumentation Engineers (SPIE) Conference Series, volume 6913 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, apr 2008. [105] E.A. Vokurka, Watson N.A., Watson Y., Thacker N.A., and Jackson A. Improved high resolution mr imaging for surface coils using automated intensity non-uniformity correction: Feasibility study in the orbit. Journal of Magnetic Resonance Imaging, 14(5):540–546, 2001. [106] Uroš Vovk, Franjo Pernuš, and Boštjan Likar. Mri intensity inhomogeneity correction by combining intensity and spatial information. Physics in Medicine and Biology, 49(17):4119–4133, 2004. [107] Uroš Vovk, Franjo Pernuš, and Boštjan Likar. A review of methods for correction of intensity inhomogeneity in mri. Medical Imaging, IEEE Transactions on, 26(3):405– 421, 2007.


xi

[108] W.M. Wells, W.E.L. Grimson, R. Kikinis, and F.A. Jolesz. Adaptive segmentationof mri data. Medical Imaging, IEEE transaction on, 15(4):429–442, August 1996. [109] Miles N Wernick and John N Aarsvold. Emission tomography: the fundamentals of PET and SPECT. Academic Press, 2004. [110] Brita Willershausen, Adrian Kasaj, Bernd Röhrig, and Benjamin Briseño Marroquin. Radiographic investigation of frequency and location of root canal curvatures in human mandibular anterior incisors in vitro. Journal of Endodontics, 34(2):152–156, 2008. [111] Enhua Wu and Youquan Liu. Emerging technology about gpgpu. In Circuits and Systems, 2008. APCCAS 2008. IEEE Asia Pacific Conference on, pages 618–622. IEEE, 2008. [112] Fu-Che Wu, Wan-Chun Ma, Rung-Huei Liang, Bing-Yu Chen, and Ming Ouhyoung. Domain connected graph: the skeleton of a closed 3d shape for animation. The Visual Computer, 22(2):117–135, 2006. [113] Weishi Xia, R.M. Lewitt, and P.R. Edholm. Fourier correction for spatially variant collimator blurring in spect. Medical Imaging, IEEE Transactions on, 14(1):100–115, mar 1995. [114] Fang Xu and K. Mueller. Accelerating popular tomographic reconstruction algorithms on commodity pc graphics hardware. Nuclear Science, IEEE Transactions on, 52(3):654–663, june 2005. [115] Fang Xu and Klaus Mueller. Gpu-acceleration of attenuation and scattering compensation in emission computed tomography. na, jul. 2007. [116] Fang Xu and Klaus Mueller. Real-time 3d computed tomographic reconstruction using commodity graphics hardware. Physics in Medicine and Biology, 52(12):3405, 2007. [117] Fang Xu, Wei Xu, Mel Jones, Bettina Keszthelyi, John Sedat, David Agard, and Klaus Mueller. On the efficiency of iterative ordered subset reconstruction algorithms for acceleration on gpus. Comput Methods Programs Biomed, 98(3):261–270, 2010. [118] Masashi Yamada, Yoshinobu Ide, Satoru Matsunaga, Hiroshi Kato, and K Nakagawa. Three-dimensional analysis of mesiobuccal root canal of japanese maxillary first molar using micro-ct. The Bulletin of Tokyo Dental College, 52(2):77, 2011.

dc_592_12 xii [119] Habib Zaidi and Bruce Hasegawa. Determination of the Attenuation Map in Emission Tomography. J Nucl Med, 44(2):291–315, 2003. [120] G.L. Zeng, Chuanyong Bai, and G.T. Gullberg. A projector/backprojector with sliceto-slice blurring for efficient three-dimensional scatter modeling. Medical Imaging, IEEE Transactions on, 18(8):722–732, aug. 1999. [121] G.L. Zeng, G.T. Gullberg, B.M.W. Tsui, and J.A. Terry. Three-dimensional iterative reconstruction algorithms with attenuation and geometric point response correction. Nuclear Science, IEEE Transactions on, 38(2):693–702, apr 1991. [122] Y. Brady M. Smith S. Zhang. Segmentation of brain mr images through a hidden markov random field model and the expectation-maximization algorithm. IEEE Transactions on Medical Imaging, 20(1):45–57, 2001.

MTA DOKTORA PÁLYÁZAT DOKTORI ÉRTEKEZÉS KÉPALKOTÓ ÉS KÉPFELDOLGOZÓ ELJÁRÁSOK HATÉKONYSÁGÁNAK NÖVELÉSE AZ ORVOSINFORMATIKÁBAN

Recommend Documents