A kutatás célkitűzései
A számítógépes grafika területén egyrészt a színterek egyre nagyobb geometriai komplexitásának hatékony kezelése miatt van szükség párhuzamos feldolgozásra, másrészt a különféle komplex optikai jelenségek (például globális illuminációs effektusok) fizikailag korrekt számítása is csak párhuzamos architektúrán implementálható valós időben. Az orvosi vizualizációs és képfeldolgozó alkalmazásokban elsősorban a nagyfelbontású 3D adatok (például CT, MRI, PET) mérete teszi szükségessé a feldolgozási feladatok adatszintű partícionálását. Moderált tárigényű színtereket illetve adathalmazokat megfelelő hatékonysággal lehet kezelni egyetlen grafikus feldolgozó egységen is, a komplexebb modellek feldolgozásának viszont gátat szab a GPU limitált lokális memóriája. Az egyik lehetőség az adatpartíciók operatív memória és lokális GPU-memória közötti mozgatása. Ekkor azonban a GPU számítási kapacitása nem használható ki teljesen, mivel a sávszélesség jelenti a szűk keresztmetszetet. A másik megoldás a GPU klaszterek alkalmazása, ahol egy-egy GPU-t rendelünk az egyes partíciókhoz, így nem kell nagy mennyiségű adatot dinamikusan mozgatni. Az egyes GPU csomópontok (node) a partíciók feldolgozásának eredményét, mint 2D képi információt továbbítják egymásnak. Ilyenkor a fő cél a kommunikációs overhead minimalizálása.
Annak ellenére, hogy a konvencionális GPU-t eredetileg hagyományos számítógépes grafikai alkalmazásokra fejlesztették ki, ma már egyre inkább használják általános célú társprocesszorként, például jel- (Fourier-transzformáció, diszkrét jelek folytonos rekonstrukciója) és képfeldolgozási (újramintavételezés, tömörítés, szegmentálás, regisztráció, tomográfiás rekonstrukció) feladatokra, fizikai szimulációra (például folyadékszimuláció vagy többszörös szóródás modellezése inhomogén anyagokban), vagy lineáris egyenletrendszerek illetve differenciálegyenletek megoldására. A GPU implementáció lehetőségét minden olyan alkalmazás esetében célszerű megvizsgálni, amely részben vagy kizárólag SIMD (Single Instruction Multiple Data) típusú processzálást végez. Abban az esetben, ha dinamikus adatszerkezetek kezelésére is szükség van, a szuperszámítógép CPU alapú számítási kapacitására lehet építeni, de ekkor is érdemes a feldolgozási feladatok CPU és GPU közötti ésszerű megosztását elemezni.
A szuperszámítógép GPU kapacitását főként 3D orvosi vizualizációs és képfeldolgozási feladatokra szeretnénk használni. Hosszabb távon a cél egy olyan keretrendszer felépítése, amely a nagy tárigényű CT és MRI adatok tárolását, feldolgozását, és vizualizációját a szerveren valósítja meg, kliens oldalról pedig csak a generált képeket kell letölteni. A kliens gépektől nem várunk el nagy számítási kapacitást, a lokális GPU-t legfeljebb a tömörített formában továbbított képek valós idejű kicsomagolására használnánk. Ennek a megoldásnak az a legfőbb előnye, hogy nem kell a kliens oldalon nagy mennyiségű 3D adatot letölteni és tárolni, így a rendszer kellőképpen gyors hálózati összeköttetéssel akár egy olyan laptopról is használható lenne, amelybe nincs nagy teljesítményű GPU integrálva.
A szerver oldalon azt érdemes megvizsgálni, hogy a képfeldolgozó és vizualizációs műveleteket miként lehet a leghatékonyabban megosztani a CPU és GPU csomópontok között. A CPU implementáció rugalmasabb és általában nagyobb numerikus precizitást biztosít. A GPU SIMD-jellegű architektúrája viszont jobban kihasználható konvolúció-szerű műveletek gyors végrehajtására. GPU implementáció esetén számolni kell az adatfeltöltés illetve kiolvasás járulékos időköltségével. Ezt a többletköltséget csak akkor lehet kompenzálni, ha a feldolgozás a GPU-n legalább egy nagyságrenddel gyorsabb mint a sokmagos CPU-n. Kutatási szempontból azt is érdemes elemezni, hogy egy tisztán CPU implementáció megfelelő optimalizálással mennyire marad el teljesítményben egy tisztán GPU implementációtól, figyelembe véve az adatmozgatási többletköltséget.
Első projekt előrehaladási jelentés (2011.09 - 2011.12)
Hatékony rekonstrukciós szűrők tervezése BCC rácson mintavételezett térfogati adatokhoz
Ha egy izotrópikusan sávkorlátozott 3D jelet a hagyományos kockarács (Cartesian Cubic lattice – CC rács) helyett egy optimális tércentrált kockarácson (Body-Centered Cubic lattice – BCC rács) mintavételezünk, akkor egységnyi térfogatból közel 30%-kal kevesebb mintát kell venni, ahhoz hogy az eredeti jelet tökéletesen rekonstruálni lehessen. A BCC rács ugyanakkor nem szeparálható, ezért az 1D rekonstrukciós szűrők előnyős tulajdonságai (interpoláló vagy kvázi-interpoláló jelleg, approximációs rend) általában nem öröklődnek a BCC rácson, ha a hagyományos szeparálható vagy gömbszimmetrikus kiterjesztéseket alkalmazzuk. A szeparálható B-spline szűrőket ezért egy olyan nem szeparálható koszinuszos ablakkal moduláltuk, amely jól illeszkedik a BCC rács geometriájára. A kapott rekonstrukciós szűrők numerikusan, és vizuálisan is jobb eredményt adnak, mint a korábban javasolt box spline szűrők. Elméleti szempontból fontos eredmény, hogy a szűrők energiája a lezáró sávban az aliasing spektrumok közötti holttérre koncentrálódik, ezért a postaliasing hatás numerikusan és vizuálisan is minimális. Egy interaktív térfogat-vizualizációs alkalmazáson demonstráltuk, hogy a koszinuszosan modulált B-spline szűrők GPU implementációja sokkal hatékonyabb, mint a box spline szűrőké.
Megvilágítástól függő opacitásmoduláción alapuló térfogat-vizualizáció
3D orvosi adatok (CT, MR) interaktív megjelenítéséhez kidolgoztunk egy új térfogat-vizualizációs modellt, amely a voxelek opacitását a megvilágítás irányától függően definiálja. Ha N szintfelületet szeretnénk megjeleníteni, akkor N különböző frekvenciájú és irányú fényforrást alkalmazunk. Minden egyes szintfelület csak egy adott frekvenciájú fényre érzékeny, tehát csak a megfelelő fényforrásból érkező fényt veri vissza. Az egyes felületi pontok opacitása akkor lesz maximális, ha a beeső fénysugár merőleges a felületre. Másrészt minél inkább érintőleges a beeső fénysugár, az opacitás annál kisebb lesz. Ennek a modellnek az alkalmazásával vizuálisan jobban elkülöníthetők a szintfelületek, mint a hagyományos térfogat-vizualizáció esetén, mivel a különböző szintfelületek hatása a képtérben jobban elkülönül. Ezáltal a belső struktúrák könnyebben értelmezhetők. Főként akkor kapunk részletgazdag képeket, ha a vetítősugarak mentén nagy mintavételi frekvenciával számolunk, ami jelentősen megnöveli a számítási költséget. A módszer egyetlen GPU alkalmazásával nem biztosít interaktív megjelenítést, ezért a későbbiekben szeretnénk GPU-klaszterre adaptálni.
3D orvosi adatok GPU-alapú szegmentálása
Zajos 3D orvosi adatok szegmentálására kidolgoztunk egy absztrakt távolság-transzformáción alapuló módszert. A felhasználó a szegmentálandó objektumot egy referenciaponttal jelöli ki. Ezután az algoritmus minden voxelhez kiszámolja annak az ösvénynek a költségét, melyen keresztül a referenciapontból a legkisebb költséggel lehet eljutni az adott voxelbe. Az ösvény mentén az összköltséget úgy definiáljuk, hogy vesszük az érintett voxeleket, és kiszámoljuk ezek sűrűségének átlagos négyzetes eltérését a referenciapont sűrűségétől. A voxeleket a hozzájuk rendelt költségek küszöbözésével tudjuk osztályozni. A módszer hasonló eredményt ad, mint a népszerű szinthalmazos (level set) szegmentálás, de GPU-n jóval hatékonyabban implementálható.
Iteratív tomográfiás rekonstrukció GPU implementációja
A CT képalkotás gyakorlatában a szűrt visszavetítésen alapuló rekonstrukció (Filtered Back Projection - FBP) a leggyakrabban használt módszer. Bár az FBP a leghatékonyabb tomográfiás rekonstrukciós algoritmus, képminőség szempontjából jóval rosszabb eredményt ad, mint a költségesebb iteratív módszerek (Algebrai Rekonstrukció – ART, Expectation Maximization - EM). Kihasználva a mai modern GPU-k számítási kapacitását az iteratív módszerek versenyképessé tehetők. Az FBP, ART és EM módszereket C++/CUDA környezetben implementáltuk és részletesen összehasonlítottuk képminőség, numerikus pontosság és hatékonyság szempontjából.
Második projekt előrehaladási jelentés (2012.01 - 2012.06)
A Marschner-Lobb tesztjel kalibrálása CC, BCC, és FCC rácsokra
A közismert Marschner-Lobb tesztjelet eredetileg hagyományos kockarácson használható folytonos rekonstrukciós szűrők vizuális összehasonlítására javasolták. Az utóbbi időben viszont alkalmazták BCC és FCC rácsokra tervezett rekonstrukciós sémák elemzésére is. Elméletben és gyakorlatban is igazoltuk, hogy ezek a korábban publikált összehasonlítások félrevezetőek, mivel a Marschner-Lobb tesztjel sajátos spektrumának köszönhetően erősen preferálja a szuboptimális FCC mintavételezést az optimális BCC mintavételezés helyett.
Az optimális 3D mintavételezés elmélete azon az előfeltételezésen alapul, hogy az eredeti jel sávkorlátozott és a spektruma jól körülhatárolható egy gömbfelülettel. A leggyakrabban használt Marschner-Lobb jel ugyanakkor nem teljesíti ezt az előfeltételt, mivel a spektruma inkább korong mint gömb alakú. Ennek következtében a spektrum legsűrűbb, átlapolódásmentes pakolását a frekvenciatartományban egy BCC rács biztosítja. Mivel az FCC rács a BCC rács duálisa a tértartományban, ezért az FCC rácson ritkábban lehet a tesztjelet mintavételezni mint a BCC rácson, amely gömbösen sávkorlátozott jelek mintavételezéséhez optimális. Ennek következtében a Marschner-Lobb jel nem alkalmas a BCC és az FCC mintavételezés korrekt összehasonlítására, mivel a korong alakú spektrumának köszönhetően eleve előnyben részesíti az FCC rácsot. A korábbi összehasonlító elemzések ezért hibásan vonják le azt a következtetést a Marschner-Lobb tesztek eredményeiből, hogy az FCC rács biztosítja általában a legkisebb prealiasing hatást.
Állításaink gyakorlati igazolásához, mindhárom rácson kiszámoltuk az ideális rekonstrukciós szűrésnek egy nagy pontosságú Lánczos szűrős közelítését, majd a mintavételi frekvencia fokozatos növelésével ekvivalens, prealiasingmentes reprezentációkat kerestünk. Ehhez szükség volt egy hatékony GPU implementációra, mivel CPU-n a számítás napokig tartott volna. Megmutattuk, hogy valóban az FCC rács biztosítja a Marschner-Lobb jel leghatékonyabb mintavételezését.
Sziluettek és belső élek hatékony renderelése geometria-processzálással
Olyan, nem-fotorealisztikus képalkotásban használható eljárást dolgoztunk ki, amely egy háromszögháló karakterisztikus körvonalait egy menetben, a GPU geometriai árnyalóegységét kihasználva jeleníti meg. A módszer sziluettek és a felület görbületét megszakító jellegzetes belső élek rajzolására is alkalmas. A megközelítés lényege, hogy a geometriai reprezentációt egy előfeldolgozási lépésben olyan információval egészítjük ki, amely minimális tárigény- és sávszélességigény-növekedés mellett, lehetővé teszi a geometriai árnyalóegység számára a rajzolandó élek azonosítását. Az előfeldolgozási lépés a szokásos háromszög-szomszédossági számításba beépíthető, így nem jelent többletköltséget. Hasonlóképpen, az élek megjelenítése közel azonos költségű mint a háromszögháló megjelenítése. A módszer elkerüli a képtérbeli élkereső módszerek diszkretizálásból adódó hibáit, így kiváló minőségű eredményt ad.
PET rekonstrukció optimalizálása a parciális térfogati effektus csökkentésével
Az úgynevezett parciális térfogati effektus (partial volume effect) csökkentésére egy olyan algoritmust dolgoztunk ki, amely morfológiai műveleteken alapul és képes azonosítani azokat a régiókat, ahol a PET rekonstrukció eredményében és egy megfelelően regisztrált és szegmentált CT reprezentációban is határfelület található. A korrekció ezekre a régiókra koncentrálódik. A módszer hatékony implementálásához kihasználtuk a GPU párhuzamos számítási kapacitását.
Harmadik projekt előrehaladási jelentés (2012.07 - 2012.12)
Anatómiai régiók automatikus detektálása
Egy olyan új eljárást dolgoztunk ki, amely MR adatokban képes automatikusan felismerni különböző anatómiai régiókat. Létrehoztunk egy referencia adatbázist, amely négy páciens teljes test MR felvételeit tartalmazza felcímkézve. Az egyes szeletek leíróvektorait Zernike-transzformációval generáltuk, amely alkalmasnak bizonyult a szeleten ábrázolt anatómiai struktúra geometriájának tömör és intuitív reprezentációjára.
Az elkészített adatbázist felhasználva egy neurális háló alapú tanulórendszert építettünk, amely képes a Zernike-transzformáltak alapján a szeletek durva osztályozására. Mivel ez önmagában nem adott kielégítő eredményt, így egy utófeldolgozási szakaszt is készítettünk a rendszerhez.
Az utófeldolgozás során a neurális háló által legmegbízhatóbban felismert szeleteket, mint referenciákat felhasználva, újracímkézzük a többi szeletet. Ehhez egy kölcsönös információ alapú összehasonlítást alkalmazunk. Az így elért eredmények pontossága 90% felett van.
Analitikus MIP grafikus hardveren
Angiográfiás CT/MR adatokról legegyszerűbben úgy lehet vetületi képet készíteni, hogy minden vetítő sugár mentén megkeressük a 3D sűrűségfüggvény maximumát (Maximum Intensity Projection - MIP), és ezt az értéket rendeljük az adott sugárhoz tartozó pixelhez. Megmutattuk, hogy a mai grafikus kártyákon MIP képeket analitikusan is ki lehet számolni. Módszerünk egyrészt kihasználja a sugarak inkrementális kiértékelését, másrészt a behuzalozott trilineáris interpolációt. A hatékony implementációnak köszönhetően a megjelenítés interaktív, és a felhasználónak nem kell külön mintavételi frekvenciát beállítani.
Átlagoló és Metropolis iteráció pozitron emissziós tomográfiában
Az iteratív pozitron emissziós tomográfiás algoritmusok a voxel tér és a LOR tér között végeznek projekciókat, amelyek matematikai értelemben nagydimenziós integrálok kiszámítását jelentik. A dimenzió nagyon magas, elméletileg végtelen is lehet, ha szóródási eseményeket is kompenzálni kell. Nagydimenziós numerikus integráláshoz a Monte Carlo módszer használható. Mivel a voxelek és a LOR-ok száma többszáz millió lehet, és a projekciók a rekonstruálandó objektumtól is függenek, az integrálokat nem lehet előreszámítani, hanem az iteratív megoldás minden lépésében újra meg kell határozni. Javasoltuk a maximum-likelihood elvű iteratív rekonstrukció egy módosítását annak érdekében, hogy a Monte Carlo módszer hibája csökkenthető legyen. Ha a Monte Carlo minták helye minden iterációban azonos, az iteráció egy torzított eredményhez konvergál. Ha viszont a mintahelyek független valószínűségi változók, akkor az iteráció vagy divergál, vagy pedig fluktuál. A célunk az iteráció pontosságának növelése volt a minták számának és így a számítási idő növelése nélkül. Először az iterációs folyamatot elemeztük, majd két módszerre tettünk javaslatot. Az átlagoló iteráció az előrevetítések eredményeit átlagolja. A Metropolis módszer pedig elutasítja azon előrevetítéseket, amelyek meghiusítanák a pontos iterációt. Demonstráltuk, hogy ezekkel a módszerekkel a számítási idő jelentősen csökkenthető. A javasolt algoritmusok beépültek a TeraTomo rendszerbe.
Negyedik projekt előrehaladási jelentés (2013.01 - 2013.06)
Koszinuszos B-spline rekonstrukció optimális BCC rácson
Optimális tércentrált kockarácson (body-centered cubic lattice – BCC rács) mintavételezett térfogati adatok folytonos rekonstrukciójához kidolgoztunk egy hatékony, GPU-ra optimalizált algoritmust. Mivel a hagyományos trilineáris interpolációt a GPU hardveresen támogatja, ezért a szűrőkernellel vett konvolúciót visszavezettük trilineáris minták súlyozott összegére. Kihasználtuk, hogy a BCC rács két átlapolódó kockarácsból áll, melyeken egymástól függetlenül lehet hagyományos B-spline rekonstrukciót végrehajtani. A kapott mintákat egy olyan koszinuszos súlyozó függvénnyel kombináltuk, amely a BCC rácson biztosítja az interpolációt, azaz mind az elsődleges mind pedig a másodlagos kockarács rácspontjaiban tökéletesen reprodukáljuk az eredeti BCC mintákat. Egy alapos összehasonlító elemzéssel megmutattuk, hogy BCC rácson a koszinuszos B-spline rekonstrukció nemcsak gyorsabban kiértékelhető, de nagyobb numerikus pontosságot és magasabb képminőséget is produkál.
Térfogati metamorfózis a Radon-transzformáltakon végzett eloszlás-interpolációval
Térfogati modellek folytonos metamorfózisához kidolgoztunk egy újszerű módszert, amely a Radon-transzformáltak alaktartó interpolációján alapul. Elsődleges célunk az volt, hogy minél több irányból számolt 1D vetületek folytonos átmenetét biztosítsuk. Ehhez először elvégeztük a Radon-transzformációt, majd a megfelelő vetületek, mint 1D sűrűségeloszlások alaktartó eloszlás-interpolációjával kiszámoltuk az átmeneti állapotok Radon-transzformáltját, végül kiértékeltük az inverz Radon-transzformációt a klasszikus szűrt visszavetítés (Filtered Back-Projection) algoritmussal. A módszer összes lépését GPU-n implementáltuk. A szuperszámítógépen történő futtatást a jelentős számítási költség indokolta.
Ötödik projekt előrehaladási jelentés (2013.07 - 2013.08)
Ebben az időszakban az anatómiai régiók automatikus detektálására kidolgozott módszerünket tökéletesítettük. Az elkészült folyóiratcikk bírálat alatt áll.
