Look-at quaternion using up vector

Question

I have a camera (in a custom 3D engine) that accepts a quaternion for the rotation transform. I have two 3D points representing a camera and an object to look at. I want to calculate the quaternion that looks from the camera to the object, while respecting the world up axis.

This question asks for the same thing without the "up" vector. All three answers result in the camera pointing in the correct direction, but rolling (as in yaw/pitch/roll; imagine leaning your head onto your ear while looking at something).

I can calculate an orthonormal basis of vectors that match the desired coordinate system by:

lookAt = normalize(target - camera)
sideaxis = cross(lookAt, worldUp)
rotatedup = cross(sideaxis, lookAt)

How can I create a quaternion from those three vectors? This question asks for the same thing...but unfortunately the only and accepted answer says ~"let's assume you don't care about roll", and then goes about ignoring the up axis. I do care about roll. I don't want to ignore the up axis.

Answer 1

A previous answer has given a valid solution using angles. This answer will present an alternative method.

The orthonormal basis vectors, renaming them F = lookAt, R = sideaxis, U = rotatedup, directly form the columns of the 3x3 rotation matrix which is equivalent to your desired quaternion:

Multiplication with a vector is equivalent to using said vector's components as the coordinates in the camera's basis.

A 3x3 rotation matrix can be converted into a quaternion without conversion to angles / use of costly trigonometric functions. Below is a numerically stable C++ snippet which does this, returning a normalized quaternion:

inline void CalculateRotation( Quaternion& q ) const {
  float trace = a[0][0] + a[1][1] + a[2][2];
  if( trace > 0 ) {
    float s = 0.5f / sqrtf(trace + 1.0f);
    q.w = 0.25f / s;
    q.x = ( a[2][1] - a[1][2] ) * s;
    q.y = ( a[0][2] - a[2][0] ) * s;
    q.z = ( a[1][0] - a[0][1] ) * s;
  } else {
    if ( a[0][0] > a[1][1] && a[0][0] > a[2][2] ) {
      float s = 2.0f * sqrtf( 1.0f + a[0][0] - a[1][1] - a[2][2]);
      q.w = (a[2][1] - a[1][2] ) / s;
      q.x = 0.25f * s;
      q.y = (a[0][1] + a[1][0] ) / s;
      q.z = (a[0][2] + a[2][0] ) / s;
    } else if (a[1][1] > a[2][2]) {
      float s = 2.0f * sqrtf( 1.0f + a[1][1] - a[0][0] - a[2][2]);
      q.w = (a[0][2] - a[2][0] ) / s;
      q.x = (a[0][1] + a[1][0] ) / s;
      q.y = 0.25f * s;
      q.z = (a[1][2] + a[2][1] ) / s;
    } else {
      float s = 2.0f * sqrtf( 1.0f + a[2][2] - a[0][0] - a[1][1] );
      q.w = (a[1][0] - a[0][1] ) / s;
      q.x = (a[0][2] + a[2][0] ) / s;
      q.y = (a[1][2] + a[2][1] ) / s;
      q.z = 0.25f * s;
    }
  }
}

Source: http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion

Converting this to suit your situation is of course just a matter of swapping the matrix elements with the corresponding vector components:

// your code from before
F = normalize(target - camera);   // lookAt
R = normalize(cross(F, worldUp)); // sideaxis
U = cross(R, F);                  // rotatedup

// note that R needed to be re-normalized
// since F and worldUp are not necessary perpendicular
// so must remove the sin(angle) factor of the cross-product
// same not true for U because dot(R, F) = 0

// adapted source
Quaternion q;
double trace = R.x + U.y + F.z;
if (trace > 0.0) {
  double s = 0.5 / sqrt(trace + 1.0);
  q.w = 0.25 / s;
  q.x = (U.z - F.y) * s;
  q.y = (F.x - R.z) * s;
  q.z = (R.y - U.x) * s;
} else {
  if (R.x > U.y && R.x > F.z) {
    double s = 2.0 * sqrt(1.0 + R.x - U.y - F.z);
    q.w = (U.z - F.y) / s;
    q.x = 0.25 * s;
    q.y = (U.x + R.y) / s;
    q.z = (F.x + R.z) / s;
  } else if (U.y > F.z) {
    double s = 2.0 * sqrt(1.0 + U.y - R.x - F.z);
    q.w = (F.x - R.z) / s;
    q.x = (U.x + R.y) / s;
    q.y = 0.25 * s;
    q.z = (F.y + U.z) / s;
  } else {
    double s = 2.0 * sqrt(1.0 + F.z - R.x - U.y);
    q.w = (R.y - U.x) / s;
    q.x = (F.x + R.z) / s;
    q.y = (F.y + U.z) / s;
    q.z = 0.25 * s;
  }
}

(And needless to say swap y and z if you're using OpenGL.)

Answer 2

Assume you initially have three ortonormal vectors: worldUp, worldFront and worldSide, and lets use your equations for lookAt, sideAxis and rotatedUp. The worldSide vector will not be necessary to achieve the result.

Break the operation in two. First, rotate around worldUp. Then rotate around sideAxis, which will now actually be parallel to the rotated worldSide.

Axis1 = worldUp
Angle1 = (see below)

Axis2 = cross(lookAt, worldUp) = sideAxis
Angle2 = (see below)

Each of these rotations correspond to a quaternion using:

Q = cos(Angle/2) + i * Axis_x * sin(Angle/2) + j * Axis_y * sin(Angle/2) + k * Axis_z * sin(Angle/2)

Multiply both Q1 and Q2 and you get the desired quaternion.

Details for the angles:

Let P(worldUp) be the projection matrix on the worldUp direction, i.e., P(worldUp).v = cos(worldUp,v).worldUp or using kets and bras, P(worldUp) = |worldUp >< worldUp|. Let I be the identity matrix.

1. Project lookAt in the plane perpendicular to worldUp and normalize it.

   tmp1 = (I - P(worldUp)).lookAt
   n1 = normalize(tmp1)

2. Angle1 = arccos(dot(worldFront,n1))

3. Angle2 = arccos(dot(lookAt,n1))

EDIT1:

Notice that there is no need to compute transcendental functions. Since the dot product of a pair of normalized vectors is the cosine of an angle and assuming that cos(t) = x, we have the trigonometric identities:

• cos(t/2) = sqrt((1 + x)/2)
• sin(t/2) = sqrt((1 - x)/2)

Answer 3

A previous answer has given a valid solution using angles. This answer will present an alternative method.

The orthonormal basis vectors, renaming them F = lookAt, R = sideaxis, U = rotatedup, directly form the columns of the 3x3 rotation matrix which is equivalent to your desired quaternion:

Multiplication with a vector is equivalent to using said vector's components as the coordinates in the camera's basis.

A 3x3 rotation matrix can be converted into a quaternion without conversion to angles / use of costly trigonometric functions. Below is a numerically stable C++ snippet which does this, returning a normalized quaternion:

inline void CalculateRotation( Quaternion& q ) const {
  float trace = a[0][0] + a[1][1] + a[2][2];
  if( trace > 0 ) {
    float s = 0.5f / sqrtf(trace + 1.0f);
    q.w = 0.25f / s;
    q.x = ( a[2][1] - a[1][2] ) * s;
    q.y = ( a[0][2] - a[2][0] ) * s;
    q.z = ( a[1][0] - a[0][1] ) * s;
  } else {
    if ( a[0][0] > a[1][1] && a[0][0] > a[2][2] ) {
      float s = 2.0f * sqrtf( 1.0f + a[0][0] - a[1][1] - a[2][2]);
      q.w = (a[2][1] - a[1][2] ) / s;
      q.x = 0.25f * s;
      q.y = (a[0][1] + a[1][0] ) / s;
      q.z = (a[0][2] + a[2][0] ) / s;
    } else if (a[1][1] > a[2][2]) {
      float s = 2.0f * sqrtf( 1.0f + a[1][1] - a[0][0] - a[2][2]);
      q.w = (a[0][2] - a[2][0] ) / s;
      q.x = (a[0][1] + a[1][0] ) / s;
      q.y = 0.25f * s;
      q.z = (a[1][2] + a[2][1] ) / s;
    } else {
      float s = 2.0f * sqrtf( 1.0f + a[2][2] - a[0][0] - a[1][1] );
      q.w = (a[1][0] - a[0][1] ) / s;
      q.x = (a[0][2] + a[2][0] ) / s;
      q.y = (a[1][2] + a[2][1] ) / s;
      q.z = 0.25f * s;
    }
  }
}

Source: http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion

Converting this to suit your situation is of course just a matter of swapping the matrix elements with the corresponding vector components:

// your code from before
F = normalize(target - camera);   // lookAt
R = normalize(cross(F, worldUp)); // sideaxis
U = cross(R, F);                  // rotatedup

// note that R needed to be re-normalized
// since F and worldUp are not necessary perpendicular
// so must remove the sin(angle) factor of the cross-product
// same not true for U because dot(R, F) = 0

// adapted source
Quaternion q;
double trace = R.x + U.y + F.z;
if (trace > 0.0) {
  double s = 0.5 / sqrt(trace + 1.0);
  q.w = 0.25 / s;
  q.x = (U.z - F.y) * s;
  q.y = (F.x - R.z) * s;
  q.z = (R.y - U.x) * s;
} else {
  if (R.x > U.y && R.x > F.z) {
    double s = 2.0 * sqrt(1.0 + R.x - U.y - F.z);
    q.w = (U.z - F.y) / s;
    q.x = 0.25 * s;
    q.y = (U.x + R.y) / s;
    q.z = (F.x + R.z) / s;
  } else if (U.y > F.z) {
    double s = 2.0 * sqrt(1.0 + U.y - R.x - F.z);
    q.w = (F.x - R.z) / s;
    q.x = (U.x + R.y) / s;
    q.y = 0.25 * s;
    q.z = (F.y + U.z) / s;
  } else {
    double s = 2.0 * sqrt(1.0 + F.z - R.x - U.y);
    q.w = (R.y - U.x) / s;
    q.x = (F.x + R.z) / s;
    q.y = (F.y + U.z) / s;
    q.z = 0.25 * s;
  }
}

(And needless to say swap y and z if you're using OpenGL.)