ˇ´ıtac ˇova ´ cvic ˇen´ı Poc ˇ ´ho modelova ´ n´ı Skola matematicke 2014
´ c, Marie Sadowsk´a, Maty´aˇs Theuer Petr Beremlijski, Rajko Cosi´
ˇ´ıtac ˇova ´ cvic ˇen´ı Poc ˇ ´ho modelova ´ n´ı Skola matematicke
´ c, Marie Sadowsk´a, Maty´aˇs Theuer Petr Beremlijski, Rajko Cosi´
Katedra aplikovan´e matematiky Fakulta elektrotechniky a informatiky ˇ - Technick´a univerzita Ostrava VSB 2014
Podˇ ekov´ an´ı Autoˇri textu dˇekuj´ı Doc. RNDr. Jiˇr´ımu Bouchalovi, Ph.D. za pozorn´e pˇreˇcten´ı a ˇradu cenn´ych pˇripom´ınek k tomuto textu.
Tato akce byla podpoˇrena z prostˇredk˚ u projektu Matematika s radost´ı (registraˇcn´ı ˇc´ıslo CZ.1.07/1.1.00/26.0042). Jeho webov´ a str´ anka je http://msr.vsb.cz.
ˇen´ı 1: Matlab – na ´stroj pro matematick´ ´n´ı Cvic e modelova Abychom se mohli vˇenovat numerick´emu ˇreˇsen´ı matematick´ych u ´ loh, potˇrebujeme vhodn´e prostˇred´ı, kter´e n´am to umoˇzn´ı. A tak jako fyzici ˇci chemikov´e maj´ı sv´e laboratoˇre, maj´ı i numeriˇct´ı matematici svou Maticovou laboratoˇr1 - Matlab. Podrobnˇe se tomuto pracovn´ımu prostˇred´ı a jeho pˇr´ıkaz˚ um vˇenuje pˇriloˇzen´y Matlabovsk´y slabik´aˇr [2] nebo tak´e u ´ vod textu [1]. My si v tomto textu uvedeme pouze struˇcn´y pˇrehled matlabovsk´ych promˇenn´ych a pˇr´ıkaz˚ u, kter´e budeme potˇrebovat. Prostˇ red´ı • help, demos, intro, who, whos, clear, size, length Promˇ enn´ e • Skal´ary • Vektory • Matice Pˇ r´ıkazy • Skal´arn´ı funkce - sin, cos, tan, cot, exp, log, abs, sqrt, round • Vektorov´e funkce a generov´an´ı vektor˚ u - max, min, sort • Maticov´e funkce a generov´an´ı matic - det, rand, ones, zeros, eye • Skal´arn´ı operace - +, −, ∗, /, b
• Maticov´e a vektorov´e operace - +, −, ∗, ´(transponov´an´ı), \ (A\v = x ⇔ Ax = v) Operace po prvc´ıch“ - .∗ , .b, ./ ” • 2D grafika (vykreslen´ı graf˚ u funkc´ı jedn´e promˇenn´e) - plot, hold on, hold off, figure • 3D grafika (vykreslen´ı graf˚ u funkc´ı dvou promˇenn´ych) - meshgrid, mesh, contour, hold on, hold off, figure ˇ ıd´ıc´ı pˇr´ıkazy - if (podm´ınˇen´y pˇr´ıkaz), for, while (pˇr´ıkazy cyklu se zn´am´ym poˇctem • R´ opakov´an´ı a podm´ınkou na zaˇc´atku) • Relace a logick´e operace - <, >, <=, >=, ==, ∼=, &, |, ∼ • Skripty a funkce - function 1
MATrix LABoratory
1
Vˇse si nyn´ı vyzkouˇs´ıme pˇri ˇreˇsen´ı n´asleduj´ıc´ıch u ´ loh. ´ Ukol 1.1 Kolik ˇclen˚ u harmonick´e ˇrady2 mus´ıme nejm´enˇe seˇc´ıst, aby tento ˇc´asteˇcn´y souˇcet ˇrady mˇel hodnotu alespoˇ n 10 (15, 20)? N ´ Ukol 1.2 Zkuste odhadnout s vyuˇzit´ım Matlabu souˇcet ˇrady ∞ X 1 n=1
1 · . n n+1 N
´ Ukol 1.3 Sestrojte grafy n´asleduj´ıc´ıch funkc´ı: • f (x) := x2 √ • f (x) := 1 − x2 • f (x) := x2 sin
1 x2
• f (x) := |x|
N
´ Ukol 1.4 Sestrojte grafy n´asleduj´ıc´ıch funkc´ı: • f (x, y) := x2 + y 2 p • f (x, y) := x2 + y 2 2
2
• f (x, y) := (x + y ) sin • f (x, y) :=
p
|xy|
1 x2 +y 2
N
2ˇ
Radou (re´ aln´ ych ˇ c´ ısel) rozum´ımePv´ yraz a1 + a2 + · · · + an + · · · =: ∞ je an ∈ R. Harmonickou naz´ yv´ame ˇradu n=1 n1 .
2
P∞
n=1
an , kde pro kaˇzd´e n ∈ N
ˇen´ı 2: Monte Carlo Cvic V tomto cviˇcen´ı se budeme vˇenovat odhad˚ um Ludolfova ˇc´ısla π. Nejprve hledejme obsah kruhu o polomˇeru r (pˇredstavme si na chv´ıli, ˇze nezn´ame pˇr´ısluˇsn´y vzorec a nic nev´ıme o ˇc´ıslu π). Prvn´ım n´apadem by mohlo b´yt zhotoven´ı v´alcov´ych n´adob s r˚ uzn´ymi polomˇery podstav (napˇr. s polomˇery 1 a 2 jednotky) a jednotkovou v´yˇskou, viz obr. 1. Objem vody, kter´y se do takov´ych n´adob vejde, je roven obsahu podstavy v´alce, tj. ob-
v =1
v =1
r=1
r=2
Obr´azek 1: V´alcov´e n´adoby sahu kruhu s polomˇerem r. Rychle si vˇsimneme, ˇze pokud zvˇetˇs´ıme polomˇer podstavy dvakr´at, zvˇetˇs´ı se objem ˇctyˇrikr´at, dalˇs´ımi experimenty m˚ uˇzeme odvodit, ˇze obsah kruhu je pˇr´ımo u ´ mˇern´y druh´e mocninˇe polomˇeru. Tak´e zjist´ıme, ˇze druhou mocninu polomˇeru kruhu mus´ıme vyn´asobit vˇzdy stejnou konstantou, abychom dostali obsah dan´eho kruhu. Tuto konstantu oznaˇc´ıme π. Existuje mnoho moˇznost´ı, jak tuto konstantu odhadnout. Snad nejjednoduˇsˇs´ım zp˚ usobem, jak stanovit interval, ve kter´em leˇz´ı π, je vepsat do kruhu o polomˇeru r ˇctverec a stejn´emu kruhu opsat jin´y ˇctverec, viz obr. 2. D´elka strany vˇetˇs´ıho
r r
Obr´azek 2: Odhad intervalu obsahuj´ıc´ı π √ ˇctverce je 2r a z Pythagorovy vˇety d´ale zjist´ıme, ˇze d´elka strany menˇs´ıho ˇctverce je 2r, tud´ıˇz obsah vˇetˇs´ıho resp. menˇs´ıho ˇctverce je 4r 2 resp. 2r 2. Jelikoˇz jsme si uˇz odvodili“, ” ˇze obsah kruhu o polomˇeru r je d´an vzorcem πr 2 , pouh´ym porovn´an´ım obsah˚ u ˇctverc˚ ua kruhu zjist´ıme, ˇze π ∈ (2, 4). Mnohem rozumnˇejˇs´ı v´ysledek z´ısk´ame, pokud budeme dan´emu kruhu vepisovat n-´ uheln´ıky a poˇc´ıtat jejich obsahy. Tuto metodu naz´yv´ame vyˇcerp´avac´ı (exhaustn´ı) a pravdˇepodobnˇe prvn´ı ji pouˇzil Eudoxos3 . Neˇz se budeme vˇenovat odhad˚ um ˇc´ısla π, pod´ıvejme se kr´atce na v´ypoˇcet obvodu kruhu. Jistˇe v´ıme, ˇze obvod kruhu je pˇr´ımo u ´ mˇern´y dvojn´asobku jeho polomˇeru. Ot´azkou je, jak´a 3
Eudoxos (410 nebo 408 pˇr. n. l. – 355 nebo 347 pˇr. n. l.) – ˇreck´ y astronom, matematik a fyzik, student Plat´ona
3
je pˇr´ısluˇsn´a konstanta u ´ mˇernosti k. Na obr´azku 3 provedeme pˇreuspoˇr´ad´an´ı kruhu na u ´ tvar, kter´y se pro zjemˇ nuj´ıc´ı se dˇelen´ı kruhu bl´ıˇz´ı obd´eln´ıku. Porovn´an´ım obsahu kruhu (πr 2 ) a obsahu vznikl´eho obd´eln´ıku (kr 2 ) lze odvodit, ˇze konstanta k je opˇet rovna π. o = 2kr r
r kr
zjemnˇen´ı dˇelen´ı
r kr k=π
Obr´azek 3: Stanoven´ı vzorce pro obvod kruhu Nyn´ı si uk´aˇzeme nˇekolik zp˚ usob˚ u, jak nal´ezt pˇribliˇznou hodnotu ˇc´ısla π. Buffonova metoda ´ D˚ usledkem ˇreˇsen´ı tzv. Buffonova4 probl´emu s jehlou je aproximace ˇc´ısla π. Uloha spoˇc´ıv´a v opakovan´em h´azen´ı jehly o d´elce ℓ na rovinu, na kter´e m´ame vyznaˇcenu s´ıt’ rovnobˇeˇzek se vzd´alenost´ı 2ℓ. Jestliˇze jehlu hod´ıme n-kr´at a x-kr´at n´am bˇehem tˇechto pokus˚ u po dopadu zkˇr´ıˇz´ı nˇekterou z rovnobˇeˇzek, pak ˇc´ıslo n x aproximuje ˇc´ıslo π. V roce 1975 Perlman a Wichura5 publikovali n´asleduj´ıc´ı v´ysledek t´ykaj´ıc´ı se pˇresnosti Buffonovy √ metody. S pravdˇepodobnost´ı 95 procent nem´a chyba aproximace hodnotu vˇetˇs´ı neˇz 5/ n. Tzn. napˇr´ıklad pro 104 pokus˚ u n´am s pravdˇepodobnost´ı 95 procent chyba nepˇrekroˇc´ı hodnotu 0,05. ´ Ukol 2.1 Implementujte Buffonovu metodu a pouˇzijte ji k aproximaci ˇc´ısla π. Porovnejte vaˇsi aproximaci se skuteˇcnou“ hodnotou ˇc´ısla π a urˇcete chybu aproximace. N ” 4 5
G. L. Buffon (1707–1788) – francouzsk´ y pˇr´ırodovˇedec M. Perlman a M. Wichura – ameriˇct´ı matematikov´e
4
Metoda Monte Carlo Monte Carlo je tˇr´ıda v´ypoˇcetn´ıch algoritm˚ u zaloˇzen´a na prov´adˇen´ı n´ahodn´ych experiment˚ u. T´eto metody se ˇcasto pouˇz´ıv´a pro simulaci fyzik´aln´ıch a matematick´ych syst´em˚ u. V´ysledkem proveden´ı velk´eho mnoˇzstv´ı experiment˚ u je obvykle pravdˇepodobnost urˇcit´eho jevu. Na z´akladˇe z´ıskan´e pravdˇepodobnosti a zn´am´ych vztah˚ u pak spoˇc´ıt´ame potˇrebn´e v´ysledky. Protoˇze metoda vyˇzaduje generov´an´ı velk´eho souboru n´ahodn´ych dat, je vhodn´e pro jej´ı implementaci pouˇzit´ı poˇc´ıtaˇce. Metod Monte Carlo se pouˇz´ıv´a v pˇr´ıpadˇe, kdy je pˇr´ıliˇs pracn´e nebo nemoˇzn´e nal´ezt pˇresn´y v´ysledek jin´ym zp˚ usobem. Jej´ı v´yhodou je jednoduch´a implementace, nev´yhodou relativnˇe mal´a pˇresnost. Metoda byla vytvoˇrena skupinou fyzik˚ u pracuj´ıc´ıch na projektu jadern´e pumy v Los Alamos, jm´eno metody bylo navrˇzeno v roce 1940 von Neumannem6 . V matematice se Monte Carlo pouˇz´ıv´a zejm´ena pro v´ypoˇcet urˇcit´ych integr´al˚ u (zejm´ena v´ıcen´asobn´ych urˇcit´ych integr´al˚ u), kter´e je obt´ıˇzn´e ˇci nemoˇzn´e vyˇc´ıslit analyticky nebo jinou vhodnou numerickou metodou. Napˇr.R obsah plochy E ohraniˇcen´e grafem funkce y = 1 x2 , osou x a pˇr´ımkou x = 1 (tj. S(E) = 0 x2 dx, viz obr. 4) je moˇzn´e metodou Monte Carlo aproximovat n´asleduj´ıc´ım zp˚ usobem. Necht’ n´aˇs program generuje n´ahodnˇe dvojice ˇc´ısel [x, y], pˇriˇcemˇz kaˇzd´e z ˇc´ısel x a y je vybr´ano nez´avisle z intervalu h0, 1i. Tuto dvojici budeme ch´apat jako souˇradnice bodu, kter´y je n´ahodnˇe zvolen ve ˇctverci h0, 1i × h0, 1i. Pravdˇepodobnost toho, ˇze bod leˇz´ı uvnitˇr zadan´eho ˇctverce, je 1. Pravdˇepodobnost toho, ˇze bod leˇz´ı uvnitˇr podmnoˇziny E zadan´eho ˇctverce, je rovna obsahu plochy E, tj. S(E). Takˇze obsah plochy, kter´a je podmnoˇzinou zvolen´eho ˇctverce, m˚ uˇzeme odhadnout jako pravdˇepodobnost, ˇze n´ahodnˇe zvolen´y bod z dan´eho ˇctverce leˇz´ı v t´eto podmnoˇzinˇe. 1
1
S(E) =
y
y
x2 Z1
x2 dx
0
0
0
0
1
0
1
x
x
Obr´azek 4: Plocha E a Monte Carlo ´ Ukol 2.2 Implementujte metodu Monte Carlo pro pˇribliˇzn´y v´ypoˇcet v´ysledky porovnejte s analytick´ym vyˇc´ıslen´ım integr´alu. 6
John von Neumann (1903–1957) – v´ yznamn´ y mad’arsk´ y matematik
5
R1 0
x2 dx. Z´ıskan´e N
Pro odhad pˇresnoti metody Monte Carlo plat´ı n´asleduj´ıc´ı√tvrzen´ı. S pravdˇepodobnost´ı 95 procent nem´a chyba aproximace hodnotu vˇetˇs´ı neˇz 1/ n. Tzn. napˇr´ıklad pro 104 pokus˚ u n´am s pravdˇepodobnost´ı 95 procent chyba nepˇrekroˇc´ı hodnotu 0,01. Chceme-li vyuˇz´ıt metodu Monte Carlo k aproximaci ˇc´ısla π, vypoˇcteme pˇribliˇznˇe touto metodou (napˇr´ıklad) integr´al Z 1√ 1 − x2 dx. 0
Snadno si uvˇedom´ıme (viz obr. 5), ˇze t´ımto zp˚ usobem z´ısk´ame aproximaci hodnoty π/4.
´ Ukol R 1 √ 2.3 Implementujte metodu Monte Carlo pro pˇribliˇzn´y v´ypoˇcet integr´alu 1 − x2 dx a pouˇzijte ji k aproximaci ˇc´ısla π. Porovnejte vaˇsi aproximaci se skuteˇcnou“ 0 ” hodnotou ˇc´ısla π a urˇcete chybu aproximace. N
1
1
y
1 − x2
y
√
Z1 √
1 − x2 dx
0
0
0
0
1
0
x
1
x
Obr´azek 5: π a Monte Carlo
Pˇripomeˇ nme jeˇstˇe, ˇze v´yˇse uveden´e pouˇzit´ı metody Monte Carlo je pouze ilustrativn´ı; takto jednoduch´y integr´al z u ´ kolu 2.3 lze velmi efektivnˇe spoˇc´ıst vhodnou numerickou kvadraturou.
6
Aproximace π pomoc´ı ˇ c´ıseln´ ych ˇ rad Posledn´ı metodou nalezen´ı aproximace ˇc´ısla π, kterou si v tomto pˇrehledu uk´aˇzeme, je vyuˇzit´ı ˇc´ıseln´ych ˇrad. K t´eto aproximaci pouˇzijeme tzv. Gregoryho7 ˇradu, kter´a je rozvojem funkce arctg x: X x3 x5 x7 x2n−1 arctg x = x − + − + · · · =: (−1)n+1 . 3 5 7 2n − 1 n=1 ∞
Tato ˇrada m´a koneˇcn´y souˇcet pro x ∈ h−1, 1i (ˇr´ık´ame, ˇze ˇrada konverguje), nav´ıc plat´ı, ˇze ˇc´ım v´ıce je |x| menˇs´ı neˇz 1, t´ım m´enˇe ˇclen˚ u ˇrady potˇrebujeme pouˇz´ıt k nahrazen´ı arctg x s uspokojivou“ pˇresnost´ı. Prvn´ı moˇznost´ı, jak aproximovat π, je tud´ıˇz nahradit arctg 1 ” ˇc´asteˇcn´ym souˇctem Gregoryho ˇrady, jelikoˇz arctg 1 = π/4. Aproximaci s rychlejˇs´ı konvergenc´ı z´ısk´ame, pokud pouˇzijeme rovnost arctg 1 = arctg (1/2) + arctg (1/3),
(1)
viz obr. 6.
arctg 1
arctg (1/2)
arctg (1/3)
Obr´azek 6: Ilustrace k odvozen´ı vztahu (1)
´ Ukol 2.4 Implementujte metodu, kter´a vyuˇzije Gregoryho ˇradu k nalezen´ı aproximace ˇc´ısla π. N
7
James Gregory (1638 – 1675) – skotsk´ y matematik a astronom
7
ˇen´ı 3: Obra ´zky jako matice Cvic V dob´ach analogov´ych fotoapar´at˚ u se vyfocen´y obraz v kaˇzd´em okamˇziku uchov´aval skuteˇcnˇe jako obraz. At’ uˇz na negativu filmu, nebo po vyvol´an´ı na fotografick´em pap´ıˇre ˇci diapozitivu. Dnes je tomu jinak. Pˇri zm´aˇcknut´ı spouˇstˇe digit´aln´ıho fotoapar´atu se aktu´aln´ı sc´ena pˇred objektivem zachyt´ı pomoc´ı sn´ımac´ıho ˇcipu a uloˇz´ı jako soubor ˇc´ısel na pamˇet’ovou kartu. I kdybychom tuto kartu rozebrali, obr´azky na n´ı neuvid´ıme. K jejich zobrazen´ı potˇrebujeme opˇet nˇejak´e digit´aln´ı zaˇr´ızen´ı, kter´e um´ı obraz uloˇzen´y v jedniˇck´ach a nul´ach pˇrev´est do viditeln´e formy. Matice obrazu Pro jednoduchost se nejprve zab´yvejme ˇcernob´ıl´ymi obr´azky - pˇresnˇeji obr´azky ve stupn´ıch ˇsed´e. V okamˇziku exponov´an´ı dojde k zachycen´ı sn´ıman´e sc´eny na ˇcip, kter´y je tvoˇren soustavou miniaturn´ıch fotodiod uspoˇr´adan´ych do ˇr´adk˚ u a sloupc˚ u. V z´avislosti na osvˇetlen´ı pˇr´ısluˇsn´e ˇc´asti ˇcipu kaˇzd´a z diod vyprodukuje urˇcit´y n´aboj, kter´y je zmˇeˇren, a jeho hodnota je zaznamen´ana ve formˇe ˇc´ısla. V pˇr´ıpadˇe JPEG obr´azku se jedn´a o cel´e ˇc´ıslo od 0 do 255,8 pˇriˇcemˇz hodnota 0 odpov´ıd´a ˇcern´e a hodnota 255 b´ıl´e. Proces rozdˇelen´ı spojit´eho obrazu na diskr´etn´ı hodnoty v jednotliv´ych oddˇelen´ych bodech se naz´yv´a kvantov´an´ı a vzorkov´an´ı.
Obr´azek 7: Pˇr´ıklad kvantov´an´ı a vzorkov´an´ı Pomineme-li ve skuteˇcnosti bin´arn´ı uloˇzen´ı obrazu, m˚ uˇzeme si jeho digit´aln´ı podobu pˇredstavit jako obecnˇe obd´eln´ıkovou tabulku ˇc´ısel, ve kter´e kaˇzd´a hodnota reprezentuje jas urˇcit´e mal´e ploˇsky (pixelu) v zachycen´em obraze. Tuto tabulku budeme oznaˇcovat pojmem matice obrazu: a1,1 a1,2 · · · a1,n a2,1 a2,2 · · · a2,n A = .. .. .. , kde ai,j ∈ {0, 1, 2, ..., 255} . . . . . . . am,1 am,2 · · · am,n Pˇr´ıklad ˇc´asti matice obrazu m˚ uˇzeme vidˇet na obr. 8. 8
Tyto hodnoty odpov´ıdaj´ı osmibitov´e reprezentaci ˇc´ısla.
8
41 38 34 35 38 62 64 58 42 66 120 137 127 111 57 59
51 41 35 39 38 68 60 46 35 77 133 134 122 126 80 81
51 71 50 50 70 63 52 43 65 101 152 142 118 107 92 101
44 64 68 67 96 49 47 39 82 150 167 142 118 113 98 98
49 46 40 55 63 50 37 47 136 177 165 137 149 119 107 101
75 47 43 33 38 40 63 124 177 179 165 146 148 127 111 94
106 78 30 32 43 72 138 184 182 171 177 154 138 124 92 59
139 133 96 75 122 163 190 177 178 170 175 157 152 126 90 49
182 190 177 182 194 203 189 169 175 184 174 164 155 140 111 78
175 174 171 177 176 186 179 170 180 181 186 161 158 137 111 89
160 169 174 168 174 176 177 172 173 189 185 158 149 124 100 62
155 163 169 173 180 188 182 177 187 198 181 154 152 124 92 46
137 156 159 171 179 181 179 172 180 188 173 156 155 143 119 86
134 174 171 172 173 173 173 167 173 176 160 151 153 144 134 126
102 155 168 174 176 181 173 168 177 172 148 145 157 142 123 116
51 97 137 164 177 178 169 160 156 146 136 147 166 162 125 103
ˇ ast matice obrazu Obr´azek 8: C´ Histogram M´ame-li obr´azek ve stupn´ıch ˇsed´e (ˇc´ıslech od 0 do 255), m˚ uˇze n´as zaj´ımat, kolikr´at se v obr´azku jednotliv´e stupnˇe (ˇc´ısla) vyskytuj´ı. Graf ˇcetnosti v´yskytu jednotliv´ych hodnot v obr´azku se naz´yv´a histogram.
Obr´azek 9: Obr´azek a jeho histogram Pˇrestoˇze obr´azek nen´ı sv´ym histogramem jednoznaˇcnˇe definov´an, m˚ uˇzeme z histogramu 9
o p˚ uvodn´ım obr´azku hodnˇe vyˇc´ıst. Tmav´e obr´azky maj´ı v histogramu velk´e hodnoty nahromadˇen´e v lev´e ˇc´asti grafu. Naopak svˇetl´e obr´azky maj´ı v histogramu velk´e hodnoty v prav´e ˇc´asti grafu. Histogram obr´azk˚ u s n´ızk´ym kontrastem vypad´a jako jeden relativnˇe u ´ zk´y ko” pec“, zat´ımco obr´azky s vysok´ym kontrastem maj´ı vˇetˇsinou dva kopce“, kaˇzd´y na opaˇcn´e ” stranˇe grafu.
Obr´azek 10: Pˇr´ıklad r˚ uzn´ych obr´azk˚ u a jejich histogram˚ u
10
Jednoduch´ eu ´pravy Jak jsme si ˇrekli, obr´azky jsou v poˇc´ıtaˇci uloˇzeny jako matice ˇc´ısel, takˇze m´ısto manipulace s barevn´ymi ploˇskami n´am staˇc´ı prov´adˇet operace s ˇc´ısly. Nyn´ı si v jednoduchosti pop´ıˇseme nˇekter´e z´akladn´ı operace s maticemi obr´azk˚ u. Pˇredpokl´adejme, ˇze m´ame obr´azek ve stupn´ıch ˇsed´e, kter´y m´a 512×512 pixel˚ u. Matice takov´eho obrazu, kterou oznaˇc´ıme Iorig , m´a 512 ˇr´adk˚ u a 512 sloupc˚ u. Dohromady je to 262144 ˇc´ısel, takˇze je jasn´e, ˇze nebudeme poˇc´ıtat s kaˇzdou hodnotou ruˇcnˇe, ale nastav´ıme nˇejak´e pravidlo, kter´e poˇc´ıtaˇci ˇrekne, co m´a s jednotliv´ymi hodnotami udˇelat. • Oˇrez obr´azku provedeme jednoduˇse tak, ˇze vybereme jen omezen´y rozsah ˇr´adk˚ ua sloupc˚ u p˚ uvodn´ı matice. Napˇr´ıklad pro v´yˇrez prav´e horn´ı ˇctvrtiny obr´azku Iorig vezmeme pouze 1. aˇz 256. ˇr´adek a 257. aˇz 512. sloupec. Dostaneme tam matici Icrop , kter´a m´a 256 ˇr´adk˚ u a 256 sloupc˚ u. • Zesvˇetlen´ı obr´azku provedeme napˇr´ıklad tak, ˇze k hodnotˇe kaˇzd´eho pixelu v matici Iorig pˇriˇcteme nˇejak´e kladn´e ˇc´ıslo, napˇr´ıklad 100. Touto operac´ı se vˇsak m˚ uˇze st´at, ˇze se nˇekter´e hodnoty ocitnou mimo rozsah 0 aˇz 255. Oprav´ıme to napˇr´ıklad tak, ˇze vˇsechny hodnoty, kter´e vyjdou vˇetˇs´ı neˇz 255, zmenˇs´ıme pr´avˇe na 255. Ve vˇsech bodech, kde doˇslo k takov´em oˇrezu hodnot ovˇsem ztr´ac´ıme obrazovou informaci. • Ztmaven´ı obr´azku m˚ uˇzeme prov´est obdobnˇe jako zesvˇetlen´ı, ale mus´ıme zkontrolovat, zda nˇekter´e hodnoty nevyˇsly z´aporn´e. Pokud ano, nastav´ıme je na 0. Dalˇs´ı moˇznost´ı je vyn´asobit vˇsechny hodnoty matice obrazu nˇejak´ym ˇc´ıslem z intervalu (0, 1). Zde nehroz´ı, ˇze bychom se dostali mimo rozsah 0 aˇz 255, ale m˚ uˇze se st´at, ˇze v´ysledkem nebudou jen cel´a ˇc´ısla, a tak je potˇreba kaˇzdou v´yslednou hodnotu zaokrouhlit na nejbliˇzˇs´ı cel´e ˇc´ıslo. • Negativ obr´azku je obraz, kter´y m´a obr´acenou reprezentaci ˇc´ısel. V klasick´em obr´azku odpov´ıd´a 0 ˇcern´e a 255 b´ıl´e hodnoty mezi tˇemito ˇc´ısly odpov´ıdaj´ı r˚ uzn´ym stupˇ n˚ um ˇsed´e od tmav´e po svˇetlou. V negativn´ım zobrazen´ı je to naopak. Abychom nemuseli pˇreprogramovat zobrazovac´ı zaˇr´ızen´ı k obr´acen´ı stupnice, m˚ uˇzeme obr´atit stupnici pˇr´ımo v naˇsem obrazu. Hodnoty 0 zmˇen´ıme na 255 a naopak. Vˇsechny hodnoty tedy vypoˇc´ıt´ame tak, ˇze p˚ uvodn´ı hodnotu odeˇcteme od ˇc´ısla 255 a dostaneme pr´avˇe negativ p˚ uvodn´ıho obr´azku, kter´y jiˇz zobraz´ıme klasicky. • Prahov´an´ı je jednoduch´a operace, kde zvol´ıme t ∈ (0, 255) (jednoduch´y pr´ah). Jednotliv´e prvky matice Iprah nastav´ıme na 0, pokud pˇr´ısluˇsn´a hodnota matice Iorig je menˇs´ı neˇz t, a na 255, pokud je pˇr´ısluˇsn´a hodnota matice Iorig je vˇetˇs´ı nebo rovna t. • Vyrovn´an´ı histogramu je metoda, kter´a slouˇz´ı k automatick´emu upraven´ı kontrastu tak, aby byla ˇcetnost jednotliv´ych hodnot rozdˇelena pˇribliˇznˇe rovnomˇernˇe. Tato metoda je automatick´a a nezohledˇ nuje, o jak´y obr´azek se jedn´a, takˇze v´ysledek m˚ uˇze p˚ usobit ponˇekud umˇele.
11
Obr´azek 11: Negativ obr´azku (vlevo), obr´azek po prahovan´ı s parametrem t = 50 (uprostˇred) a t = 150 (vpravo)
Obr´azek 12: Pˇr´ıklad vyrovn´an´ı histogramu obr´azku s n´ızk´ym kontrastem Barevn´ e obr´ azky Pro reprezentaci barevn´eho obrazu se pouˇz´ıvaj´ı matice rovnou tˇri - jedna pro kaˇzdou barevnou sloˇzku RGB9 . V´ysledn´y obraz pak vznikne sloˇzen´ım tˇechto tˇr´ı obraz˚ u. Hodnota kaˇzd´eho pixelu obrazu se tedy skl´ad´a ze tˇr´ı ˇc´ısel, a tak se tato ploˇska zobrazuje v jedn´e z 2563 moˇzn´ych barev. Pro jin´y form´at souboru m˚ uˇze matice obrazu obsahovat i jin´a ˇc´ısla neˇz cel´a ˇc´ısla od 0 do 255. M˚ uˇzeme napˇr´ıklad pouˇz´ıt re´aln´a ˇc´ısla z intervalu (0, 1), kdy ˇrekneme, ˇze ˇcern´a je 0 9
RGB je syst´em rozloˇzen´ı barevn´eho obrazu na ˇcerven´ y (R), zelen´ y (G) a modr´ y (B) kan´al.
12
a b´ıl´a odpov´ıd´a hodnotˇe 1. Pokud uloˇz´ıme hodnoty kaˇzd´eho pixelu ve form´atu 14-bitov´ych ˇc´ısel, dostaneme aˇz 163843 barev. Bˇeˇzn´e zobrazovac´ı zaˇr´ızen´ı vˇsak neum´ı takov´y rozsah barev zobrazit, a proto ve vˇetˇsinˇe pˇr´ıpad˚ u staˇc´ı ukl´adat obr´azky do osmi bit˚ u.
Obr´azek 13: JPEG obr´azek sloˇzen´y z RGB kan´al˚ u ´ Ukol 3.1 Naˇctˇete sv˚ uj obr´azek do Matlabu pomoc´ı pˇr´ıkazu imread. Zobrazte tento obr´azek pomoc´ı funkce imshow a zobrazte jeho jednotliv´e sloˇzky R, G a B. Pˇreved’te tento obr´azek do stupˇ n˚ u ˇsed´e pomoc´ı funkce rgb2gray a uloˇzte pomoc´ı funkce imwrite. N
13
´ Ukol 3.2 Proved’te libovoln´y v´yˇrez vaˇseho ˇcernob´ıl´eho obr´azku obsahuj´ıc´ı 512 × 512 pixel˚ u a uloˇzte tento oˇr´ıznut´y obr´azek. Zobrazte histogram tohoto obr´azku pomoc´ı funkce imhist. N ´ Ukol 3.3 Proved’te zesvˇetlen´ı a ztmaven´ı vaˇseho obr´azku a vytvoˇrte jeho negativ. Pro’ ved te prahov´an´ı s prahem 10, 100, 150, 200. N ´ Ukol 3.4 Proved’te zv´yˇsen´ı kontrastu vaˇseho obr´azku pomoc´ı vyrovn´an´ı histogramu. N Na z´avˇer poznamenejme, ˇze mnohem v´ıce lze o problematice zpracov´an´ı obrazu nal´ezt v textu [3].
14
ˇen´ı 4: Konvoluce Cvic V tomto cviˇcen´ı si zkus´ıme aplikovat na obr´azky diskr´etn´ı konvoluci. Diskr´etn´ı konvoluce je zobrazen´ı, do kter´eho vstupuje obr´azek a takzvan´a konvoluˇcn´ı maska. Konvoluˇcn´ı maska specifikuje, co konkr´etnˇe konvoluce s obr´azkem provede. Protoˇze se v poˇc´ıtaˇc´ıch obr´azky vyskytuj´ı nejˇcastˇeji ve formˇe matice, jej´ıˇz kaˇzd´a hodnota odpov´ıd´a jasu v dan´em pixelu, nepˇrekvap´ı n´as, kdyˇz i konvoluˇcn´ı maska bude ze stejn´eho svˇeta. Diskr´etn´ı konvoluce tedy vezme matici obr´azku a kaˇzd´emu bodu pˇriˇrad´ı novou hodnotu podle n´asleduj´ıc´ıho sch´ematu:
Obr´azek 14: Sch´ema diskr´etn´ı konvoluce Jednoduˇse ˇreˇceno: poloˇz´ıme masku na obr´azek tak, aby stˇred masky byl v bodˇe, kde konvoluci poˇc´ıt´ame, pˇren´asob´ıme kaˇzdou hodnotou v masce pˇr´ısluˇsnou hodnotu obr´azku a potom vˇse seˇcteme10 . Pomoc´ı konvoluce jsme schopni odstranit z obr´azk˚ u ˇsum nebo zv´yraznit hrany, coˇz je hojnˇe vyuˇz´ıv´ano, chceme-li, aby poˇc´ıtaˇc dok´azal s´am rozpoznat objekty na obr´azku (tˇreba pˇri zpracov´an´ı dat z pr˚ umyslov´ych kamer). ˇ ım je maska vˇetˇs´ı, t´ım vˇetˇs´ı okol´ı bodu n´am m˚ C´ uˇze zas´ahnout do v´ypoˇctu. N´azornˇe je to vidˇet v n´asleduj´ıc´ım pˇr´ıkladu, kdy byla nejdˇr´ıve pouˇzita maska 1 1 1 1 1 1 1 , · 9 1 1 1 kter´a vlastnˇe novou hodnotu spoˇc´ıt´a jako pr˚ umˇer hodnot v okoln´ıch bodech. Pot´e na stejn´y obr´azek aplikujeme masku vˇetˇs´ı, ale opˇet s hodnotami, kter´e zpr˚ umˇeruj´ı hodnoty obr´azku: 10
Konvoluci vypoˇcteme pomoc´ı pˇredpisu (f ∗ h) (x, y) =
popisuje obr´azek a funkce h masku.
15
k P
k P
r=−k s=−k
f (x + r, y + s) h (r, s), kde funkce f
1 25
·
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
.
Obr´azek 15: Lena pˇred aplikac´ı konvoluce (vlevo), po konvoluci s maskou 3 × 3 (uprostˇred) a po konvoluci s maskou 5 × 5 (vpravo) Jak si m˚ uˇzeme vˇsimnout, naˇse maska zpr˚ umˇerov´an´ım hodnot na okol´ı obraz rozmazala, a to t´ım v´ıce, ˇc´ım vˇetˇs´ı okol´ı bodu se do pr˚ umˇeru zapoˇc´ıtalo. Kdyby byla maska stejnˇe velik´a jako obr´azek, tak by cel´y v´ysledn´y obr´azek byl jen v jedn´e barvˇe. V minul´em cviˇcen´ı jsme si uk´azali, jak se v Matlabu pracuje s maticemi obrazu. Nyn´ı si zkus´ıme naprogramovat diskr´etn´ı konvoluci podle n´asleduj´ıc´ıho algoritmu. Algoritmus 1: Diskr´etn´ı konvoluce o br a zek = n a c t i o b r a z e k ( ) N,M = rozmery ( o br a zek ) maska = ma t ice ( n , n ) for i = 1 . .N for j = 1 . .M r e s t r i k c e = o br a zek ( ( i −n / 2 ) : ( i+n / 2 ) , ( j−n / 2 ) : ( j+n / 2 ) ) no vy o br a zek ( i , j ) = suma prvku ( maska . ∗ r e s t r i k c e ) end end u l o z o b r a z e k ( no vy o br a zek )
Abychom se vyhnuli pˇreˇcn´ıv´an´ı masky ,,do pr´azdna“ pˇri poˇc´ıt´an´ı hodnot u krajn´ıch bod˚ u, nebudeme konvoluci poˇc´ıtat v bodech, ve kter´ych by maska pˇreˇcn´ıvala pˇres okraj obr´azku. Vytvoˇr´ıme si vlastnˇe takov´y r´ameˇcek, na kter´em v´ypoˇcet konvoluce prostˇe vynech´ame. 16
´ Ukol 4.1 Naprogramujte diskr´etn´ı konvoluci podle n´avodu v´yˇse a na obr´azc´ıch Lena.png a bricks.jpg (m˚ uˇzete pouˇz´ıt i sv´e obr´azky z minul´eho cviˇcen´ı) otestujte n´asleduj´ıc´ı konvoluˇcn´ı masky. Obr´azky naˇc´ıtejte v odst´ınech ˇsedi. 1 1 1 • 91 · 1 1 1 . . . odstran´ı ˇsum 1 1 1 1 2 1 1 • 16 · 2 4 2 . . . Gaussova maska 1 2 1 0 1 0 • 1 −4 1 . . . Laplaceova maska 0 1 0 0 −1 0 5 −1 . . . zv´yrazn´ı hrany • −1 0 −1 0 −1 0 1 • −2 0 2 . . . zv´yrazn´ı svisl´e hrany −1 0 1 −1 −2 −1 0 0 . . . zv´yrazn´ı vodorovn´e hrany • 0 1 2 1
U posledn´ıch ˇctyˇr masek pˇriˇctˇete ke kaˇzd´emu bodu jeˇstˇe hodnotu 128, at’ obr´azky nejsou moc tmav´e. N Pro ilustraci praktick´eho vyuˇzit´ı diskr´etn´ı konvoluce pouˇzijeme na zaˇsumˇen´y obr´azek konvoluˇcn´ı masku 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
viz obr. 16. Jak je vidˇet, pouˇzili jsme stejnou masku, na kter´e jsme si ukazovali rozmaz´an´ı obr´azku. Pˇri odstraˇ nov´an´ı ˇsumu pomoc´ı diskr´etn´ı konvoluce se rozmaz´an´ı obr´azku pravdˇepodobnˇe ˇ ım vˇetˇs´ı ˇsum potˇrebujeme z obr´azku odfiltrovat, t´ım v´ıce budeme muset nevyhneme. C´ obr´azek rozmazat.
17
Obr´azek 16: Lena pˇred odstranˇen´ım ˇsumu (vlevo) a Lena po odstranˇen´ı ˇsumu (vpravo) Velmi podobn´y v´ysledek bychom obdrˇzeli i pˇri pouˇzit´ı Gaussovy masky 1 4 7 4 1 4 16 26 16 4 1 7 26 41 26 7 · , 273 4 16 26 16 4 1 4 7 4 1 viz obr. 17.
Obr´azek 17: Lena pˇred odstranˇen´ım ˇsumu (vlevo) a Lena po odstranˇen´ı ˇsumu pomoc´ı Gaussovy masky (vpravo) Dalˇs´ı informace o vyuˇzit´ı konvoluce pro zpracov´an´ı obrazu m˚ uˇze zv´ıdav´y ˇcten´aˇr z´ıskat v textu [3].
18
Reference [1] T. Kozubek, T. Brzobohat´y, V. Hapla, M. Jaroˇsov´a, A. Markopoulos: Line´arn´ı algebra s Matlabem. Text vytvoˇren´y pˇri realizaci projektu Matematika pro ” inˇzen´yry 21. stolet´ı“, Vysok´a ˇskola b´an ˇ sk´a - Technick´a univerzita Ostrava (2012). (http://mi21.vsb.cz/modul/linearni-algebra-s-matlabem) [2] K. Sigmon: Matlab Primer. University of Florida (1993). [3] E. Sojka, J. Gaura, M. Krumnikl: Matematick´e z´ aklady digit´ aln´ıho zpracov´an´ı obrazu. Text vytvoˇren´y pˇri realizaci projektu Matematika pro inˇzen´yry ” 21. stolet´ı“, Vysok´a ˇskola b´an ˇ sk´a - Technick´a univerzita Ostrava (2012). (http://mi21.vsb.cz/modul/matematicke-zaklady-digitalniho-zpracovani-obrazu)
19
Apendix A: Matice Definice A.1 Necht’ jsou d´any prvky a1,1 , a1,2 , . . . , am,n z dan´e mnoˇziny F. Matice typu (m, n) (struˇcnˇe m × n matice) je obd´eln´ıkov´a tabulka " a . . . a1,n # 1,1 .. .. .. A= , . . . am,1 . . . am,n kter´a m´a m · n prvk˚ u ai,j uspoˇr´adan´ych do m ˇr´adk˚ u riA a n sloupc˚ u sA ze j , takˇ
r1A A A = ... = sA 1 , . . . , sn , A " a # rm 1,j .. A A ri = ai,1 , . . . , ai,n , sj = . . am,j Struˇcnˇe p´ıˇseme t´eˇz A = [ai,j ].
20
