[eng]

   Po absolvování předmětu Anlýza signálů a obrazů na VUT Brno (mimochodem s výstižnou zkratkou MASO :-) jsem se nějakou dobu vrtal v problematice zpětné projekce. To je celkem jednoduchý algoritmus, kterým lze z 1D projekcí, např. rentgenových, pořízených z různých úhlů složit 2D řez objektem. Pro jeden řez to fungovalo celkem dobře, ale už tehdy mě docela nahlodávala myšlenka, jestli by šlo provádět zpětnou projekci přímo z 2D projekcí do 3D objemového modelu. Později jsem se dočetl, že se tomuto typu říká kónická zpětná projekce (cone-beam backprojection). Když jsem nedávno zahlédl na Youtube nějaká rentgenová videa rotujících objektů, tak jsem si řekl, že by se to dalo použít jako testovací data a zkusil jsem k tomu účelu napsat rekonstrukční nástroj, který je předmětem následujícího textu.

   Nejdřív jsem to zkusil v Matlabu, ale bylo to poněkud pomalé, takže jsem se rozhodl procvičit své oblíbené C/C++. Překvapivě se mi povedlo celkem rychle rozchodit obě základní metody (ray/voxel driven algorithm) a výkon byl velmi uspokojivý. Protože exportovat a načítat data do prohlížeče 3D objemových modelů 3D Slicer je poněkud neohrabané a pomalé (a navíc ten SW má řadu buggů/nepříjemných vlastností), implementoval jsem i svůj vlastní raycaster získaného modelu, což je vlastně přímá projekce, tedy přesný opak zpětné projekce a většina algoritmu tak zůstala zcela identická.

   Zpětná projekce je velmi jednoduchá metoda. V zásadě jde jen o sečtení 1D projekcí pořízených pod několika úhly do jednoho 2D obrazu. Aby byl výsledek ostrý, je před tím třeba provést filtraci (proto se metoda jmenuje filtrovaná zpětná projekce). K tomuto účelu se používá něco na způsob horní propusti. Tento filtr taky odstraní bias (střední hodnotu jasu celého řezu). Graficky znázorněný princip ukazuje obr. 1.1. Detaily v angličtině [2].

obr. 1.1 - Základní princip filtorvané zpětné projekce [2].

   3D model se touto základní technikou získává tak, že se skenovací rám nebo skenovaný objekt postupně posouvá a postupně získané řezy se pak složí do výsledného 3D objemového modelu, který lze zobrazit např. ve zmíněném 3D Sliceru (dostupný zadarmo).
   Tato metoda je ovšem poněkud pomalá a náročná na skenovací aparát. Rekonstrukci lze naštěstí provést přímo z 2D projekcí do 3D modelu, tím odpadá potřeba posuvu vzorku, stačí rotace. Této metodě se pak říká kónická zpětná projekce (cone-beam backprojection). Princip skenování je uveden na obr. 1.2.

obr. 1.2 - Proces skenování pro kónickou zpětnou projekci [4].

   Vlastní algoritmus je v tomto případě analogický: Virtuální detektor a ohnisko rentgenky se umístí v 3D prostoru přesně tak, jak byly během skenování, obraz na detektoru se "protáhne" směrem k ohnisku (podobně jako v případě 2D verze na obr. 1.1) a část vzniklého jehlanu, která protíná rekonstrukční oblast, v tomto případě krychli, se do ní jednoduše přičte. Totéž se opakuje pro všechny dostupné projekce. Nevýhoda je, že jich musí pro čistý obraz být docela dost, ale ne zase nijak extrémně (na kuře stačilo 45 projekcí na 360° rotaci).
   Praktická realizace algoritmu je možná v zásadě dvěma způsoby: paprsková a voxelová metoda (ray drivnen, voxel driven method). První metoda spočívá jednoduše ve vysílání paprsků z jednotlivých pixelů detektoru do ohniska a integraci (přičítání) intenzity aktuálního pixelu do buněk (voxelů) voxelového bufferu (3D pole), kde se vytváří výsledný objemový model. To se zopakuje pro každý pixel každě projekce. Jde tedy jen o výpočty průsečíků přímka-krychle a rasterizaci (tedy spíše "voxelizaci", když jde o 3D) paprsku do pracovního 3D pole (např. 3D Bresenhamův algoritmus). Výhodou metody je jendoduchost, také je celkem svižná a navíc je blízká fyzikální realitě skenování, takže jako vedlejší efekt proběhne váhování podle hustoty paprsků, tj. oblast blíže k ohnisku bude mít vyšší průměrnou intezitu, protože paprsky jsou blíže u sebe. Problémem metody je, že bude uspokojivě fungovat jen při velkém rozlišení projekcí (nebo malém rozlišení voxelového modelu) a velkém počtu projekcí, jinak se stane, že některé voxely budou paprsky zcela minuty nebo projde třeba jen jeden. Pak je výsledek hodně zašumněný. Řešit se to dá interpolací (upsamplingem) projekcí na výšší rozlišení, ale tím zase rostou výpočetní nároky.
   Druhá metoda jde na problém z opačného konce: začne se u voxelu a končí u pixelu projekce. Funguje to tak, že se pro každý voxel modelu vyšle virtuální paprsek z ohniska skrz tento voxel a pokud paprsek "trefí" některý pixel virtuálního detektoru, je jeho intenzita přičtena k intenzitě voxelu. Proces se opět opakuje se všemi projekcemi. Výhoda metody je, že dává čisté výsledky pro libovolné rozlišení projekcí i voxelového modelu. Nevýhoda je v tom, že je potřeba provádět váhování intenzity podle hustoty paprsků. To se celkem podepisuje na výkonu, ale lze to zanedbat, pokud bylo ohnisko umístěno dostatečně daleko od rekonstrukční oblasti, protože paprsky jsou pak téměř rovnoběžné.

   Popisovaný program je schopen načíst obrázky projekcí exportované např. ze zmíněných 3D videí, provést filtraci, zpětnou projekci a na konec zobrazit výsledek buďto po řezech nebo pomocí raycasteru (3D).

obr. 2.1 - Nástroj pro kónickou zpětnou projekci.

   Program je jako obvykle napsán v BDS 2006 Turbo C++ a byl testován ve Windoze XP/7. Je kompilován jako 32bit aplikace, takže by měl jet jak v 32bit, tak i 64bit systémech. Program používá pouze základní sadu 586, takže teoreticky jede i na starém železe, ale rozumných výsledků lze dosáhnout jen na moderních vícejádrových CPU a předevsím s rychlou pamětí! Zatímco mezi výpočetním výkonem v základních instrukcích nebude mezi nějakým 5 let starým Core 2 Duo a Intel i5 valný rozdíl, paměť udělá při té velikosti polí (stovky MB) skutečně hodně. Konkrétně jsem program ladil na Intel i5 3450 (4 jádra/4 vlákna) s 1600MHz DDR3 pamětí. Na starším Core 2 Duo s 800MHz DDR2 to v přepočtu na jedno jádro jelo asi čtvrtinovým tempem (docela překvápko). Aby program využil dnešní vícejádrové CPU, jsou pochopitelně všechny náročnější výpočty distribuované do více vláken. Poprvé jsem také zcela oddělil výpočty od VCL vlákna (GUI), takže rozhraní by mělo reagovat i během výpočtu a ten tak lze obvykle přerušit.

Download, V1.00, 14.4.2013: cone_beam_backprojection_tool_2013_04_14.zip (2,1 MB).

   Prográmek je freeware, takže může být používán i šířen jak libo, ale bez jakékoliv záruky správné funkce. Distribuce je povolena, ale jen s přiloženou složkou help, tj. tímto helpem! Autor se zříká odpovědnosti za jakoukoliv škodu způsobenou použitím tohoto programu ať už by byla jakákoliv. Použitím programu vyjadřujete souhlas s uvedenými podmínkami, v opačném případě program nepoužívejte.

   K vytvoření nového projektu je třeba někam zkopírovat INI soubor přiloženého příkladu. Lze to buď udělat ručně nebo načtením vzoru a uložením pod jiným názvem. Následně je třeba ručně vyplnit sekci [PATH] nového INI souboru (projektu). Sekce obsahuje jen dva klíče:

path - relativní nebo absolutní cesta ke složce s obrázky projekcí,
mask - wildcard řetězec pro vyfiltrování soubourů ve složce, např. "img_????.bmp".

   Obrázky jsou podporovány jen ve formátu BMP a to buďto 8-bit indexované nebo 24-bit RGB (oboje pouze ve stupních šedi). Obrázky musí být pojmenovány tak, aby po seřazení podle názvu byly v pořadí, v jakém byly skenovány a musí se jednat o projekce s pevným krokem úhlu. rozlišení všech obrázků musí být identické. Zbytek konfigurace lze provést už v programu samotném.

   První věc, co je třeba nastavit pro nový projekt, je geometri scény, tj. vzdálenost rentgenky, rozměry detektoru, úhlový krok skenování apod. Význam použitých nastavení je uveden na obr. 2.2 (červená krychle uprostřed je rekonstrukční oblast, kde vzniká objemový model).

obr. 2.2 - Geometrie scény, originální obrázek z [3].

   Ofsety du a dv jsou myšleny relativně vzhledem k rozměrům detektoru v rozsahu ±0.5. Výška detektoru se dopočítává automaticky z poměru stran obrázků projekcí. Jakékoliv změny je třeba potvrdit stiskem "Refresh" nebo klávesy Enter.

   Zde se nastavuje filtrace obrázků projekcí. První filtr v pořadí je práh a saturace v rozsahu 0 až 255, kterou lze odstranit část šumu zdrojových dat. Následuje filtrace okénkovou funkcí, kterou lze odstranit ostré přechody po obvodu projekcí. Poslední filtr je vlastní filtr zpětné projekce. Program umožňuje filtrovat vertikálně i horizontálně, ale obvykle se používá jen horizontální filtrace (u-axis). Standardní nastavení (cosine filter): f=1.0 a s=1.0. Výsledná data po všech filtracích jsou pak normalizovány na nastavenou úroveň (standardně 1024). Při nastavení "common" se všechny projekce berou jako jako jeden obrázek, na něm se provede normalizace a zase se rozseká na jednotlivé projekce. Díky tomu nedojde k vzájemné změně intenzity jednotlivých projekcí. Při nastavení "individual" se každá projekce normalizuje samostatně. To se může hodit pokud je expozice jednotlivých projekcí proměnlivá. Každou změnu je třeba potvrdit stiskem klávesy Enter (stisk Refresh neprovede přepočet filtrů!). Po přepočtu filtrů se spustí nová rekonstrukce.

   Prvním nastavením je rozlišení N v rozsahu 64 až 400. Decimace umožňuje vyřadit část projekcí pro urychlení rekonstrukce. Hodí se to při experimentálním ladění geometrie. Nastavení "Mode" umožňuje měnit metodu rekonstrukce, jak bylo popsáno v teorii. "Ray driven" metoda je obvykle rychlejší, ale hodí se maximálně pro rychlé experimentální ladění geometrie a filtrů. Kvalitu výstupu za cenu výpočetního času lze zvýšit upsamplingem. Pro finální rekonstrukci je vhodné použít metodu "Voxel driven". Ta sice zabere více času, ale dává podstatně čistší výsledky. V tomto případě lze alternativně zapnout lineární interpolaci, což se může hodit pro malé rozlišení projekcí.
   Na tomto panelu lze také nastavit počet výpočetních vláken. Standardně program nastavuje maximální počet dle CPU, ale počet lze klidně snížit. Někdy to může být i nezbytné, protože při velkých rozlišeních voxelového pole a velkém počtu vláken může být snadno překročen paměťový limit 32bit aplikací. Posledním nastavením je perioda překreslování během výpočtů.
   Rekonstruovaný objemový model může být exportován jako řezy. To lze provést tlačítkem "Save volume". Stačí nastavit cílovou složku a zadat prefix názvu, např. "obr". Program pak do složky vygeneruje obrázky jednotlivých řezů s názvy obr_000.bmp, obr_001.bmp, ... Formát je BMP, 8-bit, indexované barvy. Projekce je uložena bez post-filtrace. Soubor řezů v tomto formátu pak lze snadno načíst do 3D Sliceru.

   V tomto panelu lze nastavit orinetaci a pozici řezu v objemovém modelu. Pro potlačení šumu jsou k dispozici základní filtry (práh a gamma korekce). Klávesa Esc resetuje trackbary do defaultní polohy. Dvojklik na obrázek řezu nabídne export ve formátu BMP.

   Tlačítkem "View Projections" lze zobrazit prohlížeč obrázků projekcí. Lze zobrazit originální obrázek a jeho filtrovanou variantu.

obr. 2.3 - Prohlížeč projekcí.

   Protože používat 3D Slicer je celkem o nervy, implementoval jsem svůj vlastní raycaster. Má sice jen pár základních funkcí, ale rozhodně je jeho použití efektivnější, než po každé změně nastavení exportovat data, načítat je do 3D Sliceru a znovu ho nastavovat (ukládání projektu v 3D Sliceru v současné verzi jaksi nefunguje - vytuhne to).
   Stejně jako pro vlastní zpětnou projekci lze nastavit počet výpočetních vláken. Snažil jsem se to napsat tak, aby byl minimální overhead, takže i mezi třemi a čtyřmi vlákny při výpočetním čase 30ms je na mém PC stále znatelný rozdíl. Stisk tlačítka "Reset" obnoví nastavení prvků z INI souboru projektu. Stisk klávesy Esc na kterémkoliv trackbaru vrátí jeho defaultní hodnotu.

Fig. 2.4 - Nastavení raycasteru objemového modelu.

   Raycaster vyžaduje provést pár ručních nastavení. Ty jsou rozděleny do několika záložek a jsou uspořádány tak, jak jsou během vlastního výpočtu používány. Nejdůležitějším nastavením je práh, pomocí kterého je třeba odstranit šum okolo modelu. Dále je k dispozici nastavení "cutoff", které pro změnu umožňuje zahodit voxely přesahující nastavenou hodnotu.
   V další kartě se nastavuje zobrazovaná oblast (bounding box). Pro lepší orientaci je vhodné zobrazit si osy a "bounding box" (checkboxy vlevo dole). Zmenšení zobrazované oblasti může o něco urychlit raycasting.
   V kartě "Camera" lze nastavit úhel pohledu kamery, její polohu a měřítko. Stejné nastavení lze provádět i myšítkem v zobrazovací oblasti. Pravé myšítko slouží k posunu modelu, levé k rotaci podle základní osy Z, levé myšítko+ctrl k rotaci podle všech os, ale to je težce experimentální funkce, která zrovna moc nefunguje. Dvojklik resetuje kameru na defaultní nastavení. Pravým myšítkem lze zobrazit menu pro export (24-bit BMP), reset kamery na default nebo reset na původní nastavení dle INI souboru projektu.
   Záložka "Mode" slouží k nastavení režimu raycasteru. Základní režim jednoduše sečte všechny voxely po trase paprsku, model tak bude plně průhledný, viz obr. 2.5a. Druhý režim umožňuje podmínit začátek paprsku nějakou podmínkou, v tomto případě překročením nastaveného prahu. Po té raycasting pokračuje až do nastavené hloubky. výsledek je jen částečně průhledný model, viz obr 2.5b.

Fig. 2.5a - Plný raycasting

Fig. 2.5b - Omezený raycasting

   Poslední záložka slouží k nastavení post-filtrace. Jako obvykle je k dispozici gamma a negativní režim. Především je ale k dispozici možnost omezit úroveň saturace (bíla barva). Raycaster si standardně nastavuje úroveň bílé sám, ale lze mu vnutit vlastní hodnotu. Při zaškrtlém nastavení "Fixed" bude úroveň bíle nastavena na zadanou hodnotu i v prípadě, že skutečná úroveň bílé je nižší (to může také vést k tomu, že celý obraz bude černý, takže na toto nastavení nezapomínat).

   Zde je pár již exportovaných a nakonfigurovaných datových souborů:

Dvoupólový přepínač (naskenoval Danyk [5]), download: ct_switch.zip (3,9 MB).

Zářivkový startér (naskenoval Danyk [5]), download: ct_starter.zip (3,2 MB).

Pajšl ze zapalovače (naskenoval Danyk [5]), download: ct_lighter.zip (3,2 MB).

Pouzdro DIL8 (naskenoval Danyk [5]), download: ct_DIL8.zip (819 kB).

Budík (naskenoval Ruslan [7]), download: ct_clock.zip (10 MB).

Lebka Rysa (naskenoval Noah Spurrier [9]), download: ct_lynx_skull.zip (5,9 MB).

Kuře (naskenoval Ben Krasnow [6]), download: ct_chicken.zip (1,3 MB).

Vajíčko (naskenoval Oliver), download: ct_eier.zip (5,7 MB).

   Vzhledem k tomu, že jsem byl tentokrát líný kreslit všechno sám, použil jsem cizí obrázky z níže uvedených zdrojů. Ty také obsahují podstatně více teorie, kdyby se tím někdo chtěl zabývat:

   Zdrojová data pro ladění programu jsem získal z 3D rentgenových videí na Youtube. Našel jsem zatím tyto autory:

   Jeden příklad zdrojových videí za všechny [5]:

Last update: 25.7.2018