3 Předzpracování Obsah Předchozí kapitola Následující kapitola Titulní strana

          Předzpracování je velice důležitá součást digitálního zpracování obrazu. Ve své diplomové práci této oblasti věnuji poměrně hodně prostoru, protože na výsledcích předzpracování výrazně závisí kvalita rozkladu na tiskové barvy a ta hraje rozhodující úlohu při rozpoznávání. Předzpracování je možné rozdělit do dvou samostatných oblastí:

          Předzpracování pracuje s barevnou bitmapou a jeho výsledkem je opět barevná bitmapa. Pro účely své diplomové práce jsem použil pouze relevantní část metod předzpracování, ty jsem upravil pro dané účely a jejich výsledky demonstroval na digitalizovaných mapách.


3.1 Restaurace obrazu

          Restaurací obrazu rozumíme odstranění degradací vzniklých při převodu mapy vytištěné na papíře do digitalizované podoby. Pro převod dokumentů do digitální podoby se používají digitizéry - obvykle to bývá scanner nebo kamera. V dalším textu budu předpokládat, že mapy byly digitalizovány na stolním scanneru, druhy degradací totiž úzce souvisí s druhem použitého digitizéru.

          Digitalizaci map jsem prováděl v rozlišení 600 x 600 dpi (viz kapitola "2.1 Měřítko a orientace mapy") nezávisle na dvou stolních scannerech:

          Ve výsledných obrázcích získaných pomocí uvedených zařízení bylo patrné jednak mírné rozmazání detailů a také relativní posun RGB vrstev (barevné skvrny na okrajích objektů). Uvedené degradace se projevovaly při použití obou scannerů, ale v různé míře. Obě tyto degradace výrazným způsobem ovlivňují klasifikaci barev a proto je nezbytné je alespoň částečně odstranit. Pro zobrazení degradací a výsledků jejich odstraňování jsem použil výřez mapy [M8] velikosti 128 x 128 pixelů. Na obrázku 3.1 je zobrazen výsledek digitalizace.

Degradace obrazu způsobené digitalizací na stolním scanneru (originál).
originál
Degradace obrazu způsobené digitalizací na stolním scanneru (výsledek digitalizace).
výsledek digitalizace

Obrázek 3.1
Degradace obrazu způsobené digitalizací na stolním scanneru.

          Na tomto místě bych chtěl ještě upozornit na chybný postup při digitalizaci, kterého je třeba se vyvarovat. Po digitalizaci mapy je vždy nutné zkontrolovat nastavení jasu a kontrastu - správné nastavení lze ověřit tak, že výsledek neobsahuje žádný pixel, jehož některá ze složek RGB nabývá buď minimální nebo maximální hodnoty. V opačném případě je pravděpodobné, že skutečné intenzity byly mimo rozsah a byly upraveny (a potom jsou barvy příslušných pixelů zkreslené). Na druhou stranu je žádoucí, aby intenzity pokrývaly celou škálu hodnot (mimo již zmíněných zcela krajních hodnot) - v takovém případě dosáhneme minimální chyby způsobené kvantováním.


3.1.1 Korekce posunu RGB vrstev

          Velmi nepříjemnou, ale snadno odstranitelnou degradací je vzájemný posun RGB vrstev. Scanner používá tři nezávislé sady fotocitlivých snímačů, které mohou být mírně posunuty vůči jejich korektní poloze. Při velikém rozlišení jsou potom výsledkem barevné skvrny na okrajích objektů. U obou použitých scannerů byl při rozlišení 600 x 600 dpi posun menší než 1 pixel.

          Úloha korekce posunu RGB vrstev sestává ze dvou částí. Nejprve je nutné odhadnout parametry transformace a poté transformaci provést. Nejprve uvedu způsob provedení transformace a teprve potom odhad jejích parametrů, protože odhad využívá samotné transformace.


Transformace

          Vstupem transformace je degradovaný obraz a sada koeficientů Rx, Ry, Gx, Gy, Bx, a By. Koeficienty udávají, o kolik pixelů se mají posunout roviny R, G a B ve směru os x a y. Koeficienty jsou reálná čísla, protože posun vrstev bývá zpravidla jenom zlomek velikosti pixelu. Není nutné používat všechny koeficienty, je možné jednu rovinu pevně fixovat (zpravidla rovinu R) - potom budou koeficienty Rx a Ry nulové.

          Transformace je speciálním případem převzorkování se zachováním měřítka. Převzorkování se provádí odděleně ve vrstvách R, G a B pomocí bilineární interpolace, což při zachování měřítka plně postačuje. Ve vzorcích pro převzorkování je označen vstupní obraz jako g a výsledek transformace jako f. Transformační vzorec je uveden pouze pro vrstvu R, transformace v ostatních vrstvách probíhá shodně.

          Rovnice 3.1 Rovnice 3.1 (3.1)


Odhad parametrů transformace

          Odhad parametrů lze v zásadě provést dvojím způsobem. Buď parametry odhadneme podle výsledku digitalizace (což je možné) nebo můžeme použít odhad automatický. Nejlepších výsledků lze zpravidla dosáhnout kombinací obou zmíněných způsobů - nejprve je třeba odhadnout maximální hodnotu parametrů, poté lze aplikovat automatickou detekci parametrů a výsledek případně ještě nakonec ručně poopravit.

          Pro automatický odhad parametrů transformace je třeba nalézt oblast obrazu, v níž se vyskytují (resp. mají vyskytovat) pouze monochromatické barvy. Vlastnosti vybrané oblasti vyplynou z algoritmu pro odhad parametrů.

          Vstupem algoritmu pro automatický odhad parametrů je obdélníková (není podmínkou) oblast obrazu <xmin, xmax> x <ymin, ymax> a odhad maximálního posunu vrstev (pro všechny vrstvy společný). Výstupem algoritmu je sada koeficientů Rx, Ry, Gx, Gy, Bx a By.

          Nejprve algoritmus fixuje vrstvu R a výsledné koeficienty Rx a Ry položí nulové. Poté aplikuje transformaci na vrstvu G pro všechny možné kombinace parametrů Gx a Gy v rozsahu stanoveném odhadem maximálního posunu (posuny procházíme po zadaném kroku, obvykle postačuje 0,1 pixelu). Symbol g značí degradovaný obraz, funkce T(g.G, Gx, Gy) značí posun vrstvy G obrazu g o Gx horizontálně a Gy vertikálně. Pro každou provedenou transformaci spočte funkci CG určující krytí vrstev R a G:

          Rovnice 3.2 (3.2)

          Koeficienty Gx a Gy budou ty, pro které vyšla minimální hodnota CG. Analogický postup aplikujeme na vrstvu B obrazu g (výpočet krytí provedeme opět vůči vrstvě R). Referenční oblast je nutné vybrat tak, aby funkce CG (CB) nabývala minima pro korektní koeficienty Gx a Gy (Bx a By). Z toho důvodu je kromě již zmíněné podmínky monochromatických barev nutné, aby v oblasti bylo co nejvíce přechodů a to jak v horizontálním, tak i ve vertikálním směru. Těmto požadavkům obvykle nejlépe vyhovuje výřez názvu mapy nebo nějaké reklamy.

          Praktické výsledky korekce posunu RGB vrstev jsou demonstrovány na obrázku 3.2.

Korekce posunu RGB vrstev (výsledek digitalizace).
výsledek digitalizace
Korekce posunu RGB vrstev (korekce posunu RGB vrstev).
korekce posunu RGB vrstev

Obrázek 3.2
Korekce posunu RGB vrstev.


Další informace

          S řešením obdobné degradace jsem se v literatuře nesetkal. V knize [4] je naznačen odhad parametrů transformací a bilineární interpolace.


3.1.2 Impulsní odezva scanneru

          V této a následujících kapitolách se budu držet následující konvence pro označování obrazů:


Digitalizace obrazu

          Digitalizací obrazu budeme rozumět převod reálného obrazu (mapy), který je dvojrozměrnou funkcí, do diskrétní podoby. Při digitalizaci budeme používat výhradně obdélníkové (prakticky čtvercové) vzorkování s počátkem v [x0, y0] a obdélníkem velikosti Dx x Dy. V ideálním případě by intenzita pixelu o souřadnicích [m, n] závisela pouze na obsahu oblasti ohraničené obdélníkem pro tento pixel:

          Rovnice 3.3 (3.3)

          Při reálné digitalizaci je ale skutečnost jiná. Díky vlastnostem scanneru intenzita pixelu závisí i na pixelech v jeho okolí - jejich váhu určuje funkce h. Funkci h budeme nazývat impulsní odezvou scanneru.

          Rovnice 3.4 (3.4)

          Funkce 3.4 představuje polohově invariantní zobrazovací model - výsledek funkce h závisí pouze na relativní poloze vůči středu pixelu. Tento model je vyžadován metodami pro rekonstrukci obrazu (navíc vystihuje skutečnost), proto se obecným modelem zabývat nebudeme. Také je třeba si uvědomit, že funkce pro ideální vzorkování je speciálním případem reálné vzorkovací funkce. Bude-li mít totiž impulsní odezva h tvar obdélníkového impulsu, funkce 3.3 a 3.4 budou mít stejný výsledek.

          Rovnice 3.5 (3.5)

          Nás ale bude zajímat především diskrétní podoba reálného vzorkování. Impulsní odezva scanneru h je sice definována v intervalu <-inf, inf> x <-inf, inf>, ale nezanedbatelných hodnot nabývá prakticky pouze na omezené oblasti v blízkosti nuly velikosti <-M, M> x <-N, N>. Rovnice 3.4 přepsaná do diskrétní podoby s konečnou impulsní odezvou bude mít tvar:

          Rovnice 3.6 (3.6)

          Protože se jedná o diskrétní konvoluci, můžeme zjednodušeně psát g = f*h. Příklad interakce sousedních pixelů při digitalizaci je uveden na obrázku 3.3. V ideálním případě by měl červeně orámovaný pixel bílou barvu, v reálném případě jeho intenzita závisí na pixelech v okolí velikosti 5 x 5 pixelů (M = N = 2).

Impulsní odezva scanneru.

Obrázek 3.3
Impulsní odezva scanneru.


Odhad impulsní odezvy

          Existují metody, které ze znalosti impulsní odezvy scanneru a reálně digitalizovaného obrazu dokážou určit ideálně digitalizovaný obraz. Tento druh korekce je při zpracování barevných map velmi důležitý, protože prvním krokem před vlastním rozpoznáváním je klasifikace barev. Rozmazání obrazu a přetiskování barev má za následek vysoký počet odstínů barev a jejich obtížnou rozlišitelnost. Pro použití zmíněných metod je ale nutné znát impulsní odezvu h. Na tomto místě bych chtěl upozornit, že impulsní odezva h je při snímání na stejném zařízení a se stejnými parametry neměnná.

          První velmi jednoduchou metodou odhadu impulsní odezvy je digitalizovat izolovaný bod. Předem je nutné odhadnout velikost jádra konvoluce (obvykle postačuje velikost 7 x 7, tj. N = M =3). Dále je třeba na kvalitní tiskárně vytisknout obrázek obsahující souvislou bílou a černou plochu a izolovaný bod. Všechny tyto objekty by od sebe měly být dostatečně vzdáleny. Posledním krokem je digitalizace vytištěného obrázku. Tisk i digitalizaci je nutné provádět ve stejném rozlišení, v jakém jsou digitalizovány obrazy, které je třeba restaurovat. Struktura obrázku a výsledek digitalizace je zobrazen na následujícím obrázku:

Digitalizace izolovaného bodu.

Obrázek 3.4
Digitalizace izolovaného bodu.

          Oblast v níž se objevila impulsní odezva označím jako r, r0,0 bude ve středu této oblasti. W bude označovat intenzitu pixelu s bílou barvou, B intenzitu pixelu s černou barvou. Impulsní odezvu již snadno spočteme:

          Rovnice 3.7 (3.7)

          Výhodou této metody její jednoduchost. Potíže ale mohou vzniknout při pořízení speciálního obrázku. Izolovaný pixel vytištěný na reálné tiskárně zcela jistě nebude ideální - proto doporučuji obrázek dlaždicově zopakovat a na výtisku vybrat nejlepší (porovnání výsledků lze nejlépe provést pod mikroskopem). Druhou potíží je relativní posun hranic vytištěného pixelu a snímaného pixelu - zde doporučuji provést více digitalizací při různých polohách papíru, nebo vytisknout pixel v menším rozlišení, než je použito pro snímání.

          Jinou podobnou metodou je odhad odezvy pomocí vertikální nebo horizontální ostré hrany. Tato metoda je sice pracnější na výpočet, ale nevyžaduje zvláštní obrázek (běžné mapy obsahují vyhovujících objektů celou řadu). Další výhodou této metody je její větší stabilita.

Digitalizace ostré hrany.

Obrázek 3.5
Digitalizace ostré hrany.

          V digitalizovaném obrazu vyhledáme posloupnost intenzit odpovídající digitalizované hraně a pomocí této posloupnosti již snadno dopočítáme koeficienty pro jednorozměrnou konvoluci h0 - hk.

          Rovnice 3.8 (3.8)

          Uvedený postup doporučuji provést několikrát, výsledky zanést do grafu a z grafu odhadnout průběh impulsní odezvy (viz obrázek 3.5). Jako poslední krok je nutné rozšířit jednorozměrnou konvoluci na dvojrozměrnou. Nejjednodušším způsobem je předpokládat, že její průběh je ve všech směrech stejný a z grafu odhadnout hodnoty podle vzdálenosti od počátku. Jinou možností je určit jednorozměrnou konvoluci ve vertikálním a horizontálním směru a dvojrozměrný výsledek vypočítat jako jejich kombinaci.

          Příkladem impulsní odezvy může být výsledek této metody aplikovaný na mapu digitalizovanou na scanneru HP ScanJet IIcx v rozlišení 600 x 600 dpi.

          Rovnice 3.9 (3.9)

          Další možností je předpokládat, že impulsní odezvu lze aproximovat průběhem nějaké funkce. Podstatnou výhodou je, že místo několika desítek koeficientů stačí odhadnout pouze parametry této funkce (mnohdy bývá pouze jeden). Protože je parametrů málo, odhad lze provést zkusmo na základě výsledků restaurace. Obvykle je možné tímto způsobem dosáhnout alespoň částečného zaostření. Příkladem takových funkcí může být normální rozložení N nebo normální rozložení s nulovou střední hodnotou N0:

          Rovnice 3.10 (3.10)

          Jistě by bylo možné aplikovat i různé variace uvedených metod. Prakticky jsem otestoval všechny uvedené způsoby, pro restauraci používám odhadu pomocí ostrých hran.


Další informace

          Impulsní odezva zařízení se ve zpracování obrazu nevyskytuje pouze v souvislosti s degradacemi vzniklými při digitalizaci. Se znalostí impulsní odezvy je například možné restaurovat fotografie poškozené pohybem nebo rozostřením.

          Digitalizací obrazu se zabývá prakticky každá kniha o digitálním zpracování obrazu (např. [6] nebo [4]), kapitoly o konvolucích a impulsní odezvě zařízení jsou obsaženy v [6].


3.1.3 Zpětná konvoluce pomocí soustavy lineárních rovnic

          Nyní máme k dispozici reálně digitalizovaný obraz g a impulsní odezvu zařízení h. Dále známe závislost g = f*h. Intuitivním způsobem získání ideální digitalizace f je řešení soustavy rovnic získané pomocí definice konvoluce.


Soustava rovnic

          Soustavu rovnic získáme přímo rozepsáním výpočtu pro jednotlivé pixely z definice konvoluce. Tato soustava bude ale značně veliká, protože pro každý pixel obrazu bude obsahovat jednu rovnici. Technickou záležitostí bude vyřešit okrajové podmínky. Podle definice konvoluce závisí hodnoty výsledného pixelu na pixelech z jeho okolí a tedy přesnou definici nebude možné použít na okrajové body obrazu. Také není možné použít rovnice pouze pro ty body, které jsou dostatečně vzdáleny od okrajů, protože potom by soustava obsahovala méně rovnic, než proměnných. Zvolil jsem tedy takové řešení, kde předpokládám, že intenzita pixelů za okrajem bude shodná s intenzitou pixelů na okraji. Tato úprava může způsobovat chyby na okrajích restaurovaného obrazu.

          Potíží při aplikaci všech metod dekonvoluce je, že obraz g je zatížen aditivním šumem z důvodu kvantování intenzit (viz obrázek 3.6). Proto není žádoucí, abychom nalezli skutečné řešení soustavy rovnic - výsledek by neodpovídal očekávání. Prakticky jsem používal iterační metody vycházející z řešení f shodného s pravou stranou soustavy g (degradovaný obraz). Tyto metody postupně zpřesňovaly řešení f od degradovaného obrazu přes částečně restaurovaný obraz až po výsledek zatížený šumem.

Aditivní šum způsobený kvantováním (digitalizovaný obrázek).
digitalizovaný obrázek
Aditivní šum způsobený kvantováním (aditivní šum).
aditivní šum

Obrázek 3.6
Aditivní šum způsobený kvantováním.

          Je-li K x L velikost vstupního obrazu, potom označím P = KL počet pixelů vstupního obrazu. Matice H má velikost P x P, vektory f, g a b mají velikost P. Vektory g a f odpovídají vstupnímu a ideálně digitalizovanému obrazu zapsanému po řádcích, vektor b určuje korekci pro okrajové pixely.

          Rovnice 3.11 (3.11)

          Matice E obsahuje pouze koeficienty z impulsní odezvy h. Na každém řádku (a v každém sloupci) je nejvýše (2M+1)(2N+1) nenulových hodnot pravidelně rozmístěných v okolí hlavní diagonály. Pro impulsní odezvu velikosti 3 x 3 (M = N = 1) má E má tvar:

          Rovnice 3.12 (3.12)

          Z předchozího příkladu je zřejmé, že matice E pro obecné M a N bude obsahovat 2N+1 diagonálních pásů složených z matic Ai, matice Ai bude potom obsahovat 2M+1 diagonálních pásů. Nyní zbývá ještě definovat hodnoty vektoru b. Tento vektor kompenzuje koeficienty matice E chybějící na daném řádku s využitím nejbližšího pixelu od chybějícího pixelu gnearest.

          Rovnice 3.13 (3.13)


Řešení soustavy rovnic

          Jak jsem již zmínil, dobrých výsledků jsem dosáhl při použití iteračních metod - Jacobiho, Gauss-Seidelovy nebo Southwellovy. Všechny zmíněné metody ale kladou podmínku na diagonální dominanci matice E a tím tedy i na impulsní odezvu h:

          Rovnice 3.14 (3.14)

          Jinými slovy absolutní hodnota h0,0 musí být alespoň taková jako součet absolutních hodnot ostatních prvků h. Uvedená podmínka je postačující, není ale nutná - to znamená, že iterační metody mohou konvergovat i při jejím nesplnění. Podmínka je ale příliš omezující a reálné impulsní odezvy takové vlastnosti zpravidla nemají, možností je ale použít rozkladu h na více konvolucí a iterační metodu požít vícekrát.

          I přes omezující podmínku 3.14 jsou tyto iterační metody výhodné hlavně pro svou rychlost a dobré výsledky. Další velikou výhodou těchto metod je, že nemodifikují matici E, ta nemusí být uložena v paměti. Její prvky lze totiž snadno vypočíst:

          Rovnice 3.15 (3.15)

          Prakticky jsem aplikoval Jacobiho a Gauss-Seidelovu iterační metodu. Jacobiho iterační metoda se provádí podle následujícího postupu:

          Rovnice 3.16 (3.16)

          Při použití Jacobiho metody se vždy spočítá celý vektor řešení pomocí výsledku z minulého kroku. Podobná Gauss-Seidelova iterační metoda používá spočtené koeficienty ihned (při výpočtu kroku (k+1) se používají pro výpočet fi koeficienty fj z kroku (k) i z kroku (k+1)):

          Rovnice 3.17 (3.17)

          Gauss-Seidelova metoda je paměťově úspornější, protože vyžaduje pouze jeden vektor řešení (Jacobiho metoda vyžaduje vektor řešení z minulého kroku). Samozřejmostí je, že při praktickém výpočtu je nutné k výpočtu použít místo P pouze MN nenulových koeficientů.

          U všech uvedených iteračních metod navíc nutné testovat podmínku na ukončení iterace. K tomuto je možné použít dva způsoby: buď provést zadaný počet iterací, nebo sledovat chybu (reziduum) a iteraci ukončit až při poklesu této hodnoty pod určitou hranici. Prakticky jsem použil první metodu - lze snadněji sledovat postup zaostření. Praktické výsledky zaostření pomocí této metody jsou demonstrovány na obrázku 3.7. V tomto případě nesplňovala impulsní odezva podmínku 3.14 a při použití iteračních metod muselo být jádro konvoluce upraveno. Obrázek ukazuje postupné zaostření při několikanásobné aplikaci.

Zpětná konvoluce pomocí soustavy lineárních rovnic (digitalizovaný obrázek).
digitalizovaný obrázek
Zpětná konvoluce pomocí soustavy lineárních rovnic (1x aplikovaná dekonvoluce).
1x aplikovaná dekonvoluce
Zpětná konvoluce pomocí soustavy lineárních rovnic (2x aplikovaná dekonvoluce).
2x aplikovaná dekonvoluce

Obrázek 3.7
Zpětná konvoluce pomocí soustavy lineárních rovnic.


Další informace

          Postup výpočtu zpětné konvoluce pomocí soustavy lineárních rovnic jsem v použité literatuře nenalezl. Uvedené metody pro řešení soustav lineárních rovnic byly převzaty z [1]. Kromě popsaných metod existují i výhodnější iterační metody, které nevyžadují splnění podmínky 3.14. Protože jsem se tímto způsobem zpětné konvoluce zabýval pouze okrajově, nevyhledával jsem optimálnější iterační metody.


3.1.4 Diskrétní Fourierova transformace

          Fourierova transformace (FT) je důležitý nástroj používaný pro převod obrazu do jeho sinových a cosinových komponent. Fourierova transformace nalézá uplatnění v řadě aplikací, z nichž nedůležitější jsou analýza, komprese a filtrace obrazu. Pro účel diplomové práce jsem ji použil při restauraci obrazu, což je speciální případ filtrace.

          V oblasti digitálního zpracování obrazu nalézá uplatnění diskrétní verze Fourierovy transformace (DFT), především pak dvojrozměrná diskrétní Fourierova transformace (2D DFT), proto se spojitým případem zabývat nebudu.

          Výsledek DFT neobsahuje všechny frekvence formující obraz, ale pouze sadu frekvencí, která postačuje k přesnému popisu navzorkované verze obrazu. Počet těchto frekvencí souvisí s počtem pixelů obrazu - velikost obrazu v prostorové oblasti a ve frekvenční oblasti je stejná. 2D DFT převádí čtvercový obraz velikosti NxN pixelů na pole komplexních čísel stejné velikosti. Reálná složka výsledku udává zastoupení funkce cos pro danou frekvenci, imaginární složka zastoupení funkce sin.

          Protože jsem 2D DFT implementoval v aplikaci, kromě definice a příkladů uvedu i podrobný rozbor použitých metod výpočtu.

Prostorová a frekvenční oblast (prostorová oblast).
prostorová oblast
Prostorová a frekvenční oblast (frekvenční oblast (amplituda)).
frekvenční oblast (amplituda)
Prostorová a frekvenční oblast (frekvenční oblast (fáze)).
frekvenční oblast (fáze)

Obrázek 3.8
Prostorová a frekvenční oblast.

          Úvodem upozorním na některé implementační detaily:


Definice diskrétní Fourierovy transformace

          Dvojrozměrná diskrétní Fourierova transformace a inverzní dvojrozměrná diskrétní Fourierova transformace jsou definovány následujícím způsobem:

          Rovnice 3.18 (3.18)

          Protože výpočet DFT a inverzní DFT je analogický, v dalším textu se budu zabývat pouze výpočtem DFT.

          Z definice je zřejmé, že pro výpočet jednoho koeficientu je třeba O(N2) komplexních násobení a sčítání. Celkem je potřeba spočítat N2 koeficientů, takže celková složitost výpočtu je O(N4). Tento způsob výpočtu je použitelný pouze pro velmi malé velikosti obrazu. Protože je FT separabilní, můžeme ji rozepsat pomocí dvou jednorozměrných FT:

          Rovnice 3.19 (3.19)

          V případě, že požadujeme pouze jeden koeficient Fx,y, žádné úspory nedosáhneme. Pokud ale požadujeme všechny koeficienty, nejprve spočteme jednorozměrné DFT po řádcích a potom z výsledku spočteme jednorozměrné DFT po sloupcích. V každém případě je potřeba spočítat N2 koeficientů 1D DFT, výpočet jednoho koeficientu je možné provést v čase O(N). Celková složitost výpočtu tedy klesla z O(N4) na 2N2.O(N) = O(N3).

          V dalším odstavci ukáži, jakým způsobem je možné spočítat všechny koeficienty 1D DFT v čase O(N.log2(N)), složitost 2D DFT klesne na O(N2.log(N)).


Metody rychlého výpočtu Fourierovy transformace

          V roce 1965 Cooley a Tukey nalezli rychlý algoritmus pro výpočet DFT (podobné algoritmy byly publikovány již několikrát od roku 1903, ale svého času nenalezly ohlas). Dnes je známa celá řada modifikací tohoto algoritmu, které se souhrnně označují jako rychlá Fourierova transformace (fast Fourier transform - FFT). Rychlé algoritmy výpočtu DFT lze interpretovat jako faktorizaci jádra této transformace na součin řídkých matic, čímž se vyloučí většina zbytečných operací.

          V minulém odstavci jsem ukázal, jak spočítat 2D DFT pomocí 2N 1D DFT. Dále se budu zabývat pouze výpočtem jednorozměrných DFT. Podrobně ukážu princip výpočtu jednoho ze dvou základních algoritmů: rozkladu v časové oblasti. Pokud máme vstupní posloupnost fx sudé délky, můžeme DFT rozepsat následujícím způsobem:

          Rovnice 3.20 (3.20)

          Dále definujeme dvě posloupnosti gx a hx poloviční délky. Posloupnost gx bude obsahovat sudé prvky posloupnosti fx, posloupnost hx bude obsahovat liché prvky posloupnosti fx.

          Rovnice 3.21 (3.21)

          Obě posloupnosti gx a hx mají samozřejmě své DFT.

          Rovnice 3.22 (3.22)

          Nyní vyjádříme Fx pomocí Gx a Hx.

          Rovnice 3.23 (3.23)

          Za povšimnutí stojí, že při výpočtu Fx jsou použity hodnoty Gx a Hx s indexem x mimo rozsah <0, N/2>. Zde využijeme toho, že DFT je periodická a tedy Gx = Gx mod N/2 a Hx = Hx mod N/2.

          Princip výpočtu FFT spočívá v rekurzivním dělení posloupnosti až na posloupnosti délky 1 - potom složitost klesne z O(N2) na O(N.log(N)). DFT posloupnosti délky 1 lze snadno vypočíst:

          Rovnice 3.24 (3.24)

          Dále bych chtěl upozornit na některé implementační detaily:


Implementace výpočtu FFT

          Popis uvedený v minulém odstavci přímo dává návod, jak FFT implementovat. Pro názornost uvedu příklad na posloupnosti délky 8 znázorněný pomocí grafu. Koeficient u šipky značí komplexní násobení příslušné hodnoty, znaménko + v uzlu značí součet dvou komplexních čísel. Výsledné koeficienty F0 - F7 je nakonec ještě třeba vydělit odmocninou z N.

          Z grafu je zřejmé, že budeme potřebovat dvě pole délky N, jedno pro vstup a druhé pro výsledek na dané úrovni. Vstupem vyšší úrovně je výsledek z minulé úrovně. Na každé úrovni se provádí N komplexních sčítání a násobení, úrovní je log2(N).

FFT pro N = 8 s přirozeným uspořádáním vstupu.

Obrázek 3.9
FFT pro N = 8 s přirozeným uspořádáním vstupu.

          Kromě výše zmíněného existuje ale výhodnější způsob implementace, ve kterém jsou nároky na operační paměť poloviční. Základní myšlenkou je změna uspořádání vstupní posloupnosti - posloupnost je nyní uspořádána podle bitově převráceného indexu (např. 001 --> 100). Na každé úrovni vždy 2 vstupní hodnoty dávají za výsledek 2 výstupní hodnoty se stejnými indexy - není tedy nutné pole velikosti N pro výsledek, ale postačí pouze jedna pomocná komplexní proměnná.

FFT pro N = 8 s bitově převráceným uspořádáním vstupu.

Obrázek 3.10
FFT pro N = 8 s bitově převráceným uspořádáním vstupu.

          Při použití rozkladu ve frekvenční oblasti je vstupní posloupnost uspořádána přirozeně, výstupní posloupnost je naopak uspořádána bitově převráceně.


Další informace

          Při implementaci mi nešlo o to použít co nerychlejší algoritmus, ale implementovat algoritmus pracující v reálném čase. Největším nedostatkem použité metody je nutnost doplnění obrazu na čtverec o velikosti hrany 2m - existují algoritmy, které umožňují výpočet pro obdélníkový vstup bez doplnění. Techniky rychlého výpočtu DFT prochází dosud vývojem.

          Zběžné seznámení s principem Fourierovy transformace obsahuje většina knih o digitálním zpracování obrazu. Podrobný popis spojité a diskrétní Fourierovy transformace lze nalézt v [6], praktické výsledky předkládá interaktivní kurs [2]. Základní seznámení s algoritmy pro rychlý výpočet DFT obsahuje [7].


3.1.5 Inverzní filtrace

          Po seznámení s Fourierovou transformací můžeme přistoupit k metodám restaurace obrazu založených výpočtu ve frekvenční oblasti. Metody jsou založeny na konvolučním teorému (operátor F značí Fourierovu transformaci):

          Rovnice 3.25 (3.25)


Inverzní filtr

          Nejjednodušší z restauračních filtrů pracujících ve frekvenční oblasti je inverzní filtr. Stejně jako zpětná konvoluce pomocí soustavy lineárních rovnic využívá znalost impulsní odezvy h a lineárního zobrazovacího modelu bez šumu g = f * h. Vychází přímo z konvolučního teorému:

          Rovnice 3.26 (3.26)

          Na obrázku 3.11 jsou výsledky aplikace této metody. Je zřejmé, že výsledek restaurace neodpovídá očekávání. Důvodem je aditivní šum způsobený kvantováním (viz obrázek 3. 6). Použití této metody je omezeno pouze na uměle vytvořené případy uložené pomocí reálných čísel, které nejsou zatíženy šumem. Nicméně je tato metoda základem pro další použitelné metody.

Inverzní filtr (digitalizovaný obraz).
digitalizovaný obraz
Inverzní filtr (inverzní filtr).
inverzní filtr

Obrázek 3.11
Inverzní filtr.


Pseudo-inverzní filtr

          Důvodem, proč inverzní filtr nefunguje podle očekávání, je přítomnost aditivního šumu. Ten způsobuje zastoupení vysokých frekvencí i v obraze g, ačkoli detaily by měly být díky impulsní odezvě h rozmazány a vysoké frekvence by měly být zastoupeny pouze minimálně. Naopak v reálné impulsní odezvě nejsou vysoké frekvence příliš zastoupeny a tedy filtrace podle 3.26 zesiluje šum. Proto výsledek filtrace neodpovídá očekávání.

          Vysvětlení příčiny chyb dává návod, jak zvýraznění šumu předejít. Aby se zabránilo růstu vysokých frekvencí, zavedeme koeficient T a výpočet podle 3.26 upravíme:

          Rovnice 3.27 (3.27)

          Výsledky pseudo-inverzní filtrace pro různé hodnoty T jsou ukázány na obrázku 3.12. Na obrázku je vidět, že se sice podařilo odstranit šum, ale výsledky stále ještě nejsou dobré. Také je třeba si uvědomit, že položíme-li T = 0, získáme inverzní filtr, s roustoucím T bude klesat zastoupení vysokých frekvencí a obrázek bude stále více zkreslený.

Pseudo-inverzní filtr (digitalizovaný obraz).
digitalizovaný obraz
Pseudo-inverzní filtr (pseudo-inverzní filtr T = 5x10-3).
pseudo-inverzní filtr T = 5x10-3
Pseudo-inverzní filtr (pseudo-inverzní filtr T = 5x10-4).
pseudo-inverzní filtr T = 5x10-4

Obrázek 3.12
Pseudo-inverzní filtr.


Wienerův filtr

          Nejpoužívanější metodou pro restauraci obrazu je Wienerův filtr. Opět vychází s inverzního filtru a omezuje růst vysokých frekvencí, stejně jako pseudo-inverzní filtr, ale činí tak jiným způsobem. Zatímco pseudo-inverzní filtr pokládá koeficienty s velkým nárůstem za nulové, Wienerův filtr pouze omezuje jejich růst. Podobně jako v minulém případě je i zde zaveden parametr - SNR (signal to noise ratio).

          Rovnice 3.28 (3.28)

          Podobně jako v předchozím případě se i zde s klesajícím SNR filtr přibližuje k inverznímu filtru. Při určování SNR lze vyjít z vysokého podílu šumu (kolem 1x10-3) a podle výsledků podíl snižovat. Praktické výsledky použití Wienerova filtru pro různé koeficienty SNR jsou na obrázku 3.13.

Wienerův filtr (digitalizovaný obraz).
digitalizovaný obraz
Wienerův filtr (Wienerův filtr SNR = 1x10-6).
Wienerův filtr SNR = 1x10-6
Wienerův filtr (Wienerův filtr SNR = 1x10-7).
Wienerův filtr SNR = 1x10-7

Obrázek 3.13
Wienerův filtr.


Další informace

          Pro restauraci digitalizovaných map jsem používal převážně Wienerův filtr. Popis inverzního, pseudo-inverzního a Wienerova filtru včetně příkladů lze nalézt v interaktivním kursu [3].


3.2 Úpravy obrazu

          Zatímco v předchozí kapitole o restauraci obrazu jsme se zabývali odstraněním degradací způsobenými digitalizací za předpokladu znalosti vlastností digitalizačního systému, v této kapitole se ještě pokusíme upravit a vylepšit restaurovaný obraz tak, abychom usnadnili další fáze zpracování.


3.2.1 Definice použité oblasti

          Na běžných mapách jsou kromě vlastní mapy vytištěny i další doprovodné informace - hlavička mapy, rozdělení prostoru mezi autory mapy, reklamy apod. Pro účely rozpoznávání je nutné tyto doprovodné informace odstranit, aby zůstala pouze mapa.

          Mapa je obvykle ohraničena buďto nějakými významnými liniovými objekty nebo v případě jejich lokální absence lomenou čarou. Z těchto důvodu k vymezení oblasti s mapou stačí pouze polygon (v nejjednodušších případech to je obdélník). Automatické rozpoznání oblasti s mapou je obtížné a podle mého názoru i zbytečné - ruční definice hranice polygonální oblasti není nijak pracná.

          Pro jednoduchost je definice použité oblasti zahrnuta do předzpracování a tedy je to funkce, která převádí bitmapu na bitmapu. Proto je nutné nahradit oblasti s doprovodnými informacemi barvou pozadí. Barvu pozadí je možné zadat ručně, nebo ji vyhledat jako nejsvětlejší barvu vyskytující se v digitalizovaném obrazu (tato varianta funguje spolehlivě).

          Výsledky ořezání mapy jsou vidět na obrázku 3.14. Pro ořezání jsem použil upravený algoritmus pro vyplňování polygonu uvedený v [1].

Definice použité oblasti (mapa včetně dalších informací).
mapa včetně dalších informací
Definice použité oblasti (ořezaná mapa).
ořezaná mapa

Obrázek 3.14
Definice použité oblasti.


3.2.2 Rozšíření barevné škály

          Rozšíření barevné škály je další z jednoduchých filtrů. Je to vlastně jednoduchá korekce jasu a kontrastu digitalizovaného obrazu. Důvodem zavedení je skutečnost, že úroveň jasu a kontrastu ve výsledku digitalizace je závislá na mnoha faktorech - vlastnostmi papíru počínaje a nastavením scanneru konče. Snahou je "normovat" jas a kontrast tak, aby barvy odpovídaly očekávání.

          Vstupem filtru je obraz a dvě barvy cdark a clight, první z nich určuje nejtmavší barvu (požadovaná černá) a druhá nejsvětlejší barvu (požadovaná bílá) ve výsledném obrazu. Nejprve je potřeba ve vstupním obraze nalézt nejtmavší a nejsvětlejší barvu, ty označím jako cmin a cmax. Pro všechny pixely v obraze provedeme výpočet nové barvy. Novou barvy pixelu cenh spočteme pomocí původní barvy c lineární interpolací:

          Rovnice 3.29 (3.29)

          Výše uvedený postup platí pro monochromatické obrázky. V případě barevných obrázků je nutné postup aplikovat odděleně na vrstvy R, G a B (scanner může zkreslovat v každé barevné rovině jiným způsobem). Výsledky filtru jsou ukázány na obrázku 3.15.

Rozšíření barevné škály (před aplikací filtru).
před aplikací filtru
Rozšíření barevné škály (po aplikaci filtru).
po aplikaci filtru

Obrázek 3.15
Rozšíření barevné škály.

          Na použití tohoto filtru je ale třeba dát pozor. Je nutné si uvědomit, že pracuje s kvantovaným obrazem a výsledek je tedy zatížen chybou zaokrouhlování (při opakovaném použití se chyba může zvětšovat). Možnost alternativní úpravy jasu a kontrastu by měly umožňovat i filtry, které pracují s reálnými čísly (např. filtry založené na FT), úpravu by měly provést ještě před konverzí na celá čísla. Dále je třeba úpravu barevné škály používat v souladu s pravidly pro klasifikaci barev - viz následující kapitola. Popis úpravy jasu a kontrastu obrazu je popsán prakticky každé knize o digitálním zpracování obrazu (např. [2], [3] nebo [4]).


3.2.3 Konvoluční filtr

          Konvoluční (často nazývané také jako lineární) filtry se ve zpracování obrazu hojně využívají. Slouží především k odstraňování šumu, zvýrazňování a vyhledávání hran apod. Při předzpracování map je výhodné využít zejména filtry pro odstranění šumu a zvýraznění hran. Odstranění šumu pomocí konvolučního filtru má ale za následek rozmazání obrazu, což je nepřípustné. Pro odstranění šumu proto doporučují použít pořadové filtry.

          Výpočet filtru vychází z definice diskrétní polohově invariantní konvoluce upravené pro konečnou masku h velikosti 2N+1 x 2N+1:

          Rovnice 3.30 (3.30)

          Prakticky se používají filtry s velikostí masky 3 x 3 - 7 x 7 (N = 1, 2, 3). V následující tabulce jsou uvedeny některé nejběžnější masky velikosti 3 x 3 použitelné pro předzpracování map:

název masky maska výsledek poznámka
průměrování
(mean filter)
Tento filtr sice nejlépe odstraňuje šum, ale na druhou stranu nejvíce rozmazává přechody a detaily. Z důvodů uvedených v kapitole o rozkladu na tiskové barvy jej nedoporučuji používat.
odstranění šumu
(smooth filter, lowpass filter)
Stejně jako předchozí filtr i tento poměrně dobře odstraňuje šum, jeho výhodou oproti průměrování je vyšší míra zachování detailů. Nicméně stále komplikuje klasifikaci barev.
zvýraznění hran
(enhancement filter)
Tento filtr slouží pro zvýraznění hran (zaostření obrazu). Doporučuji ho používat ve spojení s restaurací obrazu pro omezení výskytu přechodů mezi barvami.

Tabulka 3.1
Nejpoužívanější masky pro konvoluční filtr.

          Kromě masek velikosti 3 x 3 se běžně používají i masky větší. Konvoluční filtry je třeba používat v souladu s pravidly pro klasifikaci barev. Další informace jsou uvedeny prakticky v každé knize zabývající se digitálním zpracováním obrazu (např. [2], [3] nebo [4]), masky pro konvoluční filtry jsem převzal z knihy [4].


3.2.4 Pořadový filtr

          Poslední skupinou filtrů pro předzpracování jsou filtry založené na uspořádané posloupnosti intenzit pixelů v okolí bodu velikosti 2N+1 x 2N+1. Pro každý bod obrazu sestrojím dvě posloupnosti. Posloupnost S bude obsahovat intenzity bodů z jeho okolí, posloupnost R bude setříděná posloupnost S:

          Rovnice 3.31 (3.31)

          Pořadové filtry definují hodnoty gk,l jako jednu z hodnot nebo kombinaci hodnot z posloupnosti R. Nejběžnější filtry jsou:

název prvek posloupnosti výsledek (pro N = 3) poznámka
eroze
(erosion)
r1 Tento filtr není pro předzpracování map příliš užitečný. Snižuje totiž podíl syté barvy v detailech a tím komplikuje klasifikaci barev.
medián
(median)
r|K/2| Poměrně úspěšný filtr pro odstraňování šumu. Na rozdíl od konvolučních filtrů není výsledek tolik rozmazán.
dilatace
(dilation)
rK Rozšiřování detailů má za následek zvýšení podílu sytých barev a tím usnadňuje klasifikaci barev. Na druhou stranu nesmí dojít k přílišnému rozšíření a slití blízkých objektů (např. rastrů).

Tabulka 3.2
Aplikace pořadového filtru.

          Kromě filtrů uvedených v tabulce je možné samostatně použít libovolný jiný prvek posloupnosti (tím dosáhneme mírnější eroze/dilatace) nebo kombinaci prvků. Běžně se používají průměry vybraných prvků posloupnosti k dosažení lepšího odstranění šumu, je to ale za cenu většího rozmazání.

          U ostatních druhů filtrů implementace přirozeně vyplynula z jejich definice. U tohoto druhu se zmíním o jednoduché a poměrně efektivní implementaci, kterou jsem převzal z [8]. Využívá následujících vlastností:

          Algoritmus používá pole A celých čísel velikosti počtu různých hodnot, které lze uložit do datového typu pro intenzitu (pro 1 byte je to pole velikosti 256), a proměnnou Min určující index minimálního nenulového prvku v poli A. Pole A slouží pro uložení počtu pixelů v okně rozdělených podle intenzit. Algoritmus prochází vstupní obraz po řádcích. Nejprve je nutné inicializovat pole A s proměnnou Min pro první polohu okna na řádce - ty nám umožní snadné určení k-té hodnoty v setříděném poli (budeme procházet pole A od indexu Min tak dlouho, dokud nepřejdeme k pixelů). Při přechodu na další pozici okna na řádce stačí pouze pole A s proměnnou Min aktualizovat (odstranění starých hodnot a přidání nových).

          Další informace jsou uvedeny prakticky v každé knize zabývající se digitálním zpracováním obrazu (např. [2], [3] nebo [4]).