Kvalita obrazuKontrast, rozlišení, neostrost, šum a jeho eliminace, digitalizace obrazu

(kvantování a vzorkování), aliasing a antialiasing, ekvalizace histogramu.

Mgr. David ZoulFakulta biomedicínského inženýrství ČVUT

2013

Předzpracování obrazu

Image enhancement methods:

Rekonstrukce obrazu:

• Jas, kontrast.• Ekvalizace histogramu.• Odstranění šumu. • Periodické poškození obrazu. • Detekce a zvýraznění hran a rohů.

• Rozpoznání poškození.• Invertování.• Odstranění šumu.

Kontrast obrazuV levé části obrázku je znázorněno zobrazení jednoduché struktury kruhového tvaru průměru d a optické denzity A, obklopené homogenním prostředím - pozadím - o optické denzitě B. Výraznost struktury vůči pozadí můžeme charakterizovat jako kontrast objektu Cobj = (A-B)/B.

Vlivem nedokonalého prostorového rozlišení se ostré kontury původního objektu A rozostřily a snížil se rozdíl mezi maximem v obraze A* a pozadím B* - kontrast obrazu Cimg = (A*max-B*)/B* je nižší než kontrast objektu Cobj: Cimg<Cobj.

Za předpokladu struktury kruhového tvaru a Gaussovského konvolučního rozmazání (odezvová funkce bodového zdroje PSF kamery má tvar Gaussovy křivky s pološířkou FWHM) je vztah mezi kontrastem objektu a obrazu dán exponenciálním výrazem:

    

kde d je průměr struktury. Pro velké struktury (d > 4 FWHM) se kontrast zobrazením téměř nezmění (Cimg Cobj). U struktur malých rozměrů, srovnatelných nebo menších než FWHM, je však degradace kontrastu velmi výrazná, Cimg Cobj .

2FWHMd

img objC C e

Kontrast obrazuKontrast: Rozdíl v signálu sousedních oblastí výsledného obrazu – závisí na rozdílu zeslabovacích vlastností zobrazovaného objektu: a) rozdíl tlouštěk z

b) rozdíl v µ

1) kontrast objektu2) kontrast detektoru3) kontrast monitoru4) vnímavost oka

Redukce kontrastu:

kde S je množství rozptýleného a P množství primárního záření, dopadajícího na receptor obrazu (žádoucí, aby byl podíl S/P (SPR) co nejmenší).

0

S1P

CC

0 1 zA B

A

I IC eI

Příklad 1: dokažte

platnost této rovnosti

Pomůcky pro kontrolu zobrazení u přístrojů s přímou digitalizací obrazu Měření kontrastu

RozlišeníRozlišení zobrazovacího systému lze popsat prostřednictvím odezvy na bodový impuls.Obraz bodového mpulsu – PSF (Point Spread Function) – má v ideálním případě tvar gaussovského píku, u něhož stanovujeme rozlišení standardně jako FWHM (Full Widh in Half Magnitude).Bodový impuls (Diracovu funkci) realizujeme při klasické RTG projekci tzv. dírkovou kamerou, při CT projekci tenkým kovovým drátkem kolmým k rovině řezu. V reálu má však ohnisko rentgenky, stejně jako otvor v dírkové kameře konečné rozměry, a obrazem bodového impulsu je proto právě Gaussova funkce.V důsledku rozbíhavosti svazku se při klasické RTG projekci tvar PSF v rámci detektoru obrazu mění místo od místa – gaussovskému píku odpovídá přibližně na centrální ose svazku

Rozlišení

Reálný obraz je složen z množství bodových impulsů různé amplitudy.PSF určuje neostrost výsledného obrazu.Původní ostrý obraz lze získat pomocí dekonvoluce, je-li známo jádro konvoluce – PSF.

Rekonstrukce obrazu

I f g n

Poškozený obraz „Skutečný“ obraz Point Spread Function (PSF)Předpokládejme, že je známa

Šum

DekonvoluceVztahu pro neznámou X se často říká konvoluční rovnice.Proces vyřešení této rovnice se nazývá dekonvoluce.Jestliže g není známé resp. je částečně známé, pak se procesu říká krátkozraká čili myopická dekonvoluceKorelace je operace podobná jako konvoluce:

R

I x f g x f t g t x dt

f g n I

Konvoluční Věta

Přímou aplikací této věty je dekonvoluce

gfgfF ˆˆ

1

ˆˆ ˆ ˆˆ ˆˆ

ˆˆ

I f g n

I f g n

I nfg

f F f

Rychlá metoda pomocí FFT Pro nízké frekvence je druhý člen

velký a metoda v přítomnosti šumu selhává.

Fourier quotient method.

gfgfF ˆˆ

Pro korelace platí

Wienerův Filtr

Myšlenka je následující:

ˆˆ ˆ ˆ

ˆChceme najít filtr takový, abyˆ ˆˆ ˆ ˆ ˆ

a pak ˆˆ ˆ ˆ

ˆˆ

I f g n

I f g nw

w f g n f g

w f g nf

g

Wienerův Filtr

Tento přístup je bohužel příliš optimistický. Reálnějším přístupem je:

2

ˆˆ ˆ ˆˆ

ˆˆa hledáme tak aby

ˆ ˆ byl minimální

w f g nf

gw

f f

2

2 2

Takový filtr má tvar:

ˆ ˆˆ

ˆ ˆ ˆ

f gw

f g n

Originál Poškození + šum

Rekonstrukce pomocí Wienera

Periodicke Artefakty

Originál Poškození + šum

Rekonstrukce pomoci Wienera

RozlišeníSchopnost zobrazovacího systému rozlišit blízké kontrastní objekty ve výsledném obrazu

FWHM signálu PSF (point spread function) pořízeného skrze dírkovou kameru 10 µmLSF (line spread function) signálu pořízeného štěrbinovou kamerouESF (edge spread function) odezva na ostrou hranu – používá se pro hodnocení vlivů závislých na velikosti pole (např. SPR)

Gaussova funkceBodový impuls PSF, či průmět snadněji měřitelné LSF, má charakter Gaussovy funkce 

 kde 0 je parametr určující „šířku“ funkce. Stanovení jejího Fourierova obrazu provedeme z definičního vztahu, v němž obě exponenciely sloučíme a exponenty převedeme na součet čtverce a části nezávislé na x:  

2

2x

x e

2222

2 22 2 2 2

22 2

2

2 4 24

4

1 12 2

1 1 2 2

1 .2

xx i xi x

x i x i

x

x e e dx e dx

e dx e e dx

e e dx

F

Gaussova funkcePoslední integrand nemá primitivní funkci, takže jej bylo nutno zintegrovat lebesgueovsky:

čili

Lze tedy psát

což je hledaný Fourierův obraz Gaussovy funkce. Je vidět, že tvar funkce se zachová, ale šířky originálu a obrazu jsou si (v souladu s principem neurčitosti) nepřímo úměrné.

2 2 2 2

2 2 2

2 2

2 2

2

0 0 0 0

22

00

2 2

2 2 ,2

x x y

e dx e dydx e d d

e d e

2

2x

e dx

2 2

4

2x e

F

Princip neurčitosti

FT

prostor frekvenční

FT

prostor frekvenční

Malý objekt v prostoru je velký ve frekvenční oblasti a

na opak.

2

Jestliže je velikost objektu v prostoru a je velikostobjektu ve frekvenční oblasti, pak platí

116

x yu v

x y u v

RozlišeníF (prostorová frekvence) udává se v lp/mm.Měření F lze provádět buď pomocí Siemansovy hvězdice, nebo pomocí čárových testů.

RozlišeníMTF (modulační přenosová funkce) popisuje, jakým způsobem zobrazovací systém zaznamenává objekty se zvyšujícím se F. Při vysokém F dochází k modulaci MTF.MTF lze spočítat Fourierovou transformací LSF. Pokud má LSF gaussovský průběh, rovněž i MTF má gaussovský průběh. Čím je užší LSF (a tedy lepší rozlišení), tím je širší její Fourierova transformace MTF (důsledek principu neurčitosti) a tím vyšší je rozlišitelná F.

Rozlišení

Neostrost obrazuNeostrost:

Geometrická neostrost (polostín)

kde D je nominální rozměr ohniska

Pohybová neostrost

kde

a α je úhel směru pohybu vzhledem k rovině receptoru

Materiálová neostrost – souvisí s použitím zesilujících fólií, ovlivněna jejich tloušťkou, velikostí zrn, kontaktem film-fólie.

2 2 2geom pohyb materialU U U U

1 objekt receptorgeom

ohnisko objekt

D df U D M

d

cospohybU vtM

ohnisko film

objekt film

dIMO d

ŠumAditivní model NoiseSignalSignal CleanM

NoiseSignalSignal CleanM *Modely šumu

Typy šumu

Korelované

• Elektrická interference• Zdroj-Detektor interference• Moiré patterns

• Kvantový šum na CCD kamery• Kvantizovaný šum na dig. fotkách• Neuronal noise in a retina

Nekorelované

Multiplikativní model

Příklad 2: Vypočtěte celkovou neostrost rentgenovaného objektu, který se pohybuje rychlostí 0,5 ms-1, jestliže expoziční čas činil 5 ms, vzdálenost ohnisko – objekt byla 60 cm, SID = 100 cm, a nominální velikost ohniska byla 0,3 mm. Materiálovou neostrost zanedbejte.

Šuma) Kvantový šum (poissonovské rozdělení)

b) Poměr signál – šum

c) Kvantová detekční účinnost

d) Strukturální šum – vzniká důsledkem náhodných fluktuací v počtu absorbovaných fotonů v jednotce plochy zesilující fólie filmu v důsledku nehomogenity scintilační vrstvy

e) Zrnitost filmu – vzniká v důsledku nehomogenit v rozložení např. radiofotoluminiscenčních center na povrchu filmu a jejich odezvy při skenování

f) Poměr kontrast – šum

kde tloušťka vrstvy Al = 0,2 mm

detN

detdet dop 40NSNR N QDE N

det

dop

NQDE

N

2 2

5

2

PV Al PV normCNR

Al norm

Příklad 3: Vypočtěte poměr signál – šum, jestliže na pixel detektoru dopadlo 106 fotonů a kvantová detekční účinnost detektoru je 10%.

Příklad 4: Vypočtěte poměr kontrast – šum pro záření o polotloušťce 2 mm Al a pro stejné parametry, jako v předešlé úloze.

Kontrast a šumNáhodné fluktuace detekovaného signálu, nezávisle pro každý pixel

Relativní statistické fluktuace s/n = 1/(N) jsou tím nižší, čím vyšší je počet impulsů nastřádaný v jednotlivých pixelech obrazu. Konstantní pozadí B je tedy zobrazeno jako plocha, jejíž body kolísají zhruba mezi B*±(B*), tj. sB = ±(B*). Podobně statisticky kolísají hodnoty bodů v obraze A*. Pokud jsou tyto fluktuace příliš vysoké, srovnatelné s průměrnými hodnotami rozdílu mezi A* a B*, mohou se v nich tyto rozdíly snadno "ztratit" a příslušná struktura nebude na obraze patrná. Rušivé statistické fluktuace jsou tak principiálním limitujícím faktorem rozpoznatelnosti malých a ne příliš kontrastních objektů v obraze.

Roseovo kritériumZe statistické analýzy obrazových dat plyne, že rozpoznat (a statisticky prokázat) můžeme v obraze jen takovou strukturu, jejíž kontrast C img splňuje podmínku

Cimg > 4/(B*).

Je to podmínka statistické významnosti rozdílu A*-B* informace v obraze objektu vůči okolnímu fluktuujícímu pozadí B*. Pro kvantitativní popis vlastností obrazů se zavádějí pojmy:Signál s je rozdíl v intenzitě obrazu (počtu nastřádaných impulsů) mezi vyšetřovanou strukturou a okolím. V našem případě je dán rozdílem:

s = A*max-B*.

Šum n (Noise) představuje rušivé statistické fluktuace v obraze. Pro náš případ jsou důležité fluktuace pozadí, takže šum je dán druhou odmocninou z průměrného nastřádaného počtu impulsů v obraze pozadí:

n = sB = (B*).

Shora uvedenou statistickou podmínku detekovatelnosti struktury lze pak vyjádřit tzv.

Roseovým kritériem: SNR > 4.

Vezmeme-li v úvahu vliv rozlišení i statistických fluktuací, spojením shora uvedených vztahů můžeme základní podmínku rozpoznatelnosti struktury zformulovat takto:

2FWHMd

obj *

4C eB

Pomůcky pro kontrolu zobrazení u přístrojů s přímou digitalizací obrazu Portalvision Phantom

C-D diagramProstorové rozlišení lze nejlépe popsat pomocí MTF, kontrast pomocí SNR, popř. CNR.Diagram Contrast-Detail kvalitativně popisuje vzájemný vztah těchto parametrů zobrazovacího systému.Na ose x leží škála velikostí objektu (rozlišení detailů), na ose y pak kontrast objektu.A – lepší prostorové rozlišeníB – lepší rozlišení kontrastu

Dynamický kontrast a kvantová detekční účinnost

Dynamický (expoziční) kontrast:

Kvantová detekční účinnost:

kde NPS je tzv. výkonové spektrum šumu, které se získá Fourierovoutransformací autokorelační funkce ACF

NEQ je ekvivalent šumu (číselně odpovídá Ndet)

max

min

EDRE

2 22

2out

in

PV MTF F NEQ FSNRQDE FSNR NPS F

1

L

ACF x i x x i x dxL

Eliminace šumu

Cílem je snížit rozptyl šumové funkce.Lidskému oku nejvíce vadí vysoké frekvence šumu nízkofrekvenční filtr.Informace o hranách však také leží ve vysoké frekvenci.

Metody pro eliminaci šumu 1. „Časove“ průměrování

Snímá se vícekrát a konečný obraz je průměrem všech obrazů.

`

22

Za předpokladu, že šum je 2D

zredukuje se variance šumu z na N

m clean

m clean

S S n

nS SN

Metody pro eliminaci šumu 2. Obyčejné (souřadnicové)

průměrování

121242121

161

111111111

91

typumaticí sobrazu Konvoluce

C

nebo

C

Pokud jsou hrany známé, tak podél hran používáme jiný filtr.V ostatních částech používáme souřadnicový filtr.

Výsledek konvoluce v bodě měníme pouze když překročil zadaný práh, nebo když leží v daném intervalu (Podle toho jaké vlastnosti má šum a jaké je SNR)

Metody pro eliminaci šumu 3. Průměrování podél hran

Metody pro eliminaci šumu 4. Prahování

Používají se dva typy konvoluční matice 3x3 pro všech 8 směrů (8 sousedů počítaného bodu). Metoda pracuje na okolí 5x5 bodů a to tak, že v osmi směrech od středního bodu počítá rozptyl. Vybere oblast s nejmenším rozptylem, oblast zprůměruje pomocí konvoluce a výslednou hodnotou nahradí bod uprostřed masky 5x5.

Metody pro eliminaci šumu 5. Metoda rotujicího okna

V okně provedeme seřazení dat a prostřední prvek (medián) tvoří výsledek.Výběrové okno je třeba zvolit tak, aby hrany nebyly na koncích.Algoritmus má při klasickém třídění náročnost O(n4) ale pomocí Quicksortu nebo Heapsortu nebo MergeSortu lze redukovat na O(n2 log 2n) případně na O(n log n).Existuje nová verze (FMF, Fast Median Filter) která má náročnost O(n)

Metody pro eliminaci šumu 6. Mediánový Filtr

Metody pro eliminaci šumu 7. Zobecněný Mediánový Filtr

Provedeme seřazení dat a na setříděnou posloupnost aplikujeme váhovou funkci w rovnu např.

1 0,0, 0,1,2,1,0, ,0,04

Medián

w

U šumu typu S&P se většina pixelů správně zobrazí (bez šumu) a jen některé mají falešnou hodnotu 0 nebo 255. Pro takový šum stačí udělat průměrování v bodech majících hodnotu 0 nebo 255. Filtr se zvolí typu:

Metody pro eliminaci šumu 8. Šum typu Salt and Pepper

111101111

81C

Digitalizace obrazuADC charakterizován vzorkovací frekvencí (pixelizace) a bitovou hloubkou (kvantování obrazu). S nižší vzorkovací frekvencí a bitovou hloubkou dochází ke ztrátě informací z původního analogového signálu.Výsledný digitální obraz obsahuje informaci o poloze bodů a amplitudě signálu v těchto bodech.K rekonstrukci obrazu na zobrazovací jednotce (monitoru) je zapotřebí DAC převodník.

high-res image pixelated quantizedpixelated

& quantized

Kvantování (bitová hloubka) obrazu

Kvantování (bitová hloubka) obrazuBitová hloubka určuje rozlišení kontrastu

Vzorkování (pixelizace)Průměrování přes

obdélník

),( yxfDS yxf ,

high-res image pixelated image DS

kde je vzorkovací funkcea je obdélníková funkcese šířkou

NN

N

f f s

s

N

Vzorkování (pixelizace)

Ztratí se informace se vzorkováním?Ukazuje se, že za určitého předpokladu nikolivProces vzorkování lze chápat jako násobení původního obrázku s následujcí funkcí

které říkáme „vzorkovací funkce“ a N je vzorkovací šířka

Vzorkování (pixelizace)

Výsledný obraz je pak diskrétní verzí původního obrazu.Fourierova transformace takového obrazu je konvolucí FT původního obrazu a 2D delta funkce s intervaly 1/N.FT bude pak periodická

Vzorkování (pixelizace) obrazuUrčuje prostorové rozlišení

AliasingMatice detekčních či zobrazovacích elementů je charakterizována vzorkovací šířkou L a šířkou detekčního (zobrazovaccího) elementu – obě nenulové. Dochází ke vzorkování obrazu a zprůměrování obrazu přes šířku elementu. Interval prostorových frekvencí F, které mohou být detekovány či zobrazeny, je dána Nyquistovým kritériem F ≤ 1/(2L).

Aliasing

Aliasing

AliasingDigitální receptory mají obdélníkovou LSF. Mějme obdélníkový puls popsaný funkcí

Fourierova transformace obdélníkového pulzu je tedy

Protože f(x) je sudá funkce, je poslední integrand lichou funkcí a jeho integrál je tudíž nulový. Máme tak

Označme L vzorkovací šířku detektoru (vzdálenost středu dvou sousedních detekčních elementů), F prostorovou frekvenci (lp/mm) signálu.

0 pro 0 pro E x L

f xx L

cos sin cos sinikxF k f x e dx f x kx i kx dx f x kx dx i f x kx dx

00 0

sincos cos sin 2L L

LL

E kLF k f x kx dx E kx dx kx E Lk kL

0

sin 2,

L FMTF F L E L

L F

AliasingNyquistovo kriterium:

Je li prostorová frekvence vstupního signálu vyšší než Flim, dochází k modulaci MTF, což se projeví jako splývání struktur – aliasing. Frekvence výsledného splynutého signálu je o tolik menší než F, o kolik je větší frekvence vstupního signálu oproti F.

lim1

2F

L

0

1 pro a 0 pro

pak sin 2

f x x L f x x L

LFF f k E L

LF

Aliasing

Snímek cihlové zdi pořízený ve vysokém (vlevo) a v nízkém (vpravo) rozlišení

DS

kde je vzorkovací funkcea je obdelníková funkcese šířkou

NN

N

f f s

s

N

DSˆ

ale i opak je pravdou:Ideální rekonstruční vzorkovací filtr

má v k-prostoru tvar obdelníkový1 se šírkou (odstraníme-li vysoké

frekvence které nejsme schopni detekovat) a tudíž v t-oblas

NNf F f s F

N

max

ti bude konvoluční filtr mít tvar znamé

sinc funkce:

2sin xfx

Anti-aliasing filter

Anti-aliasing

Image enhancement methods(Metody pro vylepšení obrazu)

Histogram obrazu je funkce četnosti jednotlivých intenzit. Kumulativní histogram obrazu je funkce intenzit kde pro

každou intenzitu je funkční hodnota rovna počtu bodů majících svou intenzitu menší nebo rovnu.

Ekvalizace (Transformace) Histogramu

Obecně se jedná o transformace histogramu na požadovaný tvar. K transformaci se používá kumulativní histogram.Transformace je možné dělat i lokálně (změna jasu, kontrastu).

0 127 255

012 7

25 5

Změna jasu

Změna jasu

0 127 255

012 7

25 5

g

transform mapping

Změnakontrastu

minmax

minmax

IIIIK

Histogram pro RGB obraz

Ekvalizace Histogramu

Ztráta detailů z ekvalizace

histogramů

Zpracování digitálního obrazuKorekce digitálního obrazu Iraw(x,y) (temný šum, vadné pixely):

kdeD je průměrný temný šum elementů detektoru,G je průměrná odezva elementů detektoru na homogenní signál G(x,y).

Vstupní výkonové spektrum šumu radiačního pole na povrchu detektoru

Výkonové spektrum šumu na výstupu digitálního RTG zobrazovacího zařízení

Detekční kvantová účinnost

, ,

,rawGI x y I x y D

G x y

2,in a in aa

EW u v K SNR K dE E dE

K

2

256 256

1 1 1

, , , exp 2256 256

M

out n k i j i j n i k j

m i j

x yW u v I x y S x y i u x v yM

2 ,, ,

,in

out

W u vDQE u v MTF u v

W u v