Tloušťka kortikální části kosti je důležitou metrikou nejen při diagnóze osteoporózy, kde se analyzují dlouhé kosti a těla obratlů. Její mapování má své využití například v cerebrálním kortexu, kde její změny mohou souviset například s Alzheimerovou cho-robou [19] nebo ve spodní a horní čelisti, kde se sleduje při plánování míst ukotvení ortodontických implantátů [30]. Jak již bylo zmíněno, stanovení tloušťky kortexu v diagnostických CT datech se standardním rozlišením je vzhledem k jeho relativně malým rozměrům velmi komplikované. Většina publikovaných metod je z tohoto dů-vodu vyvíjena pro sledování kortexu femuru, který je oproti vertebrálnímu kortexu několikrát tlustší.
Tímto problémem se zabýval již Prevrhal et al. [34] ve studii provedené v roce 1997, ve které byly testovány limity přesnosti měření tloušťky kortexu. Nasnímán byl dedikovaný fantom simulující variabilitu tloušťky kortexu tvořený tělesem s ekviva-lentem hustoty vody obaleným různě silnými vrstvami hliníku. Při analýze výsledků skenování na𝜇CT skeneru s izotropním rozlišením 0,7 mm bylo zjištěno, že sledová-ním intenzit na profilu vrstvy metodou prahování na polovině jasového maxima je možné dosáhnout chyby menší než 10 % tloušťky jen tehdy, pokud jsou rozměry větší nebo srovnatelné s rozlišením, konkrétně šiřkou bodové rozptylové funkce (Point spread function, PSF) skeneru (Obr 2.8).
Obr. 2.8: Závislost naměřené kortikální tloušťky a skutečné tloušťky [34].
V roce 2001 byla studie stejnými autory rozšířena na sledování vlivu planárního rozlišení a tloušťky řezu. Ačkoliv zahrnutí modelu, který tyto vlivy korigoval, do výpočtu tloušťky kortexu umožnilo zpřesnit měření například v oblasti krčku femuru, pro měření kortexu těla obratlů byla průměrná chyba měření příliš vysoká [35].
Typickým přístupem při měření tloušťky tenké vrstvy v lékařských tomografic-kých datech je detekce průchodu nulou druhé derivace podél intenzitního profilu kolmého na měřenou vrstvu [37]. Dostatečnou přesnost pro klinické účely však ta-kový přístup zajistí pouze v případě, že vrstva dosahuje tloušťky aspoň 1,5 násobku PSF skeneru (typicky srovnatelné s rozměrem voxelu).
Jedna z prvních metod umožňující měření pod touto hranicí pochází z roku 2007.
Streekstra et al. [41] využívá při měření tloušťky chrupavky v zápěstním kloubu znalost PSF skeneru změřené při snímání fantomu tvořeného tenkými hliníkovými dráty ve vodním prostředí. Prvním krokem metody je určit v oblasti zájmu po-mocí převažujícího směru gradientu směr intenzitního profilu kolmého na vrstvu chrupavky. Metoda předpokládá, že PSF skeneru má tvar izotropní trojrozměrné Gaussovy funkce (na profilu se tedy promítne jako jednorozměrná Gaussova funkce) a původní scéna má podobu směsi dokonalých schodovitých funkcí, které odpoví-dají rozhraní kortikální kost - chrupavka - měkká tkáň - chrupavka - kortikální kost. Protože obraz snímaný odpovídajícím skenerem odpovídá výsledku konvoluce původní scény s PSF skeneru, lze toto rozhraní modelovat následující analytickou rovnicí (3.1): V této rovnici𝜇𝑖(𝑖= 0..6) představuje koeficienty útlumu jednotlivých tkání po-dél profilu a 𝑔1(𝑥, 𝜎) je jednorozměrná Gaussova funkce s odchylkou 𝜎. Optimální registrací této funkce na naměřený intenzitní profil získáme kromě odhadu koefi-cientů útlumu jednotlivých materiálů také poziční hodnoty 𝑥𝑖, ze kterých můžeme snadno odečíst tloušťky jednotlivých vrstev.
Tento přístup založený na registraci modelu nedosahuje přesnosti původního pří-stupu detekce průchodu nulou druhé derivace u objektů s většími rozměry. Pro vrstvy menší než PSF skeneru však výrazně snižuje vychýlení odhadu směrem k vyšším hodnotám typické pro původní přístup, které jej v těchto intervalech činí ne-použitelným. Metoda vykazuje nevychýlený odhad až do tlouštěk vrstev o velikosti 50% PSF skeneru (Obr 2.9).
Treece et al. [44] v roce 2010 navázal na tuto metodu a aplikoval ji na problém měření tenké kortikální vrstvy femuru. Poprvé se jedná o metodu mapující kompletní povrch proximální části kosti. Tvar kosti je extrahován semi-automatickou metodou prostého prahování a narůstání oblasti. Výsledkem je obvykle hrubý tvar kosti, jehož kontura je následně ručně upravena a převedena na trojúhelníkovou síť, jejíž vrcholy slouží jako body analýzy tloušťky kortexu. Linie jasového profilu je v každém bodě
Obr. 2.9: Závislost naměřené tloušťky vrstvy a skutečné tloušťky. Měření metodou průchodu nulou druhé derivace () a metodou založenou na registraci modelu (∙) při PSF skeneru 0,49 mm [41].
extrahována jako normála na povrch trojúhelníkové sítě a její přesná manuální de-finice je tedy kritická pro následující automatickou analýzu. Jasový profil je poté metodou kubického splajnu nadvzorkován na 100 vzorků.
Model použitý při registraci na jasový profil je v zásadě stejný jako v předchozí metodě (rovnice 3.1). Protože je v případě kortexu femuru však přítomna jen jedna tenká vrstva (místo pěti vrstev v případě chrupavky), jsou zde jen tři jasové para-metry (pro trabekulární kost, kortikální kost a okolní tkáň) a dva poziční parapara-metry (počátek a konec kortexu podél profilu). K tomu je volně ponechán také parametr𝜎 vyjadřující předem neznámý rozměr PSF. Jasové parametry trabekulární kosti a okolní tkáně jsou při registraci nalezeny automaticky díky dostatečnému prodlou-žení jasového linie na obě strany od kortexu. Kromě hledaných pozičních parametrů tedy zbývá jasová hodnota kortikální kosti a parametr 𝜎. Závěrem autorů je, že nej-přesnějších odhadů je dosaženo v případě, že je parametr𝜎 ponechán volný a jasová hodnota kortikální kosti nastavena pevně na 1750 HU, popřípadě volena globálně pro každý femur zvlášť.
Přesnost metody byla stanovena pomocí snímání několika femurů na diagnos-tickém CT skeneru s anizotropním rozlišením 0,6 x 0,6 x 1,0 mm a 𝜇CT skenerem s izotropním rozlišením 0,08 mm. Tloušťky kortexu byly na 𝜇CT datech naměřeny původní metodou prahování na polovině jasového maxima, která je s takto vyso-kým rozlišením prokazatelně dostatečně přesná [34]. Přesnost odhadu je ve studii vizualizována mapováním na model segmentovaného femuru. Z výsledků metody na obrázku 2.10 je zřejmé, že ačkoliv metoda dosahuje dostatečné přesnosti ve většině
kosti, kde tloušťky kortexu dosahují větších rozměrů, chyba odhadu roste v přípa-dech, kdy se tloušťka blíží rozměrům okolo 0,5 mm. Pro tloušťky menší než 0,3 mm autoří přesnost neuvádí.
Obr. 2.10: Příklad mapování barevně kódovaného výsledku odhadu tloušťky kortexu a chyby tohoto odhadu na model femuru [44].
Tloušťka kortikální kosti těl obratlů se fyziologicky pohybuje mezi 180 a 600𝜇m [39]
a vyšších hodnot, které jsou spohelivě měřitelné výše uvedenými metodami, dosa-hují pouze v případech osteoblastických lézí nasedajících na kortex nebo formací osteofytů. Pro pozorování patologických změn způsobených osteoporózou je žádoucí v odhadu spolehlivě rozeznat změny okolo 50 𝜇m při hodnotách okolo spodní hra-nice fyziologického intervalu (viz kapitola 1.2). V době psaní této práce v literatuře neexistuje automatická metoda, kterou by bylo možné použít na odhad tloušťky kortexu těla obratlů.
3 ANALYTICKÉ METODY VYUŽITÉ V TÉTO PRÁCI
Tato práce si klade za cíl sestavení plně automatického algoritmu, který dokáže minerální denzitu a tloušťku kortexu těl obratlů automaticky určit. V této kapitole jsou popsány teoretické principy metod segmentace pomocí řezu grafem a regrese náhodným lesem, které jsou v navrženém algoritmu využity.
3.1 Segmentace pomocí řezu grafem
Graf𝐺lze matematicky chápat jako strukturu tvořenou množinou vrcholů𝑉 a mno-žinou hran 𝐸, kde prvky množiny 𝐸 jsou dvouprvkovými podmnožinami množiny vrcholů𝑉. Pro účely segmentace obrazu uvažujeme graf váhovaný (každá hrana má váhu z množiny reálných čísel, která vyjadřuje její hodnotu vůči ostatním prvkům) a orientovaný (rozlišujeme hranu z vrcholu𝑣𝑖 do𝑣𝑗 a hranu z vrcholu𝑣𝑗 do vrcholu𝑣𝑖).
Síť je speciálním případem grafu, jehož hrany lze považovat za systém potrubí, kde každá hrana může pojmout omezený maximální tok vyjádřený její váhou. Při zvyšování toku sítí pro její uzly platí, že přicházející tok se musí rovnat odcháze-jícímu toku s vyjímkou zdrojového a terminálního uzlu. Naplníme-li libovolnou síť maximálním možným tokem ze zdrojového do terminální uzlu, získáme množinu sa-turovaných hran, které společně zdrojový uzel od terminálního zcela oddělují. Tato množina hran se proto nazývá také řez grafem. Při vhodné reprezentaci obrazových dat pomocí sítě se řez odpovídající maximálnímu toku rovná optimální hranici mezi objektem a pozadím v obraze [7].
Konverze obrazových dat na síť probíhá následujícím způsobem. Každý voxel je reprezentován jedním uzlem, který je se svými sousedními voxely v 6-okolí propojen tzv.n-hranami. Váha hrany mezi sousedními voxelyp aq má hodnotu
𝐵𝑝,𝑞= exp(︂−(𝐼𝑝−𝐼𝑞)2 2𝜎2
)︂
· 1
𝑑𝑖𝑠𝑡(𝑝, 𝑞), (3.1) kde hodnoty Ip a Iq vyjadřují jasové hodnoty odpovídajících voxelů a 𝜎 je lo-kální směrodatná odchylka intenzit v obraze. Jedná se o tzv. boundary term, který exponenciálně roste se snižujícím se rozdílem sousedních voxelů.
Dále je každý z uzlů napojen tzv. t-hranami na zdrojový uzel představující po-myslný objekt v obraze a na terminální uzel představující pomyslné pozadí objektu.
Váhy těchto hran pro voxel p mají hodnoty
𝑅𝑝(„obj“) =−ln Pr(𝐼𝑝 | „obj“) (3.2)
𝑅𝑝(„bkg“) =−ln Pr(𝐼𝑝 | „bkg“) (3.3) a vyjadřují záporný logaritmus pravděpodobnosti příslušnosti daného voxelu do objektu nebo pozadí na základě jeho jasové hodnoty, tedy tzv. region term. Tyto pravděpodobnosti se obvykle odhadují z histogramu obrazových dat.
Libovolnou segmentaci nyní můžeme formulovat jako řez touto sítí a spočítat její energii jako
𝐸(𝐴) = 𝜆·𝑅(𝐴) +𝐵(𝐴), (3.4) kde A je množina hran protnuta tímto řezem. Region term𝑅(𝐴) zahrnuje všechny t-hrany spojující voxely segmentované jako objekt s terminálním uzlem, tedy po-myslným pozadím a všechny t-hrany spojující voxely segmentované jako pozadí se zdrojovým uzlem, tedy pomyslným objektem. Boundary term𝐵(𝐴) potom zahrnuje všechny n-hrany ležící na rozhraní objekt - pozadí.𝜆 je volitelný parametr umožňu-jící upřednostňovat jednu ze složek energie. Příklad řezu dvourozměrným obrázkem o devíti pixelech a hran, které do tohoto řezu náleží, je pro ilustraci uveden na obrázku 3.1.
Obr. 3.1: Příklad řezu dvourozměrným obrázkem reprezentovaným sítí. Řez znázor-něn fialově, hrany do něj zahrnuté oranžově. V případě zdrojového uzlu s a termi-nálního uzlu t odpovídá pravých šest pixelů objektu.
Energie řezu je tedy nejnižší v případě, že voxely zařazené do objektu mají níz-kou náležitost k pozadí, voxely zařazené do pozadí mají nízníz-kou náležitost k objektu a hranice mezi objektem a pozadím je vedena mezi voxely, jejichž jasové hodnoty se silně liší. Protože problém nalezení řezu s minimální energií je ekvivalentní k pro-blému nalezení maximálního toku v odpovídající síti, lze pro řešení použít efektivní
Ford-Fulkersonův algoritmus, který je schopný nalezení globálního minima této ener-gie [17]. Nalezení tohoto minimálního řezu vyjadřuje nalezení optimální segmentace v odpovídajících obrazových datech.