Úvod Jednou z intenzivně užívaných metod v oblasti výpočetní neuroanatomie je morfometrie založená na voxelech (voxel-based morphometry VBM), která je stále objektem aktivního výzkumu a diskuzí [1], [2]. Princip této metody spočívá ve zkoumání anatomických skenů z MRI voxel po voxelu s cílem automaticky vyznačit v mozku oblasti se signifikantními rozdíly mezi skupinou pacientů a kontrolní skupinou dobrovolníků. VBM je zřetězením několika algoritmů pro zpracování obrazů, přičemž samotná pružná registrace zde hraje klíčovou roli. Hlavní myšlenkou pružné registrace je nalézt lokální síly, které zdeformují plovoucí obraz tak, aby se více podobal obrazu referenčnímu. Použité nelineární transformace bývají založeny buď na hladkých bázových funkcích [3], [4], [5] nebo vychází z fyzikálních interpretací, např. z mechaniky těles [6], [7], [8]. Zatímco v prvním případě se získávají hladké málorozměrné deformace, kterými je možno potlačit pouze globální rozdíly ve tvaru mozku mezi jednotlivými subjekty, cílem druhé skupiny metod je dosáhnout v registrovaných obrazech perfektního slícování. V původní VBM [9], [10] jsou obrazy afinní transformací převedeny do stereotaktického prostoru a následně jsou zdeformovány na referenční obraz pomocí málorozměrné parametrické transformace založené na bázových funkcích užívaných v diskrétní kosinové transformaci. Koeficienty transformace jsou hledány optimalizačním algoritmem, který minimalizuje sumu čtverců rozdílů mezi intenzitami v obrazech a zároveň maximalizuje hladkost dosažených deformací. Rozdíly v intenzitách obrazů jsou postiženy pouze jedním škálovacím parametrem, a tak je tato metoda vhodná pouze pro monomodální obrazy.
V tomto článku je navržen algoritmus pružné registrace vhodný pro multimodální obrazy. Jsou zde jednak vysvětleny použité metody a dále jsou prezentovány experimentální výsledky získané z pokusů při kvantitativním hodnocení tohoto algoritmu. Metody Zde navržená pružná registrace se provádí pomocí víceúrovňového lícování podobrazů, viz obr. 1. Plovoucí obraz N je v iteračním procesu deformován tak, aby lícoval s referenčním obrazem M. Výsledné vychýlení u je pro každý voxel interpolováno z lokálních translací v obrazu N s využitím radiálních bázových funkcí (RBF). Translace podobrazů reprezentující deformační síly f jsou nalezeny v lokálních registracích, ve kterých je maximalizována symetrická regionová podobnostní metrika. Předpokládá se, že plovoucí obraz N byl v předešlém kroku transformován do souřadného systému obrazu M pomocí lineární registrace.
 Obr. 1. Schéma pružné registrace (detaily jsou popsány v textu).
Regionové symetrické podobnostní metriky V navrženém algoritmu mohou být použity různé podobnostní metriky. Regionová podobnost se vypočítá zprůměrováním bodových podobností v podobrazu (dále regionu) [11]:
 (1)
kde SW označuje podobnostní metriku regionu W se středem v bodě w a s KW překrývajícími se voxely x, ve kterých jsou vyhodnocovány bodové podobnosti S. Bodová podobnostní metrika SMI odvozená z dobře známé globální podobnostní metriky zvané vzájemná informace je definována jako [11]:
 (2)
kde pMN označuje sdruženou hustotu pravděpodobnosti intenzit a pM, pN jsou marginální hustoty pravděpodobnosti intenzit v obrazech M a N. Jiná bodová podobnostní metrika SUH je navržena v [6]:
 (3)
kde SH je bodová podobnostní metrika odvozená z globální sdružené entropie obrazů. Všechny zde definované metriky závisejí na sdružené hustotě pravděpodobnosti, která bývá odhadována ze vzájemného histogramu obrazů, který ovšem není znám, dokud nejsou obrazy přesně slícované. Sdružená hustota pravděpodobnosti proto bývá odhadována z obrazů slícovaných pomocí předchozí lineární registrace. Takto je tomu i ve [12], kde je navržena pro nelineární registraci regionová metrika založená na podmíněných pravděpodobnostech. Její definice je zde přepsána do podoby bodové podobnostní metriky:
 (4)
která je tedy dána pravděpodobností závislosti mezi danou intenzitou m v referenčním obrazu M a intenzitou n plovoucího obrazu N. Hustoty podmíněných pravděpodobností jsou získány normalizací hodnot v každém řádku vzájemného histogramu rovnoběžném s osou intenzit plovoucího obrazu N. Poslední zkoumaná podobnostní metrika je zde odvozena z (2):
(5) V každé úrovni dělení obrazů jsou hledána posunutí podobrazů v plovoucím obrazu N pomocí optimalizačního algoritmu, který maximalizuje vybranou regionovou podobnostní metriku. S využitím myšlenky symetrické registrace [13] je zde navržena symetrická regionová podobnostní metrika, která se počítá jako součet dvou metrik: pro přímou registraci podobrazu v N vzhledem k referenčnímu obrazu M a pro směr opačný.
Optimalizační algoritmus je zde navržen s cílem vyhnout se lokálním optimům. Jedná se o kombinaci algoritmu rozsáhlého vyhledávání (z angl. extensive search) a algoritmu nejstrmějšího sestupu. V prvním kroku je prostor všech možných translací prohledán s relativně velkým krokem. P nejlepších bodů je pak určeno jako počáteční body pro algoritmus nejstrmějšího sestupu. Za globální maximum se vezme maximum z P lokálních maxim nalezených nestrmějším sestupem.
Víceúrovňová deformace
Po nalezení lokálních translací všech podobrazů je spočítáno vychýlení u interpolací s využitím RBF. Vychýlení je počítáno zvlášť pro každou z D souřadnicových os:
 (6)
kde uk(x) je vychýlení bodu x obrazové mřížky podél osy k, R je radiální bázová funkce vzdálenosti ||x-wi|| mezi bodem obrazové mřížky a středem podobrazu w i. Koeficienty αi se spočítají dosazením lokálních translací f do (6) a vyřešením vzniklého systému B lineárních rovnic – opět pro každou souřadnicovou osu zvášť. Jako funkce R je zde použito Wendlandovy funkce s kompaktním nosičem, která byla úspěšně použita v [4] pro pružnou registraci založenou na význačných bodech. Její výhodné matematické vlastnosti umožňují její prostorovou dilataci a kompresi, což je důležité pro zavedení víceúrovňové strategie. Pro každou úroveň dělení podobrazů je nastavena velikost podobrazu na polovinu velikosti z předcházející úrovně. Vychýlení jsou sčítána přes všechny úrovně, čímž dochází k postupnému zpřesňování deformace. Algoritmus lze urychlit vynecháním z lícování těch oblastí, které neobsahují žádné kontury nebo povrchy. Další dělení podobrazu je prováděno jen v případě, když alespoň v jednom jeho voxelu je normalizovaný gradient obrazové intensity větší než předdefinovaný práh, viz obr. 2.
 Obr. 2: Pěti-úrovňové adaptivní dělení podobrazů.
Tkáňové pravděpodobnostní mapy Výpočet vzájemného histogramu obrazů není jediným způsobem, jak lze odhadnout sdruženou hustotu pravděpodobnosti obrazů M a N. Provádí-li se registrace ve stereotaktickém prostoru, často jsou k dispozici tkáňové pravděpodobnostní mapy reprezentující užitečnou apriorní informaci. S využitím této informace je zde počítán ještě jeden odhad sdružené hustoty pravděpodobnosti a tento je dále kombinován s obvyklým odhadem založeným na vzájemném histogramu obrazů váhovacím parametrem λ :
 (7) Takto jsou ve výsledném odhadu pMN zvýrazněny ty intenzitní páry, které reprezentují hlavní mozkové tkáně.
Během registrace jsou odhady sdružené hustoty pravděpodobnosti přepočítávány v každé iteraci. Nelze přitom použít klasické inkrementování v histogramových kontejnerech, poněvadž plovoucí obraz nemá vzhledem k dosud vypočteným vychýlením pravidelnou mřížku. Je proto použito zobecněné interpolace částečných objemů (GPVE) [14], kdy histogramové kontejnery související se všemi voxely referenčního obrazu sousedícími s vychýleným voxelem plovoucího obrazu jsou navýšeny o hodnotu interpolační jádrové funkce, kterou je v tomto případě kubický B-spline. Tohoto postupu je zde využito pro výpočet vzájemného histogramu obrazů. Při výpočtu založeném na tkáňových pravděpodobnostních mapách se interpolace provádí mezi voxely plovoucího obrazu, a tak byl tento postup adaptován do interpolace nepravidelně rozložených dat založené na váhování vzdálenostmi [15]. Experimentální výsledky Kvalita registrace navrženým algoritmem s využitím různých podobnostních metrik byla vyhodnocena s využitím 2D obrazových dat získaných z databáze simulovaných obrazů (Simulated Brain Database SBD) [16]. Velikost transverzálních řezů v SBD obrazech je 181 × 217 pixelů s velikostí pixelu 1×1 mm. Obrazy byly rozšířeny na velikost 217×217 pixelů, poněvadž čtvercový rozměr je výhodnější pro víceúrovňové dělení. Byly vygenerovány umělé deformace, a to z náhodných translací umístěných do 10 % náhodně vybraných pixelů. Tato silová pole byla vyhlazena pomocí gaussovských filtrů s náhodně volenou směrodatnou odchylkou (σ=10±5 mm). Výsledná vychýlení byla aplikována na 20 intenzitních obrazů a k nim korespondující segmentované obrazy. Průměrná počáteční chyba překryvu segmentovaných částí v původních nezkreslených obrazech a obrazech po deformaci byla 41,0%. Deformace pak byly potlačeny pomocí zde navržené pružné registrace a byl vypočten průměrný pokles chyby překryvu Δe. Pro plovoucí obrazy byly požity obrazy váhované T1, T2 a PD (hustotou protonů). Tyto obrazy byly zarušeny šumem a artefaktem označovaným jako nelinearita přenosu obrazové informace (intensity nonuniformity). Jako referenční obraz posloužil T1 váhovaný obraz bez simulovaných artefaktů. Výsledky pro různé maximálně úrovně dekompozice jsou v tab. 1. Chyba překryvu se zmenšovala až do 5. úrovně dělení, kdy byla velikost podobrazů 7x7 pixelů. Při dalším dělení podobrazů sice vzrostla ještě globální vzájemná informace, ovšem slícování obrazů hodnocené chybou překryvu a vizuální kontrolou bylo stejné nebo horší. Proto je zde 5. úroveň považována za maximální úroveň dělení podobrazů pro tento registrační algoritmus.
Tab. 1. Kvantitativní hodnocení navrženého algoritmu pomocí poklesů průměrné chyby překryvu Δe. Jsou zde prezentovány výsledky získané při použití různých podobnostních metrik, pro různé páry intenzitních obrazů a pro různé úrovně dělení podobrazů.
| Obrazy | Δe [%] | | SPC | SUH | SMI | SPMI | | 1. úroveň | | T1-T1 | 5,2 | 1,6 | 1,7 | 1,5 | | T1-T2 | 5,3 | 1,3 | 1,4 | 2,8 | | T1-PD | 4,4 | 1,5 | 1,3 | 2,2 | | 2. úroveň | | T1-T1 | 12,6 | 8,3 | 8,4 | 8,5 | | T1-T2 | 12,3 | 4,7 | 4,8 | 7,4 | | T1-PD | 10,1 | 5,6 | 6,0 | 5,8 | | 3. úroveň | | T1-T1 | 22,6 | 17,9 | 17,9 | 19,5 | | T1-T2 | 21,2 | 9,8 | 9,8 | 16,4 | | T1-PD | 17,1 | 11,4 | 12,3 | 12,2 | | 4. úroveň | | T1-T1 | 27,2 | 24,8 | 24,9 | 26,6 | | T1-T2 | 25,3 | 16,0 | 16,2 | 23,3 | | T1-PD | 20,2 | 15,6 | 16,4 | 17,4 | | 5. úroveň | | T1-T1 | 27,7 | 26,8 | 27,1 | 28,5 | | T1-T2 | 25,1 | 19,4 | 19,9 | 24,9 | | T1-PD | 20,0 | 16,9 | 17,9 | 19,2 | Závěr V tomto článku byl navržen algoritmus pro málorozměrnou registraci MRI obrazů vhodnou pro pružnou registraci obrazů na digitální atlas mozku. V experimentu s umělými deformacemi byla studována kvalita registrace při použití čtyř různých podobnostních metrik. Kvalita registrace byla kvantitativně hodnocena pomocí poklesu chyby překryvu v segmentovaných obrazech. Podobnosti podobrazů byla měřeny s využitím vzájemného histogramu obrazů a dále také s využitím tkáňových pravděpodobnostních map. Chyba překryvu byla nejmenší při použití metrik SPC a SPMI. Při implementaci algoritmu bylo použito zobecněné interpolace částečných objemů, takže nebylo nutné během registračního procesu přepočítávat deformovaný plovoucí obraz. Navržený registrační algoritmus je vhodný pro automatickou morfometrii v obrazech označovanou jako voxel-based morphometry, protože přesnost registrace lze zde řídit nastavením maximální úrovně dělení. Takto mohou být potlačeny v obrazech pouze globální tvarové rozdíly, zatímco jemná anatomická variabilita, která je důležitá pro statistické parametrické testy, může být uchována. Navíc díky použití multimodálních podobnostních metrik je možné tímto algoritmem registrovat obrazy na libovolný atlas mozku bez nutnosti další transformace intenzit v registrovaných obrazech.
Poděkovaní
Tato práce byla podpořena granty 1ET2085120511 AV ČR, 102/04/0472 a 305/04/1385 GAČR a dále výzkumným záměrem Vysokého učení technického v Brně MSM 0021630513.
Literatura | [1] | Friston, K. J., Ashburner, J.: Generative and Recognition Models for Neuroanatomy. NeuroImage, vol. 23, pp. 21–24, 2004. | | [2] | Davatzikos, Ch.: Why Voxel-based Morphometric Analysis Should Be Used with Great Caution When Characterizing Group Ddifferences. NeuroImage, vol. 23, pp. 17–20, 2004. | | [3] | Pauchard, Y., Smith M. R., Mintchev M. P.: Modeling Susceptibility Difference Artifacts Produced by Metallic Implants in Magnetic Resonance Imaging with Point-Based Thin-Plate Spline Image Registration. In Proc. 26th Conf. IEEE EMBS, 2004, pp. 1766–1769. | | [4] | Fornefett, M., Rohr, K., Stiehl, H. S.: Radial Basis Functions with Compact Support for Elastic Registration of Medical Images. Image Vision Comput., vol. 19, pp. 87–96, 2001. | | [5] | Rohlfing, T., Maurer, C. R., Bluemke D. A., Michael J. A.: Volume-Preserving Nonrigid Registration of MR Breast Images Using Free-Form Deformation With an Incompresibility Constraint. IEEE Trans. on Medical Imaging, vol. 12, pp. 730–741, 2003. | | [6] | Rogelj, P., Kovačič, S., Gee, J. C.: Point Similarity Measures for Non-rigid Registration of Multi-modal Data. Computer Vision and Image Understanding, vol. 92, pp. 112–140, 2003. | | [7] | Tang, S., Jiang, T.: Fast Nonrigid Medical Image Registration by Fluid Model. In: Proc. 6th Asian Conference on Computer Vision, 2004. | | [8] | Christensen, G. E., He, J.: Consistent Nonlinear Elastic Image Registration. In: IEEE Proc. MMBIA, 2001, pp. 37–43. | | [9] | Ashburner, J., Friston, K. J.: Voxel-Based Morphometry – The Methods. NeuroImage, vol. 11, pp. 805–821, 2000. | | [10] | Ashburner, J., Friston, K. J.: Nonlinear Spatial Normalization Using Basis Functions. Human Brain Mapping, vol. 7, pp. 254–266, 1999. | | [11] | Rogelj, P., Kovačič, S.: Point Similarity Measure Based on Mutual Information. In Biomedical Image Registration: revised papers, 2003, pp.112–121. | | [12] | Maintz, J. B. A., Meijering, E. H. W., Viergever M. A.: General Multimodal Elastic Registration based on Mutual Information. In Medical Imaging: Image Processing (Proc. of SPIE, vol. 3338), 1998, pp. 144–154. | | [13] | Rogelj, P., Kovačič, S.: Symmetric Image Registration. In Medical Imaging: Image Processing (Proc. of SPIE, vol.5032), 2003, pp. 334–343. | | [14] | Chen, H., Varshney, P. K.: Mutual Information-Based CT-MR Brain Image Registration Using Generalized Partial Volume Point Histogram Estimation. IEEE Trans. on Medical Imaging, vol. 22, pp. 1111–1119, 2003. | | [15] | Franke, R. Nielson, G.: Smooth Interpolation of Large Sets of Scattered data. International Journal for Numerical Methods in Engineering, vol. 15, pp. 1691–1704, 1980. | | [16] | Collins, D.L., Zijdenbos, A.P., Kollokian, V., Sled, J.G., Kabani, N.J., Holmes, C. J., Evans A. C.: Design and Construction of a Realistic Digital Brain Phantom. IEEE Trans. on Medical Imaging, vol. 17, pp. 463–468, 1998. | |
in English |
Česky
English
