Fundamental Matrices

TECHNICAL FIELD

The present invention relates to the field of 3D reconstruction and computer vision in general, and to the field of fundamental matrix estimation schemes in particular.

BACKGROUND ART Introduction

The fundamental matrix (FM) describes the geometric constraints that exist between two images of the same scene. These geometric constraints are called the epipolar constraint. If the pose and the intrinsic parameters of the two cameras that took the two images are known, the FM can be computed directly. Computer vision, however, solves the inverse problem by first estimating the FM using a set of corresponding points identified on the two images and then deducing the parameters of the two cameras. The determination of the FM has several applications in computer vision, photogrammetry, robotics and computer graphics. These applications include camera self- calibration, 3D scene modeling from two to or more images, image rectification, robot vision, Visual Simultaneous Localization And Mapping (V-SLAM), and image-based rendering techniques such as view interpolation, view morphing and other transfer methods.

Epipolar Geometry

The mathematical background of the FM is described next. Let / and / ' be two perspective images of the same scene taken from different viewpoints. The epipolar constraint specifies that there exists a 3 x 3 rank-2 matrix F such that if p and p' are the image coordinates of any 3D point P on respectively / and / ' , then

p ^r Fp' = 0. (1 ) where p = (p ^r l) . The matrix F is called the Fundamental Matrix (FM) and the pair of points M = (p, p') is referred to as a correspondence. Figure 1 shows an example correspondence (p, p') obtained by projecting the 3D point P on / and / ' , respectively. The quantity d _A (M) = d _A (ρ, ρ') = ρ' Fp' is referred to as the algebraic distance. Simple algebraic manipulations show that d _A (M) can be written in the following dot-product form: d _A (M) = m' f.

where m = (xx ' xy ' x yx ' yy ' y x ' y ' l) ⁷ is called the measurement vector 1 and f = ( , f ₂ . .. / ₉ ) ⁷ is the 9-vector containing the 9 entries of F in row-major order.

Background of Fundamental Matrix Estimation

In theory, seven correspondences between / and / ' are generally enough to estimate F . In practice, we are given a number n > 7 of noisy correspondences {M, = (p, , p' )}" ₌₁ · In this case, there is no a 3 x 3 rank-2 matrix that exactly satisfies the epipolar constraint (1 ) for all these noisy correspondences. So, we select the FM F that best fits these correspondences by solving an optimization problem of the following form:

n

argmin C (f ) =∑r, ² ,

subject to D(f ) = 0. where r. measures the residual of the correspondence M, as a function of the FM F(f ) , and D(f ) is the determinant of F(f ) :

D(f) = det F(f)

Different FM estimation schemes may differ in the particular definition they adopt for the residual function r _i , the way in which they enforce the rank-2 constraint, and/or the numerical scheme they use to solve the optimization problem in Eq. (2). Ideally, the residual function r should measure the Euclidean distance d in 4D space between the noisy correspondence M, and the hyper-surface S implicitly defined by the epipolar constraint (1). The Euclidean distance d _E is also referred to as the reprojection error. The rank-2 matrix ¥ _ML that minimizes is the maximum-likelihood (ML) estimate and the schemes that carry out this optimization, such as bundle adjustment (BA) schemes, are regarded as the gold standard (GS). Since they are relatively expensive, GS techniques mainly serve as a benchmark against which other approximate schemes are evaluated.

A more efficient, yet accurate alternative to d _E is Sampson distance d _s . It has the following closed-form expression:

= _dA ² _i l {l ² + l ₂ ² + l[ ² + l ' ² ) where w _Si is the Sampson weighting, I _} is the j th component of the epipolar line 1 = Fp , , and / ' is the j th component of the epipolar line Γ = F ^r p, . Sampson distance d _s ² has the advantage of being a first-order approximation of .

Previous Estimation Schemes

Numerous schemes have been proposed for estimating the FM from a set of noisy correspondences. Prominent examples are the (normalized) eight-point scheme (8P), the iteratively re-weighted eight-point scheme (W8P), the Levenberg-Marquardt (LM) iterative minimization of the Sampson distance (LMS), Hartley's LM-based algebraic distance minimization scheme (ADMS), the renormalization scheme (RNS), the fundamental numerical scheme (FNS), the heteroscedastic-errors-in-variable scheme (HEIV), the constrained FNS (CFNS), and the extended FNS (EFNS). Most of these schemes suffer from one or more of the following issues:

1. Implementation difficulty: Apart from the 8P and W8P schemes, most FM estimation schemes are relatively difficult to implement and understand for non-experts.

2. Improper incorporation of the singularity constraint: Some schemes do not take the singularity constraint into account while estimating the FM. Among these are the 8P and W8P schemes, FNS, HEIV, and RNS. An additional post-processing step is performed to correct the obtained estimate F ₃ so that it satisfies the singularity constraint. This is performed either by the singular value decomposition (SVD) correction or by iteratively correcting the estimates obtained by FNS or HEIV. Whereas SVD-correction gives inaccurate results, iterative-correction takes longer to run and makes the FM estimation process more complicated.

3. Minimization of non-geometric cost functions: Some schemes minimize the algebraic cost function to gain more speed at the expense of lower accuracy. From accuracy point of view, however, it was experimentally and mathematically found that the switch from a geometric to an algebraic cost function has a smaller effect on accuracy than that incurred by the inappropriate incorporation of the singularity constraint. 4. Slow convergence: All FM estimation schemes run several times slower than the 8P scheme. For instance, found in an FM estimation experiment involving « = 100 correspondences that FNS (followed by SVD correction), FNS (followed by iterative correction), CFNS, LMS, and EFNS take 13, 15, 25, 22, and 37 times, respectively, longer to run than the 8P scheme.

5. Lack of robustness: Apart from LM schemes, all the FM estimation schemes that take the singularity constraint into account (such as EFNS and CFNS) produce poor results with data sets contaminated with outliers and/or badly-located data. This is because they are non-robust least-squares techniques that consider all data, including outliers and badly-located ones, equally important.

DISCLOSURE OF THE INVENTION

It is an object of this invention to estimate the FM accurately and efficiently. This is achieved through a novel approach for taking the rank-2 constraint into consideration during the estimation of FMs using any of the popular 8P or W8P schemes. The method allows the minimization of either the algebraic cost function d _A (as in the 8P scheme) or the weighted algebraic cost function d _w = wd _A (as in the W8P scheme). For ease of reference, the two schemes covered under the method of this invention will be referred to as the Extended Eight-Point (E8P) scheme and the Extended Weighted Eight-Point (EW8P) scheme. Given an initial guess ¥ _k of the optimal FM estimate, the present invention replaces the nonlinear determinant function by a Taylor linear approximation at ¥ _k . A better estimate F _{4 +1} is then obtained by solving the resulting linear system of equations, and the process is repeated until convergence. This is in contrast to the 8P and the W8P schemes which totally ignore the singularity constraint.

The present invention has several advantages:

1. It is as simple to understand and implement as the traditional 8P and W8P schemes, \item The E8P scheme was experimentally found to give exactly the same results as ADMS while being almost as fast as the simplest scheme (i.e., the 8P scheme).

2. The EW8P scheme not only improves the accuracy of the W8P scheme, but also runs about twice as fast. Experimental results indicate that it gives near-optimal results (its reprojection error is within 0.1% of the best result most of the time) while being 8-16 times faster than the more complicated schemes such as Levenberg-Marquardt schemes. 3. Similar to the LM scheme and the W8P scheme, the EW8P scheme offers a very flexible structure that permits the use of geometric cost functions and/or robust weighting functions.

4. The FM estimates obtained by the E8P and the EW8P schemes have unit Frobenius norm (after applying the denormalization transformation, the Frobenius norm of the obtained FM estimate may not remain one) and perfectly satisfy the singularity constraint, eliminating the need for a post-processing step to enforce the rank-2 constraint.

BRIEF DESCRIPTION OF DRAWINGS

• Figure 1 shows a flowchart of fundamental matrix estimation according to the present invention.

• Figure 2 illustrates the epipolar geometry. The fact that the projections p and p' are coplanar with the optical centers O _c and O^. gives rise to the epipolar constraint.

BEST MODE FOR CARRYING OUT THE INVENTION

Two preferred embodiments of the present invention are described below. The first preferred embodiment is the Extended Eight-Point (E8P) scheme. The second preferred embodiment is the Extended Weighted Eight-Point (EW8P) scheme. Although the EW8P scheme is not as efficient as the E8P scheme, the EW8P is more robust and accurate.

The Extended Eight-Point Scheme

Derivation

The E8P scheme aims to minimize the same algebraic cost function of the 8P scheme but subject to the additional constraint that the estimated matrix is singular:

argminC, (f) = ∑(mf f) ² ,

f

subject to |f f = 1, (3) D(f) = 0.

Let the measurement matrix 2 M be the n x 9 matrix whose / th row is m, (defined in 1). It follows that C _A (f) = |Mf| ² = f ^r M ^r Mf = f ^r Af . We hereafter refer to the 9 x 9 positive (semi- )definite matrix A as the moment matrix. Let g _{ (i) = f | - 1 and g ₂ (f) = D(f) . We shall refer to g _{ and g ₂ as the constraint functions. It follows that Eq. (3) can be re-written as follows: argminC, (f) = |Mff ,

f

subject to g, (f) = 0, (4) g ₂ V) = 0.

Introducing Lagrange multipliers converts (4) into the following unconstrained optimization problem (using arbitrary non-zero constants to the left of Lagrange multipliers does not alter the optimization problem): argmin L (f , λ, , λ, ) = |Mf f + 2 g _x (f ) + 2^g ₂ (f ).

ΜιΛ

At the optimal point (ί,λ^,λ^ , the partial derivatives ( d _f L , d^L , and d^L ) must vanish. This leads to the following system of equations:

M ^r Mf _f £, (f) + ^d _f £ ₂ (f) = 0,

g,(f) = 0, (5)

Let f _k be a given a guess of the minimizer f . Let be the first-order Taylor series approximation of the constraint function g . at f _k . Then, _; is given as follows:

/, (f) = s, (f _t ) + ¾ (f* )(f - f* )-

The linearization / ₍ can be used to approximate the behavior of g _i at any point f , but the approximation is highly accurate when f is sufficiently close to f _k and when g _i is a low-order polynomial that is free from exponentials, trigonometric functions, or any highly nonlinear functions. Consequently, if the guess f _k is provided sufficiently close to f , we can replace each constraint function g _i in Eq. (5) by its linearization I _i to obtain:

M ^r Mf + λ _ι δ ₍ Ι _ι (ί) + λ ₂ δ _ί 1 ₂ ί(ί) = 0,

/,(f) = 0, (6)

/ ₂ (f) = 0. Substituting with the definitions of into Eq. (6) yields

o,

o,

0. where 9 _f g,(f) = 2f and d _f g ₂ (f) is given as follows:

^ ₂ ⁽ o

Define the 2x9 matrix _k = [d _f g _l (f _k ) d _t g ₂ (f _k )] ^? , the 2x1 vector λ = {λ _χ A ₂ , and the

2 x 1 vector c _k = J _k i _k - (g,(f _t ) g ₂ (f* )/ · The above system of equations can then be compactly rewritten as follows:

M ^r Mf +3 ^T _k A

(7)

= c,.

This is a linear system of 11 equations in the vector of 11 unknowns x Finally,

we can rewrite Eq. (7) in the following more elegant form:

L*x = b, (8) where the 11x11 symmetric matrix h and the 11 -vector are defined as follows: J The linear system of equations (8) can be solved for x using standard elimination techniques such as the LU decomposition. The first nine entries of x give the desired estimate f , whereas the two other entries give the less-interesting value of the Lagrange multiplier λ .

The estimate of f obtained from Eq. (8) may not exactly be equal to the minimizer f of (4) because the linearizations /, may not have approximated the behavior of the original constraint functions g _: at f exactly. Nevertheless, this estimate can be used as a better guess t _{k +1} of the minimizer f , allowing us to obtain more accurate linearizations /, of the constraints g _t and solve Eq. (8) again. Equation (8) should then be written in the following recurrent form:

This process is repeated until we find that subsequent estimates of f are sufficiently close.

Since the E8P approach relies on linear Taylor approximations near the solution, it exhibits a similar behavior to the Newton-Raphson scheme for root-finding: Once the current estimate is sufficiently close to the solution, the convergence becomes highly rapid. This is especially true in our case because the constraints being linearized are low-order polynomials that involve no exponentials or trigonometric functions.

A key advantage of this scheme is that the count of flops it makes in each iteration is independent of n . Indeed, it takes only (2d' ³ 1 3) flops to solve a linear system such as (9) every iteration, where d ' = 11 . Yet, a more efficient update procedure can be obtained if the 9 x 9 moment matrix A = M ^r M is invertible. This procedure further reduces the count of flops performed in each iteration to just (2d ² ) , where d = 9 .

The more efficient procedure can be derived as follows. If A is invertible, we can multiply both sides in Eq. (7) by A ^-1 to obtain the following formula for f _+i in terms of λ :

Substituting Eq. (10) into Eq. (7) yields:

-J _A A-'J^ = c, (1 1 )

Define the 9 x 2 matrix T _t = A ^"1 J ^T _k and the 2 x 2 symmetric matrix Ν _λ = J _k T _k - 3 _k A ^~l J ^T _k .

We can solve Eq. (1 1) by multiplying both sides by - ^ ¹ to obtain λ = -(N^c _A ) . Substituting into Eq. (10) yields the following more efficient update equation:

^ = ¾ (Ν¾ ). (12) An efficient and robust implementation should be able to switch between the two alternative solution procedures (9) and (12) based on the rank of M . When it discovers that M has rank 9, it should calculate A ^-1 once in a pre-processing step and use Eq. (12) to update the solution at each iteration. If it discovers that M does not have rank 9, it should calculate A once in a preprocessing step and use Eq. (9) to advance the solution. We describe one such implementation at the end of this section.

The E8P scheme requires an initial estimate f ₀ of the minimizer f . As with the traditional 8P scheme, we can use the right singular vector v ₉ corresponding to the smallest singular value d ₉ of the measurement matrix M . Among all unit vectors, v ₉ is the minimizer of |Mf| ² when the singularity constraint is ignored. With this particular choice, it takes the E8P scheme only 3-4 iterations to converge in most of the cases. We have also considered different choices of the initial estimate f ₀ . When we used the estimate obtained by running the 7-point scheme on a random subset of 7 correspondences, it took our technique about 5-6 iterations to converge. In all our experiments, the E8P scheme never diverged to a bad solution. This indicates that our technique converges rapidly and reliably to the correct solution whenever provided with a reasonable initial estimate f ₀ .

Implementation

We present an efficient implementation of the E8P scheme 4 in a MATLAB-like pseudocode. The implementation switches to the more efficient iterative formulation (12) if it discovers that the rank of M is 9. This will almost always be the case in practice because we are provided with n » 8 noisy correspondences. It is assumed that we are given as input a set of correspondences C = {M, })Z" , a convergence constant δ , a maximum number of iterations max, and an optional initial estimate f ₀ . The coordinates of the input correspondences can optionally be normalized in a pre-processing step to improve the numerical conditioning of the calculation. This is done by translating and scaling the points in each in each images so that their center becomes the origin and the average of the distances of the points to the origin becomes

M <— BuildMeasurementMatrix( C );

(d, V) <- SVD(M );

if no f ₀ was provided, set f ₀ — V(:, 9) ; if n > 8 AND d(9) / d(l) > 1 x 1(Γ ¹⁰

Precalculate A ^"1 <— V*diag(d(l) ^~2 ,d(2) ^~2 ,...,d(9) ^~2 )* V ⁷ ;

Use update equation (12) to advance the solution;

Stop when max iterations are exhausted or when |f _A+1 - ≤ δ ² ; else

Precalculate A = M ⁷ M ;

Use update equation (9) to advance the solution;

Stop when max iterations are exhausted or when -f _k ≤ δ ² ; r ^eturn Lnv ;

The Extended Weighted Eight-Point (EW8P) Scheme

The (traditional) weighted eight-point scheme (W8P) attempts to minimize a cost function of the following form:

subject to ]f I = 1. where each weight w _j depends on both the correspondence M, and the FM f . The scheme carries out the optimization in the following iterative fashion. At iteration k , the FM estimate f _k is used to evaluate the weight w _i for each correspondence M,.. The factor w , is then used to weight the corresponding row m ⁷ in the measurement matrix M . The optimization problem can then be written as follows:

arg - 1).

where the nx9 matrix M _t =W _A M and W _t is the nxn diagonal matrix diag(w,,w ₂ ,...,w _n ) . The W8P scheme can be made to minimize the Sampson geometric cost function by setting w, =w _Si =l/|Vt^,|. The scheme that results when this particular form of w is used is called the Sampson scheme. It can also be made to minimize a robust geometric cost function by setting w =w _Si w _Hi where w _Hi is the robust Huber weighting that helps lower the effect of badly-located correspondences on the estimated FM.

Similar to the 8P scheme, the W8P can be extended to directly incorporate the rank-2 constraint into the optimization. We refer to the resulting scheme as the extended weighted eight-point (EW8P) scheme. It differs from the E8P scheme in that the measurement matrix M has to be re- weighted in each iteration k by the weight matrix W _t , yielding M. _k = W _A M . Consequently, the update equation (9) changes into: or shortly into

'k -

The implementations of the EW8P and the 8P schemes differ in two ways. First, the EW8P scheme calculates in each iteration a set of different weights and subsequently a new moment matrix A _k = M. ⁷ _k M. _k . Second, the EW8P scheme does not use the solution procedure given in

Eq. (1 1) because it needs to compute a different inverse A^ ¹ every iteration which wastes any possible gain in efficiency. Apart from these two differences, the EW8P and the 8P schemes have quite similar implementations.

The pseudocode implementation of the EW8P scheme 6 is given below. It takes the same input as the E8P implementation in addition to the weighting function w (M, f ) .

Input: C = {M, }) " , δ , max, f ₀ , w

Steps:

M <— BuildMeasurementMatrix( C );

%Initialize f with v ₉ if no f ₀ was provided.

if f _n is NIL

(d, V) <- SVD( M );

f <- V(:, 9) ;

else

f <- f ₀ ;

k - 0 ;

H <- ^• zeros(l 1, 1 1) ;

b <- zeros(l 1, 1) ; repeat

<- f ;

\textrm{\%Calculate the weights}

w^( _w (f _oW ,M,) w(f _old ,M ₂ ) ... w(f _M ,M _n )f ;

H(l:9,l:9) -M^ *M _k ;

J - [2 * t _o ^T _id ; NablaDet(f _oM ) ^r ];

H(l : 9,10 : 11) = J ⁷ ;

H 10 : 11,1 : 9) = J ;

%Solve Hx = b

x<-H\b;

f<-x(l:9);

k <- k +1;

until k > max OR |f - f _M f≤ δ ²

return f ;

Unlike the E8P scheme, the count of flops made in each iteration of the EW8P scheme is no more independent of n . This is due to the need to calculate a set of different weights and subsequently calculate a new moment matrix A _k = M^M^ every iteration. This scheme may not, consequently, be as efficient as the E8P scheme. Practical experiments, however, show that the EW8P scheme runs many times faster than the other more complicated schemes (such as FNS, EFNS, LM). It is even faster than the less-accurate W8P scheme and Hartley's ADMS (as long as n is not very large).

While the present invention has been described with respect to certain preferred embodiments, it will be apparent to those skilled in the art that various changes and modifications may be made without departing from the scope of the invention.