Skip to content

Conversation

@nbelakovski
Copy link
Contributor

Reference issue

N/A

What does this implement/fix?

The implements 2 new methods in the Rotation class of the spatial/transform package, as_mrp and from_mrp, which allow the class to accept input and provide output in the form of Modified Rodrigues Parameters.

@tylerjereddy tylerjereddy added enhancement A new feature or improvement scipy.spatial labels Aug 5, 2020
Copy link
Contributor

@tylerjereddy tylerjereddy left a comment

Choose a reason for hiding this comment

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

Just added a few minor comments.

@nmayorov will probably have more detailed comments

I see the manuscript has been cited > 2000 times, which may help for justification of breadth of interest, but I'm not immediately familiar with this myself.

Notes
-----

.. versionadded:: 1.5.3
Copy link
Contributor

Choose a reason for hiding this comment

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

Just FYI -- this will probably have to be 1.6.0 if the feature is approved because we have already released the 1.5.x series and almost always restrict point releases to bug fixes rather than enhancements/new features.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Updated

quat[ind, 0] = 2 * cmrp[ind, 0] / (1 + ss)
quat[ind, 1] = 2 * cmrp[ind, 1] / (1 + ss)
quat[ind, 2] = 2 * cmrp[ind, 2] / (1 + ss)
quat[ind, 3] = (1 - ss) / (1 + ss)
Copy link
Contributor

Choose a reason for hiding this comment

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

minor: not sure if the compiler will optimize the repeated (1 + ss) calls or if assigning the value to i.e., a variable once might help

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I'm not sure either, I went ahead and changed ss to 1+ss

Notes
-----

.. versionadded:: 1.5.3
Copy link
Contributor

Choose a reason for hiding this comment

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

1.6.0

@nbelakovski nbelakovski marked this pull request as ready for review August 5, 2020 13:32
@nbelakovski
Copy link
Contributor Author

I see the manuscript has been cited > 2000 times, which may help for justification of breadth of interest, but I'm not immediately familiar with this myself.

I'm not sure myself what the breadth of interest is. I'd never heard of them until I started taking a course on Coursera for spacecraft attitude control. Perhaps it would be better to split this work out into it's own little Python library and then merge it into scipy if it seems popular enough?

@nmayorov
Copy link
Contributor

nmayorov commented Aug 5, 2020

MRP fits well into Rotation, it's a small incremental addition which I support. I will review on the weekend.

Copy link
Contributor

@nmayorov nmayorov left a comment

Choose a reason for hiding this comment

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

Initial review.

def from_mrp(cls, mrp):
"""Initialize from Modified Rodrigues Parameters (MRPs).

MRPs are a 3 element attitude description commonly used in the aerospace
Copy link
Contributor

Choose a reason for hiding this comment

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

I think a more substantive definition is required here.

The easiest way it to define it through axis and angle of rotation, perhaps pointing out similarity to rotation vector but with different scaling.

Copy link
Contributor

Choose a reason for hiding this comment

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

description -> representation?

Do you think it's valuable to mention its usage in the aerospace industry? I'd say all common representations are widely used in the aerospace industry. I would probably suggest removing.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I went ahead and used your provided definition. I think MRPs are not widely used outside of aerospace, so I guess I'd say the target audience for that factoid is people outside of the aero industry, but I don't think it hurts one way or another if the factoid is in the docs or not, we can leave it out.


Initialize a single rotation:

>>> r = R.from_mrp([0, 0, 0.41421356])
Copy link
Contributor

Choose a reason for hiding this comment

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

I suggest using unity MRP that corresponds to rotation of 180 degrees.

Perhaps use as_euler in degrees to show this correspondence.

>>> r.as_rotvec().shape
(2, 3)

It is also possible to have a stack of a single rotaton:
Copy link
Contributor

Choose a reason for hiding this comment

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

I think this example is not necessary.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yes, but all the other definitions have this same example.

Copy link
Contributor

Choose a reason for hiding this comment

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

The examples in Rotation are not very good, no point following them in my opinion.

cdef double ss, qnorm

for ind in range(num_rotations):
ss = 1 + _dot3(cmrp[ind, :], cmrp[ind, :])
Copy link
Contributor

Choose a reason for hiding this comment

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

Perhaps name it at least denom or denominator.

def as_mrp(self):
"""Represent as Modified Rodrigues Parameters (MRPs).

3D rotations can be represented using Modified Rodrigues Parameters, which
Copy link
Contributor

Choose a reason for hiding this comment

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

Again work on definition and make sure they are the same in from_mrp and as_mrp.

>>> r.as_mrp().shape
(3,)

Represent a stack with a single rotation:
Copy link
Contributor

Choose a reason for hiding this comment

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

Again I suggest removing this marginal example.

cdef double[:] quat

for ind in range(num_rotations):
if self._quat[ind, 3] < 0: # w > 0 to ensure 0 <= angle <= pi
Copy link
Contributor

Choose a reason for hiding this comment

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

I would probably not introduce quat variable just adjust the formula for mrps components depending on the case (using if in loop by i or before loop).



def test_from_1d_single_mrp():
mrp = [0, 0, 0.41421356]
Copy link
Contributor

Choose a reason for hiding this comment

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

Prefer using "exact" constants, like 0.25 * np.pi.

I guess unity MRP is a better testing example.



def test_from_generic_mrp():
mrp = [
Copy link
Contributor

Choose a reason for hiding this comment

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

Ugly formatting (and more cases below).

[0, 0, 0, 1]
])
assert_array_almost_equal(
Rotation.from_mrp(mrp).as_quat(),
Copy link
Contributor

Choose a reason for hiding this comment

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

Questionable formatting here as well.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Resolved, along with others. Do you guy prefer to hit the "resolve conversation" button when a comment is addressed or just leave it?

@nbelakovski
Copy link
Contributor Author

nbelakovski commented Aug 12, 2020

I've updated the PR to address comments. In the update I have

  • Made the definition more substantial and made it consistent across from_mrp and as_mrp with one small difference in as_mrp to denote that the rotation returned is <=180
  • Modified the examples to use unity MRP and Euler angles where it made sense
  • Renamed the ss variable. Renaming it to denom wasn't the best choice since it is also used in the numerator for calculating the scalar part of the quat, so I just went explicit and named it "mrp_squared_plus_1"
  • Modified tests to have more consistent formatting (not sure if I got what you were after, but it looks more consistent to me now) and use the unity mrp in some cases
  • Kept the "stack with a single rotation" examples despite the suggestion to remove it because the other representations all have the same example.
  • Minor fixes (corrected the shape issue pointed out by @nmayorov , dropped use of the word "shadow" and updated docs accordingly, removed quat in the as_mrp calculation and just used a sign variable to make sure the set of MRPs being calculated corresponded to a rotation of <=180)

"""Initialize from Modified Rodrigues Parameters (MRPs).

MRPs are a 3 element attitude description commonly used in the aerospace
industry. It is similar to the rotation vector, except instead of scaling
Copy link
Contributor

Choose a reason for hiding this comment

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

Maybe just following definition in from_rotvec is a good idea.

MRPs is a 3 dimensional vector co-directional to the axis of rotation and whose magnitude is equal to tan(theta / 4), where theta is the angle of rotation (in radians) [1].

the rotation vector by the angle of rotation, is it scaled by tan(theta/4).

MRPs have a singuarity at 360 degrees which can be easily avoided by ensuring
that the angle of rotation is always <= 180
Copy link
Contributor

@nmayorov nmayorov Aug 13, 2020

Choose a reason for hiding this comment

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

is always <= 180 -> doesn't exceed 180 degrees?

industry. It is similar to the rotation vector, except instead of scaling
the rotation vector by the angle of rotation, is it scaled by tan(theta/4).

MRPs have a singuarity at 360 degrees which can be easily avoided by ensuring
Copy link
Contributor

Choose a reason for hiding this comment

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

Remove "easily"?

-------
rotation : `Rotation` instance
Object containing the rotations represented by input rotation
vectors.
Copy link
Contributor

Choose a reason for hiding this comment

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

Not rotation vectors.

Initialize a single rotation:

>>> r = R.from_mrp([0, 0, 1])
>>> r.as_euler('xyz')
Copy link
Contributor

Choose a reason for hiding this comment

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

Add degrees=True or use as_rotvec otherwise.

>>> r.as_euler('xyz').shape
(2, 3)

It is also possible to have a stack of a single rotation:
Copy link
Contributor

Choose a reason for hiding this comment

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

I saw your consistency argument.

Not sure, I'm fine with either keeping or removing it.

sign = -1 if self._quat[ind, 3] < 0 else 1

for i in range(3):
mrps[ind, i] = sign * self._quat[ind, i] / (1 + sign * self._quat[ind, 3])
Copy link
Contributor

Choose a reason for hiding this comment

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

Moving cdef double denominator = 1 + sign * self._quat[ind, 3] outside of the loop is a reasonable option.

[1, 2, 2],
[1, -1, 0.5],
[0, 0, 0]
])
Copy link
Contributor

@nmayorov nmayorov Aug 13, 2020

Choose a reason for hiding this comment

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

Super minor, but I'd put ]) on the previous line (in all places).

@nbelakovski
Copy link
Contributor Author

New update with latest feedback:

  • Used @nmayorov's provided definition for MRPs and cleaned up the language a little
  • Set Euler examples to use degrees
  • Moved denominator calculation outside of for loop in as_mrp
  • Other minor fixes (fixed copy/paste error spotted by @nmayorov, fixed formatting in test_rotation)

"""Initialize from Modified Rodrigues Parameters (MRPs).

MRPs is a 3 dimensional vector co-directional to the axis of rotation and whose
magnitude is equal to tan(theta / 4), where theta is the angle of rotation (in radians) [1].
Copy link
Contributor

Choose a reason for hiding this comment

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

I suggest adding double ticks (``) around "tan(theta / 4)" and "theta" to render it as code/expression.

The reference should be done as [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.

Done

# to a rotation of <= 180
sign = -1 if self._quat[ind, 3] < 0 else 1

denominator = 1 + sign * self._quat[ind, 3]
Copy link
Contributor

@nmayorov nmayorov Aug 14, 2020

Choose a reason for hiding this comment

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

I suggest adding a test which will show that a returned MRPs will be correct and have magnitude <= 180 degrees, when the rotation was created with magnitude >= 180 degrees (through Euler or rotation vector). The scalar quaternion component must be negative in this case.

Reasoning - to test that this mathematical transformation does in fact work, because it's not obvious.

Added: but actually it might be considered obvious, because we just flip quaternion sign here. But still I think it won't hurt to test it.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I think that's a good idea, test added

This is an attitude description commonly used in aerospace.
@nbelakovski
Copy link
Contributor Author

My latest changes seem to be causing some Windows failures, but I'm having trouble figuring out exactly what's going on. Any help?

@tylerjereddy
Copy link
Contributor

My latest changes seem to be causing some Windows failures, but I'm having trouble figuring out exactly what's going on. Any help?

Not related to your changes, I should be able to fix shortly

@nmayorov
Copy link
Contributor

It looks very good now, I'm merging.

Thank you for the contribution!

@nmayorov nmayorov merged commit 047ab37 into scipy:master Aug 17, 2020
@nbelakovski nbelakovski deleted the mrp branch August 17, 2020 21:24
@ilayn ilayn added this to the 1.6.0 milestone Aug 17, 2020
@nbelakovski
Copy link
Contributor Author

It looks very good now, I'm merging.

Thank you for the contribution!

Great working with you, looking forward to future contributions!

@ilayn
Copy link
Member

ilayn commented Aug 17, 2020

Thanks @nbelakovski . Please have a moment to boast about this in the release notes ;)

https://github.com/scipy/scipy/wiki/Release-note-entries-for-SciPy-1.6.0

swallan pushed a commit to swallan/scipy that referenced this pull request Sep 6, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

enhancement A new feature or improvement scipy.spatial

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants