Selecting the right essential matrix

Viewed 948

I want to code myself the sfm pipeline using Matlab because I need some outputs that opencv functions don't provide. However, I'm using opencv for comparison.

The Opencv function [E,mask] = cv.findEssentialMat(points1, points2, 'CameraMatrix',K, 'Method','Ransac'); provides the essential matrix solution using Nister's fivepoint algorithm and RANSAC.

the inlier indices are found using :InliersIndices=find(mask>0);

I used this Matlab impelmentation of Nister's algorithm:

Fivepoint_algoithm_code

The call to the function is as follows:

[E_all, R_all, t_all, Eo_all] = five_point_algorithm( pts1, pts2, K, K);

The algorithm outputs up to 10 solutions of essential matrices. However, I encountered the following issues:

  1. The impelmentation stated above is only for perfect correspondances (without Ransac) and I'm providing to the algorithm 5 correspondances using InliersIndices, the outputted essential matrices (up to 10) are all different from the one returned by Opencv.
  2. All the returned essential matrices should be solutions so why when I triangulate for each one using the below function, I don't obtain the same 3D points?
  3. How to choose the right essential marix solution?

I triangulate using the function of matlab toolbox

Projection matrices:

P1=K*[eye(3) [0;0;0]];
P2=K*[R_all{i} t_all{i}];

[pts3D,rep_error] = triangulate(pts1', pts2', P1',P2');

Edit

The returned E from [E,mask] = cv.findEssentialMat(points1, points2, 'CameraMatrix',K, 'Method','Ransac');

E =

    0.0052   -0.7068    0.0104
    0.7063    0.0050   -0.0305
   -0.0113    0.0168    0.0002

For the 5-point Matlab implementation,5 random indices from inliers are taken so:

pts1 =

  736.7744  740.2372  179.2428  610.5297  706.8776
  112.2673  109.9687   45.7010   91.4371   87.8194

pts2 =

  722.3037  725.3770  150.3997  595.3550  692.5383
  111.7898  108.6624   43.6847   90.6638   86.8139

K =

  723.3631    7.9120  601.7643
   -3.8553  719.6517  182.0588
    0.0075    0.0044    1.0000

and 4 solutions are returned:

E1 =

   -0.2205    0.9436   -0.1835
    0.8612    0.2447   -0.1531
    0.4442   -0.0600   -0.0378

 E2 =

   -0.2153    0.9573    0.1626
    0.8948    0.2456   -0.3474
    0.1003    0.1348   -0.0306
E3 =

    0.0010   -0.9802   -0.0957
    0.9768    0.0026   -0.1912
    0.0960    0.1736   -0.0019
E4 =

   -0.0005   -0.9788   -0.1427
    0.9756    0.0021   -0.1658
    0.1436    0.1470   -0.0030

Edit2:

pts1 and pts2 when triangulated using the essential matrix E, R and t returned [R, t] = cv.recoverPose(E, p1, p2,'CameraMatrix',K);

X1 =

   -0.0940    0.0478   -0.4984
   -0.0963    0.0497   -0.4987
    0.3033    0.1009   -0.5202
   -0.0065    0.0636   -0.5053
   -0.0737    0.0653   -0.5011

with

R =

   -0.9977   -0.0063    0.0670
    0.0084   -0.9995    0.0305
    0.0667    0.0310    0.9973

and

t =

    0.0239
    0.0158
    0.9996

When triangulated with the Matlab code, the chosen solution is E_all{2}

R_all{2}=

   -0.8559   -0.2677    0.4425
   -0.1505    0.9475    0.2821
   -0.4948    0.1748   -0.8512

and

t_all{2}=

   -0.1040
   -0.1355
    0.9853

X2 =

    0.1087   -0.0552    0.5762
    0.1129   -0.0578    0.5836
    0.4782    0.1582   -0.8198
    0.0028   -0.0264    0.2099
    0.0716   -0.0633    0.4862

When doing

X1./X2

ans =

   -0.8644   -0.8667   -0.8650
   -0.8524   -0.8603   -0.8546
    0.6343    0.6376    0.6346
   -2.3703   -2.4065   -2.4073
   -1.0288   -1.0320   -1.0305

There is an almost constant scale factor between triangulated 3D points. However, rotation matrices are different and there is no scale factor between translations.

 t./t_all{2}=
       -0.2295
       -0.1167
        1.0145

which makes the plotted trajectory wrong

1 Answers

Answering your numbered questions:

  1. Beware that Nister's 5 point algorithm has many implementations, but most of them don't work well. Personal experience and unpublished work by colleagues show that OpenCV does not have a good implementation. The open implementation in Bundler and other working SfM pipelines work better in practice (but there is a lot of room for improvement).

  2. The 10 solutions are simply zeros of a certain polynomial equation. As far as the polynomial equation can describe the problem, these 10 solutions all zero the equation. The equation does not describe that these 10 points are real, or that the 3D points corresponding to the 5 point correspondences have to be the same for each solution, but only that there are some 3D points (for each solution) that project to the 5 points, without even considering if the 3D points are in front of the respective cameras. Moreover, there may well be two sets of 3D points and cameras that happen to generate the same images of 5 points, so you would have to weed them out with some other procedure (below).

  3. The choice of the right solution among the 10 complex solutions is usually done by many techniques:

    • Discard solutions that would lead to purely complex points or 3D points with negative depth (currently Bundler does not do this last check)
    • Discard solutions that are not physical for some other reason (you may have to do some of that yourself for your application)
    • The more usual procedure: For each remaining solution, check which one is more consistent with additional correspondences. In a real system you don't know which additional correspondences are right and which are pure trash. So run RANSAC for each of the solutions and keep the one with the most inliers. This is computationally heavy so should be used as a last resort.

You can see how Bundler does this at file 5point.c line 668:

    generate_Ematrix_hypotheses(5, r_pts_inner, l_pts_inner, &num_hyp, E);

    for (i = 0; i < num_hyp; i++) {
        int best_inlier;
        double score = 0.0;

        double E2[9], tmp[9], F[9];
        memcpy(E2, E + 9 * i, 9 * sizeof(double));
        E2[0] = -E2[0];
        E2[1] = -E2[1];
        E2[3] = -E2[3];
        E2[4] = -E2[4];
        E2[8] = -E2[8];

        matrix_transpose_product(3, 3, 3, 3, K2_inv, E2, tmp);
        matrix_product(3, 3, 3, 3, tmp, K1_inv, F);

        inliers = evaluate_Ematrix(n, r_pts, l_pts, // r_pts_norm, l_pts_norm, 
                                   thresh_norm, F, // E + 9 * i, 
                                   &best_inlier, &score);

        if (inliers > max_inliers ||
            (inliers == max_inliers && score < min_score)) {
            best = 1;
            max_inliers = inliers;
            min_score = score;
            memcpy(E_best, E + 9 * i, sizeof(double) * 9);
            r_best = r_pts_norm[best_inlier];
            l_best = l_pts_norm[best_inlier];
        }

        inliers_hyp[i] = inliers;
    }
Related