Spectral Adjustment

In my opinion, spectral adjustment is the hardest part of decomposing an affine matrix. We need to take the non-unique result of our spectral decomposition and turn it into something that's geometrically meaningful. In Shoemake's reference code, this method is called snuggle

Spectral composition leaves us with the following decomposition: \(S = UKU^{T}\). Changing the labels & reversing teh directions of along which \(S\) is diagonalized will change the values of \(U\) and \(K\), but not the product \(UKU^{T}\). Out of all the valid variations of \(U\), we need to pick the one that has the smallest rotation angle.

To do this, consider interpolating between two stretch matrices, \(S_{1}\) and \(S_{2}\). We need to choose their diagonalizing rotations, \(U_{1}\) and \(U_{2}\) so that they are as similar as possible. The transform that goes from \(U_{1}\) to \(U_{2}\) is \(U_{12} = U_{1}^{T}U_{2}\). We need to minimize the angle of rotation performed by \(U_{12}\). We know what \(U_{1}\) is, how about \(U_{2}\)? \(U_{2}\) can be one of 24 matrices. These matrices are:

6 * 8 / 2 = 24, \(U_{2}\) can be one of 24 matrices. Knowing all possible values of \(U_{2}\), we need to compute \(U_{12}\) for each potential combination. To find the one with the smallest angle, convert \(U_{12}\) to a quaternion and pick the one with the largest w component. The eigen values in \(K\) should be re-arranged to match the axis permutation of the chosen \(U_{2}\). Multiply the eigen vectors by \(U_{2}^{T}\) to find the new eigen vectors, \(U = U_{1} U_{2}^{T}\)

The tricky part of finding all the possible values for \(U_{2}\) is figuring out which of the 48 combinations are achiavable by a rotation. The interactive sample below helps visualize this. The basis on the left is identiy (think of it as a potential value for \(U_{1}\)). The middle basis shows all axis permutations. The matrix on the right is all the sign combinations of the current permutation. Use the top slider to change permutation, the bottom slider to change sign combination. Look for matrices on the right that are rotations of the matrix on the left.

Canvas support required WebGL support required

Below is a table with all axis permutations and the sign combinations that can be achieved by a rotation.

Permutation Valid Axis 1 Valid Axis 2 Valid Axis 3 Valid Axis 4
X, Y, Z +X, +Y, +Z -X, -Y, -Z +X, -Y, -Z -X, +Y, -Z
X, Z, Y -X, +Z, +Y -X, -Z, -Y +X, +Z, -Y +X, -Z, +Y
Y, X, Z -Y, +X, +Z -Y, -X, -Z +Y, +X, -Z +Y, -X, +Z
Y, Z, X +Y, +Z, +X -Y, -Z, -X +Y, -Z, -X -Y, +Z, -X
Z, X, Y +Z, +X, +Y -Z, -X, -Y +Z, -X, -Y -Z, +X, -Y
Z, Y X -Z, +Y, +X -Z, -Y, -X +Z, +Y, -X +Z, -Y, +X 
struct SpectralAdjustmentResult {
    Matrix U; // Each basis vector is an eigenvector
    Matrix K; // Contains eigenvalues on main diagonal
    Matrix Ut; // Transpose of U
}

SpectralAdjustmentResult SpectralDecompositonAdjustment(SpectralDecompositionResult input) {
  Matrix U1 = input.U;
  Matrix U1t = input.Ut;

  Matrix[24][9] m_permutations = [
    // Permutation 0: x, y, z
    [ 1, 0, 0,  0, 1, 0,  0, 0, 1 ],
    [-1,-0,-0, -0,-1,-0,  0, 0, 1 ],
    [ 1, 0, 0, -0,-1,-0, -0,-0,-1 ],
    [-1,-0,-0,  0, 1, 0, -0,-0,-1 ],
    // Permutation 1: x, z, y
    [-1,-0,-0,  0, 0, 1,  0, 1, 0 ],
    [-1,-0,-0, -0,-0,-1, -0,-1,-0 ],
    [ 1, 0, 0,  0, 0, 1, -0,-1,-0 ],
    [ 1, 0, 0, -0,-0,-1,  0, 1, 0 ],
    // Permutation 2: y, x, z
    [-0,-1,-0,  1, 0, 0,  0, 0, 1 ],
    [-0,-1,-0, -1,-0,-0, -0,-0,-1 ],
    [ 0, 1, 0,  1, 0, 0, -0,-0,-1 ],
    [ 0, 1, 0, -1,-0,-0,  0, 0, 1 ],
    // Permutation 3: y, z, x
    [ 0, 1, 0,  0, 0, 1,  1, 0, 0 ],
    [-0,-1,-0, -0,-0,-1,  1, 0, 0 ],
    [ 0, 1, 0, -0,-0,-1, -1,-0,-0 ],
    [-0,-1,-0,  0, 0, 1, -1,-0,-0 ],
     // Permutation 4: z, x, y
    [ 0, 0, 1,  1, 0, 0,  0, 1, 0 ],
    [-0,-0,-1, -1,-0,-0,  0, 1, 0 ],
    [ 0, 0, 1, -1,-0,-0, -0,-1,-0 ],
    [-0,-0,-1,  1, 0, 0, -0,-1,-0 ],
    // Permutation 5: z, y, x
    [-0,-0,-1,  0, 1, 0,  1, 0, 0 ],
    [-0,-0,-1, -0,-1,-0, -1,-0,-0 ],
    [ 0, 0, 1,  0, 1, 0, -1,-0,-0 ],
    [ 0, 0, 1, -0,-1,-0,  1, 0, 0],
  ];

  float x = input.K[0];
  float y = input.K[5];
  float z = input.K[10];

  Vector[6][3] eigen_value_permutations = [
    [ x,  y,  z], // Permutation 0
    [ x,  z,  y], // Permutation 1
    [ y,  x,  z], // Permutation 2
    [ y,  z,  x], // Permutation 3
    [ z,  x,  y], // Permutation 4
    [ z,  y,  x], // Permutation 5
  ];

  int saved_index = -1
  int saved_value = -1

  // The rotation taking U1 into U2 is U1t * U2
  for (int i = 0; i < 24; ++i) {
    Matrix U2 = m_permutations[i];
    Matrix U12 = U1t * U2;

    Quaternion QU12 = ToQuaternion(U12);
    
    // Optimize for largest w, which is smallest angle of rotation
    if (saved_index == -1 || QU12.w > saved_value) {
      saved_value = QU12.w
      saved_index = i
    }
  }

  Matrix U2t = Transpose(m_permutations[saved_index])
  var index = saved_index/4; // Integer division (floor this)

  SpectralAdjustmentResult result;
  result.U = Mul3(U1, U2t);
  result.K = Matrix(
    eigen_value_permutations[index][0], 0, 0, 0,
    0, eigen_value_permutations[index][1], 0, 0,
    0, 0, eigen_value_permutations[index][2], 0,
    0, 0, 0, 1
  );
  result.Ut = Transpose(result.U)

  return result;
}