Visão estéreo - correspondência e reconstrução -

Visão estéreo- correspondência e reconstrução -

Reconstrução da forma

Captura de movimento

Basic principle to recover position from stereo images: Triangulation

• Requires correspondence and camera calibration

Mapa de profundiade

Disparidade x Profundidade

Profundidade versus disparidade

zx ol x or

z

P

cl cr

xlxr

T

Z

f

fZ

xxT

Z

T rl

d

Tf

xx

TfZ

rl

Uso do mapa de profundidade

Cesar Palomo

Larry Zitnick

Time slice, bullet-time or frozen time effects

Correpondência por semelhança

Sum of Square Differences – SSD

ou

Correlação

Correspondênciapor vizinhança correlacionada

Semelhança de duas regiões WW(SSD – Sum of Squared Difference)

min),(),(),( 00 évusvuyxd

W

Wx

W

Wy

yvyxuxIyyxxIvus 2002001 ))(,)((),(),(

x0

y0

x0+u

y0+v

Semelhança de duas regiões WW(correlação)

W

Wx

W

Wy

yvyxuxIyyxxIvus 2002001 ))(,)((),(),(

W

Wx

W

Wy

yyxxIvus 2001 ),(),(

W

Wx

W

Wy

yvyxuxIyyxxI ))(,)((),(2 002001

W

Wx

W

Wy

yvyxuxI 2002 ))(,)((

constante

constante

W

Wx

W

Wyvu

yvyxuxIyyxxI ))(,)((),( 002001),(

max

Semelhança de duas regiões WW(Normalização)

W

Wx

W

Wyvu

yvyxuxIyyxxI ))(,)((),( 002001),(

max

222

2

11

2211

),()ˆ),((ˆ),(

)ˆ),()(ˆ),((max

IvyuxIIyxI

IvyuxIIyxIvu

Normalizando:

2

222

22

2

11

11

),( )ˆ),((

)ˆ),((

ˆ),(

)ˆ),((min

IvyuxI

IvyuxI

IyxI

IyxIvu

W

Wx

W

Wyvu

yvyxuxIyyxxI 2002001

),(

))(,)((),(min

Correspondence between points

• With characteristics

-

Correspondence problems: Oclusion

-

Correspondence problem: lack of characterists

-

Ostridge egg on a Chinese checker board

Correspondência com luz estruturada

Estéreo Ativo

Taxonomy of active range acquisition methods

Active shape acquisition

Contact

Nondestructive

Non-contact

Transmissive

Reflective

CT

Asla Sá et al, Coded Structure Light for 3D-Photograpy: an Overview, Revista de Informática Teórica e Aplicada, Volume IX, Número 2, Porto Alegre, 2002

Brian Curless. New Methods for Surface Reconstruction from Range Images. PhDDissertation. Stanford University. 1997

Non-optical

Sonar

Microwave radar

Optical Radar

Triangulation

Active stereo

Active depth from defocus

Passive

Shape from focus

Destructive

Slices

CMM

Shape from shading

Shape from silhouettes

Active

…

…

Active stereo solution

Use a light source to mark corresponding points

uncalibrated light source calibrated light source

One point at the time: long capture process.

Active stereo: capturing many points

Use of a digital projector as a structured light source

Pattern with several elements in a way where each element can be identified univocally

point coding: prone to errors stripes: more robust

Methods for light coding: temporal codification

Project, in sequence, a series of slides that code in the image a binary number.

n slides for 2n stripes. Two ilumination levels. Static scene. Code one axis.

Posdamer, J. L. Altschuler, M. D. Surface Measurement by Space-Encoded Projected Beam Systems. Comput. Graphics Image Process. 18, pp. 1-17, 1982.

slide1 slide2 codeslide3

problem:all transitions occur in the same place!

can be also111 or 001!

Código de Gray código binário

Código binário1 bit: 0 12 bits: 00 01 10 113 bits: 000 001 010 011 100 101 110 111 Código de Gray1 bit: 0 12 bits: 00 01 11 103 bits: 000 001 011 010 110 111 101 100

ordem invertida

Código de Gray código binário

Código binário Código de Gray

Robust temporal codification: Gray coding

Inokuchi, Seiji. Sato, Kosuki. Matsuda, Fumio. Range Imaging for 3D Object Recognition. Proc. Int. Conf. on Pattern Recognition, pp.806-808, 1984.

transitions occur in different places

Example of Gray coding

needs too many slides!

A 10-bit horizontal Gray-code sequence.

http://community.middlebury.edu/~schar/papers/structlight/p1.html

Color Gray coding

better yet…

reduces the number of slides by 3

(b,s)-BCSL Coding

Sá, Asla Medeiros. Medeiros, Esdras Soares. Carvalho, Paulo Cezar Pinto. Velho, Luiz. Coded Structured Light for 3D-Photography: an Overview. Revista de Informática Teórica e Aplicada, Vol. 9, No. 2, outubro 2002

20

A practical difficulty in the border detection

example with the monochrome Gray code

Edge detection Projecting positive and negative slides

is a robust way to recover edges.

51

60

2100

40410

01816

32rgb-BCSL coding

(+)

(-)

slide 1 slide 2

Recovering colored codes

BBBB

GGGG

RRRR

pruI

pruI

pruI

,

,

ii

ii ru

uI

if

if

1

0

i

i

p

p negative slide

positive slide

ambient lightreflection factorsprojected light

iii IIr

iii III

Implementação do BCSL//A função getBcslStripeCode retorna o código de transição de faixa conforme a seqüência de cores fornecida.//Observe a ordem em que as cores devem ser passadas:// Primeiro as cores da imagem 1 e depois da imagem 2// Primeiro a faixa da esquerda e depois a faixa da direita////O código das cores e das bases é conforme a tabela abaixo.

//Padrão 3_2 //base 3//1 - vermelho//2 - verde//3 - azul

//Padrão 4_2 //base 4//1 - vermelho//2 - verde//3 - azul//4 - magenta

//Padrão 6_2//base 6//1 - vermelho//2 - verde//3 - azul//4 - ciano//5 - magenta//6 - amarelo

int getBcslStripeCode(int base, int colorLeft1, int colorRight1,int colorLeft2, int colorRight2);

int getBcslStripeCode(int base, int colorLeft1, int colorRight1,int colorLeft2, int colorRight2){ int aux1, aux2,linha,coluna; colorLeft2--; colorRight2--; colorLeft1--; colorRight1--;

linha = (colorLeft1 * base) + colorLeft2; aux1 = (colorRight2 - colorLeft2); aux2 = (colorRight1 - colorLeft1); aux1 = (aux1>0)?(aux1-1):((base-1)+aux1); aux2 = (aux2>0)?(aux2-1):((base-1)+aux2); coluna = ((aux2) * (base-1)) + (aux1); switch(base){ case 3:

return matrix3_2[linha *4+coluna];break;

case 4:return matrix4_2[linha *9 +coluna];break;

case 6:return matrix6_2[linha *25 +coluna];break;

default:printf("Error: invalid BCSL base\n");return -1;

}}

int matrix3_2[4*9]={ 0, 3, 6, 9,14, 17, 19, 11,28, 34, 22, 24,26, 29, 18, 21, 1, 31, 33, 35,15, 4, 8, 13,16, 23, 32, 12,27, 5, 7, 25, 2, 10, 20, 30

};….

teoria pode ser complicadamas a implementação é muito simples!

Geometria Epipolar

Correspondência pela Geometria das Câmeras

Guido Gerig

Epipolar Geometry ctd.

Geometria Epipolar: notação

l

Ol

P

r

Or

Pl

plxcl

ycl

zcl

xcr

ycr

zcr

pr

Pr

el er

Linhaepipolar

Linhaepipolar

Example: converging cameras

Example: motion parallel with image plane

Example: forward motion

e

e’

Geometria Epipolar: relações básicas

l

P

r

xcl

ycl

zcl

xcr

ycr

zcr

)( rll

lrr eyePRP

ll

ll Z

fPp r

r

rr Z

fPp

rleye

rPlP

leyereye

MGattass

Produto Vetorialforma de lembrar

kjippp )()()( 21212121212121 xyyxxzzxyzzy

222

11121

zyx

zyx

kji

pp

kjipp22

11

22

11

22

1121 yx

yx

zx

zx

zy

zy

MGattass

Matriz do produto vetorial

xaya

zaxa

yaza

yx

xz

zy

pa

z

y

x

aa

aa

aa

xaya

zaxa

yaza

xy

xz

yz

yx

xz

zy

0

0

0

pa

0

0

0

xy

xz

yz

aa

aa

aa

a

Produto misto (revistado)

u

v

w

wv

h

wvbasedaárea

wv

wvu

altura

wvu alturabaseV

Matriz Essencial 0)( lrl

Trll PeyeeyeP

0 lrl

T

r

Tlr PeyePR

r

Tlr

rll

rll

lrr PReyePeyePRP )()(

llrl SPPeye

0

0

0

x

rly

rl

x

rlz

rl

y

rlz

rl

eyeeye

eyeeye

eyeeye

S

0llr

Tr PSRP 0l

Tr EPP

Matriz essencial

0lT

r Epprl

rl

ZZ

ff

l

eye l

P

r

eye r

Pl

pl

xcl

ycl

zcl

xcr

ycr

zcr

pr

Pr

el errleye

SRE lr

Rotação de a para b (left to right)

lrlrlr

lrlrlr

lrlrlrlr

kkjkik

kjjjij

kijiii

R

Twl

wr

lw

wr

lr RRRRR

lwlwlw

lwlwlw

lwlwlw

wrwrwr

wrwrwr

wrwrwrlr

kkjkik

kjjjij

kijiii

kkjkik

kjjjij

kijiii

R

rwrwrw

rwrwrw

rwrwrw

wlwlwl

wlwlwl

wlwlwlrl

kkjkik

kjjjij

kijiii

kkjkik

kjjjij

kijiii

R

Vetor do eye de b em a

lrwl

rl eyeeyeReye

X w

Y w

Z w

eye lxcl

ycl

zcl rleye

xcr

ycr

zcr

eye r

Matriz essencial (código C)

Matrix epiEssencialMatrix( Matrix Ra, Vector eye_a, Matrix Rb, Vector eye_b) { Matrix Rba = algMult(Rb,algTransp(Ra)); Vector eye = algMult(Ra,algSub(eye_b,eye_a); Matrix S = algVectorProductMatrix(eye); Matrix E = algMult(Rba,S);

return E;}

SRE lr

Matriz Essencial

T

l

Ol

P

r

Or

Pl

pl

xcl

ycl

zcl

xcr

ycr

zcr

pr

Pr

el er

0lT

r Epp

lr Epu

0rT

r up

r

r

f

y

x

p

c

b

a

ru 0 rrT

r cfbyaxup

Câmera para imagem

Z

Y

X

osf

osf

w

y

x

yy

xx

h

h

100

0

0

yy

im

xx

im

oZ

Y

s

fy

oZ

X

s

fx

wy

wx

y

x

h

h

im

im

)(

)(

yimy

ximx

oys

oxs

y

xp

ZY

ZX

f

Geometria Epipolar: Matriz Fundamental

lll PMp

0lT

r EPP

01 ll

Tr

Tr pEMMp

1 lT

r EMMF

Matriz fundamental0l

Tr pFp

lr pFu

Z

Y

X

osf

osf

w

v

u

yy

xx

100

0

0rrr PMp

lll pMP 1 Tr

Tr

Tr

MpP

lp lM lP

Matriz fundamental 1 lT

r EMMF

100

0

0

0

0

0

100

0

0

ylyl

l

xlxl

l

x

rly

rl

x

rlz

rl

y

rlz

rl

lrlrlr

lrlrlr

lrlrlr

T

yryr

r

xrxr

r

osf

osf

osf

osf

eyeeye

eyeeye

eyeeye

kkjkik

kjjjij

kijiii

F

333231

232221

131211

FFF

FFF

FFF

F Pode ser estimada diretamente se conhecermos pelo menos oito pares de pontos correspondentes

0lT

r pFp

Matriz Fundamental (código C)

Matrix epiFundamentalMatrix( Matrix Ma, Matrix Ra, Vector eye_a, Matrix Mb, Matrix Rb, Vector eye_b){ Matrix E = epiEssencialMatrix(Ra,eye_a,Rb,eye_b); Matrix invMa = algInv(Ma); Matrix invMbTransp = algTransp(algInv(Mb)); Matrix tmp = algMult(invMbTransp,E); Matrix F = algMult(tmp,invMa); return F;}

Estimativa direta da Matriz Fundamental

O algoritmo de 8 pontos

Estimating Fundamental Matrix

Each point correspondence can be expressed as a linear equation

0

1

1

333231

232221

131211

r

r

ll y

x

FFF

FFF

FFF

yx

01

33

32

31

23

22

21

13

12

11

F

F

F

F

F

F

F

F

F

yxyyyxyxyxxx rrlrlrllrlrl

The 8-point algorithm

0lT

r pFp

Estimating Fundamental Matrix

0

1

1

1

1

33

32

31

23

22

21

13

12

11

33

22

11111

F

F

F

F

F

F

F

F

F

xx

xx

xx

xyxxx

nr

nl

rl

rl

lllrl

The 8-point algorithm

0FA }]{[

TUDVA

F é a coluna de V correspondente ao menor valor singular

Estimating Fundamental MatrixThe 8-point algorithm

333231

232221

131211

FFF

FFF

FFF

F deveria ter posto 2!

TUDVF Seja D' = D com o menor valor singular = 0

TVUDF ''

FF 'min0'det F

The Normalized 8-point AlgorithmRichard Hartley

0lT

r Fpp

0'' llTr

T

r pFTTp

Tpp '

lTr FTTF '

:= T

1000 0 -500

0 1000 -500

0 0 1

:= Tt

1000 0 0

0 1000 0

-500 -500 1

1250000 250000 -500

250000 1250000 -500

-500 -500 1

The Normalized 8-point AlgorithmRichard Hartley

centro

escale para a distância média ficar em 2

Retificação de Imagens

Retificação

UNC-CH

Rectification ctd.

before

after

Guido Gerig

Retificação de imagens

z'

xc

yc

zc

y'

x'O O'

fz

y

x

z

y

x

c

c

c

'

'

'

pontoprincipal

Trucco e Verri

Retificação de imagens

l

Ol

P

r

Or

plxcl

ycl

zcl

xcr

ycr

zcr

pr

el erT

T

Te 1

0

1222 x

y

yx

T

T

TTe

213 eee

T

T

T

rect

3

2

1

e

e

e

R

Trucco e Verri

Retificação de imagens

T

T

T

rect

3

2

1

e

e

e

R

Trucco e Verri

1. Construa:

2. Defina:rectrrectl RRRRR ,

3. Aplique:

Pr= R(Pl - T)

'

'

'

z

y

x

llpR

'

'

'

''

z

y

x

z

flp

3. Aplique:

'

'

'

z

y

x

ddpR

'

'

'

''

z

y

x

z

fdp

Reconstrução

Reconstrução por triangulação

r

reye

l

leye

lp rp

rP

rrll pRpn

lapnc

rTb pR

rlr

rll cba eyenpRp

)( rll

lrr eyePRP

rleye

Reconstrução por triangulação

z

lr

y

lr

x

lr

zzrrllz

yyrrlly

xxrrllx

c

b

a

np

np

np

eye

eye

eye

pR

pR

pR

)(

)(

)(

npP2

ca ll l

Tlwl

lww PRPRP

rrll pRpn

rlr

rll cba eyenpRp

Outro processo de reconstrução

Z

Y

X

osf

osf

w

v

u

yy

xx

100

0

0

10100

00

00

Z

Y

X

osf

osf

w

v

u

yy

xx

Miguel Ribo, Axel Pinz, Anton L. Fuhrmann “A new Optical Tracking System for Virtual and Augmented Reality Applications”,

ll P0IKp 3

lr PTRKp

l

l

l

l

l P

p

p

p

p

3

2

1

l

r

r

r

r P

p

p

p

p

3

2

1

44

32

31

32

31

44

xrr

rr

ll

ll

x

pvppuppvppup

B

r

r

l

l

0BP l

ll

lllu

Pp

Pp3

1

013 lllll u PpPp

Reconstruction up to a Scale Factor

• Assume that intrinsic parameters of both cameras are known• Essential Matrix is known up to a scale factor (for example, estimated from the 8 point algorithm).

Steve Seitz, University of Washington

Reconstruction up to a Scale Factor

Rtk

T TT tRRtk 2 Tttk 2

22222

22222

22222

YXZYZX

ZYZXYX

ZXYXZY

TTkTTkTTk

TTkTTkTTk

TTkTTkTTk

TraceT 222222 22 tkTTTk ZYX

tk R

t

tkR

t

tk

sgnsgn

TT ttEE ˆˆˆˆ

2

2

2

ˆ1ˆˆˆˆ

ˆˆˆ1ˆˆ

ˆˆˆˆˆ1

ZZYZX

ZYYYX

ZXYXX

TTTTT

TTTTT

TTTTT

Rtk ˆsgn E

Steve Seitz, University of Washington

Reconstruction up to a Scale Factor

T

T

T

E

E

E

E

3

2

1

ˆ

ˆ

ˆ

ˆ

T

T

T

R

R

R

R

3

2

1

Let 3,2,1 ,ˆˆ itEw ii

It can be proved that

2133

1322

3211

wwwR

wwwR

wwwR

Steve Seitz, University of Washington

Reconstruction up to a Scale Factor

We have two choices of t, (t+ and t-) because of sign ambiguityand two choices of E, (E+ and E-).

This gives us four pairs of translation vectors and rotation matrices.

Steve Seitz, University of Washington

Reconstruction up to a Scale Factor

Given and E t

1. Construct the vectors w, and compute R2. Reconstruct the Z and Z’ for each point3. If the signs of Z and Z’ of the reconstructed points are

a) both negative for some point, change the sign ofand go to step 2.

b) different for some point, change the sign of each entryof and go to step 1.

c) both positive for all points, exit.

t

E

pRfRx

tRfRxfZ T

T

13

13

pfRxR

tfRxRfZ T

T

13

13

Steve Seitz, University of Washington

Proposed system: equipament

2 cameras and 1 projector(fast)

1 moving camera and 1 projector(slow)

Proposed system: 32rgb-BCSL coding

left right

positiveslide

positiveslide

negativeslide

Where is a point in the other image?

u u

One solution: (u,v) coordinates

double the number of photos!

Epipolar geometry

l

eyel

P

r

eyer

Pl

plxcl

ycl

zcl

xcr

ycr

zcr

pr

Pr

el er

EpipolarLine

EpipolarLine

Epipolar correspondence

100

0

0

0

0

0

100

0

0

ylyl

l

xlxl

l

x

rly

rl

x

rlz

rl

y

rlz

rl

lrlrlr

lrlrlr

lrlrlr

T

yryr

r

xrxr

r

osf

osf

osf

osf

eyeeye

eyeeye

eyeeye

kkjkik

kjjjij

kijiii

F

0lT

r pFp

Reconstruction by triangulation: ideia

r

reye

l

leye

lp rp

rP

rrll pRpn

lapnc

rTb pR

rlr

rll cba eyenpRp

)( rll

lrr eyePRP

rleye

Reconstruction by triangulation: algebra

z

lr

y

lr

x

lr

zzrrllz

yyrrlly

xxrrllx

c

b

a

np

np

np

eye

eye

eye

pR

pR

pR

)(

)(

)(

wpP2

ca ll l

Tlwl

lww PRPRP

rlr

rll cba eyenpRp

Captured data

Alinhamento de nuvens de pontos

Problema de alinhamento

• Dadas duas nuvens de pontos P e Q, encontrar o movimento rígido (R, t) que minimiza o erro de ajuste

onde q(pi) é o ponto em Q correspondente ao ponto pi de P.

• Dificuldade: não se sabe, a priori, qual é o ponto em Q que corresponde a pi

2||)()(||min tRppq ii

i

Algoritmo ICP (Iterative Closest Point)

• inicia com estimativa grosseira para R e t

• a cada iteração, q(pi) é escolhido como o ponto em Q mais próximo de Rpi + t

• R e t são recalculados de modo a minimizar o erro de ajuste

P. J. Besl and N. D. McKay, A Method for Registration of 3-D Shapes, IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. 14, No. 12, February 1992

Algoritmo ICP (Iterative Closest Point)

Algoritmo ICP (Iterative Closest Point)

• Também se aplica a ajuste de nuvens de pontos a modelos (por exemplo, cilindros)

Subproblema fundamental: alinhamento de pontos correspondentes

• Dados os pares de pontos correspondentes pi, qi (i = 1, ..., n), determinar R e t que minimizam

• Resolvido por Horn[88, 89]

2||)(||min tRpq ii

i

Alinhamento de pontos correspondentes

– Obter os centróides p0 e q0

– Definir pi’= pi – p0 , qi’= qi – q0

– Definir

onde

2||)(||min tRpq ii

i

zzzyzx

yzyyyx

xzxyxx

SSS

SSS

SSS

M

etc,, ,,,,iy

iixxyix

iixxx qpSqpS

Alinhamento de pontos correspondentes

– Rotação: R = M(MT M) –1/2

– Translação: t = q0– Rp0

– Observação: se a transformação é escrita na forma R(p + t), R não muda e t é simplesmente q0– p0

2||)(||min tRpq ii

i

Alinhamento de pontos correspondentes

– Pode ser incluído um fator de escala:

– Rotação: R = M(MT M) –1/2

– Escala:

– Translação: t = q0– sRp0

2||)(1

||min tRpsqs

ii

i

2,2, ||||/|||| i

ii

i pqs

Exemplo (2D)

131219

121170

208136

,

11478

11418

21820

qp

Exemplo (2D)

131219

121170

208136

,

11478

11418

21820

qp

)33.153,175(),66.148,66.38( 00 qp

0.94600.3241

0.3241-0.9460R

( = 18o)

Exemplo (2D)

131219

121170

208136

,

11478

11418

21820

qp

131219

121170

208136

,

3.1334.223

8.1137.166

87.2128.134

qtRp

Variantes do ICP

• ICP funciona melhor quando todo ponto em P corresponde a algum ponto em Q (não é o caso no alinhamento de malhas obtidas por fotografia 3D)

• Diversas estratégias para melhorar o desempenho para alinhamento de malhas com superposição parcial

Variantes do ICP

• Tomar pontos em ambas as malhas

• Usar outros critérios para casamento (textura, normal)

• Atribuir pesos às correspondências

• Rejeitar outliers• Alterar a medida de distância (ponto-

superfície no lugar de ponto-ponto)

Cylinder model

n

jpj

n

j

Tpjpj tensor

nn 00

)(1

))((1

cpcpcpM

zzzyzx

yzyyyx

xzxyxx

zyx

z

y

x

z

y

x

tensorM

1e

2e

3e

321321 ),(ˆˆˆ Meee rseigenvecto

n

jjp n 0

1pc

axis of the points pi:

covariance matrix:

centroid:

Initial cylinder adjustment

pc 2e

3e

cc

1ez c

3ˆ)( ecc aRpc

)( . ˆmin 3 pjja cpe a

tangent plane perpendicular to ê3:

first guess for cc:

first guess for zc:

Results of the initial cylinder adjustment

Projection of a point on a cylinder

0p

p

cz

ppip

Plane : 0 czpp

Axis : ctt zpp 0)(

cc

cit zz

zpp

0

ciii tt zppp 0)(

i

ip

R

pp

ppp

Given : Rc ,,, 0 zpp

Compute : pp

Results

Direct measure

http://vision.middlebury.edu/stereo/data/scenes2006/