Transcript Radiodiagnostika 5

Zpracování digitálního obrazu
Konvoluce, dekonvoluce, Wienerův filtr, Fourierova řada a Fourierova
transformace funkce, derivace obrazu – detekce a zvýraznění hran,
k-prostor.
Mgr. David Zoul
Fakulta 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.
z x    f  g x    f t g x  t dt
R
z x, y    f  g x, y     f t , s g x  t , y  s dtds
RR
z i, j    f  g i, j    f u, v g i  u, j  v 
A
B
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)
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í



f g  g f
f  g  h    f  g   h
f  g  h    f  g    f  h 
a f  g   af   g  f  ag 
Je asociativní se skalárním součinem
Tedy je bilineárním zobrazením z L x L  L
Diracova funkce: f t   t   f t 
f t   t  T   f t  T 
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á:
J ( x, y)   f  h( x, y) 


f ( x   , y   ) h(  ,  )




   
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 konvoluce
V 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í
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í
originál
3x3 průměr
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.
1 1 1
1
1 1 1
9 1 1 1


Convolution Examples: 33 Blur
1 1 1 1 1
1 1 1 1 1
1 
1 1 1 1 1

25 1 1 1 1 1


1 1 1 1 1
Convolution Examples: 55 Blur
1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1


1 1 1 1 1 1 1 1 1
1 
1 1 1 1 1 1 1 1 1
81 1 1 1 1 1 1 1 1 1


1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1


Convolution Examples: 99 Blur
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
289 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Convolution Examples: 1717 Blur
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
HranaVý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
G  f   G  f
Často nepočítame Zerocrossing points ale
maximální hodnotu po
provedení filtrace.
Laplaceův operátor ∆
Lze ho také použít pro zvýraznění hrany.
Provedeme konvoluci s filtrem
0 1 0


  1  4 1
0 1 0


Laplaceův operátor
Kirshův operátor
 5 3 3 
C   5 0 3 
 5 3 3 


Prewittové operátor
 1 0  1


C   1 0  1
 1 0  1


(Tyto filtry jsou pro svislé hrany - detekce v ose x bude
dána transponovanou maticí)
Sobelův operátor
Pro svislé hrany konvoluční filtr vypadá
následovně
 1 0 1 


C    2 0 2
 1 0 1 


Dává větší váhu středu, čímž by mělo docházet
k lepší lokalizaci hran.
Robinsonův operátor
Jako konvoluční filtr pro detekci svislé hrany se
používá
  1 1 1


C    1  2 1
  1 1 1


Convolution Examples: Original Images
  1
 2
  1
 
Convolution Examples: Vertical
Difference
 1
2  1
Convolution Examples: Horizontal
Difference
Detekce hran
Metody pro detekci hrany jsou většinou velmi
citlivé na šum. Proto je rozumné obraz vyhladit a
aplikovat filtr na šum:
1
1
C1   1
9 
1
 1
C2   2
 1

1 1
1 1
1 1
0 1
0 2 
0 1 
Konečný filtr bude vypadat takto
C3  C1  C2  ?
Všesměrovou detekci realizujeme nezávisle v 8
směrech a výsledky spojíme dohromady.
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
a
f t   0 
2


n 1

an cos  nt  

bn sin  nt 
n 1
kde an, bn jsou funkce tvořící tzv. spektrum operátoru f.
Fourierova řada
Nejprve 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í
sin  mT 
a prointegrujeme přes interval délky
T
Dostaneme rovnici
T

2

T
f  t  sin  mt  dt 
0

1
a0 sin  mt  dt 
2
0
T



b sin  mt sin  nt  dt 
n
n1
0

T
a sin  mt cos  nt  dt.
n
n1
0
Fourierova řada
Využitím vzájemné ortogonality funkcí 1, sin, cos dostaneme
T

f  t  sin  mt  dt 
bmT
2
0
Podobně postupujeme při určení koeficientu an čímž získáme vztahy
T
am 
2
T
 f t  cos  mt  dt,
0
T
bm 
2
T
 f t sin  mt  dt.
0
Pro praktické počítání obvykle vyjadřujeme Fourierovu řadu na intervalu ve
m
tvaru
Fm  x  
a0

2

k 1
kde
2 kx
2 kx 

 bk sin
 ak cos

ba
ba 

b
ak 
2
ba

f  x  cos
2 kx
dx,
ba
a
b
bk 
2
ba

a
f  x  sin
2 kx
dx.
ba
Fourierova řada
1 sine
2 sines
4 sines
8 sines
16 sines
32 sines
Gibbsův jev
Fourierova řada
Pří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 transformace
Vý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
T 2
f   
1
T

f  t  dt 
T 2
2
T

T 2
  f t  cos  n   t  dt,
n 1 T 2
Budeme-li uvažovat pouze funkce absolutně integrovatelné na celé reálné ose

 f t  dt  

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
n  ,
   
dostaneme
1


T 2


f  t  cos    t   dt 


 T 2


n1
2
T
Fourierova transformace
Výraz sumace vyjadřuje v limitě T  ∞ integrální součet a poslední rovnice
přejde ve dvojný integrál
 
f   
1

 f t  cos    t  dtd
0 
Dosadíme-li sem podle Eulerova vzorce za funkci cos, dostaneme konečný
výraz pro Fourierův integrál
 
f   
1
2

f t  e
i  t 
dtd 
 
Tento vztah se dá zapsat v symetrickém tvaru jako

1
 1
f   
2  2
 



 f t  e
 it


dt ei d 


Výraz uvnitř hranaté závorky považujeme za Fourierovu transformaci funkce a
zbylá část vztahu udává inverzní Fourierovu transformaci

F  f  t    F   
F
1
1
2
 F     f  t  

f  t  e  it dt

1
2



ei d 
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:
 
F   ,    F  f  x, y   
1
2

f  x, y  ei x  y  dxdy,
 
 
f  x, y   F
1
 f  t    21

 
F  ,   ei x  y  d  d  .
k - prostor
V 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
F k  

f  r  e2 ikr dr
V
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.
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í obrazu
Jestliž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