Openmvg: Potential issue in lambda twist p3p inclusion to OpenMVG

[mpizenberg]

@rjanvier

Hi, I'm currently porting the lambda twist p3p algorithm to Rust and so first, I'm basically rewriting solver_resection_p3p_nordberg.cpp into a compilable Rust version. While reading the gauss_newton_refineL function I think I spotted a mistake that makes this function potentially wrong.

At lines 80-82, l1, l2, l3 are recomputed using the same unchanged L vector. This means that the new residual is always the same that the old one and the condition at line 86 is therefore always false, leading to L always being updated by the new lambda L1, even in cases where the residual increases.

Let me know if I'm mistaken, I'm not very familiar with the algorithm but this stroke me as something weird.

[pmoulon]

HI @mpizenberg,

Thx you for your feedback.

Are you considering the last version that exist in the develop branch (see https://github.com/openMVG/openMVG/commit/c37c0824e101382290290137a90c3edcd3a7ee5e)

Could you make a wrapper of C++ code in Rust, so this way you could leverage existing code?

[mpizenberg]

Ah, yes it is fixed in develop indeed! sorry for the noise.

[pmoulon]

No worries. It is good to have feedback.

[mpizenberg]

I have another question related to the code if I may.

In function gauss_newton_refineL documentation, the author defines the b12, b13, b23 parameters as being the cosine of angle between corresponding bearing vectors.

In the main function computePosesNordberg the previous function is called as follows:

gauss_newton_refineL(Ls[i], a12, a13, a23, b12, b13, b23);

Where b12, b13, b23 are previously defined as follows:

double b12 = -2.0 * (f1.dot(f2));
double b13 = -2.0 * (f1.dot(f3));
double b23 = -2.0 * (f2.dot(f3));

In those lines, f1, f2, f3 are the normalized bearing vectors, so these are not exactly the cosine values. I am wondering if the refine function should be called with c12, c23, c31 instead, that are defined as:

double c31 = -0.5 * b13;
double c23 = -0.5 * b23;
double c12 = -0.5 * b12;

So is the issue with the documentation comments, the original implementation, or the rewrite? Knowing that the cosines are indeed called b12, b23, b13 in the paper, it may have been the cause of the confusion. Or maybe it is just a documentation comment issue.

[mpizenberg]

Hum, actually I should have just check the reference implementation. The original author implementation is the same that the one here. And I just tested locally with c12, ... and it gives wrong results. So it is just a documentation issue.

[mpizenberg]

This notation difference between the paper and the code for b12, ... picked my curiosity so I have tried to look for inconsistencies between the paper and the code. I found two of these, the second one potentially with consequences.

1. The coefficient c in the paper in equation (15) does not correspond to the code version. I believe it is just a parenthesis typo. The paper version is:

(a13 + a13 - a12) * w0^2 - 2 * a13 * b12 * w0 // paper notation

The code version is (cf code)

a13 + (a13 - a12) * w0^2 + a13 * b12 * w0 // code notation

The issue here is that the paper version is doing (a13 + a13 instead of a13 + (a13.
The difference with the sign and the multiplication by 2 in the second term is just normal due to the notation difference between paper and code.

1. Equation (16) in the paper is followed by "which is solved by ..." in a way that is wrong if we simply develop equation (16). The paper states that the denominator of lambda_2k squared is:

tau * (tau + b23) + 1 // paper notation

while it should be

tau * (tau - 2 * b23) + 1 // paper notation

I believe the author of the paper made a mistake in taking the code notation instead of the paper notation. What is a bit more suspect is that they continue justifying that this denominator is positive since tau > 0 and |b23| <= 1.0. However we have -2 * b23 (paper notation) and not b23, so the proof of the denominator being positive does not hold without further conditions on tau.

I thought I should raise this to your attention since lambda twist is being integrated to OpenMVG. Please let me know if I'm mistaken somewhere. If not, this means that we should also check for the sign of the denominator before continuing to compute and validate l1, l2, l3.

[rjanvier]

Thanks for pointing this to our attention, very appreciated. I feel you should report this upstream to M. Persson @midjji since it’s more about difference between the paper and the reference implementation than about integration in openMVG. Thanks a lot!

[midjji]

Hi,
good catch, this is indeed an error in the paper. It is an error in two ways, first because the reasoning is incorrect and second because it does not hold numerically in general.

The fix, which I have pushed to the reference implementation, is to add a check to discard those solutions, see below. This should also be added to the openmvg version.

if(tau1>0) {
    T tau=tau1;
    T d=a23/(tau*(b23 + tau) + T(1.0));
    if(d>0){ // this line
      T l2=std::sqrt(d);

      T l3=tau*l2;

      T l1=w0*l2 +w1*l3;
      if(l1>=0){
          Ls[valid]={l1,l2,l3};
          ++valid;
      }
    } // this line
}
if(tau2>0){
    T tau=tau2;
    T d=a23/(tau*(b23 + tau) + T(1.0));
    if(d>0){ // this line
      T l2=std::sqrt(d);

      T l3=tau*l2;

      T l1=w0*l2 +w1*l3;
      if(l1>=0){
          Ls[valid]={l1,l2,l3};
          ++valid;
      }
    } // this line
}  

[pmoulon]

Nit: Updated last post to format the code in a more readable way

[midjji]

Just to be clear, I think the statement is correct, but the proof is wrong but its also wrong numerically and for degenerate cases and only holds for limited pinhole cameras. Arxiv version is on its way.