Registrace obrazů buněk podle intenzit v obraze

MASARYKOVA UNIVERZITA Fakulta informatiky

Registrace obrazů buněk podle intenzit v obraze diplomová práce

Ondřej Daněk 20. prosince 2006

Prohlašuji, že tato práce je mým 'původním autorským dílem, které jsem vypracoval samostatně. Všechny zdroje, prameny a literaturu, které jsem při vypracování používal nebo z nich čerpal, v práci řádně cituji s uvedením úplného odkazu na příslušný zdroj.

Rád bych poděkoval vedoucímu práce RNDr. Petru Matulovi, Ph.D. za jeho cenné rady a čas, který mi věnoval.

Shrnutí práce Práce se věnuje registraci obrazů buněk pouze na základě intenzit obsažených v ob­ razech, tedy bez využití uměle vytvořených, či jinak dodaných, registračních značek. Re­ gistrací se rozumí vyhledání geometrické transformace mezi dvěma obrazy umožňující maximální slícování objektů obsažených v obrazech. Je prezentováno několik podobnost­ ních mír použitelných pro tento způsob registrace a tyto míry jsou podrobně zkoumány a porovnávány s ohledem na jejich využitelnost při registraci dat dodaných Centrem pro analýzu biomedicínských obrazů Fakulty informatiky Masarykovy univerzity v Brně. Dále jsou testovány algoritmy pro globální optimalizaci (nalezení minima nebo maxima) těchto funkcí, přičemž je kladen důraz na přesnost a rychlost zpracování. V závěrečné části práce jsou shrnuty dosažené výsledky a doporučen nejvhodnější postup registrace. Součástí práce je také vytvoření registrační knihovny implementující všechny výše uvedené algoritmy a její podrobná dokumentace.

Klíčová slova registrace obrazů, podobnostní funkce, vzájemná informace, Fourierova transformace, optimalizace

Keywords voxel-based image registration, similarity function, mutual information, Fourier transformation, optimization

Obsah 1

Úvod 1.1 Původ zkoumaných obrazů 1.2 Vznik odchylek v obrazech 1.2.1 Barevné vady mikroskopu 1.2.2 Opakované snímání (hybridizace) 1.2.3 Působení chemikálií 1.3 Cíle práce

2

Teoretická část 2.1 Registrace obrazů 2.1.1 Formální zápis 2.1.2 Metody registrace 2.2 Podobnostní míry 2.2.1 Kritérium stochastické změny znaménka (SSC) 2.2.2 Kritérium součtu absolutních hodnot rozdílů (SAVD) 2.2.3 Kritérium součtu čtverců rozdílů (SSD) 2.2.4 Korelační koeficient (CC) 2.2.5 Normalizovaný korelační koeficient (NCC) 2.2.6 Vzájemná informace (MI) 2.2.7 Zobecněná vzájemná informace (GMI) 2.2.8 Woodsův algoritmus 2.3 Prostor transformací 2.4 Stupeň interpolace 2.5 Integrální transformace 2.6 Optimalizační metody 2.6.1 Úplné vyhledávání 2.6.2 Rekurzivní vyhledávání 2.6.3 Genetický algoritmus (GA) 2.6.4 Simulované žíhání 2.6.5 Adaptivní simulované žíhání 2.6.6 Pyramidový algoritmus

3 Výběr podobnostní míry 3.1 Registrace 2D obrazů 3.1.1 Registrace uměle transformovaných dat

i

1 1 2 2 3 3 3 4 4 5 5 6 6 7 7 7 7 8 8 9 9 10 10 10 11 11 12 13 13 13 14 14 15

11

3.2 3.3

3.1.2 Registrace reálných dat 3.1.3 Vliv degradace 3.1.4 Citlivost na parametry Registrace 3D obrazů Doporučené podobnostní míry

19 31 32 35 37

4 Výběr optimalizační metody 4.1 Registrace posunů 4.1.1 Úplné prohledávání 4.1.2 Rekurzivní prohledávání 4.1.3 Genetický algoritmus 4.1.4 Simulované žíhání 4.1.5 Adaptivní simulované žíhání 4.1.6 Integrální transformace 4.2 Registrace změny měřítka 4.3 Registrace otočení 4.4 Registrace kombinací transformací 4.5 Další možnosti optimalizace

38 38 38 39 39 40 40 40 41 41 41 42

5

43 43 44

Závěr 5.1 Výsledky 5.2 Další rozšíření

A Obsah CD

47

B Registrační knihovna B.l Podobnostní funkce B.2 Optimalizační metody B.3 Logovací třída B.4 Export grafů B.5 Hlášení chyb B.6 Statické funkce

48 48 50 52 53 53 53

C Podklady pro subjektivní hodnocení

54

D Ukázky registrace dalších dat D.l První pár D.2 Druhý pár D.3 Třetí pár

55 55 56 57

Seznam obrázků 1.1 1.2

Vznik autofocus projekce Čočka s chromatickou vadou - zdroj Wikipedia

2.1 2.2

Příklad satelitních snímků k registraci Postup rekurzivního vyhledávání

4 12

3.1

Obraz jádra buňky a tentýž obraz posunutý o 15 pixelu vpravo a 7 pixelu dolů Průběh extrémů funkcí při registraci obrazů z obr. 3.1 [část 1] Průběh extrémů funkcí při registraci obrazů z obr. 3.1 [část 2] Tatáž buňka po opakovaném snímání Histogram obrazu A z obr. 3.4 Histogram obrazu B z obr. 3.4 Průběh extrémů funkcí při registraci obrazů z obr. 3.4 [část 1] Průběh extrémů funkcí při registraci obrazů z obr. 3.4 [část 2] Průběh extrémů funkcí SSC a SAVD pro obrazy 3.4 bez vyrovnání histogramů Obrazy buněk s objekty na okraji Histogram obrazu B z obr. 3.10 Průběh extrémů funkcí při registraci obrazů z obr. 3.10 [část 1] Průběh extrémů funkcí při registraci obrazů z obr. 3.10 [část 2] Výsledek registrace obrazů z obr. 3.10 pomocí MI Průběh extrému MI při registraci změny měřítka obrazů z obr. 3.10 Výsledek registrace změny měřítka u obrazů z obr. 3.10 pomocí MI Registrace dat obsahujících nelineární transformaci intenzit Registrace 3D obrazu pomocí F-NCC

15 16 17 20 20 20 21 22 22 25 25 26 27 29 30 30 32 36

3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10 3.11 3.12 3.13 3.14 3.15 3.16 3.17 3.18

C.l Srovnání registrace pomocí NCC, MI a PIU u obrazů z obr. 3.4 C.2 Srovnání registrace pomocí NCC a MI u obrazů z obr. 3.10

ni

2 3

54 54

Kapitola 1

Úvod Tato práce vznikla ve spolupráci s Centrem pro analýzu biomedicínských obrazů Fakulty informatiky Masarykovy univerzity v Brně 1 a zabývá se registrací obrazů buněk pořízených v tamních laboratořích. Registrace obrazů spadá do oblasti zpracování dat a představuje nalezení geometrické transformace umožňující maximální možné slícování objektů obsaže­ ných ve dvou registrovaných obrazech. Z několika možných metod registrace se práce zabývá registrací pouze na základě in­ tenzit obsažených v obrazech, tedy bez využití dalších dodatečných informací (více v kapi­ tole 2). V prvním kroku jsou představeny známé, dříve prezentované, podobnostní míry2 umožňující právě tento způsob registrace a ty jsou srovnávány s důrazem na jejich použi­ telnost pro registraci dodaných dat. Důležitá je především přesnost výsledku a také doba zpracování. Dalším krokem nezbytným k nalezení optimální transformace je vyhledání globálního minima (maxima) podobnostní funkce. Prostor transformací je často velmi rozsáhlý a jeho úplné prohledávání není se současným technickým vybavením možné. V práci je proto uvedeno několik optimalizačních algoritmů pro redukci velikosti prostoru transformací a nalezení globálního optima v přiměřeném čase. Opět s přihlédnutím k přesnosti výsledku.

1.1

Původ zkoumaných obrazů

Zkoumané obrazy obsahují snímky živočišných buněk. Buňka je základním stavebním ka­ menem živých organismů, její hlavní funkcí je dělení se a tvorba bílkovin (proteinů). V já­ dře buňky je obsažen soubor chromozómů (označovaný jako karyotyp). Jejich nejdůležitější částí je DNA (deoxyribonukleová kyselina). DNA je tvořena řetězem nukleotidů a jejich sekvence (geny) hrají důležitou roli při produkci proteinů. Pro jejich produkci je dále podstatná látka RNA (ribonukleová kyselina), která přenáší informaci zapsanou v DNA řetězci do místa v buňce (ribozóm), kde proběhne vlastní syntéza bílkoviny. Chceme-li tedy porozumět chování živých organismů, je třeba zkoumat prostorovou strukturu a uspořádání DNA a její vztah k vlastní funkci buňky (tj. kdy se bude produ­ kovat jaká bílkovina). K tomu se v praxi používají 2 techniky označované jako FISH a immuno-FISH. Tyto techniky umožňují fixování buňky ve zkoumaném stavu chemickým 1 2

http://www.fl.muni.cz/lom/ Lze také použít termín podobnostní funkce. V následujícím textu budou používány oba dva výrazy.

1

KAPITOLA

1. ÚVOD

m

MAX

Obrázek 1.1: Vznik autofocus projekce

činidlem a vizualizaci (obarvení) zajímavých míst v řetězci DNA (technika FISH), respek­ tive proteinů účastnících se přepisu informace z RNA a vzniku nového proteinu (technika immuno-FISH). Obarvená místa v preparovaných buňkách je potom možné snímat mik­ roskopem. Pokud je třeba pozorovat více objektů, provádí se barvení po etapách. V každé etapě jsou obarveny pouze určité objekty, po snímání je provedeno odstranění barvy a obarvení dalších objektů. Výsledná data testovaná v této práci jsou pořízena fluorescenčním konfokálním mik­ roskopem3, a jsou ve formě třírozměrných šedotonních obrazů. Z důvodu možnosti lepší vizualizace výsledků (ve formě 3D grafů) jsou pro některé výpočty a srovnání použity pouze dvourozměrné obrazy vzniklé maximalizací intenzity obrazu ve směru osy z. Tyto obrazy jsou v textu označovány jako autofocus projekce a jejich vznik ilustruje obrázek 1.1.

1.2

Vznik odchylek v obrazech

Proč je vlastně třeba získaná data registrovat? Důvodů je hned několik. Z nich nejzásad­ nější jsou tyto: • Barevné vady mikroskopu • Opakované snímání • Působení chemikálií 1.2.1

B a r e v n é v a d y mikroskopu

Barevné vady jsou někdy též nazývány jako chromatické aberace 4 . Jedná se o vadu způso­ benou závislostí ohniskové vzdálenosti čočky na vlnové délce světla. Čočky s touto vadou pak lámou světlo každé vlnové délky jinak (viz obrázek 1.2). Bohužel ani u velmi kvalit­ ních čoček nelze tento typ vady vyloučit. Snímáme-li potom dvakrát ten stejný objekt, přičemž v každém případě ho máme obarvený barvivem vyzařujícím při nasvícení světlo různé vlnové délky, dostáváme ve výsledku dva obrázky téhož objektu, které ale mohou být vůči sobě mírně posunuty, zmenšeny, případně zvětšeny. 3 4

h t t p : / / www.fi.muni.cz/lom/pages/laboratory/hardware. shtml littp://cs. wikipedia.org/wiki/Cliromaticka_aberace

1.3. CÍLE PRACE

3

Obrázek 1.2: Čočka s chromatickou vadou - zdroj Wikipedia.

1.2.2

Opakované snímání (hybridizace)

Zřejmě největší zdroj odchylek ve snímaných datech, která pro tuto práci byla k dispozici, byl způsoben opakovaným snímáním téhož (resp. v laboratoři upraveného) preparátu. Jakmile je nutné vyjmout sklíčko s preparátem z mikroskopu, nějakým způsobem s ním manipulovat a poté jej opět do mikroskopu vrátit, není již možné dosáhnout přesně stejné pozice snímaných objektů pod mikroskopem. Dochází tak k různým posunům, natočením a dalším odchylkám. 1.2.3

P ů s o b e n í chemikálií

Jak již bylo výše uvedeno, zkoumané objekty se musí nejdříve zafixovat chemickým či­ nidlem a poté dochází k jejich postupnému barvení a odbarvování. Všechny chemikálie aplikované na buňky způsobují větší či menší změny ve vnitřní struktuře, které se opět mohou ve výsledku projevit různými posuny a jinými změnami. Navíc pro přebarvení je nutné preparát z mikroskopu vyjmout a dochází tak k problémům uvedeným v předchozí sekci.

1.3

Cíle práce

Pro účely této práce byla dodána data (vždy dvojice obrazů téhož objektu - buňky) vyka­ zující menší či větší geometrické odchylky (posuny, natočení, změny měřítka). Na těchto datech byly postupně testovány podobnostní míry a optimalizační algoritmy uvedené v této práci. Výsledky byly porovnány s ohledem na přesnost, rychlost, citlivost na parametry a další faktory a byl doporučen optimální postup pro registraci těchto dat. Vzhledem k neexistenci ucelené programové knihovny využitelné pro registraci obrazů je součástí práce i implementace všech, v práci uvedených, podobnostních mír a optimali­ začních algoritmů, s možností generování grafických průběhů a zaznamenávání informací o době registrace, nastavení parametrů, zkoumaném transformačním prostoru a dalších.

Kapitola 2

Teoretická část 2.1

Registrace obrazů

Registrace (lícování) obrazů, je definována jako hledání geometrické transformace, která slícuje body a objekty jednoho obrazu s odpovídajícími body a objekty druhého obrazu. Oba obrazy přitom typicky představují pohled na stejný objekt, který ale mohl být pořízen různými pořizovacími metodami nebo v různých časech. Registrace obrazů je v současnosti nejvíce využívána pro následující účely: • Lékařství — Registrace dat pořízených kombinací několika pořizovacích senzorů. Tedy například morfologických dat pořízených počítačovou tomografií (CT), jader­ nou magnetickou rezonancí (MR) a funkčních dat pořízených jednofotonovou emisní tomografií (SPÉCT) či pozitronovou emisní tomografií (PET). Registrace dat z ně­ kolika modalit je podmínkou jejich důkladnějšího vizuálního zkoumání. • Hledání změn — Porovnávání snímků pořízených v různých časech nebo za jiných podmínek. Například registrace satelitních snímků jak ukazuje obrázek 2.1. • Porovnávání a rozpoznání objektů — Například ověření podpisu, či segmentace (čle­ nění) objektů.

Obrázek 2.1: Příklad satelitních snímků k registraci

4

2.1. REGISTRACE 2.1.1

OBRAZŮ

5

Formální zápis

Mějme dány dva obrazové soubory / a J a funkci S oceňující registraci obrazů, které jsou jí předány jako parametry. Předpokládejme dále, že optimální registrace je dosažena v minimu této funkce. Proces registrace poté můžeme definovat vztahem: ropt = a r g m m 5 ( / , r ( J ) ) ,

(2.1)

kde T vyjadřuje geometrickou transformaci, T je množina všech uvažovaných transformací (která je před výpočtem také určitým způsobem specifikována) a Tcypt je transformací v níž je dosaženo optimální registrace. Analogicky lze registraci definovat i pro funkce dosahující optimální registrace v maximu. 2.1.2

M e t o d y registrace

Podle způsobu jakým postupujeme můžeme vyčlenit 5 nejvýznamnějších metod regis­ trace [1], [2]: • Interaktivní metody • Registrace pomocí korespondencí značek • Korespondence hranic oblastí • Maximalizace globální podobnosti • Registrace využívající integrální transformace V následujícím textu budou jednotlivé metody stručně představeny. Tato práce se pak detailně zabývá posledními dvěmi jmenovanými. Interaktivní metody Interaktivní metody představují lícování obrazů pod vizuální kontrolou operátora. Uplat­ ňují se především v případech kdy by bylo obtížné použití ostatních metod. Tedy přede­ vším při silné deformaci zobrazených objektů, malé velikosti překrývající se oblastí, nízkém kontrastu apod. Registrace pomocí korespondencí značek Tento typ registrací se v anglické literatuře označuje jako point-based registrace. V duchu klasického členění registrací (např. devítidimenzionální schéma pro klasifikaci registrací v medicínském prostředí dle Maintze a Viergevera [2]) se typicky jedná o vnější registrace, kdy jsou do registrovaných dat zavedeny dodatečné informace, například ve formě značek (navrtáním dat laserem), které nám usnadní nalezení optimální transformace. Případně je rovněž možné využít jako registrační značky významná místa v obraze (rohy, křižo­ vatky, apod.). Tyto metody lze dále rozdělit na metody interaktivní (značky jsou v obraze označeny operátorem) a neinteraktivní (vyhledání značek provádí sám stroj segmentací).

KAPITOLA

6

2. TEORETICKÁ

CÁST

Korespondence hranie oblastí Při těchto metodách je základním předpokladem segmentace objektů v obrazech a sta­ novení jejich hranic. Následně dochází k aplikaci algoritmů umožňujících, na základě na­ lezených hranic, určení hledané registrační transformace. Typickými představiteli těchto algoritmů je například Houghova transformace či metoda „klobouk a hlava". Maximalizace globální podobnosti Tyto metody se v anglické literatuře označují jako voxel-based registrace. Registrace dat v tomto případě závisí na zvolení vhodné podobnostní míry a její globální optimalizaci (hledání globálního maxima nebo minima) na prostoru uvažovaných transformací. Vy­ hodnocování algoritmů tak probíhá pouze na základě intenzit voxelů či pixelů obsažených v obrazech, přičemž se využívá informace z celého datového souboru. Výhodou těchto metod je, že není vyžadována segmentace a většinou ani výrazné předzpracování obrazů. Registrace využívající integrální transformace I tyto metody lze zahrnout pod pojem voxel-based registrace. Využívají skutečnosti, že integrální transformace (například Fourierova) mají v případě translace, rotace a změny měřítka své transformanty ve frekvenční oblasti. Tedy na základě výpočtu transformace obrazů je možné efektivně určit optimální registraci. Jakým postupem toho dosáhnout je uvedeno dále v této kapitole.

2.2

Podobnostní míry

Podobnostní míry (nebo též funkce) tvoří základ voxel-based registrací. Jak vyjadřuje rov­ nice 2.1 jedná se o binární funkce, které pro dané dva obrazy vrací míru jejich podobnosti. Optimální registrace je vždy dosaženo buď v globálním maximu funkce nebo v jejím mi­ nimu. Pro účely této práce se uvažuje pouze druhý případ, tedy funkce dosahující optimální registrace v maximu jsou násobeny hodnotou —1 a tím je jejich optimalizace převedena na minimalizační úlohu. Ze známých podobnostních funkcí budeme podrobně zkoumat těchto osm: kritérium stochastické změny znaménka, kritérium součtu absolutních hodnot rozdílů, kritérium součtu čtverců rozdílů, korelační koeficient, normalizovaný korelační koeficient, vzájemnou informaci, zobecněnou vzájemnou informaci a dvě míry jejichž autorem je R. P. Woods (jedná se o Ratio-Image uniformity a Partitioned Intensity Uniformity). 2.2.1

K r i t é r i u m stochastické z m ě n y z n a m é n k a ( S S C )

Podstatou tohoto kritéria [4] je sčítání počtu průchodů nulou při výpočtu rozdílů intenzit (hodnot jasu) korespondujících obrazových bodů. Kritérium je odvozeno z neparametrické statistické teorie a obrazy jsou registrované v případě maximálního počtu průchodů nulou. Metoda předpokládá, že slícované obrazy jsou identické až na nekorelovaný šum s nulovou střední hodnotou, který má symetrické rozložení hustoty pravděpodobnosti. Tedy tato metoda uvažuje přítomnost šumu v obrazech.

2.2. PODOBNOSTNÍ

2.2.2

MIRY

7

Kritérium součtu absolutních hodnot rozdílů (SAVD)

Toto kritérium [5] je definováno vztahem: SAVD(I, J,T) = ±J2

\W - T(Jm,

(2-2)

i

kde I a J jsou registrované obrazy, T je zkoumaná transformace, I (i) představuje prvek na pozici i v obrazu I (i může být i vektor několika souřadnic pro vícedimenzionální obrazy) a T(J)(i) představuje prvek na pozici i v obrazu který vznikne aplikací transformace T na obraz J. N je celkový počet obrazových elementů v průniku obrazů. Funkce SAVD je relativně výpočetně nenáročná. Uvažujeme-li i normalizační faktor 1/N je obor funkčních hodnot nezávislý na velikosti registrovaných dat. Pro šedotónová (grey-scale) data nabývá hodnot v intervalu (0,255). Obrázky jsou registrovány pro mini­ mum funkce. 2.2.3

K r i t é r i u m s o u č t u čtverců rozdílů ( S S D )

Kritérium vychází z kritéria SAVD a je definováno vztahem: SSD(I, J,T) = J2 (I(i) - T{J){i)f.

(2.3)

i

Předpokladem je, že k optimální registraci dochází jestliže obrazy / a T(J) jsou identické. Pro tyto obrazy dostáváme hodnotu SSD rovnou nule. Pokud obrazy nejsou slícovány vychází hodnota SSD > 0, tedy optima je dosaženo pro minimum funkce. 2.2.4

Korelační koeficient ( C C )

Vzájemná korelační funkce je definována vztahem:

cc(i,j,T) = j2mnj)(i).

(2.4)

i

Výhodou korelačního koeficientu je opět jeho nízká výpočetní náročnost, naopak jeho nevýhodou je menši robustnost v porovnání s některými ostatními funkcemi (viz kapitola 3), protože funkce není normalizovaná. Optimální registrace je dosaženo v maximu funkce. 2.2.5

N o r m a l i z o v a n ý korelační koeficient ( N C C )

Normalizovaný korelační koeficient [6] eliminuje většinu nedostatků předchozí podobnostní míry, při zachování obdobné výpočetní složitosti. Je definován vztahem:

VEÍ

(/(«) - /*)

EÍ

( J M - J*)

kde /* a J* jsou průměrné hodnoty v obrazech I a J. Funkce popisuje lineární závislost dvou souborů dat a nabývá hodnot v rozmezí (—1,1). Čím více se hodnota blíží jedné tím lepší je registrace. Funkce vykazuje velmi dobré vlastnosti a proto je jednou z nejpoužíva­ nějších pro měření podobnosti obrazů.

KAPITOLA

8

2.2.6

2. TEORETICKÁ

CÁST

V z á j e m n á informace (MI)

Vzájemná informace (Mutual information) [7] je netradičním způsobem měření podobnosti dat. Odhaduje obecnou (tj. i nelineární) závislost dvou souborů dat. Registrace pomocí vzájemné informace je tedy robustnější a je vhodná i pro registraci dat získaných z různých zobrazovacích modalit. Algoritmy pro její výpočet jsou velmi obecné a jsou použitelné pro široké spektrum dat. Výpočet vychází ze statistické teorie — 2 náhodné veličiny x a y jsou statisticky nezávislé právě tehdy, když jejich simultánní hustota pravděpodobnosti je součinem jejich marginálních hustot. Tedy platí: p(x,y)

=Px(x)py(y)

a maximálně statisticky závislé jestliže platí: p(x,y) =px(x)

=Py(y).

Vzájemná informace náhodných veličin I a T(J) měří stupeň jejich závislosti KullbackLeiblerovou vzdáleností mezi sdruženou hustotou pravděpodobnosti pu(I(i),J(i)) a sou­ činem marginálních hustot pi(I(i))pj(J(i)). Vypočítá se vztahem (jednotkou jsou bity):

MW,J,T> = ^ W W W ^ ^ f f j ; » , .

(2-6)

Optimální registrace je dosaženo pro maximum funkce. Jiný způsob vypočtu vzájemné informace využívá Shannonových entropií: MI(I, J, T) = H(I) + H(J(T))

- H(I, T(J)),

(2.7)

kde H(I) a H(T(J)) jsou Shannonovy entropie veličin I a T(J) a H(I,T(J)) jejich vzájemná entropie: H{I) = - $ > / ( / ( * ) ) log 2 pj(/(i)), (2-8) i

H{T{J)) = - £ P j ( T ( J ) ( i ) ) l o g 2 í > j ( T ( J ) ( i ) ) ,

(2.9)

i

H(I,T(J))

= -J2'PiÁKý),T(J)(ý))^g2Pij(m,T(J)(i)).

(2.10)

i

2.2.7

Z o b e c n ě n á v z á j e m n á informace ( G M I )

Definice zobecněné vzájemné informace [8] je v mnohém podobná definici MI, pouze Shan­ nonovy entropie ve vztahu 2.7 jsou nahrazeny zobecněnými Rényiho entropiemi druhého řádu: GMI(I,J,T) =H{-2\l) + H{-2\j{T))-H{-2\l,T{J)), (2.11) kde

tf(2)(/) = -log 2 ][>?(/(*)),

(2.12)

i

HW(T(J))

=-log^^CZVX*)),

(2.13)

2.3. PROSTOR

TRANSFORMACI H^\I,T{J))

9 = -log2J2pÍAm,T(J)(i)).

(2.14)

i

Také GMI dosahuje optimální registrace pro maximum funkce a stejně jako MI vyka­ zuje výborné vlastnosti i při registraci dat vykazujících nelineární závislosti. 2.2.8

W o o d s ů v algoritmus

Následující dva algoritmy byly navrženy Woodsem [9], [10] především pro registraci PETPET a MR-PET objemových dat. Ukazuje se ale, že jsou využitelné i při řešení obecnějších problémů. Podobně jako vzájemná informace dokáží eliminovat nelineární transformace intenzit. Názvy algoritmů jsou uváděny v anglickém jazyce. Ratio-Image Uniformity (RIU) Algoritmus Ratio-Image Uniformity je používán především pro registraci obrazů poříze­ ných stejnými modalitami. Jeho jádrem je následující postup výpočtu:

R(l) =

WŘj

R* = ^J2R(^> *« = Jj^rrEW)-**) 2 , RIU = g .

(2-15) (2.16)

Registrace dat je dosaženo pro minimum funkce. Partitioned Intensity Uniformity Algoritmus Partitioned Intensity Uniformity navíc umožňuje registraci dat pořízených růz­ nými modalitami, tedy například MR-PET registrace. Je založený na myšlence, že jsou-li obrazy registrované, pak pro danou hodnotu voxelu (objemového elementu) souboru / je rozsah intenzit odpovídajících voxelu v souboru J minimální. Výpočet podobnosti je následující:

p{a) =

Eb bh{a, b)

ňH^ř

,

a{a) =

1

J2bh(a,b)(b-p*(a))2

mi—ňH^J)—'

PIU = YJW{a)^bhi^A)l

(2 17)

-

(2.18)

a

kde h(a, b) je vzájemný histogram obrazů I a J. Registrace dat je dosaženo pro minimum funkce.

2.3

Prostor transformací

Velikost prostoru prohledávaných transformací je klíčovým prvkem při registraci obrazů. Pokud je prostor příliš velký nemusí se podařit najít globální minimum podobnostní funkce nebo může jeho nalezení zabrat neúměrně dlouhou dobu. Převážná část této práce se bude věnovat registraci posunů (translací), tedy pro 2D obrazy dostáváme dvourozměrný prostor

KAPITOLA

10

2. TEORETICKÁ

CÁST

transformací a pro 3D obrazy třírozměrný. V závěrečných částech bude zmíněn i postup pro registraci anisotropní 1 změny měřítka (scale) a rotací.

2.4

Stupeň interpolace

Dalším faktorem ovlivňujícím výsledek je stupeň interpolace používaný při transformaci obrazů. Interpolaci stupně 0 se také označuje termínem interpolace typu nejbližší soused. Dalším stupněm je interpolace lineární, kdy se pro určení hodnoty voxelu používá inter­ polace mezi dvěma sousedními voxely. Následují interpolace kvadratická (též bilineární), kubická, atd. Práce se zabývá nejen vlivem stupně interpolace na transformaci obrazů, ale také na důsledky použití vyšších stupňů interpolace při převzorkování dat.

2.5

Integrální transformace

Klasický výpočet korelační funkce a normalizované korelační funkce se provádí v prosto­ rové oblasti. Je však možné počítat je také s použitím korelačního teorému Fourierovy transformace ve frekvenční oblasti: CC(I, J, T) = F-\F(I)F*(T(J))),

(2.19)

kde F (I) je Fourierova transformace obrazu I, F*(J) je komplexně sdružená Fourierova transformace obrazu J a F~l je inverzní Fourierova transformace. Obdobně NCCd J T) = l ' ' j y/F~^F(I

F-\F(I-P)*F*(T(J)-J*)) - I*) * F(K))F-i(F(T(J) - J*) * F{K))'

l

'

kde /* a J* jsou střední hodnoty obrazů / a J a K je obraz sestávajících ze samých intenzit s hodnotou 1. Z teorie Fourierových transformací vyplývá ještě jeden důležitý důsledek: posun v ob­ razové doméně se ve frekvenční oblasti projeví posunem fáze, přičemž amplitudy zůstanou zachovány. Můžeme tedy definovat funkci fázová korelace (PCC):

PCCd JT)- F-iflVhí^rWL)

(2 21)

FCC(l,J,l)-r (|F(/)MF(T(J))|)I 2 - 21 ) Výsledkem je delta funkce se středem v bodě určujícím posun registrovaných obrazů, tedy v místě optimální transformace. Velkou výhodou je, že při výpočtu výše uvedených funkcí integrální transformací automaticky získáme hodnoty podobnostní funkce pro celý obor transformací a není již třeba aplikovat žádné další optimalizační algoritmy. Nevýhodou je potom, že bez dalších úprav nedokážeme integrální transformací registrovat jiné transfor­ mace než posuny (translace).

2.6

Optimalizační metody

Optimalizace lze definovat [11] například jako obor zabývající se určením nejlepšího řešení jistého matematicky definovaného problému. Pro naše účely budeme optimalizací rozumět 1

Změna ve směru jednotlivých os není stejná

2.6. OPTIMALIZAČNÍ

METODY

11

nalezení globálního optima (minima či maxima) podobnostní funkce pro dané dva regis­ trované obrazy na oboru uvažovaných transformací. Je zřejmé, že pro rozsáhlý prostor transformací se jedná o velmi netriviální úlohu. Existuje celá řada optimalizačních metod, z nichž se v této práci zaměříme na úplné vyhledávání, rekurzivní vyhledávání, genetické algoritmy, simulované žíhání, adaptivní simulované žíhání a pyramidový algoritmus. Z dalších významných metod jmenujme především metody největšího spádu založené na výpočtu gradientu, kterými se ale tato práce nezabývá. Je to jednak z důvodu nemož­ nosti analyticky vyjádřit optimalizované funkce a také z důvodu možného výskytu mnoha lokálních minim v těchto funkcích, znesnadňujících použití gradientních metod. Nicméně jejich použití bude diskutováno v závěru práce. 2.6.1

Ú p l n é vyhledávání

Nejjednodušším možným postupem je úplné prohledávání prostoru transformací. Výhodou tohoto postupu je to, že nám zaručuje nalezení globálního maxima podobnostní funkce. Tuto výhodu ovšem převažuje několik zásadních nedostatků: • Výpočetní náročnost — jak už bylo několikrát zmíněno, prostor transformací může mít, a také většinou má, mnoho dimenzí (registrace v této práci pracují s prostorem transformací až do dimenze 9). Už malý rozsah parametrů v jednotlivých dimenzích nám způsobí, že nebudeme schopni použít tento postup v rozumném čase. Jako pří­ klad může uvažovat 3D obraz o rozměrech 200 x 200 x 40 voxelů a požadavek registro­ vat posuny v intervalu (—50; 50) voxelů ve směru všech os. Jednoduchým výpočtem dostáváme, že podobnostní funkci bude třeba vyhodnotit 1000 000 x pokaždé pro obraz obsahující 1600 000 voxelů. Už takovéto, na první pohled jednoduché, zadání povede na současném běžně používaném hardwaru k výpočtům trvajícím několik desítek minut, spíše však hodin, což je neakceptovatelné. • Spojité parametry — z podstaty definice úplného vyhledávání je problémem re­ gistrace spojitých parametrů. Typickým příkladem je změna měřítka, protože ta představuje nekonečně mnoho testovatelných transformací. Pro takovéto parametry je třeba zavést předem definovanou velikost kroku, po které bude algoritmus po­ stupovat. Tento „trik" nám pomůže i v případě, kdy potřebujeme touto metodou registrovat například posuny se sub-voxelovou přesností. 2.6.2

Rekurzivní vyhledávání

Rekurzivní (někdy též stupňové) vyhledávání eliminuje obě nevýhody úplného vyhledá­ vání. Daní, kterou musíme zaplatit, je ztráta přesnosti výpočtu a možnost uváznutí vyhle­ dávání v bodě lokálního minima podobnostní funkce. Postup algoritmu názorně ilustruje obrázek 2.2 V každém kroku je rozsah jednoho parametru transformace rozdělen na části a funkce se vyhodnocuje jen v předem stanoveném počtu bodů (v prvním kroku to jsou body ozna­ čené 1). Je nalezen bod s nejnižší hodnotou podobnostní funkce a v jeho okolí se rekurzivně pokračuje ve vyhledávání (kroky 2 a 3). Algoritmus zastaví buď po daném počtu kroků nebo jakmile vzdálenosti testovaných bodů klesnou pod stanovenou hranici. Je zřejmé, že tímto postupem dosáhneme razantní redukce výpočetní náročnosti a není také třeba

KAPITOLA

12

'•

t

T

^r a_

r-

2. TEORETICKÁ

CÁST

T 1 --^..

'3Z_a 3--apň-a.4-fjqUB

^Xfzfí

3 33

T

T

L

"1

1

+"

Obrázek 2.2: Postup rekurzivního vyhledávání

řešit problém spojitých parametrů. Při správném nastavení může být rekurzivní vyhledá­ vání velice efektivním postupem. Záleží ovšem zejména na tvaru podobnostní funkce, nad kterou algoritmus běží. 2.6.3

G e n e t i c k ý algoritmus ( G A )

Genetické algoritmy [3] jsou inspirovány evoluční biologií. Evoluce je ve své podstatě me­ todou, která hledá optimální řešení nad množinou velkého množství možností. Genetické algoritmy se snaží simulovat evoluční procesy ve snaze dosáhnout optima řešeného pro­ blému. Vlastní algoritmus má dvě fáze. V první fázi je náhodně vygenerován předem stanovený počet jedinců (prvků prostoru možných řešení), označovaný jako populace nebo též první generace. V dalších fázích se poté neustále opakuje proces tvorby další generace z generace předchozí, přičemž je snahou nalézt lepší řešení než v generaci předchozí. Celý proces končí po dosažení předem daného počtu generací nebo nalezení dostatečně dobrého řešení. Tvorba další generace sestává z několika opakujících se kroků. Nejprve jsou jedinci seřazeni podle hodnoty jejich ohodnocovací funkce (v terminologii GA označována jako fitness funkce, v našem případě se jedná o míru podobnosti). Následně jsou vybráni 2 jedinci (čím větší fitness, tím větší pravděpodobnost výběru), ti jsou s určitou pravděpo­ dobností (obvykle více než 50%) zkříženi a s menší pravděpodobností (obvykle do 5%) je u nich provedena mutace. Výslední jedinci přechází do nové generace. Je také možné do nové generace přímo přenést jedince s velmi dobrou hodnotou fitness a zachovat je tak pro další evoluci (tzv. elitismus). Tyto kroky jsou opakovány tak dlouho dokud není v nové

2.6. OPTIMALIZAČNÍ

METODY

13

populaci dostatečný počet jedinců. Jednotlivé prvky prostoru řešení (v našem případě prostoru transformací) jsou obvykle zakódovány do bitových řetězců, což nám usnadňuje proces křížení a mutace. Při křížení jsou bitové řetězce jedinců rozděleny v určitém místě a jejich části vyměněny. Při mutaci dochází k jednobitové změně na náhodném místě v bitovém řetězci. Testování dokazují, že genetické algoritmy jsou řešením velice robustním. Jejich velkou nevýhodou je však pomalá konvergence k vrcholu extrému funkce. Často tak nalézáme řešení, které se pouze značně blíží řešení optimálnímu. 2.6.4

Simulované žíhání

Dalším algoritmem inspirovaným přírodou a přírodními ději, konkrétně chováním termo­ dynamických systémů, kdy se změnou teploty látka mění skupenství, tj. např kov je žíhán a ochlazován je simulované žíhání [1, 19]. Je-li tavenina ochlazována pomalu, dochází k tvorbě krystalické struktury, která představuje minimální energetický stav systému. Metoda byla vyvinuta hlavně pro řešení třídy problémů, kdy je vyhledáván globální extrém v okolí mnoha extrémů lokálních. Je založena na tzv. Boltzmannově distribuci: P(E) ~ e ~ * r , kde P(E) udává pravděpodobnost, že při teplotě T přejde systém do energetického stavu E a k je Boltzmannova konstanta. Tedy i při nízké teplotě je zde malá pravděpodobnost, že systém přejde do vyššího energetického stavu. V řeči našich funkcí to lze vyjádřit tak, že algoritmus je schopný přejít z lokálního minima přes místa s vyšší energií do míst s ještě menší energetickou hladinou. To jej zásadně odlišuje například od metody rekurzivní. 2.6.5

A d a p t i v n í simulované žíhání

Metoda adaptivního simulovaného žíhání dále vylepšuje algoritmus simulovaného žíhání. Její předností je jednak teplotní schéma s rychlejším ochlazování a také přizpůsobování žíhacího předpisu zvlášť pro každý parametr transformace podle jeho citlivosti. Protože se jedná o algoritmus poměrně značně komplikovaný a jelikož jeho zkoumání není předmětem této práce, odkazuji zde pouze na internetové stránky autora (viz [12]), kde je k nalezení podrobný popis. 2.6.6

P y r a m i d o v ý algoritmus

Jak uvádí Čapek [1], velmi progresivním způsobem urychlení nalezení globálního extrému je pyramidový algoritmus. Je vybudována pyramida (série) 2D nebo 3D obrazů s obrazem v originálním rozlišení na dně pyramidy a s obrazy se snižujícím se rozlišením na vyš­ ších stupních. Nejprve se registrují obrazy na vrcholu pyramidy (s nejnižším rozlišením) a je nalezeno globální minum. V dalších patrech pyramidy se pozice globálního extrému už pouze zpřesňuje vyhledáváním v okolí extrému nalezeného v předchozím stupni. Ne­ výhodou této metody je opět riziko uváznutí v lokálním minimu, nejčastěji díky chybné registraci v některém z vyšších pater pyramidy. Tomu lze předejít vhodným nastavením algoritmu a použitím dostatečně velkého rozlišení na vrcholu pyramidy tak, aby byly za­ chovány důležité rysy obou registrovaných obrazů.

Kapitola 3

Výběr podobnostní míry V této kapitole se pokusíme srovnat jednotlivé podobnostní míry a doporučit ty nejvhod­ nější pro registraci dodaných obrazů buněk. Pokud nebude možné jednoznačné doporučení, pokusíme se alespoň formulovat podmínky, za kterých jsou funkce použitelné a dávají ale­ spoň uspokojivá řešení. K tomu, abychom to mohli uskutečnit bylo třeba provést dostatečný počet testů na velkém množství dodaných obrazů. Pro každý testovaný obraz byly zkoumány vlastnosti podobnostní míry a faktory ovlivňující její průběh a pozici globálního extrému. Největší důraz byl samozřejmě kladen na schopnost podobnostní míry dosáhnout svého globálního extrému v místě optimální registrace obrazů a také na dobu výpočtu. Další vlastností na kterou je vhodné se zaměřit, je Čapkem [1] zavedená zóna atrakce funkce1, která vyjadřuje velikost oblasti (vyjádřené v jednotkách voxelů), z jejíhož každého bodu se lze dostat k vrcholu globálního minima pomocí optimalizace založené na spádu gradientu, neboli oblasti vymezené hranicí, od které by se pomyslná kulička po ploše grafu průběhu funkce ještě skutálela k vrcholu minima. Čím větší má funkce svou ZAF, tím snadnější je lokalizace globálního extrému. Čapek ve své práci uvádí funkce pro Matlab 2 , pomocí nichž je možné odhadnout velikost ZAF. My se však spokojíme pouze s vizuálním odhadem. Pro registraci větších obrazů pokročilými optimalizačními metodami (například pyra­ midovým algoritmem), bylo nutné přistoupit k určitému předzpracování registrovaných obrazů (převzorkování na menší rozměry, změna počtu kvantizačních úrovní 3 ,...). Proto se v této části práce zaměříme i na citlivost podobnostních funkcí a jejich zkoumaných vlastností na výše uvedené změny v registrovaných datech.

3.1

Registrace 2D obrazů

Jak už bylo řečeno v úvodní kapitole, dodaná data jsou ve formě třírozměrných šedotonních obrazů. Pro prozkoumání vlastností podobnostních funkcí se však jeví jako vhodnější použít 2D obrazy, na kterých lze lépe vizualizovat průběh podobnostní funkce. Registrace byly proto primárně testovány na autofocus projekcích (viz úvodní kapitola) původních 1

V dalším textu označována jako ZAF. h t t p : / / www. mathworks.com/ 3 Počet úrovní intenzity nebo též šedých úrovní obrazu. 2

14

3.1. REGISTRACE

2D OBRAZŮ

15

dat. Pro tvorbu grafů a nalezení globálního minima v této a následující podkapitole byly použity metody úplného vyhledávání a integrálních transformací. 3.1.1

R e g i s t r a c e u m ě l e transformovaných dat

Pro počáteční analýzu podobnostních funkcí použijeme uměle transformovaná data, u kte­ rých budeme dopředu znát pozici registračního optima. To nám usnadní hodnocení funkcí. Jako registrovaný soubor byl vybrán jeden z obrazů buněčného jádra o rozměrech 149 x 149 pixelů, který byl v grafickém editoru posunut o 15 pixelů ve směru osy x a o 7 pixelů ve směru osy y, čímž vznikl druhý obraz. Navíc byl v malé míře přidán do obou obrazů Gaussův šum. Oba registrované obrazy s uměle vytvořeným posunem znázorňuje obrázek 3.1:

Obrázek 3.1: Obraz jádra buňky a tentýž obraz posunutý o 15 pixelů vpravo a 7 pixelů dolů Tyto obrazy byly postupně registrovány s použitím všech podobnostních mír uvede­ ných v teoretické části práce. Jako prostor transformací byly uvažovány posuny (translace) v intervalu (—40, 40) pixelů a to jak pro osu x, tak pro osu y. Podobnostní míra byla vyhod­ nocena pro všechny transformace s přesností na 1 pixel, to znamená, že byla vyhodnocena celkem 6400 x. Korelační koeficient a normalizovaný korelační koeficient byly, kromě standardního po­ stupu, vypočítány i pomocí integrální transformace. Jelikož integrální transformace v sobě implicitně zahrnuje předpoklad periodicity obrázku (tj. například pro obrázek o šířce 150 pixelů předpokládá na pozici 151 pixel se stejnou intenzitou jako je na pozici 1) mo­ hou se výsledky lišit od funkcí počítaných standardním postupem. Pokud chceme počítat funkce bez periodicity je třeba doplnit obrazy nulovou intenzitou až do dvojnásobku je­ jich rozměrů. Důkladnější rozbor kdy uvažovat periodicitu a kdy ne, provedeme až po prostudování grafů průběhů podobnostních funkcí. Ty zobrazují obrázky 3.2 a 3.3. Je ještě třeba poukázat na to, že všechny podobnostní míry, s výjimkou těch počítaných přes Fourierovu transformaci, byly upraveny tak, aby dosahovaly registračního optima ve svém globálním minimu. Míry počítané přes integrální transformaci dosahují registračního optima ve svém globálním maximu.

16

KAPITOLA

3. VÝBĚR PODOBNOSTNÍ

-9421 pra TX [ptal] = -14. TY [ptali =

= 4 37509 pra IX |pÍHel| = -15. TY [phel] =

Sum of squared differences; extrem = 729S20 pro TÁ [pixel] = -15, TY [pixel] = -7

524e-HII8 pro TX [pixel] = -15, TY [pixel] = -7

>Ä

f

Normalized cross correlation; extrem - -0.99414G pro TX [pixel] - -15, TY [pixel] - -7

Generalized mutual information; exlrem - -1.27354 pro TX [pixel] - -15. TY [pixel] - -7

Woods (Ratio-Image Uniformuj); extrem -

Obrázek 3.2: Průběh extrémů funkcí při registraci obrazů z obr. 3.1 [část 1]

MÍRY

3.1. REGISTRACE

2D OBRAZŮ

Fourier NCC; extrem = 0.993331 pro TX [pixel] = -15. TY [pixel] = -7

TY [pixel]

K [pixel]

17

Fourier PCC; extrem = 3.497792 pro TX [pixel] = -15, TY [pixel] = -7

TY [pixel]

TX [pixel]

Obrázek 3.3: Průběh extrémů funkcí při registraci obrazů z obr. 3.1 [část 2]

Na první pohled se ukazuje, že všechny funkce (s drobnou odchylkou u SSC) nalezly optimální registrační transformaci správně (druhý obraz je třeba posunout o 15 pixelů doleva a o 7 pixelů nahoru). Pro lepší porovnání si proto uveďme ještě dosažené výpočetní časy potřebné k registraci obrazů a nalezené pozice extrémů. Ty nám ukazuje tabulka 3.1. Výpočetní doby jsou samozřejmě velmi závislé na konkrétní implementaci, optimali­ zaci pro konkrétní typ obrazů a také na použitém hardware. Uvedné časy jsou měřeny za použití registrační knihovny, která vznikla v rámci této práce. Tato knihovna umož­ ňuje registrovat 2D a 3D obrazy s jedním barevným kanálem. Z možných transformací lze registrovat posuny, změny měřítka, otočení a jejich kombinace. Umožňuje také nastavit stupeň interpolace použitý při transformacích obrazů a převzorkování obrazů. Knihovna je snadno rozšiřitelná o další podobnostní míry a optimalizační metody a jedná se tedy o implementaci poměrně robustní. Použitým hardware byl počítač osazený procesorem Athlon 64 3000+ (běžícím na frekvenci 1,8 GHz) s 1 GB operační paměti typu RAM a operačním systémem Microsoft Windows XP. Pro snazší orientaci si rozdělme použité podobnostní míry do 5 skupin: • Skupina 1 (SSC, SAVD, SSD) - Míry založené na rozdílech intenzit. Vidíme, že SAVD a SSD se téměř neliší, doba výpočtu je u obou stejná, pouze SSD vykazuje o něco hladší průběh grafu. ZAF obou funkcí je tvořeno bezmála celou plochou grafu, tedy výsledek značně uspokojivý. SSC se nepatrně odlišuje. Na první pohled překvapí o sekundu delší doba výpočtu, která je však dána nutností testovat průchody nulou,

18

KAPITOLA

3. VÝBĚR PODOBNOSTNÍ

MÍRY

Tabulka 3.1: Doba výpočtu a pozice extrému podobnostních funkcí při registraci dat z obr. 3.1 Podobnostní míra Doba registrace Pozice extrému

ssc

8s

(-14, -7)

SAVD

7s

(-15, -7)

SSD

7s

(-15, -7)

CC

7s

(-15, -7)

NCC

9s

(-15, -7)

MI

16 s

(-15, -7)

GMI

14 s

(-15, -7)

RIU

8s

(-15, -7)

PIU

14 s

(-15, -7)

F-CC (s periodicitou)

0s

(-15, -7)

F-NCC (s periodicitou)

0s

(-15, -7)

F-PCC

0s

(-15, -7)

což znesnadňuje moderním procesorům jejich výpočetní optimalizace. Nicméně rozdíl v době výpočtu v řádu sekund rozhodně není klíčovým. Problémem u kritéria SSC je, že jako jediné nedokázalo přesně určit optimální registraci (liší se v ose x o jeden pixel). Z grafu je patrné, že SSC vykazuje dva extrémy v těsné blízkosti. Příznivé však je, že i zde ZAF tvoří téměř celou plochu grafu. • Skupina 2 (CC, NCC) - Míry založené na výpočtu korelačního koeficientu. Na první pohled je zřejmé, že pro tento typ obrazových dat vyjdou obě funkce, jak CC tak NCC, shodně. Výhodou normalizace je totiž zbavení se přímé úměry mezi množstvím bodů v obraze a maximální hodnotou podobnostní funkce a především možnost registrovat obrazy s neuniformním pozadím (většinou už nenulové pozadí představuje pro CC problém) a snížení vlivu objektů na hranici obrazu. K výhodám normalizace a problémům spojeným s existencí objektů na hranici obrazu se ještě jednou vrátíme u funkcí počítaných přes Fourierovu transformaci. V tomto případě, kde se nám žádný z uvedených jevů nevyskytuje, však poskytují obě funkce hladký průběh s velkou ZAF. • Skupina 3 (MI, GMI) - Míry založené na vzájemné informaci. MI vykazuje velmi výrazný extrém v místě optimální registrace, zatímco v okolí působí graf na první pohled mírně plošším dojmem. Je to však pouze zdání způsobené měřítkem a při bližším prozkoumání je zřejmé, že ZAF tvoří většinu prostoru grafu, tedy nevysky­ tují se žádná lokální minima. Podobně je tomu u GMI, které však není tolik ploché v okrajových částech. Co je nepříjemné, jsou doby registrace a výpočtu podobnostní míry pro jednu transformaci, které jsou více než dvakrát větší než u ostatních funkcí. Je to dáno především nutností sestavení vzájemného histogramu obou obrazů, nut­ ného pro výpočet Shannonovy či Rényiho vzájemné entropie. Pro urychlení výpočtu

3.1. REGISTRACE

2D OBRAZŮ

19

existuje několik postupů pro odhad vzájemné entropie, či hustoty pravděpodobnosti, které ve své práci uvádí Čapek [1]. Jejich nevýhodou je bohužel značné zhoršení vlast­ ností průběhu obou funkcí. V mé práci bude diskutováno pouze urychlení založené na zmenšení velikosti vzájemného histogramu, tedy na snížení počtu kvantizačních úrovní. • Skupina 4 (RIU, P I U ) - Míry navržené Woodsem. Ve shodě s literaturou je patrná (v tomto případě zvláště u PIU) menší velikost ZAF, než u ostatních funkcí, což Woodsovy algoritmy mírně znevýhodňuje. Horší jsou také doby výpočtu, které se nachází někde mezi prvními dvěma skupinami a skupinou třetí. RIU je časově méně náročné než PIU a řadí se spíše k prvním dvěma skupinám. PIU, naproti tomu, se svým časem blíží mírám založeným na vzájemné informaci. • Skupina 5 (F-CC, F-NCC, F-PCC) - Míry založené na integrální transformaci. Do této kategorie spadají míry počítající funkce ze skupiny 2 přes Fourierovu trans­ formaci (tedy F-CC a F-NCC) a fázová korelace (F-PCC). Jak bylo uvedeno v úvodu kapitoly, v obrázku 3.3 jsou zachyceny grafy pro varianty F-CC a F-NCC implicitně uvažující periodicitu obrázků. Výsledky se tak mohou odlišovat od standardních funkcí CC a NCC. F-NCC bez periodicity (tedy počítaná na obrazu s doplněním nulových intenzit až do dvojnásobku rozměrů obrazu) není v grafech uvedena, pro­ tože její výstup by byl naprosto shodný s NCC počítanou klasickou cestou. F-CC bez periodicity není uvedena ze 2 důvodů. Ten první je, že podobně jako u NCC, by její výstup byl shodný s klasickou CC. Ten druhý, závažnější, je, že využitelnost CC bez periodicity je omezena na velmi malou skupinu obrazů (objekty obklopené uniformním, nejlépe nulovým, pozadím). Pro jiné typy obrazů je už z definice CC, jako takové, patrné, že globálního extrému bude dosaženo pro nulovou transformaci, neboť v průniku obrazů bude nejvíce bodů a tedy suma součinů vyjde největší. Zdání plochého průběhu funkcí F-CC a F-NCC je opět dáno měřítkem grafů a jeho nato­ čením. Zcela jiným případem je F-PCC, která pro tyto registrované obrazy dává výborné výsledky: v místě registrace se skutečně objevuje delta funkce s výraznou amplitudou, zatímco pro všechny ostatní uvažované transformace je míra podobnosti téměř nulová. Skutečně výjimečné jsou pak doby výpočtu. Hodnota 0 v tabulce vy­ jadřuje dobu výpočtu (všechny testované transformace) výrazně nižší než 1 sekunda, a to pro všechny podobnostní míry z této skupiny. Shrneme-li výsledky, shledáme všechny podobnostní míry jako použitelné pro tento typ obrazových dat (výrazný objekt nebo objekty obklopené uniformním pozadím). Jedno­ značným vítězem jsou potom míry počítané integrální transformací, díky bezkonkurenčně nízkým dobám registrace. O něco horší jsou potom funkce z první a druhé skupiny. Algo­ ritmy ze třetí a čtvrté skupiny lze pak označit za méně vhodné. U skupiny 3 je limitujícím faktorem doba výpočtu. U skupiny 4 navíc relativně menší ZAF. Výhody těchto funkcí se však projeví při registraci „obtížnějších" dat. 3.1.2

R e g i s t r a c e reálných dat

Po otestování a srovnání základního chování podobnostních mír na uměle transformova­ ných datech se přesuneme už ke zcela reálným případům. Cílem nyní bude sledovat jak se

KAPITOLA

20

3. VÝBĚR PODOBNOSTNÍ

MÍRY

změna dat projeví na grafech funkcí. Jako registrovaný pár číslo 1, byly vybrány obrazy zobrazené na obrázku 3.4. Obrazy pocházejí ze dvou různých snímání téže buňky, která byla mezi prvním a druhým snímáním upravena v laboratoři.

Obrázek 3.4: Tatáž buňka po opakovaném snímání

Na první pohled je vidět, že oba dva obrazy znázorňují ten stejný objekt, nicméně tentokrát jsou zde již patrné faktory, které by mohly znesnadnit registraci. Jedná se o roz­ díly v histogramech obou obrazových souborů. Jak ukazují obrázky 3.5 a 3.6, tak zatímco obraz vlevo (A) využívá téměř celého spektra intenzit z intervalu (0,255), spektrum ob­ razu napravo (B) je spíše v oblasti nižších hodnot jasu. Globální charakter histogramu je zachován a jeví se proto jako vhodné vyrovnání histogramů obou obrazů. Vyrovnání bylo provedeno lineární interpolací intenzit v obou obrazech tak, aby bylo využito celého spektra.

n

I

Obrázek 3.5: Histogram obrazu A z obr. 3.4

n

i

Obrázek 3.6: Histogram obrazu B z obr. 3.4

Po vyrovnání histogramů byly na obrazech, stejně jako v předešlém případě, spočteny všechny podobnostní míry. Výsledné průběhy ukazují grafy na obrázcích 3.7 a 3.8:

3.1. REGISTRACE

2D OBRAZŮ

21

I, TY [pixel] = -4

Sum of squared differ

rem = 1 84179e^i07 pro TX [pixel] = 7. TY [pixel] =

Sum of absolute values of differ.

= 17.4169 pro TX [pixel] = 6. TY [pixel] = -2

s correlation; extrem = -1.76981 e-HľiOS pro TX [pixel] = 7. TY [pixel] =

: -0.926586 pro TX [pixel] = 7, TY [pixel] = -2

Generalized mutual information; extrem = -1.36035 pro TX [pixel] = 11, TY [pixel] = -B

Woods ÍRatio-lmege Uniformity); extrem = 0.536535 pro TX [pixel] = 8, TY [pixel] = -5

Obrázek 3.7: Průběh extrémů funkcí při registraci obrazů z obr. 3.4 [část 1]

KAPITOLA

22

3. VÝBĚR PODOBNOSTNÍ

MÍRY

Woods (Partitioned Intensity Uniformity); extrem = 0 364211 pro TX [pixel] = 7, TY [r.

í

0.2 4 [

Fourier ľCC; extrem = 0.92905 pro TX [pixel] = 7, TY [pixel] = -2

oK|pixel]=14,TY

Obrázek 3.8: Průběh extrémů funkcí při registraci obrazů z obr. 3.4 [část 2]

Lze položit otázku, zda bylo nutné před registrací provést vyrovnání histogramů. Od­ pověď je složitější - závisí na zvolené podobnostní funkci. Pro názornost byla vyzkoušena registrace bez vyrovnání histogramů. U většiny funkcí došlo pouze ke změnám zanedba­ telného charakteru (změny pozice globálního extrému do 1 pixelu). Problémy se projevily nejvíce u podobnostních funkcí SSC a SAVD. Jak ukazuje obrázek 3.9, je pro ně v tomto případě vyrovnání histogramů zásadní. Nedojde-li k vyrovnání, můžeme vidět, že průběh SSC je pro registraci nepoužitelný a pozice globálního extrému se nalézá zcela mimo oblast, kde by měla být. U SAVD rovněž pozice extrému neodpovídá skutečnosti. Stochastic sign change; f

n = -463 pra TX [pixel] = -27, TY [f

Sum of absolute

Obrázek 3.9: Průběh extrémů funkcí SSC a SAVD pro obrazy 3.4 bez vyrovnání histo­ gramů

3.1. REGISTRACE

2D OBRAZŮ

23

Naměřené hodnoty a doby výpočtu shrnuje tabulka 3.2: Tabulka 3.2: Doba výpočtu a pozice extrému při registraci obrazů z obr. 3.4 Podobnostní míra Doba registrace Pozice extrému SSC

9s

(11, -4)

SAVD

8s

(6, -2)

SSD

8s

(7, -2)

CC

8s

(7, -2)

NCC

10 s

(7, -2)

MI

18 s

(8, -4)

GMI

15 s

(11, -6)

RIU

9s

(8, -5)

PIU

16 s

(7, -4)

F-CC (s periodicitou)

0s

(7, -2)

F-NCC (s periodicitou)

0s

(7, -2)

F-NCC (bez periodicity)

0s

(7, -2)

F-PCC

0s

(14, 13)

SSC (bez vyrovnání histogramů)

8s

(-27, 6)

SAVD (bez vyrovnání histogramů)

8s

(0,0)

Tentokrát se už pozice globálního extrému jednotlivých funkcí od sebe odlišují. Pokud nebudeme uvažovat zjevně chybné hodnoty dostáváme rozpětí 5 pixelů ve směru osy x a 4 pixely ve směru osy y, což jsou poměrně vysoká čísla. Bohužel nelze objektivně určit, kde přesně by mělo registrační optimum ležet. Jedinou možností v tomto případě je subjek­ tivní hodnocení. Na základě porovnání obrazu s registrovaným protějškem (viz příloha) se mi jeví jako nejlepší transformace navržená Woodsovým algoritmem Partitoned Intensity Uniformity (PIU), tedy posun obrazu B o vektor (7, -4). Registrace s použitím SSC (i přes vyrovnání histogramů) a poněkud překvapivě i GMI, lze na základě subjektivního porov­ nání jednoznačně označit za neuspokojivé. Rovněž výsledek PCC je zjevně chybný. Proč došlo k chybné registraci se pokusím rozebrat níže. Další hodnocení je rozděleno podle skupin, jako v případě uměle transformovaných obrazů: • Skupina 1 (SSC, SAVD, SSD) - Míry založené na rozdílech intenzit. SSC ani přes vyrovnání histogramů nedokázalo registrovat obrazy správně. Tuto míru lze tedy již nyní označit za nevhodnou. Globální extrém funkcí SAVD i SSD naopak i nyní odpovídá realitě, přičemž o něco lepší vlastnosti a hladší průběh opět vykazuje SSD. Vzhledem k poměrně nízké výpočetní náročnosti lze tyto podobnostní míry pro obrazy podobného typu doporučit. • Skupina 2 (CC, NCC) - Míry založené na výpočtu korelačního koeficientu. Mírně

24

KAPITOLA

3. VÝBĚR PODOBNOSTNÍ

MÍRY

lepší výsledek tentokrát podala funkce CC (je rychlejší). Nicméně jak CC tak NCC se ukazuje jako jednoznačně vhodnější počítat přes integrální transformaci. Skupina 3 (MI, GMI) - Míry založené na vzájemné informaci. MI registrovala obrazy správně. Jeho nevýhodou je však nejvyšší výpočetní doba. Na druhé straně GMI, jak už bylo zmíněno, vyhodnotila jako nejlepší transformaci, která po vizuální kontrole neodpovídá bodu registrace obou obrazů. Navíc podobně jako MI potřebuje pro svůj výpočet více výpočetního času než ostatní podobnostní míry. • Skupina 4 (RIU, P I U ) - Míry navržené Woodsem. Pro tento konkrétní obraz

dávají obě míry jak RIU, tak PIU, velmi dobré výsledky. Co se týká přesnosti vy­ hodnocení optimální registrace PIU se umístila na prvním místě, RIU subjektivně na třetím místě spolu s S SD, CC a NCC za vzájemnou informací (MI). Je ovšem třeba opět zmínit vyšší výpočetní doby (zejména PIU) a ve srovnání s ostatními podobnostními mírami nižší ZAF. Tedy teoreticky těžší optimalizaci při složitějších registracích. • Skupina 5 (F-CC, F-NCC, F-PCC) - Míry založené na integrální transformaci. F-CC a F-NCC (s periodicitou i bez ní) dopadly podle očekávání stejně jako je­ jich protějšky počítané standardní cestou, znovu pouze s velmi výraznou časovou úsporou. Pokud bychom tedy uvažovali o registraci pomocí korelačního koeficientu či normalizovaného korelačního koeficientu je rozumné použít varianty počítané přes Fourierovu transformaci. Podmínkou je ovšem registrace pouze posunů, protože pro registraci změny měřítka či rotací nelze Fourierovu transformaci v této podobě po­ užít. Zbývá ještě rozebrat podobnostní míru PCC. Podíváme-li se na graf průběhu vidíme, že PCC je pro tento obrázek nevhodné. Otázkou je proč. Odpověď není v tomto případě snadná. Z definice PCC vyplývá, zeje nutné (narozdíl od CC a NCC) počítat ji nad obrázkem s uvažovanou periodicitou. Z toho ovšem dostáváme, že obrázky u nichž se vyskytují objekty na hranici nelze pomocí PCC registrovat. Poskládáme li totiž obrazy periodicky vedle sebe vytvoří se nám „zuby" v intenzitě na hranicích a po převodu do frekvenční oblasti Fourierovou transformací se tím objeví nechtěné frekvence, které znemožní obrazy fázovou korelací zregistrovat. Na registrovaných obrazech na obrázku 3.4 nejsou patrné žádné objekty na hranici, provedeme-li ale jejich zesvětlení zjistíme, že v levém dolním rohu vystupuje z pozadí jakási „mlha", která vede od buňky až na hranici obrazu. Domnívám se, že právě ona by mohla být důvodem rozdílných výsledků a špatného výsledku PCC. Je to však pouze domněnka. Jisté je, že PCC nelze pro tyto obrazy spolehlivě použít. Shrnutím výše uvedeného dostáváme první množinu podobnostních mír použitelných na dodaných datech. S výjimkou SSC a PCC lze použít libovolnou ze zbývajících funkcí pro nalezení vyhovující registrace. SSC a PCC není možné v žádném případě doporučit. Naopak nejvýhodněji se jeví normalizovaný korelační koeficient počítaný přes integrální transformaci, tedy funkce v grafech a tabulkách označovaná jako F-NCC. Jako další 2 obrazy vhodné k prověření vlastností jednotlivých funkcí byla vybrána dvojice obrazů s výrazným centrálním objektem, avšak velkým množstvím menších ob­ jektů ležících na hranici. Vizuálním srovnáním je navíc patrné, že oba obrazy se od sebe

3.1. REGISTRACE

2D OBRAZŮ

25

liší více, než jen posunem. Mění se zde jak rozestupy jednotlivých objektů (buněk), tak jejich tvar a částečně i pozice. Cílem bude co nejlépe registrovat centrální objekt, který pravděpodobně bude biologa zkoumajícího obrazy nejvíce zajímat. Obrazy si můžeme pro­ hlédnout na obrázku 3.10

Obrázek 3.10: Obrazy buněk s objekty na okraji

Co bude rovněž zajímavé sledovat je změna doby výpočtu registrací. Zatímco předchozí 2 registrované soubory obsahovaly obrazy o velikosti 150 x 150 pixelů, tak v tomto případě registrujeme dvojici, v níž každý z obrazů má velikost 365 x 371 pixelů, tedy celkem 135 415 obrazových bodů. Připomínám, že zatím všechny registrace probíhaly na obrazech v plném rozlišení a bude tomu tak i nyní. V předchozím případě nám pomohlo vyrovnání histogramů zlepšit výsledky registrace. Pokud se o něco podobného pokusíme i nyní nebudeme bohužel úspěšní. Globální charakter histogramu obou obrazů je totožný, pouze histogram obrazu vpravo (B) vykazuje jisté známky degradace, jak ukazuje obrázek 3.11.

Ik

I

Obrázek 3.11: Histogram obrazu B z obr. 3.10 Pokusil jsem se to napravit aplikací nízkofrekvenčního Gaussova filtru, což sice pomohlo opravit průběh histogramu, ale jak prokázaly mé testy, na výsledky registrace, které vidíme na obrázcích 3.12 a 3.13, to mělo pouze velmi minoritní vliv.

KAPITOLA

26

o TX [pixel] = -14, TY [pixel] = -2

Sum of absc

Sum ofetpjared difference

Normalized cross correlatí

3. VÝBĚR PODOBNOSTNÍ

elation; extrem = -1.22939e-HII09 pro TX [pixel] = 0. TY [pixel] = -10

= -0.815294 pro TX [pixel] = 10, TY

Generalized mutual information; extrem = -0.991982 pro TX [pixel] = 16, TY [pixel] = -10

Mutual information; extrem - -0.700466 pro TX [pixel] - 14, TY [pixel] - -14

Woods (Rat

34315 pro TX [pixel] = 8. TY [pixel] = -17

Obrázek 3.12: Průběh extrémů funkcí při registraci obrazů z obr. 3.10 [část 1]

MÍRY

3.1. REGISTRACE

2D OBRAZŮ

27

24129e^09 pro IX [pixel] = 6, TY [pixel]-

Woods (Partitioned Intensity Uniformity); extrem = 0.264513 pro TX [pixel] = 10, TY [pixel] = -17 0.5^

f

1.25-.

l045-

1 S

1.2-, 0.4-

• ^

1

= 035- •••""":"""" ^ "

I"' •Š 0.25-

I

0.2zi

8

Jß wr^Jr

1

""^^jj? •

!

•

•

~r^-^Z

"

"

'

"

^

^

^

115-, '•11.05,

^

' ' " ' ••'•'•' - ä Í E ~ r ~ = = i .

1 5

;

"~~xT

4o

Obrázek 3.13: Průběh extrémů funkcí při registraci obrazů z obr. 3.10 [část 2]

Tentokrát se podle očekávání pozice globálních extrémů jednotlivých funkcí odlišují již velmi výrazně, což nám umožní provést lepší selekci vhodných a nevhodných podob­ nostních funkcí. Naměřené doby registrací a pozice globálního extrému přehledně shrnuje tabulka 3.3. Pro porovnání nalezených transformací mezi sebou byla opět použita vizuální kontrola (viz příloha). Tentokrát z ní subjektivně jednoznačně nejlépe vyšly míry založené na vzájemné informaci (MI, GMI), kterým se podařilo nejlépe registrovat centrální ob­ jekt. Všechny ostatní míry s výjimkou SSC, CC a PCC dosáhly zhruba stejného výsledku, na kterém se ale projevuje vliv okolních objektů a registrace centrální buňky je viditelně mírně nedokonalá. Udělejme si ještě jednou hodnocení podle skupin: • Skupina 1 (SSC, SAVD, SSD) - Míry založené na rozdílech intenzit. Pro tuto dvojici registrovaných obrazů, podle předpokladů, opět nelze použít SSC. SSD a SAVD dávají výsledek částečně se blížící optimální registraci, přitom výsledek SAVD je subjektivně o trochu lepší. • Skupina 2 (CC, N C C ) - Míry založené na výpočtu korelačního koeficientu. Již bylo řečeno, že pro obrazy obsahující objekty na hranici není možné korelační koeficient (bez periodicity) použít, což se projevilo v tomto případě a CC nedokázala najít místo registrace daných obrazů. NCC podobné faktory neovlivňují. Její výsledek byl proto lepší, srovnatelný s SSD a SAVD z první skupiny, pouze výpočet byl delší. Ze zatím uvedených výsledků se může zdát, že funkce NCC není pro registrace vhodná,

28

KAPITOLA

3. VÝBĚR PODOBNOSTNÍ

MÍRY

Tabulka 3.3: Doba výpočtu a pozice extrému při registraci obrazů z obr. 3.10 Podobnostní míra Doba registrace Pozice globálního extrému

ssc

53 s

(-14, -2)

SAVD

49 s

(10, -18)

SSD

49 s

(9, -21)

CC

47 s

(0, -10)

NCC

66 s

(10, -18)

MI

85 s

(14, -14)

GMI

82 s

(16, -11)

RIU

54 s

(8, -17)

PIU

66 s

(10, -17)

F-CC (s periodicitou)

0s

(6, -18)

F-NCC (s periodicitou)

0s

(6, -18)

F-NCC (bez periodicity)

2s

(10, -18)

F-PCC

0s

(0,0)

neboť podává stejné výsledky jako SAVD a SSD, pouze v horších časech. Opak je však pravdou. Jak se pokusím rozebrat v části o vlivu degradace dat, NCC je oproti funkcím z první skupiny odolnější (robustnější) proti vlivu šumu či lineární transformace intenzit v obraze - např. NCC by dokázala registrovat předchozí obraz i bez vyrovnání histogramů. • Skupina 3 (MI, GMI) - Míry založené na vzájemné informaci. Tentokrát hodnotím výsledek MI a GMI, i přes drobné rozdíly v pozici extrému, jako nejlepší. Lépe z těchto dvou však dopadla jako vždy ML Několikrát zmiňována byla délka výpočtu. Ta je i zde zhruba dvakrát delší, než pro ostatní podobnostní míry. Je to však vyváženo kvalitou registrace (minimálně v případě MI). • Skupina 4 (RIU, PIU) - Míry navržené Woodsem. Tyto míry dosáhly podobných výsledků jako míry ze skupin 1 a 2. Projevil se zde ale zvláštní úkaz, a to významné zrychlení funkce PIU. Ta se dostala na stejnou dobu výpočtu jako NCC a částečně se i přiblížila druhé Woodsově podobnostní míře RIU. Důvod je prostý, pro rostoucí objem dat vstupuje do doby výpočtu stále větší mírou čas spotřebovaný na samotné transformace dat, které je nutné v každém kroku udělat před spuštěním výpočtu podobnostní míry. Dále je také důležité na čem závisí složitost výpočtu podobnostní míry. Například pro RIU jsou nutné minimálně dva průchody celým obrazem v prvním je spočítán podílový obraz a vypočtena průměrná hodnota, ve druhém je potom vypočtena standardní odchylka. U PIU však stačí jen jeden průchod obrazem, ve kterém je vytvořen vzájemný histogram. Další výpočet už závisí jen na velikosti tohoto histogramu. Uvažujeme-li ale šedotonní obrázky o 256 úrovních intenzity, pak je velikost vzájemného histogramu konstantní a velikost obrazů tedy už nemá

3.1. REGISTRACE

2D OBRAZŮ

29

na tuto část výpočtu vliv. Limitně ji lze tedy zanedbat a dá se předpokládat, že pro ještě větší objem dat bude PIU rychlejší než RIU i NCC. • Skupina 5 (F-CC, F-NCC, F-PCC) - Míry založené na integrální transformaci. V této skupině k žádnému překvapení nedošlo. Již po několikáté si můžeme všimnout, že budeme-li uvažovat periodicitu je výsledek F-CC i F-NCC naprosto shodný. Bez periodicity potom dostáváme u F-NCC stejný výsledek jako pro NCC, ale 33 x (!) nižší dobu registrace. PCC znovu nedokázalo obrazy registrovat. Zbývá se ještě zamyslet co způsobilo tak velké rozdíly v pozicích globálních extrémů. Jsou to jednoznačně okolní buňky, které, jak je vidět na úvodním obrázku, v žádném případě nevykazují mezi těmito dvěma obrazy pouze posun. Tím nám bohužel výsledky registrace velmi zkreslují a dokonce znemožňují použití některých jinak dobrých metod (PCC). Dle mého názoru je velice obtížné takovéto obrazy správně zregistrovat. Nejvhod­ nější by bylo nějakým procesem na bázi segmentace okolní buňky odstranit, ponechat pouze centrální objekt zájmu a ten poté registrovat. Nicméně jako nejlepší byla vybrána transformace navržená MI, tedy posun obrazu B o vektor (14, -14). Výslednou dvojici (obraz A a registrovaný obraz B) ukazuje obrá­ zek 3.14. Šedá místa vyznačují pixely, pro které z důvodu posunu neexistuje informace o intenzitě. Z obrázku vidíme, že slícování centrální buňky je poměrně dobré, avšak okolní buňky (zejména pak ty na levé straně) si stále pozičně neodpovídají.

Obrázek 3.14: Výsledek registrace obrazů z obr. 3.10 pomocí MI Z obrázku se zdá, že k lepšímu celkovému slícování by mohla pomoci změna měřítka. Pokusil jsem se proto ještě o registraci změny měřítka (s použitím podobnostní míry MI). Aby byl prostor parametrů stále pouze dvourozměrný zafixoval jsem posun na dříve vyhledané pozici, tedy (14, -14) a registroval pouze změnu měřítka v rozsahu 0,8 až 1,2 pro obě osy. K tomu jsem použil znovu úplné vyhledávání s krokem 0,05. Výsledný průběh funkce MI pak ukazuje obrázek 3.15 a zregistrovanou dvojici obrázek 3.16.

30

KAPITOLA

3. VÝBĚR PODOBNOSTNÍ

MÍRY

Mutual information; extrem = -0.769447 pro SX = 1.045, SY = 1.07

1.2

1-3

Obrázek 3.15: Průběh extrému MI při registraci změny měřítka obrazů z obr. 3.10

Obrázek 3.16: Výsledek registrace změny měřítka u obrazů z obr. 3.10 pomocí MI

Extrém této registrace nám udává (viz obrázek s průběhem 3.15) zvětšení v ose x o 4,5% a v ose y o 7%. Tímto se podařilo zvýšit hodnotu vzájemné informace z před­ chozích 0,700466 na 0,769447, což lze považovat za úspěch. Registrace změny měřítka má tedy v tomto případě smysl. Také je z obrázků vidět, že se nám podařilo dosáhnout mírně lepší registrace okolních buněk. Stále jsou však viditelné rozdíly, k jejichž smazání by bylo dle mého názoru potřeba minimálně použití nerigidní transformace, spíše však nejde

3.1. REGISTRACE

2D OBRAZŮ

31

v této podobě obraz zregistrovat tak, aby došlo ke slícování všech buněk. Teoreticky lze dosáhnout ještě lepší registrace, než té zde uvedené, tím že ponecháme volné i parametry posunu a nebudeme je fixovat. Poté už ale dostáváme příliš velký (čtyřdimenzionální) pro­ stor transformací a nelze tak použít pro nalezení optima úplné vyhledávání. Proto tuto registraci provedu až v následující kapitole o optimalizacích. Další dodaná data mají již podobný charakter, a proto srovnávat v tomto místě podob­ nostní míry na dalších obrazech nemá smysl. Další ukázky grafů průběhů podobnostních funkcí jsou uvedeny v příloze této práce a na doprovodném CD. 3.1.3

V l i v degradace

Při opakovaném snímání buněk či jiných objektů nebo při pořizování jednotlivých ob­ razů pomocí různých zobrazovacích modalit (zařízení) se bohužel mohou projevit jisté známky degradace informace obsažené v obrazech. Jeví se proto jako velmi vhodné zkou­ mat jakým způsobem reagují jednotlivé funkce při registraci obrazů postižených některým z druhů degradace. Tedy jakousi robustnost podobnostní funkce. Tato robustnost nám pak spolu s dalšími vlastnostmi jako výpočetní náročností, velikostí ZAF, apod. umožní vybrat nejvhodnější funkci. Mezi základní druhy degradace patří: • Šum (v tomto případě nejčastěji Poissonův) • Lineární transformace obrazových intenzit • Nelineární transformace obrazových intenzit Jak šum, tak lineární transformace obrazových intenzit se už v práci vyskytly. Šum se inherentně vyskytuje ve většině dodaných snímků z mikroskopu, a s lineární transformací obrazových intenzit jsme se setkali u prvního registrovaného páru obrazů (nepočítáme-li uměle transformovaná data), kde bylo třeba pro správnou funkčnost některých podobnost­ ních funkcí provést vyrovnání histogramů. S nelineární transformací obrazových intenzit jsme se zatím nesetkali. Ta se vyskytuje zejména při registraci dat z různých zobrazovacích modalit (např. PET-MR), což není náš případ, ale také, jak uvádí Čapek [1], při regis­ traci obrazů barvených různými fluoresnečními činidly. Čapek se ve své práci důkladně zabývá rozborem jednotlivých degradací. Není účelem této práce jeho experimenty opa­ kovat. Uvádím tak pouze orientační tabulku 3.4, ve které jsou uvedené podobnostní míry odolné proti jednotlivým druhům degradace. Protože v Čapkově práci se neověřuje možnost použít pro nelineární transformaci i GMI, provedl jsem alespoň krátký test, ve kterém jsem nelineárně upravil intenzity v ob­ raze B (pravém) z obr. 3.1 a zkusil data registrovat. Výsledek, který je možné vidět na obrázku 3.17, dokazuje, že pro registraci dat obsahujících nelineární transformaci intenzit je možné použít pouze 3 podobnostní míry uvedené v tabulce. Obrázek zobrazuje GMI a PIU, které správně určily registrační posun (-15, -7) a NCC a SSD, které byly (mimo jiné) neúspěšné. MI, které v obrázku uvedeno není, také dokázalo data registrovat správně. Zde se projevuje pravá síla MI (a částečně i GMI), která spočívá v jejich velké robust­ nosti a vyvažuje tak, drobný nedostatek v podobě delší doby výpočtu. V tabulce si můžeme

32

KAPITOLA

3. VÝBĚR PODOBNOSTNÍ

MÍRY

Tabulka 3.4: Funkce vhodné pro jednotlivé typy degradace dat Typ degradace Použitelné podobnostní míry Bez degradace

SSD, NCC, MI, GMI, WOODS, ...

Aditivní šum

SSC, SAVD, NCC, MI, GMI

Lineární transformace intenzit

NCC, MI, GMI, WOODS

Nelineární transformace intenzit

MI, GMI, WOODS (PIU)

Sum of squared differences; extrem = 4.46753e+007 pro TX |pixel] = 40. TY [pixel] = -40

Normalized cross correlation; extrem =-0.569194 pro TX |pixel] = 40, TY [pixel] = -40

Generalized mutual information; extrem = -2 919B7 pro TX [pixel] = -15, TV [pixel] = -7

Woods (Partitioned Intensity Uniformity); extrem = 0 pro TX [pixel| = -15, TY [pixel] = -7

TY [pixel]

-40

TX [pixel]

TY [pixel]

Obrázek 3.17: Registrace dat obsahujících nelineární transformaci intenzit

všimnout ještě jedné zajímavé informace a to, že přestože algoritmy vzájemné informace by měly být z definice náchylné na výskyt šumu (šum snižuje množství informace), je jejich použití i v těchto případech (pokud není šumu extrémní množství) vhodné. 3.1.4

Citlivost na p a r a m e t r y

U většiny podobnostních funkcí nelze, nebudeme-li uvažovat převzorkování obrazu na menší velikost, snížit dobu jejich výpočtu. Je to dáno už z jejich definice, protože do vý­ počtu u nich, kromě velikosti obrazu, nevstupují už žádné dodatečné parametry. Výjimku v tomto případě tvoří podobnostní míry MI, GMI a Woodsův algoritmus PIU. Ty totiž závisí také na počtu úrovní intenzity (kvant iz ač nich úrovní - k.ú.). Tento počet ovlivňuje velikost jednotlivých histogramů obrazů a jejich vzájemného histogramu a tím pádem i dobu výpočtu těchto funkcí. Doposud byly všechny registrace prováděny se standardním

3.1. REGISTRACE

2D OBRAZŮ

33

počtem k.ú pro šedotonní obrazy, tedy k.ú. = 256. Velikost vzájemného histogramu je v tomto případě 2562 = 65 536. Bylo by zajímavé pozorovat jakým způsobem se na době výpočtu a přesnosti registrace projeví snížení k.ú. Provedl jsem proto 2 testy. V prvním jsem znovu registroval relativně malé obrazy z obrázku 3.1 (tentokrát však bez přidaného šumu), přičemž jsem postupně snižoval počet kvantizačních úrovní. Podobně jsem postu­ poval i v případě větších obrazů z obrázku 3.10. Tabulka 3.5 uvádí výsledky registrací v prvním případě: Tabulka 3.5: Vliv počtu kvantizačních úrovní na doby výpočtu a pozici extrému při registraci obrazů z obr. 3.1 Podobnostní míra Doba registrace Pozice extrému Hodnota extrému 128 k.ú. MI

14 s

(-15, -7)

-3,91302

GMI

12 s

(-15, -7)

-3,32838

PIU

10 s

(-15, -7)

0

64 k.ú. MI

12 s

(-15, -7)

-3,24776

GMI

12 s

(-15, -7)

-2,6566

PIU

9s

(-15, -7)

0

32 k.ú. MI

12 s

(-15, -7)

-2,58546

GMI

12 s

(-15, -7)

-2,02129

PIU

9s

(-15, -7)

0

Porovnáme-li tuto tabulku s tabulkou 3.1, je zřejmé, že největší úspory v čase ve všech případech dosáhneme při přechodu z 256 k.ú. na 128 k.ú. Při dalším poklesu k.ú. není již časová úspora nijak zásadní. Ve všech případech registrace dopadla správně a byla nalezena transformace (-15, -7). Je dobré poznamenat, že obrazy byly registrovány správně i pro hodnotu 16 k.ú. a dokonce i v případě 8 k.ú.! Co však rozhodně stojí za povšimnutí je fakt, že při přechodu ze 256 na 128 k.ú. se ani v jednom případě nezměnila hodnota extrému funkce. Tato změna nám tedy umožňuje znatelné zrychlení výpočtu bez ztráty přesnosti. Při přechodu na 64 k.ú. a níže můžeme již dobře pozorovat pokles míry podobnosti, tedy extrém v grafu se stává méně výrazným. Jiným způsobem se ovšem pokles počtu k.ú. na grafech neprojevil, proto je zde neuvádím. Na grafech bychom mohli pozorovat změny zřejmě jen při současném převzorkování obrazů na nižší velikost, kdy by nám snížení počtu k.ú. způsobilo vyhlazení grafu (ověřeno pokusem). Ve své podstatě se jedná o odstranění vysokých frekvencí z obrazu. Jak dopadla registrace v druhém případě ukazuje tabulka 3.6. Výsledky kopírují ty předchozí. Srovnáme-li s tabulkou 3.3 znovu vidíme, že největší časové úspory dosáhneme při přechodu ze 256 k.ú. na 128 k.ú. Nicméně tato časová úspora je bohužel nižší, než by

KAPITOLA

34

3. VÝBĚR PODOBNOSTNÍ

MÍRY

Tabulka 3.6: Vliv počtu kvantizačních úrovní na doby výpočtu a pozici extrému při registraci obrazů z obr. 3.10 Podobnostní míra Doba registrace Pozice extrému Hodnota extrému 128 k.ú. MI

82 s

(14, -14)

-0,700466

GMI

81 s

(16, -10)

-0,991982

PIU

59 s

(10, -17)

0,260692

64 k.ú. MI

81 s

(14, -13)

-0,673046

GMI

81 s

(16, -10)

-0,981891

PIU

58 s

(10, -17)

0,256059

32 k.ú. MI

80 s

(13, -15)

-0,656482

GMI

80 s

(16, -10)

-0,936219

PIU

57 s

(10, -17)

0,25862

16 k.ú. MI

80 s

(13, -15)

-0,616348

GMI

80 s

(15, -10)

-0,85276

PIU

57 s

(10, -17)

0,252348

se dalo očekávat (pohybuje se zhruba okolo 5 sekund). Znovu se také opakuje výsledek, že při přechodu na 128 k.ú. se nemění hodnota extrému ani jeho pozice. Při dalším snižování už ovšem tentokrát pozorujeme jak pokles hodnoty extrému, tak mírné změny jeho pozice. Nejvíce je to patrné v případě MI. Shrnu-li naměřené hodnoty, tak počet 128 k.ú. se mi jeví jako vhodný kompromis. Při jeho použití dosáhneme největší časové úspory, při současném zachování přesnosti výpočtu. V teoretické části práce bylo zmíněno, že na registraci může mít vliv i stupeň použité interpolace a to jak při převzorkování obrazů na menší velikost (pyramidová metoda), tak při výpočtu jednotlivých transformací. Provedl jsem sérii testů z jejichž výsledků vyplývá, že vliv stupně interpolace na výsledek registrace skutečně existuje. Je však velmi malý. Jelikož vyšší stupně interpolace brzdí výpočet, doporučuji ve většině případů ponechat nejnižší stupeň — tedy interpolaci typu nejbližší soused. Pouze při převzorkování obrazů na nižší velikost je vhodné použít interpolaci lineární, bilineární či kubickou. Ukazuje se, že vyšší by neměla smysl. Vhodné je kombinovat ji se současným snížením počtu kvantizač­ ních úrovní, pro další urychlení výpočtu. Stupeň interpolace pro transformace doporučuji ponechat 0. Vyšší by našel opodstatnění pouze při výraznějších změnách měřítka či rota­ cích.

3.2. REGISTRACE

3.2

3D OBRAZŮ

35

Registrace 3D obrazů

Doposud byly všechny registrované obrazy v této práci dvourozměrné. Bohužel data která bude třeba v praxi registrovat budou většinou třírozměrná a bude je třeba registrovat ve směru všech tří os. Při přechodu na registraci 3D obrazů však není třeba provádět velké množství dalších srovnání. Všechny dosud zjištěné výsledky zůstávají v platnosti, pouze datové objemy se zvětšují, stejně jako prostor transformací. Více se projevují nedokonalosti některých podobnostních funkcí. Největším problémem při přechodu do třetího prostoru je růst datových objemů re­ gistrovaných obrazů a tím i značný růst výpočetní doby. Navíc nám přibývá třetí rozměr (uvažujeme-li pouze posuny) do prostoru transformací — tedy další zpomalení, protože je třeba prohledávat více kombinací. Pro názornou ukázku co se stane při přechodu z 2D do 3D jsem se pokusil registrovat obrazy z obrázku 3.4, pouze místo autofocus projekcí byly použity plné 3D obrazy. Plná velikost těchto obrazů je 159 x 153 x 60 voxelů. Smut­ ným zjištěním je, že ani při velmi znatelném omezení prostoru transformací nelze takovýto datový objem registrovat metodou úplného vyhledávání. Jednoduchým výpočtem dostá­ váme, že zachováme-li zhruba počet testovaných transformací, bude nám jedna registrace trvat okolo deseti minut. Proto jsem v tomto případě upustil od úplného vyhledávání a použil vyhledávání re­ kurzivní, které sice nezaručuje nalezení globálního extrému, ale umožnilo mi dostat se na podobné doby výpočtu jako v případě 2D obrazu. F-CC, F-NCC a F-PCC jsou samozřejmě i nadále počítány integrální transformací, tedy stejně jako doposud. Výsledky registrací shrnuje tabulka 3.7. Před výpočtem jsem opět provedl vyrovnání histogramů obou obrazů. Tabulka 3.7: Doba výpočtu a pozice extrému při registraci 3D obrazu Podobnostní míra Doba registrace Pozice globálního extrému

ssc

12 s

(9, -3, 0)

SAVD

13 s

(7, -2, -1)

SSD

13 s

(7, -2, -1)

CC

12 s

(5, -1, 0)

NCC

17 s

(7, -2, -1)

MI

21 s

(7, -3, -1)

GMI

21 s

(7, -4, -4)

RIU

17 s

(6, -2, -1)

PIU

16 s

(7, -3, -1)

F-CC (s periodicitou)

1s

(7, -2, -1)

F-NCC (s periodicitou)

2s

(7, -2, -1)

F-NCC (bez periodicity)

21 s

(7, -2, -1)

F-PCC

1s

(0, 0, 1)

Z tabulky je patrné, že podle očekávání nedošlo k žádným velkým překvapením. Zde je

KAPITOLA

36

3. VÝBĚR PODOBNOSTNÍ

MÍRY

podrobnější hodnocení: • Skupina 1 (SSC, SAVD, SSD) - Míry založené na rozdílech intenzit. Stejně jako ve 2D případě SSC nedokázalo obrazy správně registrovat. SAVD a SSD byly naopak spolehlivé. Je to však dáno tím, že tato dvojice obrazů je jednou z jednodušších. • Skupina 2 ( C C , N C C ) - Míry založené na výpočtu korelačního koeficientu. Zde už se výsledek od 2D případu mírně odlišuje. Zatímco CC bylo pro 2D obraz úspěšné, tak přidáním třetího rozměru se od optimální registrace vzdaluje. Tato podobnostní míra zde není rozhodně vhodná. NCC bylo úspěšnější a s registrací nemělo problémy. • Skupina 3 (MI, GMI) - Míry založené na vzájemné informaci. MI, jako zatím vždy, našlo optimální registraci spolehlivě. Opět je pouze mírně delší doba výpočtu. GMI naproti tomu tentokrát doporučit nelze, protože podobně jako u autofocus projekce se jeho výsledek od ostatních funkcí relativně zásadně odlišuje (o 3 pixely v ose z). Co je příčinou těchto problémů je otázkou. • Skupina 4 (RIU, P I U ) - Míry navržené Woodsem. PIU opět nelze nic vytknout. RIU je mírně horší - jak časově, tak v přesnosti. • Skupina 5 ( F - C C , F - N C C , F - P C C ) - Míry založené na integrální transformaci. S výjimkou PCC (u kterého se to čekalo) byly všechny tři funkce spolehlivé. Pouze u F-NCC bez periodicity se začíná značně projevovat nutnost registrovat 8 x (jedná se o 3D obraz) větší objem dat a stala se v tomto testu nejpomalejší ze všech. Ovšem je třeba uvážit, že v jejím případě, narozdíl třeba od klasické NCC, máme díky integrální transformaci zaručeno nalezení globálního minima. Míry bez periodicity pak byly opět nejrychlejší a i tento objem dat byly schopny registrovat pod 2 sekundy.

Founsr NCC; extrem = 0.322043 pra TY [pixel] = -2, TZ [pixel] = -1

TZ [pixel]

TY [pixel]

Obrázek 3.18: Registrace 3D obrazu pomocí F-NCC Pro zajímavost uvádím ještě průběh podobnostní míry F-NCC na těchto obrazech v osách y a z se zafixovaným posunem v ose x o 7 pixelů (obrázek 3.18). Uvádět další

3.3. DOPORUČENÉ PODOBNOSTNÍ

MÍRY

37

grafické průběhy (neliší se příliš od 2D případu) či dále rozebírat registraci 3D obrazů na tomto místě nemá smysl. Důležitější je se zde zaměřit na problém doby výpočtu a velikosti prostoru transformací. Tím se zabývá následující kapitola o výběru optimalizační metody.

3.3

Doporučené podobnostní míry

Tuto kapitolu bych rád završil doporučením podobnostní míry pro konkrétní typy dat. Shrneme-li výsledky z registrací z předchozích podkapitol, zjistíme, že ve všech případech našla dobře registrační optimum pouze funkce ML Navíc ta vychází nejlépe i z kapitol o degradacích. Bude-li tedy vyžadována co největší robustnost je volba zřejmá. Její jedinou nevýhodou je minimálně dvojnásobná doba výpočtu oproti ostatním funkcím. Proto pro registraci běžných dat, nevykazujících nelineární transformace intenzit apod., lze spíše doporučit podobnostní míru NCC, která také neměla v testech větší problémy. Pokud bude třeba registrovat pouze posuny, pak nejlépe její variantu počítanou přes Fourierovu transformaci, zajišťující extrémně rychlý výpočet. Velmi dobře ze srovnání vyšel také Woodsův algoritmus PIU, který je relativně rychlý a účinný. Navíc, podobně jako MI, je odolný proti degradaci dat. Zatím však v práci nebyla zmíněna jedna vlastnost obou Woodsových algoritmů a to jejich „nesymetričnost". Tedy PIU(A,B) se může lišit (velmi!) od PIU(B,A), podobně pro RIU. Použití PIU či RIU tak může být rizikem. Z dalších postřehů lze zmínit, že téměř stejné výsledky podávají funkce SSD a SAVD. Pouze SSD je mírně robustnější. Úplně stejné výsledky pak podaly funkce F-CC a F-NCC s uvažovanou periodicitou obrazů. Rozhodně nelze doporučit podobnostní míry SSC, CC a F-PCC.

Kapitola 4

Výběr optimalizační metody Kritickým faktorem při registraci velkých obrazů s uvažovaným velkým prostorem trans­ formací, je počet ohodnocení podobnostní funkce / nutný k nalezení globálního minima. Metody minimalizující tento počet na základě sofistikovaných algoritmů jsou označovány jako optimalizační metody. Optimalizačních metod existuje velké množství, výčet těch na které se zaměříme v této práci je uveden v teoretické části. Pro lepší hodnocení rozdělím ještě registrace podle toho, zda registrujeme pouze posuny, změny měřítka a otočení nebo i jejich kombinace. Vhodné postupy se totiž mohou pro tyto jednotlivé případy lišit. Rovněž už v této části nebudeme testovat všechny postupy na všech podobnostních mírách, ale zaměříme se pouze na ty, které budou s největší pravděpodobností v praxi použity — tedy NCC a MI.

4.1

Registrace posunů

Pro srovnání jednotlivých optimalizačních metod byly vybrány obrazy z obrázku 3.10 v plném rozlišení, což je 365 x 371 x 40 voxelů. Tyto dva obrazy byly postupně registrovány jednotlivými optimalizačními metodami, přičemž byly uvažovány následující transformace: posun v osách x a y od -30 do 30 voxelů a v ose z od -10 do 10 voxelů. Pro srovnání, extrém podobnostních funkcí ležel u NCC v bodě (8, -13, 3) a měl hodnotu -0.800439, u MI potom v bodě (11, -12, 3) s hodnotou -0.569735.

4.1.1

Úplné prohledávání

Úplné vyhledávání není optimalizační metodou v pravém slova smyslu. Na startu metody jsou nastaveny hranice prostoru parametrů a velikost kroku po kterém bude algoritmus postupovat (může se lišit pro jednotlivé dimenze). Budeme-li uvažovat nejjednodušší pří­ pad, kdy není třeba registrovat obrazy se sub-pixelovou přesností, bude nám stačit krok v podobě jednoho voxelů ve všech směrech. Po této volbě dostáváme velikost prostoru transformací rovnu číslu 72 000 (60 x 60 x 20). Dle měření trvá, v případě zvolených obrazů, testování jedné transformace podobnostní mírou NCC 0,4 sekundy. Při úplném prohledávání se tak dostáváme na hodnotu doby registrace rovnou 28800 sekundám. To je přesně 8 hodin. Pro MI je potom doba testování jedné transformace 0,48 sekundy. Celková doba registrace tedy 34560 sekund = 9,6 hodin. Je zřejmé, že obě dvě hodnoty jsou nepoužitelné.

38

4.1. REGISTRACE

POSUNŮ

39

Možností urychlení je několik. První je omezit prostor transformací. Omezení by však muselo být příliš drastické. Druhou možností je zvýšit hodnotu kroku. Tím ale dostáváme pouze hrubou registraci. Třetí možností je převzorkování obrazu na menší velikost. Ve své podstatě je to ale to stejné jako předchozí možnost. Závěr je zřejmý, úplné prohledávání nelze v tomto případě uspokojivě použít. 4.1.2

Rekurzivní prohledávání

Rekurzivní vyhledávání je naopak zcela jistě optimalizační metodou. Pro určité typy funkcí může být tato metoda, založená ve své podstatě na dělení intervalů, velice rychlá a přesná. Jelikož jak NCC, tak MI mají dle grafů z předchozí kapitoly v okolí extrémů poměrně velmi pěkný průběh (velká ZAF) bez lokálních extrémů, mohla by tato metoda být účinná. Co se týká nastavení, tak metoda má dva vstupní parametry a to prostor transformací a počet dělení rozsahu parametrů v každé dimenzi. Cím větší je počet dělení, tím snáze se odhalí globální extrém ukrytý mezi extrémy lokálními. Zkoušel jsem dvě nastavení: dělení na tři části ve všech dimenzích a na pět částí ve všech dimenzích. Při dělení na tři části byl průběh metody při použití NCC následující: Rekurzivní krok 1 (27 ohodnocení funkce) 2 (54 ohodnocení funkce) 3 (81 ohodnocení funkce) 4 (108 ohodnocení funkce) 5 (135 ohodnocení funkce)

Hodnota extrému -0.787417 -0.796471 -0.796471 -0.800011 -0.800439

Pozice globálního extrému (15, -15, 5) (7.5, -15, 2.5) (7.5, -15, 2.5) (9.375, -13.125, 3.125) (8.4375, -13.125, 2.8125)

Celkem tak bylo nutno provést 135 ohodnocení podobnostní funkce a celá registrace trvala 54 sekund. V případě dělení na pět částí byl extrém znovu úspěšně lokalizován. Tentokrát bylo potřeba provést tři rekurzivní kroky a 350 ohodnocení podobnostní míry, tedy zhruba dvakrát delší výpočet. Potvrdilo se tím, že funkce NCC má pro testované posuny dostatečně „pěkný" průběh a stačí tak použít menší počet dělení. V případě MI se výsledky opakovaly, pouze rozdíl mezi použitím tří a pěti dělících kroků se ještě více prohloubil. 4.1.3

G e n e t i c k ý algoritmus

Narozdíl od rekurzivního prohledávání má genetický algoritmus poměrně široké možnosti nastavení. Důležitá je velikost populace, počet generovaných generací, pravděpodobnost křížení a mutace a další. Pro testy jsem nastavil velikost populace na 30 jedinců, prav­ děpodobnost křížení (je implementováno dvoubodové křížení) na 95 % a pravděpodob­ nost mutace na 10%. Protože genetický algoritmus je algoritmem stochastickým, tedy při opakovaných spuštěních se stejným nastavením může dát jiná řešení, provedl jsem čtyři spuštění jak pro NCC, tak pro ML Maximální počet generací byl nastavený na šest (300 ohodnocení funkce). Pro NCC nedokázal genetický algoritmus nalézt extrém ani jednou ze čtyř pokusů a vždy skončil v jeho těsné blízkosti. Typické bylo rychlé nalezení oblasti extrému, ale poté

KAPITOLA

40

4. VÝBĚR OPTIMALIZAČNÍ

METODY

velmi pomalá konvergence k samotnému bodu extrému. Toto je bohužel obecně vlast­ nost genetický algoritmů. Pro MI byl algoritmus mírně úspěšnější a nalezl přesné řešení v jednom ze čtyř případů. V ostatních se opakovala situace jako v případě NCC. Tato relativní neúspěšnost však neznamená nevhodnost použití genetických algoritmů pro registraci. Naopak. Například pokud hrozí existence více lokálních minim, je genetický algoritmus vhodnější než rekurzivní prohledávání, protože je podstatně robustnější. Při použití pyramidové metody se například velmi hodí do vyšších pater pyramidy, kde není třeba lokalizovat extrém přesně, ale stačí určit přibližnou oblast kde se bude vyskytovat. 4.1.4

Simulované žíhání

Podobně jako genetický algoritmus je i simulované žíhání metodou stochastickou a velmi robustní, s širokou možností nastavení. Podstatná je především počáteční teplota systému a předpis žíhání. Protože testované funkce nevykazovaly nijak složitý průběh, bylo žíhání nastaveno tak, aby byly upřednostňovány spíše kroky znamenající snížení hodnoty funkce. Algoritmus byl spuštěn z bodu (0, 0, 0) a pro obě podobnostní míry byly provedeny čtyři testy. Metoda měla 100 % úspěšnost, tedy lokalizovala správně extrém ve všech osmi případech a potřebovala k tomu průměrně 200 ohodnocení podobnostní funkce. Ani simulované žíhání však není metodou samospasitelnou, jeho velká úspěšnost byla dána především hladkým průběhem funkce a malou vzdáleností extrému od bodu startu (což by však měl být případ většiny běžných registrací). Pro lepší hodnocení by bylo třeba provést větší množství důkladnějších testů, což však už překračuje rámec této práce. 4.1.5

A d a p t i v n í simulované žíhání

Adaptivní simulované žíhání (ASA) vzniklo jako vylepšení algoritmu simulovaného žíhání. Poněkud překvapivě má ve své základní verzi trochu méně nastavení, metoda si o způsobu dalšího postupu rozhoduje více sama. Nenastavuje se tak například počáteční teplota žíhání ani předpis žíhání. Znovu bylo provedeno osm testů, přičemž počet iterací algoritmu byl omezen na 300. Při použití NCC našlo ASA pozici extrému ve dvou případech a dvakrát skončilo po 300 iteracích v jeho těsné blízkosti. V případě MI bylo ASA stoprocentní a pro lokalizaci extrému potřebovalo průměrně 250 ohodnocení funkce. Úspěšnost ASA je tak velmi vysoká a je ověřeno, že při složitějších registracích (větší prostor transformací, rozdíly v citlivosti parametrů) je tato úspěšnost vyšší, než u standardního simulovaného žíhání. 4.1.6

Integrální transformace

Ve své podstatě jsou i integrální transformace optimalizační metodou. Budeme-li totiž uvažovat pouze posuny, umožňují nám jedním výpočtem získat míru podobnosti pro celý prostor transformací. Čas který takto ušetříme je velmi výrazný, jak je vidět z předchozích měření v kapitole o podobnostních mírách, tedy výhody jsou zřejmé. Pokusím se zde uvést však i jejich nevýhody, těch je hned několik: • Tímto postupem lze registrovat jen pomocí CC a NCC. Při výpočtu MI nám nepo­ může. • Nelze registrovat sub-pixelově.

4.2. REGISTRACE

ZMĚNY MĚŘÍTKA

41

• Bez dalších úprav není možné registrovat jiné transformace než posuny. • Paměťová náročnost. Potřebujeme-li počítat skutečné NCC (tj. bez periodicity ob­ razů) je integrální transformace opravdu extrémně paměťově náročná. Je totiž třeba držet v paměti několik obrazů o čtyřnásobném (respektive osminásobném v 3D pří­ padě) objemu obrazů původních. Proto už například zde testovaný obraz nebylo možné na testovacím stroji s 1 GB operační paměti registrovat.

4.2

Registrace změny měřítka

Pokud uvažujeme samotnou registraci změny měřítka se zafixovanými ostatními parametry (posun, rotace) nedostáváme prakticky žádný rozdíl od registrace posunů. Lze tedy použít všechny algoritmy testované v části o registraci posunů s podobnými výsledky. Jedinou výjimku tvoří integrální transformace, které v základní podobě registraci změny měřítka neumožňují. Jak však uvádí literatura [1], po úpravě výpočetních vztahů lze i tento postup použít. Pokud se registrují větší změny měřítka je dobré ještě jednou poznamenat, že pro lepší výsledky je dobré nastavit vyšší stupeň interpolace pro transformace obrazů.

4.3

Registrace otočení

Zde už dostáváme oproti registraci posunů několik rozdílů: • Tím hlavním je, že je třeba při registracích brát do úvahy vnitřní rozlišení obrazů (typicky v mikrometrech na pixel). • Při převzorkování obrazů (např. v pyramidové metodě) musí být dodrženy poměry stran, tedy jak šířka, tak výška a hloubka musí být změněny stejným poměrem. Při nedodržení budou výsledkem registrace nesmyslné hodnoty. • Je třeba dát velký pozor na pořadí otočení kolem jednotlivých os v implementaci, protože operace otočení v tomto smyslu není komutativní. Jinak platí zásady uvedené výše, stejně jako poznámka u registrace změny měřítka o použití integrálních transformací.

4.4

Registrace kombinací transformací

Registrace kombinací několika typů transformací je tím nejobtížnějším případem. Při sou­ časné registraci posunů s například změnou měřítka či rotací roste prostor transformací a složitost celého problému registrace zásadním způsobem. Typicky lze toto ukázat na příkladu rekurzivního vyhledávání, které se s rostoucím počtem dimenzí stává nepoužitel­ ným stejně jako úplné vyhledávání, protože složitost rekurzivního vyhledávání je na počtu dimenzí bohužel závislá exponenciálně1. Zde přichází vhod pyramidový algoritmus, který umožňuje další škálování celkové slo­ žitosti výpočtu. V horním patře (patrech) pyramidy najdeme na obraze převzorkovaném 1 Příklad: 6 parametrů, 3 dělení každého parametru = 3 6 testovaných transformací v každém rekurzivním kroku.

42

KAPITOLA

4. VÝBĚR OPTIMALIZAČNÍ

METODY

na menší velikost hrubou pozici extrému a tu poté v dalších patrech už pouze zpřesňujeme. Tohoto postupu jsem s úspěchem využil při pokusu registrovat oba obrazy s ohledem na posuny a změny měřítka. Tedy šestidimenzionální prostor transformací. Byly srovnány dvě varianty. Pro obě byla použita dvoupatrová pyramida, přičemž vrchní patro obsahovalo obraz převzorkovaný na poloviční velikost a spodní patro potom obraz v plném rozlišení. Rovněž omezení prostoru parametrů bylo pro obě varianty shodné. V horním patře se testovaly posuny v osách x a y mezi -30 a 30 voxely, v ose z mezi -10 a 10 voxely a změna měřítka mezi 0,9 a 1,2 pro všechny tři osy. V dolním patře už se pouze registrace upřesňovala v posunech ±4 voxely a změně měřítka ±0.05. Jako podobnostní míra byla vybrána MI. 1. varianta — V horním patře byl použit genetický algoritmus o velikosti populace 50 jedinců a bylo generováno 10 generací. V dolním patře byl potom použit algoritmus ASA limitovaný na 300 iterací. Po první fázi bylo nalezena hrubá pozice extrému pro posun (10, -14, 6) a změnu měřítka (1,04; 1,12; 0,95). Ta byla ve druhé fázi upřesněna metodou ASA na posun (11,2117; -12,3092; 2,60163) a změnu měřítka (1,06365; 1,08114; 0,90239). Výsledná hodnota MI byla -0.652004 (tedy značný posun). Během výpočtu byla podob­ nostní funkce ohodnocena celkem 800 x a registrace trvala 3 minuty. 2. v a r i a n t a — Zde byl použit algoritmus ASA jak pro horní, tak pro spodní patro. V horním patře byl limitován na 500 iterací, ve spodním potom opět na 300. Tento postup se ukázal být jako rychlejší, byl nalezen stejný extrém jako v předchozí variantě, byl však nalezen hned na začátku druhé fáze už po zhruba 60 iteracích. Ostatní parametry výpočtu se nezměnily. Na tomto místě lze samozřejmě navrhnout spoustu dalších postupů, které budou vhod­ nější pro ten který druh registrací. Nicméně po vizuální kontrole byla registrace obrazů velmi dobrá a i pro další obrazy dokázala tato metoda extrém úspěšně lokalizovat. Proto tento postup (zejména pak druhou variantu) doporučuji pro současnou registraci posunů a změny měřítka.

4.5

Další možnosti optimalizace

V rámci práce byla rovněž testována možnost registrace obrazů pouze na základě určitého objemového okénka. Tedy postup, kdy byla podobnostní funkce počítána jen v předem pevně stanoveném kvádru uvnitř obrazů. Cílem bylo další urychlení registrace. Toto zrych­ lení ale nebylo dostatečné a metoda přinesla výrazné zhoršení přesnosti registrací. Proto ji zde nebudu dále rozebírat.

Kapitola 5

Závěr 5.1

Výsledky

Práce se zabývala známými podobnostními mírami a jejich využitelností při registraci dodaných testovacích dat. Protože registrace objemově rozsáhlých obrazů se ukázala být časově velmi náročnou úlohou, byly rovněž vybrány některé známé optimalizační metody a byl zkoumán vliv jejich použití na rychlost a přesnost registrace. Jednotlivé podobnostní míry jsou představeny v teoretické části. V kapitole o výběru podobnostní míry je následně na příkladu uměle transformovaných dat ukázána oprávně­ nost použití těchto funkcí pro registraci dat. Na dalších, tentokrát již reálných datech, je podrobně zkoumána jejich rychlost, přesnost a citlivost na parametry. Nejlépe z tohoto srovnání vychází podobnostní míra MI (vzájemná informace), která jako jediná uspěla ve všech testech a vykazuje navíc i odolnost proti nelineární degradaci dat. Pro data nepostižená touto degradací (tedy většinu dodaných dat) je vhodná i podobnostní míra NCC (normalizovaný korelační koeficient). Tuje rovněž možné počítat rychlou Fourierovou transformací. S úspěchem jsou použitelné i oba Woodsovy algoritmy. Naopak mezi podob­ nostní míry, které ve srovnání neuspěly lze rozhodně zařadit SSC (stochastická změna znaménka), CC (korelační koeficient) a PCC (fázová korelace). U většiny podobnostních funkcí se prokázal hladký průběh v okolí globálního extrému, bez výskytů extrémů lokálních. Proto je v části o optimalizacích s úspěchem použito algo­ ritmu rekurzivního prohledávání s nízkým počtem dělení. I ostatní optimalizační metody však potvrdily oprávněnost svého zařazení do této práce a prokázaly výborné schopnosti při hledání globálního extrému. Pro registrace s širokým prostorem transformací se pak ukázala být nejvhodnější pyramidová metoda (zejména v kombinaci s algoritmem adaptiv­ ního simulovaného žíhání), kde je nejprve nalezena hrubá registrace na obraze převzorkovaném na menší velikost a následně je v jejím okolí tato registrace postupně zpřesňována. Diskutována byla také možnost registrace na základě objemových okének, která se ale ukázala být nevhodnou. Všechny postupy a výpočty prezentované v práci byly uskutečněny s pomocí registrační knihovny, která byla vyvíjena souběžně s prací. Jedná se o robustní knihovnu pro registraci šedotonních obrazů, pomocí podobnostních funkcí a optimalizačních metod uvedených v této práci. Při vývoji byl ovšem kladen důraz na snadnou rozšiřitelnost o nové funkce. Knihovna již nyní umožňuje například podrobné zaznamenávání průběhu registrace, či

43

KAPITOLA

44

5.

ZÁVĚR

vytváření grafických výstupů do programu Matlab. Podrobná dokumentace knihovny je uvedena v příloze práce.

5.2

Další rozšíření

Co se týká výhledů do budoucna, lze jistě počítat s dalším rozšiřováním registrační knihovny o nové funkce a zejména pak o nové optimalizační postupy umožňující dále zrychlit a zpřes­ nit registraci. Zajímavá by jistě byla implementace metod založených na spádu gradientu nebo úpravy pro možnost využití na víceprocesorových stanicích, či pro možnost běhu v paralelním distribuovaném prostředí. Další rozšíření by mohla směřovat do oblasti vícekanálových (barevných) obrazů.

Literatura [I] Čapek, M.: Registrace snímků z konfokálního mikroskopu [disertační práce]. ČVUT FEL, Praha, 1999 [2] Maintz, J. B. A, Viergever, M. A.: A Survey of Medical Image Registration. Medical Image Analysis, Vol. 2, No. 1, 1998 [3] Jedlička, M.: Využití genetických algoritmu pro optimalizaci výběru sanačních techno­ logií [diplomová práce]. FI MU, Brno, 2006 [4] Venot, A., Lebruchec, J., Golmard, J, Roucavrol, J.: An automated method for the normalization of scintigraphic images. Journal of Nuclear Medicine, 24, 1983 [5] Wagner, R., Galiana, H. L.: Evaluation of Three Template Matching Algorithms for Registering Images of the Eye. IEEE Trans. Biomed. Eng., Vol. 39, No. 12, December 1992 [6] Barnea, D. I., Silverman, H. F.: A Class of Algorithms for Fast Digital Image Regis­ tration. IEEE Trans. Comp., Vol. C-21, No. 2, p. 179-186, 1972 [7] Maes, F., Collignon, A., Vandermeulen, D., Marchal, G., Suetens, P.: Multimodality Image Registration by Maximization of Mutual Information. IEEE Trans. Med. Imag. Vol. 16, No. 1, 1996 [8] Pompe, B.: Measuring Statistical Dependencies in a Time Series. Journal of Statistical Physics, Vol. 73, p. 195-206, 1993 [9] Woods, R. P., Cherry, S. R., Mazziotta, J. C : Rapid Automated Algorithm for Aligning and Reslicing PET Images. Journal of Computer Assisted Tomography, Vol. 16, No. 4, 1992 [10] Woods, R. P., Mazziotta, J. C , Cherry, S. R.: MRI-PET Registration with Automated Algorithm. Journal of Computer Assisted Tomography, Vol. 17, No. 4, 1993 [II] Svobodvá Vařeková, R.: Optimalizace Dokument dostupný na URL http://ncbr.chemi.muni.cz/ nl9n/vyuka/optimalizace/ (prosinec 2006) [12] Ingber, L.: Papers Using ASA. Dokument dostupný na URL http://www.ingber.com/asa_papers.html (prosinec 2006)

45

46

[13] Eckel, B.: Myslíme v jazyku C++. Grada, 2000. [14] Frigo, Matteo and Johnson, Steven G.: FFTW User Manual. Dokument dostupný na URL http://www.rTtw.org/rTtw3_doc/ (prosinec 2006) [15] Kitware, Inc.: CMake documentation. Dokument dostupný na URL http://www.cmake.org/HTML/Documentation.html (prosinec 2006) [16] Wikipedia article: Cross-correlation. Dokument dostupný na URL http://en.wikipedia.org/wiki/Cross-correlation (prosinec 2006) [17] Wikipedia article: Mutual information. Dokument dostupný na URL http://en.wikipedia.org/wiki/Mutual_information (prosinec 2006) [18] Wikipedia article: Genetic algorithm. Dokument dostupný na URL http://en.wikipedia.org/wiki/Genetic_algorithm (prosinec 2006) [19] Wikipedia article: Simulated annealing. Dokument dostupný na URL http://en.wikipedia.org/wiki/Simulated_annealing (prosinec 2006) [20] Článek ze serveru Wikipedia: Chromozóm. Dokument dostupný na URL http://cs.wikipedia.org/wiki/Chromozom (prosinec 2006) [21] Článek ze serveru Wikipedia: Buňka. Dokument dostupný na URL http://cs.wikipedia.org/wiki/Bunka (prosinec 2006) [22] Článek ze serveru Wikipedia: Chromatická aberace. Dokument dostupný na URL http://cs.wikipedia.org/wiki/Chromaticka_aberace (prosinec 2006)

LITERATURA

Příloha A

Obsah CD K této práci je přiložen kompaktní disk obsahující: Adresář

Obsah

/data/

Některá data uváděná v této práci, výstupy registrací, výstupy do Matlabu.

/examples/

Příklady použití registrační knihovny. Zdrojové kódy i přeložené spusti­ telné soubory.

/obrazy/

Testované obrazy.

/viewer/

Prohlížeč 3D obrazů

/text/

Zdrojové texty této diplomové práce pro systém I^TßX

/voxreg/

Registrační knihovna

47

Příloha B

Registrační knihovna Zde se nachází podrobný popis registrační knihovny, jejích tříd a rozhraní. Podstatou knihovny je třída RegSimFunc, ze které jsou odvozeny všechny podobnostní míry a třída RegAlgo, která je předkem všech optimalizačních metod. Celá knihovna je postavena na knihovně ißd 1 a knihovně fftw 2 .

B.l

Podobnostní funkce Diagram tříd

RegSimFunc Abstraktní třída pro podobnostní míru. Všechny další podobnostní míry musí být jejími potomky. Mimo jiné si třída pamatuje dosud nejlepší nalezené řešení, jeho hodnotu a počet kolikrát byla volána. Tyto hodnoty tedy není třeba uchovávat jinde v optimalizačních metodách, či hlavním programu. Metody této třídy jsou následující: RegSimFunc - konstruktor. V odvozených třídách by měla být v konstruktoru nastavena proměnná m-name, udávající jméno podobnostní míry. 1

Dostupná v laboratoři http://www.fl.muni.cz/lom/ 2 http://www.fftw.org/

optické mikroskopie na Fakultě informatiky

48

Masarykovy

univerzity,

B.l.

PODOBNOSTNÍ

FUNKCE

49

ComputeFunc - privátní virtuální metoda, která musí být implementována v potomkovi. Této metodě je předán obraz A a transformovaný obraz B. Metoda počítá na těchto obrazech vlastní podobnostní míru a vrací její hodnotu. ComputeFuncVoi - tato privátní virtuální metoda slouží ke stejnému účelu jako předchozí metoda, počítá ale podobnostní míru pouze nad zvoleným objemovým okénkem. Nemusí být implementována (nyní je implementována pouze pro míry SSC, SAVD, SSD, CC a NCC). SetupFunc - tato privátní virtuální metoda je volána metodou Start před samotným zahá­ jením registrace. Zde je možné inicializovat proměnné, případně předupravit obrazy před registrací (v MI například použito pro alokování paměti pro histogramy a snížení počtu kvantizačních úrovní). Start - metoda volaná z metody Init třídy RegAlgo před započetím registrace. Je zde vynulován počet volání funkce, upraveno objemové okénko podle koeficientů převzorkování a volána metoda SetupFunc. Tato metoda by neměla být v dceřinných třídách znovu implementována. End - metoda volaná z optimalizačních algoritmů po skončení registrace. Vrací nejmenší nalezenou hodnotu podobnostní funkce a její pozici v prostoru transformací. GetName - vrací jméno podobnostní míry. LogParameters - metoda sloužící k zaznamenávání průběhů registrace. Měla by vypisovat hodnotu všech podstatných parametrů. Compute - metoda volaná z optimalizačních algoritmů pro výpočet jedné konkrétní hod­ noty podobnostní funkce. Metoda má tři parametry a to obraz A a transformovaný obraz B a danou transformaci. Tato metoda volá privátní metodu ComputeFunc (případně Com­ puteFuncVoi), zajišťuje uchování nejlepší dosažené hodnoty a také volá zpětně (callback) uživatelskou funkci, pokud je nějaká nastavená. SetCallback - nastavuje uživatelskou funkci pro zpětné volání z metody Compute. Toto využívá například třída RegGraphExp_c, která zajišťuje generování grafů. SetVOI - nastavuje objemové okénko, nad kterým bude podobnostní funkce počítána. Standardně je počítána nad celým obrazem. GetNumberOfCalls - vrací počet volání podobnostní funkce. RegFuncSSC Implementace stochastické změny znaménka. RegFuncSAVD Implementace součtu absolutních hodnot rozdílů. RegFuncSSD Implementace součtu čtverců rozdílů. RegFuncCC Implementace výpočtu korelačního koeficientu. RegFuncNCC Implementace výpočtu normalizovaného korelačního koeficientu.

PRÍLOHA B. REGISTRAČNÍ

50

KNIHOVNA

RegFuncMI Implementace výpočtu vzájemné informace. SetMaxLevel - metoda pro nastavení počtu kvantizačních úrovní. Důležité: nastavuje se počet bitů potřebných k uložení hodnoty intenzity, tedy 8 = 256 k.ú., 7 = 128 k.ú., ... RegFuncGMI Implementace výpočtu zobecněné vzájemné informace. SetMaxLevel - metoda pro nastavení počtu kvantizačních úrovní. Důležité: nastavuje se počet bitů potřebných k uložení hodnoty intenzity, tedy 8 = 256 k.ú., 7 = 128 k.ú., ... RegFuncWRIU Implementace výpočtu Woodsova algoritmu Ratio-Image Uniformity. RegFuncWPIU Implementace výpočtu Woodsova algoritmu Partitioned Intensity Uniformity. SetMaxLevel - metoda pro nastavení počtu kvantizačních úrovní. Důležité: nastavuje se počet bitů potřebných k uložení hodnoty intenzity, tedy 8 = 256 k.ú., 7 = 128 k.ú., ...

B.2

Optimalizační metody D iagram tříd

RegAlgo Abstraktní třída pro optimalizační metody. Všechny další optimalizační metody musí být jejími potomky. Mimo jiné třída zajišťuje převzorkování obrazů, výpočet doby registrace a kontrolu hodnot parametrů. Třída je šablonou, která je závislá na typu vstupních obrazů. V současné době jsou v knihovně explicitní instance pro obrazy s šedotonními osmibito­ vými intenzitami. Metody této třídy jsou následující: Setlmages - této metodě jsou předány vstupní obrazy, které se mají registrovat. SetSimilarityFunction - této metodě je předána podobnostní míra, kterou se bude regis­ trovat.

B.2. OPTIMALIZAČNÍ

METODY

51

GetSimilarityFunction - vrací odkaz na momentálně nastavenou podobnostní míru. SetlnterpolationDegree - nastavení stupně interpolace použitého při transformacích. Resample - nastavení poměru převzorkování a stupně interpolace použitého při převzorkování. Poměr převzorkování je vektor 3 hodnot mezi 0 a 1. Těmito hodnotami je vynásobena šířka, výška a hloubka obrazu, čímž jsou získány rozměry na jaké se obraz převzorkuje. Důležité je, že například při zmenšení obrazu na polovinu se podobně změní i objemové okénko v podobnostní funkci, dále krok po kterém postupuje algoritmus úplného vyhledá­ vání a další parametry! Register - touto metodou se spouští vlastní registrace. Metodě jsou předány 2 parametry: min a max. Tyto parametry mohou mít 3, 6 nebo 9 prvků a určují prostor transformací. První 3 hodnoty omezují posuny, další 3 změnu měřítka a poslední 3 rotace, všechny v pořadí podle os x, y, z. Je dobré nezapomenout, že pokud chceme registrovat například pouze rotace, musí být u 4 až 6 hodnoty (změny měřítka) min i max nastaveno na 1, nikoliv 0! GetName - vrací jméno optimalizační metody. LogParameters - metoda se stará o logování parametrů optimalizační metody. RegAlgoFS Implementace úplného vyhledávání. SetStep - nastavení kroku po kterém bude postupovat metoda úplného vyhledávání. Délka tohoto vektoru by měla být větší nebo rovna délce vektorů min a max, předaných metodě Register. ExportGraph - tato metoda se stará o export grafu do souboru čitelného programem Matlab. Jako parametr je třeba předat proměnnou typu RegGraphExp_c RegAlgoGrid Implementace rekurzivního vyhledávání. SetNumOfDivs - metodě se předává vektor, jehož prvky určují na kolik částí se budou v každém kroku dělit jednotlivé dimenze prostoru transformací. Délka tohoto vektoru by opět neměla být menší než je délka vektorů min a max. SetStop - při postupu rekurzivní metody dochází k postupnému zmenšování prostoru transformací. Vektor předaný této metodě určuje kdy má algoritmus zastavit. Algoritmus končí v okamžiku kdy je velikost prohledávaného intervalu ve všech dimenzích (max - min) menší než odpovídající hodnota ve stop vektoru. UseFast - metoda umožňuje zapnutí rychlejší varianty algoritmu. Jakmile se prohledá­ vaný interval pro některý z parametrů zmenší pod hodnotu stop vektoru, je určena finální hodnota tohoto parametru jako střed posledního intervalu. Tedy tento interval není již dále dělen a to ani pokud intervaly ostatních parametrů nedosáhly svých stop hodnot. RegAlgoGenetic Implementace genetického algoritmu. SetParams - nastavení genetického algoritmu. Metodě se předává požadovaná velikost po­ pulace, počet elitních jedinců , kteří automaticky přejdou do další generace, počet nových náhodných jedinců generovaných do 2 a další generace, pravděpodobnosti křížení a mutace a celkový počet generací.

PRÍLOHA B. REGISTRAČNÍ

52

KNIHOVNA

RegAlgo Annealing Implementace simulovaného žíhání. SetParams - nastavení parametrů simulovaného žíhání. Počet iterací, počáteční teplota, po kolika krocích dojde k snížení teploty a žíhací předpis (tedy číslo mezi 0 a 1, kterým se násobí teplota). RegAlgo ASA Implementace adaptivního simulovaného žíhání. SetParams - nastavení parametrů algoritmu ASA. Počet iterací, posun podobnostní míry (přičte se ke všem ohodnocením, ASA nefunguje pro záporné hodnoty) a další 2 parametry, které doporučuji ponechat na původních hodnotách (slouží k jemnému doladění metody ASA). RegAlgoPyramid Implementace pyramidového algoritmu. AddPhase - přidání dalšího stupně do pyramidy. Tento stupeň je tvořen optimalizační metodou a parametry min a max. Pro první stupeň (vrchol pyramidy) jsou parametry min a max absolutní, pro další stupně pak relativní vzhledem k extrému nalezenému na předchozím stupni. Flush - metoda slouží k vyprázdnění pyramidy. RegAlgo FourierCC Implementace výpočtu korelačního koeficientu Fourierovou transformací a nalezení nej­ větší hodnoty v oblasti zadané parametry min a max. Je uvažována periodicita obrázku. RegAlgo FourierNCC Implementace výpočtu normalizovaného korelačního koeficientu Fourierovou transformací a nalezení největší hodnoty v oblasti zadané parametry min a max. RegAlgoFouňerNCC - konstruktor má jeden parametr a tím je zapnutí/vypnutí paddingu. Defaultně je vy­ pnutý. Padding - metoda pro zapnutí/vypnutí paddingu (doplnění obrázků nulami do dvojnásobku velikosti). Vypnutý padding znamená, že výpočet bude probíhat s uvažova­ nou periodicitou obrázku. RegAlgo FourierPCC Implementace výpočtu fázové korelace a nalezení největší hodnoty v oblasti zadané para­ metry min a max. Je uvažována periodicita obrázku.

B.3

Logovací třída

Knihovna obsahuje globální proměnou RegLog, která slouží k podrobnému zaznamenávání průběhu registrací. Lze použít 3 důležité metody: LogToScreen - zapne výpis zaznamenávaných událostí na obrazovku. LogToFile - zapne ukládání výpisu do zvoleného souboru. LogString - metoda používaná k zápisu do registru událostí.

BA.

EXPORT

B.4

GRAFŮ

53

Export grafů

Export grafů umožňuje optimalizační metoda úplné vyhledávání a pak také všechny me­ tody založené na Fourierově transformaci. Export se provádí voláním metody ExportGraph, které je předán odkaz na proměnou typu RegGraphExp.c. Po skončení registrace je ještě třeba zavolat metodu SaveGraph. Pro správnou funkčnost je nutné, aby prostor transfor­ mací prohledávaných při registraci byl dvoudimenzionální.

B.5

Hlášení chyb

Report chyb provádí knihovna přes globální proměnou RegErr. Zda došlo při běhu k nějaké chybě zjistíme její metodou Ok. Výpis potom metodou PopError.

B.6

Statické funkce

Nejdůležitější statickou funkcí v knihovně je procedura RegTransform, která zajišťuje jednotlivé transformace obrazů (posuny, změny měřítka, rotace) a převzorkování. Tato procedura je postavena na kódu Philippe Thevenaze (kontakt je uvedený ve zdrojovém souboru), speciálně upravené pro zvýšení rychlosti registrací. Důležité je, že v současné implementaci procedura RegTransform nebere ohled na vnitřní rozlišení obrazů. Překlad knihovny je možné provést nejjednodušeji pomocí utility CMake3. Podrobně dokumentované příklady použití jsou na doprovodném CD.

3

http://www.cmake.org

Příloha C

Podklady pro subjektivní hodnocení Protože v práci byla několikrát vybrána nejvhodnější registrace na základě subjektivního hodnocení, uvádím zde příklady jak může vypadat porovnání výsledku jednotlivých po­ dobnostních funkcí. Následující dva obrázky sestávají vždy z několika obrazů, které jsou vytvořeny současným zobrazením obrazu A (červená složka) a obrazu B (zelená složka) registrovaného pomocí podobnostních mír uvedených v popisku (od leva):

Obrázek C.l: Srovnání registrace pomocí NCC, MI a PIU u obrazů z obr. 3.4

Obrázek C.2: Srovnání registrace pomocí NCC a MI u obrazů z obr. 3.10

54

Příloha D

Ukázky registrace dalších dat D.l

První pár

Sum of squared diffei

5916-HDD7 pra TA [pixel] = -6, TY [pixel] = 5

Mutual information; eitrem =-1.0561 pro TX [pixel] = -6, TY [pixel] = 4

Normalized cross correlation; e Jitrem = -0.95427 pro TX [pixel] = -6, TY [pixel] = 4

Woods (Partitioned Intensity Uniformity); extrem = 0.254562 pro TV, [pixel] = -6, TY [pix.

55

PRÍLOHA D. UKÁZKY REGISTRACE Podobnostní míra SSD NCC MI PIU

2

D o b a registrace 10 s 13 s 20 s 17 s

DALŠÍCH DAT Pozice e x t r é m u (-6, 5) (-6, 4) (-6, 4) (-6, 4)

Druhý pár

Sum oí scared differences; extrem = 2.03537e-HJ07 pra TX [pixel] = 4, TY [pixel] = -1

TX [pixel]

Mutual information; extrem = -0.715235 pro TX [pixel] = 3, TY [pixel] = -1

TY [pixel]

Podobnostní míra SSD NCC MI PIU

TX [pixel]

Normalized cross correlation; extrem - -0.93198 pro TX [pixel] = 4, TY [pixel] - -1

TY [pixel]

TX [pixel]

Woods (Partitioned Intensity Uniformity); extrem = 0.356054 pro TX [pixel] = 3, TY [pixel] = -1

TY [pixel]

TX [pixel]

D o b a registrace 10 s 14 s 21 s 19 s

Pozice e x t r é m u (4, -1) (4, -1) (3, -1) (3, -1)

D.3.

TRETIPAR

D.3

Třetí pár

equsrsd diforoncse, smrsm = 7.S32SS»«07 pra TX [pi.sl] = -40, T

57

l24pro-rX[pi«el] = 7.TY [pi»el] = -4

i^l

Mutual information; extrem = -0.773892 pro TX [pixel] = 5, TY [pixel] = -1

Podobnostní míra SSD NCC MI PIU

Woods (Partitioned Intensity Uniformity); extrem = 0.44597 pro TX [pixel] = 4. TY [pixel] = -2

D o b a registrace 8s 11 s 19 s 16 s

Pozice e x t r é m u (-40, -6) (7, -4) (5, -1) (4, -2)