Transcript Multiple View Geometry in Computer Vision

Singular Value Decomposition
Amn  UmmΣmn VnTn
 1 0  0 
0   0
2



  
Σ

0
0


n






 0 0  0 
mn
1   2     n  0
UT U  I
VT V  I
A  U1 1 V1T  U2  2 V2T   Un  n VnT
UΣ
Σ VT X
Singular Value Decomposition
A  UΣ V
T
• Homogeneous least-squares
min AX subject to X  1
• Span and null-space
S L  U1 U2 ; N L  U3U4 
S R  V1V2 ; N R  V3V4 
solution X  Vn
 1 0
0 
2
Σ
0 0

0 0
0
0
0
0
0
0
0

0
• Closest rank r approximation
~
~ T ~
A  UΣ
UΣ V   diag1, 2 ,, r , 0 ,, 0 
• Pseudo inverse
A   VΣ  U T   diag11, 21,, r1, 0 ,, 0 
Singular Value Decomposition
A  UΣ V
T
UΣ
Σ VT X
Homogeneous least-squares
min AX subject to X  1
solution X  Vn
Parameter estimation
• 2D homography
Given a set of (xi,xi’), compute H (xi’=Hxi)
• 3D to 2D camera projection
Given a set of (Xi,xi), compute P (xi=PXi)
• Fundamental matrix
Given a set of (xi,xi’), compute F (xi’TFxi=0)
• Trifocal tensor
Given a set of (xi,xi’,xi”), compute T
Number of measurements required
• At least as many independent equations
as degrees of freedom required
• Example:
x'  Hx
 x   h11
λ  y  h21
 1  h31
h12
h22
h32
h13   x 
h23   y 
h33   1 
2 independent equations / point
8 degrees of freedom
4x2≥8
Approximate solutions
• Minimal solution
4 points yield an exact solution for H
• More points
• No exact solution, because
measurements are inexact (“noise”)
• Search for “best” according to some cost
function
• Algebraic or geometric/statistical cost
Gold Standard algorithm
• Cost function that is optimal for
some assumptions
• Computational algorithm that
minimizes it is called “Gold
Standard” algorithm
• Other algorithms can then be
compared to it
Direct Linear Transformation
(DLT)
xi  Hxi  0 xi  Hxi
 h 1T x 
 T i
T
xi  xi, yi, wi  Hxi   h 2 x i 
 3T 
 yh 3T x  wh 2 T x   h x i 
i
i
 i T i
1
3T

x i  Hx i  wih x i  xih x i 


2T
1T
 xih x i  yih x i 


 0T

T

 wi x i
 yix iT

T

 wi x i
T
1



yi x i
h 
 2
T
T
0
 xix i  h   0

T
T
3


xi x i
0  h 
Ai h  0
•
Direct Linear Transformation
(DLT)
Equations are linear in h
Ai h  0
• Only 2 out of 3 are linearly independent
(indeed, 2 eq/pt)
 0T

T

w
x
 i i
 yix iT

 wix iT
0T
xix iT
yix iT  h1 



 xix iT  h 2   0
0T  h 3 
xiA1i  yiAi2  wiA3i  0
 0T
 T
 wix i
 wix iT
0T
1


h
yix iT  2 
h   0
T
 xix i  3 
h 
(only drop third row if wi’≠0)
•Holds for any homogeneous representation,
e.g. (xi’,yi’,1)
Direct Linear Transformation
(DLT)
• Solving for H
Ah  0
 A1 
A 
 2 h  0
A 3 
 
A 4 
size A is 8x9 or 12x9, but rank 8
Trivial solution is h=09T is not interesting
1-D null-space yields solution of interest
pick for example the one with h  1
Direct Linear Transformation
(DLT)
• Over-determined solution  A1 
Ah  0
A 
 2 h  0
  
 
A n 
No exact solution because of inexact measurement
i.e. “noise”
Find approximate solution
- Additional constraint needed to avoid 0, e.g.
- Ah  0 not possible, so minimize Ah
h 1
DLT algorithm
Objective
Given n≥4 2D to 2D point correspondences {xi↔xi’},
determine the 2D homography matrix H such that xi’=Hxi
Algorithm
(i) For each correspondence xi ↔xi’ compute Ai. Usually
only two first rows needed.
(ii) Assemble n 2x9 matrices Ai into a single 2nx9 matrix A
(iii) Obtain SVD of A. Solution for h is last column of V
(iv) Determine H from h
Inhomogeneous solution
Since h can only be computed up to scale, ~
pick hj=1, e.g. h9=1, and solve for 8-vector h
0
0
 xi wi '  yi wi '  wi wi ' xi yi ' yi yi ' ~   wi yi ' 
 0

h  
x w ' y w ' w w '

0
0
0
xi xi ' yi xi '
i i
i i
 i i
 wi xi ' 
Solve using Gaussian elimination (4 points) or
using linear least-squares (more than 4 points)
However, if h9=0 this approach fails
also poor results if h9 close to zero
Therefore, not recommended
Note h9=H33=0 if origin is mapped to infinity 0
l T Hx 0  0 0 1H 0  0
1
Degenerate configurations
x1
x4
x2
x3
Constraints:
Define:
Then,
x1
H?
x4
x2
x3
(case A)
xi  Hxi  0
H’?
(case B)
x1
x4
x2
x3
i=1,2,3,4
H*  x4lT
H*xi  x4 lT x i  0, i  1,2,3
H*x 4
 
 x l x   kx
T
4
4
4
H* is rank-1 matrix and thus not a homography
If H* is unique solution, then no homography mapping xi→xi’(case B)
If further solution H exist, then also αH*+βH (case A)
(2-D null-space in stead of 1-D null-space)
Solutions from lines, etc.
2D homographies from 2D lines
li  HT li
Ah  0
Minimum of 4 lines
3D Homographies (15 dof)
Minimum of 5 points or 5 planes
2D affinities (6 dof)
Minimum of 3 points or lines
Conic provides 5 constraints
Mixed configurations?
Cost functions
• Algebraic distance
• Geometric distance
• Reprojection error
• Comparison
• Geometric interpretation
• Sampson error
Algebraic distance
DLT minimizes
e  Ah
ei
Ah
residual vector
partial vector for each (xi↔xi’)
algebraic error vector
 0
 wix
d alg xi , Hxi   ei  
T
T


w
x
0
 i i
algebraic distance
T
2
2
T
i
 yix
 xix
T
i
T
i

h

2
d alg x1 , x 2   a12  a22 where a  a1, a2 , a3 T  x1  x 2
2
2



d
x
,
Hx
 alg i i   ei
i
2
 Ah  e
2
2
i
Not geometrically/statistically meaningfull, but given good
normalization it works fine and is very fast (use for initialization)
Geometric distance
x measured coordinates
xˆ
x
estimated coordinates
true coordinates
d(.,.) Euclidean distance (in image)
Error in one image
e.g. calibration pattern
2
ˆ
H  argmin d x , Hx 
H

i
i
i
Symmetric transfer error
2
-1
ˆ
H  argmin  d x , H x    d x  , Hx
i
H
i
i
2

i
i
Reprojection error
Hˆ , xˆ , xˆ    argmin d x , xˆ   d x, xˆ  
2
i
i
H, xˆ i , xˆ i
i
i
i
ˆ xˆ
subject toxˆ i  H
i
2
i
i
Reprojection error


d x, H x  d x, Hx
-1
2
2
2


ˆ
ˆ
d x, x   d x , x 
2
Comparison of geometric
and algebraic distances
Error in one image
T
T








ˆ
ˆ
ˆ
ˆ
xi  xi , yi , wi  xi  xi , yi , wi   Hx
 yiwˆ i  wi yˆ i 

A i h  ei  
 wixˆi  xiwˆ i 
 0T
 T
 wix i
 wix iT
0T
dalg xi , xˆ i    yiwˆ i  wi yˆi   wixˆi  xiwˆ i 
2
2

1


h
T

yix i  2 
h 
T
 xix i  3 
h 
2
d xi , xˆ i    yi / wi  yˆ i / wˆ i   xˆi / wˆ i  xi / wi 
 d alg xi , xˆ i  / wiwˆ i
2
2

2 1/ 2
ˆ i  h3xi , except for affinities
wi  1 typical, but not w
For affinities DLT can minimize geometric distance
Possibility for iterative algorithm
Geometric interpretation of
reprojection error
Estimating homography~fit surface ν Hto points X=(x,y,x’,y’)T in 4
xi  Hxi  0 represents 2 quadrics in 4 (quadratic in X)
ˆ
Xi  X
i
2
  xi  xˆi    yi  yˆ i   xi  xˆi    yi  yˆ i 
2
2
2
2
2
2
 d x i , xˆ i   d xi , xˆ i 
d xi , xˆ i   d xi , xˆ i   d Xi , H 
2
2
2
Analog to conic fitting
dalg x, C  x T Cx
2
d x, C
2
Sampson error
between algebraic and geometric error
2
ˆ is
ˆ that minimizes the geometric error X  X
Vector X
the closest point on the variety ν H to the measurement X
ˆ
Sampson error: 1st order approximation of X
Ah  CH X  0
CH
CH X  δ X   CH X  
δX
X
C
CH X   H δ X  0
X
ˆ X
δX  X

ˆ 0
CH X
CH
Jδ X  e with J 
X
Find the vector δ X that minimizes δX subject to Jδ X  e
Find the vector δ X that minimizes δX subject to Jδ X  e
Use Lagrange multipliers:
minimize
δX δX - 2λJδX  e  0
T
derivatives 2δX - 2λ J  0
T
T
2Jδ X  e  0
 δX  J T λ
 JJ T λ  e  0
 e
 J JJ 
 λ   JJ
 δX
ˆ  Xδ
X
X
δX
2
T
 
 δ TX δ X  e T JJ T
T 1
1
e
T 1
e
Sampson error
between algebraic and geometric error
2
ˆ is
ˆ that minimizes the geometric error X  X
Vector X
the closest point on the variety ν H to the measurement X
ˆ
Sampson error: 1st order approximation of X
Ah  CH X  0
CH
CH X  δ X   CH X  
δX
X
C
CH X   H δ X  0
X
ˆ X
δX  X

ˆ 0
CH X
Jδ X  e
Find the vector δ X that minimizes δX subject to Jδ X  e
δX
2
 
 δ δ  e JJ
T
X X
T
T 1
e
(Sampson error)
Sampson approximation
δX
2
 
 e JJ
T
T 1
e
A few points
(i)
(ii)
(iii)
For a 2D homography X=(x,y,x’,y’)
e  CH X  is the algebraic error vector
C
J  H is a 2x4 matrix,
X e.g. J    wx T h 2  yx T h3 / x  wh  yh
11
i i
i i
i 21
i 31


T
(iv) Similar to algebraic error e  e e
2
(v)
2
in fact, same as Mahalanobis distance e JJ T
Sampson error independent of linear reparametrization
(cancels out in between e and J)
  
T
T 1
e JJ
e
(vi) Must be summed for all points
(vii) Close to geometric error, but much fewer unknowns
Statistical cost function and
Maximum Likelihood Estimation
• Optimal cost function related to noise model
• Assume zero-mean isotropic Gaussian noise
(assume outliers removed)
Pr x  
1
2 πσ
2
e

 d  x, x 2 / 2 2

Error in one image
Pr xi | H   
i
1
2 πσ
log Pr xi | H   

xi ,Hx i 2 / 2 2 
e
2
d
1
2σ
2
2

 d x i , Hx i   constant
Maximum Likelihood Estimate
 d xi , Hxi 
2
Statistical cost function and
Maximum Likelihood Estimation
• Optimal cost function related to noise model
• Assume zero-mean isotropic Gaussian noise
(assume outliers removed)
Pr x  
1
2 πσ
2
e

 d  x, x 2 / 2 2

Error in both images
Prxi | H   
i
1
2 πσ
2
e
 d x i , x i 2  d

Maximum Likelihood Estimate
2
2


ˆ
ˆ




d
x
,
x

d
x
,
x
 i i
i
i
xi ,Hx i 2  / 2 2 
Mahalanobis distance
• General Gaussian case
Measurement X with covariance matrix Σ
XX
2

 X  X   1 X  X 
T
Error in two images (independent)
2
2
X  X   X  X 
Varying covariances
X
i
i
 Xi
2
i
 Xi  X i
2
i
Invariance to transforms ?
x   Hx
~
x  Tx
~
x   T x 
?
~
H  T  HT
-1
~~
~
x  Hx
~
Tx  HT x
-1 ~
x  T HT x
will result change?
for which algorithms? for which transformations?
Non-invariance of DLT
Given x i  xi and H computed by DLT,
xi  T xi , ~
xi  Txi
and ~
~
~
Does the DLT algorithm applied to x i  xi
~
-1
yield H  THT ?
Effect of change of coordinates on algebraic error


~~
~
~

ei  xi  Hxi  Txi  THT-1 Txi  T* xi  Hxi   T*ei
for similarities
0
 R
sR t 
*
T  s  T
T  


t
R
s


 0 1
~~
T
T
~
~
so A i h  e1 , e2   sR e1 , e2   s A i h

~~
~


dalg x i , Hxi   sd alg x i , Hx i

Non-invariance of DLT
Given x i  xi and H computed by DLT,
xi  T xi , ~
xi  Txi
and ~
~
~
Does the DLT algorithm applied to x i  xi
~
-1
yield H  THT ?
2

minimize d alg x i , Hxi  subject to H  1
i
~~ 2
~
 minimize d alg xi , Hx i subject to H  1
i
~~ 2
~
~

 minimize d alg x i , Hx i subject to H  1


i


Invariance of geometric error
Given x i  xi and H,
~
-1
~
~
~
~




x

T
x
,
x

T
x
,

x

x
,
and i
i
i
i
i H  T HT
i
Assume T’ is a similarity transformations

 

~
d~
xi , H~
x i  d Txi , THT-1T xi  d Txi , THxi 
 sd xi , Hxi 
Normalizing transformations
• Since DLT is not invariant,
what is a good choice of coordinates?
e.g.
• Translate centroid to origin
• Scale to a 2 average distance to the origin
• Independently on both images
Or
Tnorm
0
w / 2
w  h
  0
w  h h / 2 
 0
0
1 
1
Importance of normalization
0
x
 i
0
0  xi  yi  1
yi
1
~102 ~102 1
0
0
0
yixi
 xixi
~102
~102
1
~104
yi yi
 xi yi
~104
orders of magnitude difference!
 h1 
yi  2 
h   0

 xi  3 
h 
~102
Normalized DLT algorithm
Objective
Given n≥4 2D to 2D point correspondences {xi↔xi’},
determine the 2D homography matrix H such that xi’=Hxi
Algorithm
 xi
(i) Normalize points ~
xi  Tnorm x i , ~
xi  Tnorm
(ii) Apply DLT algorithm to ~
xi  ~
xi ,
~
(iii) Denormalize solution H  T-1 H
T
norm
norm
Iterative minimization metods
Required to minimize geometric error
(i) Often slower than DLT
(ii) Require initialization
(iii) No guaranteed convergence, local minima
(iv) Stopping criterion required
Therefore, careful implementation required:
(i) Cost function
(ii) Parameterization (minimal or not)
(iii) Cost function ( parameters )
(iv) Initialization
(v) Iterations
Parameterization
Parameters should cover complete space
and allow efficient estimation of cost
• Minimal or over-parameterized? e.g. 8 or 9
(minimal often more complex, also cost surface)
(good algorithms can deal with over-parameterization)
(sometimes also local parameterization)
• Parametrization can also be used to restrict
transformation to particular class, e.g. affine
Function specifications
(i) Measurement vector XN with covariance Σ
(ii) Set of parameters represented by vector P N
(iii) Mapping f : M →N. Range of mapping is
surface S representing allowable measurements
(iv) Cost function: squared Mahalanobis distance
X  f P    X  f P   1 X  f P 
2
T
Goal is to achieve f P  X, or get as close as
possible in terms of Mahalanobis distance
Error in one image
 d xi , Hxi 
2
f : h  Hx1 , Hx2 ,...,Hxn 
X  f h 
Symmetric transfer error


 d x i , H xi  d xi , Hx i 
i
2
-1
2

f : h  H-1x1 , H-1x2 ,...,H-1xn , Hx1, Hx2 ,...,Hxn
X  f h 
Reprojection error
 d xi , xˆ i   d xi , xˆ i 
2
2
f : h, xˆ 1 , xˆ 2 ,...,xˆ n   xˆ 1 , xˆ 2 ,...,xˆ n 
X  f h 

Initialization
• Typically, use linear solution
• If outliers, use robust algorithm
• Alternative, sample parameter space
Iteration methods
Many algorithms exist
• Newton’s method
• Levenberg-Marquardt
• Powell’s method
• Simplex method
Gold Standard algorithm
Objective
Given n≥4 2D to 2D point correspondences {xi↔xi’},
determine the Maximum Likelyhood Estimation of H
(this also implies computing optimal xi’=Hxi)
Algorithm
(i) Initialization: compute an initial estimate using
normalized DLT or RANSAC
(ii) Geometric minimization of -Either Sampson error:
● Minimize the Sampson error
● Minimize using Levenberg-Marquardt over 9 entries of h
or Gold Standard error:
● compute initial estimate for optimal {xi}
2
2
● minimize cost  d x i , xˆ i   d xi , xˆ i  over {H,x1,x2,…,xn}
● if many points, use sparse method
Robust estimation
• What if set of matches contains gross outliers?
RANSAC
Objective
Robust fit of model to data set S which contains outliers
Algorithm
(i) Randomly select a sample of s data points from S and
instantiate the model from this subset.
(ii) Determine the set of data points Si which are within a
distance threshold t of the model. The set Si is the
consensus set of samples and defines the inliers of S.
(iii) If the subset of Si is greater than some threshold T, reestimate the model using all the points in Si and terminate
(iv) If the size of Si is less than T, select a new subset and
repeat the above.
(v) After N trials the largest consensus set Si is selected, and
the model is re-estimated using all the points in the
subset Si
Distance threshold
Choose t so probability for inlier is α (e.g. 0.95)
• Often empirically
2
d
• Zero-mean Gaussian noise σ then  follows
 m2 distribution with m=codimension of model
(dimension+codimension=dimension space)
Codimension
Model
t2
1
l,F
3.84σ2
2
H,P
5.99σ2
3
T
7.81σ2
How many samples?
Choose N so that, with probability p, at least one
random sample is free from outliers. e.g. p=0.99
1  1  e 
s N
 1 p

N  log1  p / log 1  1  e
s

proportion of outliers e
s
2
3
4
5
6
7
8
5%
2
3
3
4
4
4
5
10% 20% 25% 30% 40% 50%
3
5
6
7
11
17
4
7
9
11
19
35
5
9
13
17
34
72
6
12
17
26
57
146
7
16
24
37
97
293
8
20
33
54 163 588
9
26
44
78 272 1177
Acceptable consensus set?
• Typically, terminate when inlier ratio reaches
expected ratio of inliers
T  1  en
Adaptively determining
the number of samples
e is often unknown a priori, so pick worst case,
e.g. 50%, and adapt if more inliers are found,
e.g. 80% would yield e=0.2
• N=∞, sample_count =0
• While N >sample_count repeat
•
•
•
•
Choose a sample and count the number of inliers
Set e=1-(number of inliers)/(total number of points)
Recompute N from e
N  log1 p/ log1 1 es 
Increment the sample_count by 1
• Terminate
Robust Maximum Likelyhood
Estimation
Previous MLE algorithm considers fixed set of inliers
Better, robust cost function (reclassifies)
e 2 e 2  t 2 inlier
R   ρd i  with ρe   2 2 2
i
t e  t outlier
Other robust algorithms
• RANSAC maximizes number of inliers
• LMedS minimizes median error
• Not recommended: case deletion,
iterative least-squares, etc.
Automatic computation of H
Objective
Compute homography between two images
Algorithm
(i) Interest points: Compute interest points in each image
(ii) Putative correspondences: Compute a set of interest
point matches based on some similarity measure
(iii) RANSAC robust estimation: Repeat for N samples
(a) Select 4 correspondences and compute H
(b) Calculate the distance d for each putative match
(c) Compute the number of inliers consistent with H (d<t)
Choose H with most inliers
(iv) Optimal estimation: re-estimate H from all inliers by
minimizing ML cost function with Levenberg-Marquardt
(v) Guided matching: Determine more matches using
prediction by computed H
Optionally iterate last two steps until convergence
Determine putative
correspondences
• Compare interest points
Similarity measure:
• SAD, SSD, ZNCC on small neighborhood
• If motion is limited, only consider interest points
with similar coordinates
• More advanced approaches exist, based on
invariance…
Example: robust computation
Interest points
(500/image)
Putative
correspondences (268)
Outliers (117)
Inliers (151)
Final inliers (262)
•
Maximum Likelihood Estimation
X
i
i
 Xi
2
i
 Xi  X i
2
i
•
•
DLT not invariant  normalization
Geometric minimization invariant
•
Iterative minimization
•
•
•
•
Cost function
Parameterization
Initialization
Minimization algorithm
Automatic computation of H
Objective
Compute homography between two images
Algorithm
(i) Interest points: Compute interest points in each image
(ii) Putative correspondences: Compute a set of interest
point matches based on some similarity measure
(iii) RANSAC robust estimation: Repeat for N samples
(a) Select 4 correspondences and compute H
(b) Calculate the distance d for each putative match
(c) Compute the number of inliers consistent with H (d<t)
Choose H with most inliers
(iv) Optimal estimation: re-estimate H from all inliers by
minimizing ML cost function with Levenberg-Marquardt
(v) Guided matching: Determine more matches using
prediction by computed H
Optionally iterate last two steps until convergence
Algorithm Evaluation
and Error Analysis
• Bounds on performance
• Covariance estimation
?
uncertainty
?
residual error
Algorithm evaluation
Test on real data or
test on synthetic data
x  measured coordinates
xˆ , H estimated quantities
x, H true coordinates
• Generate synthetic correspondences xi  x i
• Add Gaussian noise, yielding xi  x i
ˆ from xi  x i using algorithm
• Estimate H
maybe also xˆ i , xˆ i 
ˆ x or x  H
ˆx
• Verify how well xi  H
i
i
i
• Repeat many times (different noise, same )
• Error in one image
xi  x i
xi  xi  G(0, σ)
ˆ , then
H
Estimate
eres
1/ 2
 1
2
   d xi , xˆ i  
 2n 1

n
Note: residual error ≠ absolute measure of quality of
ˆ
H
e.g. estimation from 4 points yields eres=0
more points better results, but eres will increase
• Error in two images
xi  x i
x i  x i  G(0, σ2 ) xi  xi  G(0, σ2 )
ˆ , xˆ , xˆ  so that
Estimate H
i
i
1/ 2
1 
2
2

  d x i , xˆ i    d xi , xˆ i  
4n  1
1

n
eres
ˆ xˆ , then
xˆ i  H
i
n
Optimal estimators (MLE)
Estimate expected residual error of MLE,
Other algorithms can then be judged to this standard
f : M → N (parameter space to measurement space)
X N P  M
f P   X
f
X
P
SM
M
N
dimension of submanifold SM = #essential parameters

X  G X, Nσ2
ˆ  MLE(X)
X

Assume SM locally planar around X
projection of isotropic Gaussian distribution on N
with total variance N2 onto a subspace of
dimension s is an isotropic Gaussian distribution
with total variance s2
X
n
X
X
SM
N measurements (independent Gaussian noise 2)
model with d essential parameters
(use s=d and s=(N-d))
(i) RMS residual error for ML estimator
ˆ  X / N
eres  E  X


2
1/ 2
  1  d / N 
1/ 2
(ii) RMS estimation error for ML estimator
2

ˆ
eest  E X  X / N 


1/ 2
  d / N 
1/ 2
X
n
X
X
SM
Error in one image
d  8 and N  2n
eres   1  d / N 
1/ 2
eest   d / N 
1/ 2
1/ 2


eres   1 4 / n
eest   4 / n
1/ 2
Error in two images
d  8  2n and N  4n
 n4
eres   

 2n 1/ 2
n4
eest   

 2n 
1/ 2
Covariance of estimated model
• Previous question: how close is the error
to smallest possible error?
• Independent of point configuration
• Real question: how close is estimated
model to real model?
• Dependent on point configuration
(e.g. 4 points close to a line)
Forward propagation of covariance
Let v be a random vector in M with mean v and
covariance matrix , and suppose that f: M
→N is an affine mapping defined by
f(v)=f(v)+A(v-v). Then f(v) is a random variable
with mean f(v) and covariance matrix AAT.
Note: does not assume A is a square matrix
Example:
 

x  G 0,12 , y  G 0,22

x  f ( x, y)  3x  2 y  7
x   f (0,0)  7
1 0


0
4


A  3 2
1 0 3
AA  3 2
 25



0 4 2
T
std ( x)  5
Example:
 

x  G 0,12 , y  G 0,22

x  3 x  2 y
y  3x  2 y
1 0


0
4


3 2 
A

3

2


3 2  1 0 3 3   25  7
AA  







3

2
0
4
2

2

7
25



 

T
Exx  25; Eyy  25; Exy  7;
Non-linear propagation of covariance
f (v)  f ( v)  J(v  v)
f i
J ij 
v j
Let v be a random vector in M with mean v and
covariance matrix , and suppose that f: M →N
differentiable in the neighborhood of v. Then, up to
a first order approximation, f(v) is a random
variable with mean f(v) and covariance matrix JJT,
where J is the Jacobian matrix evaluated at v
Note: good approximation if f close to linear
within variability of v
Example:
  0  1 0  

x  ( x, y )  G  , 


0
0
4






x  f ( x, y)  x 2  3x  2 y  5
T

x     P( x, y) f ( x, y)dxdy


1 ( x 2  y 2 / 4) / 2 2
P ( x, y ) 
e
8π
 2    P( x, y) f ( x, y)  x 2 dxdy

x  5  σ2
σ2x  25σ2  2σ4
x  5
1   3 
2
2
2
σ x  σ 3  2

25
σ
  2
4

 
J  3  2
Example:
f ( x, y)  ax2  bxy  cy2  dx  ey  f
mean  aσ2x  cσ2y  f
variance 2a2σ4x  bσ2x σ2y  2c2σ2y  dσ2x  eσ2y
est.mean  f
est. variance dσ2x  eσ2y
Backward propagation of covariance
f : M → N
ˆ N Pˆ  M
X
1 ˆ
ˆ
f
X  Pˆ
ηX X


f 1 h X  Pˆ
f -1
X
h
P
X
Backward propagation of covariance
f
1
h X  Pˆ
assume f is affine
what about f -1oh ?
ˆ
minimize: X  X

X
f P   X
f P  f P   JP-P 

f
P
X
h
-1


 X  f Pˆ  X  X   J Pˆ  P
 



solution: P
ˆ  P  J T 1J -1 J T 1 X  X
f
1



ˆ
P  f 1 X Pˆ  f 1 X



T 1 -1 T 1
ˆ
h X  P  J  J J  X  X   f 1 X
hX 
Backward propagation of covariance
f
1
h X  Pˆ

1


X
-1 T
A  J  J J  1
T
f
P
X
h
-1
1

-1 T
1
1
 
1

-1 T
 f 1 h  J  J J   J  J J 
T
 
 J  J 
T
1
-1 T
T
1
-1
T
T

1
 J  J J   J J  J
T

-T
1

T
Backward propagation of covariance
f
1
h X  Pˆ
f
P

If f is affine, then  P  J T  x1J
X
h
-1
X

-1
non-linear case, obtain first order
approximations by using Jacobian
Over-parameterization
In this case f is not one-to-one and rank J<M

1
x

-1
so  P  J  J can not hold
T
e.g. scale ambiguity  infinite variance!
However, if constraints are imposed, then ok.
A

1
x
P  J  J
T

A

1
x

1
 A A J  JA AT
T T
Invert d xd in stead of MxM
Over-parameterization
When constraint surface is locally orthogonal
to the null space of J

1
x
P  J  J
T


(pseudo-inverse)
e.g. usual constraint ||P||=1
||P||=1
Example: error in one image
ˆ from the data
(i) Estimate the transformation H
(ii) Compute Jacobian J f  X / h , evaluated at hˆ
(iii) The covariance matrix of the estimated hˆ is

T

1
given by   J  J
h

x
x iT
1 ~
J i  xi / h  
wi  0

0
~
xT
i
 xi~
x iT 

 yi~
x iT 

T 1 
 h  J  J    J i i J i 
 i


T
1
x



Example: error in both images
J  A | B
separate in homography and point parameters
T 1
T 1

A

A
A
 x B
T 1
x
J  x J   T 1
T 1 
 B  x A B  x B 
Using covariance matrix in point transfer
Error in one image
 x  J h  J
T
h h
Error in two images
x  J h h JTh  J x x JTx
(if h and x independent, i.e. new points)
Example:
=1 pixel =0.5cm
(Crimisi’97)
Example:
=1 pixel =0.5cm
(Crimisi’97)
Example:
(Crimisi’97)
Monte Carlo estimation of
covariance
• To be used when previous assumptions do
not hold (e.g. non-flat within variance) or
to complicate to compute.
• Simple and general, but expensive
• Generate samples according to assumed
noise distribution, carry out computations,
observe distribution of result