POČÍTAČOVÁ FYZIKA II
Rudolf Hrach
PF UJEP Ústí nad Labem 2003
2
Obsah 1 Úvod
3
2 Počítačová grafika a vizualizace 2.1 Grafický hardware . . . . . . . . . . . . . 2.2 Grafický software . . . . . . . . . . . . . . 2.2.1 Typy grafického software . . . . . . 2.2.2 Souřadné systémy . . . . . . . . . . 2.3 Vybrané algoritmy počítačové grafiky . . . 2.3.1 Společné problémy rastrové grafiky 2.3.2 Kreslení úsečky . . . . . . . . . . . 2.3.3 Kreslení elipsy . . . . . . . . . . . . 2.3.4 Vyplňování oblasti . . . . . . . . . 2.3.5 Kreslení prostorových objektů . . . 2.4 Vizualizace . . . . . . . . . . . . . . . . . 2.5 Shrnutí . . . . . . . . . . . . . . . . . . . . 2.6 Problémy k řešení . . . . . . . . . . . . . . 3 Zpracování obrazu 3.1 Formulace problému . . . . . . . . . . . . 3.2 Metody zpracování obrazu nízké úrovně . . 3.2.1 Digitalizace . . . . . . . . . . . . . 3.2.2 Geometrické transformace . . . . . 3.2.3 Filtrace . . . . . . . . . . . . . . . 3.2.4 Binarizace . . . . . . . . . . . . . . 3.2.5 Rozpoznávání objektů . . . . . . . 3.3 Metody zpracování obrazu vysoké úrovně . 3.3.1 Integrální charakteristiky . . . . . . 3.3.2 Informace o jednotlivých objektech 3.3.3 Informace o rozložení objektů . . . 3.4 Další problémy z oblasti zpracování obrazu 3.5 Shrnutí . . . . . . . . . . . . . . . . . . . . 3.6 Problémy k řešení . . . . . . . . . . . . . .
1
. . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . .
6 6 8 8 10 11 12 15 17 19 23 25 29 30
. . . . . . . . . . . . . .
31 31 33 33 34 34 43 45 47 48 49 52 57 58 58
OBSAH 4 Integrální transformace 4.1 Fourierovy řady . . . . . . . . . . 4.2 Fourierova transformace . . . . . 4.3 Diskrétní Fourierova transformace 4.4 Rychlá Fourierova transformace . 4.5 Další integrální transformace . . . 4.6 Shrnutí . . . . . . . . . . . . . . . 4.7 Problémy k řešení . . . . . . . . .
2
. . . . . . .
60 61 62 63 63 64 65 65
5 Další směry počítačové fyziky 5.1 Řízení experimentu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Symbolické manipulace . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
66 66 68
6 Závěr
69
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
Kapitola 1 ÚVOD Úvodem druhé části studijního textu „Počítačová fyzika IIÿ si zopakujeme to, co jsme uvedli v jeho první části. Počítačová fyzika patří mezi nové a prudce se rozvíjející směry vědeckého bádání. Název, jenž vznikl spojením dvou pojmů – fyzika a počítač – je výstižným popisem jejího zaměření. Na rozdíl od jiných směrů vědecké práce, kde charakterizujícím faktorem je objekt studia, v případě počítačové fyziky je rozhodující metodika řešení fyzikálního problému. Tím dochází k omezování již dlouho probíhajícího a ne vždy vítaného trendu stále hlubší specializace ve vědecké práci. Z hlediska metodického se vědecká práce ve fyzice již dlouho dělí na dva základní směry – na experimentální a teoretický. Experimentální fyzika produkuje data, dává podněty pro vytváření teorií a umožňuje tyto teorie testovat. Teoretická fyzika interpretuje a zobecňuje experimentální poznatky a navrhuje další experimenty sloužící jednak k prohloubení znalostí ve studovaných oblastech a jednak dávající podněty pro studium oblastí nových. Tato symbióza obou přístupů velmi úspěšně probíhá po mnoho staletí, resp. tisíciletí. V posledních několika desetiletích však začíná být doplňována právě počítačovou fyzikou. Ta se snaží popsanou dvoučlennou spolupráci experimentální a teoretické fyziky změnit na rovnocennou spolupráci tří partnerů: experimentální, teoretické a počítačové fyziky. Vedle již popsaných vazeb historických směrů tak přibývají další vazby. Experimentální fyzika generuje data využívaná jak teoretickou tak i počítačovou fyzikou, počítačová fyzika naopak pro experimentální fyziku provádí analýzu těchto dat, zabezpečuje řízení experimentálních zařízení, vytváří modely reálných procesů a na jejich základě analyzuje experimentální data a navrhuje nové experimenty. Podobná vazba se vytváří i mezi počítačovou a teoretickou fyzikou: teoretická fyzika poskytuje rovnice pro řešení metodami počítačové fyziky a naopak interpretuje výsledky počítačovými metodami získané, počítačová fyzika pak poskytuje teoretické fyzice pomoc při provádění rozsáhlých výpočtů a podobně jako experimentální fyzika dává prostřednictvím tzv. počítačových experimentů podněty pro vytváření nových teorií a pomáhá při jejich testování. Takto nastíněná spolupráce tří rovnocenných metodik řešení fyzikálních problémů je v současné době pro počítačovou fyziku zatím ještě poněkud ambiciozním plánem, s rozvojem hardwarových a softwarových prostředků se však stává stále reálnější. 3
KAPITOLA 1. ÚVOD
4
Co to vlastně je počítačová fyzika? Přesná definice tohoto vědního oboru je velmi obtížná, ne-li nemožná. Stručně řečeno, jsou to takové postupy při řešení fyzikálních problémů, při kterých počítač hraje podstatnou roli. Důležité je slovo „podstatnouÿ, protože v současné době najdeme počítač v každé laboratoři, kde se užívá při drobnějších výpočtech, psaní textů (publikací, dizertačních prací), sběru dat, atd. Do počítačové fyziky určitě patří např. počítačové simulace nebo symbolické manipulace, částečně i automatizace experimentu, apod. S rozvojem hardware se vytvářejí i nové oblasti počítačové fyziky jako např. práce v počítačových sítích nebo studium paralelních procesů. Počítačová fyzika jednak poskytuje výstup výsledků do jednotlivých fyzikálních disciplin a jednak vytváří nové metody práce, čímž dále počítačovou fyziku jako disciplínu rozvíjí. Hranice oboru jsou z tohoto důvodu mlhavé a neustále se měnící. Pro potřeby tohoto studijního textu byla počítačová fyzika rozdělena do následujících základních směrů, kterým se budeme postupně věnovat: • Počítačové modelování • Počítačová grafika a vizualizace • Zpracování obrazu • Integrální transformace Z těchto oblastí má v počítačové fyzice výsadní postavení počítačové modelování, což je technika používaná fyziky různých specializací nejčastěji, neboť jim poskytuje často neocenitelnou pomoc při experimentálním i teoretickém studiu nejrůznějších fyzikálních jevů a procesů. Je také nejvíce rozpracována a dělí se na řadu dílčích metodik, které jsou často známé i mimo rámec počítačové fyziky - např. metoda Monte Carlo. Proto se i my budeme této oblasti počítačové fyziky věnovat podrobněji a bude tvořit těžiště prvního dílu tohoto studijního textu. Ostatní tři vyjmenované směry počítačové fyziky budou tvořit náplň druhého dílu studijního textu. Jelikož vyjmenovanými základními směry počítačová fyzika zdaleka nekončí, v závěru druhého dílu se stručně zmíníme i o dalších metodikách počítačové fyziky, které buď patří do počítačové fyziky jen okrajově (např. Řízení experimentu) a nebo jsou tak rozsáhlé a obtížné, že překračují rámec tohoto studijního textu (Symbolické manipulace, Teorie waveletových funkcí, atd.). I když počítačová fyzika patří mezi nejmladší vědecké disciplíny, již se rozvinula do takové šíře a dosáhla tolika výsledků, že se začíná vnitřně dělit. Stále častěji se hovoří o • klasické počítačové fyzice • moderní počítačové fyzice. Pod pojem moderní počítačová fyzika se zahrnují nové progresívní směry, které se objevily v posledních letech buď přímo v rámci počítačové fyziky a nebo byly jako aplikace převzaty z jiných vědeckých oblastí (nejčastěji z matematiky a informatiky). Do moderní počítačové fyziky v současné době patří neuronové sítě (resp. jejich aplikace ve fyzice), fuzzy logika a
KAPITOLA 1. ÚVOD
5
evoluční programování (genetické algoritmy), přičemž tento seznam je nadále otevřený. Z tohoto hlediska předložený studijní text se věnuje klasické počítačové fyzice, pouze v závěru druhého dílu budou velmi stručně nastíněny problémy řešené postupy moderní počítačové fyziky. Závěrem ještě jednu poznámku ke členění textu. Základní informace v textu obsažené jsou určené pro bakalářskou úroveň studia fyziky. Pokud to ale bude vhodné, bude výklad doplněn tak, aby se zájemce (nebo případně posluchač magisterského studia) dozvěděl o příslušné problematice poněkud více. Rozšiřující informace v tomto studijním textu budou od základního výkladu odděleny graficky – budou psány drobnějším písmem a odsazeny od okrajů stránky.
K účelu doplnění a rozšíření znalostí získaných při přednáškách slouží též seznam literatury uvedený na konci obou dílů studijního textu. Naneštěstí neexistuje dostatečné množství literatury k jednotlivým partiím počítačové fyziky v češtině, proto je většina uváděných knih v angličtině a v knihovnách jsou jen v omezeném počtu. Pro bakalářské studium však čtení rozšiřující literatury není nezbytné (i když samozřejmě je užitečné). Každá delší kapitola obou dílů studijního textu bude zakončena sekcí Shrnutí, v níž budou uvedeny nejdůležitější poznatky v příslušné kapitole popsané. Tyto poznatky budou též vyžadovány u zkoušek – dílčích i státních. Dále budou v některých kapitolách uvedeny příklady určené k procvičování látky v kapitole popsané. Některé z těchto příkladů budou již vyřešeny, většina však bude určena k samostudiu, případně k procvičení v rámci organizované výuky.
Kapitola 2 POČÍTAČOVÁ GRAFIKA A VIZUALIZACE Mezi základními směry práce v počítačové fyzice existují dva, které místo s numerickými daty pracují s obrazovou informací. Jedním z těchto směrů je počítačová grafika spolu s vizualizací velkých souborů dat a druhým směrem zpracování obrazů. Při počítačové grafice vycházíme z numerických údajů, které jsou uloženy v počítači jako výsledek předchozích operací, např. počítačového modelování nebo zpracování experimentálních dat, a my tuto informaci chceme získat jako výstup v grafické podobě. Důvodem preferování grafické formy prezentování vědeckých výsledků je, že lidské smysly jsou přednostně orientované na vnímání vnějšího světa pomocí obrazů (člověk získává převážnou většinu podnětů z okolí pomocí zraku). Proto též jeden dobře navržený graf nebo jiný typ grafického výstupu výsledků nám dá v přehledné podobě mnohem více informací o chování studovaných systémů než podrobné tabulky s výsledky měření nebo modelování. Naproti tomu při zpracování obrazů jsou obrazová data vstupní informací pro další výpočty. Tato situace je běžná např. v astronomii, kde výsledkem pozorování jsou fotografie astronomických objektů, ale i v jiných oblastech fyziky, kde se soustavně pracuje s obrazy z optických nebo elektronových mikroskopů. Tyto obrazy je třeba zdigitalizovat, přenést do počítače a v něm analyzovat a to je právě úkol směru počítačové fyziky zvaného zpracování obrazů. Počítačové grafice a vizualizaci se budeme věnovat v této kapitole studijního textu, zatímco podrobné informace o metodách zpracování obrazu budou tvořit obsah kapitoly příští.
2.1
Grafický hardware
K tomu, abychom mohli získat grafický výstup z našich výpočtů, simulací nebo měření, musíme mít počítač vybavený grafickým výstupním zařízením. Takových zařízení existuje celá řada a my je budeme z hlediska formy získávaných výstupů dělit na dvě skupiny, na
6
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
7
– zařízení poskytující výstupní informaci v trvalé formě, již lze následně využít v dizertačních pracích, vědeckých publikacích, apod., a – zařízení poskytující výstup na obrazovce (existují ovšem postupy, jak i tyto dočasné výstupy učinit trvalými). Mezi zařízení první skupiny patří různé typy zapisovačů (digitální souřadnicové zapisovače, analogové zapisovače a fotozapisovače) a v poslední době především tiskárny (mozaikové, inkoustové a laserové). O kvalitě grafického výstupu a oblastech použití rozhodují parametry těchto zařízení, především jejich rozlišení a možnost barevného tisku. V případě zapisovačů s krokovými motory je typické prostorové rozlišení velikosti kolem 0, 1 mm a rozměr výsledného obrazu se může pohybovat od několika centimetrů (u stolních zařízení určených pro vytváření obrázků pro vědecké publikace) až po několik metrů (u specializovaných zapisovačů určených pro kreslení kartografických a meteorologických map, apod.). U tiskáren charakterizujeme kvalitu jejich výstupu pomocí parametru dpi (dot per inch), udávajícího kolik bodů je tiskárna schopna vytisknout na úsečku délky 25, 2 mm. V současné době bývá minimálního hodnota tohoto parametru 300 dpi nebo 600 dpi a kvalitní tiskárny umožňují tisk i s rozlišením 1200 dpi a větším. Přepočítáme-li toto rozlišení např. na stránku A4, uvidíme, že jsme schopni vytisknout obrázek obsahující řádově 104 bodů v každém směru, což výrazně překračuje možnosti běžných zapisovačů a proto se tiskárna stává základním výstupním zařízením. Navíc většina tiskáren umí různými mechanismy měnit velikost tištěného bodu, což dále jejich možnosti rozšiřuje. Druhá skupina výstupních zařízení je tvořena grafickými displeji, poskytujícími výstup na obrazovce. Konstrukčně existují tři typy displejů – alfanumerické, umožňující výstup pouze textové informace a proto pro počítačovou grafiku nevhodné, rastrové a vektorové. Nejběžnější jsou displeje rastrové, s kterými se setkáváme u mikropočítačů (a které již dávno známe u televizorů). U těchto displejů jsou obrazové body, tzv. pixely, rozloženy v pravidelné síti (obvykle čtvercové) a o výsledném grafickém rozlišení rozhoduje počet pixelů ve směru vodorovném a svislém. V oblasti PC mikropočítačů je obvyklá stupnice možných grafických rozlišení rastrových monitorů začínající v současné době u rozlišení 640 × 480 (tzv. rozlišení VGA) a pokračující až po 1600 × 1200, specializované rastrové displeje umožňují ovšem pracovat i s rozlišením podstatně vyšším. Pokud se podíváme podrobněji na řadu možných grafických rozlišení u grafických karet a monitorů u PC mikropočítačů, tj. 640×480, 800×600, 1024×768, 1280×1024 a 1600 × 1200, vidíme, že s výjimkou rozlišení 1280 × 1024 jsou vodorovné a svislé rozměry vzniklých obrazů v poměru 4 : 3. Při použití rozlišení 1280 × 1024 mohou proto nastat problémy s přenositelností grafických programů.
Vedle rastrových displejů existují i tzv. displeje vektorové. Při práci s nimi zadáváme počáteční a koncový bod úsečky, kterou je třeba vykreslit, a elektronový paprsek na monitoru tyto body propojí spojitou čarou (na rozdíl od obecně schodovité čáry u displejů
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
8
Obrázek 2.1: Vektorová (vlevo) a rastrová (vpravo) grafika. rastrových) – viz obr. 2.1. Vektorové displeje patří spíše do profesionální oblasti a proto se s nimi na školách běžně nesetkáme. V terminologii rastrových a vektorových zařízení můžeme popsat i grafická výstupní zařízení první skupiny. Podle tohoto dělení tiskárny a digitální souřadnicové zapisovače jsou zařízeními rastrovými, zatímco starší analogové zapisovače a fotozapisovače patří mezi zařízení vektorová. Při programování grafických výstupů je třeba pečlivě rozlišovat, do které skupiny dané výstupní zařízení patří, protože použité algoritmy se též budou dělit na rastrové a vektorové. Při běžné práci však budeme používat rastrová zařízení, ať již displeje nebo tiskárny, a proto i popisovaný software bude odpovídat tomuto typu zařízení. Proto též nebudeme uvádět, že se jedná o algoritmy pro rastrová zařízení a pouze v případě odchylek pro zařízení vektorová na to zvláště upozorníme. Vedle obvyklých displejů a tiskáren se můžeme setkat i s dalšími grafickými výstupními zařízeními. Z nichž největší perspektivu zřejmě mají zařízení, která budou umožňovat výstup trojrozměrné obrazové informace. Rozpracováno je několik principů – 3D monitory, stereoskopické brýle nebo zařízení pracující na principu holografie. Tato zařízení jsou velmi vhodná především pro vizualizaci velkých souborů dat, statických i dynamických.
2.2 2.2.1
Grafický software Typy grafického software
V krátké historii vývoje software pro práci s grafikou bylo realizováno několik různých možností (o konkrétních algoritmech pro grafický výstup budeme mluvit v následující sekci). • Klasické programovací jazyky (FORTRAN, Pascal, atd.) neobsahovaly žádné prostředky pro grafickou práci. Jedinou možností proto byla tzv. pseudografika, kdy se využil výstup alfanumerických znaků k velmi hrubému výstupu grafické informace.
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
9
• Jazyky, které vznikly poněkud později (jako BASIC) obsahovaly specifické grafické příkazy zaměřený pouze na určitý typ hardware a takto vytvořené programy proto neumožňovaly přenositelnost na jinou hardwarovou základnu. • Velké výpočetní systémy, pokud byly určené i pro grafickou práci, obsahovaly tzv. grafické podsystémy tvořené specializovaným hardwarem i softwarem. Jednalo se obvykle o špičková grafická zařízení, která svými možnostmi výrazně překračovala i současný standard tvořený PC mikropočítači – např. již před 20-30 lety umožňovala práci s rastrovými zařízeními s rozlišením 4096 × 4096 a lepším. Vazba na klasické programovací jazyky se realizovala pomocí techniky podprogramů, tj. grafický podsystém se uživateli jevil jako specializovaný grafický program, kde jednotlivé činnosti bylo možno vyvolat formalizovaným voláním jednotlivých podprogramů. Tento postup umožňoval do běžných programovacích jazyků zahrnout i velmi silné grafické programovací prostředky – např. jedním příkazem s řadou volitelných parametrů bylo možno nakreslit osy grafu, navrhnout jejich dělení a popsat je. Pokračovateli tohoto směru vývoje jsou současné graficky orientované pracovní stanice. • Na běžných pracovištím se pro grafické výstupy v současné době využívají standardní PC mikropočítače. Zde se po stránce hardwarové opíráme o standardní monitory a tiskárny, resp. zapisovače. Pro jejich ovládání se rozšiřuje soubor instrukcí programovacího jazyka (např. místo jazyka Pascal se používal jazyk TurboPascal, apod.), případně jsou do překladače zahrnuty moduly pro grafickou práci (Compaq FORTRAN, atd.). Úroveň práce s takovým softwarem se velmi liší podle použitého programovacího jazyka – na jedné straně můžete mít pouze základní příkazy typu „nakresli bodÿ, „zvol barvuÿ, „ukaž barvu zvoleného boduÿ a „vymaž obrazovkuÿ (např. PUTPIXEL, SETCOLOR, GETPIXEL, CLEARSCREEN) a na druhé straně sadu silných programových prostředků pro vytváření kvalitních grafických výstupů. Pokud bychom chtěli ovládat grafické možnosti našeho počítače co nejdokonaleji, bylo by vhodné seznámit se i s organizací ukládání grafických dat do paměti mikropočítače a případně umět naprogramovat některé operace až na úrovni strojového kódu mikroprocesoru, resp. assembleru. To však překračuje rámec tohoto studijního textu.
• Práce se specializovaným softwarem, buď přímo zaměřeným na grafiku (např. program ORIGIN) nebo softwarem, kde práce s grafikou je součástí většího celku (např. systém MATLAB). V poslední době stoupá počet uživatelů, kteří dávají přednost poslední možnosti, neboť jim umožňuje získávat skutečně kvalitní grafické výstupy s nejmenší námahou. V tomto případě program produkující data, ať již jako výstup počítačového modelování nebo zpracování výsledků měření, tato data ukládá jako soubory na vnější médium (pevný disk, CD-ROM, apod.) a teprve po ukončení tohoto výpočtu je zavolán specializovaný program, který vytvořená data načte a zobrazí.
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
10
V tomto případě se počítačový stejně jako teoretický nebo experimentální fyzik stává pouhým uživatelem, který vůbec nemusí umět programovat a hlavně nemusí detailně znát použité grafické algoritmy. Proto též Počítačová grafika jako jeden ze směrů počítačové fyziky ztrácí na důležitosti. K tomu, aby uživatel volil správné metody z nabídky, kterou specializovaný software zobrazí, je však vhodné znát aspoň základní principy nejdůležitějších grafickým algoritmů. Pro kvalifikovanější uživatele, kterými by měli být počítačoví fyzici magisterské i bakalářské úrovně, se jako optimální jeví kombinace dvou postupů: – v průběhu hlavního výpočtu (simulace nebo analýzy experimentálních dat) naprogramovat zjednodušený grafický výstup na obrazovku s cílem nejprve ladění programu a později sledování jeho běhu, a – data vygenerovaná hlavním programem a uložená na médiu vykreslit dodatečně kvalitním specializovaným softwarem. Při tomto postupu je vhodné používat pro vytváření hlavního programu programovací jazyk umožňující aspoň základní grafické operace (což jsou v současné době překladače všech hlavních programovacích jazyků používaných na PC mikropočítačích) a pak mít v laboratoři patřičný graficky orientovaný software. Protože však obvykle grafický výstup hlavní výpočet zpomaluje, je vhodné mít v programu možnost volitelného zapnutí nebo vypnutí grafického výstupu. I při tomto systému práce je ovšem vhodné znát základní principy nejdůležitějších grafických algoritmů.
2.2.2
Souřadné systémy
I když grafických programů je velký počet, postupně se ustálilo pouze několik základních systémů, v nichž jsou uváděny souřadnice kreslených grafických objektů. Důvodem je, že jediný typ souřadného systému by výrazně zjednodušil přejímání programů a algoritmů. Ve skutečnosti takové systémy existují tři: • Reálné souřadnice (anglicky „world coordinatesÿ) Při zobrazování pomocí těchto souřadnic jsou konvence zcela stejné jako při kreslení grafů v matematice – osa x směřuje zleva doprava, osa y zdola nahoru a obě souřadnice pracují s racionálními čísly. Souřadnice proto přesně odpovídají zobrazovanému útvaru. Při tomto způsobu práce je uživatel zcela odizolován od grafického hardwaru, na kterém bude výsledný obraz prezentován, a proto tento systém zaručuje nejlepší přenositelnost programů. • Souřadnice zařízení (angl. „device coordinatesÿ) Tyto souřadnice jsou celočíselné. Osa x směřuje zleva doprava, osa y však shora dolů. Počátek (0, 0) je umístěn v levém horním rohu. Velikost kroku na obou osách je volena tak, aby odpovídala konkrétnímu grafickému zařízení, obvykle displeji nebo zapisovači. Např. pro rozlišení obrazovky VGA má pravý dolní bod souřadnice (639, 479). Práce s těmito souřadnicemi je nejjednodušší, přenositelnost je ale minimální.
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
11
• Normalizované souřadnice (angl. „normalized device coordinatesÿ) Obvykle odpovídají souřadnicím zařízení, tj. počátek soustavy souřadné je umístěn vlevo nahoře a osa y je orientována směrem dolů. Pracují však s idealizovaným grafickým zařízením s mimořádně velkým počtem kroků, např. v rozmezí 0 až 32737. Jedná se o kompromis mezi předchozími systémy i s kompromisními výhodami a nevýhodami. Každý z těchto souřadných systémů má své výhody a nevýhody. Při práci s reálnými souřadnicemi se nejjednodušším způsobem do programu přenášejí data z matematického modelu. V souřadnicích zařízení se velmi jednoduše naprogramuje grafický výstup na obrazovku, neboť se pracuje přímo s hodnotami pixelů. Problém však nastává tehdy, když chceme změnit grafické rozlišení, to pak znamená všechny číselné hodnoty v programu změnit. Proto se často jako kompromis sice pracuje se souřadnicemi zařízení, maximální hodnoty na obou osách se však uvádějí ve formě konstant xM ax a yM ax a délkové údaje se udávají v zaokrouhlených zlomcích těchto konstant. Změna grafického rozlišení se pak provede prostou změnou hodnot konstant xM ax a yM ax. Jak již bylo výše uvedeno, v tomto případě je třeba dát pozor na grafické rozlišení 1280 × 1024, které nezachovává poměr 4 : 3. Pokud by toto rozlišení bylo výchozí a nebo by se na něj grafický výstup převáděl, dojde k deformaci kružnic na elipsy a k dalším zkreslením.
Normalizované přístrojové souřadnice mají výhodu v tom, že jsou přístrojově nezávislé. Pokud se grafický výstup naprogramuje v nich, převod na skutečné přístrojové souřadnice se snadno provede vydělením vhodnou konstantou a vzhledem k velkým hodnotám koncových bodů normalizovaných souřadnic bude zaokrouhlovací chyba při dělení dostatečně malá. Tyto souřadnice se proto nejčastěji používají k přenosu obrazu mezi jednotlivými grafickými zařízeními. Pokud sami grafický program vytváříme, je volba vhodného souřadného systému plně v naší pravomoci – pouze se nedoporučuje zavádět další typy souřadnic než tři výše uvedené. Pokud však budeme pracovat s komerčním grafickým programem, je nutno zjistit, jaký souřadný systém používá a podřídit se mu.
2.3
Vybrané algoritmy počítačové grafiky
V této sekci se seznámíme se základními principy nejdůležitějších grafických algoritmů. Tento výčet nebude zdaleka úplný a též výběr algoritmů bude podřízen oblasti, v nichž čtenáři nejčastěji budou vlastní programy vytvářet a komerční programy používat. Úplnější informace lze nalézt v literatuře, např. [1], [2], [3]. S jedinou výjimkou má popis algoritmů uvedený v této kapitole význam pouze informativní, aby čtenář mohl správně používat knihovní programy, které při své budoucí práci bude pravděpodobně využívat.
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
2.3.1
12
Společné problémy rastrové grafiky
Ať již budeme na počítači vykreslovat úsečky, kružnice nebo složitější křivky, narazíme při použití rastrového výstupního zařízení na společné problémy, proto si je souhrnně uvedeme úvodem (s některými problémy bychom se setkali i při použití vektorových výstupních zařízení). Prohlédněme si obr. 2.2. Na tomto obrázku vidíme dva problémy současně – při kreslení úsečky s jiným sklonem než 0, ±45 a ±90 stupňů se použitý rastr projeví schodovitostí výsledné čáry, a – výraznost nakreslené čáry závisí na jejím sklonu.
Obrázek 2.2: Zobrazení úseček různého sklonu v rastrové grafice. První problém se zmírňuje při použití vyššího grafického rozlišení, kdy „schodyÿ jsou menší a zvětšující tloušťka čáry je lépe kompenzuje (viz dále). Druhý problém se s přechodem k vyššímu grafickému rozlišení nezmění a pokud nám vadí, musíme jej redukovat v rámci technických možností našeho výstupního zařízení. Např. při výstupu na displeji (obrazovce nebo LCD) budeme body ležící blíže u sebe vykreslovat vhodným stupněm šedi místo bíle (vzdálenost bodů √ ležících na vodorovné úsečce a úsečce vedené pod úhlem 45 stupňu je v poměru 1 : 2, proto pixely na vodorovné nebo svislé úsečce „rozsvítímeÿ pouze na 70 % maximálního jasu). Při práci na tiskárně často nemáme možnost měnit barvu bodu, můžeme ale někdy měnit v určitých mezích jeho velikost a tak dosáhneme stejného efektu.
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
13
Další problémy se projeví, pokud budeme kreslit přerušovanou křivku – čárkovanou, čerchovanou, apod. Nejjednodušší je programovat přerušovanou kresbu čáry pomocí osy x, vždy pro určité hodnoty x křivku kreslit a pro určité kresbu přerušit. Pokud ale budeme kreslit přímky různého sklonu nebo křivku, jejíž sklon se mění, dostaneme výsledky podobné jako jsou znázorněny na obr. 2.3.
Obrázek 2.3: Kreslení přerušované křivky v rastrové grafice. Na tomtéž obrázku vidíme ještě jeden problém. V obrázku jsou bílými kroužky znázorněny počáteční a koncové body křivek. Je zřejmé, že mechanické přepínání kreslení podle x-ové souřadnice pixelu většinou povede na situaci, kdy křivka neobsahuje koncový bod. Profesionální software tuto situaci řeší tak, že první a poslední úsečka čárkované křivky má jinou délku než vnitřní úsečky tak, aby koncový bod byl do křivky zahrnut. V případě napsání vlastního programu pro kreslení přerušovaných křivek musí uživatel tuto situaci zvládnout sám. Dalším z problémů, o kterém je nutno se zmínit, se týká tloušťky čar. Dříve, kdy se pracovalo s malým rozlišením (typicky VGA), měly všechny křivky tloušťku jednoho pixelu a žádné problémy tohoto typu nenastávaly. Současný hardware však umožňuje pracovat s rozlišením podstatně vyšším – řádu 1600 × 1200 a větším – a z důvodu omezení schodovitosti rastrové grafiky je vhodné toto rozlišení plně využít. V závislosti na použitém software, např. ve FORTRANu, je možno i na obrazovce pracovat s vyšším rozlišením, než je maximální rozlišení použité grafické karty, a displej chápat pouze jako okno do výsledného obrazu. To lze dokumentovat i v našem studijním textu, kde obrázky 2.3 až 2.5 byly vytvořeny s rozlišením 2000 × 1000 pixelů. Též v kapitole věnované zpracování obrazu jsou převzaty některé ilustrace z vědeckých publikací a je možno srovnat, jak vypadal grafický výstup před deseti lety a nyní (porovnejte např. obr. 3.13 až 3.15 navzájem i s dalšími obrázky v tomto studijním textu). Při velkém rozlišení však čára tloušťky jeden pixel je příliš slabá a je třeba ji odpovídajícím způsobem rozšířit – obvykle na 2-3 pixely pro osy grafu a 3-5 pixelů pro vlastní
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
14
Obrázek 2.4: Kreslení čar větší tloušťky (zde 4 pixely) v rastrové grafice. křivky v grafu. Při tom ale pozorujeme další jev – efektivní tloušťka čáry opět závisí na jejím sklonu, viz obr. 2.4. Na druhé straně je v tomto obrázku vidět, že zvětšením tloušťky čáry skutečně omezujeme její schodovitost. Toto chování je zřejmé z přímky, která prochází obrázkem úhlopříčně. Na obrázku obr. 2.4 je ale vidět ještě jeden problém, který přináší zvýšení tloušťky křivek, a to jejich zakončení. U křivek vodorovných a svislých je použité „pravoúhléÿ řešení vhodné. Naproti tomu u šikmých čar by mělo být jejich zakončení buď kolmé na směr čáry a nebo zaoblené. Profesionální grafické systémy tyto a další možnosti nabízejí, zatímco uživatel samostatně vytvářející grafický výstup se musí o odpovídající řešení postarat ve svém programu sám.
Obrázek 2.5: Správné kreslení křivky tloušťky 5 pixelů v rastrové grafice.
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
15
Na obr. 2.5 je znázorněno správné řešení, kdy tloušťka křivky je automaticky regulována v závislosti na jejím sklonu a též zakončení křivky je zaoblené. Proti reálné bitmapě jsou tento i předchozí obrázky 10× zvětšené, proto ve skutečnosti budou mít tyto křivky rozumnou tloušťku i hladkost. Vektorová grafika musí řešit podobné problémy, ale jinými technickými prostředky. Zapisovač má tloušťku čáry pevně danou (a malou), tlustou křivku proto vytvoříme soustavou rovnoběžných čar, které též musíme odpovídajícím způsobem zakončit.
2.3.2
Kreslení úsečky
Jako základní grafický prvek si uvedeme úsečku, protože většina složitějších obrazů kresbu úsečky ve svém algoritmu využívají. Úsečku můžeme kreslit pod libovolným úhlem, jak to problém vyžaduje. Pokud však máme možnost volby úhlu čáry, což nastává např. při vyplňování větších rovinných útvarů barvou, dáváme přednost vodorovnému směru, neboť kresba vodorovných čar probíhá výrazně rychleji než při jiném sklonu. Zdůvodnění leží v organizaci ukládání grafických dat do paměti mikropočítače. Úsečku budeme mít zadanou souřadnicemi počátečního (x1 , y1 ) a koncového bodu (x2 , y2 ). Z těchto souřadnic můžeme jednoduše spočítat směrnici úsečky a a úsek na ose y, b, a úsečku popsat parametricky y = ax + b. (2.1) Nejjednodušší algoritmus, který nás napadne, by mohl mít tvar: – začni od hodnoty x1 – zvětši x o jedničku a podle vztahu (2.1) spočítej hodnotu y – postup opakuj tak dlouho, pokud se x nebude rovnat x2 . Pokud bychom takový algoritmus skutečně naprogramovali, dočkali bychom se dvou nepříjemných překvapení – pro některé kombinace okrajových bodů by kresba úsečky byla velmi nepřesná, protože úsečka by obsahovala příliš málo bodů (obr. 2.6), a program by byl též zbytečně pomalý. Zdůvodnění prvního problému je v použití rastrové grafiky. Zatímco u vektorové grafiky by byl počáteční a koncový bod propojen plnou čarou, v případě rastrové grafiky se kreslí pouze body, které jsou v odpovídajících uzlech rastru. Pokud by směrnice a ležela v rozmezí h0, 1i, nakreslená úsečka by sice nebyla zcela hladká, ale její body by byly u sebe tak těsně jak použité rozlišení dovolí. Problémy nastávají při směrnici větší než jedna (i v absolutní hodnotě). Před použitím parametrického vyjádření pro kreslenou úsečku je proto nejprve třeba určit hodnotu směrnice a se vztahem (2.1) pracovat pouze tehdy, když kreslená čára bude mít směr v úhlovém rozmezí h315, 45i stupňů a nebo v rozmezí h135, 225i stupňů, tj. když směrnice a bude ležet v rozmezí h−1, 1i. V opačném případě vztah (2.1) převedeme na alternativní vyjádření b 1 (2.2) x= ·y− , a a
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
16
Obrázek 2.6: Kreslení úseček různých sklonů při použití řídící osy x. kde nová směrnice a1 již bude výše uvedenou podmínku splňovat. Znamená to ovšem též, že osou, jejíž hodnotu budeme postupně o jednotku zvyšovat, se stane osa y – této ose říkáme řídící osa. Podobnou transformaci musíme v případě rastrové grafiky provést vždy, když se bude jednat o kreslení čáry s příliš velkým sklonem. Pro zrychlení algoritmu pro kreslení úsečky je možno použít více postupů. Za jeden z nejsilnějších je považován tzv. Bresenhamův algoritmus, který si výjimečně v tomto studijním textu uvedeme celý [1]: procedure Usecka_Bresenham (x1, y1, x2, y2: integer); var deltax, deltay, konst1, konst2: integer; x, y, x_koncove, predikce: integer; begin if x1 > x2 then begin x := x2; y := y2; x_koncove := x1; end else begin x := x1; y := y1; x_koncove := x2; end; deltax := x_koncove - x; deltay := abs(y1 - y2); konst1 := 2 * deltay; konst2 := 2 * (deltay - deltax);
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
17
predikce := 2 * deltay - deltax; Put_pixel (x,y); while x < x_koncove do begin x := x + 1; if predikce < 0 then predikce := predikce + konst1 else begin y := y + 1; predikce := predikce + konst2 end; Put_pixel (x,y) end end; Tento algoritmus z roku 1965 [4] pracuje v celočíselné aritmetice a je optimalizován na rychlost. Popis jeho činnosti lze nalézt v práci [1] nebo [2].
2.3.3
Kreslení elipsy
Při kreslení kružnice, nebo obecněji elipsy, nastává situace podobná jako při kreslení úseček. Je třeba zvolit algoritmus jak odolný proti chybám rastrové grafiky tak dostatečně rychlý.
Obrázek 2.7: Kreslení kružnice v parametrickém vyjádření s řídící osou x. Prvním problémem je zvolit vhodné parametrické vyjádření odpovídající vztahům (2.1) a (2.2). V případě kružnice a elipsy se ale situace komplikuje v tom, že směrnice kreslené
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
18
čáry se mění. Z tohoto důvodu je nevhodné obvyklé parametrické vyjádření kružnice založené na vztahu (x − x0 )2 + (y − y0 )2 = R2 , kde (x0 , y0 ) jsou souřadnice počátku a R poloměr kružnice. Z tohoto vztahu sice jednoduše dostaneme parametrické vyjádření kružnice q
y = y0 ±
R2 − (x − x0 )2 ,
(2.3)
při jeho použití však bude kružnice v částech přiléhajících k ose x tvořena příliš malým počtem bodů – viz obr. 2.7. Stejná situace by nastala i pro elipsu. Situaci by trochu zlepšilo, kdybychom v průběhu kreslení kružnice přepínali postupně mezi řídícími osami x a y podle sklonu kreslené křivky. Nejlepším řešením však je přejít na parametrické vyjádření kružnice (nebo elipsy) založené na středovém úhlu ϕ, tj. pracovat v polárních souřadnicích. V tomto případě bude mít parametrické vyjádření kružnice tvar x = x0 + R cos ϕ y = y0 + R sin ϕ
(2.4)
a všechny body kružnice budou v jejím obraze rozloženy ekvidistantně. Jejich vzdálenost si můžeme (v rámci rastru daného grafického rozlišení) volit pomocí kroku úhlu ϕ – viz obr. 2.8.
Obrázek 2.8: Kreslení kružnice v parametrickém vyjádření s úhlem ϕ. Výpočet můžeme ještě zefektivnit, když si uvědomíme, že body na kružnici jsou rozloženy symetricky podle všech poloos. Pro kružnici proto stačí spočítat souřadnice jen jedné osminy potřebných bodů a ostatní získat postupnými permutacemi mezi hodnotami x, y, −x, −y – na obr. 2.8 jsou dopočítané body znázorněné bílými kroužky. Podobná situace platí i v případě elipsy s tím, že má nižší symetrii a proto musíme spočítat jednu čtvrtinu potřebných bodů.
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
19
Počet bodů, které při kreslení kružnice (nebo elipsy) musíme spočítat, bude záviset jednak na technice kreslení a dále na velikosti kružnice. Při obvyklém postupu spočítáme přiměřený počet bodů, které budeme chápat jako počáteční a koncové body úseček, kterými je spojíme. Pokud však budeme chtít nakreslit kružnici co nejpřesněji, spočítáme tolik bodů, aby to odpovídalo rastru displeje. Dále platí, že čím větší bude poloměr kružnice, tím více bodů musíme pro kreslení získat – např. na obr. 2.8 jsme pracovali se středovým úhlem 3 stupně. Naopak, pokud chceme kreslit velmi malé kroužky jako jsou např. značky v grafech, obvykle nám místo kružnice stačí nakreslit šesti- nebo osmiúhelník. Maximální efektivita kreslení kružnic a elips se v rastrové grafice docílí tak, že pro kreslení bodů určených parametrickým vyjádřením typu (2.4) se opět použije Bresenhamův algoritmus v modifikaci pro kreslení elipsy [1]. V případě použití vektorové grafiky je situace poněkud odlišná, neboť máme k disposici pouze možnost kreslení úseček. Pomocí vztahu (2.4) proto spočítáme dostatečný počet bodů, které pak úsečkami spojíme, tj. vždy nahrazujeme kružnici mnohoúhelníkem s velkým počtem vrcholů. Pokud chceme docílit vysoké přesnosti kreslení, musíme zvážit, zda požadujeme mnohoúhelník vepsaný, opsaný nebo kompromis mezi těmito dvěma případy (což je nejobvyklejší).
2.3.4
Vyplňování oblasti
Při vyplňování určité oblasti barvou je základní otázkou, jak bude zadána hranice této oblasti. Rozeznáváme dva základní případy: – hranice zadaná křivkami, které jsou popsány analyticky (obvykle úsečkami, částmi kuželoseček, apod.), a – hranice zadaná aspoň částečně křivkami, které umíme charakterizovat pouze pomocí souřadnic pixelů které je vytvářejí (např. čáry kreslené „od rukyÿ). V prvním případě problém převádíme na kreslení úseček a využíváme přitom faktu, že nejrychleji jsou kresleny úsečky vodorovné. Algoritmus, zvaný metoda řádkového rozkladu, můžeme stručně popsat takto: – najdi nejvyše položený bod vyplňované oblasti a urči jeho y-ovou souřadnici ymax – tímto bodem veď myšlenou přímku rovnoběžnou s osou x – spočti průsečíky této přímky s hranicemi oblasti (výsledkem mohou být jeden nebo dva body) – pokud dostaneš jeden bod, vyplň jej požadovanou barvou – pokud dostaneš dva body, považuj je za počáteční a koncový bod vodorovné úsečky a tuto úsečku požadovanou barvou nakresli – zmenši y o jedničku a postup opakuj – výpočet ukonči tehdy, až myšlená přímka nebude mít žádný průsečík s hranicí oblasti.
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
20
Obrázek 2.9: Vyplňování oblastí zadaných analyticky popsanými hranicemi. Pro oblast jednoduchého tvaru tento algoritmus povede k správnému cíli. Pro oblast složitějšího tvaru však může nastat případ, že budeme mít více než dva průsečíky (na obr. 2.9 vpravo budeme pro určité hodnoty y dostávat tři nebo čtyři průsečíky a pro poslední y dokonce několik desítek „průsečíkůÿ). Pak musíme v programu ošetřit, aby byla barvou vyplněna správná část oblasti. Detaily je možno opět nalézt v práci [1] nebo [3]. Při aplikaci metody řádkového rozkladu musíme dát velký pozor na formulaci problému, tj. na skutečný tvar vyplňované oblasti. Pokud se hraniční čáry nikde neprotínají, jako je tomu na obr. 2.9, je problém formulován jednoznačně. Pokud však hraniční křivky rozdělují oblast na více podoblastí, je třeba přesně určit, které její části máme ve skutečnosti za úkol vyplnit, a tomuto úkolu výše popsaný algoritmus přizpůsobit.
Druhým základním případem je vyplňování oblasti, která je definována pomocí seznamu hraničních pixelů, tj. kde hranice jsou nakreslené a ne zadané geometricky. Tuto skupinu algoritmů můžeme stručně charakterizovat těmito kroky: – prověř, zda je zadaná hranice úplná – najdi jeden vnitřní bod oblasti – z tohoto bodu obarvuj sousední body tak dlouho, pokud ve všech směrech nedojdeš k hranici. Prvním problémem, který podle tohoto schématu musíme vyřešit, je otestování úplnosti zavedení hranice, tj. její „nepropustnostÿ. Na čtvercové síti pixelů rozeznáváme dva případy – každý pixel může mít čtyři sousedy (pak se jedná o tzv. 4-spojitou hranici a 4-spojité vyplňování [1]) a nebo osm sousedů (8-spojité vyplňování) – viz obr. 2.10. Nejjednodušší algoritmus je nazývaný semínkové vyplňování (název pochází z představy, že původní vnitřní bod oblasti představuje rostlinu plevele, který se každým rokem vysemeňuje a šíří do okolí, tj. každým rokem vytváří 4 nebo 8 nových rostlin). Algoritmus má tvar:
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
21
Obrázek 2.10: Dvě varianty vytvoření neprostupné hranice – hranice 4-spojitá (vlevo) a hranice 8-spojitá (vpravo). – vezmi původně zadaný vnitřní bod oblasti a obarvi jej – urči jeho nejbližší sousedy (4 nebo 8 podle typu úlohy) – jeden z nově určených bodů vezmi jako nový výchozí bod a souřadnice ostatních (3 nebo 7 bodů) ulož na zásobník – otestuj, zda nový výchozí bod je vnitřním bodem oblasti nebo již bodem hranice – pokud je tento bod vnitřním bodem oblasti, obarvi jej a jdi na krok 2 – pokud je tento bod bodem hraničním, vezmi ze zásobníku souřadnice posledního úloženého bodu a jdi na krok 2 – pokud je zásobník prázdný, ukonči činnost. Tento algoritmus je velmi jednoduchý a za pomoci rekurentních vztahů je dobře programovatelný. Má však základní slabinu – velmi neefektivně pracuje s pamětí, takže je v praxi použitelný jen pro vyplňování velmi malých oblastí. Důvodem je, že nebere v úvahu již vybarvené body a velmi snadno proto dojde k přeplnění zásobníku. V praxi se proto výhradně používá jeho modifikace nazvaná řádkové semínkové vyplňování. Rozdíl proti původnímu semínkovému vyplňování je v tom, že ve vodorovném směru se „plevelÿ v jednom kroku šíří neomezeně, resp. až na hranici oblasti. Pak na zásobník ukládáme pouze jeden bod v jednom ze svislých směrů, zatímco v druhém směru pokračujeme ve vyplňování. Dále na zásobník případně ukládáme body charakterizující část oblasti, která je od hlavní oblasti oddělena hranicí – viz spodní přímka na obr. 2.9 vpravo. Díky těmto modifikacím původního semínkového vyplňování je algoritmus mnohem efektivnější, pracuje s vodorovnými úsečkami a se zásobníkem zachází mnohem úsporněji. Detaily jsou opět uvedeny v práci [1]. I při řešení problematiky vyplňování oblasti narážíme na problémy s rastrovou grafikou, tj. na problémy o kterých jsme hovořili v sekci 2.3.1. Jedná se o to, že při použití dostatečně jemného rastru bude sice vnitřek oblasti vyplněn barvou homogenně, schodovitost se ale projeví na okrajích oblasti, na jejích hranách.
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
22
Obrázek 2.11: Dvě techniky vyhlazování hran. Řešení tohoto problému závisí na možnostech použitého hardwaru. Dvě obvyklá řešení jsou znázorněna na obr. 2.11. Vlevo je uveden případ, kdy můžeme pracovat s více stupni šedi, což je obvyklé u displejů nebo barevných tiskáren, vpravo případ vhodný pro monochromatické tiskárny. Na příkladu z obrázku jsme předpokládali možnost zavedení vedle bílé a černé barvy ještě tří odstínů šedi, případně místo jednoho většího černého bodu tisknout čtyři černé body menší. To umožnilo zjemňovat hranu po 25 % a tak např. místo rozlišení 600 dpi kreslit hranu s efektivním rozlišením 1200 dpi. Pro srovnání si můžeme uvést i obrázek kreslený v původním rozlišení – viz obr. 2.12.
Obrázek 2.12: Kresba z obr. 2.11 bez dodatečného zjemnění. Výše uvedené postupy vyhlazování hran jsou založené na čistě geometrických úvahách – kolik procent plochy pixelu patří do vybarvované oblasti, takovým stupněm šedi onen pixel nakreslíme. Ve vývoji jsou ale zcela odlišné algoritmy založené na vlastnostech lidského zraku, resp. na jeho nedokonalosti, podobné různým optickým klamům a hříčkám. První dosažené výsledky nasvědčují tomu, že vhodně navržená hrana, i když z hlediska geometrického méně přesná, se může jevit pozorovateli mnohem dokonalejší než matematicky přesné řešení.
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
23
Poznámka – tento postup na jemnější vykreslování hran větších ploch lze s úspěchem využít i pro vyhlazování okrajů tlustých čar (sekce 2.3.1). Postupy pro vyhlazování hran jsou standardní součástí řady systémů – ať již laserových tiskáren tak i operačních systémů. Na obr. 2.13 je příklad, jak zjemňování hran písmen řeší operační systém Windows, případně pod ním spouštěné programy jako Painbrush (Malování), Compaq FORTRAN, apod.
Obrázek 2.13: Výstup textu na displeji nebo tiskárně, nahoře – v původním rozlišení, dole – s vyhlazováním hran založeném na technice více stupňů šedi.
2.3.5
Kreslení prostorových objektů
Hlavním úkolem počítačové grafiky na fyzikálním pracovišti je presentovat v přehledné formě výsledky měření nebo počítačového modelování. Obvykle pracujeme s grafy funkcí jedné proměnné, jako jsou např. volt-ampérové charakteristiky, optická spektra, frekvenční závislosti, apod. Pro tyto úkoly bychom již měli mít dostatečné znalosti a navíc v této oblasti existuje velké množství více či méně dostupného softwaru. Občas však narazíme na úkol složitější, když námi studovaná závislost je ovlivňována dalším fyzikálním parametrem. V tomto případě s jedním grafem znázorňujícím funkci jedné proměnné nevystačíme. Řešením je místo jednoho grafu shrnout naše výsledky do celé série obrázků pro různé hodnoty parametru. Podstatně přehlednější ale je, když všechna data vykreslíme do stejného obrázku. Lze sice pracovat i „plošnýmÿ grafem obsahujícím několik křivek odpovídajících jednotlivým hodnotám parametru, vhodnější však je připravit třírozměrný graf. Jeho příklad je uveden na obr. 2.14. Dalšími problémy, které jsou obvyklejší spíše v matematice, je vykreslování průběhu funkcí více proměnných či obecně zobrazování prostorových objektů. Pro tyto účely lze v literatuře nalézt větší množství algoritmů a existují i vhodné softwarové balíky, které programy tohoto typu obsahují. Avšak při vykreslování souborů dat
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
24
Obrázek 2.14: Maxwellovo rozdělení rychlostí elektronů N (v relativních jednotkách, tj. normalizováno na jednotkovou plochu) v závislosti na jejich rychlosti v a teplotě T . ve fyzikální laboratoři občas míváme specializované požadavky na formát grafického výstupu (pro dizertace, vědecké publikace, apod.), které komerční software nedokáže splnit. Je proto vhodné znát 1-2 algoritmy pro vykreslování složitějších obrazů, abychom si v případě potřeby mohli potřebný graf připravit sami. Když se podíváme na výše uvedený obrázek, vidíme, že problém spočívá v tom, že existují oblasti grafu, které jsou pro nás „neviditelnéÿ, tj. které nechceme vykreslovat, abychom získali dojem prostorovosti a obrázek byl přehlednější. Relativně jednoduchý algoritmus, který takové požadavky splňuje, se nazývá malířův algoritmus. Jeho podstata spočívá v tom, že jednotlivé křivky kreslíme odzadu a chápeme je jako plochy ohraničené ze spodní strany osou x. Tyto plochy vyplníme neprůhlednou (bílou) barvou a tak zakryjeme čáry v pozadí, které nemají být vidět. Pokud chceme nějaké čáry ponechat i v zakryté části obrazu, vykreslíme je až nakonec. Nevýhoda tohoto postupu spočívá v nižší efektivitě práce, neboť vykreslujeme i čáry, které budou později smazány. Výhodou naopak je jednoduchost připravovaného programu, což běžný uživatel ocení. Na obr. 2.14 jsme takto zakryli části křivek Maxwellova rozdělení a šikmé přímky v rovině (v, T ) odpovídající jednotlivým hodnotám rychlosti. Naopak přímky odpovídající dělení teplotní osy jsme chtěli v obrázku ponechat, proto byly kresleny dodatečně. Pro jednoduché úkoly odpovídající běžným požadavkům fyziků je malířův algoritmus plně vyhovující. Pokud však chceme přesně zobrazit průběhy funkcí dvou proměnných
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
25
apod., použijeme raději silnější algoritmus. Ten se nazývá metoda plovoucího horizontu. Metoda plovoucího horizontu vykresluje na rozdíl od malířova algoritmu obraz odpředu a při dalším kreslení se berou v úvahu části obrazu, které jsou již zakryté. K tomuto účelu se používají dva horizonty, horní a dolní, což jsou celočíselné vektory s délkou, která přibližně odpovídá šířce obrazu (v pixelech). Základní kroky algoritmu metody plovoucího horizontu jsou následující: – na počátku naplň horní horizont hodnotami menšími a spodní horizont hodnotami většími než jakákoliv souřadnice budoucího obrazu – při vykreslování první křivky ukládej její y-ové souřadnice a vertikální umístění osy x v obraze do příslušných horizontů – větší hodnotu do horního a menší do spodního, jako index buňky příslušného vektoru vezmi x-ovou souřadnici vykreslovaného bodu křivky – tím do vektorů-horizontů bude uložena informace o té části obrazu, která je již zakryta první křivkou, resp. plochou ohraničenou první křivkou a osou x – při kreslení druhé a dalších křivek se nejprve podívej do obou horizontů a příslušný pixel vykresli pouze tehdy, pokud nepadne do zakryté části obrazu – v případě, že nový bod bude vykreslen, oprav hodnotu příslušného prvku vektoruhorizontu – při vykreslování dalších křivek ber v úvahu posunutí počátku této křivky – index příslušného prvku horizontu změň o rozdíl x-ových souřadnic, rozdíl y-ových souřadnic připočti k hodnotě, kterou do toho prvku ukládáš – po vykreslení všech křivek grafu v případě potřeby dokresli do zakryté části obrazu informaci, kterou tam chceš mít. Tato varianta metody plovoucího horizontu odpovídá situaci vykreslené v obr. 2.14, kde řezy jsou rovnoběžné s osou x. To též znamená, že při kresbě křivky je řídící osou (viz sekce 2.3.2) osa x. To se záporně projevilo u křivky Maxwellova rozdělení pro teplotu 1000 K. V experimentální fyzice, kdy úkolem je vykreslit graf funkce naměřené pro předem daný počet x-ových souřadnic, další informace pro kreslení nemáme k disposici. Řešením v tomto případě je propojit naměřené body úsečkami. V matematice nebo teoretické fyzice, kde úkolem často bývá nakreslení funkcí dvou proměnných složitého průběhu, je však vhodné algoritmus metody plovoucího horizontu rozšířit tak, aby mohl pracovat s oběma řídícími osami. Informace lze nalézt např. v práci [1].
2.4
Vizualizace
Počítačová grafika je součástí počítačové fyziky od samého začátku. V poslední době však se od ní počala oddělovat nová oblast – vizualizace velkých souborů dat. Důvodem je, že existují oblasti fyziky i dalších vědních disciplín, které produkují obrovské množství informace jež je nutno přehledně presentovat a to někdy dokonce velmi rychle.
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
26
Příkladem může být fyzika vysokoteplotního plazmatu. Při experimentech na zařízeních typu tokamak (výbojová komora ve tvaru toroidu s magnetickým polem, které udržuje pulzně generované plazma ve středu komory, a doplněná diagnostikou vznikajícího horkého plazmatu) je při jednom pulzu získáváno řádově 100 MB dat. Tato data je třeba on-line zpracovat (resp. předzpracovat, protože výsledné vyhodnocení se provádí dodatečně) a výsledky předat experimentátorům, aby mohli modifikovat parametry dalšího pulzu, přičemž perioda pulzů je řádu 10 za hodinu. Takové parametry kladou mimořádné požadavky jak na použitý hardware tak i na software určený pro předzpracování naměřených dat a na presentování výsledků.
V současné době již existuje velké množství technik použitelných k tomuto účelu. Jejich podrobný popis přesahuje zaměření našeho studijního textu, neboť jsou obsahem specializovaných přednášek v magisterském a doktorském studiu. Přehledně však můžeme uvést, že při presentování velkých souborů dat tak, aby je uživatel co nejlépe vstřebal, se používají moderní techniky založené na vytváření dynamické obrazové informace (jednotlivé obrázky vytvořené postupy počítačové grafiky se řadí do sekvencí, které se promítají velkou rychlostí tak že vytvářejí jakýsi film), pracuje se s barvou, využívají se postupy pseudo-prostorové presentace popsané v sekci 2.3.5 a v poslední době se začíná pracovat s pravou prostorovou presentací, ať již založenou na technice steroskopie (polarizační brýle, přepínání obrazů pro obě oči, atd.) nebo pomocí holografie. Jako příklad relativně jednoduché dynamické vizualizace, jež je v silách uživatele s běžnými znalostmi, si uvedeme postup použitý v práci [5]. Úkolem bylo studovat procesy, které se odehrávají na rozhraní nízkoteplotního plazmatu a do něho vnořené pevné látky. Plazma se skládá (kromě jiného) ze směsi záporných elektronů a kladně, případně i záporně, nabitých iontů. Tyto nosiče náboje se velmi odlišují svou hmotností, což vede na zajímavé efekty jak ve statickém tak zejména dynamickém režimu. Pokud na vnořenou pevnou látku (kov) přiložíme vůči plazmatu určité předpětí, budou se náboje jednoho znaménka přitahovat a druhého odpuzovat a to tak dlouho, dokud se v okolí pevné látky nevytvoří přechodová vrstva, která vliv nabitého povrchu odstíní. Tento proces probíhá i ve statickém režimu a vede na určitý ustálený stav soustavy plazma-pevná látka (jedná se o jakýsi ekvivalent známějšího p-n přechodu v polovodičích). V dynamickém režimu se pak výrazně projeví rozdílné hmotnosti nosičů náboje, což právě bylo cílem našeho studia. Z hlediska výpočetního se při studiu této problematiky využívají postupy popsané v prvním dílu našeho studijního textu, tj. různé techniky počítačového modelování. Konkrétně, pro řešení dynamiky procesů probíhajících při interakci plazma-pevná látka se jako optimální jeví hybridní částicový postup kombinující deterministickou metodu molekulární dynamiky a stochastickou metodu Monte Carlo. Je samozřejmě možné dynamiku procesu, při kterém na pevnou elektrodu vkládáme vůči plazmatu proměnné (skokové nebo střídavé) předpětí, popsat soustavou jednotlivých grafů pro různé čásové okamžiky, různá předpětí a různé složení plazmatu. Mnohem přehlednější a na zahraniční konferenci, kde byla práce [5] presentována,
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
27
i mnohem efektivnější bylo připravit velký počet jednotlivých časových „snímkůÿ a seskupit je do ucelené presentace – v našem případě měla 10 000 „okénekÿ, což vytvořilo několikaminutový film. Vliv rozdílného složení plazmatu na dynamiku tohoto procesu jsme ukázali tak, že na snímcích jsou současně znázorněny dva různé procesy a při jejich souběhu je nejlépe demonstrován vliv tohoto složení – viz obr. 2.15 (snížená kvalita zobrazení je způsobena tím, že bylo sejmuto během běhu vytvořeného „filmuÿ).
Obrázek 2.15: Znázornění dynamiky interakce plazma-pevná látka pro dvě složení plazmatu směsi argon/kyslík obsahujícího různé koncentrace záporných atomárních iontů kyslíku.
Druhý příklad prakticky použité vizualizace pochází z fyziky tenkých vrstev. Při růstu tenkých vrstev různých materiálů technikou vakuového napařování dopadají na podložku jednotlivé atomy, po podložce migrují a postupně vytvářejí shluky atomů. O tvaru těchto shluků rozhodují síly které působí mezi atomy nanášené látky a podložky, resp. mezi jednotlivými atomy nanášené látky. Z hlediska výpočetního se jedná o techniku částicového počítačového modelování, konkrétně metodu molekulární dynamiky. Jelikož síly působící na atomární úrovni jsou krátkodosahové a závisejí na typech konkrétních atomech, je potřeba k jejich určení použít buď experimentálně stanovené potenciály a nebo data získat z kvantově-mechanických výpočtů.
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
28
Pro analýzu vzniklých útvarů se používají různé techniky zpracování obrazu, o kterých budeme mluvit v další kapitole. Tyto postupy dokáží z fotografií tenkých vrstev získaných v elektronovém mikroskopu vyhodnotit koncentraci objektů a vytvořit statistiku velikostí jednotlivých objektů i jejich rozložení po podložce. Pokud však chceme získat názornou představu o průběhu studovaných procesů, je vhodné opět použít vizualizační postupy. Na následujících dvou obrázcích převzatých z práce [6] si ukážeme průběh počáteční fáze růstu tenkých vrstev – obr. 2.16, a fáze slévání objektů – obr. 2.17. V publikaci byly použity tyto dva statické obrazy, v případě potřeby je opět možno vytvořit celou sekvenci obrazů a demonstrovat i dynamiku celého procesu.
Obrázek 2.16: Znázornění počáteční fáze růstu kovové vrstvy na dielektrické podložce.
Obrázek 2.17: Znázornění koalescence dvou kovových ostrůvků na dielektrické podložce. Obrazy uvedeného typu si byl schopen zkušený uživatel naprogramovat i sám, vhodnější však je využít nějaký profesionální produkt – zde konkrétně jsme pracovali se systémem MATLAB. Třetí příklad použití vizualizace v praxi pochází z fyziky pevných látek. Zde je za určitých podmínek možno získat obrazy s více stupni šedi (obr. 2.18). Interpretace těchto stupňů šedi se v různých aplikacích liší – může to být skutečná barva, ale také výškový profil nebo nějaká fyzikální charakteristika povrchu (výstupní práce, koeficient sekundární emise, apod.).
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
29
Obrázek 2.18: Mnohoúrovňové obrazy v 256 stupních šedi. Naším úkolem je obraz zvýraznit tak, abychom mohli snadno odlišit dva podobné ale ne stejné obrazy a abychom v následném zpracování mohli vyhodnotit studovanou fyzikální vlastnost. Na zvážení uživatele je, zda obraz ponechá v původních stupních šedi, převede jej do umělé barevné škály a nebo zvýrazní nějakým jiným způsobem. Na obr. 2.19 je ukázána technika převodu víceúrovňového obrazu na obraz binární, ale s použitím vrstevnic, které následně umožňují snadnější zpracování různými gradientními metodami. Převzato z práce [7].
Obrázek 2.19: Informace z obr. 2.18 převedená do binární podoby.
2.5
Shrnutí
Počítačovou grafiku stačí zvládnout pouze na uživatelské úrovni, protože existuje velké množství dostatečně kvalitních grafických programů, jež můžeme pro vykreslení vlastních dat použít. Podstatné pojmy, které by měl čtenář tohoto studijního textu znát, jsou: • Rozdíly mezi rastrovou a vektorovou grafikou a jak se tyto rozdíly promítnou do grafických algoritmů
KAPITOLA 2. POČÍTAČOVÁ GRAFIKA A VIZUALIZACE
30
• Problémy při kreslení čar v rastrové grafice • Principy algoritmů pro: – kreslení úsečky – kreslení kružnice – vyplňování oblastí • Algoritmus pro kreslení grafu funkce dvou proměnných.
2.6
Problémy k řešení
Problémy určené k procvičování látky popsané v této kapitole buď formou samostudia nebo v rámci organizované výuky: 1. Napište program pro výpočet směrnice a volbu řídící osy při kreslení úseček pod různým sklonem. Vstupní data: souřadnice počátečního bodu (x1 , y1 ) a koncového bodu úsečky (x2 , y2 ). 2. Na základě výpisu programu analyzujte Bresenhamův algoritmus pro kreslení úsečky. 3. Napište program pro kreslení kružnice o zadaném středu a poloměru. 4. Napište program pro kreslení elipsy o zadaném středu a velikosti obou poloos. 5. Napište program pro vykreslení funkce dvou proměnných f (x, y).
Kapitola 3 ZPRACOVÁNÍ OBRAZU Jak jsme si uvedli dříve, při počítačové grafice vycházíme z numerických údajů, které jsou uloženy v počítači jako výsledek předchozích operací, a my tuto informaci chceme získat jako výstup v grafické podobě. Naproti tomu při zpracování obrazů jsou obrazová data vstupní informací pro další výpočty. Zpracování obrazů patří mezi velmi důležité směry počítačové fyziky, neboť obrazová informace je základní v mnoha oblastech fyziky i v dalších oborech vědy a techniky – biologii, medicíně, strojírenství, atd. Základem jsou různé techniky zpracování obrazu (anglicky se tato disciplína nazývá buď „Image processingÿ nebo „Image analysisÿ). Od této základní disciplíny se odvozují další techniky jako např. počítačová tomografie. Na rozdíl od počítačové grafiky, kde již byly vytvořeny kvalitní komerční programy či soubory programů a uživatel proto potřebuje pouze základní znalosti, v oblasti zpracování obrazu je situace zcela jiná. I zde sice existují společné problémy pro uživatele z nejrůznějších oblastí vědy a proto pro jejich řešení i zde vznikl kvalitní software na komerční úrovni (to se týká zejména oblasti předzpracování obrazu, např. odstraňování šumu). V dalších krocích se však potřeby různých uživatelů natolik liší, že komerční programové soubory nemohou na všechny tyto požadavky reagovat a uživatel si proto musí specializované programy napsat sám. Z tohoto důvodu patří metodika zpracování obrazu mezi důležité směry počítačové fyziky, který se navíc stále ještě vyvíjí. Tato důležitost je samozřejmě hodnocena různě uživateli z různých oblastí vědy a techniky.
3.1
Formulace problému
Při zpracování obrazu vycházíme ze základní obrazové informace, která je nejčastěji výsledkem pozorování – fotografie z dalekohledů při sledování astronomických objektů, fotografie z optických mikroskopů v biologii, medicíně, metalurgii nebo fyzice pevných látek, mikrofotografie z transmisních elektronových mikroskopů ve fyzice tenkých vrstev nebo metalurgii, obrazy z STM nebo AFM mikroskopů ve fyzice povrchů, apod. Dříve byly vstupní obrazy získávány pouze ve formě fotografií, v současnosti stále více do vědy pronikají techniky on-line přenášení obrazu – např. v digitálním mikroskopu je umístěn CCD čip, který přímo 31
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
32
poskytuje elektrický signál. Podobně v astronomii bývá dalekohled přímo propojen se CCD kamerou. Tyto moderní techniky však modifikují pouze první krok obrazové analýzy, zatímco další řetězec operací ovlivněn není. Cílem následné analýzy těchto obrazů je získání fyzikální, biologické a jiné informace, která je v těchto obrazech skryta. Úplná analýza obrazu proto musí obsahovat tyto kroky: 1. Digitalizace 2. Analýza nízké úrovně (a) Geometrické transformace (b) Filtrace (c) Binarizace (d) Rozpoznávání objektů 3. Analýza vysoké úrovně (a) Integrální informace (b) Informace o jednotlivých objektech (c) Informace o rozložení objektů 4. Získávání odborných informací. Při konkrétní práci v dané oblasti vědy a techniky není váha všech uvedených kroků stejná a některé z nich mohou případně zcela chybět. Dříve se kladl velký důraz na etapu digitalizace, protože ta ovlivňovala výsledek celé obrazové analýzy. S pokrokem výpočetní a měřící techniky tato etapa ustupuje do pozadí a v případě použití dalekohledů nebo mikroskopů s CCD prvkem se neuplatní vůbec. Analýza obrazu nízké úrovně se musí používat při studiu obrazové informace získávané v libovolné oblasti vědy. Proto též zde se nejvíce uplatňují připravené knihovní programy, buď specializované a nebo tvořící součást větších programových souborů – např. Image Processing Toolbox v systému MATLAB. V této etapě zpracování obrazu je typické, že základní operace se provádějí pomocí komerčního software a pouze určité specifické operace si uživatel naprogramuje sám. Podle požadavků jednotlivých oblastí vědy však některé kroky analýzy nízké úrovně mohou chybět. To se týká především binarizace, kde je převáděn obraz z více úrovní šedi na binární a přitom se část informace ztrácí. Existují aplikace, kdy základní hledaná informace je skryta právě v různých stupních šedi nebo v barvě a proto se etapa binarizace přeskakuje. Při analýze vysoké úrovně je rozdíl požadavků různých uživatelů již tak velký, že pouze výjimečně existují dostatečně silné programové soubory, které uživateli vyhovují. Proto zde se předpokládá, že podíl vlastní programátorské práce uživatele je již značný. V této etapě je též typické, že různí uživatelé užívají pouze vybrané metody ze souboru technik analýzy vysoké úrovně – např. pro lékaře vyšetřujícího složení krve pacienta je zcela bezpředmětné
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
33
využívat metody popisující rozložení objektů v obraze a klíčové pro něj jsou metody analyzující tvar jednotlivých krvinek. Naproti tomu ve fyzice tenkých vrstev jsou mimořádně podstatné informace o rozložení objektů po povrchu, neboť z tohoto rozložení lze usuzovat na fyzikální procesy probíhající při nukleaci a následném růstu vrstev. Část metod z oblasti zpracování obrazu lze vyvíjet a aplikovat počítačovými fyziky, zatímco u dalších metod je nezbytná přítomnost uživatele z konkrétní oblasti vědy nebo techniky. Podíl kvalifikované práce daného uživatele-fyzika, biologa, lékaře, atd. narůstá s postupem této analýzy a poslední etapa „4. Získávání odborných informacíÿ je již plně v rukou odborníků na příslušnou vědní disciplínu. V následujících sekcích se seznámíme s jednotlivými etapami analýzy obrazu nízké i vysoké úrovně. Podrobnější informace mohou čtenáři nalézt například v [2] a [8].
3.2 3.2.1
Metody zpracování obrazu nízké úrovně Digitalizace
Cílem kroku „Digitalizaceÿ je přenést studovaný obraz do počítače. Konkrétní průběh tohoto kroku závisí na analyzovaném obrazu a použité technice. V minulosti se pracovalo především s televizní technikou, která snímala fotografie na papírovém médiu, takže získaná digitální verze fotografie byla zatížena chybami použité optiky. V současné době se při digitalizaci fotografií používají především skenery, které tímto problémem netrpí. Pokud se však zpracovávají obrazy on-line, často se opět používá optická cesta přenosu, ať již digitální fotoaparáty nebo specializované kamery, a otázka chyb optického zobrazení se opět musí brát v úvahu. Důležitou otázkou je prostorové rozlišení při digitalizaci. V období používání televizních kamer pro digitalizaci fotografií bývalo typické rozlišení výsledného obrazu 256 × 256 až 512 × 512. V současnosti při používání digitálních fotoaparátů a kamer není problém získat prostorové rozlišení zdigitalizované fotografie v řádu 2000 × 2000 a více. A pokud použijeme nejobvyklejší cestu: fotografie + skener, můžeme dosáhnout i prostorového rozlišení podstatně většího. Například při zpracování fotografie velikosti A4, tj. přibližně 8” × 12”, dostaneme při rozlišení skeneru 600 dpi výsledné rozlišení fotografie 4800 × 7200 a při rozlišení 1200 dpi dokonce 9600 × 14400. Jaké rozlišení tedy máme při zpracování konkrétní fotografie obsahující fyzikální informaci zvolit? V úvahu musíme brát dva protichůdné požadavky: – Velikost nejmenších detailů, které potřebujeme analyzovat. Otázkou je, co potřebujeme na fotografii studovat. Pokud nám stačí informace, že v daném místě fotografie nějaký objekt existuje, stačí pracovat s rastrem zdigitalizované fotografie polovičním než je velikost nejmenšího objektu. Pokud však nás bude zajímat i informace o tvaru objektů, měl by krok digitalizace být aspoň o jeden řád menší než je velikost nejmenšího studovaného objektu.
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
34
– Požadavky na paměť počítače. Budeme-li předpokládat (což je při studiu fyzikálních objektů typické), že se pracuje s 256 úrovněmi šedi, zaujímá v paměti počítače 1 pixel velikost 1 bytu. Pak fotografie velikosti A4 zaujímá při rozlišení 600 dpi velikost téměř 35 MB a při rozlišení 1200 dpi více než 120 MB. Pokud si uvědomíme, že při zpracování digitalizované fotografie potřebujeme pracovní prostor o velikosti dvou- až osminásobku velikosti fotografie v závislosti na použité metodě zpracování, jedná se o extrémně velký požadavek na paměť počítače. Kromě těchto odhadů je třeba brát v úvahu fakt, že s růstem velikosti digitalizované fotografie roste i doba výpočtu, lineárně nebo i kvadraticky v závislosti na použitém algoritmu analýzy obrazu. Nesmíme též zapomenout na to, jakou technikou byl obraz připraven, abychom nepracovali s větším rozlišením než je informační hodnota studovaného obrazu. Proto je třeba rozlišení digitalizace pečlivě zvážit a nepracovat při zbytečně velkém rozlišení. Z tohoto důvodu se podle charakteru studovaných objektů obvykle nepoužívá maximální rozlišení skeneru, ale pracuje se s rozlišením v rozmezí 200 dpi až 600 dpi. Výsledkem etapy digitalizace bude matice uložená v paměti počítače a obsahující v jednotlivých prvcích hodnoty pixelů původního obrazu. Tuto matici budeme označovat A a pokud nebude řečeno jinak, její prvky aij budou mít velikost INTEGER(1). V případě on-line snímání obrazu se etapa digitalizace prakticky přeskakuje, ale obraz elektronicky přenesený do paměti počítače opět označíme stejným způsobem a v dalším zpracování proto již nebudeme zdroj dat rozlišovat.
3.2.2
Geometrické transformace
Tímto pojmem rozumíme změnu obsahu matice A, tj. digitalizované fotografie, s cílem odstranit chyby, které vznikly ve fázi snímání obrazu nebo jeho digitalizace. Typickou chybou ve fázi snímání je deformace obrazu vlivem nevhodné geometrie, např. snímání povrchu Země z družice obíhající Zemi v malé vzdálenosti nebo fotografování povrchu vzorku v laboratoři ze šikmého směru. Též sem patří korekce vad optického řetězce použitého buď při pořizování původní fotografie nebo při její digitalizaci. Jedná se o klasické chyby zobrazení známé z optiky. Metoda nápravy je ve všech případech stejná. Je třeba nalézt transformační funkci, kterou se matice A vynásobí s cílem získat původní nezkreslený obraz. Tuto transformační funkci buď vypočítáme ze známé geometrie problému a nebo ji najdeme empiricky. Pokud totiž budeme celou sérii obrazů snímat stejným způsobem (např. fotografování stínítka obsahujícího obrazovou informaci kamerou, která nemůže být umístěna kolmo proti stínítku), nasnímáme stejným způsobem testovací obraz tvořený pravoúhlou mřížkou a z deformace známého obrazu mřížky dopočítáme transformační funkci.
3.2.3
Filtrace
Pojmem „filtraceÿ rozumíme při zpracování obrazu změnu hodnoty konkrétního pixelu obrazu, tj. prvku aij digitalizované matice, na základě hodnot okolních pixelů. Tímto způ-
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
35
sobem se budeme snažit zmenšit chyby, které již ve snímaném obrazu byly, případně které se do něho dostaly ve fázi snímání. Typickým příkladem druhého typu chyb je změna obrazového signálu přenášeného kosmickými sondami při průzkumu Sluneční soustavy během cesty tohoto signálu k Zemi. Existuje celá řada počítačových technik, které se při filtraci používají. Tyto techniky se liší podle toho, jaký je charakter chyb a jakou metodu jejich redukce volíme (a částečně též podle toho, zda potlačujeme šum v signálu dvourozměrném nebo jednorozměrném). V teorii chyb známe celou řadu chyb (v experimentální fyzice většinou hovoříme o šumu, který se skládá s měřeným signálem). Z nich nejobvyklejší je Gaussovský šum, kde výsledný signál je tvořen původním signálem na nějž je superponován šum, přičemž četnost odchylek klesá v závislosti na amplitudě šumové složky podle Gaussovy rozdělovací funkce. V této dvouparametrické funkci N (µ, σ) je první parametr nulový a druhý určuje intenzitu šumu. Tento typ šumu je v přírodě nejběžnější, jak to odpovídá centrální limitní větě počtu pravděpodobnosti. Druhým nejčastějším typem šumu je tzv. šum pepř a sůl, kde původní hodnota signálu je zcela nahrazena hodnotou šumovou, která nabývá extrémních hodnot, tj. černé nebo bílé. Název tohoto typu šumu je právě odvozen od vzhledu zašuměného obrazu. Filtrace je nejdůležitější nízkoúrovňovou operací, která se aplikuje jak při zpracování obrazů tak především při zpracování jednorozměrných signálů, tj. výsledků fyzikálních měření. Bylo proto vyvinuto velké množství metod na potlačování šumu, které můžeme rozdělit do několika skupin: – klasická filtrace – morfologické operátory – fourierovská filtrace. Klasické metody filtrace Klasické filtrační metody byly původně vyvinuty pro potlačování šumu ve fyzikální laboratoři, nejčastěji při měření volt-ampérových charakteristik nebo optických spekter. Po převedení z jednorozměrné formulace na dvourozměrnou je můžeme použít i pro potlačování šumu v obrazech. Těchto metod je známo několik desítek a my si zde uvedeme ty nejrozšířenější: – Metoda koherentní sumace Tato metoda je založena na skutečnosti, že šum je s časem náhodně proměnný. Sečtemeli proto několik zašuměných signálů a zprůměrujeme je, užitečný signál se nezmění zatímco šum poklesne. Zmenšení šumu je přibližně rovné druhé odmocnině z počtu zprůměrovaných obrazů. Výhodou této metody je, že jako jediná potlačuje pouze šum a nedeformuje užitečný signál. Při jejím použití však musí být splněna podmínka koherentnosti, tj. že užitečný signál se během celého průběhu operace odšumování nesmí měnit.
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
36
– Metoda obyčejného průměrování V experimentální fyzice se tato metoda nazývá metodou klouzavého průměru. V našem dvourozměrném případě si zvolíme malou matici nazývanou maska nebo strukturní element, např. velikosti 5 × 5 – obr. 3.1. Prostřední bod této masky si označíme. Masku budeme postupně přikládat na jednotlivé pixely obrazu, tj. prvky matice A (tuto operaci nemůžeme provést na krajních řádcích a sloupcích matice). Spočítáme průměr ze všech hodnot prvků matice, které leží pod maskou, a touto hodnotou nahradíme původní prvek, který leží pod označeným bodem matice. Tuto operaci postupně opakujeme se všemi vnitřními prvky matice A. Přitom je nezbytné, abychom masku postupně kladli na původní matici A, zatímco změněné hodnoty budeme ukládat do nové matice A0 .
Obrázek 3.1: Symetrické masky velikosti 3 × 3, 5 × 5 a 7 × 7 pro filtraci obrazu. Tato metoda potlačuje šum úměrně velikosti masky, pro masku 3 × 3 klesne šum v obraze přibližně třikrát, atd. Současně však mění i užitečný signál. V jednorozměrném případě metoda klouzavého průměru snižuje velikost maxim a zaplňuje minima, což v případě zpracování obrazu odpovídá rozmazávání hran a zmenšování kontrastu celého obrazu. Vedle této základní slabiny se u metody obyčejného průměrování (stejně jako při původní metodě klouzavého průměru) projeví další drobný nedostatek – jelikož masku typu obr. 3.1 nedokážeme přiložit až k okraji matice A, bude několik prvních a posledních řádků a sloupců obrazu neodfiltrováno. – Metoda mediánového průměrování Další metody mají za cíl zachovat schopnost metody obyčejného průměrování potlačovat šum, ale přitom omezit její nežádoucí vliv na hrany objektů v obraze. Metoda mediánového průměrování to dělá pomocí statistické funkce medián. Začátek této metody je stejný jako v předchozím případě. Zvolíme masku (např. typu znázorněného na obr. 3.1), postupně ji klademe na všechny vnitřní prvky odšumované matice A a po zpracování hodnot jednotlivých prvků masky výslednou hodnotou nahradíme pixel obrazu pod označeným prvkem masky. Rozdíl je pouze v technice zpracování hodnot prvků masky. Místo obyčejného průměrování jako v předchozím případě, nyní tyto hod-
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
37
noty seřadíme podle velikosti a jako výsledek si vezmeme jejich medián (z prvního dílu studijního textu připomínáme, že medián z lichého počtu členů posloupnosti je číslo, které leží po jejím uspořádání přesně uprostřed). Tato metoda dostatečně potlačuje šum (poněkud méně než metoda obyčejného průměrování, ale dostatečně) a přitom podstatně méně rozmazává hrany objektů v obraze. – Metoda rotující masky Ještě silnější technikou je metoda rotující masky. Při ní musíme zvolit masku asymetrickou, např. masku velikosti 3 × 3, na které označíme místo středu jeden prvek v rohu masky – viz obr. 3.2. Tato asymetrie nám umožní na stejný prvek matice A uložit masku čtyřmi různými způsoby. Použití této metody je standardní – masku postupně pokládáme na všechny vnitřní prvky matice A, spočítáme novou hodnotu opravovaného prvku matice a tu uložíme do výsledné matice A0 . Novou hodnotu dostaneme tak, že masku aplikujeme všemi čtyřmi možnými způsoby, obyčejným průměrováním určíme čtyři možné korigované hodnoty a z nich použijeme tu, která se nejméně liší od od původní hodnoty opravovaného prvku matice A.
Obrázek 3.2: Masky pro filtraci obrazu metodou rotující masky. Tato metoda opět dostatečně potlačuje šum (stejně jako původní metoda obyčejného průměrování) a velmi málo rozmazává hrany. V jejím algoritmu je ale zabudován předpoklad, že šum má Gaussovskou povahu, tj. tato metoda je nepoužitelná pro odstraňování šumu typu pepř a sůl. – Další metody Existuje ještě celá řada dalších metod, které se od výše uvedených liší v tom, že se v nich používá vážené průměrování, což znamená, že při vytváření nové hodnoty pixelu se prvkům masky dávají různé váhy. Tyto váhy se volí tak, aby metoda byla optimalizována pro určitý typ problémů. Tímto způsobem lze masku přizpůsobit předpokládanému typu šumu a dokonce lze vytvořit masky, které hrany nerozmazávají ale naopak zvýrazňují. Extrémním případem jsou filtry, které odstraní z obrazu jakýkoliv neproměnný signál (což při zpracování obrazu znamená smazání ploch vyplněných jednou barvou) a ponechají pouze signál skokově se měnící, tj. hrany. Takový obraz se potom změní v perokresbu, kde na bílém pozadí jsou nakresleny hranice objektů. V matematice takové operaci odpovídá diferenciální operátor gradient. Více informací o těchto metodách najdete např. v [2], případně knihovně systému MATLAB – Image Processing Toolbox.
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
38
Ve všech výše uvedených metodách můžeme účinnost filtru regulovat. Standardní metodou regulace je volba velikosti masky, tj. 3 × 3, 5 × 5, 7 × 7 nebo i větší. Jelikož ale při zvětšování masky s rostoucí účinnosti potlačování šumu rostou i škodlivé účinky na hrany objektů, je vhodné při práci mít celý proces pod kontrolou. K tomu se lépe hodí druhá metoda, při které se pracuje s malou maskou např. 3 × 3 a celý proces se podle potřeby opakuje. Jelikož jakékoliv zásahy do obrazu jsou nebezpečné a hrozí nevratným poškozením analyzované obrazové informace, je vhodné mít originální matici A uloženou, pracovat na její kopií a teprve, když jsme si jisti výsledkem, originál nahradit změněnou maticí. Tato zásada platí nejen pro filtraci, ale i pro další morfologické operace. Další zásadou pro práci s obrazy, platnou opět nejen pro operace fitrace, je spoléhat více na zrak než na různá kritéria nastavená v programu. Je proto vhodné výsledky jednotlivých dílčích operací znázorňovat na obrazovce spolu s původní digitalizovaným obrazem a při rozhodování o dalším postupu se řídit vlastním zrakem. To platí například i o volbě stupně filtrace zmíněné výše, kdy proces opakované aplikace masky je vhodné naprogramovat interaktivně. Filtrace pomocí morfologických operátorů Další, méně často používanou ale efektivní, skupinou metod filtrace jsou tzv. morfologické operátory. Tyto operátory jsou dílčím výsledkem silné matematické teorie vzniklé v 70. letech zvané matematická morfologie (např. [9]). Tato teorie je využívána především při vytváření metod zpracování obrazů vysoké úrovně, proto se s ní seznámíme později. Morfologické operátory byly původně navrženy pro práci s již binarizovanými obrazy, tj. obrazy, kde objekty jsou znázorněny černou barvou a pozadí bílou. Později sice byla jejich definice rozšířena i na obrazy obsahující více úrovní šedi, my však zde budeme používat původní definici morfologických operátorů. Rozlišujeme více druhů morfologických operátorů: – Morfologické operátory 1. stupně 1. Dilatace 2. Eroze – Morfologické operátory 2. stupně 1. Otevření 2. Uzavření 3. Podmíněné operátory 4. . . . Při všech aplikacích morfologických operátorů opět využíváme strukturní element neboli masku. Na rozdíl od masek v klasické filtraci může mít tento strukturní element libovolný tvar a velikost, nejmenší strukturní element má proto rozměry 2 × 1. Jeden prvek
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
39
strukturního elementu si opět označíme. Při všech operacích postupujeme podobně jako při práci s maskou v metodách klasické filtrace – strukturní element postupně přikládáme na jednotlivé prvky původní matice A, hodnotu prvku této matice ležící pod označeným prvkem strukturního elementu změníme, změnu zapíšeme do druhé matice A0 a postup postupně opakujeme se všemi prvky matice A na které můžeme strukturní element položit. Při operaci dilatace označíme prvek matice A ležící pod označeným prvkem strukturního elementu černou barvou tehdy, když aspoň jeden prvek strukturního elementu leží na objektu (tj. na černém prvku matice A). Výsledkem operace dilatace proto bude zvětšení objektů úměrné velikosti strukturního elementu a poloze označeného prvku v něm.
Obrázek 3.3: Význam morfologického operátoru dilatace. Vlevo původní obraz, vpravo obraz po dilataci s použitím uvedeného strukturního elementu. Přidané pixely jsou označeny dvojitým okrajem.
Obrázek 3.4: Operace dilatace provedená na digitalizované matici stejné jako na obr. 3.3 s použitím většího strukturního elementu. Přidané pixely jsou opět označeny.
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
40
Na obr. 3.3 je znázorněn výsledek dilatace provedené pomocí malého strukturního elementu do obrázku zakresleného – je to nejmenší strukturní element, který lze vytvořit. Pro srovnání byla provedena stejná operace se stejným obrazem, ale jiným strukturním elementem. Výsledek je uveden na obr. 3.4. Při operaci eroze vyplníme prvek matice A ležící pod označeným prvkem strukturního elementu bílou barvou (tj. barvou pozadí) tehdy, když aspoň jeden prvek strukturního elementu leží mimo objekt (tj. na bílém prvku matice A). Výsledkem operace eroze bude zmenšování objektů a některé další efekty.
Obrázek 3.5: Význam morfologického operátoru eroze. Vlevo původní obraz, vpravo obraz po erozi. Přidané pixely jsou označeny bílými kroužky. Oba morfologické operátory prvního stupně by mohly najít při zpracování obrazu své uplatnění: – Dilatace se používá tehdy, obsahují-li objekty dutiny, které chceme zaplnit, případně chceme-li odstranit bílé body na objektech způsobené šumem typu pepř a sůl. Na obr. 3.3 vidíme, že došlo k zaplnění jen malých dutin, zatímco na obr. 3.4 se podařilo zaplnit i dutinu větší. Závěr z tohoto pozorování je, že velikost strukturního elementu musíme volit úměrnou velikosti maximální dutiny, kterou chceme odstranit. – Eroze se používá při opačných problémech – odstranění černých bodů šumu typu pepř a sůl na pozadí, resp. odstranění jiných malých poruch či neúplných objektů. Opět velikost strukturního elementu určuje, jak velkou poruchu na pozadí dokážeme odstranit. Vedle těchto pozitivních funkcí mají ale oba operátory základní nevýhodu v tom, že mění užitečnou informaci v obrazu, tj. zvětšují resp. zmenšují objekty, případně mění i jejich počet – eroze odstraňuje velmi malé objekty a dilatace může dva objekty spojit v jeden. Obě tyto procesy jsou přítomné v obr. 3.4 a 3.5. Jelikož tyto nevýhody při praktickém použití převažují, používají se proto v převážné míře morfologické operátory druhého stupně.
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
41
Morfologický operátor druhého stupně otevření je definován jako uspořádaná dvojice operátorů eroze + dilatace, která pracuje se stejným strukturním elementem. Při aplikaci tohoto operátoru se provedou změny obrazu stejné jako při původním operátoru eroze, tj. odstraní se poruchy z pozadí, a následnou operací dilatace se velikost objektů vrátí na původní hodnotu. Název operátoru pochází z typické aplikace, kdy na astronomické fotografii se dvě blízké hvězdy vlivem procesů ve fotografické emulzi propojí krčky, které ve skutečnosti neexistující. Při správném použití tohoto operátoru se krčky odstraní a obrazy hvězd zůstanou nedotčené (a poruchy na pozadí budou též odstraněny) – viz obr. 3.6.
Obrázek 3.6: Morfologický operátor otevření. Vlevo původní obraz, uprostřed obraz po provedení eroze a vpravo výsledný obraz. Byl použit symetrický strukturní element 3 × 3. Morfologický operátor druhého stupně uzavření je definován jako uspořádaná dvojice operátorů dilatace + eroze, která opět pracuje se stejným strukturním elementem. Při aplikaci tohoto operátoru se nejprve provede operace dilatace, tj. zaplní se dutiny a další bílá místa v objektech, a pak se pomocí eroze uměle zvětšené objekty vrátí přibližně do původní velikosti (obr. 3.7). Název operátoru pochází z jeho základní funkce – uzavírat dutiny v objektech. Volbu velikosti a tvaru strukturního elementu provádíme s ohledem na upravovaný obraz, tj. element musí být větší než největší porucha, kterou má v obrazu odstranit. Musí však být menší než nejmenší objekt nebo jeho část, která má být zachována. Morfologické operátory otevření a uzavření nejsou vůči sobě zcela inverzní, i když by se to podle jejich definice zdálo (pokud by tomu tak bylo, neměly by žádný praktický význam). Některé změny v obrazu jsou totiž nevratné a druhý operátor je již nedokáže napravit. To se při operaci otevření týká nejen poruch na pozadí, které se po odstranění již znovu nemohou vytvořit, ale i některých povrchových útvarů na objektech – srovnejte např. původní a rekonstruované obrazy na obr. 3.6 (kde ke změnám došlo) a na obr. 3.7 (kde změny nenastaly). Největší síla morfologických operátorů spočívá v možnosti zavést podmíněné morfologické operátory. Tyto operátory si můžeme definovat tak, aby přesně řešily náš konkrétní problém. Vezměme si například úlohu odstranit z fotografie všechny objekty, které
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
42
Obrázek 3.7: Morfologický operátor uzavření. Vlevo původní obraz, uprostřed obraz po provedení eroze a vpravo výsledný obraz. Byl použit symetrický strukturní element 3 × 3. jsou přerušeny okrajem fotografie. Tato úloha je typická při analýze velikostí objektů, neboť neúplné objekty by nám vytvářenou statistiku zkreslovaly. Úlohu sice můžeme vyřešit i jinak, napsání vhodného programu by však bylo náročné. Za pomocí podmíněného morfologického operátoru se úloha vyřeší velmi jednoduše. Stačí, když si dodefinujeme operátor eroze, který se aplikuje pouze na objekty dotýkající se okraje fotografie a tato eroze se opakuje tak dlouho, dokud se celý objekt neodstraní. Dodatečně byly morfologické operátory dilatace a eroze (a tím samozřejmě i všechny odvozené operátory) předefinovány tak, aby mohly pracovat i s nebinarizovaným obrazem obsahujícím libovolný počet úrovní šedi. Operátor dilatace byl zaveden tak, že se strukturní element postupně přikládá na všechny dostupné prvky matice A, najde se maximum ze všech hodnot pod strukturním elementem a tímto maximem se nahradí hodnota prvku matice A na nějž ukazuje označený prvek strukturního elementu. Operátor eroze pracuje obdobně, pouze místo maxima se počítá minimum ze všech hodnot pixelů pod strukturním elementem. Tyto nové definice přecházejí pro binarizovaný obraz v definice původní a tím je zajištěna kontinuita metod.
Filtrace pomocí integrálních transformací Třetí a v určitých ohledech nejsilnější skupina technik filtrace obrazů je založena na integrálních transformacích, především na použití Fourierovy transformace. Základní myšlenka spočívá v tom, že se provede přímá Fourierova transformace studované soustavy, ve vzniklém Fourierově obraze se odstraní oblasti odpovídající šumu a po provedení zpětné Fourierovy transformace dostaneme soustavu podstatně méně šumem zatíženou.
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
43
V případě jednorozměrných funkcí bývá typický šum vysokofrekvenční. Převod tohoto původního zašuměného signálu do jeho fourierovského obrazu znamená rozklad ve frekvencích. Pokud předpokládáme, že užitečný signál je převážně nízkofrekvenční, stačí (vhodným způsobem – viz kapitolu 4) potlačit vysokofrekvenční část fourierovského obrazu a tak šum výrazně zeslabit. Při práci s obrazovou informací místo nízkých a vysokých frekvencí se bude jednat o velké a malé objekty. Pokud v dvourozměrném fourierovském obrazu původní obrazové informace opět potlačíme vysokofrekvenční složku, potlačíme ve výsledném obraze malé objekty na úkor větších. Budeme-li předpokládat, že velké objekty jsou skutečné a malé objekty představují poruchy, použití tohoto postupu provede požadovanou filtraci. Integrální transformace nám poskytují ještě mnohem více možností. Pokud známe konkrétní charakteristiku šumu nebo obecněji chyb v obraze které chceme odstranit, přímá část integrální transformace nám umožní jednotlivé složky obrazu dobře odseparovat. Pak můžeme najít chybné prvky, odstranit resp. opravit je a po provedení zpětné transformace dostaneme mnohem kvalitnější obraz. Tyto postupy již překračují rámec pouhé filtrace a patří do speciální partie počítačové fyziky zvané Fourierovská optika.
3.2.4
Binarizace
Dalším úkolem, který v rámci analýzy obrazu nízké úrovně (někdy se používá i pojem předzpracování) musíme provést, je převod obrazu s více úrovněmi šedi (nebo obecně obraz barevný) na obraz obsahující pouze bílou a černou barvu. Tomuto postupu se říká binarizace nebo též prahování. Pro nás to znamená, že matici A, jejíž prvky aij nabývaly hodnot z rozmezí h0, 255i nahradíme maticí s prvky obsahující pouze dvě hodnoty: 0 pro pozadí a 255 (nebo 1) pro objekty. Zda budeme pro objekty používat označení 255 nebo 1 záleží na tom, zda budeme i nadále s maticí A pracovat jako s typem INTEGER(1) a nebo zda ji budeme chtít převést do binární podoby. Tím sice zmenšíme velikost matice na jednu osminu, ale zpomalíme další výpočty, neboť před použitím budeme muset bitové hodnoty jednotlivých pixelů zpět převádět na hodnoty bytové. Standardní postup spočívá v tom, že z původní matice A si vytvoříme histogram, který na ose x obsahuje čísla v rozmezí 0 až 255 označující stupně šedi (pokud pracujeme s 256 úrovněmi šedi v obrazu, jinak si dělení osy patřičně upravíme) a na ose y odpovídající zastoupení jednotlivých odstínů v obrazu. Kvalitní kontrastní obraz vytvoří histogram tvořený dvěma jasně oddělenými maximy. Pak stačí všechny stupně šedi odpovídající levému maximu nahradit hodnotou 0 a stupně odpovídající pravému maximu hodnotou 255, resp. 1, a úloha je vyřešena. Ve většině případů však tak kvalitní obrazy nemáme a dostaneme proto histogram obsahující sice dvě maxima, ale minimum mezi nimi nenabývá nulové hodnoty. Pak je třeba v oblasti tohoto minima vhodně zvolit hodnotu, tzv. práh, která určí hranici od které se budou pixely obrazu počítat k pozadí nebo k objektům. Na obr. 3.8 jsou schématicky uvedeny nejběžnější případy, s kterými se při provádění převodu obrazu do binární podoby můžeme setkat. Je však řeba si uvědomit, že uvedené
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
44
Obrázek 3.8: Použití histogramu četnosti výskytu šedi při binarizaci obrazu. grafy jsou idealizované a že v praxi budou jednotlivé křivky dosti zašuměné, což analýzu dále zkomplikuje. Některé stupně šedi v obrazu mohou zcela chybět, takže histogram stupňů šedi může obsahovat několik minim s hloubkou až k nule, což může zmást některé algoritmy na hledání prahu. Význam jednotlivých křivek v obr. 3.8 je následující: – Případ a: Ideální situace, kdy všechny objekty mají barvu výrazně tmavší než pozadí, takže obě maxima jsou jasně oddělena. Hodnotu prahu volíme kdekoliv v oblasti mezi oběma maximy. Proces binarizace nebude zatížen žádnou chybou. – Případ b: Reálná situace, kdy objekty jsou na obraze zřetelně vidět a rozdíl jasu objektů a pozadí je dosti značný. V obrazu však již existuje přechodová oblast mezi objekty a pozadím, takže proces binarizace nebude zcela přesný. Minimum na histogramu četnosti stupňů šedi je však dosti hluboké, chyba binarizace proto nebude příliš velká.
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
45
– Případ c: Situace, kde analýza bude zatížena značnou chybou. V histogramu sice ještě pozorujeme dvě maxima, ale přechod mezi nimi je nevýrazný, což může v obrazu odpovídat široké přechodové vrstvě mezi objektem a pozadím. – Případ d: Existují dva důvody, proč histogram stupňů šedi vykazuje podobné chování. Velmi často jsou přechody mezi objekty a pozadím tak nevýrazné, že nemá smysl další analýzu fotografie provádět. Existuje však i možnost, že objekty jsou proti pozadí dosti výrazné a přechodová oblast úzká, ale celá fotografie obsahuje prostorový gradient zčernání. Důsledkem může být, že v jedné části obrazu jsou objekty světlejší než v jiné části obrazu pozadí a přitom lokální kontrast objektů vůči pozadí je zachován. V tomto případě lze binarizaci s úspěchem provést, je však nutno nejdříve obraz upravit. Starší postupy předpokládaly rozdělit obraz na několik částí a v každé provést binarizaci samostatně s jiným prahem. Novější postupy využívají možnosti vyrovnání úrovně jasu pozadí. Těmito postupy se často podaří převést situaci na dobře analyzovatelnou (odpovídající případu b). Též situace uvedena v případu c může mít podobnou příčinu a je vhodné se pokusit o stejnou metodu nápravy. Jelikož v řadě případů bývá rozhodování o poloze prahu dosti obtížné, bývá operace binarizace často zdrojem velkých chyb pro celou následnou obrazovou analýzu. Jsou proto navrhována různá sofistikovaná kritéria, jak správně práh zvolit. I když se v některých programech používají různě složité algoritmy pro optimální navržení hodnoty prahu, i zde je vhodné preferovat interaktivní způsob práce a definitivní rozhodnutí ponechat na lidském zraku. Konkrétně by program řešící tuto etapu zpracování obrazu mohl pracovat tak, že v jedné části obrazovky znázorní původní obraz a v druhé vytvořený histogram četností stupňů šedi. Automaticky nebo volbou uživatele se nastaví určitá hodnota prahu a do původního obrazu se dokreslí hranice objektů. Pokud není uživatel spokojen, hodnotu prahu mění tak dlouho, až hranice objektů odpovídají jeho představám. V některých případech může mít uživatel k disposici pomocnou informaci, která s analýzou obrazu přímo nesouvisí. Pokud ji ale na obrazovou převede, může mu to pomoci volit velikost prahu správným způsobem. Například ve fyzice tenkých vrstev může být touto dodatečnou informací elektrický odpor studované struktury nebo množství látky, které studovanou strukturu vytváří.
3.2.5
Rozpoznávání objektů
Poslední etapou analýzy obrazu nízké úrovně je rozpoznávání objektů. Jedná se o problém podobný mnohem známějšímu problému rozpoznáváni písma. Prozatím pro nás obraz reprezentovala matice A s prvky odpovídajícími pixelům obrazu černé nebo bílé barvy (zde předpokládáme, že již proběhla etapa binarizace). Cílem následujícího postupu bude získat informaci o počtu objektů, z kterých se obraz skládá a ke každému objektu přiřadit seznam pixelu tvořících tento objekt. Taková informace již bude postačující pro následnou analýzu vysoké úrovně.
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
46
Při prozkoumávání matice A za účelem identifikace jednotlivých objektů budeme předpokládat, že objekt je tvořený malými čtvercovými prvky, které se dotýkají bočními stranami, tj. kontakt pouze přes úhlopříčku (v jednom bodě) nebudeme za vazbu považovat. Nejjednodušší algoritmus je založen na postupném procházení matice A(i, j) po řádcích a po sloupcích odleva doprava a od shora dolů. – Na počátku si zavedeme jednu buňku paměti jako počitadlo počtu objektů a tuto buňku vynulujeme. – Pokud při procházení matice narazíme na prvek objektu aij mající nenulovou hodnotu, podíváme se na jeho sousedy jimiž jsme již prošli, tj. na dva prvky ai−1,j a ai,j−1 . – Pokud je hodnota obou těchto prvků nulová, víme (resp. se to v tomto okamžiku domníváme), že prvek aij představuje prvek nového objektu. Přidáme proto do počitadla jedničku a vzniklou hodnotu přiřadíme prvku aij . Současně zapíšeme souřadnice (i, j) tohoto prvku mezi pixely nového objektu. – Pokud některý z prvků ai−1,j nebo ai,j−1 obsahuje nenulovou hodnotu, interpretujeme ji jako pořadové číslo objektu a náš prvek aij přiřadíme ke stejnému objektu, tj. jeho číslo uložíme do prvku aij a současně souřadnice tohoto prvku uložíme mezi souřadnice pixelů patřičného objektu. – Postup opakujeme s dalším prvkem matice A. Při ukončení výpočtu máme v počítadle uložen celkový počet objektů a současně známe souřadnice pixelů všech objektů.
Obrázek 3.9: Základní algoritmus rozpoznávání objektů. Pixely, které dávájí nesprávnou informaci o objektu do kterého patří, jsou označeny.
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
47
Ve výjimečných případech může mít výpočet takto hladký průběh. Pokud je však tvar objektů poněkud složitější, dojde občas ke konfliktu, který by výše popsaný algoritmus nedokázal vyřešit. Konflikt bude spočívat v tom, že náš prvek aij bude mít dva nenulové sousedy, ale s rozdílnými hodnotami. To znamená, že jsme založili omylem dva objekty, které jsou ve skutečnosti částmi objektu jediného. Na obr. 3.9 je to ukázáno pro některé pixely, pro které po dotazu nahoru a vlevo byly nesprávně založeny nové objekty. Správný postup v tomto případě je, přiřadit do prvku aij číslo nižšího objektu a vyšší objekt zrušit. Úplné zrušení omylem založeného objektu je však dosti namáhavé, je nutno se postupně v matici A vracet zpět, všechny prvky se špatným číslem objektu přečíslovat, přesunout souřadnice pixelů ke správnému objektu a upravit též počitadlo počtu objektů. Tyto operace se promítnou do složitosti algoritmu. Pokud budeme pracovat se čtvercovou maticí A o hraně N , bude složitost algoritmu O(N 4 ). Poznámka – ve skutečnosti nám většina chybně označených pixelů v obr. 3.9 příliš nevadí, neboť na chybu se přijde ihned v následujícím kroku a oprava je proto okamžitá. Výjimkou je pouze jediný bod ležící v prvním řádku pravé části velkého objektu, který svou chybu přenese na dalších několik desítek bodů ve čtyřech řádcích a oprava proto bude pomalá. Z důvodu silné funkční závislosti složitosti navrženého algoritmu na velikosti matice N je tento algoritmus použitelný pouze pro velmi malé matice. Dalším problémem uvedeného algoritmu je, že z důvodu vracení se a opravování hodnot prvků matice musíme mít neustále celou matici přístupnou, přičemž existují problémy, při kterých matici postupně generujeme a po zpracování příslušného řádku opět zapomínáme. Nutnost ukládat celou matici do paměti počítače limituje maximální velikost takto generované matice. Pro práci s většími maticemi, které navíc nemusejí být po zpracování dále k disposici, byly vytvořeny silnější algoritmy se složitostí až O(N 2 · ln N ) [10]. Základní myšlenka těchto silných algoritmů spočívá v tom, že se při kofliktu nevracíme zpět a objekty a jejich pixely nepřečíslováváme, pouze si nějakým způsobem poznamenáme, které objekty ve skutečnosti tvoří jeden celek. Pak pokračujeme v procházení a analýze matice A, takže doba výpočtu bude úměrná ploše matice N 2 . Teprve po proběhnutí této fáze spočítáme skutečný počet objektů v obrazu, tj. vyloučíme nadbytečné objekty tvořící jeden celek, a složíme dohromady pixely těchto objektů tvořících celek. Režie uvedených operací je řádu ln N . Jednotlivé algoritmy se liší v metodě, jak si informaci o neúplných objektech ukládají a jak je pak efektivně skládají dohromady.
3.3
Metody zpracování obrazu vysoké úrovně
Po ukončení analýzy nízké úrovně budeme mít k dispozici předzpracovaný obraz, tj. matici A obsahující obraz s redukovaným šumem, počitadlo s celkovým počtem objektů a ke každému objektu soubor souřadnic všech jeho pixelů. Obraz v matici A bude převedený do černobílé podoby, pouze ve výjimečném případě se odstíny šedi nebo barva zachovají.
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
48
Při následné analýze vysoké úrovně se budeme snažit získat z předzpracovaného obrazu informaci tam uloženou. Cílem těchto postupů bude najít vhodná kritéria, která nám dovolí od sebe odlišit i velmi podobné obrazy. Tato kritéria jsou obvykle formulována tak, že obraz je popsán nějakou formalizovanou funkční závislostí. Těmto závislostem říkáme morfologické charakteristiky. Většina z nich je založena na teorii matematické morfologie [9], [15], první charakteristiky však byly formulovány již v 30. letech minulého století. I když systém morfologických charakteristik je v současné době již dobře propracován a s úspěchem se používá při řešení nejrůznějších fyzikálních problémů, vývoj v této oblasti stále pokračuje. Práce s morfologickými charakteristikami ve tvaru více či méně složitých funkčních závislostí je v některých případech dosti komplikovaná a není vždy jednoduché nalézt přímou vazbu mezi touto charakteristikou a fyzikálním parametrem studovaného systému (teplotou, tlakem, elektrickým polem, apod.) který příslušné chování analyzovaného obrazu vyvolává. Proto se nové práce v tomto oboru soustřeďují na nalezení tzv. příznaků, což jsou speciální morfologické charakteristiky ve tvaru číselných hodnot místo funkcí. Některé příznaky jsou získávány z obrazu přímo, jiné jsou druhotně odvozovány z morfologických charakteristik.
Pokud chceme informaci obsaženou ve studovaných obrazech nějak kvantitativně charakterizovat, snažíme se popsat určité typické vlastnosti vlastnosti těchto obrazů. Tyto vlastnosti můžeme rozdělit do několika skupin: 1. Integrální informace 2. Informace o jednotlivých objektech 3. Informace o rozložení objektů v rovině 4. Informace o rozložení objektů v prostoru 5. Informace získané z obrazů zachovávajících stupně šedi V následujících paragrafech se seznámíme podrobněji s prvními třemi skupinami charakteristik. Čtvrtá skupina morfologických charakteristik se snaží z dvourozměrného obrazu aspoň částečně rekonstruovat původní třírozměrnou informaci (něco podobného se již delší dobu prakticky používá jako počítačová tomografie). Pátá skupina charakteristik se zabývá informací, která se ztrácí po provedení operace binarizace. Tyto poslední dvě skupiny metod popisu obrazu jsou zatím ve vývoji a jejich výklad překračujeme rámec tohoto základního kurzu přednášek.
3.3.1
Integrální charakteristiky
Integrální charakteristiky popisující obraz jako celek, proto přinášejí nejméně detailní informaci, pro některé účely však plně postačující. Dvě základní morfologické charakteristiky z této skupiny metod jsou:
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
49
– Počet objektů (Koncentrace) Jedná se přímo o dílčí výsledek operace Rozpoznávání objektů – je to obsah počitadla počtu objektů po provedení této operace. Pokud však známe měřítko obrazu, je vhodné tuto charakteristiku normalizovat a prostý počet objektů nahradit objektivnější veličinou koncentrace objektů. – Stupeň pokrytí Tato charakteristika představuje relativní část plochy obrazu zaplněnou objekty. Charakteristika označovaná θ nabývá hodnot z intervalu h0, 1i. Určí se velmi jednoduše z matice A jako relativní počet pixelů (prvků matice) obsazených objekty. Samotná jediná integrální charakteristika nemá příliš velkou vypovídací hodnotu. Existují však oblasti fyziky, kde získáváme sekvence obrazů odpovídající časovému vývoji studované soustavy a pak např. časový vývoj koncentrace objektů může přinést cenné fyzikální poznatky – viz např. výsledky studia růstu tenkých vrstev znázorněné na obr. 3.10
Obrázek 3.10: Závislost koncentrace kovových ostrůvků vytvářených na dielektrické podložce N na parametrech přípravy vrstvy – tloušťce vrstvy s (vyjádřené v monovrstvách) a napařovací rychlosti R. Převzato z práce [11]. Problém studia růstu tenkých vrstev odpovídá situaci znázorněné na obr. 2.16 a 2.17, ukazující fáze nukleace, růstu jednotlivých objektů a jejich slévání. Kombinace těchto tří fází vytváří pozorované změny koncentrace objektů. Srovnáním těchto obrázků a současného obr. 3.10 vidíme, jak různé postupy počítačové fyziky nám ze stejného studovaného jevu mohou přinést rozdílnou informaci, přesně podle potřeb uživatele.
3.3.2
Informace o jednotlivých objektech
Další skupina morfologických charakteristik popisuje vlastnosti jednotlivých objektů. Přitom je důležité si uvědomit základní předpoklad, že při uváděných technikách popisu obrazu
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
50
žádný objekt není důležitější než ostatní. Kvůli objektivitě je proto nutné danou vlastnost studovat u všech objektů a výsledek statisticky vyhodnotit. Zde se konkrétně používá technika histogramů neboli rozdělení četností výskytu dané vlastnosti v souboru objektů. Při dále popsané analýze vyjmeme ze souboru objektů ty objekty, kde daná vlastnost je narušena a statistiku by nám zkreslovaly. Obvykle se jedná o objekty ležící na okraji obrazu tak, že část objektu již není pro analýzu k dispozici. Ve skupině charakteristik popisujících jednotlivé objekty máme na výběr více metod a v konkrétním případě volíme jednu nebo několik z nich tak, aby co nejlépe charakterizovaly objekty ve studovaném obrazu a aby byly co nejvíce citlivé i na malé změny ve studovaném obrazu.
Obrázek 3.11: Závislost rozdělení poloměrů kruhových objektů na parametrech přípravy vrstvy – tloušťce vrstvy s (v monovrstvách) a napařovací rychlosti R [11]. – Rozdělení poloměrů Pokud mají objekty tvar kruhů (resp. u třírozměrných objektů jejich průměty mají tvar kruhu), máme v této skupině k dispozici jedinou morfologickou charakteristiku. Spočteme poloměr každého objektu a sestavíme z nich rozdělovací funkci – příklad je znázorněn na obr. 3.11. K nalezení poloměru objektu slouží celkový počet pixelů tvořících daný objekt (s využitím informace o kroku digitalizace, tj. rastru matice A).
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
51
– Rozdělení efektivních poloměrů Pokud mají objekty nepravidelný tvar, je třeba k jejich popisu použít dvě charakteristiky současně – první bude rozdělení efektivních poloměrů a druhou popis odchylek tvarů objektů od přesně kruhového tvaru, zde použijeme popis pomocí tvarového faktoru. Při určení rozdělení efektivních poloměrů postupujeme tak jako by objekty byly přesně kruhové, tj. spočítáme počet pixelů objektu (normalizovaný podle rastru matice) a po2 ložíme jej rovný πRef . – Rozdělení tvarových faktorů Tvarový faktor má za úkol statisticky charakterizovat odchylky tvaru daného objektu od kruhu. Zavádíme jej obvykle tak, aby pro dokonalý kruh nabýval tvarový faktor hodnoty jedna a s růstem deformace objektu aby jeho hodnota klesala k nule. Tomuto popisu vyhovuje definice tvarového faktoru F F (značka pochází z anglického „form factorÿ) F F = 4π ·
S , O2
(3.1)
kde S je plocha objektu a O jeho obvod. Hodnoty S a O se musejí určit z počtu pixelů příslušejících k danému objektu (toto určení, zejména u malých objektů, nemusí být vždy jednoduché). – Rozdělení velkých a malých poloos a rozdělení směrů (efektivních) elips Pokud není vhodné objekt charakterizovat kruhovým tvarem, ať již přesně nebo přibližně, hledá se jiný tvar, jemuž by objekty analyzovaného obrazu vyhovovaly. Relativně často bývá tímto tvarem elipsa. Postup práce je obdobný jako práce s kruhem, pouze místo jednoho parametru-poloměru máme parametry tři – délku velké a malé poloosy a orientaci elipsy v prostoru. V případě přesně eliptických objektů proto konstruujeme tři histogramy, v případě přibližných tvarů navíc ještě zavádíme příslušně modifikované tvarové faktory (3.1). – ... Pro složitější tvary objektů postupujeme podobně. Vždy se snažíme nalézt nějaký idealizovaný tvar, uděláme histogram jeho typického parametru a případně i odchylek od tohoto ideálního tvaru. Příkladem může být popis malých krystalických objektů, které mají v průmětu obvykle tvar pravidelných mnohoúhelníků. – ... Pro ještě složitější tvary objektů se pracuje s kritérii také složitějšími. Kvantitativně se popisuje konkávnost nebo konvexnost objektů, počet vrcholů, rozdělení vrcholových úhlů, apod. Každý uživatel má možnost si pro svůj konkrétní případ dodefinovat kritérium podle své potřeby s využitím výše popsaných zásad. Na obr. 3.11 jsou znázorněna rozdělení poloměrů pro dva odlišné obrazy. Jedná se opět o případ analýzy růstu tenkých vrstev kovů na dielektrických podložkách převzatý podobně jako předchozí obrázek z práce [11].
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
3.3.3
52
Informace o rozložení objektů
Třetí skupinou morfologických charakteristik je popis rozložení objektů v obrazu. Toto rozložení může být skutečné (pokud obraz odpovídá dvourozměrnému rozložení objektů v přírodě) nebo jen zdánlivé (pokud dvourozměrný obraz vznikl jako průmět prostorového rozložení objektů – např. astronomická fotografie nebo mikrofotografie tenké kompozitní vrstvy). V některých oblastech vědy se jedná o nejdůležitější skupinu morfologických charakteristik, především v astronomii, metalurgii, fyzice tenkých vrstev, apod. Proto tato skupina morfologických charakteristik je nejrozsáhlejší a též první morfologické charakteristiky vůbec byly navrženy astronomy pro charakterizování rozložení hvězd a dalších vesmírných objektů. Naopak existují oblasti vědy, kde se tato skupina morfologických metod vůbec neuplatní – a to tehdy, kdy objekty jsou pohyblivé a obraz představuje pouze jejich náhodné uskupení. Jako příklad si můžeme uvést oblast lékařskou, kde se často počítačově zpracovávají fotografie krve získané v optickém mikroskopu a pomocí koncentrace a tvarů krvinek a dalších objektů se provádí diagnóza. Proto jsou zde základní morfologickou skupinou metody popisujících tvary objektů (a integrální charakteristiky), zatímco rozložení objektů v obrazu je zcela nezajímavé. Metody určené pro popis rozložení objektů po podložce mají ve skutečnosti kořeny dva. Jedním je již zmíněná astronomie, kde je obraz tvořen převážně bodovými objekty. Druhá skupina metod je spojená s materiálovým výzkumem jenž byl podnětem pro vytvoření metody matematické morfologie [9], kde naopak jsou typické objekty plošné. Proto též morfologické charakteristiky přinášející informaci o rozložení objektů v ploše můžeme rozdělit na dvě třídy: – metody popisující rozložení bodových objektů, a – metody popisující rozložení plošných objektů. Některé algoritmy jsou použitelné pouze pro jeden typ objektů, ať již bodových nebo plošných, jiné metody jsou univerzální a můžeme je použít přímo nebo po určité modifikaci pro oba typy objektů.
Uveďme si nejpodstatnější morfologické charakteristiky vhodné pro popis rozložení objektů v obraze: – Radiální distribuční funkce – Rozdělení nejbližších sousedů – Kovariance Nyní si tyto metody probereme detailně. Další metody, které překračují rozsah tohoto studijního textu, budou součástí výuky ve vyšších ročnících.
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
53
– Radiální distribuční funkce Jedná se o metodu navrženou původně pro bodové objekty. Tuto metodu můžeme zobecnit i na plošné objekty za předpokladu, že jsou mnohem menší než vzdálenost mezi nimi – v tomto případě objekty můžeme přibližně nahradit jejich těžišti (která opět vypočítáme ze seznamu pixelů tvořících daný objekt).
Obrázek 3.12: K zavedení radiální distribuční funkce. Radiální distribuční funkci (někdy se též nazývá párová distribuční funkce) P (r) určíme způsobem znázorněným na obr. 3.12: – Vezmeme nějaký objekt jako počátek soustavy souřadné a vedeme z něho osu r. – Na této ose zvolíme interval šířky ∆r, hr, r + ∆ri. – Vytvoříme mezikruží a spočítáme počet objektů v něm ležících, ∆n. – Postup opakujeme pro další r měnící se od nuly až po nějaké rmax . – Radiální distribuční funkce se spočítá podle vztahu ∆n . 1 P (r) = · , n0 2πr∆r
(3.2)
kde n0 je koncentrace objektů. Tento vztah je dostatečně přesný pro velká r, poblíž objektu musíme plochu mezikruží určit přesněji. Při použití výše uvedeného postupu by byla získaná radiální distribuční funkce P (r) příliš zašuměná a šum by se zvětšoval poblíž objektu, neboť tam plocha mezikruží klesá. Řešením je postup opakovat pro větší počet objektů ze střední části obrazu a následným průměrováním získaných funkcí šum redukovat.
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
54
Radiální distribuční funkce má význam relativní lokální koncentrace ve vzdálenosti r od daného objektu. Pro zcela náhodně rozložené bodové objekty bychom měli dostat konstantní hodnotu rovnou jedné. Nerovnoměrnost v rozložení objektů se projeví ve zmenšení (oblast, kde je objektů méně) nebo zvětšení (pokud je lokálně objektů více) této konstantní hodnoty. Podle této úvahy bychom mohli pozorovat následující typy funkcí P (r) (obr. 3.13 – převzato z práce [12]):
Obrázek 3.13: Základní typy radiálních distribučních funkcí RDF .
– Zašuměná konstanta (rovná jedné) Tato situace odpovídá náhodnému rozložení bodových objektů. – Skokový vzrůst funkce P (r) od nuly k jedné Tato situace odpovídá náhodnému rozložení plošných objektů. Interval na ose r, kde je funkce P (r) nulová, je roven dvojnásobku poloměru objektů.
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
55
– Pozvolný vzrůst funkce P (r) od nuly k jedné Tento průběh značí nějaký fyzikální proces, který brání těsnému přiblížení objektů. – Vzrůst funkce P (r) s jedním nebo více zákmity Zatímco všechny předchozí průběhy označovaly náhodné rozložení objektů v ploše, zákmity na funkci P (r) napovídají, že objekty nejsou rozložené náhodně. Intenzita zákmitu nebo zákmitů určuje stupeň uspořádanosti studované soustavy, poloha maxim a minim nám dává informaci o prostorovém uspořádání soustavy. – Rozdělení nejbližších sousedů Algoritmus pro výpočet rozdělení nejbližších sousedů je podobný výpočtu radiální distribuční funkce. V obou případech musíme počítat vzdálenosti mezi různými objekty. V případě výpočtu radiální distribuční funkce z těchto vzdáleností sumováním podle určitého klíče dostaneme rozdělení P (r). Rozdělení nejbližších sousedů se spočítá tak, že pro každý objekt najdeme minimum ze vzdáleností k ostatním objektům a z těchto minim se vytvoří histogram. Příklady rozdělení nejbližších sousedů (anglicky „Distribution of Nearest Neighboursÿ, odtud zkratka DN N ) jsou uvedeny spolu s radiálními distribučními funkcemi RDF na obr. 3.13. S nepatrně zvýšenou námahou lze současně vytvořit i rozdělení 2. nejbližších sousedů, pak 3., atd. Soubor těchto histogramů má téměř stejnou informační hodnotu jako radiální distribuční funkce. Pro plošné objekty existují dvě modifikace metody rozdělení nejbližších sousedů. Pro ně (na rozdíl od bodových objektů) lze zvolit, zda vzdálenosti budeme určovat mezi středy (těžišti) objektů nebo mezi jejich okraji. – Kovariance Kovariance je metoda plošná, kterou nelze pro bodové objekty použít. Metoda má více formulací, my si zde uvedeme formulaci maticovou. Kovarianční funkce C(h) se spočítá ze vztahu 1 XX C(h) = · A(i, j) · A(i + h, j) . (3.3) K i j Zde A(i, j) je matice obsahující digitalizovaný obraz a A(i + h, j) je tatáž matice posunutá ve směru osy x o nějaký úsek h. Toto posunutí můžeme provádět v libovolném směru, z praktických důvodů se však volí pouze kladné nebo záporné posuvy podél obou os. Sčítání ve vztahu (3.3) se provádí přes oba indexy matice a to pouze v rozsahu, kde se matice překrývají. Konstanta K je určena k normalizaci výsledku a její velikost je rovna velikosti překryvu obou matic. V kovarianční funkci C(h) (tuto funkci dostaneme, když vztah (3.3) aplikujeme pro řadu hodnot posunutí h) je obsaženo větší množství morfologické informace: – hodnota kovariance v počátku, C(0), udává stupeň pokrytí θ, – stupeň pokrytí též můžeme získat z hodnot C(h) pro velká posunutí jako
q
C(∞),
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
56
– směrnice kovarianční funkce poblíž počátku, −(dC(h)/dh)h=0 , vypovídá o středním poloměru objektů, – pokud je průběh kovarianční funkce monotonní, vypovídá to o rovnoměrném rozložení objektů v obraze, a – pokud jsou na kovarianční funkci vidět maxima a minima, nasvědčuje to nerovnoměrnému rozložení objektů a z poloh a velikostí těchto extrémů můžeme usuzovat na charakter nerovnoměrností. Velkou výhodou metody kovariance je, že posuvy lze provádět v různých směrech a tak můžeme popsat směrovou závislost v obrazu, pokud tam existuje. Předchozí metody, kde se přes různé směry sčítalo, tuto možnost neměly. Na obr. 3.14 je schématicky znázorněn soubor osmi obrazů (odpovídající kovovým vrstvám s postupně rostoucí tloušťkou) a příslušné kovariační funkce C(h).
Obrázek 3.14: Určení kovariance souboru objektů. Převzato z práce [13].
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
57
Na obr. 3.15 je uveden příklad komplexnější analýzy souboru objektů – obrázek obsahuje radiální distribuční funkci RDF , rozdělení nejbližších sousedů DN N a kovarianci C(h).
Obrázek 3.15: Kombinovaná analýza obrazu s použitím více morfologických metod – radiální distribuční funkce, rozdělení nejbližších sousedů a kovariance. Převzato z práce [14]. Kromě výše uvedených morfologických metod pro popis rozložení objektů v obraze existuje ještě velké množství metod dalších. Za významnější lze považovat: – kontaktní distribuční funkci – popis pomocí Wignerových-Seitzových buněk (a příbuzného Voronoiovského dláždění) – metodu „Chord-length distributionÿ, a především – metodu „Quadrat Countsÿ. Tato metoda je důležitá tím, že výstupem z ní je jednočíselná hodnota, takže tato morfologická metoda je současně i příznakem. Rozboru všech morfologických metod, těchto i uvedených výše, a hledání příznaků v nich bude věnována samostatná přednáška ve vyšším ročníku studia. Některé zkušenosti s využitím těchto morfologických metod při jejich aplikaci na fyziku tenkých vrstev jsou popsány v [15], [16].
3.4
Další problémy z oblasti zpracování obrazu
Další oblastí, kde se v současnosti ještě intezívně pracuje a stále se hledají nové morfologické metody a jejich příznaky, je snaha o získání informací o rozložení objektů v prostoru z jejich
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
58
pouze dvourozměrných obrazů. K tomuto problému existují dva přístupy. První, který representuje především počítačová tomografie, je založen na tom, že informaci ztracenou převodem třírozměrného souboru objektů do jeho plošného obrazu doplníme tak, že těchto obrazů získáme větší počet a to z různých směrů. Pak je úplná rekonstrukce ztracené informace aspoň teoreticky možná, i když při praktické realizaci nastávají určité problémy s výpočetní náročností procesu rekonstrukce. Druhý přístup je založený na práci s jediným průmětem. Tato situace nastává v případě, že z nějakých důvodů nelze soustavu objektů snímat z různých stran – např. v astronomii, kde studované objekty jsou příliš daleko od Země, nebo ve fyzice tenkých vrstev, kde studované objekty jsou příliš malé a celá geometrie problému je dvourozměrná. Zde úplná rekonstrukce možná není, ale za pomoci aparátu matematické statistiky a s případným využitím různé doplňkové informace lze aspoň získat rozdělovací funkce rozložení objektů a jejich velikostní charakteristiky. Pro kompozitní vrstvy je současný stav poznání této problematiky uveřejněn v práci [17].
3.5
Shrnutí
V této kapitole je těžiště druhého dílu studijního textu z Počítačové fyziky. Pozornost je třeba věnovat především: • Přehledu jednotlivých kroků probíhajících při úplné analýze obrazu. • Z metod zpracování obrazu nízké úrovně: – metodám filtrace – metodám rozpoznávání objektů. • Z metod zpracování obrazu vysoké úrovně: – metodách charakterizujícím velikosti a tvary jednotlivých objektů – metodám určeným pro popis rozložení objektů na ploše.
3.6
Problémy k řešení
Problémy určené k procvičování látky popsané v této kapitole: 1. Seznamte se s činností metod filtrace obrazu tak, jak jsou presentovány v demonstrační části knihovny Image Processing Toolbox systému MATLAB. 2. Zkuste navrhnout silnější algoritmus pro rozpoznávání objektů s nižší závislostí na velikosti matice N než je standardní O(N 4 ). Nápověda: pro zlepšení efektivity algoritmu je třeba potlačit vracení se s cílem opravit chybně označené pixely.
KAPITOLA 3. ZPRACOVÁNÍ OBRAZU
59
3. Nagenerujte náhodně 1000 bodů rozmístěných v matici o rozměrech 1000 × 1000 a spočítejte jejich radiální distribuční funkci. 4. Udělejte totéž pro středy náhodně generovaných kružnic o konstantním poloměru, přičemž ohlídejte aby se kružnice neprotínaly (pokud by k tomu došlo, danou kružnici zrušte a nagenerujte místo ní jinou). Počet kružnic volte opět 1000. Výpočet opakujte pro několik různě velkých poloměrů kružnice – pozor, pro příliš velký poloměr se všechny kružnice do matice nevejdou. Srovnávejte změny radiální distribuční funkce s rostoucím poloměrem kružnic. Nápověda: tím, že nedovolujeme překrývání kružnic, zvětšujeme uspořádanost souboru bodů ≡ středů kružnic. 5. Spočítejte pro soubory bodů nagenerovaných v předchozím problému jejich kovarianci. Sledujte jak se mění charakter kovarianční funkce v závislosti na rostoucí uspořádanosti souboru.
Kapitola 4 INTEGRÁLNÍ TRANSFORMACE Integrální transformace jsou jednou z obtížných partií matematiky, kde se tyto transformace vyvíjejí a studují se jejich vlastnosti. Důležité místo ale mají i ve fyzice, kde mají rozsáhlé využití při zpracování výsledků měření. Jsou základem několika vědních disciplín jako je harmonická analýza nebo fourierovská optika. Existuje více typů integrálních transformací a většina z nich je ve fyzice prakticky využitelná. Výsadní místo mezi nimi má Fourierova transformace. Tato integrální transformace je známa a používána již od 19. století, její mimořádna popularita však nastala až po roce 1965, kdy byla objevena tzv. rychlá Fourierova transformace [18]. Základní myšlenkou všech typů integrálních transformací je, že obousměrným způsobem propojují dvě veličiny – obvykle čas a frekvenci. Tak zvanou přímou transformací převedeme časově závislé děje v přírodních vědách na jejich obrazy obsahující frekvenční rozklad studovaného procesu. Inverzní proces, tak zvaná zpětná transformace, opět převede frekvenční závislosti na závislosti časové. V principu existuje dvojí využití těchto procesů. Pokud nás přímo zajímá frekvenční spektrum zkoumaného přirodního jevu, provedeme pouze přímou integrální transformaci. Takto postupoval tvůrce této metody baron Fourier, který použil po něm pojmenovanou transformaci k řešení parciální diferenciální rovnice popisující nestacionární vedení tepla. Jelikož řešení této rovnice nemohl nalézt, provedl její transformaci a na tomto transformovaném obrazu se mu již podařilo popsat vlastnosti studovaného systému a posoudit dynamiku jeho chování. Univerzálnější ale je postup, při kterém provedeme přímou transformaci rovnic popisující náš systém. Převážná většina rovnic se tímto postupem zjednoduší a mohou být vyřešeny. Získaný výsledek se pak zpětnou transformací převede do obvyklých (časových) proměnných, přičemž tento výsledek je i řešením původního problému. Tuto původní myšlenku využití přímé a zpětné integrální transformace při řešení matematických úloh je možno zobecnit. Po převedení našeho problému do jeho obrazu můžeme tento obraz upravit a tyto úpravy se po zpětné transformaci promítnou i do původního systému. V předchozí kapitole jsme se již zmínili, že tanto postup je silnou metodou filtrace zašuměných obrazů. Pokud obraz v němž je zakomponovaný šum přetransformujeme, dokážeme lehce užitečný signál a šumovou složku od sebe odlišit, šum pak potlačíme a po 60
KAPITOLA 4. INTEGRÁLNÍ TRANSFORMACE
61
provedení zpětné transformace získáme podstatně kvalitnější obraz. Zápis integrálních transformací a především jejich praktická aplikace jsou dosti složité a překračují rámec základního kurzu, proto si v této kapitole uvedeme pouze základní vztahy a naznačíme, kterým směrem se vývoj ubírá. Pokud bude třeba přímou a zpětnou transformaci prakticky provádět, doporučuji se opřít o knihovní programy, kterých existuje dostatečný počet a nesnažit se programovat si celý výpočet sám. To platí zejména při práci s rychlou Fourierovou transformací a ještě obtížnější je situace při zpracovávání dvourozměrné obrazové informace. Knihovní programy pro přímou i zpětnou Fourierovu transformaci v jedné i dvou dimenzích jsou součástí řady programových balíků. Mezi nejdostupnější patří zejména tzv. „numerické receptyÿ, tj. soubor knih Numerical Recipes v různých programovacích jazycích – existují jazykové mutace pro FORTRAN 77 a 90, C, C++ a Pascal. Citace knihy Numerical Recipes in Fortran 77 je např. [19]. Detailnější informace o teorii těchto transformací lze získat v [20]. V následujících sekcích budeme pracovat se skalárními funkcemi jedné proměnné. Tyto funkce mohou být reálné i komplexní. Pro vektorové funkce a pro funkce více proměnných lze uvedené vztahy jednoduše zobecnit.
4.1
Fourierovy řady
Mějme periodickou obecně komplexní funkci x(t) s periodou T > 0. Potom Fourierovou řadou k této funkci nazveme součet ξ(t) =
∞ X
Ck e
2πjkt T
,
−∞ < t < ∞ .
(4.1)
k=−∞
Zde i v dalším textu budeme symbolem j označovat imaginární jednotku. Pro koeficienty Ck platí −2πjkt 1 Z t0 +T Ck = x(t) e T dt , k = 0, ±1, ±2, . . . , (4.2) T t0 kde t0 je libovolný bod a o funkci x(t) předpokládáme, že platí x(t) = x(t + ν · T ) ,
ν = 0, ±1, ±2, . . . .
Pokud je funkce x(t) reálná, pro koeficienty Ck platí Im C0 = 0 C−k = Ck∗ ,
k = 1, 2, . . . ,
kde symbolem ∗ označujeme komplexně sdruženou hodnotu. Důležitou vlastností z hlediska praktických aplikací je, že s rostoucím indexem k klesají absolutní hodnoty koeficientů Ck k nule. Tento pokles je dosti rychlý: pro funkci x(t) spojitou včetně svých derivací až do řádu α − 1 klesají koeficienty Ck k nule jako | k |α+1 . I pokud je funkce x(t) nespojitá, klesají koeficienty k nule jako | k |.
KAPITOLA 4. INTEGRÁLNÍ TRANSFORMACE
62
Tato rychlost konvergence se projeví na přesnosti výpočtu. My při použití Fourierových řad nahrazujeme původní funkci x(t) jejím rozvojem, v kterém bereme jenom konečný počet členů a chyba je dána velikostí členů vynechaných.
4.2
Fourierova transformace
Na rozdíl od Fourierovy řady obsahující diskrétní koeficienty Ck , pracuje Fourierova transformace se spojitým spektrem. Zavádíme přímou a zpětnou Fourierovu transformaci. Přímá transformace je definována vztahem X(ω) =
Z ∞ −∞
x(t) e−jωt dt .
(4.3)
Zpětná transformace je definována vztahem 1 Z∞ x(t) = X(ω) e+jωt dω . 2π −∞
(4.4)
Pro funkce X(ω) platí podobná tvrzení jako pro koeficienty Ck Fourierovy řady, především že tyto funkce se zmenšují se zvyšující se frekvencí ω lim X(ω) = 0 .
ω→∞
Pokud je funkce x(t) nenulová pouze na konečném intervalu, lze si ji představit jako periodicky prodlouženou do +∞ i do −∞. Potom lze mluvit o její Fourierově řadě. Pro koeficienty této řady dostaneme vztahy, které propojují Fourierovu transformaci a Fourierovu řadu ! 1 2πk Ck = X . T T Základní vlastnosti Fourierovy transformace, které se projeví při při jejím použití ve fyzice, jsou – linearita – konvoluce – dekonvoluce – autokorelace – derivace. Fyzikální význam Fourierovy transformace a Fourierovy řady je v tom, že pokud původní funkce x(t) představovala časovou závislost nějaké veličiny, vyjadřují Fourierovská transformace nebo řada spektrální rozklad této veličiny do prostoru frekvencí. Absolutní hodnota transformované veličiny určuje, jak silně je příslušná frekvence ve spektru zastoupena, fázi této frekvence lze spočítat z poměru reálné a imaginární části té transformované veličiny.
KAPITOLA 4. INTEGRÁLNÍ TRANSFORMACE
63
Vztahy (4.1) a (4.4) ukazují, že průběh každé fyzikální veličiny můžeme vyjádřit jako superpozici nekonečně mnoha harmonických funkcí. Po provedení přímé transformace můžeme každou složku tohoto rozkladu vyšetřovat samostatně. Po provedení požadovaných změn můžeme pomocí zpětné transformace opět poskládat jednotlivé harmonické funkce dohromady a vytvořit původní nebo modifikovanou fyzikální veličinu. Limitní chování koeficientů Ck a funkcí X(ω) je odrazem faktu, že žádný děj v přírodě nemůže probíhat s nekonečnou rychlostí (tj. frekvencí), neboť by to znamenalo nekonečnou energii potřebnou k její změně.
4.3
Diskrétní Fourierova transformace
Jak je známo již z numerické matematiky, počítače nejsou vhodné pro práci se spojitými veličinami. Pro práci s Fourierovou transformací při diskrétním vyjádření a ekvidistantním dělení nezávisle proměnných t a ω se místo původní transformace (4.3) a (4.4) zavádí diskrétní Fourierova transformace. K tomuto účelu zavedeme následující značení. Vezmeme posloupnost N obecně komplexních čísel xi , i = 0, 1, . . . , N − 1. Obvykle je tato posloupnost získána ekvidistantním navzorkováním původní funkce x(t) na intervalu h0, T i. Transformací této posloupnosti xi získáme opět obecně komplexní posloupnost Xk =
N −1 X
xi e
−2πjik N
,
k = 0, 1, . . . , N − 1 .
(4.5)
i=0
Toto je přímá diskrétní Fourierova transformace. Zpětná diskrétní Fourierova transformace je definována vztahem xi =
N −1 X
Xk e
2πjik N
,
i = 0, 1, . . . , N − 1 .
(4.6)
k=0
Pokud je původní funkce x(t) a tím i posloupnost xi reálná, platí Xk = XN∗ −k ,
k = 1, 2, . . . , N − 1 .
Význam posloupnosti Xk je obdobný významu posloupnosti Ck a funkcí X(ω), neboli představuje frekvenční rozklad původní funkce x(t) do prostoru frekvencí. Tam můžeme provést plánované operace a zpětnou transformací (4.6) dostat modifikovanou funkci x(t). Při praktické aplikaci transformací (4.5) a (4.6) narážíme na řadu technických problémů. Proto je třeba s fourierovskými transformaci postupovat obezřetně, jinak se nám v rekonstruované funkci objeví falešné signály, které nám výsledek zkreslí.
4.4
Rychlá Fourierova transformace
Rychlá Fourierova transformace nepředstavuje další soustavu vzorců, ale jedná se o metodiku efektivního počítání diskrétní Fourierovy transformace (4.5), (4.6). Pokud bychom
KAPITOLA 4. INTEGRÁLNÍ TRANSFORMACE
64
počítali diskrétní transformaci přímo, výpočet by byl pro větší délku vstupní posloupnosti N časově náročný. Pokud bychom navíc prováděli transformaci dvourozměrnou, což je případ zpracování obrazů ve fourierovské optice, doba výpočtu by se prodloužila na hodiny až desítky hodin. Prakticky použitelnou se tedy Fourierova transformace stala tehdy, když se podařilo nalézt efektivní algoritmus jejího výpočtu. Princip rychlé Fourierovy transformace je velmi jednoduchý. Jestliže počet prvků posloupnosti N lze rozložit na součin dvou celých čísel K · L, potom lze i prováděné operace rozložit do dvou skupin o K a L prvcích. K výpočtu jednoho prvku potom stačí K + L operací a celá transformace vyžaduje přibližně N · (K + L) operací, což je méně než při přímém zpracování posloupnosti o N prvcích. Výpočet je efektivní tehdy, pokud volíme délku posloupnosti N jako mocninu čísla 2, N = 2m . Potom, jestliže výše uvedený rozklad čísla N na dvě čísla K a L provedeme, K a L budou nižší mocniny čísla 2. Jestliže celý postup budeme opakovat tak dlouho, až budou všechny dílčí posloupnosti mít délku 2, bude platit, že výpočetní náročnost celé transformace bude jen N · m = N · ln 2 . Pokud to srovnáme s výpočetní náročností původní diskrétní Fourierovy transformace, jež je N 2 , je úspora času zřejmá. Při zpracování obrazu, zejména když budeme pracovat s větším rozlišením, se může jednat o rozdíl několika řádů, tj. sekundy místo desítek hodin. Je tedy zřejmé, že teprve objevení algoritmů (konkrétních modifikací popsané základní myšlenky totiž existuje více) rychlé Fourierovy transformace umožnilo prudký rozvoj aplikací integrálních transformací ve fyzice a dalších vědních disciplínách.
4.5
Další integrální transformace
I když je Fourierova transformace v praxi nejvíce rozšířena, existuje větší počet dalších integrálních transformací. Pro některé z nich byly též nalezeny rychlé algoritmy jejich výpočtu a mohou proto Fourierově transformaci v efektivitě práce úspěšně konkurovat. O jejich nasazení pro řešení konkrétních fyzikálních problémů proto rozhodují další kritéria. Zaprvé existují reálné integrální transformace. Určitá slabina fourierovských transformací je v tom, že i když vstupní funkce x(t) je reálná, její obraz X(ω) je komplexní a totéž platí pro diskrétní posloupnost Xk . Protože se musí v tomto případě manipulovat s reálnou a imaginární části fourierovských obrazů, je výpočet relativně pomalejší a především náročný na ukládání dat. V případě jednorozměrných signálů toto omezení obvykle nepociťujeme, avšak při zpracování obrazové informace se často jedná až o stovky MB. Při použití reálné integrální transformace je i obrazová posloupnost reálná a tato omezení odpadají. Existuje několik variant reálných integrálních transformaci. Společné pro ně je, že místo rozkladu do harmonických funkcí ve tvaru komplexních exponenciel pracují s funkcemi sin, cos a s funkcemi od nich odvozenými, např. sin + cos. Posledně uvedená transformace se nazývá Hartleyova transformace a její výhoda je kromě reálnosti i v tom, že vztah pro přímou i zpětnou transformaci je identický, což ulehčuje přípravu kódu.
KAPITOLA 4. INTEGRÁLNÍ TRANSFORMACE
65
Zvláštní postavení v posledních letech získala tzv. waveletová transformace. Jedná se o mimořádně silný aparát, jenž může být použit pro řešení velkého množství nejrůznějších problémů. Teorie této transformace je příliš složitá i pro zjednodušený popis. Velmi stručně ji lze charakterizovat tak, že místo aby se provedl úplný integrál nebo součet přes celý definiční obor časů (u přímých transformací) nebo frekvencí (u zpětných transformací), provádí se skládání přes obě složky současně, ale jen v úzkých intervalech. Tím se podstatně zvýší vypovídací hodnota transformovaných signálů. Jako analogii můžeme použít na jedné straně převedení nějakého hudebního díla do frekvenčního prostoru, kdy se převážná část informace a celá krása skladby nenávratně ztratí, což je případ Fourierovy a dalších příbuzných transformací. Na druhé straně nám waveletová transformace představuje zápis této skladby v notách. I zde k určité redukci informace došlo, ale to co zůstalo může být použito k novému nahrání této skladby.
4.6
Shrnutí
Integrální transformace jsou velmi silnou partií počítačové fyziky a mají mnoho aplikací – od potlačování šumu v experimentální fyzice, přes zpracování obrazu až po provádění vícenásobných derivací a nebo efektivního řešení parciálních diferenciálních rovnic. Jsou však prakticky jen obtížně aplikovatelné a odpovídají teorie překračuje zaměření tohoto kurzu. Proto jako partie k zapamatování jsou uvedeny pouze základy, a to: • Definiční vztahy Fourierovy řady. • Definiční vztahy přímé a zpětné Fourierovy transformace. • Diskrétní vztahy přímé a zpětné diskrétní Fourierovy transformace. • Základní myšlenka algoritmů rychlé Fourierovy transformace.
4.7
Problémy k řešení
Jelikož se jedná o příliš obtížnou partii, nejsou uváděny žádné problémy k řešení v rámci studia. Místo toho by se měl čtenář seznámit s konkrétními algoritmy pro diskrétní Fourierovu transformaci v knihovně programů, kterou má k dispozici. Dobře je tato partie zpracována např. v Numerical Recipes [19] a nebo v systému MATLAB, kde je možno si prohlédnout i připravené demonstrační programy.
Kapitola 5 DALŠÍ SMĚRY POČÍTAČOVÉ FYZIKY Ve dvou dílech tohoto studijního textu jsme se seznámili s většinou důležitých partií počítačové fyziky, a to v prvním díle s různými technikami počítačového modelování, jež tvoří jádro počítačové fyziky, a v tomto druhém díle s počítačovou graikou, zpracováním obrazu a integrálními transformacemi.
5.1
Řízení experimentu
Do výčtu relativně často používaných partií počítačové fyziky chybí ještě dvě. Za prvé se jedná o řízení experimentu. Tato partie leží na rozhraní počítačové fyziky a elektroniky. Ve skutečnosti se jedná o dva okruhy problémů: – Jednodušší úrovní je sběr a zpracování dat, což představuje problém získat data ze studovaného problému z experimentální aparatury a přenést je do počítače, kde už je bude možno analyzovat postupy počítačové fyziky. – Vyšší úrovní je automatizace nebo též kybernetizace experimentu. Ta se skládá z první části sběru dat, ale v tomto případě jsou data po zpracování v počítači použita k ovládání experimentální aparatury a tak k modifikaci dalšího průběhu sběru dat. Oba okruhy obsahují část patřící plně do elektroniky, resp. přístrojové techniky – problém snímání různých fyzikálních veličin a jejich převod na veličiny elektrické (snímače), transformace elektrických veličin do digitálních signálů v normě vhodné pro přenos do počítače (analogově-digitální převodníky), problém přenosu do počítače (sběrnice), a v případě kybernetizace navíc ještě problém převodu digitálních signálů na jiné elektrické veličiny (digitálně-analogové převodníky) a vlastní ovládání experimentální aparatury (regulátory, krokové motory, atd.). Pak jsou tam části patřící plně do počítačové fyziky – vyhodnocení naměřených dat včetně potlačování šumu, zpracování obrazu, atd. Hraničním okruhem patřícím do obou vědních disciplín je plánování celého procesu řízení experimentu, tj. vypracování schématu řízení. Zde rozeznáváme několik úrovní: 66
KAPITOLA 5. DALŠÍ SMĚRY POČÍTAČOVÉ FYZIKY
67
– systém jednoho uživatele, – systém více uživatelů, – hierarchický systém, – použití speciálních procesorů. V minulosti se současně používaly všechny čtyři úrovně v závislosti na složitosti experimentálního zařízení (určitě je rozdíl mezi experimentální aparaturou v laboratoři školy a velkým vědeckých zařízením typu CERNu nebo termojaderného výzkumu na tokamacích). V současnosti masový nástup relativně výkonných mikropočítačů rozdíl mezi těmito úrovněmi řízení poněkud setřel. Stručně lze rozdíly mezi jednotlivými schématy řízení charakterizovat tím, kde jsou data pocházející z experimentu zpracovávána. V případě systému jednoho uživatele je ke každého experimentu připojen jeden počítač, který vykonává veškerou řídící i výpočetní práci (to je v současnosti nejběžnější případ při využití mikropočítačů typu PC). Systém více uživatelů byl založen na silnější řídící technice, která dokázala obsloužit několik experimentů současně. Výhodou byla možnost centralizovat finanční prostředky a takový počítač vybavit kvalitními perifériemi. Nevýhodou byla možnost vzájemného ovlivnění průběhu experimentů na různých pracovištích. Hierarchický systém předpokládá zapojení více počítačů – první přímo komunikuje s experimentálním zařízením, předzpracovává získaná data a ta pak k výslednému zpracování posílá do nadřazeného počítače nebo počítačů. Určité rysy tohoto systému můžeme pozorovat i u současného systému jednoho uživatele, neboť mikropočítače ve fyzikální laboratoři bývají připojeny k Internetu a ten v případě potřeby může přenos dat na silnější výpočetní techniku zprostředkovat. Zvláštním případem je použití speciálních procesorů. V předchozích třech případech byly jako řídící a výpočetní technika použity univerzální počítače různé výkonnosti, které kromě řízení experimentu mohly provádět i nejrůznější další činnosti. Toto řešení, i když je v mnohých ohledech výhodné (cena, dostupnost, atd.), může mít slabinu v omezené výkonnosti. Pokud je součástí zpracování dat z experimentu provádění větší množství dosti složitých výpočetních operací, bývá někdy výhodnější použít specializovaný mikroprocesorový systém, který sice ztrácí univerzálnost ale zato daný typ operací dokáže provádět mnohem efektivněji. Příkladem jsou např. spektrometry vybavené specializovaným procesorem pro provádění Fourierovy transformace na hardwarové úrovni.
Do tohoto okruhu znalostí patří i používání mikroprocesorové techniky při řízení včetně základů programování mikroprocesorů. Všechny tyto znalosti opět překračují základní kurs a budou jim věnovány specializované přednášky.
KAPITOLA 5. DALŠÍ SMĚRY POČÍTAČOVÉ FYZIKY
5.2
68
Symbolické manipulace
Jedná se partii počítačové fyziky, svým zadáním blížší spíše klasické než numerické matematice. Cílem je například spočítat integrál z nějaké funkce, ale na rozdíl od numerické matematiky, kde výstupem výpočtu je konkrétní číslo při zadaných konkrétních mezích integrace, zde požadujeme nalezení primitivní funkce k námi zadané funkci. Podobně tomu je i v dalších matematických oblastech, kde od počítače požadujeme odpověď takového typu, jako se hledá při cvičeních v matematické analýze nebo algebře. K tomuto účelu byly vypracovány velké specializované komerční programy jako je MATHEMATICA, REDUCE nebo MAPLE a fyzik se stává pouhým jejich uživatelem. Podobně jako pro experimentální fyziky je důležitou partií řízení experimentu, existují oblasti fyziky, kdy symbolické manipulace jsou považovány za jednu z nejpodstatnějších partií počítačové fyziky – typickým příkladem je teoretická fyzika.
Kapitola 6 ZÁVĚR Předložený studijní text se věnoval pouze některým partiím počítačové fyziky, přičemž některé probíral hlouběji a jiné jen na informativní úrovni. Výběr podstatných a méně podstatných partií byl podřízen dvěma kritériím: – použitelnost pro další studium včetně vypracování dizertačních prací a pro budoucí zaměstnání, – neexistence nebo nedostatečná existence vhodných komerčních programových souborů a z toho vyplývající nutnost vytváření vlastních programů. Podle obou kritérii se jako nejpodstatnější části počítačové fyziky jeví techniky počítačového modelování následované zpracováním obrazu, a proto obě tyto partie tvoří jádra obou dílů studijního textu. Tento studijní text se zabývá pouze okruhy té části počítačové fyziky, které se začíná říkat klasická počítačová fyzika. Vedle ní se v současnosti rodí zcela nová oblast nazývaná moderní počítačová fyzika, tvořená partiemi jako neuronové sítě. fuzzy logika nebo evoluční programování a genetické algoritmy, resp. použitím tohoto aparátu pro řešení fyzikálních problémů. Z metodik zmíněných v našem studijním textu se začíná mezi tyto moderní partie počítat i waveletová transformace, původně jedna o technik integrálních transformací. I když i v klasické počítačové fyzice vznikají stále nové a silnější algoritmy, v oblasti moderní počítačové fyziky tento vývoj probíhá s ještě větší intenzitou. Proto též některé moderní partie počítačové fyziky budou postupně zapojovány do výuky vyšších ročníků studia počítačové fyziky a počítačových metod ve vědě a technice. Důvodem tohoto postupu je, že dílčí znalosti z těchto moderních partií se mohou vhodně kombinovat s již standardními postupy počítačové fyziky a tak pomoci efektivněji řešit nejrůznější fyzikální problémy. Jedním z typických příkladů je použití neuronových sítí při zpracování obrazu, prováděné na našem pracovišti [21]. Podobně se jako silný prostředek při řešení zpracování obrazů ve fyzice projevila i waveletová transformace.
69
Literatura [1] Žára J., Limpouch A., Beneš B., Werner T.: Počítačová grafika – principy a algoritmy, Grada, Praha 1992. [2] Žára J., Beneš B., Felkel P.: Moderní počítačová grafika, Computer Press, Praha 1998. [3] Rogers D.F.: Procedural Elements for Computer Graphics, McGraw-Hill Book Company, New York 1989. [4] Bresenham J.E.: Algorithm for Computer Control of a Digital Plotter, IBM Systems Journal 4 (1965), p. 25. [5] Hrach R., Hrachová V., Vicher M., Šimek J., Sedlák D.: Study of Plasma-Solid Interaction in Electronegative Plasma, Le Vide 307 (2003) – Numéro Speciaux, p. 103. [6] Šimek J., Hrach R.: Image Analysis in Physics – Sensitivity Analysis of Morphological Methods and their Application in Thin Film Physics, Thin Solid Films (2004) – v tisku. [7] Hrach R., Novotný D., Novák S., Pavlík J.: Multilevel Morphological Analysis of Continuous Films and Surfaces, Thin Solid Films 433 (2003), p.133. [8] Pratt W.K.: Digital Image Processing, John Wiley and Sons, New York 1991. [9] Serra J.: Image Analysis and Mathematical Morphology, Academic Press, London 1982. [10] Stauffer D., Aharony A.: Introduction to Percolation Theory, Taylor and Francis Ltd., London 1991. [11] Hrach R., Šimek J., Kostern M.: Study of Thin Film Growth by means of Computer Simulation and Image Analysis, Vacuum 67 (2002), p. 229. [12] Hrach R., Novotný D., Novák S., Pavlík J., Hrachová V.: Radial Distribution Function and Nearest-Neighbour Distribution in Spatial Analysis of 2D Objects, Proc. 6th Intern. Conf. „Physics Computing 94ÿ, Lugano 1994, Switzerland, p. 333. [13] Hrach R., Novotný D., Novák S., Pavlík J.: Morphological Analysis of Continuous Metal Films, Vacuum 57 (2000), No 3, p. 259. 70
LITERATURA
71
[14] Hrach R., Novotný D., Novák S., Pavlík J.: Analysis of Methods for the Study of Semicontinuous and Continuous Metal Film Morphology, Proc. 7th Joint Vac. Conf., Debrecen 1997, Hungary, p. 237. [15] Ripley B.D.: Spatial Statistics, John Wiley Sons Inc., New York 1981. [16] Hrach R., Novotný D., Vicher M.: Popis stupně uspořádanosti soustav s využitím klasických morfologických metod, ERGO 3 (2001), No 3, p. 71. [17] Novák S., Hrach R.: Rekonstrukce třírozměrných struktur z jejich obrazů, ERGO 4 (2002), No 3, p. 41. [18] Coley J.W., Tookey J.W.: An Algorithm for Machine Calculation of Complex Fourier Series, Math. of Computation 19 (1965), p. 297. [19] Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P.: Numerical Recipes in Fortran 77, Second Edition, The Art of Scientific Computing, Cambridge Univ. Press, Cambridge 2002. [20] Čížek V.: Diskrétní Fourierova transformace a její použití, SNTL, Praha 1981. [21] Malý M., Hrach R., Novotný D.: Analýza rozložení objektů s využitím neuronových sítí, ERGO 3 (2001), No 3, p. 83.