153YZOD Zpracování obrazových dat 2006 - 9. cvičení

Z GeoWikiCZ

[ Zpracování obrazových dat ]

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ů

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.

Obr č.1: Rozdělení frekvencí ve Fourierově spektru

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.

Obr č.2: Reálná část Fourierova spektra

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
Obr č.3: Obrazová vrstva s vypočtenou amplitudou

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).

Obr č.4: Menu modulu v.digit

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

Obr č.5: Příklad filtru použitého jako masky pro inverzní Fourierovu transformaci

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())'
Obr č.6: Vizualizace tm1_amp s aktivovanou maskou

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
Obr č.7: První pásmo LandSat5-TM po odstranění periodického šumu

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 passband 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.

Tab č.1: Frekvenční filtry
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
Obr č.8: Standardní filtry pro zvýraznění obrazu pomocí Fourierovy transformace. Tmavé plochy nejsou maskovány (hodnota "1"), světlé naopak maskovány jsou (hodnota "0")

Kruh či prstenec lze velmi jednoduše vytvořit pomocí modulu r.circle. Střed kruhu se určí pomocí parametru coordinate, poloměr pomocí minmax (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
Obr č.9: Postup při tvorbě high pass filtru a aktivace masky
Obr č.10: Zvýrazněné první pásmo LandSat5-TM (low & high pass filtr)


[ Zpracování obrazových dat ]