153YZOD Zpracování obrazových dat 2006 - 9. cvičení
Fourierova transformace
osnona
V tomto kurzu si ukážeme aplikaci Fourierovy transformace v GRASSu s důrazem na praktické využití této metody – odstranění periodického šumu či spektrální filtraci obrazu.
seznam příkazů
- g.region
- i.fft
- d.erase
- d.rast
- r.mapcalc
- r.digit
- v.digit
- v.build
- v.to.rast
- i.ifft
- g.remove
- r.colors
- d.where
- d.measure
- r.circle
teorie
Na rozdíl od doposud zmíněných metod, které operují čistě v prostorovém souřadnicovém systému, Fourierova transfomace je založena na převodu dat z geometrického znázornění (tj. prostorového souřadnicového systému) do domény frekvenční (souřadnicového systému frekvenčního) – komplexního znázornění složkových četností výskytu. Během této transformace je tak obraz rozložen na frekvenční komponenty (diskrétními DH je proložena spojitá funkce, jde o funkce sinus a kosinus s různými amplitudami, frekvencemi), které jsou uloženy jako komplexní čísla. V případě dvoudimenzionální Fourierovy transformace jsou výsledkem dvě obrazové vrstvy – reálná a imaginarní část.
Tyto četnosti výskytu jsou rozloženy podél obou os a lze je znázornit ve dvourozměrném rozptylogramu (tzv. Fourierově spektru). Nulový bod ("DC" – direct current) má četnost 0,0 (viz obr. č.1). Se vzrůstající vzdáleností od nulového bodu stoupájí četnosti výskytu (jsou kolem DC uspořány symetricky). Jevy, které mají v původním obraze horizontální trend se vyskytují ve vertikálních komponentách Fourierova spektra a naopak.
Tuto metodu lze využít např. pro identifikaci a eliminaci periodického šumu obsaženého v družicovém snímku. V tomto případě je nejprve provedena Fourierova transformace, posléze aplikována maska (šum se zobrazí ve Fourierově spektru ve tvaru pruhů) a data zpětně transformována do prostorové domény (inverzní Fourierova transformace).
FFT
Aplikace Fourierovy transformace (Fourierova transformace v diskrétním tvaru – Fast Fourier Transformation – FFT) je sama o sobě v GRASSu velmi jednoduchá:
#nastavit aktuální region # GRASS > g.region rast=tm1 # #Fourierova transformace # GRASS > i.fft in=tm1 real=tm1_r imag=tm1_i
Následně můžeme vrstvy (reálnou tm_r a imaginární část tm_i) vzniklé po Fourierově transformaci zobrazit v GRASS monitoru:
#nastavit region a zobrazit výsledné obrazové vrstvy # GRASS > g.region rast=tm1_r GRASS > d.erase GRASS > d.rast tm1_r GRASS > d.rast tm1_i
Za poznámku stojí fakt, že data po Fourierově transformaci nejsou uložena ve standardním rastrovém formátu GRASSu. Navíc lze tato data zpětně transformovat do prostorové domény pouze za předpokladu, že nebyla modifikována žádným z modulů GRASSu.
filtrace šumu pomocí FFT
Pro identifikaci šumu je použita reálná čast nebo doplňková obrazová vrstva obsahující vypočtenou amplitudu (vytvořená pouze pro účel vizualizace):
#pomocná vrstva s vypočtenou amplitudou # GRASS > r.mapcalc 'tm1_amp = sqrt (tm1_r^2 + tm1_i^2)' GRASS > d.rast tm1_amp
Periodický šum je identifikovatelný jako pruhy ortogonální vůči směru poruch v originálních datech (často totožného se směrem letu nosiče skeneru). Tyto pruhy se tedy maskují a tím pádem neovlivní výsledek inverzní Fourierovy transformace. Masku (v podstatě binární filtr) lze vytvořit přímo pomocí r.digit (zásadní nevýhodou modulu je nemožnost zoomování či posunu během vytváření nové rastrové vrstvy) či nepřímo – vytvoření masky ve vektorovém tvaru (v.digit, oproti r.digit je nesrovnatelně konformnější) a posléze rasterizace masky pomocí v.to.rast.
postup s r.digit
Po spuštění modulu
GRASS > r.digit
se objeví v konzoli menu
Please choose one of the following A define an area C define a circle L define a line Q quit (and create map)
Zvolíme "A" (tj. napíšeme do konzole a stikneme ENTER)
Buttons: Left: where am i Middle: mark point Right: done
a pomocí myši (střední tlačítko) obtáhneme vybraný pruh, tím určíme hranice nově vytvářené plochy, plochu uzavřeme pravým tlačítkem myši. Následně budeme vyzváni k zadání kategorie plochy, zadáme "1", textový popisek můžeme vynechat, či napsat např. "aktivni". Na otázku
Look ok? (y/n)
odpovíme "y". Tímto způsobem označíme všechny viditelné pruhy (detekovaný periodický šum). Práci s modulem ukončíme "Q" a zadáme název výsledné rastrové vrstvy (např. filtr).
postup s v.digit
Při spuštění modulu vytvoříme novou vektorovou vrstvu a na pozadí vykreslíme rastrovou vrstvu s vypočtenou amplitudou.
GRASS > v.digit -n filtr bg="d.rast tm1_amp"
Modul se ovládá pomocí speciálního menu (viz obr. č.4), nás budou zajímat funkce "digitize new boundary" (1), "digitize new centroid" (2), "zoom in/out" (4) a "exit" (6).
Nejprve si zvětšíme detail pruhu, který chceme vektorizovat, klikneme na tlačítko "digitize new boundary". V menu je nutno změnit mode na "Manual entry" (všechny hranice ploch obdrží kategorii "1") a vypnout "Insert new record to table" (nebude zakládat žádnou atributovou tabulku). Poté plochu obtáhneme, sledujeme pouze aktuální funkci tlačítek myši v menu. Plocha se pravděpodobně neuzavře, v tomto případě si zvětšíme dané místo a opěrný bod (vertex) hranice přesuneme pomocí funkce "Move vertex" (3). Barevné nastavení modulu lze změnit v "Open settings" (5). Takto zvektorizujeme všechny identifikovatelné pruhy (ortogonální na směr letu nosiče). Poté každou plochu označíme centroidem (tzv. štítkovým bodem) nezbytně nutným pro chystanou rasterizaci této vektorové vrstvy. Mode opět zvolíme "Manual entry", číslo kategorie změníme na "1" a vypneme "Insert new record to table". Práci s modulem ukončíme "Exit".
Sestavíme-li si topologii vektorové mapy, počet ploch (number of areas) by měl souhlasit s počtem centroidů (number of centroids).
#sestavení topologie vektorové vrstvy # GRASS > v.build filtr Building topology ... 45 primitives registered Building areas: 100% 8 areas built 8 isles built Attaching islands: 100% Attaching centroids: 100% Topology was built. Number of nodes : 45 Number of primitives: 45 Number of points : 0 Number of lines : 0 Number of boundaries: 37 Number of centroids : 8 Number of areas : 8 Number of isles : 8
Nyní provedeme rasterizaci této vektorové vrsvy, k tomu nám poslouží modul v.to.rast.
#nutno přenastavit aktivní region # GRASS > g.region rast=tm1_amp # #rasterizace na základě kategorií ploch # GRASS > v.to.rast in=filtr out=filtr use=cat
Masku vytvoříme jako doplněk k rastrové vrstvě filtr pomocí r.mapcalc:
#vytvoření masky # GRASS > r.mapcalc 'MASK=if(isnull(filtr),1,null())'
Nyní aplikuje inverzní Fourierovu transformaci, výslednou vrstvu nazveme např. tm1_fft:
#inverzní Fourierova transformace # GRASS > i.ifft real=tm1_r imag=tm1_i output=tm1_fft # #odstranění (deaktivace) masky # GRASS > g.remove rast=MASK # #nastavení tabulky barev na "grey.eq" # GRASS > g.region rast=tm1 GRASS > r.colors map=tm1_fft color=grey.eq GRASS > d.rast tm1_fft
spekrální filtrace pomocí FFT
V jednom z předchozích kurzů jsme již definovali filtry s nízkou a vysokou propustností, nyní se k této problematice vrátíme. Rozlišujeme tři základní typy frekvenčních filtrů: low pass, high pass a band pass (viz tab č.1 a obr. č.8). Tyto filtry lze velmi dobře použít pro modifikaci texturálních struktur či zvýraznění obrazu jako takového, zlepšení kontrastu obrazu. Geometrický tvar těchto filtrů je na rozdíl od filtru pro detekci šumu kruhového tvaru.
| filter | popis | aplikace |
|---|---|---|
| high pass | odstraňuje nízké frekvence | detekce hran |
| low pass | odstraňuje vysoké frekvence | zvýraznění ploch |
| band pass | odstraňuje nízké a vysoké frekvence | ostření malých struktur |
| cut pass | odstraňuje střední frekvence | rozostření malých struktur |
Kruh či prstenec lze velmi jednoduše vytvořit pomocí modulu r.circle. Střed kruhu se určí pomocí parametru coordinate, poloměr pomocí min a max (to podle toho zda jde o kruh či prstenec). Přepínač -b vytvoří navíc binární vrstvu.
Je tedy nejprve nutno určit souřadnice DC bodu (pomocí modulu d.where) a posléze poloměr(y) kruhu, resp. prstence. Příklad pro high pass filtr:
#zjištění středu a poloměru kruhu # GRASS > g.region rast=tm1_amp GRASS > d.rast tm1_amp # #přepínač -1 má za následek ukončení modulu po výběru prvního bodu #zaroveň výstup modulu přesměrujeme do textového souboru # GRASS > d.where -1 > dc-xy.txt # #aktuální funkce tlačítek myši je uvedena v konzoli # GRASS > d.measure GRASS > g.region rast=tm1_amp # #vytvoření rastrové vrstvy obsahující definovaný kruh # GRASS > r.circle -b out=ff coord=-769058,-988168 max=1550
Alternativně lze využít souřadnice uložené v textovém souboru a poněkud si zkomplikovat zápis, každopadně tím získáme obdiv našich přátel;-)
#alternativní zápis
#
GRASS > r.circle -b out=ff coord=`cat dc-xy.txt | awk '{OFS=","; print $1,$2}'` max=1550
Masku vytvoříme jako doplněk k rastrové vrstvě ff:
#aktivace masky # GRASS > r.mapcalc 'MASK=if(isnull(ff),1,null())' # #zobrazení vrstvy s aktivovanou maskou # GRASS > d.rast tm1_amp # #aplikace inverzní Fourierovy transformace # GRASS > i.ifft real=tm1_r imag=tm1_i out=tm1_hp # #deaktivace masky, nastavení tabulky barev na "grey.eq" # GRASS > g.remove rast=MASK GRASS > g.region rast=tm1 GRASS > r.colors map=tm1_hp col=grey.eq GRASS > d.rast tm1_hp
