Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Code error in Geometry? #13

Closed
dondonddon opened this issue Feb 13, 2022 · 8 comments
Closed

Code error in Geometry? #13

dondonddon opened this issue Feb 13, 2022 · 8 comments

Comments

@dondonddon
Copy link

Is line 53 in Geometry.cpp supposed to say: theta = 1.0e-5 ??

I am trying to match up your code with the algorithm description in: https://www.cl.cam.ac.uk/techreports/UCAM-CL-TR-696.html

Also, you have a class called AngularVelocity, but I think the values are supposed to be gyro rates * integration time ? Otherwise, how do you account for the delta-T in order to integrate the angular rates?

@tomstewart89
Copy link
Owner

tomstewart89 commented Feb 13, 2022

Hey @dondonddon, thanks for your interest in my library. Let me take a run at your questions:

Is line 53 in Geometry.cpp supposed to say: theta = 1.0e-5 ??

That if-statement is meant to avoid divide-by-zero errors that would occur on line 56 when AngularVelocitys that are close to zero are passed to exp. In those cases, with theta = 1.0 on line 53, so3_skew will be close to zero and the terms involving theta on line 57 will evaluate to ~zero. With theta = 1e-5 then so3_skew won't be close to zero but the same terms on line 57 will evaluate to ~zero anyway since sin(1e-5)~= 0 and 1 - cos(1e-5) ~=0.

So as far as I can tell both 1.0 and 1e-5 will both produce more or less the same result. Are you having some trouble with that function? Is there something I'm missing there?

Also, you have a class called AngularVelocity, but I think the values are supposed to be gyro rates * integration time ? Otherwise, how do you account for the delta-T in order to integrate the angular rates?

Angular velocity represents the amount by which a rigid body will rotate per unit time (just as linear velocity represents the distance a body will travel per unit time). If you want to calculate how much a body will rotate in some integration time then yeah you'll need to scale accordingly.

For what it's worth, if you have a rigid body starting at some initial rotation R_t which rotates with a constant angular velocity w over some time interval dt then the code to integrate this velocity looks like this:

Rotation R_t;
AngularVelocity w;
double dt;

Rotation R_t_plus_dt = exp(w * dt) * R_t; // rotation of some rigid body at time t + dt

Hope that helps!

@dondonddon
Copy link
Author

dondonddon commented Feb 14, 2022 via email

@tomstewart89
Copy link
Owner

Hi Don,

Looking at the case for small theta, I think the goal was really to use a small angle approximation

That's a sensible guess, but to be honest I just wanted to avoid the potential divide-by-zero error on line 56 by setting theta to an arbitrary (non-zero) value. Since so3_skew would be a 3x3 matrix whose elements are almost zero in those cases anyway, the logic was that multiplying it by sin / cos of an arbitrary constant wouldn't really affect the output of the exp function.

On reflection though, I can see how that can give quite inaccurate results when integrating rates from a gyro. How would you feel about the following implementation? Here we explicitly return an identity rotation when the magnitude of theta is below 1e-5. Depending the sensitivity and sampling rate of your gyro that 1e-5 might cause some of your readings to be ignored and if that's the case we can change the underlying element type from float to double and then decrease that to say 1e-7.

Rotation exp(const AngularVelocity& w)
{
    auto theta = Norm((BLA::Matrix<3>)w);

    if (fabs(theta) < 1e-5)
    {
        return BLA::Identity<3>();
    }
    else
    {
        auto so3_skew = skew(w / theta);
        return BLA::Identity<3>() + so3_skew * sin(theta) + so3_skew * so3_skew * (1 - cos(theta));
    }
}

Incidentally, regarding this point:

I still ponder why you used the label so3_skew

The so there refers to the special orthogonal group; a lie group which represents the set of all 3D rotations (hence the 3). If you want to read up on that in the context of robotics have a look at chapter 3 of this book, it's a good read!

@dondonddon
Copy link
Author

dondonddon commented Feb 18, 2022 via email

@tomstewart89
Copy link
Owner

Hey @dondonddon,

Perhaps it is better to rename those methods: firstAngle(), secondAngle(), thirdAngle() – and leave the remaining code as is?

That sounds reasonable, I went ahead and made those changes (I went with just first second and third though). I also pushed the changes relating to handling small angular velocities in exp. In the end I just dropped the 1e-5 epsilon and compared theta to 0. Since we're only looking to avoid divide-by-zeros here I think the floating point equality comparison is ok.

You can check out both of those changes here: https://github.com/tomstewart89/Geometry/releases/tag/2.2

And with that I might close this issue. Thanks again for your feedback, I'm glad you're getting some use out of this library (give it a star if you haven't already!). If you have any other suggestions feel free to open up another issue.

@dondonddon
Copy link
Author

dondonddon commented Feb 23, 2022 via email

@dondonddon
Copy link
Author

dondonddon commented Feb 23, 2022 via email

@dondonddon
Copy link
Author

dondonddon commented Feb 23, 2022 via email

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants