ImathMatrixAlgo.h

Go to the documentation of this file.

//
// SPDX-License-Identifier: BSD-3-Clause
// Copyright Contributors to the OpenEXR Project.
//

//
//
// Functions operating on Matrix22, Matrix33, and Matrix44 types
//
// This file also defines a few predefined constant matrices.
//

#ifndef INCLUDED_IMATHMATRIXALGO_H
#define INCLUDED_IMATHMATRIXALGO_H

#include "ImathEuler.h"
#include "ImathExport.h"
#include "ImathMatrix.h"
#include "ImathNamespace.h"
#include "ImathQuat.h"
#include "ImathVec.h"
#include <math.h>

IMATH_INTERNAL_NAMESPACE_HEADER_ENTER

//------------------
// Identity matrices
//------------------

/// M22f identity matrix
IMATH_EXPORT_CONST M22f identity22f;
/// M33f identity matrix
IMATH_EXPORT_CONST M33f identity33f;
/// M44f identity matrix
IMATH_EXPORT_CONST M44f identity44f;
/// M22d identity matrix
IMATH_EXPORT_CONST M22d identity22d;
/// M33d identity matrix
IMATH_EXPORT_CONST M33d identity33d;
/// M44d identity matrix
IMATH_EXPORT_CONST M44d identity44d;

//----------------------------------------------------------------------
// Extract scale, shear, rotation, and translation values from a matrix:
//
// Notes:
//
// This implementation follows the technique described in the paper by
// Spencer W. Thomas in the Graphics Gems II article: "Decomposing a
// Matrix into Simple Transformations", p. 320.
//
// - Some of the functions below have an optional exc parameter
//   that determines the functions' behavior when the matrix'
//   scaling is very close to zero:
//
//   If exc is true, the functions throw a std::domain_error exception.
//
//   If exc is false:
//
//      extractScaling (m, s)            returns false, s is invalid
//      sansScaling (m)                  returns m
//      removeScaling (m)                returns false, m is unchanged
//      sansScalingAndShear (m)          returns m
//      removeScalingAndShear (m)        returns false, m is unchanged
//      extractAndRemoveScalingAndShear (m, s, h)
//                                       returns false, m is unchanged,
//                                                      (sh) are invalid
//      checkForZeroScaleInRow ()        returns false
//      extractSHRT (m, s, h, r, t)      returns false, (shrt) are invalid
//
// - Functions extractEuler(), extractEulerXYZ() and extractEulerZYX()
//   assume that the matrix does not include shear or non-uniform scaling,
//   but they do not examine the matrix to verify this assumption.
//   Matrices with shear or non-uniform scaling are likely to produce
//   meaningless results.  Therefore, you should use the
//   removeScalingAndShear() routine, if necessary, prior to calling
//   extractEuler...() .
//
// - All functions assume that the matrix does not include perspective
//   transformation(s), but they do not examine the matrix to verify
//   this assumption.  Matrices with perspective transformations are
//   likely to produce meaningless results.
//
//----------------------------------------------------------------------

//
// Declarations for 4x4 matrix.
//

/// Extract the scaling component of the given 4x4 matrix.
///
/// @param[in] mat The input matrix
/// @param[out] scl The extracted scale, i.e. the output value
/// @param[in] exc If true, throw an exception if the scaling in `mat` is very close to zero.
/// @return True if the scale could be extracted, false if the matrix is degenerate.
template <class T>
bool extractScaling (const Matrix44<T>& mat, Vec3<T>& scl, bool exc = true);

/// Return the given 4x4 matrix with scaling removed.
///
/// @param[in] mat The input matrix
/// @param[in] exc If true, throw an exception if the scaling in `mat`
template <class T>
Matrix44<T> sansScaling (const Matrix44<T>& mat, bool exc = true);

/// Remove scaling from the given 4x4 matrix in place.  Return true if the
/// scale could be successfully extracted, false if the matrix is
/// degenerate.
//
/// @param[in] mat The matrix to operate on
/// @param[in] exc If true, throw an exception if the scaling in `mat` is very close to zero.
/// @return True if the scale could be extracted, false if the matrix is degenerate.
template <class T> bool removeScaling (Matrix44<T>& mat, bool exc = true);

/// Extract the scaling and shear components of the given 4x4 matrix.
/// Return true if the scale could be successfully extracted, false if
/// the matrix is degenerate.
///
/// @param[in] mat The input matrix
/// @param[out] scl The extracted scale
/// @param[out] shr The extracted shear
/// @param[in] exc If true, throw an exception if the scaling in `mat` is very close to zero.
/// @return True if the scale could be extracted, false if the matrix is degenerate.
template <class T>
bool extractScalingAndShear (
    const Matrix44<T>& mat, Vec3<T>& scl, Vec3<T>& shr, bool exc = true);

/// Return the given 4x4 matrix with scaling and shear removed.
///
/// @param[in] mat The input matrix
/// @param[in] exc If true, throw an exception if the scaling in `mat` is very close to zero.
template <class T>
Matrix44<T> sansScalingAndShear (const Matrix44<T>& mat, bool exc = true);

/// Extract scaling and shear from the given 4x4 matrix in-place.
///
/// @param[in,out] result The output matrix
/// @param[in] mat The return value if `result` is degenerate
/// @param[in] exc If true, throw an exception if the scaling in `mat` is very close to zero.
template <class T>
void sansScalingAndShear (
    Matrix44<T>& result, const Matrix44<T>& mat, bool exc = true);

/// Remove scaling and shear from the given 4x4 matrix in place.
//
/// @param[in,out] mat The matrix to operate on
/// @param[in] exc If true, throw an exception if the scaling in `mat` is very close to zero.
/// @return True if the scale could be extracted, false if the matrix is degenerate.
template <class T>
bool removeScalingAndShear (Matrix44<T>& mat, bool exc = true);

/// Remove scaling and shear from the given 4x4 matrix in place, returning
/// the extracted values.
//
/// @param[in,out] mat The matrix to operate on
/// @param[out] scl The extracted scale
/// @param[out] shr The extracted shear
/// @param[in] exc If true, throw an exception if the scaling in `mat` is very close to zero.
/// @return True if the scale could be extracted, false if the matrix is degenerate.
template <class T>
bool extractAndRemoveScalingAndShear (
    Matrix44<T>& mat, Vec3<T>& scl, Vec3<T>& shr, bool exc = true);

/// Extract the rotation from the given 4x4 matrix in the form of XYZ
/// euler angles.
///
/// @param[in] mat The input matrix
/// @param[out] rot The extracted XYZ euler angle vector
template <class T> void extractEulerXYZ (const Matrix44<T>& mat, Vec3<T>& rot);

/// Extract the rotation from the given 4x4 matrix in the form of ZYX
/// euler angles.
///
/// @param[in] mat The input matrix
/// @param[out] rot The extracted ZYX euler angle vector
template <class T> void extractEulerZYX (const Matrix44<T>& mat, Vec3<T>& rot);

/// Extract the rotation from the given 4x4 matrix in the form of a quaternion.
///
/// @param[in] mat The input matrix
/// @return The extracted quaternion
template <class T> Quat<T> extractQuat (const Matrix44<T>& mat);

/// Extract the scaling, shear, rotation, and translation components
/// of the given 4x4 matrix. The values are such that:
///
///     M = S * H * R * T
///
/// @param[in] mat The input matrix
/// @param[out] s The extracted scale
/// @param[out] h The extracted shear
/// @param[out] r The extracted rotation
/// @param[out] t The extracted translation
/// @param[in] exc If true, throw an exception if the scaling in `mat` is very close to zero.
/// @param[in] rOrder The order with which to extract the rotation
/// @return True if the values could be extracted, false if the matrix is degenerate.
template <class T>
bool extractSHRT (
    const Matrix44<T>&       mat,
    Vec3<T>&                 s,
    Vec3<T>&                 h,
    Vec3<T>&                 r,
    Vec3<T>&                 t,
    bool                     exc /*= true*/,
    typename Euler<T>::Order rOrder);

/// Extract the scaling, shear, rotation, and translation components
/// of the given 4x4 matrix.
///
/// @param[in] mat The input matrix
/// @param[out] s The extracted scale
/// @param[out] h The extracted shear
/// @param[out] r The extracted rotation, in XYZ euler angles
/// @param[out] t The extracted translation
/// @param[in] exc If true, throw an exception if the scaling in `mat` is very close to zero.
/// @return True if the values could be extracted, false if the matrix is degenerate.
template <class T>
bool extractSHRT (
    const Matrix44<T>& mat,
    Vec3<T>&           s,
    Vec3<T>&           h,
    Vec3<T>&           r,
    Vec3<T>&           t,
    bool               exc = true);

/// Extract the scaling, shear, rotation, and translation components
/// of the given 4x4 matrix.
///
/// @param[in] mat The input matrix
/// @param[out] s The extracted scale
/// @param[out] h The extracted shear
/// @param[out] r The extracted rotation, in Euler angles
/// @param[out] t The extracted translation
/// @param[in] exc If true, throw an exception if the scaling in `mat` is very close to zero.
/// @return True if the values could be extracted, false if the matrix is degenerate.
template <class T>
bool extractSHRT (
    const Matrix44<T>& mat,
    Vec3<T>&           s,
    Vec3<T>&           h,
    Euler<T>&          r,
    Vec3<T>&           t,
    bool               exc = true);

/// Return true if the given scale can be removed from the given row
/// matrix, false if `scl` is small enough that the operation would
/// overflow. If `exc` is true, throw an exception on overflow.
template <class T>
bool checkForZeroScaleInRow (const T& scl, const Vec3<T>& row, bool exc = true);

/// Return the 4x4 outer product two 4-vectors
template <class T>
Matrix44<T> outerProduct (const Vec4<T>& a, const Vec4<T>& b);

///
/// Return a 4x4 matrix that rotates the vector `fromDirection` to `toDirection`
///
template <class T>
Matrix44<T>
rotationMatrix (const Vec3<T>& fromDirection, const Vec3<T>& toDirection);

///
/// Return a 4x4 matrix that rotates the `fromDir` vector
/// so that it points towards `toDir1.  You may also
/// specify that you want the up vector to be pointing
/// in a certain direction 1upDir`.
template <class T>
Matrix44<T> rotationMatrixWithUpDir (
    const Vec3<T>& fromDir, const Vec3<T>& toDir, const Vec3<T>& upDir);

///
/// Construct a 4x4 matrix that rotates the z-axis so that it points
/// towards `targetDir`.  You must also specify that you want the up
/// vector to be pointing in a certain direction `upDir`.
///
/// Notes: The following degenerate cases are handled:
/// (a) when the directions given by `toDir` and `upDir`
/// are parallel or opposite (the direction vectors must have a non-zero cross product);
/// (b) when any of the given direction vectors have zero length
///
/// @param[out] result The output matrix
/// @param[in] targetDir The target direction vector
/// @param[in] upDir The up direction vector
template <class T>
void
alignZAxisWithTargetDir (Matrix44<T>& result, Vec3<T> targetDir, Vec3<T> upDir);

/// Compute an orthonormal direct 4x4 frame from a position, an x axis
/// direction and a normal to the y axis. If the x axis and normal are
/// perpendicular, then the normal will have the same direction as the
/// z axis.
///
/// @param[in] p The position of the frame
/// @param[in] xDir The x axis direction of the frame
/// @param[in] normal A normal to the y axis of the frame
/// @return The orthonormal frame
template <class T>
Matrix44<T> computeLocalFrame (
    const Vec3<T>& p, const Vec3<T>& xDir, const Vec3<T>& normal);

/// Add a translate/rotate/scale offset to a 4x4 input frame
/// and put it in another frame of reference
///
/// @param[in] inMat Input frame
/// @param[in] tOffset Translation offset
/// @param[in] rOffset Rotation offset in degrees
/// @param[in] sOffset Scale offset
/// @param[in] ref Frame of reference
/// @return The offsetted frame
template <class T>
Matrix44<T> addOffset (
    const Matrix44<T>& inMat,
    const Vec3<T>&     tOffset,
    const Vec3<T>&     rOffset,
    const Vec3<T>&     sOffset,
    const Vec3<T>&     ref);

/// Compute 4x4 translate/rotate/scale matrix from `A` with the
/// rotate/scale of `B`.
///
/// @param[in] keepRotateA If true, keep rotate from matrix `A`, use `B` otherwise
/// @param[in] keepScaleA If true, keep scale  from matrix `A`, use `B` otherwise
/// @param[in] A Matrix A
/// @param[in] B Matrix B
/// @return Matrix `A` with tweaked rotation/scale
template <class T>
Matrix44<T> computeRSMatrix (
    bool               keepRotateA,
    bool               keepScaleA,
    const Matrix44<T>& A,
    const Matrix44<T>& B);

//
// Declarations for 3x3 matrix.
//

/// Extract the scaling component of the given 3x3 matrix.
///
/// @param[in] mat The input matrix
/// @param[out] scl The extracted scale, i.e. the output value
/// @param[in] exc If true, throw an exception if the scaling in `mat` is very close to zero.
/// @return True if the scale could be extracted, false if the matrix is degenerate.
template <class T>
bool extractScaling (const Matrix33<T>& mat, Vec2<T>& scl, bool exc = true);

/// Return the given 3x3 matrix with scaling removed.
///
/// @param[in] mat The input matrix
/// @param[in] exc If true, throw an exception if the scaling in `mat`
template <class T>
Matrix33<T> sansScaling (const Matrix33<T>& mat, bool exc = true);

/// Remove scaling from the given 3x3 matrix in place.  Return true if the
/// scale could be successfully extracted, false if the matrix is
/// degenerate.
//
/// @param[in] mat The matrix to operate on
/// @param[in] exc If true, throw an exception if the scaling in `mat` is very close to zero.
/// @return True if the scale could be extracted, false if the matrix is degenerate.
template <class T> bool removeScaling (Matrix33<T>& mat, bool exc = true);

/// Extract the scaling and shear components of the given 3x3 matrix.
/// Return true if the scale could be successfully extracted, false if
/// the matrix is degenerate.
///
/// @param[in] mat The input matrix
/// @param[out] scl The extracted scale
/// @param[out] shr The extracted shear
/// @param[in] exc If true, throw an exception if the scaling in `mat` is very close to zero.
/// @return True if the scale could be extracted, false if the matrix is degenerate.
template <class T>
bool extractScalingAndShear (
    const Matrix33<T>& mat, Vec2<T>& scl, T& shr, bool exc = true);

/// Return the given 3x3 matrix with scaling and shear removed.
///
/// @param[in] mat The input matrix
/// @param[in] exc If true, throw an exception if the scaling in `mat` is very close to zero.
template <class T>
Matrix33<T> sansScalingAndShear (const Matrix33<T>& mat, bool exc = true);

/// Remove scaling and shear from the given 3x3e matrix in place.
//
/// @param[in,out] mat The matrix to operate on
/// @param[in] exc If true, throw an exception if the scaling in `mat` is very close to zero.
/// @return True if the scale could be extracted, false if the matrix is degenerate.
template <class T>
bool removeScalingAndShear (Matrix33<T>& mat, bool exc = true);

/// Remove scaling and shear from the given 3x3 matrix in place, returning
/// the extracted values.
//
/// @param[in,out] mat The matrix to operate on
/// @param[out] scl The extracted scale
/// @param[out] shr The extracted shear
/// @param[in] exc If true, throw an exception if the scaling in `mat` is very close to zero.
/// @return True if the scale could be extracted, false if the matrix is degenerate.
template <class T>
bool extractAndRemoveScalingAndShear (
    Matrix33<T>& mat, Vec2<T>& scl, T& shr, bool exc = true);

/// Extract the rotation from the given 2x2 matrix
///
/// @param[in] mat The input matrix
/// @param[out] rot The extracted rotation value
template <class T> void extractEuler (const Matrix22<T>& mat, T& rot);

/// Extract the rotation from the given 3x3 matrix
///
/// @param[in] mat The input matrix
/// @param[out] rot The extracted rotation value
template <class T> void extractEuler (const Matrix33<T>& mat, T& rot);

/// Extract the scaling, shear, rotation, and translation components
/// of the given 3x3 matrix. The values are such that:
///
///     M = S * H * R * T
///
/// @param[in] mat The input matrix
/// @param[out] s The extracted scale
/// @param[out] h The extracted shear
/// @param[out] r The extracted rotation
/// @param[out] t The extracted translation
/// @param[in] exc If true, throw an exception if the scaling in `mat` is very close to zero.
/// @return True if the values could be extracted, false if the matrix is degenerate.
template <class T>
bool extractSHRT (
    const Matrix33<T>& mat,
    Vec2<T>&           s,
    T&                 h,
    T&                 r,
    Vec2<T>&           t,
    bool               exc = true);

/// Return true if the given scale can be removed from the given row
/// matrix, false if `scl` is small enough that the operation would
/// overflow. If `exc` is true, throw an exception on overflow.
template <class T>
bool checkForZeroScaleInRow (const T& scl, const Vec2<T>& row, bool exc = true);

/// Return the 3xe outer product two 3-vectors
template <class T>
Matrix33<T> outerProduct (const Vec3<T>& a, const Vec3<T>& b);

//------------------------------
// Implementation for 4x4 Matrix
//------------------------------

template <class T>
bool
extractScaling (const Matrix44<T>& mat, Vec3<T>& scl, bool exc)
{
    Vec3<T>     shr;
    Matrix44<T> M (mat);

    if (!extractAndRemoveScalingAndShear (M, scl, shr, exc)) return false;

    return true;
}

template <class T>
Matrix44<T>
sansScaling (const Matrix44<T>& mat, bool exc)
{
    Vec3<T> scl;
    Vec3<T> shr;
    Vec3<T> rot;
    Vec3<T> tran;

    if (!extractSHRT (mat, scl, shr, rot, tran, exc)) return mat;

    Matrix44<T> M;

    M.translate (tran);
    M.rotate (rot);
    M.shear (shr);

    return M;
}

template <class T>
bool
removeScaling (Matrix44<T>& mat, bool exc)
{
    Vec3<T> scl;
    Vec3<T> shr;
    Vec3<T> rot;
    Vec3<T> tran;

    if (!extractSHRT (mat, scl, shr, rot, tran, exc)) return false;

    mat.makeIdentity ();
    mat.translate (tran);
    mat.rotate (rot);
    mat.shear (shr);

    return true;
}

template <class T>
bool
extractScalingAndShear (
    const Matrix44<T>& mat, Vec3<T>& scl, Vec3<T>& shr, bool exc)
{
    Matrix44<T> M (mat);

    if (!extractAndRemoveScalingAndShear (M, scl, shr, exc)) return false;

    return true;
}

template <class T>
Matrix44<T>
sansScalingAndShear (const Matrix44<T>& mat, bool exc)
{
    Vec3<T>     scl;
    Vec3<T>     shr;
    Matrix44<T> M (mat);

    if (!extractAndRemoveScalingAndShear (M, scl, shr, exc)) return mat;

    return M;
}

template <class T>
void
sansScalingAndShear (Matrix44<T>& result, const Matrix44<T>& mat, bool exc)
{
    Vec3<T> scl;
    Vec3<T> shr;

    if (!extractAndRemoveScalingAndShear (result, scl, shr, exc)) result = mat;
}

template <class T>
bool
removeScalingAndShear (Matrix44<T>& mat, bool exc)
{
    Vec3<T> scl;
    Vec3<T> shr;

    if (!extractAndRemoveScalingAndShear (mat, scl, shr, exc)) return false;

    return true;
}

template <class T>
bool
extractAndRemoveScalingAndShear (
    Matrix44<T>& mat, Vec3<T>& scl, Vec3<T>& shr, bool exc)
{
    //
    // This implementation follows the technique described in the paper by
    // Spencer W. Thomas in the Graphics Gems II article: "Decomposing a
    // Matrix into Simple Transformations", p. 320.
    //

    Vec3<T> row[3];

    row[0] = Vec3<T> (mat[0][0], mat[0][1], mat[0][2]);
    row[1] = Vec3<T> (mat[1][0], mat[1][1], mat[1][2]);
    row[2] = Vec3<T> (mat[2][0], mat[2][1], mat[2][2]);

    T maxVal = 0;
    for (int i = 0; i < 3; i++)
        for (int j = 0; j < 3; j++)
            if (IMATH_INTERNAL_NAMESPACE::abs (row[i][j]) > maxVal)
                maxVal = IMATH_INTERNAL_NAMESPACE::abs (row[i][j]);

    //
    // We normalize the 3x3 matrix here.
    // It was noticed that this can improve numerical stability significantly,
    // especially when many of the upper 3x3 matrix's coefficients are very
    // close to zero; we correct for this step at the end by multiplying the
    // scaling factors by maxVal at the end (shear and rotation are not
    // affected by the normalization).

    if (maxVal != 0)
    {
        for (int i = 0; i < 3; i++)
            if (!checkForZeroScaleInRow (maxVal, row[i], exc))
                return false;
            else
                row[i] /= maxVal;
    }

    // Compute X scale factor.
    scl.x = row[0].length ();
    if (!checkForZeroScaleInRow (scl.x, row[0], exc)) return false;

    // Normalize first row.
    row[0] /= scl.x;

    // An XY shear factor will shear the X coord. as the Y coord. changes.
    // There are 6 combinations (XY, XZ, YZ, YX, ZX, ZY), although we only
    // extract the first 3 because we can effect the last 3 by shearing in
    // XY, XZ, YZ combined rotations and scales.
    //
    // shear matrix <   1,  YX,  ZX,  0,
    //                 XY,   1,  ZY,  0,
    //                 XZ,  YZ,   1,  0,
    //                  0,   0,   0,  1 >

    // Compute XY shear factor and make 2nd row orthogonal to 1st.
    shr[0] = row[0].dot (row[1]);
    row[1] -= shr[0] * row[0];

    // Now, compute Y scale.
    scl.y = row[1].length ();
    if (!checkForZeroScaleInRow (scl.y, row[1], exc)) return false;

    // Normalize 2nd row and correct the XY shear factor for Y scaling.
    row[1] /= scl.y;
    shr[0] /= scl.y;

    // Compute XZ and YZ shears, orthogonalize 3rd row.
    shr[1] = row[0].dot (row[2]);
    row[2] -= shr[1] * row[0];
    shr[2] = row[1].dot (row[2]);
    row[2] -= shr[2] * row[1];

    // Next, get Z scale.
    scl.z = row[2].length ();
    if (!checkForZeroScaleInRow (scl.z, row[2], exc)) return false;

    // Normalize 3rd row and correct the XZ and YZ shear factors for Z scaling.
    row[2] /= scl.z;
    shr[1] /= scl.z;
    shr[2] /= scl.z;

    // At this point, the upper 3x3 matrix in mat is orthonormal.
    // Check for a coordinate system flip. If the determinant
    // is less than zero, then negate the matrix and the scaling factors.
    if (row[0].dot (row[1].cross (row[2])) < 0)
        for (int i = 0; i < 3; i++)
        {
            scl[i] *= -1;
            row[i] *= -1;
        }

    // Copy over the orthonormal rows into the returned matrix.
    // The upper 3x3 matrix in mat is now a rotation matrix.
    for (int i = 0; i < 3; i++)
    {
        mat[i][0] = row[i][0];
        mat[i][1] = row[i][1];
        mat[i][2] = row[i][2];
    }

    // Correct the scaling factors for the normalization step that we
    // performed above; shear and rotation are not affected by the
    // normalization.
    scl *= maxVal;

    return true;
}

template <class T>
void
extractEulerXYZ (const Matrix44<T>& mat, Vec3<T>& rot)
{
    //
    // Normalize the local x, y and z axes to remove scaling.
    //

    Vec3<T> i (mat[0][0], mat[0][1], mat[0][2]);
    Vec3<T> j (mat[1][0], mat[1][1], mat[1][2]);
    Vec3<T> k (mat[2][0], mat[2][1], mat[2][2]);

    i.normalize ();
    j.normalize ();
    k.normalize ();

    Matrix44<T> M (
        i[0],
        i[1],
        i[2],
        0,
        j[0],
        j[1],
        j[2],
        0,
        k[0],
        k[1],
        k[2],
        0,
        0,
        0,
        0,
        1);

    //
    // Extract the first angle, rot.x.
    //

    rot.x = std::atan2 (M[1][2], M[2][2]);

    //
    // Remove the rot.x rotation from M, so that the remaining
    // rotation, N, is only around two axes, and gimbal lock
    // cannot occur.
    //

    Matrix44<T> N;
    N.rotate (Vec3<T> (-rot.x, 0, 0));
    N = N * M;

    //
    // Extract the other two angles, rot.y and rot.z, from N.
    //

    T cy  = std::sqrt (N[0][0] * N[0][0] + N[0][1] * N[0][1]);
    rot.y = std::atan2 (-N[0][2], cy);
    rot.z = std::atan2 (-N[1][0], N[1][1]);
}

template <class T>
void
extractEulerZYX (const Matrix44<T>& mat, Vec3<T>& rot)
{
    //
    // Normalize the local x, y and z axes to remove scaling.
    //

    Vec3<T> i (mat[0][0], mat[0][1], mat[0][2]);
    Vec3<T> j (mat[1][0], mat[1][1], mat[1][2]);
    Vec3<T> k (mat[2][0], mat[2][1], mat[2][2]);

    i.normalize ();
    j.normalize ();
    k.normalize ();

    Matrix44<T> M (
        i[0],
        i[1],
        i[2],
        0,
        j[0],
        j[1],
        j[2],
        0,
        k[0],
        k[1],
        k[2],
        0,
        0,
        0,
        0,
        1);

    //
    // Extract the first angle, rot.x.
    //

    rot.x = -std::atan2 (M[1][0], M[0][0]);

    //
    // Remove the x rotation from M, so that the remaining
    // rotation, N, is only around two axes, and gimbal lock
    // cannot occur.
    //

    Matrix44<T> N;
    N.rotate (Vec3<T> (0, 0, -rot.x));
    N = N * M;

    //
    // Extract the other two angles, rot.y and rot.z, from N.
    //

    T cy  = std::sqrt (N[2][2] * N[2][2] + N[2][1] * N[2][1]);
    rot.y = -std::atan2 (-N[2][0], cy);
    rot.z = -std::atan2 (-N[1][2], N[1][1]);
}

template <class T>
Quat<T>
extractQuat (const Matrix44<T>& mat)
{
    T       tr, s;
    T       q[4];
    int     i, j, k;
    Quat<T> quat;

    int nxt[3] = {1, 2, 0};
    tr         = mat[0][0] + mat[1][1] + mat[2][2];

    // check the diagonal
    if (tr > 0.0)
    {
        s      = std::sqrt (tr + T (1.0));
        quat.r = s / T (2.0);
        s      = T (0.5) / s;

        quat.v.x = (mat[1][2] - mat[2][1]) * s;
        quat.v.y = (mat[2][0] - mat[0][2]) * s;
        quat.v.z = (mat[0][1] - mat[1][0]) * s;
    }
    else
    {
        // diagonal is negative
        i = 0;
        if (mat[1][1] > mat[0][0]) i = 1;
        if (mat[2][2] > mat[i][i]) i = 2;

        j = nxt[i];
        k = nxt[j];
        s = std::sqrt ((mat[i][i] - (mat[j][j] + mat[k][k])) + T (1.0));

        q[i] = s * T (0.5);
        if (s != T (0.0)) s = T (0.5) / s;

        q[3] = (mat[j][k] - mat[k][j]) * s;
        q[j] = (mat[i][j] + mat[j][i]) * s;
        q[k] = (mat[i][k] + mat[k][i]) * s;

        quat.v.x = q[0];
        quat.v.y = q[1];
        quat.v.z = q[2];
        quat.r   = q[3];
    }

    return quat;
}

template <class T>
bool
extractSHRT (
    const Matrix44<T>&       mat,
    Vec3<T>&                 s,
    Vec3<T>&                 h,
    Vec3<T>&                 r,
    Vec3<T>&                 t,
    bool                     exc /* = true */,
    typename Euler<T>::Order rOrder /* = Euler<T>::XYZ */)
{
    Matrix44<T> rot;

    rot = mat;
    if (!extractAndRemoveScalingAndShear (rot, s, h, exc)) return false;

    extractEulerXYZ (rot, r);

    t.x = mat[3][0];
    t.y = mat[3][1];
    t.z = mat[3][2];

    if (rOrder != Euler<T>::XYZ)
    {
        Euler<T> eXYZ (r, Euler<T>::XYZ);
        Euler<T> e (eXYZ, rOrder);
        r = e.toXYZVector ();
    }

    return true;
}

template <class T>
bool
extractSHRT (
    const Matrix44<T>& mat,
    Vec3<T>&           s,
    Vec3<T>&           h,
    Vec3<T>&           r,
    Vec3<T>&           t,
    bool               exc)
{
    return extractSHRT (mat, s, h, r, t, exc, Euler<T>::XYZ);
}

template <class T>
bool
extractSHRT (
    const Matrix44<T>& mat,
    Vec3<T>&           s,
    Vec3<T>&           h,
    Euler<T>&          r,
    Vec3<T>&           t,
    bool               exc /* = true */)
{
    return extractSHRT (mat, s, h, r, t, exc, r.order ());
}

template <class T>
bool
checkForZeroScaleInRow (const T& scl, const Vec3<T>& row, bool exc /* = true */)
{
    for (int i = 0; i < 3; i++)
    {
        if ((abs (scl) < 1 &&
             abs (row[i]) >= std::numeric_limits<T>::max () * abs (scl)))
        {
            if (exc)
                throw std::domain_error ("Cannot remove zero scaling "
                                         "from matrix.");
            else
                return false;
        }
    }

    return true;
}

template <class T>
Matrix44<T>
outerProduct (const Vec4<T>& a, const Vec4<T>& b)
{
    return Matrix44<T> (
        a.x * b.x,
        a.x * b.y,
        a.x * b.z,
        a.x * b.w,
        a.y * b.x,
        a.y * b.y,
        a.y * b.z,
        a.x * b.w,
        a.z * b.x,
        a.z * b.y,
        a.z * b.z,
        a.x * b.w,
        a.w * b.x,
        a.w * b.y,
        a.w * b.z,
        a.w * b.w);
}

template <class T>
Matrix44<T>
rotationMatrix (const Vec3<T>& from, const Vec3<T>& to)
{
    Quat<T> q;
    q.setRotation (from, to);
    return q.toMatrix44 ();
}

template <class T>
Matrix44<T>
rotationMatrixWithUpDir (
    const Vec3<T>& fromDir, const Vec3<T>& toDir, const Vec3<T>& upDir)
{
    //
    // The goal is to obtain a rotation matrix that takes
    // "fromDir" to "toDir".  We do this in two steps and
    // compose the resulting rotation matrices;
    //    (a) rotate "fromDir" into the z-axis
    //    (b) rotate the z-axis into "toDir"
    //

    // The from direction must be non-zero; but we allow zero to and up dirs.
    if (fromDir.length () == 0)
        return Matrix44<T> ();

    else
    {
        Matrix44<T> zAxis2FromDir (UNINITIALIZED);
        alignZAxisWithTargetDir (zAxis2FromDir, fromDir, Vec3<T> (0, 1, 0));

        Matrix44<T> fromDir2zAxis = zAxis2FromDir.transposed ();

        Matrix44<T> zAxis2ToDir (UNINITIALIZED);
        alignZAxisWithTargetDir (zAxis2ToDir, toDir, upDir);

        return fromDir2zAxis * zAxis2ToDir;
    }
}

template <class T>
void
alignZAxisWithTargetDir (Matrix44<T>& result, Vec3<T> targetDir, Vec3<T> upDir)
{
    //
    // Ensure that the target direction is non-zero.
    //

    if (targetDir.length () == 0) targetDir = Vec3<T> (0, 0, 1);

    //
    // Ensure that the up direction is non-zero.
    //

    if (upDir.length () == 0) upDir = Vec3<T> (0, 1, 0);

    //
    // Check for degeneracies.  If the upDir and targetDir are parallel
    // or opposite, then compute a new, arbitrary up direction that is
    // not parallel or opposite to the targetDir.
    //

    if (upDir.cross (targetDir).length () == 0)
    {
        upDir = targetDir.cross (Vec3<T> (1, 0, 0));
        if (upDir.length () == 0) upDir = targetDir.cross (Vec3<T> (0, 0, 1));
    }

    //
    // Compute the x-, y-, and z-axis vectors of the new coordinate system.
    //

    Vec3<T> targetPerpDir = upDir.cross (targetDir);
    Vec3<T> targetUpDir   = targetDir.cross (targetPerpDir);

    //
    // Rotate the x-axis into targetPerpDir (row 0),
    // rotate the y-axis into targetUpDir   (row 1),
    // rotate the z-axis into targetDir     (row 2).
    //

    Vec3<T> row[3];
    row[0] = targetPerpDir.normalized ();
    row[1] = targetUpDir.normalized ();
    row[2] = targetDir.normalized ();

    result.x[0][0] = row[0][0];
    result.x[0][1] = row[0][1];
    result.x[0][2] = row[0][2];
    result.x[0][3] = (T) 0;

    result.x[1][0] = row[1][0];
    result.x[1][1] = row[1][1];
    result.x[1][2] = row[1][2];
    result.x[1][3] = (T) 0;

    result.x[2][0] = row[2][0];
    result.x[2][1] = row[2][1];
    result.x[2][2] = row[2][2];
    result.x[2][3] = (T) 0;

    result.x[3][0] = (T) 0;
    result.x[3][1] = (T) 0;
    result.x[3][2] = (T) 0;
    result.x[3][3] = (T) 1;
}

// Compute an orthonormal direct frame from : a position, an x axis direction and a normal to the y axis
// If the x axis and normal are perpendicular, then the normal will have the same direction as the z axis.
// Inputs are :
//     -the position of the frame
//     -the x axis direction of the frame
//     -a normal to the y axis of the frame
// Return is the orthonormal frame
template <class T>
Matrix44<T>
computeLocalFrame (const Vec3<T>& p, const Vec3<T>& xDir, const Vec3<T>& normal)
{
    Vec3<T> _xDir (xDir);
    Vec3<T> x = _xDir.normalize ();
    Vec3<T> y = (normal % x).normalize ();
    Vec3<T> z = (x % y).normalize ();

    Matrix44<T> L;
    L[0][0] = x[0];
    L[0][1] = x[1];
    L[0][2] = x[2];
    L[0][3] = 0.0;

    L[1][0] = y[0];
    L[1][1] = y[1];
    L[1][2] = y[2];
    L[1][3] = 0.0;

    L[2][0] = z[0];
    L[2][1] = z[1];
    L[2][2] = z[2];
    L[2][3] = 0.0;

    L[3][0] = p[0];
    L[3][1] = p[1];
    L[3][2] = p[2];
    L[3][3] = 1.0;

    return L;
}

/// Add a translate/rotate/scale offset to an input frame and put it
/// in another frame of reference.
///
/// @param inMat input frame
/// @param tOffset translate offset
/// @param rOffset rotate offset in degrees
/// @param sOffset scale offset
/// @param ref Frame of reference
/// @return The offsetted frame
template <class T>
Matrix44<T>
addOffset (
    const Matrix44<T>& inMat,
    const Vec3<T>&     tOffset,
    const Vec3<T>&     rOffset,
    const Vec3<T>&     sOffset,
    const Matrix44<T>& ref)
{
    Matrix44<T> O;

    Vec3<T> _rOffset (rOffset);
    _rOffset *= T(M_PI / 180.0);
    O.rotate (_rOffset);

    O[3][0] = tOffset[0];
    O[3][1] = tOffset[1];
    O[3][2] = tOffset[2];

    Matrix44<T> S;
    S.scale (sOffset);

    Matrix44<T> X = S * O * inMat * ref;

    return X;
}

// Compute Translate/Rotate/Scale matrix from matrix A with the Rotate/Scale of Matrix B
// Inputs are :
//      -keepRotateA : if true keep rotate from matrix A, use B otherwise
//      -keepScaleA  : if true keep scale  from matrix A, use B otherwise
//      -Matrix A
//      -Matrix B
// Return Matrix A with tweaked rotation/scale
template <class T>
Matrix44<T>
computeRSMatrix (
    bool               keepRotateA,
    bool               keepScaleA,
    const Matrix44<T>& A,
    const Matrix44<T>& B)
{
    Vec3<T> as, ah, ar, at;
    if (!extractSHRT (A, as, ah, ar, at))
        throw std::domain_error ("degenerate A matrix in computeRSMatrix");

    Vec3<T> bs, bh, br, bt;
    if (!extractSHRT (B, bs, bh, br, bt))
        throw std::domain_error ("degenerate B matrix in computeRSMatrix");

    if (!keepRotateA) ar = br;

    if (!keepScaleA) as = bs;

    Matrix44<T> mat;
    mat.makeIdentity ();
    mat.translate (at);
    mat.rotate (ar);
    mat.scale (as);

    return mat;
}

//-----------------------------------------------------------------------------
// Implementation for 3x3 Matrix
//------------------------------

template <class T>
bool
extractScaling (const Matrix33<T>& mat, Vec2<T>& scl, bool exc)
{
    T           shr;
    Matrix33<T> M (mat);

    if (!extractAndRemoveScalingAndShear (M, scl, shr, exc)) return false;

    return true;
}

template <class T>
Matrix33<T>
sansScaling (const Matrix33<T>& mat, bool exc)
{
    Vec2<T> scl;
    T       shr;
    T       rot;
    Vec2<T> tran;

    if (!extractSHRT (mat, scl, shr, rot, tran, exc)) return mat;

    Matrix33<T> M;

    M.translate (tran);
    M.rotate (rot);
    M.shear (shr);

    return M;
}

template <class T>
bool
removeScaling (Matrix33<T>& mat, bool exc)
{
    Vec2<T> scl;
    T       shr;
    T       rot;
    Vec2<T> tran;

    if (!extractSHRT (mat, scl, shr, rot, tran, exc)) return false;

    mat.makeIdentity ();
    mat.translate (tran);
    mat.rotate (rot);
    mat.shear (shr);

    return true;
}

template <class T>
bool
extractScalingAndShear (const Matrix33<T>& mat, Vec2<T>& scl, T& shr, bool exc)
{
    Matrix33<T> M (mat);

    if (!extractAndRemoveScalingAndShear (M, scl, shr, exc)) return false;

    return true;
}

template <class T>
Matrix33<T>
sansScalingAndShear (const Matrix33<T>& mat, bool exc)
{
    Vec2<T>     scl;
    T           shr;
    Matrix33<T> M (mat);

    if (!extractAndRemoveScalingAndShear (M, scl, shr, exc)) return mat;

    return M;
}

template <class T>
bool
removeScalingAndShear (Matrix33<T>& mat, bool exc)
{
    Vec2<T> scl;
    T       shr;

    if (!extractAndRemoveScalingAndShear (mat, scl, shr, exc)) return false;

    return true;
}

template <class T>
bool
extractAndRemoveScalingAndShear (
    Matrix33<T>& mat, Vec2<T>& scl, T& shr, bool exc)
{
    Vec2<T> row[2];

    row[0] = Vec2<T> (mat[0][0], mat[0][1]);
    row[1] = Vec2<T> (mat[1][0], mat[1][1]);

    T maxVal = 0;
    for (int i = 0; i < 2; i++)
        for (int j = 0; j < 2; j++)
            if (IMATH_INTERNAL_NAMESPACE::abs (mat[i][j]) > maxVal)
                maxVal = IMATH_INTERNAL_NAMESPACE::abs (mat[i][j]);

    //
    // We normalize the 2x2 matrix here.
    // It was noticed that this can improve numerical stability significantly,
    // especially when many of the upper 2x2 matrix's coefficients are very
    // close to zero; we correct for this step at the end by multiplying the
    // scaling factors by maxVal at the end (shear and rotation are not
    // affected by the normalization).

    if (maxVal != 0)
    {
        for (int i = 0; i < 2; i++)
            if (!checkForZeroScaleInRow (maxVal, row[i], exc))
                return false;
            else
                row[i] /= maxVal;
    }

    // Compute X scale factor.
    scl.x = row[0].length ();
    if (!checkForZeroScaleInRow (scl.x, row[0], exc)) return false;

    // Normalize first row.
    row[0] /= scl.x;

    // An XY shear factor will shear the X coord. as the Y coord. changes.
    // There are 2 combinations (XY, YX), although we only extract the XY
    // shear factor because we can effect the an YX shear factor by
    // shearing in XY combined with rotations and scales.
    //
    // shear matrix <   1,  YX,  0,
    //                 XY,   1,  0,
    //                  0,   0,  1 >

    // Compute XY shear factor and make 2nd row orthogonal to 1st.
    shr = row[0].dot (row[1]);
    row[1] -= shr * row[0];

    // Now, compute Y scale.
    scl.y = row[1].length ();
    if (!checkForZeroScaleInRow (scl.y, row[1], exc)) return false;

    // Normalize 2nd row and correct the XY shear factor for Y scaling.
    row[1] /= scl.y;
    shr /= scl.y;

    // At this point, the upper 2x2 matrix in mat is orthonormal.
    // Check for a coordinate system flip. If the determinant
    // is -1, then flip the rotation matrix and adjust the scale(Y)
    // and shear(XY) factors to compensate.
    if (row[0].x * row[1].y - row[0].y * row[1].x < 0)
    {
        row[1].x *= -1;
        row[1].y *= -1;
        scl.y *= -1;
        shr *= -1;
    }

    // Copy over the orthonormal rows into the returned matrix.
    // The upper 2x2 matrix in mat is now a rotation matrix.
    for (int i = 0; i < 2; i++)
    {
        mat[i][0] = row[i].x;
        mat[i][1] = row[i].y;
    }

    scl *= maxVal;

    return true;
}

template <class T>
void
extractEuler (const Matrix22<T>& mat, T& rot)
{
    //
    // Normalize the local x and y axes to remove scaling.
    //

    Vec2<T> i (mat[0][0], mat[0][1]);
    Vec2<T> j (mat[1][0], mat[1][1]);

    i.normalize ();
    j.normalize ();

    //
    // Extract the angle, rot.
    //

    rot = -std::atan2 (j[0], i[0]);
}

template <class T>
void
extractEuler (const Matrix33<T>& mat, T& rot)
{
    //
    // Normalize the local x and y axes to remove scaling.
    //

    Vec2<T> i (mat[0][0], mat[0][1]);
    Vec2<T> j (mat[1][0], mat[1][1]);

    i.normalize ();
    j.normalize ();

    //
    // Extract the angle, rot.
    //

    rot = -std::atan2 (j[0], i[0]);
}

template <class T>
bool
extractSHRT (
    const Matrix33<T>& mat, Vec2<T>& s, T& h, T& r, Vec2<T>& t, bool exc)
{
    Matrix33<T> rot;

    rot = mat;
    if (!extractAndRemoveScalingAndShear (rot, s, h, exc)) return false;

    extractEuler (rot, r);

    t.x = mat[2][0];
    t.y = mat[2][1];

    return true;
}

/// @cond Doxygen_Suppress
template <class T>
bool
checkForZeroScaleInRow (const T& scl, const Vec2<T>& row, bool exc /* = true */)
{
    for (int i = 0; i < 2; i++)
    {
        if ((abs (scl) < 1 &&
             abs (row[i]) >= std::numeric_limits<T>::max () * abs (scl)))
        {
            if (exc)
                throw std::domain_error (
                    "Cannot remove zero scaling from matrix.");
            else
                return false;
        }
    }

    return true;
}
/// @endcond

template <class T>
Matrix33<T>
outerProduct (const Vec3<T>& a, const Vec3<T>& b)
{
    return Matrix33<T> (
        a.x * b.x,
        a.x * b.y,
        a.x * b.z,
        a.y * b.x,
        a.y * b.y,
        a.y * b.z,
        a.z * b.x,
        a.z * b.y,
        a.z * b.z);
}

/// Computes the translation and rotation that brings the 'from' points
/// as close as possible to the 'to' points under the Frobenius norm.
/// To be more specific, let x be the matrix of 'from' points and y be
/// the matrix of 'to' points, we want to find the matrix A of the form
///    [ R t ]
///    [ 0 1 ]
/// that minimizes
///     || (A*x - y)^T * W * (A*x - y) ||_F
/// If doScaling is true, then a uniform scale is allowed also.
/// @param A From points
/// @param B To points
/// @param weights Per-point weights
/// @param numPoints The number of points in `A`, `B`, and `weights` (must be equal)
/// @param doScaling If true, include a scaling transformation
/// @return The procrustes transformation
template <typename T>
M44d procrustesRotationAndTranslation (
    const Vec3<T>* A,
    const Vec3<T>* B,
    const T*       weights,
    const size_t   numPoints,
    const bool     doScaling = false);

/// Computes the translation and rotation that brings the 'from' points
/// as close as possible to the 'to' points under the Frobenius norm.
/// To be more specific, let x be the matrix of 'from' points and y be
/// the matrix of 'to' points, we want to find the matrix A of the form
///    [ R t ]
///    [ 0 1 ]
/// that minimizes
///     || (A*x - y)^T * W * (A*x - y) ||_F
/// If doScaling is true, then a uniform scale is allowed also.
/// @param A From points
/// @param B To points
/// @param numPoints The number of points in `A` and `B` (must be equal)
/// @param doScaling If true, include a scaling transformation
/// @return The procrustes transformation
template <typename T>
M44d procrustesRotationAndTranslation (
    const Vec3<T>* A,
    const Vec3<T>* B,
    const size_t   numPoints,
    const bool     doScaling = false);

/// Compute the SVD of a 3x3 matrix using Jacobi transformations.  This method
/// should be quite accurate (competitive with LAPACK) even for poorly
/// conditioned matrices, and because it has been written specifically for the
/// 3x3/4x4 case it is much faster than calling out to LAPACK.
///
/// The SVD of a 3x3/4x4 matrix A is defined as follows:
///     A = U * S * V^T
/// where S is the diagonal matrix of singular values and both U and V are
/// orthonormal.  By convention, the entries S are all positive and sorted from
/// the largest to the smallest.  However, some uses of this function may
/// require that the matrix U*V^T have positive determinant; in this case, we
/// may make the smallest singular value negative to ensure that this is
/// satisfied.
///
/// Currently only available for single- and double-precision matrices.
template <typename T>
void jacobiSVD (
    const Matrix33<T>& A,
    Matrix33<T>&       U,
    Vec3<T>&           S,
    Matrix33<T>&       V,
    const T            tol = std::numeric_limits<T>::epsilon (),
    const bool         forcePositiveDeterminant = false);

/// Compute the SVD of a 3x3 matrix using Jacobi transformations.  This method
/// should be quite accurate (competitive with LAPACK) even for poorly
/// conditioned matrices, and because it has been written specifically for the
/// 3x3/4x4 case it is much faster than calling out to LAPACK.
///
/// The SVD of a 3x3/4x4 matrix A is defined as follows:
///     A = U * S * V^T
/// where S is the diagonal matrix of singular values and both U and V are
/// orthonormal.  By convention, the entries S are all positive and sorted from
/// the largest to the smallest.  However, some uses of this function may
/// require that the matrix U*V^T have positive determinant; in this case, we
/// may make the smallest singular value negative to ensure that this is
/// satisfied.
///
/// Currently only available for single- and double-precision matrices.
template <typename T>
void jacobiSVD (
    const Matrix44<T>& A,
    Matrix44<T>&       U,
    Vec4<T>&           S,
    Matrix44<T>&       V,
    const T            tol = std::numeric_limits<T>::epsilon (),
    const bool         forcePositiveDeterminant = false);

/// Compute the eigenvalues (S) and the eigenvectors (V) of a real
/// symmetric matrix using Jacobi transformation, using a given
/// tolerance `tol`.
///
/// Jacobi transformation of a 3x3/4x4 matrix A outputs S and V:
///     A = V * S * V^T
/// where V is orthonormal and S is the diagonal matrix of eigenvalues.
/// Input matrix A must be symmetric. A is also modified during
/// the computation so that upper diagonal entries of A become zero.
template <typename T>
void
jacobiEigenSolver (Matrix33<T>& A, Vec3<T>& S, Matrix33<T>& V, const T tol);

/// Compute the eigenvalues (S) and the eigenvectors (V) of
/// a real symmetric matrix using Jacobi transformation.
///
/// Jacobi transformation of a 3x3/4x4 matrix A outputs S and V:
///     A = V * S * V^T
/// where V is orthonormal and S is the diagonal matrix of eigenvalues.
/// Input matrix A must be symmetric. A is also modified during
/// the computation so that upper diagonal entries of A become zero.
template <typename T>
inline void
jacobiEigenSolver (Matrix33<T>& A, Vec3<T>& S, Matrix33<T>& V)
{
    jacobiEigenSolver (A, S, V, std::numeric_limits<T>::epsilon ());
}

/// Compute the eigenvalues (S) and the eigenvectors (V) of a real
/// symmetric matrix using Jacobi transformation, using a given
/// tolerance `tol`.
///
/// Jacobi transformation of a 3x3/4x4 matrix A outputs S and V:
///     A = V * S * V^T
/// where V is orthonormal and S is the diagonal matrix of eigenvalues.
/// Input matrix A must be symmetric. A is also modified during
/// the computation so that upper diagonal entries of A become zero.
template <typename T>
void
jacobiEigenSolver (Matrix44<T>& A, Vec4<T>& S, Matrix44<T>& V, const T tol);

/// Compute the eigenvalues (S) and the eigenvectors (V) of
/// a real symmetric matrix using Jacobi transformation.
///
/// Jacobi transformation of a 3x3/4x4 matrix A outputs S and V:
///     A = V * S * V^T
/// where V is orthonormal and S is the diagonal matrix of eigenvalues.
/// Input matrix A must be symmetric. A is also modified during
/// the computation so that upper diagonal entries of A become zero.
template <typename T>
inline void
jacobiEigenSolver (Matrix44<T>& A, Vec4<T>& S, Matrix44<T>& V)
{
    jacobiEigenSolver (A, S, V, std::numeric_limits<T>::epsilon ());
}

/// Compute a eigenvector corresponding to the abs max eigenvalue
/// of a real symmetric matrix using Jacobi transformation.
template <typename TM, typename TV> void maxEigenVector (TM& A, TV& S);

/// Compute a eigenvector corresponding to the abs min eigenvalue
/// of a real symmetric matrix using Jacobi transformation.
template <typename TM, typename TV> void minEigenVector (TM& A, TV& S);

IMATH_INTERNAL_NAMESPACE_HEADER_EXIT

#endif // INCLUDED_IMATHMATRIXALGO_H