Rekonstrukció

Tankönyv Fizikusoknak » A röntgendiagnosztika alapjai » Kúpsugaras röntgentomográfia ( Cone-beam CT ) » Rekonstrukció

A kúpsugaras CT nyersadatainak rekonstruálására jelenleg a szűrt visszavetítéses algoritmus (Filtered Backprojection – FBP) a legelterjedtebb módszer. Az FBP legfontosabb jellemzője, hogy ún. transzformációs vagy analitikai eljárás: azon a feltételezésen alapul, hogy a szkennelés során begyűjtött vetületi vagy árnyékképek a rekonstruálni kívánt lineáris gyengítési együtthatók Radon transzformációja. A Radon transzformáció analitikai inverze direkt megoldását adja a rekonstrukciós problémának. Ha a rendelkezésre álló projekciók ideálisak, akkor az FBP rekonstrukció azonos az ideális rekonstrukcióval. A tökéletes rekonstrukcióhoz a következő feltételeknek kell teljesülni :

a képalkotást végző nyalábnak monokromatikusnak kell lennie
szórt fotonok nem, csak ún. primer fotonok vesznek részt a vetületi képek létrejöttében
a mérés során az objektum egésze és részei sem mozdulnak el
a projekciók gyűjtése során a vizsgált objektum egyik részlete sem lóg ki a képmezőből (Field of View –FOV)
végtelen számú projekció áll rendelkezésre
a mért jelnek zajmentesnek kell lennie
a mért jelnek egyenesen arányosnak kell lennie a detektált fotonok számával (linearitás), a detektor kép egyetlen fotont és végtelen fotont is képes érzékelni (dinamikus tartomány)
a detektorelemeknek tulajdonságaikban azonosnak kell lenniük

A Feldkamp-Davis-Kress -féle rekonstrukció (FDK) elmélete monokromatikus röntgen nyalábra épül, aminek a homogén közegben végbemenő intenzitáscsökkenését a Lambert-Beer törvény írja le:

$y_i=I_0 \cdot e^{- \sum_{j=1}^{N}(l_{ij}\mu_j)}$

amiről a korábbi fejezetekben már részletesebben is írtunk. Az FDK típusú rekonstrukció alapjairól részletek a „Practical cone-beam algorithm” cikkben találhatók.

Vetületi képek normalizálása

Az FBP első alaplépése a nyers vetületi kép pixeleinek normalizálása:

$P'_{\theta, i}=-ln(\frac{I_{\theta,i}}{I_0})$

ahol $I_{\theta,i}$ a $\theta$ szöghöz tartozó vetületi kép (projekció) i. pixelében mért intenzitás, I_0 a forrásból kijövő intenzitás. Az I_0 -at egy kalibrációs lépés során minden energiára, és expozíciós időre lemérjük, úgy hogy a röntgenforrás és a detektor között semmilyen tárgy sem található.

Képek szűrése

Az FBP következő lépése a vetületi képek frekvencia tartománybeli szűrése. Hagyományosan az ún. RAMP, vagy más néven RamLak szűrőt alkalmaznak. Ez eliminálja a kép alacsonyfrekvenciás összetevőit, de kiemeli a nagyfrekvenciásokat. Gyakorlatban ez a szűrő emeli a kép zajtartalmát, mivel a diszkrét inverz Fourier transzformáció rosszul kezeli a jel nagy változásait, ugyanakkor a kép sokszor nem tartalmaz a rekonstrukció számára hasznos információkat a magas frekvenciás összetevőiben. Ezért ezt a szűrőt érdemes kombinálni ún. ablakozó szűrőkkel, amik a nagyfrekvenciás komponensek elnyomására alakítottak ki. Az így kapott szűrőket az ablakozó függvény típusával szokás elnevezni (8. ábra). Minél erősebb szűrést használunk a magas frekvenciákra, annál homogénebb, jobb jel/zaj arányú képeket kapunk, de ezzel együtt sajnos a rekonstruált képek térbeli felontását is csökkentjük (9. ábra), ezért mindig a megfelelő szűrőt kell kiválasztani. Ezeket a szűröket alapvetően három csoportba lehet kategorizálni.
Nagyfelbontású szűrők:

RamLak: $H(\omega)=\left | \omega \right |$
Butterworth: $H(\omega)=\left | \omega \right | \cdot \sqrt{\frac{1}{1+\omega^{2n}}}$
Shepp-Logan: $H(\omega)=\left | \omega \right | \cdot \left | \frac{sin \omega}{\pi \omega} \right |$

A legtöbb esetben javasolt a Shepp-Logan használata, mert ez 10%-kal csökkenti a képen a zajt, de a felbontást csupán 4%-kal rontja, a RamLak-hoz képest.
A második csoportba kerülhetnek:

Cosine: $H(\omega)=\left | \omega \right | \cdot cos\left ( \frac{\pi\omega}{\omega_{max}-1}-\frac{\pi}{2} \right )$
Hamming: $H(\omega)=\left | \omega \right | \cdot \left ( 0.54+0.46 \cdot \frac{2\pi\omega}{\omega_{max}-1} \right )$
Hann: $H(\omega)=\left | \omega \right | \cdot 0.5 \cdot \left ( 1 + \frac{2\pi\omega}{\omega_{max}-1} \right )$

Ezek már jelentős mértékben csökkentik a zajt, ugyanakkor a felbontást is csökkentik. Ha fontosabb a jó jel/zaj arány, mint a felbontás, érdemes ezeket használni.

Több ún. kisfelbontású szűrő is létezik, de ezek már jelentősen rontják a felbontást, ugyanakkor a BlackMan speciális esetekben használhatónak bizonyult (pl. különböző kalibrációs protokollokban).

8. ábra. Szűrők karakterisztikái. Vízszintes tengelyen a frekvencia, függőleges tengelyen a hozzá tartozó erősítés tartozik, mindkét esetben a [0..1] tartományba normalizálva

9.ábra, A felbontás és zaj kapcsolata különböző szűrők esetén. Vízszintes tengelyen a zaj mértéke látható (homogén víz Hounsfield értékének szórása- a kisebb érték a jobb), függőleges tengelyen a felbontás olvasható le (MTF- a nagyobb érték a jobb)

A szűrők finomhangolásánál lehetőség van ún. vágási frekvencia (cutoff) megadására is. Ezt általában százalékban lehet megadni, a megadott érték feletti komponenseket a szűrő ki fogja nullázni. Egy 50%-os cutoff érték tehát azt jelenti, hogy a frekvenciakomponensek felső felét szeretnénk kiszűrni, ami természetesen jót tesz a jel/zaj aránynak, de a felbontást pontosan a felére fogja csökkenteni.

Visszavetítés

Az FBP során a szűrt képeket vissza kell vetíteni a térfogatba (backprojection). A visszavetítő operátorok típusai: pixel alapú (pixel-driven ), távolság alapú (distance-driven ) és voxel alapú (voxel driven ).
Pixel alapú visszavetítés során a pixel középpontjából indulunk el a fókuszpont felé. A sugár mentén minden egyes térelemre (voxel) kiszámítjuk, hogy az adott pixelben kiolvasható érték mennyire járul hozzá az adott voxelhez.

10. ábra. Pixel alapú visszavetítés

A pixel értékek már normalizáltak ( $P'_{\theta,i}=-ln(\frac{I_{\theta,i}}{I_0})$ ), és a szűrés utáni értéket jelölje: $P_{\theta,i}$ . Az inkrementális voxel érték egy $\theta$ szöghöz tartozó projekcióra a következő képlettel számítható:

$V_{\theta,j}=\sum_{i=0}^{N} (P_{\theta,i} \cdot \frac{V_i}{V_\nu})$

ahol V_i a voxel és sugár-piramis metszetének térfogata, $V_\nu$ pedig a voxel térfogata, i minden esetben a detektor egy pixelét indexeli. A sugárpiramist úgy kell érteni, hogy a detektor egy pixel négyzet alakú, és ezt kötjük össze a röntgencső fókuszpontjával. A 10. ábran kék és sárga területek jelzik a metsző területeket. Ha a fenti számítás az összes projekciós képre kiszámítjuk, akkor ezek összege megadja a rekonstruált voxel értékét. Az ismertetett képlet, bár fizikailag pontos rendkívül számításigényes, ezért sok esetben csak vonal menti integrálokat számítanak ki, vagyis minden pixelből egy sugarat indítanak, és ez esetben csak kocka (vagy általánosabb esetben téglatest), és vonal metszését kell számítani. A pixel alapú megközelítésnek másik hátránya, hogy párhuzamos architektúrákon nem hatékony, mert a kimenő adatokon (voxelek) írási ütközés lehet a memóriában (hiszen több sugár is metszheti ugyanazt a voxelt), ami az állandó szinkronizációs műveletek miatt lassabb.

A távolság alapú visszavetítés egy viszonylag új megközelítés. A módszer során a voxelek és a detektorpixelek vízszintes és függőleges határait egy közös tengelyre vetítjük a forgástengely és arra merőleges síkjában (11. ábra). Ezután kiszámítjuk normalizált átlapolásokat mindkét síkban.

11. ábra, Távolság alapú visszavetítés

Ekkor a voxelek járulékának kiszámítása a pixelekhez:

$p=\underbrace{\frac{t}{cos\alpha cos\gamma}\frac{O_1}{w_1}\frac{O_2}{w_2}}_{l_{ij}}\mu$

Ahol a $\frac{t}{cos\alpha cos\gamma}$ a forrás és a detektorpixel által meghatározott sugár két, a detektorsíkkal párhuzamos voxelsík között áthaladó hossza. Mivel ezzel a módszerrel az előre és visszavetítés is hasonlóan megvalósítható, elsősorban iteratív rekonstrukciók esetén használják.

Voxel alapú megközelítésnél minden voxel-re függetlenül számítjuk ki a voxel közepén átmenő, fókuszból induló sugár metszéspontját a detektorral. Ez az elrendezés igen kedvező a párhuzamos architektúráknak, így a grafikus processzoroknak is (GPU). A bonyolultabb része az algoritmusnak az úgynevezett footprint (lenyomat) függvény kiszámítása minden voxelre. A lenyomatot úgy kell elképzelni, hogy a voxelt, mint kockát levetítjük a detektor síkjára, a függvény értéke annak a vonalnak a hossza lesz (12. ábra), amekkora távolságot megtesz a röntgensugár az adott pontban:

12. ábra. Voxel alapú visszavetítés

A $\theta$ szöggel jellemzett projekcióhoz tartozó inkrementális voxel érték meghatározható:

$V_{\theta,j}=\frac{\int (P_{\theta}(t) \cdot f_{\theta,j}(t)) dt}{\int f_{\theta,j}(t) dt}$

képlettel, ahol j a voxel indexe, t a detektor egy pixele, $P_\theta (t)$ függvény a logaritmikusan normált projekció értéke a t pontban, $f_{\theta,j} (t)$ a „Footprint” függvény. A voxel értéket megkapjuk, ha az összes projekcióhoz tartozó inkrementális értéket összeadjuk. A fenti integrál kiszámítása szintén számításigényes, ezért általában ezt is közelíteni szokták. Egy jó közelítési megoldás az ún. Monte-Carlo integrálás. Ennek a lényege, hogy az integrálszámítást a voxelen belüli egyenletes (kvázi-véletlen) eloszlású pontokban vett mintavételezéssel közelítjük (13. ábra). Ekkor a voxel értéke a következőképpen számítható:

$V_{\theta,j}=\frac{\sum_{k=0}^{S}P_{\theta,k}}{S}$

ahol S a mintapontok száma, k a mintapont indexe. Általánosságban elmondható, hogy minél több mintát használunk, annál jobban közelítjük az eredeti integrálszámítást. Ennek egyik lehetséges szélsőértéke lehet az is, ha csak egy mintapontunk van, sokszor használják ezt is, mivel így nagyon gyors algoritmushoz jutunk. A többszörös mintavételezéssel jelentősen csökkenthetjük a zajt, hogy közben a felbontás nem csökken .

13. ábra. Többszörös mintavételezés a voxelen belül

Hounsfield korrekció

A rekonstrukció befejező lépése a kiszámított voxel értékek normalizációja, a Hounsfield korrekció. A CT-készülék ugyanis a lineáris sugárgyengítési együtthatót méri, de a könnyebb kezelhetőség kedvéért a vízre vonatkoztatott relatív sugárgyengítést jelzi ki Hounsfield egységekben (rövidítése: angolul HU, magyarul HE).

$HE=1000 \cdot \frac{\mu-\mu_{viz}}{\mu_{viz}-\mu_{levego}}$

ahol
µ a tetszőleges anyag elnyelése
µ_víz a víz elnyelése
µ_levegő a levegő elnyelése

Hardver architektúrák

Bár a szűrt visszavetítés az egyik leggyorsabb rekonstrukciós algoritmus, a minőség és a felbontás növelésével ez esetben is nagyon nagy számítási kapacitásra, vagy hosszú rekonstrukciós időre van szükség. Egy manapság átlagosnak mondható felbontású és méretű térfogat (pl.: 1024³ voxel) előállításához egy modern, több magos CPU-n is órákra lenne szükség. Korábban több számítógépből álló fürtök segítségével sikerült az időt leszorítani, de ezek a rendszerek igen drágák voltak, nagyon sok áramot fogyasztottak, és nagyon sok helyet foglaltak. Alternatívaként célhardverekkel oldották meg a rekonstrukciót, pl.: FPGA alapú kártyák használatával. Ennek a módszernek a hátránya, hogy a célhardverek fejlesztési ideje nagy és költséges, valamint nehezen fejleszthetők tovább, ha új algoritmus beépítésére lenne szükség. Szerencsére a modern számítógépek elengedhetetlen összetevője a grafikus kártya (GPU), ami eredetileg játékok és 3D alkalmazások (pl. CAD szoftverek) futtatására lett kitalálva. A GPU-k a fejlődésük során néhány éve eljutottak egy olyan fázisba, amikor is alkalmassá váltak általános algoritmusok futtatására. Ekkor indult meg a GPU-k komoly használata igen számításigényes algoritmusok gyorsítására . Az olajipar, időjárás előrejelzés, pénzügyi alkalmazások és még számos egyéb terület után az orvosi képfeldolgozás is a GPU-k felé fordult. Egy modern GPU (Nvidia GTX-580) számítási kapacitása ~1500 GFLOPS (1500 milliárd lebegőpontos művelet másodpercenként), szemben az egyik jelenleg leggyorsabb processzor (Intel 980X) 65 GFLOPS teljesítményével. Ez a gyakorlatban azt jelenti, hogy komplett számítógép fürtök helyettesíthetők egyetlen közönséges PC-vel.

Iteratív rekonstrukciós módszerek

Mint már korábban megfogalmaztuk, a műtermékek a rekonstrukciós algoritmus által használt matematikai modell és a tényleges mérőrendszer közötti eltérések következtében lépnek fel. Ezért pontosabb matematikai modellt alkalmazó algoritmus használatával csökkenthetők vagy megszüntethetők a műtermékek.
Az iteratív algoritmusoknak két nagy csoportja vani: az ún. algebrai (Algebraic Iterative Reconstruction –ART) és a statisztikai (Statistical Iterative Reconstruction– SR) algoritmusok (14. ábra). Mindkettőről elmondható, hogy jelentősen számolásigényesebb, így a re-konstrukció hosszabb ideig tart, mint a FBP algoritmus esetén, azonban jobb minőségű rekonstrukciót tesznek lehetővé és kevésbé érzékenyek a mérőrendszer hibaforrásaira, pl. a zajra vagy a vetületi nyers adathalmaz hiányosságaira.

14. ábra. A rekonstrukciós algoritmusok típusai

A rekonstrukciós probléma algebrai megközelítése esetén a fizikai folyamatra felírható egyenletrendszerek megoldása a cél. Azonban ennél a módszernél az egyenletrendszerben az ismeretlenek száma megegyezik a rekonstruálandó térfogat voxeleinek számával, ami adott esetben több milliárd is lehet, amit csak iteratív módszerekkel lehet megoldani. Ezt a módszert ezért elsősorban a PET és SPECT készülékeknél alkalmazzák, ahol lényegesen kisebb adathalmazt kell rekonstruálni.

Az iteratív algoritmusok másik nagy csoportja, melyről a legtöbb publikáció olvasható, a statisztikai iteratív algoritmusok. Az alkalmazott matematikai modell figyelembe veszi a zaj statisztikai természetét a vetületi képekben, ami által a rekonstruált képekben csökken a zaj. Az algoritmus lényegéből adódóan kevésbé érzékeny a vetületi képekben lévő hiányosságokra, mint az analitikai inverzen alapuló módszerek. Mivel az inverz transzformációnak nem szükséges meghatározni az explicit formáját, a képalkotásban alkalmazott geometriák széles skálája esetén használhatók. Ezen felül a statisztikai rekonstrukciós algoritmusokkal csökkenthetők vagy megszüntethetők a foton szórásból, a sugárkeményedésből, a fém jelenlétéből adódó műtermékek, mivel beépíthető az algoritmusba a röntgenfoton–anyag kölcsönhatás precíz modellje. Végül a fenti tulajdonságok következtében kisebb dózis mellett jobb képminőséget eredményeznek.

Az egyetlen valódi hátránya ma még a 2-3 nagyságrenddel nagyobb számolásigény, amivel arányosan nő a rekonstrukcióhoz szükséges idő is. Azonban a rekonstrukció felgyorsítására több irányba is folynak kutatások. Módosított algoritmusokat fejlesztenek ki, az algoritmusokat speciális hardver kiegészítőkre implementálják (pl. GPU) illetve több számítógép hálózatba kapcsolásával párhuzamos számítást valósítanak meg.
Az iteratív rekonstrukciós algoritmusok általános sémája a következő (15. ábra):

15. ábra. Az iteratív rekonstrukció általános sémája

Az iterációt megelőző lépésben egy nulladik, kiindulási CT-képet kell generálni. Ez lehet egy üres mátrix, de fel lehet használni az FBP által kiszámolt képet is. Az aktuális (kezdetben az előbb említett nulladik) CT-képből a CT-szkennelés modelljét felhasználva kiszámítja az algoritmus a CT-képhez tartozó vetületi képeket. Az így kapott vetületi képeket összehasonlítja a mért, valós vetületi képekkel. Ez annak függvényében fog egyezni vagy eltérni, hogy mennyire közelíti meg a CT-kép a valóságot. Az összehasonlítás után kiszámítja a vetületi képek hibáját, amiből egy transzformáció segítségével előállítható a CT-kép hibája. Ezt ismerve pedig korrigálni lehet az előző lépés CT-képét, és újra indul a ciklus.

Export

Medical Imaging

Site Language: English