rotations.py quatProduct doesn't work for (4,n) inputs, maybe.

[kevindlewis23]

If rotations.py quatProduct on (4,n) inputs is supposed to be element-wise composition, it does not work. When running, it element-wise (seperately on (4,1) inputs), it matches running on (4,n) inputs on the first output quaternion only.

Can be reproduced with the following code:

import numpy as np
from hexrd import rotations
def test_reproduce():
    n = 3
    # Generate random quaternions
    q1 = np.random.rand(4, n) * 2 - 1
    q2 = np.random.rand(4, n) * 2 - 1

    # Fix to unit rotation quaternions (this works and is tested)
    q1 = rotations.fixQuat(q1)
    q2 = rotations.fixQuat(q2)
    print("q1:")
    print(q1)
    print("q2:")
    print(q2)

    # Product using fixQuat()
    qp = rotations.quatProduct(q1, q2)
    print("qp:")
    print(qp)

    # Elementwise product using fixQuat()
    # Written as a regular loop for readability, can be done inline
    qp_el = None
    for i in range(n):
        a = np.atleast_2d(q1[:, i]).T
        b = np.atleast_2d(q2[:, i]).T
        prod = rotations.quatProduct(a, b)
        if qp_el is None:
            qp_el = prod
        else:
            qp_el = np.hstack((qp_el, prod))
    qp_el = np.array(qp_el)
    print("qp element-wise:")
    print(qp_el)

Example output is:

q1:
[[ 0.6639988   0.45485282  0.83234881]
 [-0.20387788 -0.44343047 -0.28927565]
 [ 0.55001834 -0.70874605 -0.24522852]
 [ 0.46370166  0.30685072 -0.40420048]]
q2:
[[ 0.40840558  0.72671763  0.36315318]
 [-0.30968977 -0.04341781 -0.40514601]
 [ 0.81406359 -0.34699685 -0.69307048]
 [-0.2731256   0.591261   -0.47289511]]
qp:
[[ 0.11305919  0.16252659  0.23675927]
 [-0.23880823 -0.26005055 -0.11490376]
 [-0.96445574  0.66319071  0.62761961]
 [-0.00365778 -0.6827422   0.73269079]]
qp element-wise:
[[ 0.11305919  0.11606473  0.17603383]
 [-0.23880823  0.02941978  0.27810212]
 [-0.96445574  0.9217511   0.69289508]
 [-0.00365778 -0.36883379  0.64153544]]

[kevindlewis23]

Can be easily fixed by just using scipy multiplication, I've already written it in with testcases. But I want to make sure this is actually supposed to do element-wise multiplication.

[donald-e-boyce]

You're right. And by the way, I recently made open source the polycrystal package (in the HEXRD project), which has some overlap with hexrd, and in the quaternion operations in particular. They were done independently but both based on the matlab ODFPF package.

So here is a script that confirms your script.

import numpy as np

from polycrystal.orientations import quaternions

q1 = np.array([
    [ 0.6639988,  0.45485282,  0.83234881],
    [-0.20387788, -0.44343047, -0.28927565],
    [ 0.55001834, -0.70874605, -0.24522852],
    [ 0.46370166,  0.30685072, -0.40420048]
]).T
q2 = np.array([
    [ 0.40840558,  0.72671763,  0.36315318],
    [-0.30968977, -0.04341781, -0.40514601],
    [ 0.81406359, -0.34699685, -0.69307048],
    [-0.2731256,   0.591261,   -0.47289511]
]).T

qp = quaternions.multiply(q1, q2)
print(qp.T)

The result is:

[[-0.11305919 -0.11606473 -0.17603384]
 [-0.81660523 -0.65457512 -0.60644623]
 [ 0.56587988 -0.42403039 -0.63896874]
 [ 0.01238877  0.61502734 -0.43926532]]

Notice that the results are the negative of hexrd version. That is because q and -q are equivalent quaternions that both refer to the same rotation. In the original ODFPF, we chose to return the one with nonegative first value, for some reason which I don't remember exactly, but both values are correct.

This should definitely be fixed. Maybe using the scipy version is simplest. I believe the ordering of the quaternion components is different though, so you have to pay attention to that.