decomposeMatrix static method

void decomposeMatrix(List<double> transformMatrix, MatrixDecompositionContext ctx)

Implementation

static void decomposeMatrix(
    List<double> transformMatrix, MatrixDecompositionContext ctx) {
  // output values
  final perspective = ctx.perspective;
  final quaternion = ctx.quaternion;
  final scale = ctx.scale;
  final skew = ctx.skew;
  final translation = ctx.translation;
  final rotationDegrees = ctx.rotationDegrees;

  // create normalized, 2d array matrix
  // and normalized 1d array perspectiveMatrix with redefined 4th column
  if (_isZero(transformMatrix[15])) {
    return;
  }
  var matrix = <List<double>>[
    [0, 0, 0, 0],
    [0, 0, 0, 0],
    [0, 0, 0, 0],
    [0, 0, 0, 0]
  ];
  var perspectiveMatrix = List<double>.filled(16, 0);
  for (var i = 0; i < 4; i++) {
    for (var j = 0; j < 4; j++) {
      var value = transformMatrix[(i * 4) + j] / transformMatrix[15];
      matrix[i][j] = value;
      perspectiveMatrix[(i * 4) + j] = j == 3 ? 0 : value;
    }
  }
  perspectiveMatrix[15] = 1;

  // test for singularity of upper 3x3 part of the perspective matrix
  if (_isZero(determinant(perspectiveMatrix))) {
    return;
  }

  // isolate perspective
  if (!_isZero(matrix[0][3]) ||
      !_isZero(matrix[1][3]) ||
      !_isZero(matrix[2][3])) {
    // rightHandSide is the right hand side of the equation.
    // rightHandSide is a vector, or point in 3d space relative to the origin.
    var rightHandSide = <double>[
      matrix[0][3],
      matrix[1][3],
      matrix[2][3],
      matrix[3][3]
    ];

    // Solve the equation by inverting perspectiveMatrix and multiplying
    // rightHandSide by the inverse.
    var inversePerspectiveMatrix = inverse(perspectiveMatrix);
    var transposedInversePerspectiveMatrix =
        transpose(inversePerspectiveMatrix);
    multiplyVectorByMatrix(
        rightHandSide, transposedInversePerspectiveMatrix, perspective);
  } else {
    // no perspective
    perspective[0] = perspective[1] = perspective[2] = 0;
    perspective[3] = 1;
  }

  // translation is simple
  for (var i = 0; i < 3; i++) {
    translation[i] = matrix[3][i];
  }

  // Now get scale and shear.
  // 'row' is a 3 element array of 3 component vectors
  var row = <List<double>>[
    [0, 0, 0],
    [0, 0, 0],
    [0, 0, 0]
  ];
  for (var i = 0; i < 3; i++) {
    row[i][0] = matrix[i][0];
    row[i][1] = matrix[i][1];
    row[i][2] = matrix[i][2];
  }

  // Compute X scale factor and normalize first row.
  scale[0] = v3Length(row[0]);
  row[0] = v3Normalize(row[0], scale[0]);

  // Compute XY shear factor and make 2nd row orthogonal to 1st.
  skew[0] = v3Dot(row[0], row[1]);
  row[1] = v3Combine(row[1], row[0], 1.0, -skew[0]);

  // Compute XY shear factor and make 2nd row orthogonal to 1st.
  skew[0] = v3Dot(row[0], row[1]);
  row[1] = v3Combine(row[1], row[0], 1.0, -skew[0]);

  // Now, compute Y scale and normalize 2nd row.
  scale[1] = v3Length(row[1]);
  row[1] = v3Normalize(row[1], scale[1]);
  skew[0] /= scale[1];

  // Compute XZ and YZ shears, orthogonalize 3rd row
  skew[1] = v3Dot(row[0], row[2]);
  row[2] = v3Combine(row[2], row[0], 1.0, -skew[1]);
  skew[2] = v3Dot(row[1], row[2]);
  row[2] = v3Combine(row[2], row[1], 1.0, -skew[2]);

  // Next, get Z scale and normalize 3rd row.
  scale[2] = v3Length(row[2]);
  row[2] = v3Normalize(row[2], scale[2]);
  skew[1] /= scale[2];
  skew[2] /= scale[2];

  // At this point, the matrix (in rows) is orthonormal.
  // Check for a coordinate system flip.  If the determinant
  // is -1, then negate the matrix and the scaling factors.
  var pdum3 = v3Cross(row[1], row[2]);
  if (v3Dot(row[0], pdum3) < 0) {
    for (var i = 0; i < 3; i++) {
      scale[i] *= -1;
      row[i][0] *= -1;
      row[i][1] *= -1;
      row[i][2] *= -1;
    }
  }

  // Now, get the rotations out
  quaternion[0] = 0.5 * sqrt(max(1 + row[0][0] - row[1][1] - row[2][2], 0));
  quaternion[1] = 0.5 * sqrt(max(1 - row[0][0] + row[1][1] - row[2][2], 0));
  quaternion[2] = 0.5 * sqrt(max(1 - row[0][0] - row[1][1] + row[2][2], 0));
  quaternion[3] = 0.5 * sqrt(max(1 + row[0][0] + row[1][1] + row[2][2], 0));

  if (row[2][1] > row[1][2]) {
    quaternion[0] = -quaternion[0];
  }
  if (row[0][2] > row[2][0]) {
    quaternion[1] = -quaternion[1];
  }
  if (row[1][0] > row[0][1]) {
    quaternion[2] = -quaternion[2];
  }

  // correct for occasional, weird Euler synonyms for 2d rotation

  if (quaternion[0] < 0.001 &&
      quaternion[0] >= 0 &&
      quaternion[1] < 0.001 &&
      quaternion[1] >= 0) {
    // this is a 2d rotation on the z-axis
    rotationDegrees[0] = rotationDegrees[1] = 0;
    rotationDegrees[2] =
        roundTo3Places(atan2(row[0][1], row[0][0]) * 180 / pi);
  } else {
    quaternionToDegreesXYZ(quaternion, rotationDegrees);
  }
}