ENH: Extend Rotation.align_vectors() to allow an infinite weight, as well as minimum-distance rotations for single vectors #17542
Conversation
scottshambaugh
force-pushed
the
rotation_align_constrain
branch
4 times, most recently
from
735053b
to
a222852
Compare
scottshambaugh
force-pushed
the
rotation_align_constrain
branch
17 times, most recently
from
c440862
to
70d01a5
Compare
scottshambaugh
changed the title
scottshambaugh
force-pushed
the
rotation_align_constrain
branch
from
b5b6f2c
to
5b348ce
Compare
scottshambaugh
force-pushed
the
rotation_align_constrain
branch
2 times, most recently
from
2d31329
to
7dae5f3
Compare
scottshambaugh
force-pushed
the
rotation_align_constrain
branch
10 times, most recently
from
74997cf
to
e4aa873
Compare
|
I sketched out an example and now understand the calculation for the secondary rotation angle. Added some comments in the code to explain it, and I think this is in quite a good state now (barring the one rssd edge case in the comment above). |
| # Calculate vector components for the angle calculation. The | ||
| # angle phi to rotate a single 2D vector C to align to 2D | ||
| # vector A in the xy plane can be found with the equation | ||
| # phi = atan2(cross(C, A), dot(C, A)). |
There was a problem hiding this comment.
Yes, it's true. And if you have many vectors the same formula can be applied, but with sums over all vectors (possibly weighted). It can be derived rigorously (in 2D) by finding a stationary point of a function of one argument (the angle phi).
The trick in 3D is to consider projections of vectors into the plane. Let's denote the normal unit vector as u (in our case is a_pri). Let Cp and Ap be the orthogonal projections. Because they still live in 3D we have to project the cross product onto u, by using (Cp x Ap) . u. For dot product we don't have to change anything and simply Cp . Ap.
Now is we substitute Ap = A - (A . u) u and Cp = C - (C . u) u (remove the component along the u to get the orthogonal projection) we will get the formulas you've used here. They are just an alternative form to get the angle phi without explicitly computing Cp and Ap (by the way I think the code where Cp and Ap are computed would be also fine).
scottshambaugh
force-pushed
the
rotation_align_constrain
branch
6 times, most recently
from
da9e2ff
to
e8fa19f
Compare
scottshambaugh
force-pushed
the
rotation_align_constrain
branch
from
57a398b
to
fcd4714
Compare
|
Happy to merge this pull request. I've updated the release notes. Check if you agree with the description and whether you think the new minimal rotation should be mentioned in "Backwards incompatible changes" (I haven't done it). Ping me if you'll have other pull requests to consider. |
|
Thank you for the collaboration on this, I think the result is much better for it! I have no other pull requests at this time. The release notes look good. I would not consider this backwards incompatible behavior, since if you tried to align shape (3,) vectors previously you would get an error, and if you tried to align shape (1, 3) vectors you get the below warning. So the previous behavior had no uniqueness guarantees. |
These are good points, thank you! |
Reference issue
Closes #17462
Closes #16880
What does this implement/fix?
This implements the "align-constrain" algorithm for defining a rotation. The user passes in primary and secondary vectors in reference frames A and B. The algorithm calculates a rotation from B to A that exactly aligns the primary vectors, and then rotates about that axis to most closely align the constrained secondary vectors.
This algorithm is equivalent to
align_vectors()with sets of two vectors, an infinite weight on the first pair of vectors, and a finite nonzero weight on the second pair of vectors. There is a test included here to verify that internal consistency.Additional information
I would appreciate some feedback on how to handle the case where primary and secondary vectors are parallel. Currently I have it defined such that if the secondary vector is parallel to the primary to within an angular toleranceatol, the rotation is unconstrained and so the vector [1, 0, 0] is arbitrarily chosen as the secondary (followed by [0, 1, 0] if that is also parallel). Another option would be to error out if the vectors are parallel. A third option would be to error out if the user passes inatol = Noneand choose the secondary axis otherwise.Update: The functionality is modified so that in the instance of parallel pairs, only the primary vectors are aligned via a rotation that has no change in angle about the primary axis.