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

Allow rotation of Vector3d object by Rotation object. #453

Open
argerlt opened this issue May 19, 2023 · 3 comments
Open

Allow rotation of Vector3d object by Rotation object. #453

argerlt opened this issue May 19, 2023 · 3 comments
Labels
documentation-examples Consider making into a documentation example

Comments

@argerlt
Copy link
Contributor

argerlt commented May 19, 2023

I am writing a virtual diffractometer using Orix, and it would be useful to allow rotation of a vector by a rotation, ideally by doing:

import numpy as np
from orix.vector.vector3d import Vector3d
from orix.quaternion import Rotation

# would also be helpful to have a Vector3d.random command
v = Vector3d(np.random.rand(10,4,3)).unit
r = Rotation.random(10,4)

# it would be nice to first change Vector3d.rotate so this command does whats expected
r_new = v.rotate(r)
# and would be better if it was called like this, using left hand rotation rules
r_new_2 = r*v

Three points about this.

First, Vector3d.rotate currently uses Vector3d and angle pairs to do a rotation. What is the best way to deal with this? should rotate accept either rotation or vector/angle pairs, or should users be expected to cast vector/angles to rotations prior to the call?

Second, this loosely ties in with the points about #139 and #249, which can be broadly summarized as "what are the best universal rules for how any two orix objects use __mul__, and how should we implement those rules in the code" I think this use case is simple:

  1. follow @pc494 suggestion to always follow numpy broacasting rules.
  2. only allow left-rotation multiplication
  3. ignore symmetry by default, but give option to include if desired.
    But I want to know if anyone (especially @hakonanes , @pc494 , and @harripj ) disagree

Third, there are two ways to do the rotation, and I think the numpy-quaternion.rotate_vectors notes sum them up perfectly:

    Notes
    -----
    For simplicity, this function converts the input quaternion(s) to matrix form,
    and rotates the input vector(s) by the usual matrix multiplication.  As noted
    above, if each input quaternion is only used to rotate a single vector, this is
    not the most efficient approach.  The simple formula shown above is faster than
    this function, though it should be noted that the most efficient approach (in
    terms of operation counts) is to use the formula

      v' = v + 2 * r x (s * v + r x v) / m

    where x represents the cross product, s and r are the scalar and vector parts
    of the quaternion, respectively, and m is the sum of the squares of the
    components of the quaternion.  If you are looping over a very large number of
    quaternions, and just rotating a single vector each time, you might want to
    implement that alternative algorithm using numba (or something that doesn't use
    python).

so, we could do the math the easy way or the efficient way.

Thoughts?

@hakonanes
Copy link
Member

hakonanes commented May 19, 2023

A virtual diffractometer sounds cool, please let us know when you make it available so we can list in among related projects in the docs.

r_new_2 = r*v

This is the canonical way of rotating vectors in orix. In this example, one single vector is rotated by many orientations to return many vectors. Is this what you're after? If so it is already available? I realise that it is not covered specifically in the API reference of Vector3d, we should add this to the class docstring example...

I never use Vector3d.rotate() myself. As the method docstring states, it's a convenience method to use specifically rotation axes and angles to rotate a single vector. I'd much rather remove this method altogether than to allow Vector3d.rotate(Rotation)... What do you think? If we cover the canonical way of rotating a vector in the docstring, I think this should be fine.

  1. Currently, only left-rotation is possible, as your example above.
  2. We have Miller for handling vectors with respect to a crystal's unit cell and its symmetry. These are covered in great detail in the Crystal directions tutorial. Is the functionality you need covered by this class?

Yes, we should add a Vector3d.random().

@argerlt
Copy link
Contributor Author

argerlt commented May 19, 2023

Ah beans, you're right. Sorry, I should have tried this before making an Issue. I saw Vector3d.rotate was a function and didn't even think to try the obvious choice.
Looking in the code, Rotation is even using the efficient computation method instead of the lazy way discussed in numpy-quaternion. hats off to whoever wrote that code (@harripj i think?)

Also, quick test shows it broadcasts correctly for multiple vectors as well, so even better:

from orix.vector.vector3d import Vector3d
from orix.quaternion import Rotation
import numpy as np

v = Rotation([1, 1, 1, 0]).axis
v_repeat = Vector3d(np.tile(v.data, 2).reshape(2, 3))
v_rand = Rotation.random(2).axis

r = Rotation([0, 2, 3, 1])
r_repeat = Rotation(np.tile(r.data, 2).reshape(2, 4))
r_rand = Rotation.random(2)

print(r*v)
print(r*v_repeat)
print(r*v_rand)

print(r_repeat*v_repeat)
print(r_repeat*v_rand)

print(v*r)

As for the rest of your comments:

I'd much rather remove this method altogether ...

I agree, Vector3d.rotate() could be removed and replaced with Rotation.from_axes_angles(axis,angle)*v or just r*v as appropriate in other functions. I could do that fix this weekend and update the API as well.

We have Miller for handling vectors with respect to a crystal's unit cell and its symmetry.

That's a much better way to handle things I think. thank your for alerting me to this.

Yes, we should add a Vector3d.random().

I can add that as well, small addition.

@hakonanes
Copy link
Member

I saw Vector3d.rotate was a function and didn't even think to try the obvious choice.

I'm sure many people do this, which is understandable since we don't show on the Vector3d top reference page how to rotate a vector, aside from this method! We should definitely add this to the class docstring example.

The initial rotation of vector by a quaternion was written by the orix originator I believe, Ben Martineau. Rotation was performed using v' = q * v * q^-1 (according to this formula). @harripj replaced this "hard-coded" implementation with the current numpy-quaternion one.

Removing Vector3d.rotate() means we need to encapsulate rotation in parentheses if we want to chain operations... Also, using operators means we cannot control the outcome with any parameters (say reduction to any "fundamental zone"). For these reasons, I'm less sure of removing this method. I'd leave it for now.

Yes, please go ahead and add Vector3d.random()!

@hakonanes hakonanes added documentation Relates to the documentation documentation-examples Consider making into a documentation example and removed documentation Relates to the documentation labels Jul 26, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
documentation-examples Consider making into a documentation example
Projects
None yet
Development

No branches or pull requests

2 participants