Add nD rotation matrix generation

[kne42]

Description

Add function for the generation of a matrix to rotate a vector to coincide with another.

Checklist

[It's fine to submit PRs which are a work in progress! But before they are merged, all PRs should provide:]

• Clean style in the spirit of PEP8
• Docstrings for all functions
• Gallery example in ./doc/examples (new features only)
• Unit tests

[For detailed information on these and other aspects see scikit-image contribution guidelines]

References

• taken from Expand Gabor Kernel Generation to nD #2739
• Zhelezov OI. 2018. One Modification which Increases Performance of N-Dimensional rotation Matrix Generation Algorithm. International Journal of Chemistry, Mathematics, and Physics 2(2): pp. 13-18. https://dx.doi.org/10.22161/ijcmp.2.2.1
• ter Haar Romeny BM. 2003. Front-End Vision and Multi-Scale Image Analysis. Computation Imaging and Vision Vol. 27. Springer: Dordrecht, Netherlands: pp. 37-51. https://doi.org/10.1007/978-1-4020-8840-7_3
• https://en.wikipedia.org/wiki/Rotation_matrix#In_three_dimensions

For reviewers

(Don't remove the checklist below.)

• Check that the PR title is short, concise, and will make sense 1 year
  later.
• Check that new functions are imported in corresponding __init__.py.
• Check that new features, API changes, and deprecations are mentioned in
  doc/release/release_dev.rst.

[pep8speaks]

Hello @kne42! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

• In the file skimage/transform/tests/test_orientation.py:

Line 7:80: E501 line too long (80 > 79 characters)
Line 23:1: E302 expected 2 blank lines, found 1
Line 52:1: E302 expected 2 blank lines, found 1

Comment last updated at 2020-01-23 21:58:18 UTC

[codecov-io]

Codecov Report

Merging #3125 into master will decrease coverage by 1.78%.
The diff coverage is 99.12%.

@@            Coverage Diff            @@
##           master   #3125      +/-   ##
=========================================
- Coverage   87.89%   86.1%   -1.79%     
=========================================
  Files         325     340      +15     
  Lines       27490   27477      -13     
=========================================
- Hits        24163   23660     -503     
- Misses       3327    3817     +490

Impacted Files	Coverage Δ
skimage/transform/__init__.py	100% <100%> (ø)	⬆️
skimage/transform/tests/test_orientation.py	100% <100%> (ø)
skimage/transform/_orientation.py	98.38% <98.38%> (ø)
skimage/io/_plugins/gdal_plugin.py	0% <0%> (-62.5%)	⬇️
skimage/io/_plugins/gtk_plugin.py	0% <0%> (-30.56%)	⬇️
skimage/io/_plugins/qt_plugin.py	0% <0%> (-21.12%)	⬇️
skimage/io/_plugins/simpleitk_plugin.py	63.63% <0%> (-18.19%)	⬇️
skimage/io/_plugins/imread_plugin.py	69.23% <0%> (-15.39%)	⬇️
skimage/io/_plugins/q_histogram.py	0% <0%> (-15%)	⬇️
skimage/__init__.py	53.22% <0%> (-14.64%)	⬇️
... and 126 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 6059644...c43990b. Read the comment docs.

[emmanuelle]

About 3D rotations, did you see that a GSoC student is working on this for Scipy? https://mail.python.org/pipermail/scipy-dev/2018-May/022831.html
@kne42 @jni

[jni]

I did not see that! @kne42 the blog linked there looks like a good resource!

Thanks @emmanuelle!

[kne42]

@emmanuelle @jni @soupault @stefanv This PR is (finally) ready for review! :)

[jni]

@kne42 Wow! I will do a proper review tomorrow as part of the scikit-image Paris sprint, but on first glance: 😍😍😍

[stefanv]

This looks great, thanks @kne42. Questions: what is the situation around rotations in SciPy after their GSoC? How would this typically be used in image processing, i.e. can we come up with a good gallery example?

[kne42]

Sorry for the late reply, I keep forgetting about this 😅

So briefly, there are two rotation-centric utilities in SciPy: an nD angle-plane rotation in ndimage and the new Rotation object from GSoC. This new object offers conversion between angle-axis, vector-vector, matrix, quarternion, and (I think) euler angle rotations in 3 dimensions. Internally, everything is stored as a quaternion.

This proposed addition is an nD implementation of a vector-vector -> rotation matrix conversion. E.g. generate the matrix to rotate one vector to point in the same direction as the other.

This is probably most useful internally as a building block for nD algorithms that require rotation, but I’m sure users will be able to find direct use cases. Since I’m largely a beginner when it comes to image processing, the first example that comes to my mind would be with a rotation variant feature detector that outputs features in a plottable format. In this case, one could first rotate the image using the rotation matrix to the desired direction, then rotate the output features using the matrix transpose and plot them on the original image.

[alexdesiqueira]

Hey @kne42, thank you very much for your contribution! I am reading the code and for now I have some small changes to ask.

[kne42]

this is ready for re-review @alexdesiqueira

[alexdesiqueira]

Thank you for your work @kne42, and sorry for not checking this before.
@scikit-image/core (especially the Guardian of the Dimensions, @jni) LGTM right now! Where are we on this one?

[jni]

Damn, this is great work, @kne42! It's crazy what a bit of mulling over time will do, this went from being quite arcane to me to being crystal clear! I've made some very minor suggestions that should be easy to fix, if you can get them done soon I can push that green button! 🎉

[carterbox]

I previously implemented an ND rotation algorithm for polytope, so I'd like to point out that there is more than one rotation that aligns start and end vectors. For N > 2, not only are there multiple rotations, but these rotations do not all leave the space in the same orientation. Therefore, I'm not sure that this API allows sufficient control over the returned rotation matrix for the desired application.

As an example, if you want to rotate a pencil from pointing up to pointing down, you can flip the pencil way from you or rotate it to the right. These two rotations align the start vector with the end vector, so they could be valid matrices returned by compute_rotation_matrix. However, the pencil ends in different orientations, so which one will you get?

[rfezzani]

LGTM, despite some comments. Thank you for this contribution @kne42 ;-)

[rfezzani]

Suggested change

	unit_vector : (N, ) array
	unit_vector : (N, ) ndarray

[kne42]

changed to array-like

[rfezzani]

Suggested change

	src : (N, ) array
	src : (N, ) ndarray

[kne42]

changed to array-like

[rfezzani]

Suggested change

	dst : (N, ) array
	dst : (N, ) ndarray

[kne42]

changed to array-like

[kne42]

This is ready for re-review.

@carterbox Honestly, I'm not a mathy person so idk. I just took an algorithm from some paper I kinda understood and implemented it here.

[jni]

@carterbox thank you so much for chiming in! nD rotations have been a headache for us for a long time. 😅 Looking at polytope might be very useful for us.

@kne42 does the paper not discuss this problem at all?

[kne42]

@jni Well here's how I understood @carterbox's comment:

|
v
-->

These vectors are not aligned nor face the same direction

-->
<--

These vectors are aligned but do not face the same direction

-->
-->

These vectors are aligned and face the same direction.

---

The paper does say that it makes both vectors face the same direction (see section II):

for any vectors of the same dimension X and Y, the matrix M = compute_rotation_matrix(X, Y) can be computed such that the vector Y2 = M @ X has the same norm as X and the same direction as Y

so assuming that the algorithm is right (which I do not have the skills to verify) and I implemented it correctly, yes, it should handle direction correctly

---

sidenote: I'm pretty sure using the term "coincide" is wrong (as I have been putting it in some of the docs). what is the correct mathematical term for when a vector has the same direction as another?

[jni]

@kne42 no, the comment is more subtle than that, and you can't illustrate it with "point" vectors. Think now about rotating not a vector, but a human, 180 degrees. If you flip forward, your head will point towards the ground, but you will be facing backwards. If you flip right, your head will still point down, but you will be facing forward. So there is actually an infinite number of rotation matrices that will make your head point down, but the algorithm only returns one.

[kne42]

@jni Ohhhhhhhhhhh I see. I just skimmed over the paper again and I do not see anything mentioning this.

@carterbox what are your suggestions? I can implement a different algorithm or maybe someone who understands math can help me out...

[carterbox]

@jni, Yes! That's what I meant. Thanks for clarifying.

@kne42, My suggestions are to address the following edge cases:

1. Check that none of the vectors src, dst are the zero-vector because you can't align a vector with the zero-vector. It looks like your implementation would fail during normalization due to a division by zero. Not sure if you want explicitly check for zero-vectors to prevent that error.
2. Check that your implementation returns the identity matrix when src and dst are parallel.
3. Warn when src, dst are anti-parallel. In this case, there are many solutions, but you only return one. Alternatively, as in polytope, you can define an interface which takes a dst vector half-way along the desired rotation. This approach removes the ambiguity about how to rotate by 180 degrees, but it also might require the user to do more work when most of the time they are not rotating by 180 degrees? You could also just explain this alternative approach in the warning. Choose the approach that best fits the application here.

[kne42]

@jni @carterbox How often do you think this would be an issue for someone? Does this warrant a warning over including something about this in the notes instead?

@carterbox So would the instructions be:

M_halfway = compute_rotation_matrix(X, Y_halfway)
M = M_halfway @ M_halfway

?

[carterbox]

From the python docs:

warnings.warn() in library code if the issue is avoidable and the client application should be modified to eliminate the warning

I think a warning is appropriate because I'm not sure that the half-angle approach is compatible with homogeneous coordinates. If homogeneous_coords=True, then the space would also be translated/scaled twice? I'm actually a little rusty when it comes to homogeneous coordinates. I just remember that the order of operations matters when you mix translations and rotations. 🤷‍♂️

[jni]

@carterbox the thing is, if I read the implementation correctly, it is almost never the shortest rotation: all rotations go via (1,0, ..., 0). I don't know what happens when you try to rotate (1, 0, 0) to (-1, 0, 0), but that would be an interesting test case, @kne42! =)

[carterbox]

@jni I don't understand what you mean by "never the shortest rotation". All of the rotations are expressed as matrices of the same size, so they are all 'equally short' (computationally speaking).

[carterbox]

I didn't mean to kill the momentum of this pull request! 😅 I just wanted to point out that there should be some extra checking of the inputs and because the solution to the rotation problem is non-unique when the start and end vectors are parallel.

[jni]

@carterbox 😂 sometimes, complete discussion can have this effect. =D

The problem is that the issue you identified is in no way unique to anti-parallel rotations. For example, to point my head to the right, I can simply roll right, then my head will point right and I will face forward, or I can pitch forward and then roll (my reference frame) / yaw (external reference frame) right. Then I will be facing down, but head pointing right as before.

Unless I'm mistaken, for any rotation, there is an infinite number of ways to get there, and they will all leave the object/space in a different orientation. This PR always rotates from the source vector to the 0th axis, then to the target vector, which is different from going from source to target directly.

[carterbox]

@kne42 @jni Let's back up to my earliest question about the desired application. This function claims to align two vectors in an ND space. If the final orientation of the space doesn't matter, and it only matters that the vectors are aligned. Then, just add some documentation for people not familiar with the non-uniqueness of this problem, and that's it. The function does what it claims. Done.

If the purpose of this function is to provide a general API for rotating ND spaces, then have a look at the docstring I wrote for polytope. In those docs, I note that for N>3 dimensions and when using only two vectors, it is only possible to unambiguously define simple rotations (rotations whose motion can be projected into a single plane). Because the polytope API only promises simple rotations, using the half-rotation end-vector is enough to avoid the only ambiguous case for defining a plane (the case when the start and end vectors are co-linear). For simple rotations, the rotation matrix generating algorithm should be generating a rotation matrix in which all motion can be projected into a single plane of rotation. (I feel like there may be an easy test for this which involves row reduction or principal component analysis.)

If the user has a fully defined mapping from one orientation to another, a rotation function isn't even necessary because one can trivially construct the reorientation matrix by mapping the basis vectors the starting space to the ending space.