Efficient way to apply mirror effect on quaternion rotation?
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.
- If the scaling is negative (mirrored) along that 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
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)) thenq_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)) thenq_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)) thenq_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)) thenq_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)) thenq_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)) thenq_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)) thenq_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)) thenq_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)) thenq_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;
}