[Multiview] Add lambda twist p3p solver #1500
Conversation
- for now the code is licensed under GPL 3.0
|
Thank you guys. I will check why we have a compilation error for CLANG and VSTUDIO ;-) and make a code review. |
- few style modifications (spaces etc) - Add back and complete docstrings from authors
|
Hi,
It is crucial that the 1 root cubic solver is used. This impacts both speed
and accuracy substantially compared to using any standard 3 root cubic
solver.
The difference is caused by the numeric instability of the common higher
multiplicity roots etc.
The same thing is true for the specialized eigenvalue decomp. It cannot be
replaced by a generic component without both speed and accuracy loss.
The second order root solver isn't important.
With regards to unnecessary allocation I did mention that std vector should
not be used, and it isn't in the original. Prealloced STD vectors shouldn't
be used either because they still break alot of compiler optimisations.
Cheers
//Mikael
…On Sun, 14 Apr 2019, 21:11 Romain Janvier, ***@***.***> wrote:
***@***.**** commented on this pull request.
------------------------------
In src/openMVG/multiview/solver_resection_p3p_nordberg.hpp
<#1500 (comment)>:
> + if(v < 0.0){
+ r1 = r2 = 0.5 * b;
+ return false;
+ }
+ double y = std::sqrt(v);
+ if(b < 0.0){
+ r1 = 0.5 * (-b +y);
+ r2 = 0.5 * (-b -y);
+ }else{
+ r1= 2.0 * c / (-b + y);
+ r2= 2.0 * c / (-b - y);
+ }
+ return true;
+};
+
+static double cubick(const double & b, const double & c, const double & d){
yes I reread the paper after 6 month, it seems to one of the advantage of
the method over others p3p solvers.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
<#1500 (comment)>, or mute
the thread
<https://github.com/notifications/unsubscribe-auth/AFMMJR08Ly497KxTR8jfBVgR9OVr7OD3ks5vg31GgaJpZM4cetzb>
.
|
|
The first own, the vector.
…On Sun, 14 Apr 2019, 22:42 Romain Janvier, ***@***.***> wrote:
***@***.**** commented on this pull request.
------------------------------
In src/openMVG/multiview/solver_resection_p3p_nordberg.cpp
<#1500 (comment)>:
> +
+ Mat3 A;
+ A << A00, A01, A02,
+ A01, A11, A12,
+ A02, A12, A22;
+
+ // get sorted eigenvalues and eigenvectors given that one should be zero...
+ Mat3 V;
+ Vec3 L;
+
+ eigwithknown0(A, V, L);
+
+ double v = std::sqrt(std::max(0.0, -L(1) / L(0)));
+
+ int valid = 0;
+ std::array<Vec3, 4> Ls;
@midjji <https://github.com/midjji> or maybe the use of the std::array
here?
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
<#1500 (review)>,
or mute the thread
<https://github.com/notifications/unsubscribe-auth/AFMMJeOtz3hDWjEylK0CagpgApx_oJmmks5vg5KvgaJpZM4cetzb>
.
|
|
|
Took a look at the 3 root cubic solver linked,
While cos formulation isn't the best numerically in general.
I think you also have specific numerical issues for q~0 in the std::arcos,
this is especially bad since a cubic solver should never throw a NAN. the
main problem is that the accuracy of p^3 isn't guaranteed to be caught by
the numeric check on p and the threshold probably needs to be increased.
The input must also have been explicitly normalized, which is unwise, it
requires the user to properly account for the special cases which is ...
unlikely. Performing that division is also not required, and doing it in
advance actually reduces numerical accuracy compared to keeping it in the
formulas.
Throw a few billion random polynomials at it and id wager you'll get a few
nans and quite a few bad solutions near roots of higher multiplicity. Might
want to explicitly test bad cases like (x-0.99)(x-1 +-1e-8) too.
There are specializations many cases you'll want to add for robustness, but
I only know of a few,
but most importantly, the roots would gain several bits of accuracy from
refinement.
Not sure if any of this is important, but at least it might be good to be
aware.
…On Sun, 14 Apr 2019 at 22:52, Romain Janvier ***@***.***> wrote:
The first own, the vector.
@midjji <https://github.com/midjji>
Ok, I kept the interface we have for Ke's solver. I will stick to your
implementation then.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
<#1500 (comment)>, or mute
the thread
<https://github.com/notifications/unsubscribe-auth/AFMMJTnFwSpq5mM1Lssbj6Z17A1N3sAvks5vg5UBgaJpZM4cetzb>
.
|
|
@midjji Thank you for your feedback. I agree that each solver should use the best numerical precision and sometimes it forces to uses a specific root solver that could be optimized for speed and accuracy. We will keep your specific root solver. I was trying to see if we could factorize some code. We are glad to have your feedback. Regarding vector vs array. In order to keep some genericity due to the fact that some solver can return multiple solutions, we need vector at the end. So std::array could be used inside the solver and then only a single solution will be returned in your case. |
|
@pmoulon it returns up to 4 solutions like any other p3p solver. so the question is : should we keep the |
|
I think you should add #include in solver_resection_p3p_nordberg.cpp to compile in visual studio |
Thanks @RamadanAhmed have you any idea of the specific missing include? |
|
I agree that the Interface compatability forces the use of the vector, it's
pretty DN... convenient too. The performance increase will be minor as long
as the vector is prealloced to 4. Romain is right that it will return
multiple valid solutions, generally fewer total but more distinct ones than
the other solvers too.
Is the use of rotation matrixes important btw? I only used that for fair
speed comparisons with ke and kneip, the unit quaternion version is
slightly better.
…On Tue, 16 Apr 2019, 11:22 Romain Janvier, ***@***.***> wrote:
@pmoulon <https://github.com/pmoulon> it returns up to 4 solution like
any other p3p solver.
so the question is : should we keep the std::vector (for sake of
consistency with the other p3p interfaces in openMVG) or use arrays (to be
truly fair with the original code and in order to maximize speed) in
computePosesNordberg?
Even if i think the speed is not an issue for us, since openMVG does not
target (primarly) real time applications, I think it's not a big dealhere
for openMVG if we stick with the original code, since the call of
computePosesNordberg is wrapped inside another function anyway (
P3PSolver_Nordberg::solve) and it hides the fact we use preallocated
arrays vs. dynamic vectors.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
<#1500 (comment)>, or mute
the thread
<https://github.com/notifications/unsubscribe-auth/AFMMJR3RPJa6l1ArHBacPuxcXJoRB_DGks5vhZZSgaJpZM4cetzb>
.
|
"#include array" sorry I forget to mention :) |
|
Ok if I understand well we keep the vector in the `computePosesNordberg``? |
|
@rjanvier I agree with your comment #1500 (comment)
Congrats team we are very close to having another minimal solver in Develop to be tested! 🎉 |
|
I have also tested it and it works much better than P3P_Ke, I think it should be the main solver for p3p |
|
Thank you @RamadanAhmed for testing the WIP feature. Appreciated!
We just need to test its behavior for pinhole and 360 images. @rjanvier Seems like a minor concern is still there to the variable fx here https://app.codacy.com/app/pmoulon/openMVG/pullRequest?prid=3389059 (I can fix it) |
|
@pmoulon I think allocating fx/fpx into the loop could impact the perfomance... AFAIK it really depends on the compiler so I haven't fixed it on purpose. |
|
It could, but in practice dont worry about it for recent compilers in this
specific case.
As you might have noticed, I chased every last cycle and every last bit of
accuracy.
And I wanted to do so in a way which is cross architecture friendly, so no
explicit asm.
The latter means a bunch of strange choices in places and alot of godbolt.
…On Wed, 17 Apr 2019 at 19:15, Romain Janvier ***@***.***> wrote:
@pmoulon <https://github.com/pmoulon> I think allocating fx/fpx into the
loop could impact the perfomance... AFAIK it really depends on the compiler
so I haven't fixed it on purpose.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
<#1500 (comment)>, or mute
the thread
<https://github.com/notifications/unsubscribe-auth/AFMMJfI01Qwd_WYa54d2_lmoaAK7NwUXks5vh1aYgaJpZM4cetzb>
.
|
|
360 images? hmm might not work for that, whats the projection model?
Let me know if it does or not, if it doesnt that be great, since then I can
fix that.
On Thu, 18 Apr 2019 at 14:05, Mikael Persson <mikael.p.persson@gmail.com>
wrote:
…
It could, but in practice dont worry about it for recent compilers in this
specific case.
As you might have noticed, I chased every last cycle and every last bit of
accuracy.
And I wanted to do so in a way which is cross architecture friendly, so no
explicit asm.
The latter means a bunch of strange choices in places and alot of godbolt.
On Wed, 17 Apr 2019 at 19:15, Romain Janvier ***@***.***>
wrote:
> @pmoulon <https://github.com/pmoulon> I think allocating fx/fpx into the
> loop could impact the perfomance... AFAIK it really depends on the compiler
> so I haven't fixed it on purpose.
>
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub
> <#1500 (comment)>,
> or mute the thread
> <https://github.com/notifications/unsubscribe-auth/AFMMJfI01Qwd_WYa54d2_lmoaAK7NwUXks5vh1aYgaJpZM4cetzb>
> .
>
|
|
@midjji it's an equirectangular projection model. I did a test yesterday with 360 images, lambda twist works at least as well as Ke's solver (regarding behavior). We use bearing vectors as input (not 2 coordinates in sensor's frame) and it seems you use a bearing vector formulation too in the solver, so it's general enough but maybe I'm wrong and I'm missing some important points here, please let me know. thanks a lot! |
|
Ok my bad, it seems that the spherical model violate the assumption that all points must have the same signed distance w.r.t to the camera center. Am I right? |
|
please forget my last comment, my synthetic test was wrong. |
|
I'm unclear on if the 360 case worked or not?
…On Fri, 19 Apr 2019, 18:49 Romain Janvier, ***@***.***> wrote:
***@***.**** commented on this pull request.
------------------------------
In src/openMVG/multiview/solver_resection_p3p_nordberg.cpp
<#1500 (comment)>:
> +bool computePosesNordberg(
+ const Mat &bearing_vectors,
+ const Mat &X,
+ std::vector<std::tuple<Mat3, Vec3>> &rotation_translation_solutions)
+{
+ // Extraction of 3D points vectors
+ Vec3 P1 = X.col(0);
+ Vec3 P2 = X.col(1);
+ Vec3 P3 = X.col(2);
+
+ // Extraction of feature vectors
+ Vec3 f1 = bearing_vectors.col(0);
+ Vec3 f2 = bearing_vectors.col(1);
+ Vec3 f3 = bearing_vectors.col(2);
+
+ f1.normalize();
thanks
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
<#1500 (comment)>, or mute
the thread
<https://github.com/notifications/unsubscribe-auth/ABJQYJMKZF2MZDTAL6YGIHDPRHZZHANCNFSM4HD23TNQ>
.
|
|
I left in explicit L2 normalization, because its one of those things which
while required, might accidentally work in many cases, causing annoyring to
debug errors.
…On Mon, 22 Apr 2019 at 16:43, Romain Janvier ***@***.***> wrote:
***@***.**** commented on this pull request.
------------------------------
In src/openMVG/multiview/solver_resection_p3p_nordberg.cpp
<#1500 (comment)>:
> +bool computePosesNordberg(
+ const Mat &bearing_vectors,
+ const Mat &X,
+ std::vector<std::tuple<Mat3, Vec3>> &rotation_translation_solutions)
+{
+ // Extraction of 3D points vectors
+ Vec3 P1 = X.col(0);
+ Vec3 P2 = X.col(1);
+ Vec3 P3 = X.col(2);
+
+ // Extraction of feature vectors
+ Vec3 f1 = bearing_vectors.col(0);
+ Vec3 f2 = bearing_vectors.col(1);
+ Vec3 f3 = bearing_vectors.col(2);
+
+ f1.normalize();
@pmoulon <https://github.com/pmoulon> in fact they are not, at least in
pinhole camera model the operator() return return Kinv_ *
points.colwise().homogeneous(); so the point is expressed in homogenous
coordinates and it's not an unit vector (will be v(somescalar,
somescalar,1)).
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
<#1500 (comment)>, or mute
the thread
<https://github.com/notifications/unsubscribe-auth/ABJQYJI55CVWNBKUBTDAMSLPRXFJJANCNFSM4HD23TNQ>
.
|
|
@midjji 360 camera works on my experiments (both real cases and synthetics). I'm not 100% sure though because I do not have a strong mathematical background (I came from an "humanities" background) and I'm not able to extract by myself the "proof" from your code/paper. |
- few style modifications (spaces etc) - Add back and complete docstrings from author - fix some code formating issues thanks @pmoulon - move functions into the cpp
- add resection::SolverType::P3P_NORDBERG_ECCV18 to the SfM_Localizer - set resection::SolverType::P3P_NORDBERG_ECCV18 as the default p3p solver for sequential SfM v1, v2 and the image localization client.
|
@RamadanAhmed This feature is ready to test, feel free to try it in the develop branch and provide any feedback. |
cc Mikael Persson (@midjji)
This PR contains the code the of the P3P method exposed in Persson M., Nordberg K., Lambda Twist: An Accurate Fast Robust Perspective Three Point (P3P) Solver, ECCV 2018.
I had just adapted the code present in the original repository (https://github.com/midjji/lambdatwist-p3p) to Eigen/openMVG. I kept some comments and maybe some duplicated methods w.r.t what we already have in openMVG. I think further refinements, if needed, are out of my league and maybe @midjji and @pmoulon you could now further discuss about how to best integrate lambda twist in openMVG.
According to discussion we had with Mikael, the code is now MPL-2 licensed (instead of the original GPL-3.0 scheme) and the method is named in honor of Klas Nordberg.
The implementation is tested and functional. I made some extensive real life tests, using it for my main solver in several dataset (one was about 7790 HR images of Château de Chambord ;)). It's very stable and plays nicely with AC-RANSAC. I think it should be default p3p solver in openMVG.
Thank you so much to Klas Nordberg and Mikael Persson for making this code public under a permissive license.