Efficient way to apply mirror effect on quaternion rotation?

Question

Quaternions represent rotations - they don't include information about scaling or mirroring. However it is still possible to mirror the effect of a rotation.

Consider a mirroring on the x-y-plane (we can also call it a mirroring along the z-axis). A rotation around the x-axis mirrored on the x-y-plane would be negated. Likewise with a rotation around the y axis. However, a rotation around the z-axis would be left unchanged.

Another example: 90º rotation around axis (1,1,1) mirrored in the x-y plane would give -90º rotation around (1,1,-1). To aid the intuition, if you can visualize a depiction of the axis and a circular arrow indicating the rotation, then mirroring that visualization indicates what the new rotation should be.

I have found a way to calculate this mirroring of the rotation, like this:

• Get the angle-axis representation of the quaternion.
• For each of the axes x, y, and z.
  • If the scaling is negative (mirrored) along that axis:
    • Negate both angle and axis.
• Get the updated quaternion from the modified angle and axis.

This only supports mirroring along the primary axes, x, y, and z, since that's all I need. It works for arbitrary rotations though.

However, the conversions from quaternion to angle-axis and back from angle-axis to quaternion are expensive. I'm wondering if there's a way to do the conversion directly on the quaternion itself, but my comprehension of quaternion math is not sufficient to get anywhere myself.

(Posted on StackOverflow rather than math-related forums due to the importance of a computationally efficient method.)

Answer 1

I just spent quite some time on figuring out a clear answer to this question, so I am posting it here for the record.

Introduction

As was noted in other answers, a mirror effect cannot be represented as a rotation. However, given a rotation R1to2 from a coordinate frame C1 to a coordinate frame C2, we may be interested in efficiently computing the equivalent rotation when applying the same mirror effect to C1 and C2 (e.g. I was facing the problem of converting an input quaternion, given in a left-handed coordinate frame, into the quaternion representing the same rotation but in a right-handed coordinate frame).

In terms of rotation matrices, this can be thought of as follows:

R_mirroredC1_to_mirroredC2 = M_mirrorC2 * R_C1_to_C2 * M_mirrorC1

Here, both R_C1_to_C2 and R_mirroredC1_to_mirroredC2 represent valid rotations, so when dealing with quaternions, how do you efficiently compute q_mirroredC1_to_mirroredC2 from q_C1_to_C2?

Solution

The following assumes that q_C1_to_C2=[w,x,y,z]:

• if C1 and C2 are mirrored along the X-axis (i.e. M_mirrorC1=M_mirrorC2=diag_3x3(-1,1,1)) then q_mirroredC1_to_mirroredC2=[w,x,-y,-z]
• if C1 and C2 are mirrored along the Y-axis (i.e. M_mirrorC1=M_mirrorC2=diag_3x3(1,-1,1)) then q_mirroredC1_to_mirroredC2=[w,-x,y,-z]
• if C1 and C2 are mirrored along the Z-axis (i.e. M_mirrorC1=M_mirrorC2=diag_3x3(1,1,-1)) then q_mirroredC1_to_mirroredC2=[w,-x,-y,z]

When considering different mirrored axes for the C1 and C2, we have the following:

• if C1 is mirrored along the X-axis and C2 along the Y-axis (i.e. M_mirrorC1=diag_3x3(-1,1,1) & M_mirrorC2=diag_3x3(1,-1,1)) then q_mirroredC1_to_mirroredC2=[z,y,x,w]

• if C1 is mirrored along the X-axis and C2 along the Z-axis (i.e. M_mirrorC1=diag_3x3(-1,1,1) & M_mirrorC2=diag_3x3(1,1,-1)) then q_mirroredC1_to_mirroredC2=[-y,z,-w,x]

• if C1 is mirrored along the Y-axis and C2 along the X-axis (i.e. M_mirrorC1=diag_3x3(1,-1,1) & M_mirrorC2=diag_3x3(-1,1,1)) then q_mirroredC1_to_mirroredC2=[z,-y,-x,w]

• if C1 is mirrored along the Y-axis and C2 along the Z-axis (i.e. M_mirrorC1=diag_3x3(1,-1,1) & M_mirrorC2=diag_3x3(1,1,-1)) then q_mirroredC1_to_mirroredC2=[x,w,z,y]

• if C1 is mirrored along the Z-axis and C2 along the X-axis (i.e. M_mirrorC1=diag_3x3(1,1,-1) & M_mirrorC2=diag_3x3(-1,1,1)) then q_mirroredC1_to_mirroredC2=[y,z,w,x]

• if C1 is mirrored along the Z-axis and C2 along the Y-axis (i.e. M_mirrorC1=diag_3x3(1,1,-1) & M_mirrorC2=diag_3x3(1,-1,1)) then q_mirroredC1_to_mirroredC2=[x,w,-z,-y]

Test program

Here is a small c++ program based on OpenCV to test all this:

#include <opencv2/opencv.hpp>
#define CST_PI 3.1415926535897932384626433832795

// Random rotation matrix uniformly sampled from SO3 (see "Fast random rotation matrices" by J.Arvo)
cv::Matx<double,3,3> get_random_rotmat()
{
    double theta1 = 2*CST_PI*cv::randu<double>();
    double theta2 = 2*CST_PI*cv::randu<double>();
    double x3 = cv::randu<double>();
    cv::Matx<double,3,3> R(std::cos(theta1),std::sin(theta1),0,-std::sin(theta1),std::cos(theta1),0,0,0,1);
    cv::Matx<double,3,1> v(std::cos(theta2)*std::sqrt(x3),std::sin(theta2)*std::sqrt(x3),std::sqrt(1-x3));
    return -1*(cv::Matx<double,3,3>::eye()-2*v*v.t())*R;
}

cv::Matx<double,4,1> rotmat2quatwxyz(const cv::Matx<double,3,3> &R)
{
    // Implementation from Ceres 1.10
    const double trace = R(0,0) + R(1,1) + R(2,2);
    cv::Matx<double,4,1> quat_wxyz;
    if (trace >= 0.0) {
        double t = sqrt(trace + 1.0);
        quat_wxyz(0) = 0.5 * t;
        t = 0.5 / t;
        quat_wxyz(1) = (R(2,1) - R(1,2)) * t;
        quat_wxyz(2) = (R(0,2) - R(2,0)) * t;
        quat_wxyz(3) = (R(1,0) - R(0,1)) * t;
    } else {
        int i = 0;
        if (R(1, 1) > R(0, 0))
            i = 1;
        if (R(2, 2) > R(i, i))
            i = 2;

        const int j = (i + 1) % 3;
        const int k = (j + 1) % 3;
        double t = sqrt(R(i, i) - R(j, j) - R(k, k) + 1.0);
        quat_wxyz(i + 1) = 0.5 * t;
        t = 0.5 / t;
        quat_wxyz(0)     = (R(k,j) - R(j,k)) * t;
        quat_wxyz(j + 1) = (R(j,i) + R(i,j)) * t;
        quat_wxyz(k + 1) = (R(k,i) + R(i,k)) * t;
    }
    // Check that the w element is positive
    if(quat_wxyz(0)<0)
        quat_wxyz *= -1;    // quat and -quat represent the same rotation, but to make quaternion comparison easier, we always use the one with positive w
    return quat_wxyz;
}

cv::Matx<double,4,1> apply_quaternion_trick(const unsigned int item_permuts[4], const int sign_flips[4], const cv::Matx<double,4,1>& quat_wxyz)
{
    // Flip the sign of the x and z components
    cv::Matx<double,4,1> quat_flipped(sign_flips[0]*quat_wxyz(item_permuts[0]),sign_flips[1]*quat_wxyz(item_permuts[1]),sign_flips[2]*quat_wxyz(item_permuts[2]),sign_flips[3]*quat_wxyz(item_permuts[3]));
    // Check that the w element is positive
    if(quat_flipped(0)<0)
        quat_flipped *= -1; // quat and -quat represent the same rotation, but to make quaternion comparison easier, we always use the one with positive w
    return quat_flipped;
}

void detect_quaternion_trick(const cv::Matx<double,4,1> &quat_regular, const cv::Matx<double,4,1> &quat_flipped, unsigned int item_permuts[4], int sign_flips[4])
{
    if(abs(quat_regular(0))==abs(quat_flipped(0))) {
        item_permuts[0]=0;
        sign_flips[0] = (quat_regular(0)/quat_flipped(0)>0 ? 1 : -1);
    }
    else if(abs(quat_regular(0))==abs(quat_flipped(1))) {
        item_permuts[1]=0;
        sign_flips[1] = (quat_regular(0)/quat_flipped(1)>0 ? 1 : -1);
    }
    else if(abs(quat_regular(0))==abs(quat_flipped(2))) {
        item_permuts[2]=0;
        sign_flips[2] = (quat_regular(0)/quat_flipped(2)>0 ? 1 : -1);
    }
    else if(abs(quat_regular(0))==abs(quat_flipped(3))) {
        item_permuts[3]=0;
        sign_flips[3] = (quat_regular(0)/quat_flipped(3)>0 ? 1 : -1);
    }
    if(abs(quat_regular(1))==abs(quat_flipped(0))) {
        item_permuts[0]=1;
        sign_flips[0] = (quat_regular(1)/quat_flipped(0)>0 ? 1 : -1);
    }
    else if(abs(quat_regular(1))==abs(quat_flipped(1))) {
        item_permuts[1]=1;
        sign_flips[1] = (quat_regular(1)/quat_flipped(1)>0 ? 1 : -1);
    }
    else if(abs(quat_regular(1))==abs(quat_flipped(2))) {
        item_permuts[2]=1;
        sign_flips[2] = (quat_regular(1)/quat_flipped(2)>0 ? 1 : -1);
    }
    else if(abs(quat_regular(1))==abs(quat_flipped(3))) {
        item_permuts[3]=1;
        sign_flips[3] = (quat_regular(1)/quat_flipped(3)>0 ? 1 : -1);
    }
    if(abs(quat_regular(2))==abs(quat_flipped(0))) {
        item_permuts[0]=2;
        sign_flips[0] = (quat_regular(2)/quat_flipped(0)>0 ? 1 : -1);
    }
    else if(abs(quat_regular(2))==abs(quat_flipped(1))) {
        item_permuts[1]=2;
        sign_flips[1] = (quat_regular(2)/quat_flipped(1)>0 ? 1 : -1);
    }
    else if(abs(quat_regular(2))==abs(quat_flipped(2))) {
        item_permuts[2]=2;
        sign_flips[2] = (quat_regular(2)/quat_flipped(2)>0 ? 1 : -1);
    }
    else if(abs(quat_regular(2))==abs(quat_flipped(3))) {
        item_permuts[3]=2;
        sign_flips[3] = (quat_regular(2)/quat_flipped(3)>0 ? 1 : -1);
    }
    if(abs(quat_regular(3))==abs(quat_flipped(0))) {
        item_permuts[0]=3;
        sign_flips[0] = (quat_regular(3)/quat_flipped(0)>0 ? 1 : -1);
    }
    else if(abs(quat_regular(3))==abs(quat_flipped(1))) {
        item_permuts[1]=3;
        sign_flips[1] = (quat_regular(3)/quat_flipped(1)>0 ? 1 : -1);
    }
    else if(abs(quat_regular(3))==abs(quat_flipped(2))) {
        item_permuts[2]=3;
        sign_flips[2] = (quat_regular(3)/quat_flipped(2)>0 ? 1 : -1);
    }
    else if(abs(quat_regular(3))==abs(quat_flipped(3))) {
        item_permuts[3]=3;
        sign_flips[3] = (quat_regular(3)/quat_flipped(3)>0 ? 1 : -1);
    }
}

int main(int argc, char **argv)
{
    cv::Matx<double,3,3> M_xflip(-1,0,0,0,1,0,0,0,1);
    cv::Matx<double,3,3> M_yflip(1,0,0,0,-1,0,0,0,1);
    cv::Matx<double,3,3> M_zflip(1,0,0,0,1,0,0,0,-1);

    // Let the user choose the configuration
    char im,om;
    std::cout << "Enter the axis (x,y,z) along which input ref is flipped:" << std::endl;
    std::cin >> im;
    std::cout << "Enter the axis (x,y,z) along which output ref is flipped:" << std::endl;
    std::cin >> om;
    cv::Matx<double,3,3> M_iflip,M_oflip;
    if(im=='x') M_iflip=M_xflip;
    else if(im=='y') M_iflip=M_yflip;
    else if(im=='z') M_iflip=M_zflip;
    if(om=='x') M_oflip=M_xflip;
    else if(om=='y') M_oflip=M_yflip;
    else if(om=='z') M_oflip=M_zflip;

    // Generate random quaternions until we find one where no two elements are equal
    cv::Matx<double,3,3> R;
    cv::Matx<double,4,1> quat_regular,quat_flipped;
    do {
        R = get_random_rotmat();
        quat_regular = rotmat2quatwxyz(R);
    } while(quat_regular(0)==quat_regular(1) || quat_regular(0)==quat_regular(2) || quat_regular(0)==quat_regular(3) ||
            quat_regular(1)==quat_regular(2) || quat_regular(1)==quat_regular(3) ||
            quat_regular(2)==quat_regular(3));

    // Determine and display the appropriate quaternion trick
    quat_flipped = rotmat2quatwxyz(M_oflip*R*M_iflip);
    unsigned int item_permuts[4]={0,1,2,3};
    int sign_flips[4]={1,1,1,1};
    detect_quaternion_trick(quat_regular,quat_flipped,item_permuts,sign_flips);
    char str_quat[4]={'w','x','y','z'};
    std::cout << std::endl << "When iref is flipped along the " << im << "-axis and oref along the " << om << "-axis:" << std::endl;
    std::cout << "resulting_quat=[" << (sign_flips[0]>0?"":"-") << str_quat[item_permuts[0]] << ","
                                    << (sign_flips[1]>0?"":"-") << str_quat[item_permuts[1]] << ","
                                    << (sign_flips[2]>0?"":"-") << str_quat[item_permuts[2]] << ","
                                    << (sign_flips[3]>0?"":"-") << str_quat[item_permuts[3]] << "], where initial_quat=[w,x,y,z]" << std::endl;

    // Test this trick on several random rotation matrices
    unsigned int n_errors = 0, n_tests = 10000;
    std::cout << std::endl << "Performing " << n_tests << " tests on random rotation matrices:" << std::endl;
    for(unsigned int i=0; i<n_tests; ++i) {

        // Get a random rotation matrix and the corresponding quaternion
        cv::Matx<double,3,3> R = get_random_rotmat();
        cv::Matx<double,4,1> quat_regular = rotmat2quatwxyz(R);
        // Get the quaternion corresponding to the flipped coordinate frames, via the sign trick and via computation on rotation matrices
        cv::Matx<double,4,1> quat_tricked = apply_quaternion_trick(item_permuts,sign_flips,quat_regular);
        cv::Matx<double,4,1> quat_flipped = rotmat2quatwxyz(M_oflip*R*M_iflip);
        // Check that both results are identical
        if(cv::norm(quat_tricked-quat_flipped,cv::NORM_INF)>1e-6) {
            std::cout << "Error (idx=" << i << ")!"
                      << "\n   quat_regular=" << quat_regular.t()
                      << "\n   quat_tricked=" << quat_tricked.t()
                      << "\n   quat_flipped=" << quat_flipped.t() << std::endl;
            ++n_errors;
        }
    }
    std::cout << n_errors << " errors on " << n_tests << " tests." << std::endl;
    system("pause");
    return 0;
}

Answer 2

There is little bit easier and programmer oriented way to think about this. Assume that you want to reverse the z axis (i.e. flip z axis to -z) in your coordinate system. Now think of quaternion as orientation vector in terms of roll, pitch and yaw. When you flip z axis, notice that sign of roll and pitch is inverted but sign for yaw remains same.

Now you can find the net effect on quaternion using following code for converting Euler angles to quaternion (I'd put this code to Wikipedia):

static Quaterniond toQuaternion(double pitch, double roll, double yaw)
{
    Quaterniond q;
    double t0 = std::cos(yaw * 0.5f);
    double t1 = std::sin(yaw * 0.5f);
    double t2 = std::cos(roll * 0.5f);
    double t3 = std::sin(roll * 0.5f);
    double t4 = std::cos(pitch * 0.5f);
    double t5 = std::sin(pitch * 0.5f);

    q.w() = t0 * t2 * t4 + t1 * t3 * t5;
    q.x() = t0 * t3 * t4 - t1 * t2 * t5;
    q.y() = t0 * t2 * t5 + t1 * t3 * t4;
    q.z() = t1 * t2 * t4 - t0 * t3 * t5;
    return q;
}

Using basic trigonometry, sin(-x) = -sin(x) and cos(-x) = cos(x). Applyieng this to above code you can see that sign for t3 and t5 will flip. This will cause sign of x and y to flip.

So when you invert the z-axis,

Q'(w, x, y, z) = Q(w, -x, -y, z)

Similarly you can figure out any other combinations of axis reversal and find impact on quaternion.

PS: In case if anyone is wondering why anyone would ever need this... I needed above to transform quaternion coming from MavLink/Pixhawk system which controls drone. The source system uses NED coordinate system but usual 3D environments like Unreal uses NEU coordinate system which requires transforming z axis to -z to use the quaternion correctly.

Answer 3

I did some further analysis, and it appears the effect of a quaternion (w, x, y, z) can have it's effect mirrored like this:

• Mirror effect of rotation along x axis by flipping y and z elements of the quaternion.
• Mirror effect of rotation along y axis by flipping x and z elements of the quaternion.
• Mirror effect of rotation along z axis by flipping x and y elements of the quaternion.

The w element of the quaternion never needs to be touched.

Unfortunately I still don't understand quaternions well enough to be able to explain why this works, but I derived it from implementations of converting to and from axis-angle format, and after implementing this solution, it works just as well as my original one in all tests of it I have performed.

Answer 4

We can examine the set of all rotations and reflections in 3D this is called the Orthogonal group O(3). It can be though of as the set of orthogonal matrices with determinant +1 or -1. All rotations have determinant +1 and pure reflections have determinate -1. There is another member of O(3) the inversion in a point (x,y,z)->(-x,-y,-z) this has det -1 in 3D and we will come to this later. If we combine two transformations in the group you multiply their determinants. Hence two rotations combined give another rotation (+1 * +1 = +1), a rotation combined with a reflection give a reflection (+1 * -1 = -1) and two reflections combined give a rotation (-1 * -1 = +1).

We can restrict the O(3) to just those with determinant +1 to form the Special Orthogonal Group SO(3). This just contains the rotations.

Now the set of unit quaternions is the double cover of SO(3) that means that two unit quaternions correspond to each rotation. To be precise if a+b i+c j+d k is a unit quaternions then a-b i-c j-d k represents the same rotation, you can think of this as a rotation by ø around the vector (b,c,d) being the same as a rotation by -ø around the vector (-b,-c,-d).

Note that all the unit quaternions have determinant +1, so there is none which correspond to a pure reflection. This is why you cannot use quaternions to represent reflections.

What you might be able to do is use the inversion. Now a reflection followed by an inversion is a rotation. For example reflect in x=0 and invert, is the same as reflecting in the y=0 and reflecting in the z=0. This is the same as 180º rotation around the x-axis. You could do the same procedure for any reflection.

---

We can define a plane through the origin by using it normal vector n = (a,b,c). A reflection of a vector v(x,y,z) in that plane is given by

v - 2 (v . n ) / ( n . n) n

= (x,y,z) - 2 (a x+b y+c z) / (a^2+b^2+c^2) (a,b,c)

In particular the x-y plane has normal (0,0,1) so a reflection is

(x,y,z) - 2 z (0,0,1) = (x,y,-z)

Quaternions and spatial rotation has a nice formula for a quaternion from the axis angle formula.

p = cos(ø/2) + (x i + y j + z k) sin(ø/2)

This is a quaternion W + X i + Y j + Z k with W=cos(ø/2), X = x sin(ø/2), Y = y sin(ø/2), Z = z sin(ø/2)

Changing the direction of rotation will flip the sin of the half angle but leave the cos unchanged, giving

p' = cos(ø/2) - (x i + y j + z k) sin(ø/2)

Now if we consider reflecting the corresponding vector in x-y plane giving

q = cos(ø/2) + (x i + y j - z k) sin(ø/2)

we might want to change the direction of rotation giving

q' = cos(ø/2) + (- x i - y j + z k) sin(ø/2)

= W - X i - Y j + Z k

which I think corresponds to your answer.

We can generalise this to reflection in a general plane with unit length normal (a,b,c). Let d be the dot product (a,b,c).(x,y,z). The refection of (x,y,z) is

(x,y,z) - 2 d (a,b,c) = (x - 2 d a, y - 2 d b, z - 2 d c)

the rotation quaternion of this

q = cos(ø/2) - ((x - 2 d a) i + ((y - 2 d b) j + (z - 2 d c) k) sin(ø/2)

q = cos(ø/2) - (x i + y j + z k) sin(ø/2) + 2 d sin(ø/2) (a i + b j + c k)

= W - X i - Y j - Z k + 2 d (X,Y,Z).(a,b,c) (a i + b j + c k)

Answer 5

Note that mirroring is not a rotation, so generally you can't bake it into a quaternion (I might very well have misunderstood your question, though). The 3x3 component of the mirroring transformation matrix is

M = I-2(n*nT) 

where I is an identity 3x3 matrix, n is the mirror plane's normal represented as a 3x1 matrix, and nT is n as a 1x3 matrix (so n*nT is a 3x(1x1)x3=3x3 matrix).

Now, if the quaternion q you want to 'mirror' is the last transformation, the last transformation on the other side would be just M*q (again, this would be a general 3x3 matrix, not generally representable as a quaternion)

Answer 6

For anyone who gets here by a web-search and is looking for the math, then:

---

Reflection

To reflecting point 'p' through plane ax+by+cz=0, using quaternions:

n = 0+(a,b,c)
p = 0+(x,y,z)

where 'n' is a unit bivector (or pure quaternion if you prefer)

p' = npn

then p' is the reflect point.

If you compose with a second reflection 'm':

p' = mnpnm = (mn)p(mn)^*

is a rotation.

Rotations and reflections compose as expected.

---

Uniform scaling

Since scalar products commute and can be factor out then if we have either a rotation by unit quaternion 'Q' or a reflection by unit bivector 'b' (or any combination of) multiplying either by some non-zero scale value 's' results in a uniform scaling of s^2. And since (sqrt(s0)*sqrt(s1))^2 = s0*s1, these uniform scaling value compose as expected.

However this point is probably of no interest since in code we want to be able to assume unit magnitude values to reduce the runtime complexity.