Stereoseende i realtid

(1)

Institutionen för systemteknik

Department of Electrical Engineering

Examensarbete

Stereoseende i realtid

Examensarbete utfört i Bildbehandling vid Tekniska högskolan i Linköping

av Lars Arvidsson LITH-ISY-EX--07/3944--SE

Linköping 2007

Department of Electrical Engineering Linköpings tekniska högskola

Linköpings universitet Linköpings universitet

(2)

(3)

Stereoseende i realtid

Examensarbete utfört i Bildbehandling

vid Tekniska högskolan i Linköping

av

Lars Arvidsson LITH-ISY-EX--07/3944--SE

Handledare: Johan Skoglund

isy, Linköpings universitet Examinator: Michael Felsberg

isy, Linköpings universitet

(4)

(5)

Avdelning, Institution

Division, Department

Computer Vision Laboratory Department of Electrical Engineering Linköpings universitet

SE-581 83 Linköping, Sweden

Datum Date 2007-06-08 Språk Language Svenska/Swedish Engelska/English Rapporttyp Report category Licentiatavhandling Examensarbete C-uppsats D-uppsats Övrig rapport

URL för elektronisk version

http://www.isy.liu.se http://www.ep.liu.se/2007/3944 ISBN — ISRN LITH-ISY-EX--07/3944--SE

Serietitel och serienummer

Title of series, numbering

ISSN

—

Titel

Title

Stereoseende i realtid Stereo Vision in Real-Time

Författare

Author

Lars Arvidsson

Sammanfattning

Abstract

In this thesis, two real-time stereo methods have been implemented and evaluated. The first one is based on blockmatching and the second one is based on local phase. The goal was to be able to run the algorithms at real-time and examine which one is best. The blockmatching method performed better than the phase based method, both in speed and accuracy. SIMD operations (Single Instruction Multiple Data) have been used in the processor giving a speed boost by a factor of two.

Nyckelord

Keywords stereo, disparitet, disparity, real-time, realtid, SIMD, fas, phase, blockmatchning, blockmatching

(6)

(7)

Abstract

In this thesis, two real-time stereo methods have been implemented and evaluated. The first one is based on blockmatching and the second one is based on local phase. The goal was to be able to run the algorithms at real-time and examine which one is best. The blockmatching method performed better than the phase based method, both in speed and accuracy. SIMD operations (Single Instruction Multiple Data) have been used in the processor giving a speed boost by a factor of two.

Sammanfattning

I det här exjobbet har två stereometoder för realtidstillämpningar implementerats och utvärderats. Den ena bygger på blockmatchning och den andra på lokal fas. Målet var att kunna köra metoderna i realtid och undersöka vilken av dem som fungerar bäst. Blockmatchningsmetoden gav gott resultat medan den fasbasera-de fungerafasbasera-de sämre, båfasbasera-de vad gäller hastighet och precision. SIMD-operationer (Single Instruction Multiple Data) användes hos processorn vilket resulterade en i fördubbling av prestandan.

(8)

(9)

Tack

Jag vill tacka all personal på Computer Vision Laboratory som varit med och hjälpt till på alla sätt. Speciellt vill jag tacka min handlare Johan Skoglund och min examinator Michael Felsberg.

Jag vill också tacka min familj som alltid stöttar mig.

(10)

(11)

Innehåll

1 Inledning 1 1.1 Bakgrund . . . 1 1.2 Stereoproblemet . . . 1 1.3 Systemöversikt . . . 2 2 Teori 5 2.1 Nålhålskamera . . . 5 2.2 Triangulering . . . 6 2.3 Kamerakalibrering . . . 6 2.3.1 Sammanfattning . . . 7

2.3.2 Initiering av de intrinsiska och extrinsiska parametrarna . . 8

2.3.3 Linsdistorsion . . . 9 2.3.4 Optimering av kalibreringsparametrarna . . . 10 2.4 Rektifiering . . . 10 2.5 Stereoalgoritmer . . . 12 2.6 Blockmatchning . . . 12 2.6.1 Försteg . . . 12 2.6.2 Disparitetsberäkning . . . 12 2.6.3 Robusthet . . . 13 2.6.4 Subpixel-interpolation . . . 15 2.7 Fasbaserad stereo . . . 16 2.7.1 Skattning av fasskillnad . . . 18 2.7.2 Skattning av frekvens . . . 18 2.7.3 Säkerhetsvärden . . . 21 3 Implementation 25 3.1 Utrustning . . . 25 3.2 Blockmatchning . . . 27 3.2.1 Försteg . . . 27 3.2.2 Disparitetsberäkning . . . 28 3.2.3 SIMD-optimering . . . 30 3.3 Fasbaserad stereo . . . 30 ix

(12)

4 Utvärdering 33

4.1 Kvalitet och hastighet . . . 34

4.2 Resultat . . . 34

4.3 Metodegenskaper . . . 37

4.3.1 Parametrar hos den fasbaserade metoden. . . 39

5 Diskussion 43 5.1 Slutsats . . . 43

5.2 Svårigheter och problem . . . 43

5.3 Framtiden . . . 44

(13)

Kapitel 1

Inledning

1.1 Bakgrund

Det här examensarbetet är en del av forskningsprojektet COSPAL (COgniti-ve Systems using Perception-Action Learning) på bildbehandlingsavdelningen på ISY. COSPAL-projektet har som mål att undersöka designprinciper och arkitek-turer för tekniska system som kan lära in grundläggande kognitiva egenskaper på ett liknande sätt som vi människor gör. I projektet vill man utveckla djupskatt-ning från stereo med denna lärande metod. För att kunna jämföra dem med de befintliga metoderna vill man ha dessa implementerade i sitt system.

Det här exjobbet har som mål att implementera och utvärdera två metoder för att skatta djup från stereobilder i realtid. Den ena metoden är blockmatchning och baseras på [12]. Den andra metoden skattar djup utifrån fas och baseras på [13]. Målet har varit att skapa ett helt fungerande system som hämtar in bilder från två firewire-kameror, utför beräkningar och visar den resulterande djupkartan.

Båda metoderna har implementerats i C++ och utvärderats med testbilder från Middlebury [9], egna testbilder och verkliga bilder från två firewire-kameror. Målet har varit att få de så snabba som möjligt med god noggrannhet. Riktmärket har varit att kunna generera en djupkarta fem gånger per sekund för bilder med storleken 800x600 pixlar.

1.2 Stereoproblemet

När vi tittar med våra ögon får vi två olika bilder från vänster och höger öga. Bilderna liknar varandra men är förskjutna i sidled. Objekt som befinner sig nära ögonen är förskjutna mer än objekt långt borta. Genom att använda kameror kan man uppskatta djupet till objekt i bilderna genom att försöka skatta dessa förskjutningar som kallas dispariteter. Se figur 1.1 för ett exempel.

Det finns många tillämpningar där man har nytta av att veta djupet till objekt man ser. T ex så kan en robot lokalisera sig i sin omgivning med hjälp av två kameror. Forskning har pågått i många år inom området och det finns en uppsjö

(14)

2 Inledning

(a) Vänster bild. (b) Höger bild.

(c) Korrekt disparitet.

Figur 1.1. Bildparet tsukuba med tillhörande korrekta disparitet [9]. Objekt som befin-ner sig nära kameran är ljusare i djupkartan än objekt längre bort.

litteratur och artiklar. I takt med att tekniken utvecklats och datorerna blivit snabbare har intresset för beräkningar i realtid ökat. Sedan Zang presenterade sin enkla metod för kamerakalibrering har tillgängligheten ökat betydligt för gemene man. Det är i praktiken möjligt att få ett väl fungerande stereosystem för en väldigt liten penning.

1.3 Systemöversikt

Systemet för stereoseende visas i figur 1.2. Först hämtas bilderna in från kameror-na och rektifieras. Rektifieringen utförs mha av kalibreringsdata som skattats i ett kalibreringssteg som utförs offline. Sedan följer disparitetsberäkningen med anting-en blockmatchningsmetodanting-en eller danting-en fasbaserade metodanting-en. Till sist trianguleras en 3D-modell fram mha den slutliga disparitetskartan och kalibreringsdata.

(15)

1.3 Systemöversikt 3

Rektifiering

Kalibrering (offline)

Disparitetsskattning Triangulering 3D-data höger bild

vänster bild

djupkarta

(16)

(17)

Kapitel 2

Teori

2.1 Nålhålskamera

En kamera modelleras ofta som en nålhålskamera med ett optiskt centrum, C, och ett bildplan Π. Fokalavståndet f är det vinkelräta avståndet mellan bildplanet och fokalpunkten. En godtycklig punkt i världen, W , projiceras i punkten M i bildplanet som är skärningspunkten mellan vektorn−−→CW och bildplanet Π.

C

M

W Π

f

Figur 2.1. Modell för nålhålskameran.

Projektionen från 3D-världen till 2D-bilden beskrivs av de extrinsiska och in-trinsiska parametrarna [10]. De exin-trinsiska parametrarna talar om hur världsko-ordinaterna transformeras till kamerakoordinater och består av en rotation, R, och en translation, t. De intrinsiska parametrarna beskriver hur en punkt i ka-merakoordinater projiceras på bildplanet och består av skalfaktorerna α och β, skevningen c och principalpunkten (centrumpunkten i bilden) (u0, v0). Vi kan

be-skriva projektionen från punkten W = (X, Y, Z) till punkten M = (u, v) med hjälp 5

(18)

6 Teori

av homogena koordinater enligt

s   u v 1  = AR t     X Y Z 1     (2.1) där A =   α c u0 0 β v0 0 0 1  , R =   r11 r12 r13 r21 r22 r23 r31 r32 r33  , t =   t1 t2 t3   (2.2) och s är en skalfaktor.

2.2 Triangulering

För att kunna bestämma positionen för punkten W utifrån dess bildprojektion krävs i allmänhet två kameror. Punkten W bestäms då av skärningen mellan vektorn−−−→ClMl och

−−−→

CrMr. Se figur 2.2.

Vi behöver alltså känna till positionerna för vänster och höger bildplan, Πl

och Πr, fokalpunkterna fl och fr och punkterna Ml och Mr för att kunna göra

en fullständig triangulering. För att få fram en tät djupkarta, dvs beräkna djup för varje punkt i bilden, behöver vi känna till ett punktpar (Ml, Mr) för varje

enskild punkt i scenen vi är intresserad av. Det är här den stora utmaningen i stereoproblemet ligger – att hitta dessa punktpar.

C

_l

M

_l

Π

l

f

_l

C

r

M

_r

Π

r

f

r

W

Figur 2.2. Triangulering med två nålhålskamereror.

2.3 Kamerakalibrering

Att skatta de extrinsiska och intrinsiska parametrarna för en kamera kallas för kalibrering. Det finns flera sätt att göra det på och i det här projektet har Matlab-verktyget Caltech Calibration Toolkit [1] använts. Den bygger på en metod som

(19)

2.3 Kamerakalibrering 7

Zhang [14] utvecklat och går ut på att man fotograferar ett plant rutmönster i en mängd olika positioner och sedan parar ihop kalibreringspunkterna i bilden med motsvarande punkter i rutmönstret. Givet dessa korresponderande punkter kan man sedan skatta de sökta parametrarna. Metoden är enkel att utföra och ger gott resultat. Nedan sammanfattas de olika stegen i kalibreringen och följs sedan av en beskrivning av de intressantaste delarna.

2.3.1 Sammanfattning

Zhangs kamerakalibreringsmetod:

1. Skriv ut ett rutmönster och fäst det på ett plant underlag.

2. Fotografera planet i några olika lägen. Ca 20 stycken rekommenderas. 3. Detektera kalibreringspunkterna i bilderna. Se figur 2.3.

4. Skatta de intrinsiska och extrinsiska parametrarna med en metod på sluten form.

5. Inför linsdistorsion i kameramodellen och förfina skattningen genom att mi-nimera det totala projektionsfelet.

I Caltech Calibration Toolkit används en modifierad variant av Zhangs metod där man även tar hänsyn till den tangetiella linsdistorsionen och inte bara den radiella.

(20)

8 Teori

2.3.2 Initiering av de intrinsiska och extrinsiska

paramet-rarna

En punkt M = (X, Y, Z) i rutmönstret projiceras i punkten m = (u, v) i bildplanet enligt (2.1). Vi kan välja att placera rutmönstrets koordinatsystem så att alla punkter i rutmönstret ligger i XY-planet (Z = 0). Då får vi att

s   u v 1  = A   r11 r12 r13 t1 r21 r22 r23 t2 r31 r32 r33 t3       X Y 0 1     = A   r11 r12 t1 r21 r22 t2 r31 r32 t3     X Y 1   (2.3) Vi sätter H =   h11 h12 h13 h21 h22 h23 h31 h32 h33  = A   r11 r12 t1 r21 r22 t2 r31 r32 t3   (2.4)

vilket ger oss att

s   u v 1  =   h11 h12 h13 h21 h22 h23 h31 h32 h33     X Y 1   (2.5)

Genom att eliminera s i (2.5) får vi att

X Y 1 0 0 0 −Xu −Y u −u

0 0 0 X Y 1 −Xv −Y v −v               h11 h12 h13 h21 h22 h23 h31 h32 h33               = 0 (2.6)

Vi har alltså två ekvationer och nio obekanta.

För att skatta h11. . . h33måste vi alltså ha fler ekvationer. Med n

kalibrerings-punkter får vi 2 × n ekvationer som vi staplar på varandra. Vi sätter

L =        X1 Y2 1 0 0 0 −X1u1 −Y1u1 −u1 0 0 0 X1 Y1 1 −X1v1 −Y1v1 −v1 .. . ... ... ... ... ... ... ... ... Xn Yn 1 0 0 0 −Xnun −Ynun −un 0 0 0 Xn Yn 1 −Xnvn −Ynvn −vn        (2.7)

(21)

2.3 Kamerakalibrering 9 och h =               h11 h12 h13 h21 h22 h23 h31 h32 h33               (2.8) vilket ger Lh = 0 (2.9)

Lösningen till ekvationen ges av egenvektorn till LTL som hör till dess lägsta egenvärdet.

Utan att gå in i för mycket detaljer så kan vi härifrån lösa ut A, R och t från H med en teknik som Zang presenterar i sin rapport [14].

2.3.3 Linsdistorsion

Hittills har linsdistorsion inte tagits hänsyn till i kameramodellen. Det visar sig dock nödvändigt för att kunna göra en noggrann kalibrering.

Ur geometrin i figur 2.4 får vi att ˜u ˜ v =f z x y (2.10) där ˜u och ˜v är bildkoordinaterna.

För att översätta (˜u, ˜v) till pixelkoordinater (u0, v0) behövs pixeltäthetskoeffi-cienterna Du, Dv och skalfaktorn su. Då blir

u0 v0 =Dusuu˜ Dv˜v +u0 v0 (2.11) där (u0, v0) är centrumpunkten i bilden. (x, y, z) ˆ z ˆ u, ˆx (0, 0, f ) bildplan (˜u, ˜v)

(22)

10 Teori

I Calib Calibration Toolkit modelleras linsdistorsionen av en radiell och en tangentiell del [5]. Den radiella beskrivs av

δur δvr = ˜u(k1r 2_{+ k} 2r4+ k3r6) ˜ v(k1r2+ k2r4+ k3r6) (2.12)

där k1, k2 och k3är koefficienter för den radiella distorsion och r =p ˜u2+ ˜v2).

Det tangentiella felet modelleras som δut δvt =2p1u˜˜v + p2(r 2_{+ 2˜}_u2₎ p1(r2+ 2˜v2) + 2p2u˜˜v (2.13) där p1 och p2är koefficienterna för den tangentiella distorsionen.

Den modifierade projiceringen från kamerakoordinter till pixelkoordinter blir efter att tagit hänsyn till linsdistorsionen

u v =Dusu(˜u + δur+ δut) Dv(˜v + δvr+ δvt) +u0 v0 (2.14)

2.3.4 Optimering av kalibreringsparametrarna

Efter att linsdistorsion införts i vår kameramodell kan vi göra en ny bättre skatt-ning av kalibreringsparametrarna genom att minimera det totala projektionsfelet. Vi låter mij vara den uppmätta kalibreringspunkten j i bild i och Mj dess

tillhörande kalibreringspunkt i rutmönstret. Vidare är ˆmij den beräknade

projek-tionen av Mj i bildplanet. Givet n olika bilder med m stycken kalibreringspunkter

i varje bild vill vi minimimera följande uttryck

n X i=1 m X j=1 kmij− ˆm(A, Ri, ti, Mj, k, p)k2 (2.15)

där k = [k1k2k3]T och p = [p1 p2]T är parametrarna för de radiealla respektiva

tangetiella linsdistorsionerna. Minimeringen är ett olinjärt problem som exempel-vis löses med Levenberg-Marquardt-algoritmen.

2.4 Rektifiering

De flesta stereoalgoritmer förutsätter att de inkommande bildparen är rektifierade. Att bilderna är rektifierade innebär att de är vridna så att deras epilolarlinjer är horisontella till varandra. Detta får till följd att de korresponderande punkterna befinner sig på samma horisontella linje i de båda bilderna. Se figur 2.5 för en förklarande figur.

Resultaten från användningen av Calib Calibration Toolkit anger hur bilderna ska vridas så att bilderna blir rektifierade.

(23)

2.4 Rektifiering 11 Cl Πl Πl,rektifierad Cr Πr Πr,rektiferad W

(24)

12 Teori

2.5 Stereoalgoritmer

Metoderna som löser problemet att hitta korresponderande punkter i vänster- och högerbild kan grovt delas in i areabaserade och egenkapasbaserade.

De areabaserade metoderna löser korresponderingsproblemet för varje enskild pixel i bilden. Vanligtvis söker man var omgivningen till pixeln matchar bäst till ett block i den andra bilden.

De egenskapsbaserade försöker först hitta egenskaper så som hörn, kanter och objekt i den ena bilden och sedan hitta motsvarande punkter i den andra. Resul-tatet från dessa metoder är en disparitetskarta med värden i just dessa punkter.

Det finns några problem som är gemensamma för alla stereoalgoritmer. Ett objekt som syns i den ena bilden kan vara skymd i den andra. Det resulterar i att dessa punkter kopplas ihop med fel punkter. Ett annat problem kan vara att linsen i kameran inför fel i bilden, vilket oftast syns i kanten på bilden. Bilderna själva kan också ställa till problem genom att vissa ytor kan sakna struktur och därmed blir svåra att hantera. Likaså stöter man på svårigheter om bilden innehåller ett mönster som upprepar sig.

I efterföljande avsnitt beskrivs de två metoder som implementerats.

2.6 Blockmatchning

Blockmatchningsmetoden består av fyra steg: 1. Försteg

2. Disparitetsskattning 3. Robusthet

4. Subpixel-interpolation

2.6.1 Försteg

I försteget normaliseras de rektifierade bilderna genom att subtrahera medelvärdet av intensitetsvärdena i en liten omgivning till varje pixel. Det görs för att kompen-sera för olika ljusförhållanden som kan förekomma mellan höger och vänster bild. I försteget beräknas också variansen på intensitetsvärdena i samma omgivning. Det ger ett mått på hur väl området är texturerat.

2.6.2 Disparitetsberäkning

Blockmatchningsmetoden försöker bestämma dispariteten för en pixel genom att hitta bästa matchningen mellan block i vänster bild och block i höger bild. I och med att bilderna är rektifierade räcker det att söka i horisontell led.

Vi har att den vänstra bilden, L, är referensbilden, R, högerbilden och d ∈ [0...dmax] dispariteten. Vi börjar i punkt L(u, v) i vänster bild och söker efter

(25)

2.6 Blockmatchning 13

[R(u − dmax, v) . . . R(u, v)] i höger bild. Dispariteten ges av det d som minimerar

kostnaden . Se figur 2.6.

L(u, v) Referensbild

Vänster bild, L Höger bild, R

R(u,v)

R(u + dmax, v)

L(u, v) Referensbild

Vänster bild, L Höger bild, R

R(u - dmax, v) Startläge

Slutläge

Figur 2.6. Blockmatchning från vänster.

Det finns några olika kostnadsfunktioner föreslagna i litteraturen [8]. Bland annat SSD (Sum of Square Differences) och SAD (Sum of Absolute Differences). I praktiken är skillnaden liten. SAD har fördelen att den går snabbare att beräk-na [11]. SSD däremot används vid beräkningen av kovariandmatrisen varför den ibland används som kostnadsfunktion om den ändå räknas ut.

Givet en vänsterbild L, en högerbild R och en blockstorlek på (2n + 1 × 2n + 1) beräknas SAD för ett block centrerat i (u, v) i vänsterbilden med dispariteten d som SAD(u, v, d) = n X j=−n n X i=−n

|L(u + j, v + i) − R(u − d + j, v + i)| (2.16) Denna summa beräknas för ett antal dispariteter från noll upp till dmax som

avgör söklängden. Dispariteten som svarat mot lägst SAD är troligtvis den sanna dispariteten. Genom att välja n och dmaxkan man styra slutresultatet på det sätt

man vill. I stora drag ger större värden på parametrarna bättre resultat till priset av längre beräkningstid. För varje pixel sparar vi de fyra lägsta värdena på SAD och vilka dispariteter de svarar mot. De har vi sedan nytta av i eftersteget.

2.6.3 Robusthet

Ett vanligt test man kan göra för att upptäcka felaktiga matchningar är att utföra ett vänster-höger-test. Det går till så att man beräknar dispariteten med vänster

(26)

14 Teori

bild som referens och sedan beräknar dispariteten fast nu med höger bild som referens. De punkter där dispariteten skiljer mellan höger- och vänsterbild marke-ras som felaktiga matchningar. Genom att utföra detta test får man ett säkrare slutresultat till priset av högre beräkningsbörda.

För att hålla beräkningsbördan så låg som möjligt har även en annan metod testats för att detektera felaktiga matchningar. Den kommer från [12] och har fördelen att man endast behöver beräkna dispariteten åt ett håll. Denna metod består av två test och kräver att vi hittar det lägsta SAD-värdet och tre lokala minima. Det ena testet är ett skarphetstest och det andra ett tydlighetstest där de fyra SAD-värdena jämförs. Se figur 2.7.

algorithm will reject the previous match and accept the new one. This implies that, although the proposed approach relies on a direct matching phase only, it allows for recovering from possible previous matching errors.

The capability of our algorithm to recover from previous errors as long as better matches are found during the search is also shown inFig. 2, which plots as a function of d[ ½0; …; dmax" all the scores between the point Rðx; yÞ of right image and the points in the reference image ½Lðx2dmax; yÞ; …; Lðx; yÞ" that are allowed to establish a correspondence with R_{ðx; yÞ (see}Fig. 3).

Recalling the arcs drawn inFig. 1, we can notice that smaller d values correspond to scores computed more recently while greater d values to scores computed earlier. Considering again the two matches (1) L_{ðx þ}a2 dmax; yÞ Q Rðx; yÞ and (2) Lðx þb2 dmax; yÞ Q Rðx; yÞ; the algorithm will discard the old one, (1) since the new one, (2) has a better score with Rðx; yÞ: Moreover, if we find a new ‘collision’ when analysing the successive points of the left image; (3) L_{ðx þ}g2 dmax; yÞ Q Rðx; yÞ; the score of this new match will be compared with that associated with the current best match for Rðx; yÞ; so as to retain only one single match. That is, as in Fig. 2, since 1_{ðx þ}g2 dmax; x; yÞ # 1ðx þb2 dmax; x; yÞ; Lðx þb2 dmax; yÞ Q Rðx; yÞ will be discarded and the current match for Rðx; yÞ set to Lðx þg2 dmax; yÞ Q Rðx; yÞ:

Fig. 4shows a geometric interpretation of the proposed matching approach: as long as the left-to-right matching process proceeds, the algorithm disambiguates between all the matches that imply locating 3D points lying on the same

line of sight originating from a point of the right image. These sort of matches, referred to as ‘colliding matches’ for a given point of the right image Rðx; yÞ; can be caused by the ambiguity of photometric cues as well as by occlusions. In

Fig. 4the former situation is represented by the incorrect match betweenaand x; the latter by the incorrect match betweenband x:

To conclude this section, we highlight one major difference between the SMP and BM approaches with the aid ofFig. 5that, for a given scanline, shows on the x-axis the x-coordinates of the points of the left image that are

Fig. 2. Scores associated with point Rðx; yÞ:

Fig. 3. Potential matches in the left image for Rðx; yÞ:

Fig. 4. Geometric interpretation.

Fig. 5. The search paths of bidirectional matching. L. Di Stefano et al. / Image and Vision Computing 22 (2004) 983–1005 985

Figur 2.7. SAD-värden för alla disparitetsvärden för en pixel i vänsterbilden.

Skarphetstestet ger ett mått på hur väl de tre beräknade dispariteterna för en pixel klumpar ihop sig. Värdet ∆d blir lägre för en skarpare disparitet och ges av

∆d =

3

X

i=1

|di− dmin| (2.17)

där dmin är dispariteten med lägst SAD och di är dispariteterna för de tre lokala

miniman. De tre lokala miniman är inte de absolut lägsta pga implementerings-tekniska skäl. Men de är tillräckligt bra för att användas i praktiken.

Tydlighetstestet avgör hur väl den lägsta av de tre SAD-värdena skiljer sig från de andra. Ju högre värde på ∆SAD desto troligare är det att disparitetsskattning-en är korrekt. Tydlighetsvärddisparitetsskattning-en ges av

∆SAD =

3

X

i=1

SADi− SADmin (2.18)

där SADmin är det minsta SAD och SADi är de tre lokala miniman.

Genom att tröskla ∆d och ∆SAD med lämpliga värden kan man avgöra huruvi-da man ska klassa dispariteterna som felaktiga matchningar eller inte. Vilka trös-kelvärden man ska välja får man pröva sig fram till.

(27)

2.6.4 Subpixel-interpolation

En nackdel med blockmatchningsmetoden är att den inte ger subpixel-noggrannhet. Utan det skulle systematiska fel kunna förekomma där den beräknade dispariteten tex är 5 pixlar, men den verkliga är 4,75. Då får vi ett fel på 0,25 pixlar. Det finns i huvudsak två sätt man kan lösa detta på:

• Interpolera fram subpixelblock som används i blockmatchningen.

• Anpassa en felfunktion till SAD-värdena kring den skattade dispariteten. Sök sedan minimum av den felfunktionen.

Den första metoden ger gott resultat men är långsam. Den andra beskrivs i Johan Skoglunds LIC-avhandling [11] och ger även den gott resultat men är snabbare. Eftersom hastighet är högsta prioritet har den senare valts.

Vi har den skattade dispariteten d samt SAD(d − 1), SAD(d) och SAD(d + 1) enligt figur 2.8. Dispariteten beräknas då enligt:

k1= SAD(d − 1) − SAD(d) (2.19) k2= SAD(d + 1) − SAD(d) (2.20) dsubpixel= SAD(d − 1) − SAD(d + 1) 2max(k1, k2) + d (2.21) !1 !0.8 !0.6 !0.4 !0.2 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 3.5 4 disparity SAD SAD(d ! 1) SAD(d + 1) SAD(d) d!subpixel Figur 2.8. Subpixel-interpolation

(28)

16 Teori

2.7 Fasbaserad stereo

Den fasbaserade metoden är utgår från [13] och är inte lika intuitiv och lättför-stådd som blockmatchningsmetoden. Idén är att man vill jämföra små områden i vänster och höger bild med hjälp av den lokala fasen. Man antar att bilder-na i små omgivningar är intrinsiskt endimensionella och har en entydig fas och frekvens. Genom att skatta fasen och frekvensen lokalt i båda bilderna kan man beräkna dispariteten.

Vi antar att en lokal omgivning i en bild kan beskrivas med följande funktion:

f (u, v) = A + cos(ωu + φ) (2.22)

där A är DC-nivån, ω vinkelfrekvensen i horisontell led och φ fasen.

Vi antar vidare att den lokala omgivningen kring en punkt i vänsterbilden är identisk i samma lokala omgivning i höger bild förutom en förskjutning i sidled. Dvs de har samma frekvens, ω, men olika fas, φL och φR. Se figur 2.9 och 2.10.

Relationen mellan de lokala omgivningarna i vänster och höger bild beskriver vi som fL(u, v) = fR(u − d, v) (2.23) där fL(u, v) = A + cos(ωu + φL) fR(u, v) = A + cos(ωu + φR) (2.24)

Genom att kombinera (2.22) och (2.23) får vi

cos(ωu + φL) = cos(ω(u − d) + φR) (2.25)

Om (2.25) ska gälla för alla u i den lokala omgivningen måste det gälla att

φL= φR− ωd (2.26) Ur (2.26) får vi slutligen att d = φR− φL ω = ∆φ ω (2.27)

Observera här att den beräknade dispariteten, d, inte behöver vara ett heltal även om ∆φ och ω är det. Det innebär att vi får subpixelnoggrannhet helt gratis till skillnad mot vad som var fallet i blockmatchningsmetoden.

Det som behöver skattas i den fasbaserade metoden är alltså: • Fasskillnaden mellan vänster och höger bild, ∆φ

(29)

2.7 Fasbaserad stereo 17

Figur 2.9. Exempel på en lokal omgivning som är lika i vänster och höger bild men förskjutna i sidled i relation till varandra.

Position 1 !2 Intensity ! 0 " "/2 #" #"/2 Phase Position ! 1 !2 $% $!

Figur 2.10. Genom att skatta den lokal fasen och frekvensen kan dispariteten bestäm-mas. Här är bilderna i figur 2.9 sedda i horisontell led. Den heldragna linjen representerar högerbilden och den streckade representerar vänsterbilden. (Bild från [3]).

(30)

18 Teori

2.7.1 Skattning av fasskillnad

Fasskillnaden skattar vi med ett kvadraturfilter. Ett kvadraturfilter är ett filter som är den analytiska signalen till ett bandpassfilter. Filtreringen kan utföras på två olika sätt i huvudsak; med multiplikation i Fouriersdomänen eller faltning i signaldomänen. I denna implementation har filtrering i Fourierdomänen undvikits. Det finns två anledningar till det:

• Filtrering i Fourierdomänen är cirkulärsymmetrisk och inför därmed artifak-ter i kanartifak-terna på bilderna.

• Fourier-transformeringingen kräver att vi Fouriertransformerar våra bilder vilket är beräkningskrävande.

Filtersvaren för vänster och höger bild ges av qL och qR och fås som

qL= BLeiφL

qR= BReiφR

(2.28)

där BL och BR är amplituderna på filtersvaren och φL och φR respektive fas.

Multiplicerar vi qL med konjugatet till qR får vi

qLqR= BLBRei(φL−φR)= BLBRei∆φ (2.29)

Efter det kan vi lösa ut ∆φ som

∆φ = arg(qLqR) (2.30)

I och med att vi söker fasskillnaden i horisontell led används ett endimensionellt kvadraturfilter. I figur 2.11 ser vi ett kvadraturfilter i både signal- och frekvens-domänen. I frekvensdomänen är filtret reellt och har en amplitud som är noll för negativa frekvenser och har en nollskild amplitud för positiva frekvenser. Detta innebär att filtret i signaldomänen består både av en reell och imaginär del där dessutom den reella delen är jämn och den imaginära delen är udda.

När den reella delen verkar på en signal kan det ses som ett linjefilter. Svaret från filtreringen blir stort och positivt när filtret verkar på en ljus linje och stort och negativt för en mörk linje. Filtrerar vi något konstant blir svaret noll. På motsvarande sätt kan den imaginära delen kan ses som ett kantfilter. Se figur 2.12 för att se vilka filtersvar som fås för olika typer av insignaler.

2.7.2 Skattning av frekvens

I och med att dispariteten, d, mäts i pixlar så uttrycker vi frekvensen, ω, i fasför-ändring per pixel i horisontell led. Vi beräknar då ω som

ω = arg(q(u + 1, v)q(u, v)) (2.31)

(31)

2.7 Fasbaserad stereo 19 !20 !15 !10 !5 0 5 10 15 20 !0.25 !0.2 !0.15 !0.1 !0.05 0 0.05 0.1 0.15 0.2 0.25 realdel imaginärdel

(a) Kvadraturfiltret i signaldomänen.

!3 !2 !1 0 1 2 3 0 0.2 0.4 0.6 0.8 1

(b) Absolutbeloppet av kvadraturfiltret i frekvensdomänen.

(32)

20 Teori 20 40 60 80 100 120 140 20 40 60 80 100 120 140 0 0.5 1 20 40 60 80 100 120 140 !3.14 !1.57 0 1.57 3.14

Figur 2.12. Fasexempel. (Bild från [3]).

!

(33)

För att få ett så säkert estimat av ω som möjligt tittar vi på fasdifferensen åt vänster och höger i båda bilderna. På det sättet får vi fyra skattningar som vägs samman till en. De fyra olika differenserna är

fL−= qL(u − 1, v)qL(u, v)

fL+= qL(u, v)qL(u + 1, v)

fR−= qR(u − 1, v)qR(u, v)

fR+= qR(u, v)qR(u + 1, v)

(2.32)

och de vägs samman som

ω = arg(fL−+ fL++ fR−+ fR+) (2.33)

2.7.3 Säkerhetsvärden

Till varje punkt i bilden har vi nu skattat den lokala fasen och frekvensen, men vi är inte säkra på hur väl det stämmer överens med verkligheten. För att råda bot på denna ovisshet beräknar vi några olika säkerhetsvärden.

Säkerhetsvärde C1

Amplituden på kvadraturfiltersvaren ger oss det första säkerhetsvärdet C1. Det

ger ett mått på hur väl signalen passar filtret. C1 ges av

C1=pkqLqRk _2kq LqRk kqLk2+ kqRk2 γ (2.34) där γ styr hur mycket amplitudskillnaden mellan qLoch qRska bidra med. Enligt

[13] är 4 ett bra värde på γ som fåtts fram efter prövning. Säkerhetsvärde C2

Om bredden på kvadraturfiltret är för smalt i förhållande till den sanna dispa-riteten fås en felaktig skattning av fasen. Ju närmre den skattade frekvensen är nära sitt min- och maxvärde dvs ±π desto troligare är det att vi gjort en felaktig skattning. Därför bör vi undertrycka sådana värden. C2 ger oss ett mått på pass

säkra vi är på att filtret passar signalen och ges av C2= cos2 arg(qLqR) 2 (2.35) Säkerhetsvärde C3

Det sista säkerhetsvärdet försöker ge ett mått på hur lika de skattade frekvenserna i (2.32) är. Är den skattade frekvensen noll eller negativ vet vi säkert att den är felskattad eftersom den ska vara positiv enligt förutsättningarna. C3får vi som

C3= ( 1 4 _f L− kfL−k+ fL+ kfL+k+ fR− kfR−k+ fR+ kfR+k då ω > 0 0 f.ö. (2.36)

(34)

22 Teori

Slutligen multiplicerar vi ihop alla säkerhetsvärden till ett gemensamt C så att vi får

C = C1C2C3 (2.37)

Skalpyramid

För att kunna göra en säker skattning av dispariteten i olika typer av omgivningar måste man anpassa kvadraturfiltret och bilderna till varandra. Man kan göra det på två sätt:

• Variera storleken på filtret. • Variera storleken på bilderna.

Det senare har valt eftersom bildstorleken och därmed beräkningsbördan mins-kar för varje nivå i skalpyramiden. Dessutom kan man behålla ett och samma filter genom alla beräkningar. Först skattas dispariteten i den grövsta skalan. Sedan går man nedåt i skalpyramiden och skattar dispariteten för nästa finare nivå osv. Se figur 2.14 för en förklarande bild.

Spatial Cons. Shift D i s p a r i t y e s t i m a t i o n Shift Shift Shift Accu-mulator Spatial Cons. Edge detect. Left image Right image Edge detect.

Figur 2.14. Systemskiss för fasbaserad stereo. (Bild från [13]).

Lokalt skifte

Efter att dispariteten skattats för en viss nivå i skalpyramiden utförs ett lokalt skifte på bilderna ett steg ner i pyramiden. Vänster och höger bild skiftas in mot

(35)

varandra så att dispariteten minskas. Givet vänster bild L, höger bild R och den skattade dispariteten d så förklarar vi det matematisk som

L+(u, v) = L(u + 0.5d, v)

R+(u, v) = R(u − 0.5d, v)

(2.38) där index + betyder bilden ett steg djupare i pyramiden.

Kantdetektering

Det kan finnas disparitetsskattningar i bilden som är orimliga, så kallalade outliers. Dessa punkter befinnner sig sällan på kanter eller hörn. Därför skattas dispariteten även på kantbilder för att få ett säkrare estimat. Kantbilderna ges av styrkan på kadraturfiltersvaren som redan beräknats.

Dispariteten och säkerhetsvärdena skattas för kantbilderna och vägs sedan sam-man med de samma för gråskalebilderna. Vi beräknar den samsam-manvägda dispari-teten som

d = Cgdg+ Cede Cg+ Ce

(2.39) där Cg och Ce är säkerheten för gråskalebilden respektive kantbilden och dg och

deär den skattade dispariteten för gråskalebilden respektive kantbilden. Det

sam-manvägda säkerhetsvärdet beräknas som C = Cge iarg(pg )₂ _{+ C} eei arg(pe) 2 (2.40) där pg= qgLqgR och pe= qeLqeR. Disparitetsackkumulator

För varje nivå i skalpyramiden ackumuleras dispariteten och säkerheten till den nuvarande som dnytt= dgammalt+ d (2.41) Cnytt= p Cgammalt+ √ C 2 !2 (2.42) Robusthet

Vi införde kantdetektorer för att öka chansen till en korrekt disparitetsskattning, men det visar sig inte räcka till. Det kan finnas punkter som har väldigt låg säkerhet och orimlig disparitet. När vi utför det lokala skiftet inför dessa punkter fel nedåt i hierarkin som är svåra att bli av med. Ett sätt att reducera deras påverkan är att ersätta skattningen utifrån en liten omgivning kring punkten. Ett ny disparitet vägs samman där punkter nära och med hög säkerhet viktas högt. Ett filter som passar bra för denna tillämpning är ett gaussfilter. Den har en hög bula på mitten och avtar fort åt sidorna. Filtret kallar vi h som har följande utseende:

h(u) = 1 σ√2πe

−u2

(36)

24 Teori

där R är storleken på filtret och σ är standardavvikelsen.

Vi faltar filtret, h, med både säkerheten, C, och dispariteten, d.

m = h ∗ C (2.44)

v = h ∗ (Cd) (2.45)

Den nya säkerheten, Cnytt, och dispariteten, dnytt, får vi som

Cnytt= m (2.46)

dnytt=

v

m (2.47)

Denna operation kallas för normaliserad faltning och innebär att om filtret är pla-cerat vid en punkt med hög säkerhet så är disparitetsskattningen oberörd, men är säkerheten låg ändras skattningen till att gå mot medelvärdet av dess omgivning.

(37)

Kapitel 3

Implementation

3.1 Utrustning

Hårdvara

Eftersom målet var att kunna köra stereometoderna i realtid har snabbheten varit den kritiska punkten. Det finns flera olika sätt att få hög hastighet i sina program. Till och börja med gäller det att ha rätt förutsättningar i form av hårdvara. I det-ta projekt har en PC med en AMD XP 3800+ processor och 2 GB RAM-minne använts. Processorn stödjer en mängd SIMD-instruktioner (Single Instruction Mul-tiple Data) som kan användas för att göra flera beräkningar parallellt. Bilderna har hämtats in med hjälp av två AVT Marlin firewire-kameror.

Utvecklingsmiljö

En förutsättning för att kunna utnyttja hårdvaran till fullo är att använda rätt programmeringsspråk. Som bildbehandlare ligger ofta Matlab nära till hands ef-tersom det är lätt att arbete med och har många färdiga funktioner man kan tänkas vilja använda. När det gäller realtidstillämpningar är det dock för lång-samt. C++ har istället använts eftersom det är mycket snabbare och kan utnyttja SIMD-instruktionerna hos processorn. Man kan kombinera smidigheten hos Mat-lab med kraften hos C++ genom att utnyttja det så kallade MEX-gränssnittet. Det gör det möjligt att anropa kod skriven i C++ direkt i Matlab.

För att inte behöva skriva alla möjliga funktioner själv kan det vara värt att använda färdiga kodbibliotek. Ett lämpligt sådant är OpenCV (Open Computer Vision Library) [2] som är ett OpenSource-projekt startat av Intel som tillhanda-håller en mängd bildbehandlingsfunktioner skrivna i C.

Oavsett hur snabb hårdvara man än har så går det ofta att öka hastighe-ten ytterligare genom att optimera sin kod. Ett sätt är att aktivera de inbyggda funktionerna i kompilatorn som försöker snabba upp programmen. Eftersom den metoden är helt gratis kan det alltid vara värt att göra ett försök och se om det blir någon förbättring. Parametrar som använts i gcc-kompilatorn är

(38)

26 Implementation

-O3 -march=prescott -mfpmath=sse

och försöker optimera koden så att den går så snabbt det bara går. • -O3 - aktiverar den kraftigaste optimeringen.

• -march=prescott - kompilerar för just den valda processorn.

• -mfpmath=sse - använder SSE-instruktioner automatiskt när det är möjligt. Ett lite kostsammare sätt att optimera sina program är att försöka minimera de redundanta beräkningar på egen hand. Det kan vara arbetsamt men kan ge stora hastighetsökningar.

Slutligen så kan man försöka utnyttja parallellismen som kan finnas i meto-derna och beräkna flera saker samtidigt. Det finns flera sätt att göra detta på.

• Använda SIMD-instruktionerna i processorn.

• Utnyttja flera kärnor i processorn om den har stöd för det. • Göra beräkningarna med ett modernt grafikkort [4].

I det här examensarbetet har endast SIMD-instruktioner används för att pa-rallellisera beräkningarna.

Profiler

För att kunna avgöra var man ska lägga ner mest energi på att optimera sin kod är det en god idé att använda en profiler. En profiler samlar statistik när ett program körs och visar hur lång tid som gått åt i de funktioner som anropats. Verktyget

gprof har används i detta fall och finns gratis till de flesta platformar.

Bildformat

Enbart gråskalebilder har använts. Förmodligen skulle färgbilder kunna ge bättre resultat då mer information finna tillgängligt, men eftersom hög hastighet priorite-rats över god noggrannhet har det valts bort. I blockmatchningsmetoden används 8-bitars heltalsbilder och i den fasbaserade metoden 32-bitars flyttalsbilder. Det senare beror på att OpenCV kräver samma indata-typ som utdata-typ i sina funk-tioner. Annars är 8-bitars bilder att föredra eftersom ursprungsbilderna lagras i det formatet och går snabbare att räkna med.

SIMD

SIMD betyder Single Instruction Multiple Data och är en uppsättning instruk-tioner hos processorn som kan göra en och samma beräkning på flera tal i taget. Intel introducerade tekniken i sina x86-processorer år 1997 med namnet MMX. Åtta register adderades till processorn (MM0-MM7) som var och ett består av 64 bitar. Det innebär att man för 8-bitars tal kan göra 8 beräkningar åt gången, för 16-bitars tal 4 beräkningar i taget osv.

(39)

Den stora begränsningen hos MMX är att instruktionerna bara stödjer heltal. SSE, SSE2, SSE3, SSSE3 och SSE4 är uppföljare till MMX och stöder alla flyttal.

Det finns i huvudsak tre sätt att använda SIMD-intruktioner.

• Kompilatorn autogenererar kod som använder SIMD-instruktionerna. Intels egna C/C++ kompilator stödjer detta liksom de senaste versionerna av gcc. • Instruktionerna anropas manuellt med assemblerkod. Ett tidsödande sätt

som ger användaren mest kontroll.

• Intruktionerna anropas manuellt med C-lik kod. En variant som gör det möjligt att skriva vanlig C kod som i själva verket är assembler-kod. Nedan visas ett exempel på hur samma beräkningar anropas på de tre olika sätten. • Autogenerad kod: unsigned char a, b; b += a; • Assembler paddb %mm0 %mm1 • C-lik kod __m64 a, b; b = _mm_add_pi8(a, b);

Det sistnämnda användes i detta projekt efter som det upplevdes som lättast att använda.

3.2 Blockmatchning

Beräkningen av SAD är den tyngsta delen i blockmatchningen. Tre italienare har utformat en metod som gör detta effektivt [12]. Genom att minimera de redun-danta beräkningarna kan man minska beräkningsbördan markant. För att öka hastigheten ytterligare har SIMD instruktioner använts vilket utnyttjar den pa-rallellism som finns i metoden.

3.2.1 Försteg

I det inledande försteget beräknas medelvärdet och variansen för både vänster och höger bild. Medelvärdet för ett N2= (2n + 1) × (2n + 1) stort block i vänsterbilden beräknas som µ(u, v) = 1 N2 n X i=−n n X j=−n L(u + j, v + i) = 1 N2S1(u, v) (3.1)

(40)

och variansen kan beräknas som σ(u, v) = 1 N2 n X i=−n n X j=−n L2(u + j, v + i) − µ2(u, v) = 1 N2S2(u, v) − µ 2_{(u, v) (3.2)}

Eftersom medelvärdet och variansen beräknas för alla block i bilden i en sekven-tiell följd, kan man optimera (3.1) och (3.2) genom att utnyttja tidigare resultat. S1(u, v + 1) kan då skrivas som

S1(u, v + 1) = S1(u, v) + US1(u, v + 1) (3.3)

där US1(u, v + 1) = n X j=−n (L(u + j, v + n + 1) − L(u + j, v − n)) (3.4)

Vi kan optimera beräkningen ytterligare genom att skriva US1(u, v + 1) som

US1(u, v + 1) = US1(u − 1, v + 1)+

+ (L(u + n, v + n + 1) − L(u + n, v − n))− − (L(u − n − 1, v + n + 1) − L(x − n − 1, v − n))

(3.5)

En sista optimering kan införas och resulterar i att vi kan skriva US1(u, v + 1) som

US1(u, v +1) = US1(u−1, v +1)+(L(u+n, v +n+1)−L(u+n, v −n))−T1(˜u) (3.6)

där

T1(˜u) = (L(u − n − 1, v + n + 1) − L(u − n − 1, v − n)) (3.7)

och ˜u = u mod (2n + 1)

På samma sätt kan vi beräkna S2(u, v + 1) som

S2(u, v + 1) = S2(u, v) + US2(u, v + 1) (3.8)

där US2(u, v + 1) = US2(u − 1, v + 1) + (L 2_{(u + n, v + n + 1) − L}2_{(u + n, v − n)) − T} 2(˜u) (3.9) och T2(˜u) = (L2(u − n − 1, v + n + 1) − L2(u − n − 1, v − n)) (3.10) med ˜u = u mod (2n + 1)

3.2.2 Disparitetsberäkning

Precis som i försteget kan vi reducera de redundanta beräkningarna genom att ut-nyttja tidigare beräkningar. Vi har ett (2n+1)×(2n+1) stort block i vänsterbilden

(41)

L centrerat i punkten (u, v) och ett lika stort block i högerbilden R centrerat i punkten (u − d, v). SAD-värdet mellan blocken beräknas med (2.16) och ges av

SAD(u, v, d) = n X j=−n n X i=−n

|L(u + j, v + i) − R(u − d + j, v + i)| (3.11)

Man kan nu se att SAD(u + 1, v, d) kan fås från (3.11) som

SAD(u + 1, v, d) = SAD(u, v, d) + U (u, v + 1, d) (3.12) där U (u, v + 1, d) = n X j=−n |L(u + j, v + n + 1) − R(u − d + j, v + n + 1)| − − n X j=−n |L(u + j, v − n) − R(u − d − j, v − n)| (3.13)

U representerar skillnaden mellan SAD-värdena för raden längst upp och längst ner på det matchande blocket. Även detta U kan beräknas lite smartare som

Den här effektiva tekniken gör att beräkningsbördan inte ökar med blockenstorle-ken.

(42)

3.2.3 SIMD-optimering

Genom att utnyttja SIMD-instruktionerna beräknas SAD för åtta dispariteter i taget. Detta ger i sin tur begränsningen att metoden bara stödjer maxdispariteter med en multipel av 8. Exempel 3.1 visar hur absolutebeloppen mellan åtta pixlar i vänster och höger bild beräknas.

Exempel 3.1 // Skapar MMX-variablerna __m64 leftImage; __m64 rightImage; __m64 diff1; __m64 diff2; __m64 sad;

// Hämtar 8 pixlar från vänsterbilden leftImage = _mm_set1_pi8(left[row][col]); // Hämtar 8 pixlar från högerbilden

rightImage = _mm_set1_pi8(right[row][col-disp]); // Beräknar leftImage-rightImage, negativt svar ger 0. diff1 = _mm_subs_pu8(rightImage, leftImage);

// Beräknar rightImage-leftImage, negativt svar ger 0. diff2 = _mm_subs_pu8(leftImage, rightImage);

// Lägger ihop de båda ovanstående resultat för att // få ett korrekt absolutbelopp

sad = _mm_or_si64(diff1, diff2);

Ju färre bitar varje tal representeras av desto fler beräkningar kan göras pa-rallellt. Genom att tillåta en maximal blockstorlek på 16 × 16 räcker det att spara SAD för varje block som ett 16-bitars heltal. Det gör att vi kan jämföra fyra SAD åt gången när vi söker efter det lägsta SAD. Vidare så jämförs inte alla SAD med varandra utan som i figur 3.1. Detta innebär att vi garanterat hittar minimum och tre andra värden som är små.

3.3 Fasbaserad stereo

Implementationen av den fasbaserade stereometoden har implementerats i C++ med hjälp av OpenCV. Tyvärr har inga SIMD-operationer använts då det tids-mässigt inte fanns utrymme för det.

(43)

6

1

4

3

8

5

2

9

5

1

4

2

1

4

1 pseudominima

minimum

SAD-värden

Figur 3.1. Exempel på hur åtta SAD-värden jämförs. De tre pseudominima syns i den streckade rutan.

(44)

(45)

Kapitel 4

Utvärdering

Kvaliteten och beräkningshastighet hos de två stereometoderna har testats på någ-ra olika typer av bilder. Panåg-rametnåg-rarna i de olika metoderna har varienåg-rats för att se hur de påverkat slutresultatet. Ingen uttömmande sökning har gjorts som prö-var alla möjliga kombinationer av parametrar. Istället har parametrarna prö-varierats kvalitativt kring de värden som ger bäst resultat enligt litteraturen.

I blockmatchningsmetoden där vänster-höger-test används har följande para-metrar varierats:

1. Blockstorlek. Från 3 × 3 upp till 15 × 15. 2. Maxdisparitet. Från 8 till 24.

3. Beräkningar med och utan SIMD-instruktioner.

I blockmatchningsmetoden där dispariteten bara skattats åt ett håll har föl-jande parametrar varierats:

1. Blockstorlek. Från 3 × 3 upp till 15 × 15. 2. Maxdisparitet. Från 8 till 24.

3. Beräkningar med och utan SIMD-instruktioner.

4. Tröskelvärdet för texturparametern som avgör huruvida ett område har för lite textur och ska förkastas i disparitetsberäkningen.

5. Tröskelvärdet för skarphetstestet. 6. Tröskelvärdet för tydlighetstestet. I den fabaserade metoden varierades:

1. Bandbredd och centerfrekvens på kvadraturfiltret. 2. Djupet på skalpyramiden.

(46)

34 Utvärdering

4.1 Kvalitet och hastighet

Kvaliteten har mätts med det felmått som används på Middlebury-webbsidan [9] och hastigheten har mätts med kommandot time när metoderna körts 100 gånger i följd.

Felmåttet B anger andelen felaktigt skattade pixlar och beräknas som B = 1

N X

(u,v)

(|dC(u, v) − dT(u, v)| > δd) (4.1)

där δd är feltoleransen för dispariteten och är satt till 1. Det innebär att dispariter

som avviker mer än 1 pixel ifrån den sanna bilden klassas som felaktig. B har beräknats för följande områden i bilden:

• Hela bilden förutom kanten runt bilden där vi förmodar att disparitetsskatt-ningen är väldigt osäker pga kanteffekter i stereometoderna. Felet betecknas här med BA. Se figur 4.1(a).

• Områden som syns i både vänster och höger bild och därmed inte är skymda i någon av bilderna. Felet betecknas här med BS. Se figur 4.1(b).

• Områden där variationen av dispariteten är större än en viss tröskel. Felet betecknas här med BD. Se figur 4.1(c).

.

Bilderna som använts i detta test är två rektifierade foton från Middlebury-sidan [9], tsukuba och venus med storlekarna 384 × 288 och 434 × 383, 450 där den sanna dispariteten är känd. Se figur 4.2 och 4.3.

4.2 Resultat

Resultaten från utvärderingen redovisas nedan i tabell 4.1 och med en bild per metod och testbild. Parametrarna har valts genom prövning för att få ett så litet fel som möjligt. Detta garanterar inte bästa resultat men ger en hint om vad metoderna ger för resultat.

bildpar metod BA BS BD fps tsukuba blockmatchning 16.8 18.8 43.8 11.2 tsukuba fasbaserad 23.8 25.5 56.1 1.45 venus blockmatchning 11.5 12.5 41.0 6.14 venus fasbaserad 23.2 24.8 58.8 0.79 Tabell 4.1. Testresultat

(47)

4.2 Resultat 35

(a) Hela bilden. (b) Synliga områden.

(c) Områden med stor variation av dispariteten.

Figur 4.1. De olika områdena i bilden tsukuba där de två felmåtten beräknas. Endast de vita pixlarna tas med i beräkningarna av felen.

(48)

36 Utvärdering

(49)

4.3 Metodegenskaper 37

(a) Blockmatchning. (b) Fasbaserad.

Figur 4.4. Resultat med bildparet Tsukuba.

4.3 Metodegenskaper

I förra avsnittet redovisades endast ett resultat för varje metod och testbild. Men eftersom båda metoderna har en uppsättning parametrar man kan skruva på kan resultaten förstås variera. Därför ska vi i detta avsnitt undersöka parametrarnas påverkan på slutresultatet.

Parametrar hos blockmatchningsmetoden Här finns tre parametrar att skruva på:

• Blockstorlek • Maxdisparitet

• Beräkningar med eller utan SIMD-instruktioner.

Blockstorlek Ett stort block ger en säkrare matchning än ett mindre då fler pixlar är med i SAD-beräkningen. Dock har stor blockstorlek nackdelen att dispa-riteten smetas ut. Se figur 4.6. Och i och med att metoden är implementerad på det sett den är gjord så inverkar inte blockstorleken på beräkningsbördan.

(50)

38 Utvärdering

(a) Blockmatchning. (b) Fasbaserad.

Figur 4.5. Resultat med bildparet Venus.

(a) Litet block. (b) Stort block. Figur 4.6. Disparitetsskattningar med olika blockstorlekar.

(51)

Maxdisparitet Till skillnad från blockstorleken så inverkar maxdispariteten i hög grad på beräkningsbördan. I figur 4.7 visas hur beräkningstiden ökar med ökande maxdisparitet. 0 10 20 30 40 50 60 70 80 2 4 6 8 10 12 14 maxdisp fps

Figur 4.7. Maxdisparitet plottad mot fps (frames per second).

SIMD-operationer Att använda SIMD-operationer gav en hastighetsökning med en faktor två.

4.3.1 Parametrar hos den fasbaserade metoden.

Här finns det i huvudsak två parametrar att skruva på. • Bandbredd och centerfrekvens på kvadraturfiltret. • Djupet på skalpyramiden.

Kvadraturfiltret I figur 4.8 visas disparitetsskattningen av bilden Tsukuba för två olika kvaddraturfilter. Som vi ser ger resultatet i figur 4.8(b) mer brus än resultatet i figur 4.8(d). Det beror på att det första filtret släpper igenom mer höga frekvenser än det andra.

Djupet på skalpyramiden Djupet på skalpyramiden avgör vilken maxdispa-ritet som kan skattas. För att illustrera det används bildparet linear som visas i figur 4.9. Resultatet med djup 1, dvs ingen skalpyramid alls ses i figur 4.10(b). Som vi ser så förmår inte metoden skatta dispariteter över 5 pixlar. Ökas djupet

(52)

40 Utvärdering !3 !2 !1 0 1 2 3 0 0.2 0.4 0.6 0.8 1

(a) Absolutbeloppet av kvadratur-filtret i frekvensdomänen. (b) Resultatbild övre. !3 !2 !1 0 1 2 3 0 0.2 0.4 0.6 0.8 1 (c) Absolutbeloppet av kvadratur-filtret i frekvensdomänen. (d) Resultatbild undre.

(53)

Figur 4.9. Linear - syntetiskt skapat bildpar med känd disparitet.

kan större dispariteter skattas vilket ses i figur 4.10(c) och 4.10(d). Slutsatsen vi kan dra är att det är nödvändigt med en skalpyramid för att metoden ska kunna hantera stora dispariteter.

(54)

42 Utvärdering

(a) Korrekt disparitet. (b) Skaldjup 1.

(c) Skaldjup 2. (d) Skaldjup 3.

Figur 4.10. Exempel på hur skaldjupet inverkar på disparitetsskattningen med den fasbaserade metoden.

(55)

Kapitel 5

Diskussion

5.1 Slutsats

I inledningen till denna rapport beskrevs syftet med detta arbete. Målet var att implementera och utvärdera två olika stereoalgoritmer i realtid och avgöra vilken som är bäst. Som vi sett i resultatbilderna skiljer sig metoderna åt en del i både kvalitet och prestanda. Överlag presterar blockmatchningsmetoden bättre än den fasbaserade. Dock ska man inte lita för mycket på dessa resultat då mest tid och energi har lagts ned på den förstnämnda. En bättre implementation av den fasbaserade metoden skulle krävas för att göra en mer rättvis jämförelse, men i detta projekt presterar blockmatchningsmetiden överlägset bäst. Riktmärket som presenterades som var fem bilder per sekund på en 800x600 bild nåddes inte riktigt, men ligger inte långt borta.

5.2 Svårigheter och problem

Teorin

En stor skillnad mellan metoderna är deras komplexitet. Teorin för blockmatch-ningsmetoden är lätt att förstå vilket leder till att man snabbt kan få igång en implementation som fungerar tillfredställande. Den fasbaserade metoden kräver däremot djupare matematisk kunskap och kan ta längre tid att behärska. Meto-den är också bökigare att implementera då delresultaten i de olika stegen kan vara relativt svårverifierbara. I litteraturen är blockmatchningsmetoden mer förekom-mande än den fasbaserade vilket talar för att den är mer utbredd. Det i sin tur gör att man lättare kan hitta tips och idéer till hur man kan utveckla sin metod.

C++

Implementeringen skedde i två steg. Det första var att få fungerande kod i Matlab och det andra att konvertera den till C++. Anledningen till det var att i första

(56)

44 Diskussion

skedet slippa tänka på programmeringstekniska detaljer som C++ medför och i stället fokusera på matematiken. Detta arbetssätt visade sig vara lyckat.

5.3 Framtiden

Ämnet stereoseende är väldigt stort och denna rapport har bara skrapat lite på ytan. De bästa stereometoderna som finns idag presterar mycket bättre än meto-derna som implementerats i detta projekt, både när det gäller noggrannhet och prestanda. Med hjälp av grafikkort kan man idag göra 128 beräkningar parallellt [6]. Det kan jämföras med åtta som gjordes i vår implementation. Enligt [6] ger en grafikkortsimplementation en prestandaökning på mellan 5 och 20 gånger i vanliga fall jämfört med en implementation med endast en CPU. I vissa fås en ökning med upp emot 100 gånger.

Hur ser framtiden ut? En trend som redan startat inom datorvärlden är att vi går mot system som består av fler och fler beräkningsenheter. Det gör att parallella stereometoder blir allt viktigare. I med att utvecklingen går så fort på detta område ska det bli spännande att se vart vi kommer vara om fem år.

(57)

Litteraturförteckning

[1] Jean-Yves Bouguet. Camera calibration toolbox for matlab. http://www. vision.caltech.edu/bouguetj/calib_doc, 2007.

[2] Intel Corporation. Opencv. http://sourceforge.net/projects/ opencvlibrary, 2007.

[3] Gösta H. Granlund and Hans Knutsson. Signal Processing for Computer

Vision. Kluwer Academic Publishers, 1995. ISBN 0-7923-9530-1.

[4] Johan Hedborg. Gpgpu: Bildbehandling på grafikkort. Master’s thesis, Lin-köpings Universitet, Linköping, Sweden, 2006. LITH-MAI-EX–06/12-SE. [5] Janne Heikkila and Olli Silven. A four-step camera calibration procedure

with implicit image correction. In CVPR ’97: Proceedings of the 1997

Confe-rence on Computer Vision and Pattern Recognition (CVPR ’97), page 1106,

Washington, DC, USA, 1997. IEEE Computer Society.

[6] David Luebke and Greg Humphreys. How gpus work. Computer, 40(2):96– 100, 2007.

[7] Karsten Mühlmann, Dennis Maier, Jürgen Hesser, and Reinhard Männer. Calculating dense disparity maps from color stereo images, an efficient imple-mentation. In IEEE Conference on Computer Vision and Pattern Recognition, Kauai Marriott, Hawaii, December 2001.

[8] Daniel Scharstein and Richard Szeliski. A taxonomy and evaluation of dense two-frame stereo correspondence algorithms. Int. J. Comput. Vision, 47(1-3):7–42, 2002.

[9] Daniel Scharstein and Richard Szeliski. Middlebury stereo vision page. http: //cat.middlebury.edu/stereo, 2007.

[10] Maria Magnusson Seger. Lite om kamerageometri. 5 oktober 2005.

[11] Johan Skoglund. Robust real-time estimation of region displacements in video sequences. Lic. thesis, ISY, Linköping, Sweden, 2007. 978-91-85715-86-2. [12] L. Di Stefano, M. Marchionni, and S. Mattoccia. A pc-based real-time stereo

vision system. MG&V, 13(3):197–220, 2004. 45

(58)

46 Litteraturförteckning

[13] Carl-Johan Westelius. Focus of Attention and Gaze Control for Robot

Vi-sion. PhD thesis, Linköpings Universitet, SE-581 83 Linköping, Sweden, 1995.

Dissertation No 379, ISBN 91-7871-530-X.

[14] Zhengyou Zhang. Flexible camera calibration by viewing a plane from unk-nown orientations. In ICCV, pages 666–673, 1999.

(59)