153ZODH / 9. cvičení

Z GeoWikiCZ

Fourierova transformace

Osnova

Aplikace Fourierovy transformace s důrazem na využití této metody při odstranění periodického šumu či spektrální filtraci obrazu.

Seznam příkazů

Teorie

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 stoupají četnosti výskytu (jsou kolem DC uspořádá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):

 # nastavit výpočetní region
 #
 g.region rast=tm1
 #
 # Fourierova transformace
 #
 i.fft input_image=tm1 real_image=tm1_r imaginary_image=tm1_i

 # nastavit region a zobrazit výsledné obrazové vrstvy 
 #
 g.region rast=tm1_r
 d.erase -e
 d.rast map=tm1_r
 d.rast map=tm1_i
Obr č.2: Reálná část Fourierova spektra

Filtrace šumu pomocí FFT

Pro identifikaci šumu je použita reálná část nebo doplňková obrazová vrstva obsahující vypočtenou amplitudu (vytvořená pouze pro účel vizualizace):

 # pomocná vrstva s vypočtenou amplitudou
 #
 r.mapcalc 'tm1_amp = sqrt (tm1_r^2 + tm1_i^2)'
 d.rast map=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 formátu (v.digit či pomocí digitalizačního modulu wxGUI) a rasterizace masky pomocí v.to.rast. Následuje postup pro wxGUI.

Do seznamu vrstev přidáme rastrovou vrstvu 'tm1_amp', která bude použita jako podklad pro vektorizaci pruhů. Digitalizační modul wxGUI spustíme z mapového okna "Tools->Digitize". Novou vektorovou vrstvu vytvoříme zvolením "Select vector map->Create new vector map". Vrstvu nazveme 'filtr' a vytvoříme ji bez atributové tabulky. V nastavení digitalizačního nástroje zrušíme volbu "Attributes->Digitize new feature->Add new record into table". Dále nastavíme výchozí kategorii na "1" ("Attributes->Digitize new feature->Category mode" nastavit na "Manual entry"). Nyní přikročíme k vlastní digitalizaci.

Nejprve si zvětšíme detail pruhu, který chceme vektorizovat, z nástrojové lišty digitalizačního nástroje zvolíme "Digitize new boundary". Plochu digitalizujeme (levé tlačítko myši pro výběr dalšího lomové bodu, pravé tlačítko myši pro ukončení digitalizace). Tímto způsobem zvektorizujeme všechny identifikovatelné pruhy (ortogonální na směr letu nosiče). Poté každou plochu označíme centroidem nezbytně nutným pro rasterizaci této vektorové vrstvy. Editaci vektorové vrstvy ukončíme z kontextového menu vrstvy "Stop editing".

Obr č.4: Digitalizace pruhů v prostředí wxGUI

Provedeme rasterizaci vektorové vrstvy 'filtr' pomocí modulu v.to.rast.

 # nutno přenastavit výpočetní region
 #
 g.region rast=tm1_amp
 #
 # rasterizace na základě kategorií ploch
 #
 v.to.rast input=filtr output=filtr use=cat type=area
Obr č.5: Příklad filtru použitého jako masky pro inverzní Fourierovu transformaci

Masku vytvoříme jako doplněk k rastrové vrstvě pomocí r.mask:

 # vytvoření masky
 #
 r.mask -i input=filtr maskcats=1
Obr č.6: Vizualizace tm1_amp s aktivovanou maskou

Nyní aplikuje inverzní Fourierovu transformaci.

 # inverzní Fourierova transformace
 #
 i.ifft real_image=tm1_r imaginary_image=tm1_i output_image=tm1_fft
 #
 # odstranění (deaktivace) masky
 #
 r.mask -r
 #
 # nastavení tabulky barev na "grey.eq"
 #
 g.region rast=tm1_fft
 r.colors map=tm1_fft color=grey.eq
 d.rast map=tm1_fft
Obr č.7: První pásmo LandSat5-TM po odstranění periodického šumu

Spekrální filtrace pomocí FFT

Vrátíme se k problematice filtrů s nízkou a vysokou propustností (více 5. cvičení). 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.

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í min a max (jde-li o prstenec). Přepínač -b vytvoří navíc binární vrstvu (hodnoty 0,1).

Je tedy nejprve nutno určit souřadnice DC bodu a posléze poloměr(y) kruhu, resp. prstence. Příklad pro high pass filtr:

 # zjištění středu a poloměru kruhu
 #
 g.region rast=tm1_amp
 d.rast map=tm1_amp
 #
 # vytvoření rastrové vrstvy obsahující definovaný kruh
 #
 r.circle -b output=ff coordinate=-769058,-988168 max=1550

Masku vytvoříme jako doplněk k rastrové vrstvě 'ff':

 # aktivace masky
 #
 r.mask -i input=ff maskcats=1
 #
 # zobrazení vrstvy s aktivovanou maskou
 #
 d.rast map=tm1_amp
 #
 # aplikace inverzní Fourierovy transformace
 #
 i.ifft real_image=tm1_r imaginary_image=tm1_i output_image=tm1_hp
 #
 # deaktivace masky, nastavení tabulky barev na "grey.eq"
 #
 r.mask -r
 g.region rast=tm1_hp
 r.colors map=tm1_hp color=grey.eq
 d.rast map=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)