153ZODH / 9. cvičení
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.
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
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
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".
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
Masku vytvoříme jako doplněk k rastrové vrstvě pomocí r.mask:
# vytvoření masky
#
r.mask -i input=filtr maskcats=1
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
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.
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 (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