Z´ aklady zpracov´ an´ı obrazu Tom´aˇs Mikolov, FIT VUT Brno V tomto cviˇcen´ı si uk´aˇzeme z´akladn´ı techniky pouˇz´ıvan´e pro digit´aln´ı zpracov´an´ı obrazu. Pro jednoduchost budeme pracovat s obr´azky ve stupn´ıch ˇsedi - zpracov´an´ı barevn´ ych RGB obr´azk˚ u prob´ıh´a analogicky, mus´ıme vˇsak vˇsechny operace prov´est pro kaˇzdou barevnou sloˇzku zvl´aˇsˇt.
1
Naˇ cten´ı a vykreslen´ı obr´ azku I=imread(’lena.gif’); imshow(I); Nezapom´ınejte stˇredn´ıky za pˇr´ıkazy, kter´e vrac´ı bitmapy (matice)!
1.1
Z´ akladn´ı operace s obrazem
Velikost obr´azku m˚ uˇzeme zjistit jednoduˇse pˇr´ıkazem size(I), kter´ y n´am vr´at´ı velikost matice. Pro v´ ypis vˇsech nastaven´ ych promˇenn´ ych lze pouˇz´ıt pˇr´ıkaz whos, kde uvid´ıme i datov´ y typ pˇriˇrazen´ y promˇenn´ ym n´aˇs obr´azek bude m´ıt typ uint8. Otoˇcen´ı obrazu: Irot=imrotate(I, 30, ’bilinear’); imshow(Irot); Uloˇzen´ı obrazu: imwrite (I, ’lenacopy.gif’);
´ Ukol • zjistˇete velikost obr´azku lena.gif pˇred a po otoˇcen´ı o 45 stupˇ n˚ u
2
Z´ akladn´ı techniky zpracov´ an´ı obrazu
Mnoho jednoduch´ ych transformac´ı obrazu m˚ uˇzeme realizovat pomoc´ı pr˚ uchodu pˇres celou matici obrazu, kde postupnˇe mˇen´ıme hodnoty jednotliv´ ych prvk˚ u podle nˇejak´e pˇredem dan´e funkce. Tyto operace lze ve velk´e vˇetˇsinˇe realizovat jedn´ım pˇr´ıkazem v Matlabu. Protoˇze vˇsak smyslem tohoto cviˇcen´ı je uk´azat, jak vˇeci funguj´ı, uk´aˇzeme si ”ruˇcn´ı” pˇr´ıstup, kter´ y je sice v´ yraznˇe m´enˇe efektivn´ı, ale v´ıce n´azorn´ y. Pokud se pod´ıv´ame na matici I, kterou jsme vytvoˇrili v pˇredch´azej´ıc´ım pˇr´ıkladu naˇcten´ım obr´azku lena.gif, zjist´ıme, ˇze se n´am v n´ı objevuj´ı hodnoty v rozsahu 0-255. Samotn´a hodnota urˇcuje intenzitu b´ıl´e barvy na dan´e pozici v obraze (0 odpov´ıd´a ˇcern´e, 255 b´ıl´e barvˇe). N´azornˇe to ukazuje funkce colorbar. Pozice [0, 0] odpov´ıd´a lev´emu horn´ımu rohu - m˚ uˇzeme si vyzkouˇset jednoduch´ y experiment: Icopy=I; Icopy(1:200, 1:400)=zeros(200, 400); imshow(Icopy); colorbar;
1
2.1
Inverze barev
Pomˇernˇe zn´amou transformac´ı obrazu je inverze barev: for (x=1:512) for (y=1:512) I2(x,y)=255-I(x,y); end; end; imshow(I2); V Matlabu by staˇcilo napsat: I2 = 255 - I; imshow(I2); V tomto cviˇcen´ı V´am ale schv´alnˇe budeme vˇetˇsinu operac´ı presentovat pomoc´ı cykl˚ u (i kdyˇz je to “Matlabovˇe neoptim´aln´ı”), abyste si dok´azali pˇredstavit, jak tyto operace implementovat v C, Javˇe, atd.
2.2
Jednoduch´ a detekce hran
Pro dalˇs´ı zpracov´an´ı obrazu je ˇcasto vhodn´e detekovat hrany. Zat´ım si uk´aˇzeme velmi jednoduch´ y pˇr´ıstup, kdy zobrazujeme rozd´ıl hodnot sousedn´ıch pixel˚ u n´asoben´ y nˇejakou konstantou (v naˇsem pˇr´ıpadˇe 10) pro vetˇs´ı kontrast: for (x=1:511) for (y=1:512) I2(x,y)=(I(x,y)/2 - I(x+1,y)/2)*10; end; end; imshow(I2); Protoˇze na hran´ach je rozd´ıl hodnot sousedn´ıch pixel˚ u mnohem vˇetˇs´ı neˇz na jednobarevn´ ych ploch´ach, zobraz´ı se n´am hrany jako b´ıl´e ˇc´ary.
2.3
Thresholding
Mezi dalˇs´ı klasick´e funkce patˇr´ı thresholding (ˇcesky prahov´an´ı). Nejprve si urˇc´ıme hodnotu prahu a pot´e pixely obrazu rozdˇel´ıme podle t´eto hodnoty do dvou tˇr´ıd. prah=120; for (x=1:512) for (y=1:512) if (I(x,y)<prah) I2(x,y)=0; else I2(x,y)=255; end; end; end; imshow(I2);
´ Ukol • vyzkouˇsejte r˚ uzn´e hodnoty pro prahovou hodnotu; jak´a z nich je nejlepˇs´ı a lze ji spoˇc´ıtat?
2
3
Histogram
Histogram je funkce, kter´a n´am urˇcuje poˇcet pixel˚ u dan´e barvy. V naˇsem pˇr´ıpadˇe m´ame 256 r˚ uzn´ ych barev (0-255) - pokud tedy spoˇc´ıt´ame, kolikr´at se kter´a barva objevuje v naˇsem obr´azku a graf vykresl´ıme, dostaneme histogram: imhist(I); “Ruˇcn´ı” v´ ypoˇcet by vypadal takto: Pocet=zeros(256); for (x=1:512) for (y=1:512) Pocet(I(x,y))=Pocet(I(x,y))+1; end; end; plot(Pocet); Z histogramu je patrn´e, ˇze nˇekter´e hodnoty barev nebyly v naˇsem obr´azku pouˇzity v˚ ubec, zat´ımco jin´e byly pouˇzity mnohokr´at. M˚ uˇzeme tedy vˇsechny barvy rozt´ahnout tak, aby histogram opravdu vyuˇz´ıval vˇsechny hodnoty 0-255: I2=histeq(I); subplot(221); imshow(I); title (’original’); subplot(222); imhist(I); subplot(223); imshow(I2); title (’equalized’); subplot(224); imhist(I2); Po tomto pˇr´ıkladu subploty opˇet vypnˇete pomoc´ı subplot(111); nebo close all
4
Line´ arn´ı filtry
Jak uˇz jsme si uk´azali v pˇredchoz´ım pˇr´ıkladu, m˚ uˇzeme detekovat hrany tak, ˇze od sebe odeˇcteme hodnoty sousedn´ıch pixel˚ u. Opaˇcn´a funkce, tedy seˇcen´ı hodnot sousedn´ıch pixel˚ u, n´am zp˚ usob´ı rozmaz´an´ı obrazu. Funkce, kter´e pracuj´ı s line´arn´ı kombinac´ı okol´ı pixelu pro vyprodukov´an´ı nov´e hodnoty, lze realizovat pomoc´ı line´arn´ıch filtr˚ u.
4.1
Blur (rozmaz´ an´ı obrazu)
Nejprve si uk´aˇzeme filtrov´an´ı pomoc´ı pˇr´ıkazu imfilter: h = ones(5,5) / 25 I2 = imfilter(I,h); imshow(I2); Matice h n´am zde pˇredstavuje filtr. V tomto pˇr´ıpadˇe tedy posˇc´ıt´a vˇsechny sousedn´ı pixely do vzd´alenosti 2 a vyprodukuje pr˚ umˇernou hodnotu. Abychom l´epe pochopili, jak vˇse funguje, uk´aˇzeme si ruˇcn´ı v´ ypoˇcet: h=ones(5,5)/25 I3=zeros(512,512); for (x=3:509) for (y=3:509) for (x2=1:5) for (y2=1:5) I3(x, y)=I3(x, y)+h(x2, y2)*double(I(x+x2-3, y+y2-3)); end; end; end; end; imshow(uint8(I3)); Proch´az´ıme tedy cel´ y obraz a na kaˇzd´e pozici vypoˇc´ıt´ame novou hodnotu v obrazu tak, ˇze vyn´asob´ıme j´adro filtru a odpov´ıdaj´ıc´ı ˇca´st obrazu. 3
4.2
Zaostˇ ren´ı obrazu
Opaˇcnou funkc´ı k rozmaz´an´ı je zaostˇren´ı obrazu. Tohoto efektu dos´ahneme tak, ˇze od hodnoty pixelu odeˇcteme hodnoty sousedn´ıch pixel˚ u. h2=[-1 -1 -1; -1 9 -1; -1 -1 -1] I2 = imfilter(I,h2); imshow(I2);
4.3
Detekce hran
Jednoduchou detekci hran jsme si jiˇz uk´azali - odpov´ıdal by j´ı filtr [5 -5]. Probl´emem takov´eho pˇr´ıstupu je znaˇcn´e zaˇsumˇen´ı, protoˇze nen´ı uvaˇzov´ano v˚ ubec okol´ı pixelu. D´ale takov´ y detektor hran nen´ı schopen rozpoznat horizont´aln´ı hrany. Abychom probl´emy odstranili, pouˇzijeme vˇetˇs´ı filtr a to dvakr´at: jednou pro detekci vertik´aln´ıch hran a jednou pro detekci horizont´aln´ıch hran. h=[1 0 -1; 2 0 -2; 1 0 -1] I2 = imfilter(I,h); h2=h’ I3 = imfilter(I,h2); imshow(I2/2+I3/2);
5
Komprese ` a la JPEG pomoc´ı DCT
Zˇrejmˇe kaˇzd´ y uˇzivatel poˇc´ıtaˇce se uˇz v dneˇsn´ı dobˇe setkal s JPEG kompres´ı obrazu. V tomto pˇr´ıkladu si uk´aˇzeme jej´ı podstatu. Nejprve rozdˇel´ıme obraz na bloky 8x8 pixel˚ u a na ty aplikujeme 2D diskr´etn´ı kosinovou transformaci (DCT). Tato transformace n´am umoˇzn´ı oddˇelit sloˇzky obrazu s r˚ uznou frekvenc´ı. Pot´e ”zahod´ıme” vysokofrekvenˇcn´ı sloˇzky obrazu, d´ıky ˇcemuˇz dos´ahneme komprese. N´aslednˇe zpˇetnou transformac´ı dostaneme opˇet bloky 8x8 pixel˚ u, ze kter´ ych zrekonstruujeme obraz. Aˇckoliv je patrn´a ztr´ata informace ve vysok´ ych frekvenc´ıch (hrany jsou rozmazan´e), kvalita obrazu je poˇra´d pomˇernˇe dobr´a vzhledem k tomu, ˇze jsme zahodili 58 z 64 koeficient˚ u = 90.625% informace. Id = im2double(I); T = dctmtx(8) dct = @(x)T * x * T’; % aplikovat na kazdy 8x8 blok B = blkproc(Id,[8 8],dct); mask = [1 1 1 0 0 0 1 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 B2 = blkproc(B,[8 8],@(x)mask.* invdct = @(x)T’ * x * T; I2 = blkproc(B2,[8 8],invdct); imshow(I), figure, imshow(I2)
0 0 0 0 0 0 0 0 x);
% definuje matici DCT bazi % definuje funkci, kteou budeme % % % %
prima DCT tady to udelame maska, ktera vybere jen prvnich par koeficientu ...
0 0 0 0 0 0 0 0]; % vsechny bloky 8x8 se vymaskuji % definice zpetne DCT % kterou zde aplikujeme.
Ve skuteˇcnosti je JPEG komprese samozˇrejmˇe komplikovanˇejˇs´ı - m´enˇe podstatn´a informace nen´ı zahozena, ale uloˇzena s menˇs´ı pˇresnost´ı (na menˇs´ı poˇcet bit˚ u). To sam´e plat´ı o barevn´ ych sloˇzk´ach obrazu, protoˇze lidsk´e oko je na r˚ uzn´e barvy jinak citliv´e. Pˇr´ıklad z http: // www. mathworks. com/ access/ helpdesk/ help/ toolbox/ images/ f21-16366. html 4
´ Ukoly 1. pod´ıvejte se na ˇr´adky a sloupce matice T, kter´a obsahuje DCT b´aze, napˇr. takto: plot plot plot plot
(1:8, (1:8, (1:8, (1:8,
T(1:3,:)’) T’) T(:,1:3)) T)
% zobrazi prvni 3 radky % zobrazi vsechny radky % zobrazi prvni 3 sloupce % zobrazi vsechny sloupce
Komentujte, kter´e frekvence kter´a b´aze “hled´a”. 2. zkuste ruˇcnˇe zkonstruovat n´asleduj´ıc´ı obr´azky: • • • • •
cel´ y ˇcern´ y. cel´ y b´ıl´ y. s vodorovn´ ymi prouˇzky. se svil´ ymi prouˇzky. “ˇsachovnici” po jednoitliv´ ych pixelech.
Obr´azky proˇzeˇ nte uveden´ ym algoritmem a • Zobrazte a komentujte obsah matice B, zobrazte napˇr. B (1:8, 1:8) B2 (1:8, 1:8) • Komentujte, jak kvalitn´ı je v´ ystup. ˇ sen´ı pomoc´ı funkce blockproc je velmi efektivn´ı, ale nevid´ıme pˇresnˇe, co se dˇeje. Zkuste pˇr´ıklad 3. Reˇ pˇrepsat pomoc´ı cykl˚ u a “ruˇcnˇe” zpracovat kaˇzd´ y blok 8x8.
5
6
Chyba v obraze
Prov´ad´ıme-li pokusy s obrazem jako napˇr´ıklad v pˇredchoz´ım pˇr´ıkladu, je uˇziteˇcn´e umˇet spoˇc´ıtat, jak´e chyby jsme se dopustili. To m˚ uˇzeme prov´est napˇr´ıklad tak, ˇze spoˇc´ıt´ame pr˚ umˇernou chybu na jeden pixel: noise=0; I2=im2uint8(I2); for (x=1:512) for (y=1:512) noise=noise+double(abs(I(x,y)-I2(x,y))); end; end; noise=noise/512/512
7
Reference
http://www.mathworks.com/access/helpdesk/help/toolbox/images/f0-3373.html http://www.fit.vutbr.cz/study/courses/ISS/public/pred/2d/aux.m
8
ˇ sen´ı Reˇ
Ruˇcn´ı zpracov´an´ı DCT:
6
Id = im2double(I); % prima dct T = dctmtx(8) B = zeros (512,512); for x=1:8:511, for y=1:8:511, vysek = Id (x:(x+7),y:(y+7)); dctvyseku = T * vysek * T’; B(x:(x+7),y:(y+7)) = dctvyseku; end end % maskovani mask = [1 1 1 0 0 0 0 0 % maska, ktera vybere jen 1 1 0 0 0 0 0 0 % prvnich par koeficientu ... 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]; B2 = zeros (512,512); for x=1:8:511, for y=1:8:511, vysek = B (x:(x+7),y:(y+7)); vymaskovano = vysek .* mask; B2(x:(x+7),y:(y+7)) = vymaskovano; end end % zpetna DCT I2 = zeros (512,512); for x=1:8:511, for y=1:8:511, vysek = B2 (x:(x+7),y:(y+7)); zpet = T’ * vysek * T;; I2(x:(x+7),y:(y+7)) = zpet; end end imshow(I), figure, imshow(I2)
7