[SofaKernel] FIX bug in toEulerVector

[raffaellatrivisonne]

---

This PR:

• [x] builds with SUCCESS for all platforms on the CI.
• [x] does not generate new warnings.
• [x] does not generate new unit test failures.
• [x] does not generate new scene test failures.
• [x] does not break API compatibility.
• [x] is more than 1 week old (or has fast-merge label).

Reviewers will merge only if all these checks are true.

[sofabot]

Thank you for your pull request! Someone will soon check it and start the builds.

[guparan]

Thanks for your PR @raffaellatrivisonne :-) Could you add a small description explaining the error and your fix please? [ci-build]

[damienmarchal]

[ci-build]

[guparan]

@raffaellatrivisonne could you fix the build failure please?

[damienmarchal]

[ci-build]

[raffaellatrivisonne]

The function asin is defined in [-1,1]. The fix prevents NAN when the argument is slightly >1 due to numerical errors (1,000000000000001). Hope it's clear enough.

[guparan]

From SOFA-dev meeting: Hi @raffaellatrivisonne, is there a chance that the argument y is rounded to 1 or -1 due to anything else than a very small numerical error? If yes, we should add an epsilon to ensure the error is acceptable.

[untereiner]

I think I get your point but it adds an another parameter that changes nothing in practice. Here the point is to cut everything outside the bounds. With your suggestions, who will test the sensibility of this epsilon parameter ? It is not acceptable to have a number outside these bounds since asin is not defined at all

[guparan]

asin might be undefined, it is still better to have it crashing instead of working were it should not - again, if y contains a "big" error (greater than 10e-15 upon @raffaellatrivisonne comment). The question is: could y contain such a "big" error? If no, let's just assume that and round this way. If yes or don't know, let's add a tiny security with something like

const double epsilon = 1e-15;
Real y = Real(2.)*(q[3]*q[1] - q[2]*q[0]);
if( std::abs( double(y) ) - 1.0 > epsilon )
{
    msg_error("Quat") << "Unexpectedly out of bounds argument for asin: " << y << msgendl;
}

or

const double epsilon = 1e-15;
Real y = Real(2.)*(q[3]*q[1] - q[2]*q[0]);
if( std::abs( double(y) ) - 1.0 > epsilon )
{
    Real force_round = std::max( Real(-1.0), std::min(Real(1.0), y) );
    msg_warning("Quat") << "Unexpectedly out of bounds argument for asin: " << y
                        << "Force rounding to " << force_round << msgendl;
    y = force_round;
}

[untereiner]

Why 1e-15 ? why not 1e-11 ? I guess Raphaella's number was just an example. I completely disagree with this idea. Introducing a new random chosen very locally defined epsilon is really a bad idea. But if it a command from the consortium guy's to close this pr I have to bow

[guparan]

There is no command and I am far to pretend knowing enough how things are done in SOFA to claim my proposal is what has to be done. It is just a proposal discussed during last SOFA-dev meeting.

[untereiner]

Well, the message sounded like "you have to add an epsilon value" 🙄 . I let @raffaellatrivisonne do her job here for now on... good luck

[raffaellatrivisonne]

Hi everybody, sorry for being late with my answer. @guparan, as soon as the argument of asin is greater than 1, you have an invalid operation. It doesn't matter how "big" is your error, it can be huge or infinitesimal. In this case, it is due to numerical errors, that's why I said 10e-15. If it is the proper way, I can set an epsilon, but then how ? As @untereiner says, why 1e-15 instead of 1e-11 ? How do you set a parameter on numerical errors ?

[fjourdes]

Just passing by, I do not want to raise a flame war but actuall asin method has some documentation, espcially when it comes to domain error: http://www.cplusplus.com/reference/cmath/asin/ Just quoting;

If a domain error occurs:

• And math_errhandling has MATH_ERRNO set: the global variable errno is set to EDOM.
• And math_errhandling has MATH_ERREXCEPT set: FE_INVALID is raised.

Then looking at http://www.cplusplus.com/reference/cmath/math_errhandling/ The default behavior for math_errhandling is MATH_ERRNO, so as the doc suggest you can just check for the errno (thread-local) global variable value, and if it is set to EDOM after asin is called, then you can throw whatever error message you want.

This is probably not relevant here, since I presume the checks are there because in theory when a quaternion is normalized the value of Real(2.)*(q[3]*q[1] - q[2]*q[0]) should always belong to the range [-1;1]. The only reason it might not be is for some numerical drifting issues (?) Provided this assumption is correct you are indeed totally allowed to clamp the values there.

[guparan]

Thanks for the explanations @raffaellatrivisonne and @fjourdes, let's merge then ! :-)

[damienmarchal]

[ci-build]

[damienmarchal]

Do someone knows why do tests are failing ? (In QuaterTest)

[raffaellatrivisonne]

I'm taking a look.

[raffaellatrivisonne]

I'm trying to debug to see where it fails, but it will take a little bit more than expected. Meanwhile I just want to highlight that coming back to the un-fixed version (the one without my commit) the test doesn't fail YET the toEulerVector is doing an invalid operation (the bug that my commit is supposed to fix).

[guparan]

So if I understood well we have 2 problems here:

• q1 and p1 should not be nan
• EXPECT_EQ(p0,p1); should not pass with nan values for p1

[raffaellatrivisonne]

The inital problem is the vector (that in my screenshot I call "q0 euler") which has its second component which is NaN. It comes from the unfixed toEulerVector "q0 euler"= q0.Quater::toEulerVector() So, it doesn't surprises me that q1 and p1 are NaN, As they derive from: Quater q1 = Quater::createQuaterFromEuler(q0.Quater::toEulerVector());

    sofa::defaulttype::Vec<3,double> p1 = q1.Quater<double>::rotate(p);

But definetively EXPECT_EQ(p0,p1) should not pass with nan values for p1

[fjourdes]

maybe I am stating something already well known, but with c++11 there are some built in functions that can help to test floating point arithmetic. So with the current implementation adding something to the EulerAngle test like

for(std::size_t i=0; i<q0.size(); ++i) // same goes for q1
{
    ASSERT_FALSE(std::is_nan(q0[i]));
}

[raffaellatrivisonne]

Hi everybody, I'm trying to work on this PR, but I'm quite busy with my PhD in this moment and I don't think I will be able to finish it within a short delay. As suggested by Francois, I added the error message in Quater_test.cpp. Now is failing, as with the old-code (without my commit) NaN values may appear. In Quater.inl (function toEulerVector) I went back to the old code commenting the modifications I made with my commit. This way, if someone else takes the hands on this PR, he will better know what to do.

[guparan]

Hey, I just cleaned this oooooold PR. Should be working fine :tada: