Bildbaserad estimering av rörelse för
reducering av rörelseartefakter
Mats Jogbäck
2006-03-31
Bildbaserad estimering av rörelse för reducering av
rörelseartefakter
Examensarbete utfört i Bildbehandling
vid Tekniska högskolan i Linköping
av
Mats Jogbäck
LiTH-IMT/MI20-EX–06/429–SE
Handledare: Björn Svensson
imt, Linköpings universitet
Examinator: Hans Knutsson
imt, Linköpings universitet
Abstract
Before reconstructing a three dimensional volume from an MR brain imaging se-quence there is a need for aligning each slice, due to unavoidable movement of the patient during the scanning. This procedure is known as image registration and the method used primarily today is based on a selected slice being the reference slice and then registrating the neighbouring slices, which are assumed to be of minimal deviation.
The purpose of this thesis is to use another method commonly used in computer vision - to estimate the motion from a regular videosequence, by tracking markers indicating movement. The aim is to create a robust estimation of the movement of the head, which in turn can be used to create a more accurate alignment and volume.
Sammanfattning
För att kunna rekonstruera en tredimensionell volym av en hjärna avbildad med hjälp av magnetresonanstomografi (MRI) behöver man korrigera varje snittbild i förhållande till varandra, beroende på oundvikliga rörelser hos den röntgade pati-enten. Detta förfarande kallas bildregistrering och idag använder man sig primärt av en metod där en bild utses till referensbild och därefter anpassas närliggande bilder, som antas ha en minimal avvikelse, till referensen.
Syftet med detta examensarbete är att använda en annan metod vanligen ut-nyttjad inom datorseende för att estimera ett rörelsefält utifrån en vanlig vide-osekvens, genom att följa markörer som indikerar rörelse. Målet är att skapa en robust estimering av huvudets rörelse, som då kan användas för att skapa en mer noggrann korrigering och därmed också en bättre rekonstruktion.
Tack
Jag vill tacka följande personer; Hans Knutsson, huvudansvarig för Medicinsk Informatik (MI) vid Institutionen för Medicinsk Teknik (IMT), för att jag fick möjligheten att göra mitt examensarbete vid institutionen, och Björn Svensson, doktorand vid IMT, för all handledning och hjälp jag har fått under tiden mitt arbete har pågått. Jag skulle även vilja tacka Maria Engström, Johan Kilhberg och Kerstin Gustavsson på CMIV för möjligheten att få en demonstration av en magnetkamera. Därtill vill jag tacka Marcel Warntjes på CMIV för information om övervakningskameran som används i samband med magnetkameraundersök-ningarna.
Innehåll
1 Inledning 1 1.1 Bakgrund . . . 1 1.2 Tidigare arbete . . . 2 1.3 Problembeskrivning . . . 2 1.4 Arbetets upplägg . . . 2 1.5 Rapportens upplägg . . . 4 2 Teori 5 2.1 MRI . . . 5 2.1.1 Rekonstruktion . . . 62.2 Teknisk data om videokameran . . . 7
2.3 Bildregistrering . . . 7 2.3.1 Medicinsk bildregistrering . . . 8 2.4 Datorseende . . . 9 2.5 Projektiv geometri . . . 9 2.5.1 Koordinatsystem . . . 12 2.6 Epipolär geometri . . . 12 2.6.1 Punktmatchning . . . 14 2.6.2 Fundamentala matrisen . . . 14 2.6.3 Extrahera förflyttningsinformation . . . 14 3 Metod 17 3.1 Litteraturstudie . . . 17 3.2 Simulering . . . 17 3.2.1 Projektiv geometri . . . 17 3.2.2 Epipolär geometri . . . 18 3.2.3 Punktmatchning . . . 19 3.2.4 Lineär estimering av F . . . 21 3.2.5 Robust estimering av F . . . 22 3.2.6 Extrahera förflyttningsinformation . . . 24 4 Resultat 27 4.1 Mätningsresultat . . . 27 4.1.1 Gradientbaserad registrering . . . 27 4.1.2 Detektering av outliers . . . 28 xi
4.1.3 Estimering av F . . . 29 5 Diskussion 41 5.1 Slutsatser . . . 41 5.2 Framtida arbete . . . 42 5.2.1 Punktmatchning . . . 42 5.2.2 Epipolär geometri . . . 42 5.2.3 Estimering av F . . . 42 5.2.4 Modellering av rörelse . . . 43 5.2.5 Praktiska aspekter . . . 43 Litteraturförteckning 45
Kapitel 1
Inledning
I detta kapitel beskrivs bakgrund, inklusive tidigare arbete utfört inom området reducering av rörelseartefakter i MR, samt problembeskrivning och en planering för examensarbetet. Målsättning och en översikt av rapporten ges också här.
1.1
Bakgrund
Inom medicinsk teknik är kravet på noggrannhet mycket hög, särskilt hanteringen av data i olika former. Man eftersträvar att undvika alla typer av förenklingar, interpoleringar och avrundningar för att bibehålla validiteten hos den insamlade informationen. Därför är det viktigt att ständigt uppfinna och utvärdera nya sätt att behandla data på ett korrekt sätt.
När det gäller MRI1 finns det ett stort antal artefakter man vill undvika
-flertalet beroende på tekniska inställningar och omständigheter, men även patien-torsakade problem kan resultera i oanvändbara data.
En oundviklig företeelse är att patienten rör på sig under undersökningen vil-ket orsakar en inbördes förskjutning mellan snittbilderna. Även en ytterst liten rörelse kan ge en avvikelse från föregående snittbild. Detta blir till ett problem när man sedan ska rekonstruera en tredimensionell volym av snittbilderna - eftersom objekt ligger på olika spatiella positioner under skanningen blir rekonstruktionen inkorrekt om den inte tar till hänsyn detta faktum. För att undvika detta krävs ett extra steg mellan skanningen och rekonstruktionen där man korrigerar snitt-bilderna så att deras spatiella positioner överensstämmer.
Den metod man i huvudsak använder sig av idag bygger på att en snittbild väljs till referensbild, därpå korrigeras närliggande snittbilder därefter under antagandet att avvikelsen bilderna emellan är minimal. Detta förfarande kallas
bildregistre-ring2 och beskrivs mer i 2.3.
Artifakter i MR-bilder ofrivilligt orsakade av patienten är av naturliga skäl oundvikliga, då rörelser såsom andning, hjärtrörelser, blodflöde, ögonrörelser,
svälj-1
Magnetic Resonance Imaging 2
Image registration
2 Inledning
ningar och andra rörelser har en inverkan på de resulterande bilderna. Symptom på rörelseartifakter kan vara att de behandlade objekten blir suddiga eller att man får så kallade släpeffekter3. På grund av detta använder man sig av olika metoder
för att reducera dessa artifakter, metoder som är skilda åt beroende på vilken slags rörelse man försöker kompensera för. Det kan röra sig om enklare metoder, till exempel att inhämta data då rörelsen i fråga är minimal, eller att använda sig av kortare inhämtningsintervall, till mer komplicerade som k-space-segmentation eller medelvärdesbildande typer.
1.2
Tidigare arbete
När det gäller huvudrörelser behövs ett delmoment innan man rekonstruerar det skannade objektet, som tar till hänsyn den rigida rörelse som objektet har ge-nomgått under skanningen. Metoden beskriven i [16] estimerar parametrarna för en okänd rigid rörelse genom att minimera en energifunktion som baseras på att största energin i en ideal MR-bild ligger inom det skannade objektets gränser. En-ligt Aghaeizadeh et. al. [16] är tekniken tillräckEn-ligt robust för att klara av bilder med brus, men är begränsad till 2D-rörelser. Andra kompensationsmetoder som på likartat sätt postumt används i liknande situationer beskrivs i 2.3.
1.3
Problembeskrivning
Problemet med de kompenserande metoderna ovan nämnda är just att de är kom-penserande; man försöker inte ta reda på den underliggande fysiska rörelse som är orsaken till de individuella snittbildernas spatiella skillnad, utan istället försöker man att kompensera alla på referensbilden följande bilder. Detta är baserat på ett felaktigt antagande att bilderna är identiska medan de i verkligheten skiljer sig åt. Metoden bygger således på en tvådimensionell grund och ger ett tvådimensionellt korrigeringsresultat, medan den faktiska förflyttningen är tredimensionell, se figur 1.2.
Idén om att använda metoder avsedda för datorseende grundar sig på att för-söka estimera den underliggande fysiska rörelsen på ett korrekt sätt och sedan använda denna information för att rekonstruera snittbilderna därefter. Metoder vanligen använda inom datorseende begagnar sig ofta av just denna typ av estime-ring för att kunna rekonstruera godtyckliga scener. Tanken är att man ska kunna estimera rörelsen från videosekvensen från en kamera som används i samband med magnetkameraundersökningen.
1.4
Arbetets upplägg
Grundtanken med detta arbete är att skapa en simulerad miljö där man genom att följa markörer kan bestämma ett rörelsefält av en kropp som har utsatts för rigida rörelser. Utifrån detta rörelsefält skall man sedan kunna bestämma parametrar
3
1.4 Arbetets upplägg 3
Figur 1.1. Snittbilderna av ett mänskligt huvud ses till vänster, den högra bilden visar
den bild som kommer registreras enligt förflyttningsfältet markerat med pilar.
Figur 1.2. Snittbilderna har här utsatts för roterande rörelser och en enkel registrering
i bildplanet kommer bara vara tillräcklig för att få en godtagbar bild om rörelsen är liten. Om man vet hur själva rörelsen ser ut kan man istället göra en bättre registrering.
4 Inledning
för transformationerna och utnyttja dessa för en korrigering av de projicerade tvådimensionella bilderna.
Avgränsningar är att själva korrigeringen ej kommer genomföras. Transforma-tionerna är begränsade till rigida, det vill säga translation och rotation. Ingen deformation är därmed tillåten, punkterna kan antas vara distribuerade över en kropp som rör sig enhetligt.
1.5
Rapportens upplägg
I kapitel 2 ges teoretisk information som används som grund för nästkommande kapitel. I kapitel 3 redovisas hur arbetets gång har fortlöpt och tillvägagångssätt för detta. Rapporten avslutas med ett kapitel (4) där arbetets resultat visas samt ett därpå följande kapitel (5) innehållandes en diskussion om olika aspekter av arbetet.
Kapitel 2
Teori
I detta kapitel ges den teoretiska bakgrund som är vital för att förstå kommande kapitel. Hur magnetresonanstomografi (MRI) fungerar, den tekniska informatio-nen om videokameran, hur koordinatsystemet är uppställt i arbetet, på vilket vis projektiv geometri är relaterat till arbetet och den epipolära geometrins koppling till detta, är ämnen som tas upp.
2.1
MRI
Magnetresonanstomografi, förkortat MRI eller MRT, är en diagnostisk teknik som använder sig av magnetfält för att inhämta bildinformation från en patient. Den baseras på det fundamentala fenomenet kärnmagnetisk resonans1 vilket innebär
att vissa atomkärnor innehar ett svagt magnetiskt moment orsakat av elementar-partiklarnas spinn, och i påverkan av ett starkt magnetfält kan flera atomkärnors moment ställa in sig efter detta, antingen i riktning med eller emot. I detta till-stånd kan de påverkade atomerna exciteras, det vill säga absorbera energi om de utsätts för korta radiovågspulser (elektromagnetisk energi), så kallade RF-pulser2,
med samma resonansfrekvens, energi som därefter avges när elektronerna i atom-kärnorna återgår till sitt jämviktsläge. Denna emitterade energi detekteras och konverteras till bilder.
Patienten placeras i en så kallad magnetkamera som genererar ett mycket starkt magnetfält, ofta tiotusentals gånger starkare än jordens magnetfält. Den typ av atomkärnor man använder sig av är väteatomer, eftersom de finns naturligt i stora mängder i kroppen (ca. 63%), i form av vatten och i fettvävnader. Intensiteten hos ett detekterat organ beror därför på hur mycket väte som finns däri. MR är följdaktligen bäst lämpat för att undersöka mjukdelar, såsom hjärna, ryggmärg, hjärta, lever, mjälte. Eftersom tekniken är icke-invasiv och inte använder sig av någon form av röntgenstrålning kan en patient utsättas för upprepade undersök-ningar utan att ta fysisk skada.
1
Nuclear magnetic resonance 2
Radio Frequency
6 Teori
2.1.1
Rekonstruktion
För att kunna avbilda och skilja olika skannade områden åt krävs att de utsätts respektive för unika magnetfält - ursprungligen utsätts de för ett generellt mag-netfält vilket orsakar en allmän riktningsjustering för de nukleära magmag-netfälten, magnetiseringsvektorerna. Det man använder för att åstadkomma detta är gradi-entmagnetfält som mäter variationen för magnetfältet i varje riktning, Gx, Gy och
Gz. De olika gradientfälten appliceras vid olika tidpunkter - Gz samtidigt som
RF-pulsen för att bestämma vilken vinkel bildplanet ska ha emot z-axeln. Bild-planets koordinater [x, y] kodas efter att RF-pulsen stängts av, i x-led betecknas förfarandet frekvenskodning då en unik resonansfrekvens för varje magnetiserings-vektor bestäms (se ekvation 2.1-2.2). Det generella magnetfältet som befinner sig i magnetvektorns isocenter [x, y, z] = [0, 0, 0] betecknas B0, resonansfrekvensen där
n0och det nukleära magnetfältets gyromagnetiska kvot g.
n = g(Bo+ xGx) = no+ gxGx (2.1)
x = n − no
gGx (2.2)
I y-led kallas tillvägagångssättet faskodning och gradientfältet för detta, Gy,
appliceras efter att frekvenskodningsgradienten stängts av, så att alla magnetise-ringsvektorer har en identisk resonansfrekvens. Utifrån Gy kan fasvinkeln Φ emot
en referensaxel för respektive magnetiseringsvektor bestämmas.
Den råkodade bilden har då information i form av frekvens-fas vilket i MR-sammanhang brukar betecknas k-space. Bilden läses in i datorn där man sedan genomför rekonstruktionen. En vanlig metod är att utföra en invers
Fouriertrans-form, som transformerar data från frekvensdomänen, Fourierdomänen, till
spati-aldomänen.
Från den multidimensionella Fouriertransformen, ekvation 2.3, kan förklaras hur relationen med MR-signalen ser ut.
X(u) = Z . . . ∞ Z −∞ f (x)e−iuTx dx (2.3)
där frekvenser i varje led betecknas u= [ux uy uz]T och position x= [x y z].
Jämför det med en MR-signal S(k) bildad enligt ekvation 2.4. S(k) = ∞ Z −∞ ∞ Z −∞ I(x, y)e−in0
e−i(gxGx+gyGy)dxdy (2.4)
där I(x, y) är bilden och gxGx och gyGy är frekvenserna i respektive led. Sätt
f (x) = I(x, y)e−in0och jämför detta tillsammans med den andra termen e−i(gxGx+gyGy)
med ekvation 2.3 och den därifrån utvecklade termen e−iuTx
= e−i(uxx+uyy). Om
uxx = (gGx)x och för y på liknande vis gäller kan S(k) tolkas som en
2.2 Teknisk data om videokameran 7
Att rekonstruera med hjälp av Fouriertransform kan dock innebära olika pro-blem, till exempel så är transformerna definierade över en (eller flera) typer av grid vilket kan orsaka problem om man försöker använda standardverktyget fast fourier
transform(fft) för att beräkna dessa. Ett annat exempel som orsakar svårigheter
är om patienten rör på sig vilket resulterar i ett inhomogent magnetfält som or-sakar störning hos frekvenserna. Matematiskt innebär det att en extra term som varken beror på x eller y ingår i magnetfältet vilket betyder att S(k) inte längre är en Fouriertransform. Ett försök att rekonstruera med en invers Fouriertransform kommer i det här läget att ge ett felaktigt resultat.
Då kunskapsområdet MRI är för stort för att till fullo behandlas här rekom-menderas [5] där ämnet behandlas på signalbehandlingsnivå.
2.2
Teknisk data om videokameran
Den videokamera som kan användas i samband med en MR-undersökning är först och främst tänkt för patientövervakning i den mån detta behövs. Tanken bakom det här arbetet är att i en framtida version använda videosekvensen från kame-ran och beräkna förflyttningsfält utifrån detta. Kamekame-ran är av webbkameratyp, dvs. med relativt låg upplösning (384x287 (CCIR)) och med en svartvit CMOS-bildsensor. Både exponering och kontrast-ökning/minskning3 sköts automatiskt.
Under minimikravet på 0.5 lux i belysning kan ett SNR4 på 52 dB uppnås på
videosignalen. Videochippet sitter i ett kamerahus som är placerat 16 cm ovan-för patienten vid användning och kameran har en 4.3 mm lins. Vid användning monteras kameran i sin medföljande ställning på britsen så att den överblickar patienten.
2.3
Bildregistrering
Bildregistrering innebär att man geometriskt sett korrigerar två eller fler bilder av samma scen tagna vid olika tidpunkter5, från olika vyer6och/eller med olika
sen-sorer7så att de överensstämmer. Detta är ett viktigt steg i all slags bildbehandling
där man hämtar information från kombinationer av flera bilder, ett område som ständigt växer.
Metodiken för majoriteten av registreringsmetoder kan sammanfattas i följande fyra steg:
1. Hitta särdrag. Utmärkande drag i bilden såsom punkter, linjer, hörn, re-gioner osv. detekteras automatiskt eller manuellt och lagras på lämpligt vis. Utmärkande drag behöver inte vara direkt visuella, de kan till exempel utgöras av alla pixlar i en bild, beroende på registreringsmetoden.
3 Gain control 4 Signal-to-Noise Ratio 5 Multitemporal analysis 6 Multiview analysis 7 Multimodal analysis
8 Teori
2. Detaljmatchning. Detaljerna ovan detekterade i referensbilden och den be-handlade bilden matchas så att deras korresponderande relation bestäms. 3. Estimering av transformationsmodell. Parametrar för hur själva
registre-ringen skall gå till väga modelleras.
4. Estimering av modellens parametrar. Parametrar för transformationsmodel-len estimeras, man söker förflyttningsfältet mellan bilder.
5. Registrering och omsampling. Den behandlade bilden korrigeras, bildvärden beräknas genom någon form av interpoleringsteknik om så är nödvändigt. Registreringsmetoden man väljer att använda bör anpassas till det använd-ningsområde som är aktuellt, således kan implementeringsteknikerna för de olika stegen vara vitt skilda. Dock finns det ett antal krav metoden bör uppfylla - de-taljdetektionsalgoritmen får inte vara känslig för en eventuell förändring av bilder (till exempel i avseende om kontrast, ljusstyrka osv.), idealt sett ska den kunna detektera samma detaljer i alla slags scenprojektioner oberoende av deras natur. Detaljmatchningssteget måste vara tillräckligt robust för att kunna hantera fel-matchningar8 och brus. Parameterestimeringen bör ta i beaktning hur scenen är
uppbyggd och de villkor detta ger. En mycket viktig del är också att besluta om vilka skillnader som skall behandlas i korrigeringen, ett beslut grundat på vad för slags information man vill ta fram i slutänden. Dessutom måste man väga in hur noggrannt resultat man vill ha jämfört med beräkningskraften som tas i anspråk, om detta är en begränsning.
För en mer detaljerad beskrivning av hur de förstnämnda stegen kan genomfö-ras, exempel på olika slags tillvägagångssätt samt jämförelser däremellan rekom-menderas [15].
2.3.1
Medicinsk bildregistrering
Enligt Maintz et. al. [9] kan de medicinska bildregistreringsmetoderna delas in i två övergripande kategorier: extrinsiska, baserat på att utomstående objekt används i bilden (någon form av markörer som placeras på patientens behandlade område) och intrinsiska, baserat på att använda den befintliga bildinformation som genererats av patienten.
För en överblick av metoder för medicinsk bildregistrering rekommenderas [9]. De metoder som används huvudsakligen är av typen optiskt flöde9 och kan delas
in ytterligare i ett antal underkategorier [1]:
• Gradientbaserade metoder. Rörelseflödet för varje pixel beräknas genom att använda första eller andra ordningens spatiotemporala derivator. Skillnaden mellan bilderna antas enbart vara en förflyttning, ej till exempel intensi-tetsvariation, och en lokal omgivning i bilden ska kunna approximeras med en Taylor-utveckling. Metoderna av denna typ är inte lämpliga när brus i bilderna existerar eller om aliasing förekommer i bilderna. Eftersom flödet beräknas på pixelnivå är dessa metoder dessutom inte särskilt noggranna.
8 Outliers 9
2.4 Datorseende 9
• Regionbaserad matchning. Flödet beräknas här genom att titta på två på varandra följande bilder och hitta korresponderande objektspixlars position. Man försöker hitta bästa matchningen genom att maximera ett likhetsmått, till exempel genom att minimera ett avståndsmått som SSD10.
• Energibaserade metoder. Genom att använda frekvensfilter anpassade för olika hastigheter kan man beräkna energin utifrån bilder. Metoderna går också under namnet frekvensbaserade metoder på grund av filtren.
• Fasbaserade metoder. Även här används filter inställda för olika fasriktning-ar, till exempel kvadraturfilter. Fördelen med kvadraturfilterfasen är att den är invariant mot signalenergin. Bilder där kontrast och medelnivå varierar med tiden, eller omgivningar med lokala toppar och dalar utgör inga problem för dessa metoder.
2.4
Datorseende
Grunden till datorseende är att utifrån bildinformation få ett datorsystem att tol-ka sin omgivning och använda detta i det för systemets utformade ändamål. Hur denna bildinformation kan se ut och hur tolkningsmetoden kan te sig är vitt skilt och forskningsområdet kring datorseende är mycket stort. Användningsområden är också dessa av ett stort antal och växer ständigt - självstyrande robotar, fö-rarlösa flygplan, industriella avsyningssystem är exempel på applikationsområden där datorseende används på olika sätt.
Metodiken i ett datorseendesystem är att först bestämma sig för hur bildinfor-mationen skall tas in samt diskretiseras och lagras. Därpå ställs en kameramodell upp där användningsmiljön tas i beaktning och lämpliga koordinatsystem därefter, se 2.5 och 2.5.1 för mer information om hur sådant har ställts upp i detta arbete. När detta är gjort används någon form av algoritm för att extrahera särdrag vars uppgift är att hitta intressanta detaljer i bilderna. De resulterande punkterna ut-tagna från olika vyer jämförs sedan för att få matchande punkter, detta kan bland annat göras genom att använda geometriska villkor såsom epipolär geometri, se 2.6.
2.5
Projektiv geometri
Euklidisk geometri betraktas av gemene man som bekant då dess egenskaper in-nebär en inneboende metrik för all geometri - sidor har längd, linjer mätt mot skärande linjer har vinklar och linjer liggandes i plan som aldrig möts betraktas som parallella. Dessa egenskaper ändras heller aldrig då transformationer som translation och rotation införs. Men för att korrekt modellera en kamera be-hövs en annan typ av geometri, eftersom kamerans process gör att de euklidiska egenskaperna inte längre gäller - längder och vinklar behöver nödvändigtvis inte
10
10 Teori
längre bevaras och parallella linjer kan mycket väl mötas i bilderna. Här kommer
projektiv geometri in.
Tanken bakom projektiv geometri härstammar från konstvärlden, närmare be-stämt renässansen då man började att representera människor på olika avstånd genom att måla dessa i olika storlekar. Den enklaste formen av modell för att demonstrera detta är hålkameran11, se figur 2.1. Utgående från fokalpunkten F
kan det hela ses som en mappning från objektsrymden till den projektiva rymden, R37→ R2.
Figur 2.1. En modell av en hålkamera. I figuren ses den projicerade avbildningen
på bildplanet C där de parallella linjerna i R3
konvergerar i en så kallad idealpunkt v, förutsatt att linjerna fortsätter mot “oändligheten”.
Eftersom man inte kan beskriva ett bildplans punkter i samma koordinatsystem som det avbildade objektets koordinatsystem på grund av ovanstående skäl behö-ver man transformera objektskoordinaterna till bildplanskoordinater med hjälp av
projektiva transformationer.
Transformeringen från den euklidiska rymden Rn, objektsrymden, till den
pro-jektiva rymden Pninnebär att man här kan representera punkter vid oändligheten,
Pn är därigenom en förlängning av Rn. Ett enkelt sätt att visa detta är genom användandet av homogena koordinater. Införandet av homogena koordinater in-nebär att varje punkt ges en tredje koordinat, [ˆx ˆy w] ǫ R3, en virtuell dimension.
Om w = 0 tyder det på att punkten befinner sig vid oändligheten, i Pn en
ideal-punkt, vilket är en distinkt egenskap i Pn- punkter som befinner sig i oändligheten
mappas till godtyckliga punkter. Då w 6= 0 kan punktens homogena koordinater divideras med w för att få de kartesiska koordinaterna för x och y, [ˆx
w ˆ y
w 1], vilket
kan ses som en mappning från 3D till planet w=1, R3 7→ R2 : w = 1. Vilket
följer att en punkt med givna kartesiska koordinater [x y] kan skrivas i homogena koordinater som [x y 1]. Detta gäller analogt för Rn.
Homogena koordinater kan även ses som ett matematiskt trick för att få ett enhetligt sätt att genomföra transformationer lineärt. Rent beräkningsmässigt
11
2.5 Projektiv geometri 11
utförs skalnings- och rotationstransformationer genom matrismultiplikationer me-dan translationer görs genom addition. Då man inför homogena koordinater kan samtliga operationer istället utföras som matrismultiplikationer. En translation P = [x y]T + (T
x Ty) kan exempelvis genomföras enligt ekvation 2.5.
Kombina-tioner av flera transformaKombina-tioner kan göras genom att multiplicera de respektive transformationsmatriserna och använda den resulterande matrisen i en operation.
P = T p = 1 0 Tx 0 1 Ty 0 0 1 x y 1 (2.5)
För att applicera projektiv geometri då man skapar en bild från ett objekt brukar man modellera objektsvärlden som en tredimensionell projektiv rymd P3
och bildplansmodellen som en tvådimensionell sådan, P2. Således blir processen
en mappning P37→ P2, vilket brukar benämnas en kameramodell. Bektrakta en
punkt i P3 skriven i homogena koordinater [X Y Z W ]T. Det går att visa [4]
att den sista koordinaten är irrelevant och att motsvarande bildpunkt är punkten [X Y Z]T i P2också skriven i homogena koordinater. Mappningen kan således ses
som en mappning av homogena 3D-koordinater, representerade i form av en 3 × 4-matris P med blockstrukturen P = [I3×3| O3] där I3×3är en enhetsmatris och O3
en nollvektor, förutsatt att projektionscenter ligger i origo [0 0 0 1]T. Den mest
generella bildprojektionen, vilken bland annat tillåter ett annat projektionscenter, kan representeras i form av en 3×4-matris av rang 3. Matrisen P mappar punkter P37→ P2enligt ekvation 2.6 och brukar betecknas kameramatrisen.
x y w = P X Y Z W (2.6)
P kan skrivas om på formen P = AR[I | − ˜C] (se finite projective camera [4]), där R är en 3 × 3-rotationsmatris som beskriver hur objektskoordinater ro-teras till kamerakoordinater, ˜C är projektionscenter i objektskoordinater och A,
kamerakaliberingsmatrisen, är definierad enligt 2.7.
A = αx s x0 αy y0 1 (2.7)
där αx och αy är skalningsfaktorer för pixelkoordinater i förhållande till
bildko-ordinater, det vill säga fokallängd i pixlar längs respektive axel. Detta gäller i synnerhet kameror som har icke-kvadratiska pixlar, såsom CCD-kameror. Para-metern s är en så kallad skjuvningsparameter vilken för de flesta kameror är noll. De sista två parametrarna (x0, y0) är kameracentrum beskriven efter
pixeldimen-sioner. Parametrarna för A beskriver en kameras intrinsiska sådana, medan R och ˜C beskriver de extrinsiska. Totalt sett har en kameramodell av ovan nämnd typ 11 frihetsgrader, 11 okända parametrar. För en hålkamera reduceras antalet frihetsgrader till 9, eftersom att parametrar för linsen inte gäller.
12 Teori
Utgående från två vyer med tillhörande kameramatriser P och P′ och med
antagandet att objektskoordinatsystemet är associerat med den första vyn ges matriserna enligt följande:
P = [A | 0] P′
= [AR | At] (2.8)
där R och t återigen är en rotationsmatris och en translationsvektor som beskriver förflyttningen vyerna emellan.
2.5.1
Koordinatsystem
Koordinatsystem i denna rapport är definierade enligt följande; ett objektskoor-dinatsystem med sfäriska koordinater, där punkt P = [r θ ϕ]T transformeras till
kartesiska koordinater enligt ekvation 2.9. X Y Z = r sin θ cos ϕ r sin θ sin ϕ r cos θ (2.9)
Objektspunkterna [X Y Z]T projiceras därefter ner till ett bildplan och
bild-planskoordinaterna [x y]T enligt ekvation 2.10, se figur 2.2.
x y = −Zf X Y (2.10) där fokallängden f definieras enligt ekvation 2.11. nx betecknar bildbredden i
pixlar och αxsynfältsvinkeln12 i samma led.
f = nx 2 tanαx
2
(2.11) Rotationsmatriserna som används i rapporten är definierade enligt ekvation 2.12, translationsvektorerna trivialt. Rx= 1 0 0 0 cos α sin α 0 − sin α cos α Ry = cos α 0 − sin α 0 1 0 sin α 0 cos α Rz= cos α sin α 0 − sin α cos α 0 0 0 1 (2.12)
2.6
Epipolär geometri
Den grundläggande tanken för epipolär geometri i två matchande vyer är att för varje punkt i den ena vyn finns en definierad relation till den korresponderande punkten i den andra vyn. Genom att utnyttja detta faktum kan man teore-tiskt sett beräkna rymdkoordinaterna för de matchande punkterna och därigenom rekonstruera scenen. Epipolär geometri bygger på projektiv geometri, eftersom
12
2.6 Epipolär geometri 13 P p x y Z X Y f
Figur 2.2. Bilden illustrerar den projektiva rymden där objektspunkten P mappas till
bildpunkten p.
informationskällan är bilder som bygger på att en tredimensionell struktur har projicerats ner till bildplan.
P C1 C2 p 1 p2 c 1 c2 e1 e2 l 1 l2
Figur 2.3. En geometrisk tolkning av den epipolära geometrin. Den epipolära
geo-metrin beskriver relationen mellan de båda projicerade bildpunkterna p1 och p2 från
objektspunkten P .
Betrakta en given punkt P som är synlig i båda vyerna definierade av bild-planen C1 och C2, enligt figur 2.3. Punkten projiceras på bildplanen, till p1 och
p2. Epipolärplanet definieras som planet beskrivet av P och bildplanens
optis-ka centrum, c1 och c2. Linjen mellan c1, c2, baslinjen för epipolärplanet, brukar
benämnas just den epipolära baslinjen. Dess skärningspunkter med bildplanen C1, C2 kallas epipoler, i figur 2.3 betecknade e1 och e2. Epipolärlinjerna l1 och
l2definieras av epipolärplanets intersektion med bildplanen, linjerna passerar
ge-nom respektive punkter (p1, e1) och (p2, e2). I enlighet med denna geometriska
14 Teori
motsvarande punkt i den andra bilden ligga längs epipolärlinjen. Bakgrunden till epipolär geometri kommer från Kruppa’s ekvationer, [7] rekommenderas för både en geometrisk tolkning och en algebraisk härledning av dessa.
2.6.1
Punktmatchning
En av hörnstenarna i en implementering av en fortsättning baserat på detta arbete kommer att vara punktmatchningen - är den bristfällig kommer eventuella fel fortplanta sig i efterföljande steg. Man bör därför sträva efter att använda sig utav en så robust algoritm som möjligt.
2.6.2
Fundamentala matrisen
Den fundamentala matrisen F är en algebraisk representation av epipolär geometri, den intrinsiska projektiva geometri mellan två vyer finns kodad häri. Från sektion 2.6 fastställdes att för varje punkt p1i den ena vyn fanns en motsvarande epipolär
linje l2 i den andra vyn, vilken längs den korresponderande punkten p2 låg. Den
fundamentala matrisen beskriver denna avbildning.
En geometrisk tolkning av F är att den representerar en mappning från det projektiva planet P2 av första bilden (bildplanet C
1 i figur 2.3) till knippet av
epipolära linjer genom epipolen e2, en mappning enligt P27→ P. För en algebraisk
derivering av F utgående från kameramatriserna P och P′
för respektive vyer refereras till [4].
Ett lineärt villkor, vilket brukar betecknas epipolarvillkoret, vilket F lyder under är att för alla par av korresponderande punkter p1 ↔ p2 i två vyer gäller
2.13, vilket följer från det faktum att den korresponderande epipolära linjen för p1
är l2= F p1och på motsvarande sätt l1= F p2 för p2.
pT2F p1= 0 (2.13)
F är en 3 × 3-matris av rang 2, med 7 frihetsgrader. En vanlig 3 × 3-matris har 9 element, men skalningen för F är inte signifikant vilket gör att den har 8 frihetsgrader. Detta innebär att ett bildpar med en inbördes förflyttning på ett visst avstånd från kameran kommer att få ett likadant F som ett bildpar på ett annat avstånd, så länge som förhållandet mellan de respektive förflyttningarna och kameraavstånden är lika stort för båda bildparen. F lyder dessutom under bivillkoret att den är singulär, rang(F)=2 → det F = 0, vilket reducerar frihetsgra-derna ytterligare ett steg. Är F definierat för ett par av kameramatriser (P, P′)
är FT den motsvarande fundamentala matrisen för den omvända ordningen på
kameraparet (P′
, P ).
2.6.3
Extrahera förflyttningsinformation
Tanken är att utgående från den fundamentala matrisen F utnyttja de parametrar den beskriver för att få ut den förflyttningsinformation som skiljer vyerna åt. Enligt [4] kan den essentiella matrisen E, vilket är den fundamentala matrisens motsvarande variant för kalibrerade kameror, skrivas utifrån en kombination av
2.6 Epipolär geometri 15
en känd rotationsmatris R och translationsvektor t = [tx ty tz]T enligt ekvation
2.14:
E = t × R = R[T ]X (2.14)
där [T ]X definieras enligt ekv. 2.15:
[T ]X = 0 −tz ty tz 0 −tx −ty tx 0 (2.15)
Från E kan även R och t extraheras [4].
Relationen mellan den essentiella matrisen E och den fundamentala matrisen F är enligt ekvation 2.16:
E = ATF A (2.16)
Med detta i åtanke bör man kunna beräkna rotation och translation två vyer emellan om en fundamental matris F och en kalibreringsmatris A (se sektion 2.5) är kända.
Kapitel 3
Metod
Detta kapitel beskriver hur det faktiska examensarbetet gick till väga, med en förstudie av litteratur och därpå en simuleringsfas där själva implementationen av de valda teknikerna anges.
3.1
Litteraturstudie
Arbetet inleddes med en litteraturstudie av relevanta artiklar om bildregistrering för att få en övergripande bild av hur detta fungerar och vilka metoder som ex-isterar [9] [15], en genomgång ges i 2.3. Därefter studerades projektiv geometri [10] [11], epipolär geometri och dess användning inom datorseende [2] [3] [4] [6] [8] [14]. Teorin bakom dessa områden ges i 2.5, 2.6 samt 2.4.
3.2
Simulering
För att kunna testa de metoder som bestämts i fråga om arbetets upplägg använ-des Matlab som simuleringsmiljö. Matlab är en ypperlig miljö att arbeta i vad gäller bildbehandling då programmeringsspråket är matris- och vektororienterat. Därtill hör att det innehåller ett stort antal fördefinierade optimerade matematiska funktioner, många speciellt anpassade för linjär algebra. Förutom detta finns ett antal verktygslådor1som fungerar som specialiserade tilläggsfunktioner baserade
på programmeringsspråket i Matlab. Programmet har även ett stort stöd för visu-alisering av data och kombinationen av den interaktiva miljön och programmering medger ett lättanvändligt gränssnitt.
3.2.1
Projektiv geometri
Projektiv geometri, som beskrivet i 2.5, beskriver relationen mellan ett objekt definierat i Rn och dess projektiva motsvarighet i Rn−1, det vill säga Rn7→ Rn−1.
1
Toolboxes
18 Metod
Eftersom arbetet är tänkt för vår spatiella rymd är simuleringerna gjorda över projektionen R37→ R2.
För att kontrollera att beräkningarna från objektskoordinater till bildplansko-ordinater stämmer enligt ekvation 2.10 ställs ett system upp där en punkt P = [r θ ϕ]T placeras i R3och därefter utförs en förflyttning enligt P′
= t + RP , där t betecknar en translation och R en rotation. Koordinaterna för de projicerade bild-punkterna p och p′
beräknas enligt ekvation 2.10. Fokallängden f bestäms genom ekvation 2.11, där nysätts till en godtycklig skalning vilken betecknar bildstorleken
och αy sätts till synfältet för bilden. Kameran sätts i position [0 0 0]T.
Systemet testas genom att förflyttningen utförs ett led i taget och bara av en typ i taget (till exempel en translation i x-led) och därefter beräknas kvoten mel-lan den förflyttade objektspunkten P′ och dess motsvarande beräknade bildpunkt
p′. Denna kvot jämförs sedan med den beräknade fokallängden f. Resultatet
redovisas i tabell 3.1. Test P’ p’ Kvot f 1. T = [0.01 0 0] [0.2203, 0.0961] [4.6251, 2.0189] -41.9977 40 2. T = [0 0.01 0] [0.1396, 0.2006] [2.9116, 4.1823] -41.7057 40 3. T = [0 0 0.01] [0.0754, -0.2154] [1.5979, -4.5641] -42.3729 40 4. Rx= π/180 [-0.2489, 0.0245] [-5.0317, 0.4947] -40.4337 40 5. Ry = π/180 [0.1123, -0.1808] [2.3533, -3.7868] -41.8967 40 6. Rz= π/180 [0.1096, -0.2213] [2.2355, -4.5138] -40.7938 40
Tabell 3.1. Test för projektiv geometri
Från resultatet kan ses att skillnaden mellan kvoten och den beräknade fo-kallängden f först och främst är tecknet på respektive, detta beroende på hur den optiska axeln, i detta fall z-axeln, är definierad - pekandes inåt eller utåt från bilden. Av tabellen ses, om man bortser från tecknet, att det finns en viss skillnad mellan de två intressanta värdena. Detta kommer ifrån det faktum att f beräknas på subpixelnivå medan kvoten beräknas på pixelnivå, vilket innebär att kvoten mäts rent manuellt efter pixelkoordinater medan f beräknas rent matematiskt.
3.2.2
Epipolär geometri
Den epipolära geometrin, beskriven i sektion 2.6, och i synnerhet den fundamenta-la matrisen, beskriven i sektion 2.6.2, beräknas här utifrån två dataset av punkter. Ett fördefinierat n antal punkter P distribueras över en sfärisk yta såsom tidigare enligt en slumpmässig normalfördelning, och utsätts sedan för en transformation till P′
= RP + t där R är en rotation och t en translation. Respektive punkter projiceras ner till ett bildplan enligt ekvation 2.10 och bildpunkterna p och p′
. För att få en mer realistisk bild av hur de matchande punkterna i slutänden kommer se ut utsätts de transformerade punkterna p′
för ett gaussiskt brus i vissa tester. Således kommer varje matchande punktpar inte entydigt tyda på en enskild trans-formation utan en liten felvarians uppstår vilket är mer realistiskt tänkbart än idealfallet.
3.2 Simulering 19
Den fundamentala matrisen F beräknas för hela punktseten p, p′ genom att
använda den normaliserade 8-punktsalgoritmen beskriven i sektion 3.2.4.
3.2.3
Punktmatchning
Det huvudsakliga fokuset i det här arbetet har legat på att få en fungerande ro-bust estimering av den fundamentala matrisen och de metoder som krävs för att kunna extrahera rörelse därigenom. Prioriteten för punktmatchning har inte varit lika hög på grund av att de metoder som förmodligen kommer användas (se ka-pitel 5 för en diskussion kring detta) redan finns utvecklade och implementerade vid institutionen. Men för att få en förståelse för hur bildregistreringsmetoder av typen optiskt flöde fungerar implementerades och testades en gradientbase-rad metod, vilken redovisas här. Fördelen med en ggradientbase-radientbasegradientbase-rad metod är att den går att använda på godtyckliga bilder, man slipper till exempel att designa specialanpassade filter.
Metoden för förflyttningsfältet v(x) : R2 7→ R2, x = [x, y] beskrivet nedan
bygger på två antaganden:
1. Rörelsen två bilder I1ochI2emellan kan lokalt beskrivas som en förflyttning
∆x. Detta innebär att det enda som skiljer bilderna åt är en förflyttning, exempelvis förändras inte intensitet i punkter som motsvarar varandra.
I1(x) = I(x, t) = I(x − ∆x, t + 1) ≈ I2(x + v(x)) (3.1)
2. Bilden kan lokalt beskrivas som ett lutande plan. Vilket betyder att en lokal omgivning i bilden kan approximeras med en Taylor-utveckling av första ordningen.
I(x, t) = I(x, t) − ∇ITv + ∆
t (3.2)
där ∆t= I(x, t + 1) − I(x, t).
Under detta söks sedan förflyttningsfältet v(x) som minimerar skillnaden mel-lan bild I1 och I2, det vill säga det v(x) som minimerar ǫ givet av ekvation 3.3.
ǫ2= ||I2(x + v(x)) − I1(x)||2 (3.3)
Från antagande 2 kan ekvation 3.2 förenklas till ekvation 3.4, vilken beskriver det optiska flödet v.
∇ITv − ∆t= 0 (3.4)
Eftersom det kommer finnas många lösningar för v i varje bildpunkt behöver man en modell för hur förflyttningsfältet bör se ut och därefter söka den lösning som stämmer så bra som möjligt för alla bildpunkter.
Förflyttningsfältet v(x) modelleras uttryckt i parametervektorn p, modellen kan då skrivas som v(x) = B(x)p där B(x) är basmatrisen för p. Olika modeller förekommer, men ett affint förflyttningsfält (translation, rotation, skalning) mo-delleras med en parametervektor p i sex dimensioner. Förflyttningsfältet ser då ut som följer:
20 Metod v(x) = p1 p2 + p3 p4 p5 p6 x y = 1 0 x y 0 0 0 1 0 0 x y | {z } B(x) p (3.5)
Lösningen fås genom att hitta parametervektorn p som minimerar ǫ i följande uttryck:
ǫ2=X
i
(∇I(xi)TB(xi)p − ∆t(xi))2 (3.6)
Derivering av ǫ med avseende på p ger ekvation 3.7. ∂ǫ2 ∂p = 2 X i BT i ∇Ii(∇IiTBip − ∆ti) = 0 (3.7) En omskrivning ger X i BTi ∇Ii∇IiTBi | {z } A p =X i BiT∇Ii∆ti | {z } h (3.8)
vilken gör att lösningen fås som popt= A−1h.
En mer kompakt lösning kan också användas. I termer av I1, I2, B och p kan
ekvationen för optiskt flöde (3.4) skrivas om enligt ekvation 3.9. ∇IT 2B(x) | {z } ˜ A p ≈ I1(x) − I2(x) | {z } b (3.9)
Matrisen ˜A och vektorn b bestäms därefter enligt följande:
˜ A = ∇I2(xT1)B(x1) ∇I2(xT2)B(x2) .. . ∇I2(xTn)B(xn) b = I1(x1) − I2(x1) I1(x2) − I2(x2) .. . I1(xn) − I2(xn)
Detta ger ekvation 3.10 vilken känns igen från ovan. popt= ( ˜ATA˜ | {z } =A )−1A˜Tb |{z} =h (3.10) Tester för den implementerade gradientbaserade registreringsmetoden redovi-sas och kommenteras i kapitel 4.
3.2 Simulering 21
3.2.4
Lineär estimering av F
En lineär estimering av F utgående från 2.13 kan beräknas genom 8-punktsalgoritmen, från början utvecklad av Longuet-Higgins [6], vilken först användes till att beräkna den essentiella matrisen, den fundamentala matrisens motsvarighet för kalibrerade kameror. Samma algoritm gäller även för F och okalibrerade kameror. Som nam-net tyder på krävs minst 8 korresponderande punktpar för att kunna bestämma F för två vyer.
Om elementen i F numreras enligt 3.11 och de matchande punktseten p1 =
[u v 1]T och p
2= [u′v′ 1]T kan F skrivas om i form av en 9 element lång vektor f
och därefter kan ekvation 2.13 skrivas om på följande vis för n korresponderande punkter, ekvation 3.12. F = f11 f12 f13 f21 f22 f23 f31 f32 f33 (3.11) u′ 1u1 u′1v1 u′1 v ′ 1u1 v′1v1 v′1 u1 v1 1 u′ 2u2 u′2v2 u′2 v ′ 2u2 v′2v2 v′2 u2 v2 1 .. . ... ... ... ... ... ... ... ... u′ nun u′nvn u′n v ′ nun vn′vn vn′ un vn 1 | {z } M f11 f12 f13 f21 f22 f23 f31 f32 f33 | {z } f = 0 (3.12)
När matrisen M har bildats kan en fundamental matris ˆF estimeras genom att beräkna nollrummet ur M genom SVD2. Utifall att ˆF inte uppfyller rang( ˆF )=2,
då den har fler än 2 singulära egenvärden större än noll, kan man påtvinga singu-lariteten genom att helt enkelt sätta det minsta singulära egenvärdet till noll.
8-punktsalgoritmen har ofta kritiserats för dess känslighet för brus, ett flertal mer komplicerade algoritmer har därmed föreslagits. Hartley visade dock [3] att 8-punktsalgoritmen i sin enkelhet kan uttökas för att få ett mer robust resultat - genom en normalisering av punkterna som ett förebyggande steg innan esti-meringen av F sker kan bruskänsligheten reduceras betydligt. Problemet, enligt Hartley, bygger på att en typisk punkt p = [x y w]T = [100 100 1]T - p har (x,
y)-koordinater som är mycket större än w, vilket kommer ha en inverkan på estimatet ˆ
F . Hartley’s förslag är därför att genomföra en isotropisk skalning vilken innebär att koordinaterna har medelpunkten i origo och att medelavståndet till origo är √
2. Ett antal andra metoder finns också, varav åtta finns beskrivna i [4] där dessa också jämförs.
2
22 Metod
3.2.5
Robust estimering av F
För att få en mer robust estimering av F används den metod som beskrivs i [14], baserad på en Least Median of Squares-teknik (LMedS). Motiveringen till att just den här metoden valdes är dels för att artikeln i fråga är en välciterad sådan samt att Zhang är ett relativt välkänt namn i forskningen kring epipolära geometri och den fundamentala matrisen. Metoden innebär att man kan få en estimering som trots så kallade outliers fortfarande ger ett tillräckligt användbart resultat, något som den lineära metoden inte klarar av. Enligt Zhang et. al. [14] finns det outliers orsakade av två olika skäl:
• Felaktiga spatialpositioner. Punkter som inte följer det modellerade normal-fördelade bruset utan har en större avvikelse.
• Felmatchningar. Punkter som matchas till helt fel punkter. En möjlig an-ledning till dessa är att punkterna endast finns i en av vyerna.
En jämförelse mellan LMedS-metoden och närbesläktade så kallade M-estimatorer visar att de senare inte är robusta mot outliers av typen felmatchningar. En annan populär teknik är RANSAC vilken, likt LMedS, kan användas för att robust esti-mera en parameter. RANSAC kräver dock en initieringsparameter för felvariansen medan LMedS klarar sig utan denna. Nackdelen med LMedS är att estimeringen misslyckas då antalet outliers överstiger 50%, men för det aktuella arbetet moti-veras användandet av LMedS med att fördelarna överväger nackdelarna.
Tanken bakom LMedS-metoden är att få en robust grund för att kunna mi-nimera ett meningsfullt fysikaliskt mått. I [14] föreslås att mimi-nimera avståndet mellan punkter och deras motsvarande epipolära linjer, vilket i sin tur ger en ro-bust estimering av F . Definitionen för det euklidiska avståndet mellan en punkt p2 och dess epipolära linje F p1 ges av ekvation 3.13.
d(p2, F p1) = | p T 2F p1| p (F p1)21+ (F p1)22 (3.13) där (F p1)i är komponent i av vektorn F p1.
LMedS-metoden estimerar genom att lösa det icke-lineära minimeringsproble-met min med
|{z}
i
r2
i, där ri är residualen, det vill säga skillnaden mellan det
observe-rade värdet ai och dess estimerade värde ˆai. r2i definieras enligt ekvation 3.14.
ri2= d2(p2i, F p1i) + d2(p1i, FTp2i) (3.14)
Detta innebär att estimatorn måste få det minsta värde för medianen av de kvadrerade residualerna för hela datasetet. För att estimera F behövs, som nämnts tidigare, minst 8 punkter men en bild innehåller typiskt sett betydligt fler punkter än så. Eftersom punktrymden då kan anses vara för stor för att användas i sin helhet används en Monte Carlo-teknik för att slumpvist välja ut m subsampel om p = 8 punkter. Antalet m subsampel definieras enligt ekvation 3.15, för en härledning se [12].
3.2 Simulering 23
m = log(1 − P )
log(1 − (1 − ǫ)p) (3.15)
P är sannolikheten att minst ett subsampel är utan outliers, ǫ andelen outliers i hela datasetet och p antalet punkter i varje subsampel.
En estimering gjord efter att ha dragit m subsampel slumpvist kan dock ge en instabil lösning om inte ytterligare krav införs, till exempel om punkterna som genereras i ett subsampel ligger väldigt nära varandra, då de utvalda punkterna inte är en rättvis representation av hela punktrymden. Därför används den
för-delningsteknik3 som också beskrivs i [14]. Den baseras på att punktrymden delas
in i b × b hinkar definierade enligt min, max-koordinaterna. I varje hink placeras ett antal punkter tills alla är utplacerade, de hinkar som saknar punkter exklude-ras. För varje återstående hink beräknas sannolikheten baserad på antal punkter i hinken. Sannolikheterna tas i beaktning i nästa steg då 8 hinkar slumpvist väljs ut och därpå en slumpvis punkt ur respektive hink. Zhang et. al. [14] använder b = p = 8 i sin implementering vilket är lämpligt även för detta arbete.
För varje subsampel m beräknas sedan den fundamentala matrisen Fj lineärt
enligt sektion 3.2.4 och därefter medianen av kvadrerade residualer Mj för hela
datasetet enligt ekvation 3.16 och utifrån detta behålls det Mj som är minimalt.
Mj = medi=1,2...n[d2(p2i, Fjp1i) + d2(p1i, FjTp2i)] (3.16)
För att exkludera de punkter som räknas som outliers utförs en viktningspro-cess som grundar sig i estimeringen av den robusta standardavvikelsen4[14]), vars
ekvation 3.17 innehåller tal som empiriskt bestämts enligt [12]. ˆ
σ = 1.4826[1 + 5 n − p]
p
Mj (3.17)
där Mj är den minimala medianen bestämd ovan. Viktningen wi sätt därefter
enligt följande: wi = 1 om r2 i ≤ (2.5ˆσ)2 0 annars (3.18)
Den slutgiltiga fundamentala matrisen ˆF estimeras genom att lösa det vikta-de icke-lineära minstakvadratproblemet vikta-definierat enligt ekvation 3.19. Punkter detekterade som outliers viktas bort och estimatet är robust.
minX
i
wir2i (3.19)
Rent implementationsmässigt genomförs optimeringen genom att man para-metriserar det F man utgår ifrån (vilket inte har någon betydelse) enligt 3.20 och därefter optimeras dessa parametrar efter ekvation 3.19.
3
Regularly random selection method 4
24 Metod F = b a −ay − bx −d −c cy + dx dy′ − bx′ cy′ − ax′ −cyy′ − dy′ x + ayx′ + bxx′ (3.20)
3.2.6
Extrahera förflyttningsinformation
Som beskrivet i sektion 2.6.3 kan förflyttningsparametrar extraheras ur den es-sentiella matrisen E. Enligt Hartley [4] fås två lösningar för R och två lösningar för t, ekvation 3.21. Att det blir två par lösningar beror på ambiguiteten att om E = RTT
tgäller så gör även E = ˜RTTtdet, där ˜R = SR och S betecknar en 180
graders rotation kring axeln definierad av t [11]. Viktigt att notera är att enbart riktningen för translationen fås, inte magnituden.
R1 = U W VT
R2 = U WTVT
t1 = aT
t2 = −aT
(3.21) där U och V fås ur en SVD av E, a är tredje kolumnvektorn i U och W är definierad enligt 3.22: W = 0 −1 0 1 0 0 0 0 1 (3.22)
För att beräkna E krävs dock, som nämnts i sektion 2.6.3, att man förutom en fundamental matris F även har en kalibreringsmatris A definierad för samma vyer som F är definierad för. I en artikel från Lourakis et. al. [7] utgår man från den fundamentala matrisen F för att utföra en så kallad självkalibrering där kamerans intrinsiska kalibreringsparametrar bestäms. Metoden bygger på en SVD av F och förenklade varianter av Kruppa’s ekvationer och är en rent algebraisk sådan. Förutom att genom experiment verifiera metodens effektivitet visas också i [7] att de resulterande kalibreringsmatriserna är tillräckliga för att beräkna rörelse i 3D. Fördelen med att utgå från den fundamentala matrisen är att man inte behöver någon övrig information vilket är det huvudsakliga skälet till att denna metod valdes, nackdelen är att estimeringen av den fundamentala matrisen är mycket vital för att man ska få ett bra resultat.
Låt K beteckna den symmetriska matrisen AAT, vilken kommer estimeras och
ur den beräknas sedan A. Låt sedan ρi(SF,K)
φi(SF,K), i = 1...3 vara tre ration definierade
enligt ekvation 3.23: r2vT 1Kv1 uT 2Ku2 = rsv T 1Kv2 −uT 2Ku1 = s 2vT 2Kv2 uT 1Ku1 (3.23) där SF = [r s uT1 uT2 v1T v2T] i vilken skalärerna r och s är egenvärdena för matrisen
F FT, u
1, u2, v1, och v2 är kolumnvektorerna för U- respektive V -matriserna som
3.2 Simulering 25 K = α2 x+ x20 x0y0 x0 x0y0 α2y+ y02 y0 x0 y0 1 (3.24)
där αx, och αy, vilket nämnts tidigare, är fokallängder och (x0, y0) är
principal-punkten.
För att hitta en initial lösning till K utgås från ett ekvationsystem definierat enligt ekvation 3.25:
ρ1(SF, K)φ2(SF, K) − φ1(SF, K)ρ2(SF, K) = 0
ρ1(SF, K)φ3(SF, K) − φ1(SF, K)ρ3(SF, K) = 0 (3.25)
Genom att approximera principalpunkten till bildens mittpunkt och sätta skjuv-ningsparametern till 0 (se sektion 2.5) reduceras antalet okända i ekvationssyste-met till två - de båda fokallängderna. Som mest kan systeekvationssyste-met ha fyra lösningar av vilka de som gör att matrisen K inte är reell och positivt definit förkastas.
När en initial lösning ˜K är funnen följer en icke-lineär optimering. Låt πij(SF, K)
beteckna residualen definierad av ρi(SF,K)
φi(SF,K)−
ρj(SF,K)
φj(SF,K)och låt σ
2
πij(SF, K) vara dess
varians. Variansen approximeras enligt ekvation 3.26: σπ2ij(SF, K) = ∂πij(SF, K) ∂SF ΛSF ∂πij(SF, K)T ∂SF (3.26) där ΛSF är kovariansmatrisen tillhörandes SF. För en vägledning till hur denna
räknas ut mer i detalj hänvisas till [7].
Tanken är sedan att använda varianserna σ2
πij(SF, K) för att vikta residualerna
πij(SF, K) baserat på osäkerheten hos dem. Det slutliga estimatet av K beräknas
sedan som lösningen av ett icke-lineärt minstakvadratproblem, ekvation 3.27:
K = arg min ˜ K N X i=1 π2 12(SF, ˜K) σ2 π12(SF, ˜K) + π 2 13(SF, ˜K) σ2 π13(SF, ˜K) + π 2 23(SF, ˜K) σ2 π23(SF, ˜K) (3.27) Efter att K har estimerats extraheras A därur genom tre steg; först beräknas A−T genom en Cholesky-faktorisering av K−T, därefter transponeras den och till
Kapitel 4
Resultat
Resultatet av implementationen av arbetet är en simulering utförd i Matlab. En robust fundamental matris F beräknas från ett antal punkter och deras trans-formerade motsvarigheter. Dessutom har outlierdetekteringen av metoden testats separat. Förutom den robusta estimeringsalgoritmen har även en gradientbaserad registreringsmetod implementerats och testats.
I detta kapitel redovisas tester genomförda för de olika implementerade algo-ritmerna och de resulterande mätningsvärdena samt kommentarer till dessa.
4.1
Mätningsresultat
En viktig del i att utvärdera implementerade metoder är att genomföra tester och mätningar vars resultat kan påvisa trender som man kan dra slutsatser utifrån. Testerna bör vara utformade så att resultaten har en hög validitet och reliabilitet, dvs. att den information som är väsentlig och intressant verkligen är den som mäts, samt att upprepade försök visar på samma resultat. Det finns många olika metoder för hur man skall gå till väga när man ställer upp hur testerna skall gå till, det som har använts i det här arbetet har varit av principen “test-retest” vilket innebär att man utför upprepade försök med samma instrument och förutsättningar och därefter noterar mätvärden.
4.1.1
Gradientbaserad registrering
I testerna för den gradientbaserade registreringsmetoden (se sektion 3.2.3 för en teoretisk genomgång) användes en ytmodell av ett ansikte draperat med en tex-tur designad för att underlätta registreringsprocessen. För att visa algoritmens styrka gjordes först ett test där den transformerade bilden utsattes för en uniform skalning och därpå en rotation i bildplanet (kring z-axeln), se figur 4.1-4.3. Ef-tersom de båda bilderna innehåller samma information och förflyttningsfältet är parametriserat efter en affin modell har algoritmen inga problem med att hitta en optimal lösning och registrera bilden därefter.
28 Resultat
Däremot kommer en transformation som ses i figur 4.4-4.6 orsaka problem. Här är transformationen en rotation kring y-axeln vilket gör att den vänstra de-len av ansiktet skyms, även den högra dede-len av ansiktet kommer skilja sig från referensbilden. Eftersom bildinformationen inte längre är densamma i de båda bil-derna misslyckas algoritmen med att hitta en optimal lösning och en registrering baserad på detta kommer bli felaktig. Detta är den stora nackdelen med registre-ringsmetoder av typen optiskt flöde - eftersom de bygger på att bildinformation mellan de aktuella bilderna inte skiljer sig åt utan förutsätter att skillnaden är en transformation kommer de misslyckas i de fall detta inte gäller.
4.1.2
Detektering av outliers
För att få en bild över hur väl den robusta algoritmen för att estimera F presterar genomfördes inbördes oberoende test med 10% resp. 20% outliers införda. Algo-ritmen beräknar och markerar detekterade outliers enligt figur 4.7. Observera att termen “outliers” här används för att beteckna de outliers av typ 2 som beskrivs i sektion 3.2.5, det vill säga punkter som matchas till helt fel punkter.
Resultatet av 100 utförda tester visas i figur 4.8-4.9.
Av figurerna att döma synes att algoritmen klarar sig relativt bra i fallet med 10% outliers, med bara fåtalet missade outliers, medan den i fallet med 20% klarar sig lite sämre - ett inte helt oväntat resultat. I de flesta fall då fler än riktiga outliers detekteras missas sällan de riktiga. De gånger den markerar mindre än antalet riktiga outliers är få och då missandes med enstaka punkter. Att döma av helheten ser algoritmen ut att kunna klara sig så bra som Zhang et. al. [14] vill påstå, även om detta test inte bör betraktas som uttömmande. Hela poängen med algoritmen är att få bort så många outliers som möjligt, vilken den i stor mån klarar. Att sedan “falska” outliers markeras är mindre viktigt, resultatet blir att dessa punkter inte inkluderas i den fortsatta estimeringen av F men de punkter som återstår bör vara tillräckliga för att få en robust estimering. Som nämnt i 3.2.5 fallerar metoden om över hälften av punkterna är outliers, därför bör inte heller resultat då antalet detekterade outliers överstiger denna nivå användas. Man ser i testerna ovan att detta sker stundtals och dessa resultat bör alltså förkastas.
Det bör påpekas att outliers detekteras genom att avgöra om avståndet för de korresponderande punkterna till deras motsvarande epipolärlinjer faller utan-för en viss tröskel (se sektion 3.2.5 och ekvation 3.18 i synnerhet) och eftersom de outliers som införs slumpgenereras fram inom gränserna för den ursprungliga punktrymden kan det hända att en outlier ibland sammanfaller med den ersatta punktens faktiska position, vilket gör att den inte kommer detekteras. Omvänt kan det applicerade gaussiska bruset orsaka att “riktiga” punkter detekteras som outliers, beroende på hur hög standardavvikelse bruset har. Dessutom bygger al-goritmen på att man i förväg ställer in en toleransnivå ǫ för hur många outliers punktrymden förväntas innehålla (se ekvation 3.15) och därmed hur många sub-sampel m som används - sätter man nivån så att den ligger över den förväntade klarar sig algoritmen utan större problem.
4.1 Mätningsresultat 29
4.1.3
Estimering av F
I tester av estimering av F jämförs den lineära estimeringsmetoden med den ro-busta estimeringsmetoden. Felmåttet som används är medelvärdet av avstånden för punkter till deras respektive epipolära linje för hela punktsetet och matriserna i fråga, se ekvation 3.13. Innan detta sker skalas dock matriserna om efter en sann fundamental matris ˜F eftersom jämförelsen annars inte är meningsfull. Den sanna matrisen kan beräknas, om rörelsen är känd, utifrån ekvation 2.16 där ka-libreringsmatrisen A sätts enligt ekvation 2.7. Skalningen sker därefter genom att minimera ǫ med avseende på skalfaktorn s, definierat enligt 4.1:
ǫ = i=2,j=2 X i=1,j=1 (s2F ij− ˜Fij)2+ k=2,l=3 X k=1,l=3 (sFkl− ˜Fkl)2+ m=3,n=2 X m=3,n=1 (sFmn− ˜Fmn)2+(F33− ˜F33)2 (4.1) Rörelsen som används i följande tester är en translation i x-led enligt t = [0.1 0.0 0.0]T och standardavvikelsen σ
n för det gaussiska bruset är i relation till
förflyttningen.
Först i raden av tester är ett där den lineära estimeringsmetoden jämförs med och utan normaliserade punkter, figur 4.10. Normaliseringsprocessen ska innebära att estimeringen blir mindre bruskänslig vilket testresultatet intygar.
Därefter redovisas sex olika tester med olika förutsättningar, figur 4.13-4.17, av den robusta estimeringsmetoden jämfört med den lineära metoden:
• Figur 4.11 visar hur de båda metoderna klarar av 10% outliers. • Figur 4.12 visar ett liknande test som ovan men med 20% outliers.
• Figur 4.13 visar hur metoderna påverkas när ett gaussiskt brus med σn = 0.1
är adderat.
• Figur 4.14 visar ett liknande test som ovan men med σn = 0.5.
• Figur 4.15 visar en kombination av både 40% outliers och en brusnivå med σn= 0.5.
• Figur 4.16 visar test för ökande antal punkter med 40% outliers.
• Figur 4.17 visar test för ökande antal punkter och en kombination av 40% outliers och en brusnivå med σn = 0.5.
Sammanfattningsvis kan man säga att det som är mest påfallande och oroande är att den robusta metoden visar tecken på instabilitet. Det uppstår så kallade avvikelser i mer eller mindre alla fall, men de verkar öka då ett gaussiskt brus är adderat. Med hjälp av viktningsprocessen i den robusta estimeringsalgoritmen indikeras dessa avvikelser då alla punkter viktas bort, vilket innebär att avstånden är för stora för just den optimerade lösningen av F . Bortser man från avvikelserna presterar den robusta estimeringsmetoden generellt sett bättre än den lineära.
30 Resultat
Figur 4.1. Bilden visar referensbilden (vänster) och den transformerade bilden (höger)
4.1 Mätningsresultat 31
Figur 4.2. Här ses metoden efter 15 iterationer - i övre bilden ses i ordning
referens-bilden, den registrerade bilden samt skillnaden mellan de två. I den nedre bilden visas referensbilden samt den slutliga registrerade bilden efter 15 iterationer.
32 Resultat
Figur 4.3. Metoden har här itererat 30 gånger och man kan se att skillnaden mellan
4.1 Mätningsresultat 33
Figur 4.4. Referensbilden (vänster) och den transformerade bilden (höger) vilken den
34 Resultat
Figur 4.5. Metoden efter 15 iterationer - i övre bilden ses återigen referensbilden, den
registrerade bilden samt skillnaden mellan de två. I den nedre bilden visas referensbilden samt den slutliga registrerade bilden efter 15 iterationer.
4.1 Mätningsresultat 35
Figur 4.6. Här är bilderna efter 30 iterationer. Skillnaden mellan referensbilden och
den registrerade bilden har inte minskat nämnvärt om man jämför med bilderna efter 15 iterationer, vilket tyder på att metoden inte finner någon bättre lösning. Notera det ljusare partiet i den vänstra delen av ansiktet hos skillnadsbilden - eftersom det är just denna del som skiljer bilderna åt kommer den vara störst orsak till problem för metoden.
36 Resultat −5 0 5 −6 −4 −2 0 2 4 6 Image points p −6 −4 −2 0 2 4 6 −6 −4 −2 0 2 4 6
Image points pprim, 10 detected outliers
Figur 4.7. Bilderna visar idealfallet för 10 detekterade outliers. I den vänstra bilden
finns 100 bildpunkter p beräknade utifrån de utslumpade punkterna P . I den högra bilden har deras motsvarande translaterade bildpunkter pprim beräknats och detekterade outliers har ringats in med röda ringar, införda outliers är röda punkter.
0 10 20 30 40 50 60 70 80 90 0 5 10 15 20 25
Test for 10% outliers, 100 samples
Detected outliers Amount of detection 0 0.5 1 1.5 2 2.5 3 3.5 4 0 10 20 30 40 50 60
Test for 10% outliers, 100 samples
Missed outliers
Amount of misses
Figur 4.8. Bilderna visar data för 100 genomförda tester med 10% outliers. Den vänstra
bilden är ett histogram över hur många outliers som har detekterats och den högra de outliers som missats i samma tester.
10 20 30 40 50 60 70 80 90 100 0 1 2 3 4 5 6 7 8 9 10
Test for 20% outliers, 100 samples
Detected outliers Amount of detection 0 1 2 3 4 5 0 5 10 15 20 25 30 35
Test for 20% outliers, 100 samples
Missed outliers
Amount of misses
Figur 4.9. Här visas data för 100 genomförda tester med 20% outliers. Som ovan visas
till vänster ett histogram över detekterade outliers och till höger missade outliers. Som synes blir felfrekvensen högre ju fler outliers som införs.
4.1 Mätningsresultat 37 0 10 20 30 40 50 0 1 2 3 4 5 6
Test for difference in error, σn=0.5
Sample
Error
No normalization Normalization
Figur 4.10. Figuren visar skillnaden mellan att använda en normaliseringsprocess av
punkterna (hel röd linje) eller inte (streckad blå linje). Brusnivån är konstant σn= 0.5
och 50 tester har utförts. Som synes är versionen med normalisering stabilare än den utan. 0 10 20 30 40 50 0 0.02 0.04 0.06 0.08 0.1 0.12
Test for difference in error, 10% outliers
Sample
Error
Linear estimation Robust estimation
Figur 4.11. Här ses ett test med 10% outliers, den blå streckade linjen visar resultat för
den lineära estimeringen medan den röda heldragna linjen visar den robusta estimeringen. Som synes presterar den robusta estimeringsmetoden bättre än den andra.
38 Resultat 0 10 20 30 40 50 0 0.2 0.4 0.6 0.8 1 1.2 1.4
Test for difference in error, 20% outliers
Sample
Error
Linear estimation Robust estimation
Figur 4.12. Här är ett test med 20% outliers. Även här presterar den robusta
estime-ringsmetoden bättre än den lineära.
0 10 20 30 40 50 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
Test for difference in error, σn=0.1
Sample
Error
Linear estimation Robust estimation
Figur 4.13. Figuren visar ett test med en konstant brusnivå σn = 0.1. Här kan man
se att den robusta metoden får problem vid vissa tillfällen. Generellt sett är de robusta lösningarna bättre men de olika avvikelserna är tecken på att något inte står helt rätt till. Det som även är noterbart här är att felen generellt sett är lägre i jämförelse med testerna med outliers - detta innebär att en outlier av typen felmatchning kommer ha större inverkan på estimatet än en outlier som uppstår på grund av brus.
4.1 Mätningsresultat 39 0 10 20 30 40 50 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Test for difference in error, σn=0.5
Sample
Error
Linear estimation Robust estimation
Figur 4.14. Här har brusnivån höjts till σn = 0.5 och också här uppstår två större
avvikelser och ett antal mindre för de robusta lösningarna.
0 10 20 30 40 50 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Sample Error
Test for difference in error, 40% outliers, σn=0.5 Linear estimation Robust estimation
Figur 4.15. Detta test är en kombination av 40% outliers och en brusnivå σn = 0.5.
Återigen kan man se avvikelser för den robusta metoden medan den i övrigt presterar bättre än den lineära metoden.
