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

Implemented direct conversion to/from euler angles in any sequence #24358

Merged
merged 15 commits into from
Dec 13, 2022

Conversation

evbernardes
Copy link
Contributor

@evbernardes evbernardes commented Dec 8, 2022

Brief description of what is fixed or changed

I recently published an article about a direct formula for the conversion between a quaternion variable to Euler angles in any sequence, which can be read here (Open Access).

Compared to either having 12 separate formulas (or 24, taking into account both intrinsic and extrinsic rotations) or using the well known quaternion-to-matrix-to-euler method (used for SciPy, for example), this has 3 main advantages:

  1. Numerically, it is up to 30 times faster than the previous quaternion-to-matrix-to-euler method (used for SciPy, for example).
  2. It is a lot simpler to implement, debug and maintain than both methods.
  3. It provides the simple formulas that, I imagine, can be used for theorerical work, which is why I think it could be useful in Sympy

Other comments

The axis sequence is defined by a string of length 3, which can represent either extrinsic (when the string is all uppercase) or intrinsic rotations (when all lowercase). This choice was for consistency with the Rotation class in SciPy, since it could be beneficial for the whole Python scientific community to have this consistency.

I also implemented a test that compares transforms a quaternion into euler angles then back to a quaternion, to show it gives the same result. This is done for all possible 24 seq cases.

Since the publication, this method has been merged into the main branch of SciPy (see pull request here).

Release Notes

  • algebras
    • Added direct conversion between Euler angles and quaternions.

@evbernardes evbernardes changed the title Implemented methods for direct conversion to and from euler angles Implemented direct conversion to/from euler angles in any sequence Dec 8, 2022
Comment on lines 48 to 52
axes = set(['x', 'y', 'z'])
if not (i in axes and j in axes and k in axes):
raise ValueError("Expected axes from `seq` to be from "
"['x', 'y', 'z'] or ['X', 'Y', 'Z'], "
"got {}".format(seq))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
axes = set(['x', 'y', 'z'])
if not (i in axes and j in axes and k in axes):
raise ValueError("Expected axes from `seq` to be from "
"['x', 'y', 'z'] or ['X', 'Y', 'Z'], "
"got {}".format(seq))
bad = set(seq) - set('xyzXYZ')
if bad:
raise ValueError("Expected axes from `seq` to be from "
"['x', 'y', 'z'] or ['X', 'Y', 'Z'], "
"got {}".format(''.join(bad)))

Comment on lines 19 to 24
if axis == 'x':
return [1, 0, 0]
if axis == 'y':
return [0, 1, 0]
if axis == 'z':
return [0, 0, 1]
Copy link
Member

@smichr smichr Dec 9, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if axis == 'x':
return [1, 0, 0]
if axis == 'y':
return [0, 1, 0]
if axis == 'z':
return [0, 0, 1]
return [1 if i == axis else 0 for i in 'xyz']

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice one-liner! Wouldn't this be harder for others to read?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The simple form seems pretty clear to me.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The simple form seems pretty clear to me.

I changed the code and put the suggested one-liner instead of calling the functions. This I can avoid adding two more function definitions to the code :)

Comment on lines 27 to 29
if axis == 'x': return 1
if axis == 'y': return 2
if axis == 'z': return 3
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if axis == 'x': return 1
if axis == 'y': return 2
if axis == 'z': return 3
return 'xyz'.index(axis) + 1

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as for _elementary_axis! Since they are oneliners, I wonder if it isn't better to just put them on their calls instead of defining these functions.

I think I'm gonna do that instead, to lower the diff.

@sympy sympy deleted a comment from sympy-bot Dec 9, 2022
@smichr
Copy link
Member

smichr commented Dec 9, 2022

Looks interesting (though I do not work in this area).

I made a few minor suggestions and added a "release notes" comment which should satisfy the bot. Check the test results to see how some tests are expected to be formatted. You should be able to copy/paste the desired output into the SymPy file.

@evbernardes
Copy link
Contributor Author

evbernardes commented Dec 9, 2022

Looks interesting (though I do not work in this area).

Thanks! :)

I made a few minor suggestions and added a "release notes" comment which should satisfy the bot. Check the test results to see how some tests are expected to be formatted. You should be able to copy/paste the desired output into the SymPy file.

I integrated the suggestions and added a line to .mailmap, and fixed the doctest!

About the "release notes", where can I see it? Apparently the bot deleted the comment!

@sympy-bot
Copy link

sympy-bot commented Dec 9, 2022

Hi, I am the SymPy bot (version not found!). I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

  • algebras

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.12.

Click here to see the pull request description that was parsed.
#### Brief description of what is fixed or changed
I recently published an article about a direct formula for the conversion between a quaternion variable to Euler angles in any sequence, which can be read [here](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0276302) (Open Access).

Compared to either having 12 separate formulas (or 24, taking into account both intrinsic and extrinsic rotations) or using the well known quaternion-to-matrix-to-euler method (used for SciPy, for example), this has 3 main advantages: 

1. Numerically, it is up to 30 times faster than the previous quaternion-to-matrix-to-euler method (used for SciPy, for example).
2. It is a lot simpler to implement, debug and maintain than both methods.
3. It provides the simple formulas that, I imagine, can be used for theorerical work, which is why I think it could be useful  in Sympy

#### Other comments

The axis sequence is defined by a string of length 3, which can represent either extrinsic (when the string is all uppercase) or intrinsic rotations (when all lowercase). This choice was for consistency with the `Rotation` class in SciPy, since it could be beneficial for the whole Python scientific community to have this consistency.

I also implemented a test that compares transforms a quaternion into euler angles then back to a quaternion, to show it gives the same result. This is done for all possible 24 seq cases.

Since the publication, this method has been merged into the main branch of SciPy [(see pull request here)](https://github.com/scipy/scipy/pull/17392).

#### Release Notes

<!-- BEGIN RELEASE NOTES -->
* algebras
  * Added direct conversion between Euler angles and quaternions.
<!-- END RELEASE NOTES -->

@github-actions
Copy link

github-actions bot commented Dec 9, 2022

Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with + are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with - are speedups.

Significantly changed benchmark results (PR vs master)

Significantly changed benchmark results (master vs previous release)

       before           after         ratio
     [41d90958]       [bd7bbeaa]
     <sympy-1.11.1^0>                 
-     1.17±0.03ms         774±20μs     0.66  solve.TimeSparseSystem.time_linear_eq_to_matrix(10)
-      3.65±0.3ms      1.45±0.05ms     0.40  solve.TimeSparseSystem.time_linear_eq_to_matrix(20)
-      6.91±0.1ms      2.18±0.06ms     0.32  solve.TimeSparseSystem.time_linear_eq_to_matrix(30)

Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).

@evbernardes
Copy link
Contributor Author

evbernardes commented Dec 9, 2022

Alright, all tests are passing now :) I'm just gonna make a last minor edit.

We can see that the benchmark from github actions shows the speed has changed, which doesn't surprised me because I implemented a test that symbolic shows that applying to_euler and then from_euler gives back the same result using simplify.

I think it's worth it to leave it like this.

EDIT.: Ok I just realized I misinterpreted the message from the actions bot. This is a speedup from the master.

sympy/algebras/quaternion.py Outdated Show resolved Hide resolved
sympy/algebras/quaternion.py Outdated Show resolved Hide resolved
sympy/algebras/quaternion.py Outdated Show resolved Hide resolved
sympy/algebras/quaternion.py Outdated Show resolved Hide resolved
sympy/algebras/quaternion.py Outdated Show resolved Hide resolved
sympy/algebras/quaternion.py Outdated Show resolved Hide resolved
@smichr
Copy link
Member

smichr commented Dec 9, 2022

ping @NikhilPappu

@evbernardes
Copy link
Contributor Author

@smichr thanks for reviewing the code! I'll be modifying it next Monday :)

evbernardes and others added 2 commits December 12, 2022 18:30
As suggested by @smichr

Co-authored-by: Christopher Smith <smichr@gmail.com>
As suggested by @smichr

Co-authored-by: Christopher Smith <smichr@gmail.com>
evbernardes and others added 3 commits December 12, 2022 18:31
As suggested by @smichr

Co-authored-by: Christopher Smith <smichr@gmail.com>
@evbernardes
Copy link
Contributor Author

@smichr and @sylee957, I integrated your last suggestions!

@@ -84,6 +112,158 @@ def d(self):
def real_field(self):
return self._real_field

@classmethod
def from_euler(cls, angles, seq):
"""Returns quaternion equivalent to Euler angles represented same in
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
"""Returns quaternion equivalent to Euler angles represented same in
"""Returns quaternion equivalent to Euler angles represented in

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oops, that's indeed a typo! I meant to say Returns quaternion equivalent to rotation represented by the Euler angles.

return trigsimp(qi * qj * qk)

def to_euler(self, seq):
"""Returns Euler angles representing same in the sequence given by
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same issue as above -- "representing same". I am not sure what is the same, i.e. it's unclear.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as last one, I meant "same rotation". Fixed both!

@smichr
Copy link
Member

smichr commented Dec 13, 2022

Thanks. This seems like a "sliced bread" commit -- a nice addition for those that use it.

@smichr smichr merged commit db1d91d into sympy:master Dec 13, 2022
@evbernardes
Copy link
Contributor Author

Thanks. This seems like a "sliced bread" commit -- a nice addition for those that use it.

Awesome! Thank you for merging it :)

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

Successfully merging this pull request may close these issues.

None yet

5 participants