vector icon indicating copy to clipboard operation
vector copied to clipboard

`scipy.spatial.transform.Rotation` Integration

Open Elyasnz opened this issue 1 year ago • 3 comments

Vector Version

1.1.1.post1

Python Version

3.10.9

OS / Environment

OS: KAli Linux 2023.3 Environment: scipy==1.11.1

Describe the bug

First of all thanks for the awesome work on the library

there are some errors with Rotation.apply

these are some ways to reproduce

from scipy.spatial.transform import Rotation
from vector import VectorObject3D, VectorNumpy3D

v_object = VectorObject3D(x=1, y=0, z=0)
v_array = VectorNumpy3D({"x": [1], "y": [0], "z": [0]})
r = Rotation.from_quat((0, 0, 1, 1))

# _r = r.as_quat()
# print(v.rotate_quaternion(_r[3], _r[0], _r[1], _r[2]))
# Returns: VectorObject3D(x=0.0, y=0.9999999999999998, z=0.0)

# print(r.apply(v_object))  # Raises: IndexError(tuple index out of range)
print(r.apply(v_object.__array__().tolist()))
# Returns: [0. 1. 0.]

# print(r.apply(v_array))  # Raises: ValueError(Expected input of shape (3,) or (P, 3), got (1,).)
print(r.apply(v_array.tolist()))
# Returns: [[0. 1. 0.]]

I think it would be a good idea to add rotate method to Vector that accepts the Rotation instance as input and applies the rotation on Vector

from scipy.spatial.transform import Rotation

def rotate(self, rot: Rotation):
    pass

Any additional but relevant log output

No response

Elyasnz avatar Oct 11 '23 16:10 Elyasnz

H, @Elyasnz, thanks for opening an issue! This looks more like a feature request (an out of scope one). I am not sure if there are any plans to work on this, but I'll tag @jpivarski to check if this is something that we would want in vector.

Saransh-cpp avatar May 16 '24 15:05 Saransh-cpp

This is a bit like adding NEP-18 overloads for Python functions. It would, however, be the only interaction that Vector has with SciPy, so it would seem to stand out oddly.

It wouldn't be difficult to add this feature. Vector does not yet have a rotate method, so there's nothing to interfere with, though we might want to use this name in the future with a different interface. Since this is Python, the dispatching by rotate's signature would have to be manual, inside the rotate function, if we ever decide to do that.

scipy.spatial.transform.Rotation can be expressed 6 different ways.

>>> r = Rotation.from_quat([0, 0, np.sin(np.pi/4), np.cos(np.pi/4)])
>>> print("\n".join(name for name in dir(r) if name.startswith("as_")))
as_davenport
as_euler
as_matrix
as_mrp
as_quat
as_rotvec

and a Vector can accept 3 of these:

>>> v = vector.obj(x=1, y=2, z=3)
>>> v.transform3D         # takes a 3x3 matrix
<bound method Spatial.transform3D of VectorObject3D(x=1, y=2, z=3)>
>>> v.rotate_euler        # takes Euler angles
<bound method Spatial.rotate_euler of VectorObject3D(x=1, y=2, z=3)>
>>> v.rotate_quaternion   # takes quaternions
<bound method Spatial.rotate_quaternion of VectorObject3D(x=1, y=2, z=3)>

(I don't know if either Davenport or MRP, Modified Rodrigues Parameters, correspond to the nautical yaw-pitch-roll system in v.rotate_nautical, for a fourth.)

For instance, using a 3×3 matrix:

>>> m = r.as_matrix()
>>> v.transform3D({x + y: m[i, j] for i, x in enumerate("xyz") for j, y in enumerate("xyz")})
VectorObject3D(x=-1.9999999999999998, y=1.0000000000000004, z=3.0)

Vector is looking for an object with named fields, like the dict above, because arrays of vectors will need to broadcast the transformation to all elements. We need to make sure that a matrix (or Euler angles, etc.) gets broadcast as one matrix for every vector, rather than interpreting the matrix as an array of scalars. This seems to be working, though:

>>> vs = vector.VectorNumpy3D(
...     np.random.normal(0, 1, (10000, 3)).view([("x", np.float64), ("y", np.float64), ("z", np.float64)])
... )
>>> vs.transform3D({x + y: m[i, j] for i, x in enumerate("xyz") for j, y in enumerate("xyz")})
VectorNumpy3D([[( 0.08454731, -0.31750909,  0.26364239)],
               [( 2.22012436, -0.3274842 , -0.95178527)],
               [( 0.06190065,  0.98422746, -2.1653507 )],
               ...,
               [(-0.03083845, -0.67911075, -0.07757298)],
               [(-0.42084917,  1.05904514,  0.03699855)],
               [( 1.65960551, -0.20146602, -1.4939052 )]],
              dtype=[('x', '<f8'), ('y', '<f8'), ('z', '<f8')])

So I don't see any blockers to doing this, other than

  1. filling the namespace with more interfaces to maintain,
  2. explaining why we use SciPy for this one method and no others,
  3. deciding whether to convert all SciPy rotations into matrices or if there's some numerical-stability reason to use other methods sometimes.

It's not a high-priority item, and if @Elyasnz needs this capability, the idioms that I've written above will work.

jpivarski avatar May 16 '24 17:05 jpivarski

Hi @jpivarski thanks for the detailed response the solutions you provided are all so good and also I DO agree that it is not high-priority but it would be awesome to have the feature inside scipy

Elyasnz avatar May 18 '24 06:05 Elyasnz