153ZODH / 9. cvičení: Porovnání verzí

Z GeoWikiCZ
m 155zddp
 
(Není zobrazeno 28 mezilehlých verzí od stejného uživatele.)
Řádek 1: Řádek 1:
<div style="font-size: 13pt; font-weight: bold">Fourierova transformace</div>
{{Zastaralé|155ZDDP}}
----
{{Cvičení|153ZODH|9|Fourierova transformace}}
== Osnova ==


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


< [[153YZOD Zpracování obrazových dat|Stránky předmětu]] {{bullet}} [[153YZOD Zpracování obrazových dat - cvičení 8|Předchozí cvičení]] {{bullet}} [[153YZOD Zpracování obrazových dat - cvičení 10|Další cvičení]]
=== Seznam příkazů ===


__TOC__
* {{GrassPrikaz|g.region}}
== Osnova ==
* {{GrassPrikaz|i.fft}}
* {{GrassPrikaz|d.rast}}
* {{GrassPrikaz|r.digit}}
* {{GrassPrikaz|v.digit}}
* {{GrassPrikaz|v.to.rast}}
* {{GrassPrikaz|r.mask}}
* {{GrassPrikaz|i.ifft}}
* {{GrassPrikaz|r.colors}}
* {{GrassPrikaz|r.circle}}


== Teorie ==  
== Teorie ==  
Řádek 23: Řádek 32:
Aplikace Fourierovy transformace (Fourierova transformace v diskrétním tvaru - ''Fast Fourier Transformation'' - FFT):
Aplikace Fourierovy transformace (Fourierova transformace v diskrétním tvaru - ''Fast Fourier Transformation'' - FFT):


<source lang="bash">
  # nastavit výpočetní region
  # nastavit výpočetní region
  #
  #
Řádek 29: Řádek 39:
  # Fourierova transformace
  # Fourierova transformace
  #
  #
  i.fft input=tm1 real=tm1_r imag=tm1_i
  i.fft input_image=tm1 real_image=tm1_r imaginary_image=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  
  # nastavit region a zobrazit výsledné obrazové vrstvy  
Řádek 37: Řádek 45:
  g.region rast=tm1_r
  g.region rast=tm1_r
  d.erase -e
  d.erase -e
  d.rast tm1_r
  d.rast map=tm1_r
  d.rast tm1_i
  d.rast map=tm1_i
</source>


[[Soubor:ZOD-cv8-tm1_r.png|frame|center|Obr č.2: Reálná část Fourierova spektra]]
[[Soubor:ZOD-cv8-tm1_r.png|frame|center|Obr č.2: Reálná část Fourierova spektra]]
Řádek 46: Řádek 55:
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):
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):


<source lang="bash">
  # pomocná vrstva s vypočtenou amplitudou
  # pomocná vrstva s vypočtenou amplitudou
  #
  #
  r.mapcalc 'tm1_amp = sqrt (tm1_r^2 + tm1_i^2)'
  r.mapcalc 'tm1_amp = sqrt (tm1_r^2 + tm1_i^2)'
  d.rast tm1_amp
  d.rast map=tm1_amp
</source>


[[Soubor:ZOD-cv8-tm1_amp.png|center|frame|Obr č.3: Obrazová vrstva s&nbsp;vypočtenou amplitudou]]
[[Soubor:ZOD-cv8-tm1_amp.png|center|frame|Obr č.3: Obrazová vrstva s&nbsp;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í {{GrassPrikaz|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 ({{GrassPrikaz|v.digit}} či pomocí digitalizačního modulu {{GrassPrikaz|wxGUI}} a rasterizace masky pomocí {{GrassPrikaz|v.to.rast}}. Následuje postup pro wxGUI.
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í {{GrassPrikaz|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 ({{GrassPrikaz|v.digit}} či pomocí digitalizačního modulu {{GrassPrikaz|wxGUI}}) a rasterizace masky pomocí {{GrassPrikaz|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  
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  
Řádek 60: Řádek 71:
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".
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' ({{GRASSPrikaz|v.to.rast}}).
[[Soubor:ZOD-cv8-wxdigit.png|thumb|640px|center|Obr č.4: Digitalizace pruhů v prostředí wxGUI]]


Provedeme rasterizaci vektorové vrstvy 'filtr' pomocí modulu {{GrassPrikaz|v.to.rast}}.
<source lang="bash">
  # nutno přenastavit výpočetní region
  # nutno přenastavit výpočetní region
  #
  #
Řádek 68: Řádek 82:
  # rasterizace na základě kategorií ploch
  # rasterizace na základě kategorií ploch
  #
  #
  v.to.rast in=filtr out=filtr use=cat
  v.to.rast input=filtr output=filtr use=cat type=area
</source>


[[Soubor:ZOD-cv8-filtr.png|frame|center|Obr č.5: Příklad filtru použitého jako masky pro inverzní Fourierovu transformaci]]
[[Soubor:ZOD-cv8-filtr.png|frame|center|Obr č.5: Příklad filtru použitého jako masky pro inverzní Fourierovu transformaci]]


Masku vytvoříme jako doplněk k rastrové vrstvě  pomocí <tt>r.mapcalc</tt>:
Masku vytvoříme jako doplněk k rastrové vrstvě  pomocí {{GrassPrikaz|r.mask}}:


<source lang="bash">
  # vytvoření masky
  # vytvoření masky
  #
  #
  r.mapcalc 'MASK=if(isnull(filtr),1,null())'
  r.mask -i input=filtr maskcats=1
</source>


[[Soubor:ZOD-cv8-tm1_amp_maska.png|frame|center|Obr č.6: Vizualizace tm1_amp s aktivovanou maskou]]
[[Soubor:ZOD-cv8-tm1_amp_maska.png|frame|center|Obr č.6: Vizualizace tm1_amp s aktivovanou maskou]]
Řádek 82: Řádek 99:
Nyní aplikuje inverzní Fourierovu transformaci.
Nyní aplikuje inverzní Fourierovu transformaci.


<source lang="bash">
  # inverzní Fourierova transformace
  # inverzní Fourierova transformace
  #
  #
  i.ifft real=tm1_r imag=tm1_i output=tm1_fft
  i.ifft real_image=tm1_r imaginary_image=tm1_i output_image=tm1_fft
  #
  #
  # odstranění (deaktivace) masky
  # odstranění (deaktivace) masky
  #
  #
  g.remove rast=MASK
  r.mask -r
  #
  #
  # nastavení tabulky barev na "grey.eq"
  # nastavení tabulky barev na "grey.eq"
  #
  #
  g.region rast=tm1
  g.region rast=tm1_fft
  GRASS&nbsp;>&nbsp;r.colors map=tm1_fft color=grey.eq
  r.colors map=tm1_fft color=grey.eq
  d.rast tm1_fft
  d.rast map=tm1_fft
</source>


[[Soubor:ZOD-cv8-tm1_sum.png|frame|center|Obr č.7: První pásmo LandSat5-TM po odstranění periodického šumu]]
[[Soubor:ZOD-cv8-tm1_sum.png|frame|center|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 {{153YZODCv|5|Filtrace_obrazu|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.
<center>
{|class="border"
|+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
|-
|}
</center>
[[Soubor:ZOD-cv8-x_pass.png|frame|center|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 {{GrassPrikaz|r.circle}}. Střed kruhu se určí pomocí parametru <tt>coordinate</tt>, poloměr pomocí <tt>min</tt> a <tt>max</tt> (jde-li o prstenec). Přepínač <tt>-b</tt> 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:
<source lang="bash">
# 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
</source>
Masku vytvoříme jako doplněk k rastrové vrstvě 'ff':
<source lang="bash">
# 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
</source>
[[Soubor:ZOD-cv8-hp_postup.png|center|frame|Obr č.9: Postup při tvorbě high pass filtru a aktivace masky]]


----
[[Soubor:ZOD-cv8-tm1_lhp.png|center|frame|Obr č.10: Zvýrazněné první pásmo LandSat5-TM (low & high pass filtr)]]
< [[153YZOD Zpracování obrazových dat|Stránky předmětu]] {{bullet}} [[153YZOD Zpracování obrazových dat - cvičení 8|Předchozí cvičení]] {{bullet}} [[153YZOD Zpracování obrazových dat - cvičení 10|Další cvičení]]


{{ZOD}}
{{ZOD}}
{{GRASS}}

Aktuální verze z 5. 9. 2014, 18:48

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)