153ZODH / 6. cvičení

← předchozí cvičení • stránky předmětu • další cvičení →

Detekce hran, ostřící operátory, prahování

Osnova

Lokální zvýraznění obrazu, filtry určené pro detekci hran a ostřící operátory.

Seznam příkazů

• r.mfilter
• r.colors
• r.mapcalc
• d.rgb
• g.copy
• g.remove

Detekce hran

Hranu lze definovat jako diskontinuitu obrazové funkce (tj. signifikatní změnu digitálních hodnot). Hrana je tedy oblast v obraze, kde se výrazně nebo skokově mění digitální hodnota. Hrany lze rozdělit na tři základní typy:

• střechová hrana (roof edge) - světlejší liniový útvar na tmavším pozadí (např. polní cesta, betonová silnice)
• příkopová hrana (ditch edge) - tmavší liniový útvar na světlejším pozadí (např. vodní tok)
• stupňová hrana (step edge) - rozhraní mezi tmavějším a světlejším objektem (např. rozhraní pole-les)

Filtry (vysokofrekvenční) určené k detekci hran jsou založeny na principu plovoucího okna. Tím se vyhledávají lokální hrany, v další fázi je nutno rozhodnout, zda jde o skutečné hrany či pouhé vnitřní nespojitosti obrazu. Zde se obvykle používá metoda práhování. Nakonec se provádí ztenčení hran a konečná filtrace.

Sobelův filtr

Tento filtr nejdříve operuje ve směru sever-jih (N-S) a následně ve směru východ-západ (E-W). Definice pro r.mfilter:

TITLE Sobel (3x3)
MATRIX 3
-1 0 1
-2 0 2
-1 0 1
DIVISOR 1
TYPE P

MATRIX 3
1 2 1
0 0 0
-1 -2 -1
DIVISOR 1
TYPE P

 g.region rast=tm1
 r.mfilter input=tm1 output=tm1_sobel filter=sobel.txt
 r.colors -e map=tm1_sobel color=grey
 d.rast map=tm1_sobel

Obr. č.1: Vybraný detail prvního pásma LandSat-TM5, vlevo originální data, vpravo po aplikaci Sobelova filtru

Dále můžeme provést práhování obrazu

• pokud je DH kladná, ulož do výstupní vrstvy hodnotu "1", v opačném případě "žádná data" (null data)

v zápisu pro r.mapcalc vypadá podmínka if takto:

 r.mapcalc 'tm1_hrana=if (tm1_sobel > 0, 1, null())'
 d.rast map=tm1_hrana

Nakonec můžeme provést ztenčení detekovaných hran pomocí modulu r.thin:

 # ztenčení hran
 #
 r.thin input=tm1_hrana output=tm1_hrana1
 d.rast map=tm1_hrana1

Obr. č.2: Detekované hrany po binarizaci (vlevo) a jejich ztenčení (vpravo)

 # barevná syntéza ve skutečných barvách
 #
 d.rgb red=tm3 green=tm2 blue=tm1
 #
 # zobrazení rastrové vrstvy se ztenčenými hranami (pozor na přepínač -o: transparentní mód)
 #
 d.rast -o map=tm1_hrana1

Obr. č.3: Prezentace výsledků Sobelova filtru, na pozadí barevná syntéza ve skutečných barvách

Robertsův gradient

Tento operátor je definován vztahem

$p(i,j)=|f(i,j)-f(i+1,j+1)|+|f(i+1,j)-f(i,j+1)|\,$

v zápisu pro r.mapcalc:

 # aplikace Robertsova grafientu
 #
 r.mapcalc 'tm1_roberts = abs (tm1[0,0]-tm1[1,1]) + abs (tm1[1,0]-tm1[0,1])'
 r.colors map=tm1_roberts color=grey.eq

Prewittův operátor

Definice pro r.mfilter:

TITLE Prewittuv operator (3x3)
MATRIX 3
-1 -1 -1
0 0 0
1 1 1
DIVISOR 1
TYPE P

MATRIX 3
-1 0 1
-1 0 1
-1 0 1
DIVISOR 1
TYPE P

 r.mfilter input=tm1 output=tm1_prewitt filter=prewitt.txt
 r.colors map=tm1_prewitt color=grey.eq

Obr č.4: První kanál LandSat-TM5, Robertsův a Prewittův operátor

Práhování

Práhováním se rozumí převod DH do množiny o malém počtu prvků, zpravidla do množiny {0,1} (v tomto případě mluvíme o binarizaci obrazových dat). V jistém ohledu lze práhování považovat za segmentaci či velmi jednoduchý způsob klasifikace obrazu.

Této metody se často používá k vytvoření masky, kde oblasti s hodnotou "0" (nebo "žádná data") jsou maskovány a do dalších analýz vstupují pouze oblasti s hodnotou masky "1" (či zcela obecně s nenulovou hodnotou).

Nalezení práhu se většinou děje "empericky", jeho automatické určení je poměrně obtížné. Rozlišujeme dva základní typy:

• jednoduchý práh

p(i,j) = 1: f(i,j) >= q
p(i,j) = 0: f(i,j) < q

Příklad - vytvoření masky pevniny:

Čtvrté pásmo LandSat-TM5 (IR) lze použít pro velmi jednoduché určení vodních ploch:

DH <  40 vodní plocha
DH >= 40 pevnina

 g.region rast=tm4
 #
 # vytvoření rastrové mapy vodní plochy na základě jednoduchého práhu
 #
 r.mapcalc 'voda = if(tm4 < 40, 1, null())'
 #
 # nastavení vhodné tabulky barev (1 -> modrá barva)
 #

Vizualizace:

 # zobrazit RGB syntézu ve skutečných barvách a posléze i rastrovou mapu vodních ploch
 #
 d.rgb red=tm3 green=tm2 blue=tm1
 d.rast -o map=voda
 #
 # zobrazit vektorovou vrstvu s hranicemi vodních ploch (kromě Nechranické nádrže)
 # barva hranic červená
 #
 d.vect map=voda type=boundary color=red

Obr č.5: Maskování vodních ploch určených práhováním tm4, barevná syntéza ve skutečných barvách a vektorová vrstva vodních ploch

• vícenásobný práh

p(i,j) = 0: f(i,j) <= q1
p(i,j) = 1: q1 < f(i,j) <= q2
...
p(i,j) = n: qn1 < f(i,j) < qnn

Příklad: Během čtvrtého cvičení byla provedena jednoduchá reklasifikace NDVI podle následujícího reklasifikačního pravidla:

-100 thru 5   = 1 bez vegetace, vodni plochy
 5   thru 35  = 2 plochy s minimalni vegetaci
 35  thru 87  = 3 plochy pokryte vegetaci

Podobnou reklasifikaci lze provést přímo pomocí r.mapcalc a vícenásobného práhu (navíc se obejdeme bez nutnosti vytvářet celočíselnou mapu pro r.reclass).

 # klasifikace NDVI
 #
 r.mapcalc 'r_ndvi1 = eval( \
 a=if(ndvi >= -1.0 && ndvi < 0.05, 1, 0), \
 b=if(ndvi >= 0.05 && ndvi < 0.35, 2, 0), \
 c=if(ndvi >= 0.35 && ndvi < 0.87, 4, 0),a+b+c)'

Poznámka: Během výpočtu se přechodně vytvoří tři rastrové vrstvy 'a', 'b' a 'c', výsledná vrstva vznikne součtem těchto tří vrstev.

Ostření obrazu

Ostřící operátory zvyšují lokální maxima a snižují lokální minima. Vycházejí z první a druhé derivace, kterou lze v diskrétním obraze přibližně provádět pomocí diferencí.

Jako příklad můžeme uvést La Placeův modifikovaný filtr:

TITLE modifikovany La Placeuv (3x3)
MATRIX 3
0 -1 0
-1 4 -1
0 -1 0
DIVISOR 1
TYPE P

 g.region rast=tm1
 r.mfilter input=tm1 output=tm1_mlp filter=mlp.txt 
 #
 # nastavit tabulku barev na vyrovnané odstíny šedi
 #
 r.colors map=tm1_mlp color=grey.eq

Obr č.6: První pásmo LandSat-TM5 - původní data, po aplikaci modifikovaného La Placeova filtru

Alternativně lze použít operátor definovaný vztahem:

$p(i,j)=f(i,j)+f(i,j)-{\frac {1}{9}}{\bigg (}\sum _{k=-1}^{1}\sum _{l=-1}^{1}f(i+k,j+l)\,.\,k_{1}{\bigg )}+k_{2}$

Rovnice pro r.mapcalc (k₁ = 1; k₂ = 2):

 r.mapcalc 'tm1_o2 = tm1[0,0] + tm1[0,0] - \
 (tm1[-1,-1]+tm1[-1,0]+tm1[-1,1] + \
 tm1[0,-1]+tm1[0,0]+tm1[0,1] + \
 tm1[1,-1]+tm1[1,0]+tm1[1,1] * 1)/9 + 2'
 #
 # nastavit tabulku barev
 #
 r.colors map=tm1_o2 color=grey.eq

Obr č.7: První pásmo LandSat-TM5 - původní data, po aplikaci ostřícího operátoru