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