Zpracování digitálního obrazuKonvoluce, dekonvoluce, Wienerův filtr, Fourierova řada a Fourierova transformace funkce, derivace obrazu – detekce a zvýraznění hran,
k-prostor.
Mgr. David ZoulFakulta biomedicínského inženýrství ČVUT
2013
Zpracování digitálního obrazu
Zpracování obrazu nepřidá žádnou informaci, která v původním obraze nebyla přítomna – pouze může zvýraznit již obsaženou informaci.
Možnost dodatečné úpravy obrazu je největší výhodou oproti analogovému obrazu, jako je třeba rentgenologický film.
Konvoluce
Pohybujeme se v prostoru L(A) pro danou A omezenou.
Nechť f a g jsou funkce z L(A). Definujeme funkce z následujícím způsobem.
Funkce f představuje původní obrazovou informaci, funkce z zpracovanou (konvolvovanou) obrazovou informaci a funkce g je tzv. jádro konvoluce.
Funkci z říkáme konvoluce funkce f a g. Funkce z bude také patřit k L(A)
vjuigvufjigfjiz
dtdssytxgstfyxgfyxz
dttxgtfxgfxz
A B
R R
R
,,,,
,,,,
Konvoluce
Konvoluce dvou signálů: obdélníkového pulsu a impulsní charakteristiky RC článku. Výsledek je stejný jako odezva RC článku na stejný puls.
Některé vlastnosti
Konvoluce je komutativní
Konvoluce je asociativiní
Konvoluce je distributivní
agfgafgfa
hfgfhgf
hgfhgf
fggf
Je asociativní se skalárním součinem Tedy je bilineárním zobrazením z L x L L Diracova funkce:
tfttf TtfTttf
Konvoluce
Konvoluce je velmi často používaná operace nejen ve zpracování obrazu ale i ve fyzice, dozimetrii, spektrometrii, teorii pravděpodobnosti.
V aplikacích budeme samozřejmě používat její diskrétní verzi.
Při zpracování obrazu bude navíc dvourozměrná:
),(),(),(),( hyxfyxhfyxJ
Proces konvoluce
V praxi se pracuje se čtvercovými maticemi obrazových bodů (pixelů). Přitom hodnoty několika bodů matice původního obrazu ovlivňují hodntotu jediného (prostředního) bodu ve výsledném obrazu – konvoluce.Výsledný obraz se získá zobrazením příslušné čtvercové matice s vybraným prostředním bodem, přes tzv. filtr, čili masku, tvořící jádro konvoluce.
Proces konvoluceV případě diskrétní konvoluce lze jádro chápat jako tabulku (konvoluční maska), kterou položíme na příslušné místo obrazu. Každý pixel překrytý tabulkou vynásobíme koeficientem v příslušné buňce tabulky a provedeme součet všech těchto hodnot. Tím dostaneme jeden nový pixel.
Vyhlazení
origináloriginál 3x3 průměr3x3 průměr
Jedná se o konvoluci s jádrem tvořeným filtrem s tzv. dolní propustností (low pass filter). Všechna čísla v masce stejná a kladná – dochází ke zprůměrování okolních hodnot a tedy redukci šumu, kontrastu, ostrosti i rozlišení
Convolution Examples: Original Images
Nový pixel, který vypočteme po aplikaci na jedno místo v původním obraze, tedy bude průměrem z devíti okolních pixelů. Neudělali jsme totiž nic jiného, než že jsme sečetli hodnoty 9 pixelů a vydělili 9. Pokud aplikujeme konvoluci na celý obraz, pak dostaneme rozostřený obraz. Pokud použijeme větší konvoluční masku 5×5 s koeficienty 1/25, pak bude obraz rozostřen více.
Convolution Examples: 33 Blur
111111111
9
1
Convolution Examples: 55 Blur
1111111111111111111111111
25
1
Convolution Examples: 99 Blur
111111111111111111111111111111111111111111111111111111111111111111111111111111111
811
Convolution Examples: 1717 Blur
1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
289
1
Derivace obrazu
Koeficienty uvnitř konvoluční masky udávají vliv hodnoty pixelu pod nimi. Lze tak nadefinovat velké množství operací, např. derivaci obrazu (u diskrétního obrazu mluvíme o tzv. odhadu derivace), neboli zvýraznění hran.
Pokud hranu definujeme jako velkou změnu jasové funkce, bude v místě hrany velká hodnota derivace jasové funkce. Maximální hodnota derivace bude ve směru kolmo na hranu. Kvůli jednoduššímu výpočtu se ale hrany detekují jen ve dvou, resp. ve čtyřech směrech. Velká skupina metod na detekci hran aproximuje tuto derivaci pomocí konvoluce s vhodným jádrem. Nejjednodušší taková jsou (-1, 0, 1) a (-1, 0, 1)T.
Detekce a zvýraznění hran
HranaVýrazná změna intenzity.
Lidské oko se podle hran významně orientuje.
Plánovací systém detekuje hrany při automatickém konturovaní struktur, nebo automatchingu.
Typy hrany:
Detekce vs. Zvyraznění hran
Detekce Zvýraznění
Zvýraznění hran
Použití filtru s tzv. horní propustností (high pass filter), obsahující v masce kladná i záporná čísla.
Zvýrazní se rozhraní, zvýší se šum, zhorší se rozlišení nízkokontrastních objektů.
Inverzí jádra můžeme vždy provést dekonvoluci, čímž z konvolvovaného obrazu získáme opět původní.
Harmonizovaný obraz získáme odečtením původního a vyhlazeného obrazu – zvýrazní se pouze rozhraní
Zvýraznění hran
Laplaceův operátor ∆
Hledáme body kde druhá derivace je nulová.
Zero-crossing points Přechody mezi kladnou a zápornou hodnotou, v těch místech zaznamenáme hranu.
Marrova Varianta: vyhladíme obraz se širokým Gaussovským filtrem
Často nepočítame Zero-crossing points ale maximální hodnotu po provedení filtrace.
fGfG
Laplaceův operátor ∆Lze ho také použít pro zvýraznění hrany.
Provedeme konvoluci s filtrem
010
141
010
Laplaceův operátor
Kirshův operátor
Prewittové operátor
(Tyto filtry jsou pro svislé hrany - detekce v ose x bude dána transponovanou maticí)
5 3 3
5 0 3
5 3 3
C
101
101
101
C
Sobelův operátor
Pro svislé hrany konvoluční filtr vypadá následovně
Dává větší váhu středu, čímž by mělo docházet k lepší lokalizaci hran.
Jako konvoluční filtr pro detekci svislé hrany se používá
101
202
101
C
Robinsonův operátor
111
121
111
C
Convolution Examples: Original Images
Convolution Examples: Vertical Difference
121
Convolution Examples: Horizontal Difference
121
Detekce hranMetody pro detekci hrany jsou většinou velmi citlivé na šum. Proto je rozumné obraz vyhladit a aplikovat filtr na šum:
Všesměrovou detekci realizujeme nezávisle v 8 směrech a výsledky spojíme dohromady.
1
2
3 1 2
1 1 11
1 1 19
1 1 1
1 0 1
2 0 2
1 0 1
Konečný filtr bude vypadat takto
?
C
C
C C C
Příklad praktického využití – analýza obrazu
Zvýraznění hran použitím derivace obrazu a následného podmíněného formátování
Příklad výstupu – analýza funkce clon lineárního urychlovače automatickým změřením velikosti 3 různých polí na ozářeném
filmu importovaném do Excelu
Příklad 1: Stáhněte si textový soubor s názvem „data“, zobrazte v Excelu jeho obsah s pomocí podmíněného formátování ve stupních šedi. Zviditelněte všechny hrany pomocí konvoluce s Robinsonovým operátorem. Zvýrazněte hrany za pomoci konvoluce s Laplaceovým operátorem.
Praktické provedení ozáření malých polí filmu rostoucí a opět klesající dávkou
Zobrazení gradientů dávky
Jiný příklad – Winston-Lutzův test stereotaktické radioterapie (radiochirurgie) mozku
(průměty tužkového svazku z 5 různých polí, analyzované v Excelu)
Video
Grafy dávkové distrubuce – rovina xz a yz (vlevo) a automatická analýza polohy objektu v kruhovém poli v týchž
rovinách (vpravo)
Grafy dávkové distrubuce – rovina xz a yz (vlevo) a automatická analýza polohy objektu v kruhovém poli v týchž
rovinách (vpravo)
Výsledky analýzy odchylek ozařovaného objektu od centrální osy svazku pro 5 různých úhlů gantry a stolu
Fourierova řada
Jean Baptiste Joseph Fourier (1768 – 1830)
Nejjednodušší odvození Fourierovy transformace vychází z tzv. Fourierovy řady periodické funkce, jejíž motivaci lze nalézt ve skládání anizochronních harmonických kmitů téhož směru s takovými frekvencemi, aby výsledná funkce mohla být periodická, tedy T = nTn, kde n je celé císlo. Funkce daná touto superpozicí bude mít tvar
kde an, bn jsou funkce tvořící tzv. spektrum operátoru f.
0
1 1
cos sin2 n n
n n
af t a n t b n t
Fourierova řadaNejprve budeme uvažovat funkci periodickou na intervalu a budeme předpokládat platnost výše uvedeného rozvoje pro nějakou kombinaci koeficientuů an, bn. Obě strany rovnosti vynásobíme funkcí
a prointegrujeme přes interval délky Dostaneme rovnici
2T
sin m T
0
0 0
1 0
1 0
1sin sin
2
sin sin
sin cos .
T T
T
n
n
T
n
n
f t m t dt a m t dt
b m t n t dt
a m t n t dt
Fourierova řadaVyužitím vzájemné ortogonality funkcí 1, sin, cos dostaneme
Podobně postupujeme při určení koeficientu an čímž získáme vztahy
Pro praktické počítání obvykle vyjadřujeme Fourierovu řadu na intervalu ve tvaru
kde
0
sin2
T
mb Tf t m t dt
0
0
2cos ,
2sin .
T
m
T
m
a f t m t dtT
b f t m t dtT
0
1
2 2cos sin
2
m
m k k
k
a kx kxF x a b
b a b a
2 2cos ,
2 2sin .
b
k
a
b
k
a
kxa f x dx
b a b a
kxb f x dx
b a b a
Fourierova řada
1 sine 2 sines 4 sines
8 sines 16 sines 32 sines
Gibbsův jev
Fourierova řadaPříklad 2: Sestrojte Fourierovu řadu padesátého stupně následujících signálů jednotkové amplitudy: a) Jednotkové obdélníkové pulsy,b) Rovnoramenné pilovité pulsy,c) Cykloida jednotkového poloměru.
Fourierova transformaceVýraz pro Fourierovu transformaci můžeme odvodit z Fourierovy řady provedením limitního procesu , tedy zvolením nekonečné periody, čímž umožníme využití této metody i pro funkce, kterénejsou periodické. Dosadíme-li do Fourierovy řady vzorce pro koeficienty am, bm, pak využitím trigonometrického vztahu cos( - ) = cos cos + sin sin , dostaneme
Budeme-li uvažovat pouze funkce absolutně integrovatelné na celé reálné ose
pak první člen bude mít v limitě pro T ∞ nulovou hodnotu. Ve druhém členu máme aritmetickou posloupnost s konstantní diferencí. Označíme-li
dostaneme
2 2
12 2
1 2cos ,
T T
nT T
f f t dt f t n t dtT T
f t dt
,
2
n
T
2
1 2
1cos
T
n T
f t t dt
Fourierova transformaceVýraz sumace vyjadřuje v limitě T ∞ integrální součet a poslední rovnice přejde ve dvojný integrál
Dosadíme-li sem podle Eulerova vzorce za funkci cos, dostaneme konečný výraz pro Fourierův integrál
Tento vztah se dá zapsat v symetrickém tvaru jako
Výraz uvnitř hranaté závorky považujeme za Fourierovu transformaci funkce a zbylá část vztahu udává inverzní Fourierovu transformaci
0
1cosf f t t dtd
1
2i tf f t e dtd
1 1
2 2i t if f t e dt e d
1
2i tf t F f t e dt
F
1 1
2iF f t e d
F
Vícerozměrné zobecnění
Dvourozměrnou Fourierovu transformaci můžeme definovat v bázi z funkcí exp[−i(kx + ly)] tak, aby zůstaly zachovány vlastnosti platné pro jednoduchou transformaci. Definujeme tedy:
1
1, , , ,
2
1, , .
2
i x y
i x y
F f x y f x y e dxdy
f x y f t F e d d
F
F
k - prostorV prostorové oblasti, obvyklém eukleidovském prostoru (r-prostoru), je obraz zobrazované veličiny f popsán distribuční funkcí, neboli polem, f(x,y,z). Ve vektorovém zápisu, zavedením prostorového vektoru r, je tato funkce F(r). Obecnou Fourierovou transformací vzniká nová distribuční funkce
kde k = (k1, k2, k3) je vlnový vektor.Integruje se přes prostorovou oblast V. Distribuční funkce je definována v novém lineárním 3-rozměrném vektorovém prostoru. Prostorová f(k) i frekvenční distribuční funkce nesou tutéž informaci a souvisejí spolu přímou a inverzní Fourierovou transformací.Z matematického hlediska tedy z běžného metrického eukleidovského r-prostoru Fourierovou transformací vzniká nový "frekvenční" prostor, označovaný někdy jako k-prostor (k-space).
Název vznikl podle toho, že po Fourierově transformaci je novou nezávisle proměnnou "vlnový" vektor k (obecně komplexní). Abstraktní k-prostor je v jistém smyslu "reciproční" k obvyklému fyzikálnímu r-prostoru.
2 i
V
f e d krk r rF
Vlastnosti k – prostoru k-prostor nese úplnou informaci o obrazu zakódovanou ve frekvenční oblasti
Vysoké frekvence jsou zásadní pro kontrast obrazu, chybí však ostrost kontur
Nízké frekvence nesou informaci o konturách, chybí však kontrast
Periodické poškození obrazuJestliže je poškození periodické, bude se jasně projevovat ve Fourierově prostoru.
Transformujeme poškozený obraz pomocí FT a sledujeme symetrické píky mimo střed v k-prostoru.
Odstraníme tyto frekvence a aplikujeme inverzní FT.
x
