using quaternion.quaternion on a numpy array

[calocedrus]

Hello,
Along my way to find a method to apply the same rotation to an array of 3D points (x,y,z coordinates), I've experimented with the constructor quaternion.quaternion(0,*one_by_3_np_array).
In order to be able to directly construct a "pure" quaternion (no real part) from my array of 3D points, I thought I could do:

# define my test XYZ array
myxyz = np.random.rand(100,3)
# make it into a 100 rows quaternion array
myxyz_quat = quaternion.quaternion(0,*myxyz) 

But found out that obviously this doesn't work because unpacking my array will result in more than the 4 floats expected by the constructor.
I can define a function to build my quaternion line by line but doubt this would be efficient.
The solution I've came up with is to use quaternion.as_quat_array:

def npXYZary_to_quaternion (npary):
    tmpary = np.zeros((npary.shape[0],4))
    tmpary[:,1:] = npary
    qnpary = quaternion.as_quat_array(tmpary)
    return qnpary

But could it be possible to have the quaternion constructor take a Nx3 (Nx4 as well) numpy array ?
Ludovic

[moble]

I can make the constructor take three arguments instead of four (and I've done so in 53d8322). But this won't solve your problem, because the quaternion constructor has never taken arrays of any shape. In particular, quaternion.quaternion(0, *myxyz) just makes no sense at all — it's just a single 0, followed by some multi-dimensional array of vectors.

The solution you've come up with is the best possible one, if you need to convert to quaternions. The reason is that if you've created your data in Nx3 shape, you will necessarily have to copy the data into Nx4 shape. This isn't a limitation of my code; this is a requirement of numpy and how computers work. If possible, it's a better idea to just create your data in Nx4 shape to begin with. And if not, there's just no getting around the fact that the data will have to be copied.

However, if you're rotating multiple vectors with the same quaternion, it's fundamentally more efficient to convert the quaternion to a matrix, and then do the usual matrix multiplication. (Again, this not a limitation of this code; it's just how the math works out.) So in your case, I would do something like this:

import numpy as np
import quaternion

# Some quaternion I've made up: rotation by 0.1 radians about z axis
q = np.exp(quaternion.z * 0.1 / 2)

# Convert quaternion to matrix form
m = quaternion.as_rotation_matrix(q)

# Your vector data
myxyz = np.random.rand(100,3)

# Your rotated vector data
# (this is just matrix multiplication for each of the 100 vectors)
myXYZ = np.dot(m, myxyz.T).T

If you want to check that the result is right, you can do something like this:

>>> print((q * quaternion.quaternion(*myxyz[0]) * q.inverse()).vec)
[ 0.3589965   0.12589602  0.82435895]
>>> print(myXYZ[0])
[ 0.3589965   0.12589602  0.82435895]

You can also do multiple rotations at the same time. The code above works in the same way if you use (for example)

q = np.exp([quaternion.z * 0.1 / 2, quaternion.x * 0.2 / 2])

However, the final result might need a little more reshaping.

[calocedrus]

Thank you for the code sample and feedback, you could add this as a usage example in your readme. Your comment at the bottom of that readme influenced me: I wanted to avoid using a rotation matrix. It's now obvious that it's more efficient to use one in my case!
This is how I was doing, given my array of xyz, defined as a Nx3 numpy array (I have no choice on its input format):

import numpy as np
import quaternion
from math import sin, cos, sqrt

d2r = np.pi/180

#-------------------------------------------
### define the relevant functions
#-------------------------------------------
# makes a rotation quaternion from a rotation axis and an angle
# exponential notation (as in your reply) would be more compact and maybe faster
def qrot_axis_angle(rotaxis,rotangle):
	# angle in radian
	qrotangle = d2r*rotangle/2.0
	# normalize the rotation axis
	magrotaxis = np.sqrt(rotaxis.dot(rotaxis))
	nrotaxis = rotaxis/magrotaxis
	q_rot = quaternion.quaternion(0,*nrotaxis)
	qrot = cos(qrotangle) + sin(qrotangle)*q_rot
	return qrot

# convert a numpy array to pure quaternion
# input array must float, of size Nx3
def npXYZary_to_quaternion (npary):
	tmpary = np.zeros((npary.shape[0],4))
	tmpary[:,1:] = npary
	qnpary = quaternion.as_quat_array(tmpary)
	return qnpary

# return an array of rotated XYZ points or vectors
def rot_ptary(npary,rotquat):
	tmp_qary = npXYZary_to_quaternion(npary)
	rotd_tmp_qary = rotquat * tmp_qary * rotquat.conj()
	rotd_ary = quaternion.as_float_array(rotd_tmp_qary)
	return rotd_ary[:,1:]

Then:

# example of input array
xyz_ary = np.random.rand(100,3)
#define a quaternion corresponding to a 45 deg rotation about the x-axis
qrot = qrot_axis_angle(np.array([1,0,0]),45)
# rotate 
rotd_xyz_ary = rot_ptary(xyz_ary,qrot)

(my code formatting doesn't look as good as yours, especially the comments!)
PS:

if you've created your data in Nx3 shape, you will necessarily have to copy the data into Nx4 shape. This isn't a limitation of my code [...]

My original question wasn't aimed at suggesting any limitation of your code.

[moble]

Your comment at the bottom of that readme influenced me: I wanted to avoid using a rotation matrix.

I'm not sure I understand what you're talking about. The only comment I see is about Euler angles. Just to be clear, Euler angles should absolutely never be used — ever. Period. Never.

Rotation matrices are completely separate things. They should usually be avoided, because multiplying matrices is much less efficient than multiplying quaternions, and after a few multiplications the errors will multiply, giving you "rotation" matrices that aren't really rotations. But if you're starting with a quaternion and converting to a matrix, there's not a lot of room for error to creep in, so they're fine. The one time when you actually should use a matrix is as in your case, where you want to rotate multiple vectors.

Also, in case you missed it, I'll also point out the new function quaternion.rotate_vectors, which does the rotations and handles the various possible dimensions nicely, with no need to copy or reshape data. So your code could look as simple as

import numpy as np
import quaternion

xyz_ary = np.random.rand(100,3)
qrot = np.exp(quaternion.quaternion(1, 0, 0) * 45 * (np.pi/180) / 2)
rotd_xyz_ary = quaternion.rotate_vectors(qrot, xyz_ary)

[calocedrus]

A moment of ridicule. I've been away from rotations in space for too many years and automatically associated Euler angles to rotation matrices, extrapolating the former to the later; the bottom note of your readme is about Euler angle indeed. I had read your commit to add a function to rotate vector.

[Nirkan]

Hi there ,
I applied these methods to rotate an ellipsoid in 3D(3 dimensional array ) on an astronomical data. But when i try to view the image in 3D I cannot see that any rotation has occurred . I do see my numpy array gets changed but no rotation . I wanted to rotate 3D numpy using quaternion. Is there other method how I can use quaternion.

Niraj

[moble]

@Nirajkandpal I'd be happy to help get you started, but you haven't included enough information. Also, this is issue is not the right place to ask the question. I'm going to lock this conversation. Please open a new issue to ask your question again, with more information.

You'll need to include a minimal working example for me to understand what you've actually tried, what's happening, and how it's different from what you expected. For example, look at this code snippet I've copied from above:

import numpy as np
import quaternion

# Some quaternion I've made up: rotation by 0.1 radians about z axis
q = np.exp(quaternion.z * 0.1 / 2)

# Convert quaternion to matrix form
m = quaternion.as_rotation_matrix(q)

# Your vector data
myxyz = np.random.rand(100,3)

# Your rotated vector data
# (this is just matrix multiplication for each of the 100 vectors)
myXYZ = np.dot(m, myxyz.T).T

Here, q is the quaternion you use to do the rotation, and myxyz is the array of 3D data. Does that not do what you want?