imathmatrixalgo.h

image converter source code

///////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
// Digital Ltd. LLC
// 
// All rights reserved.
// 
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
// *       Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
// *       Redistributions in binary form must reproduce the above
// copyright notice, this list of conditions and the following disclaimer
// in the documentation and/or other materials provided with the
// distribution.
// *       Neither the name of Industrial Light & Magic nor the names of
// its contributors may be used to endorse or promote products derived
// from this software without specific prior written permission. 
// 
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
///////////////////////////////////////////////////////////////////////////


#ifndef INCLUDED_IMATHMATRIXALGO_H
#define INCLUDED_IMATHMATRIXALGO_H

//-------------------------------------------------------------------------
//
//      This file contains algorithms applied to or in conjunction with
//	transformation matrices (Imath::Matrix33 and Imath::Matrix44).
//	The assumption made is that these functions are called much less
//	often than the basic point functions or these functions require
//	more support classes.
//
//	This file also defines a few predefined constant matrices.
//
//-------------------------------------------------------------------------

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


#ifdef OPENEXR_DLL
    #ifdef IMATH_EXPORTS
        #define IMATH_EXPORT_CONST extern __declspec(dllexport)
    #else
	#define IMATH_EXPORT_CONST extern __declspec(dllimport)
    #endif
#else
    #define IMATH_EXPORT_CONST extern const
#endif


namespace Imath {

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

IMATH_EXPORT_CONST M33f identity33f;
IMATH_EXPORT_CONST M44f identity44f;
IMATH_EXPORT_CONST M33d identity33d;
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 an Imath::ZeroScale 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.
//

template <class T>  bool        extractScaling 
                                            (const Matrix44<T> &mat,
					     Vec3<T> &scl,
					     bool exc = true);
  
template <class T>  Matrix44<T> sansScaling (const Matrix44<T> &mat, 
					     bool exc = true);

template <class T>  bool        removeScaling 
                                            (Matrix44<T> &mat, 
					     bool exc = true);

template <class T>  bool        extractScalingAndShear 
                                            (const Matrix44<T> &mat,
					     Vec3<T> &scl,
					     Vec3<T> &shr,
					     bool exc = true);
  
template <class T>  Matrix44<T> sansScalingAndShear 
                                            (const Matrix44<T> &mat, 
					     bool exc = true);

template <class T>  bool        removeScalingAndShear 
                                            (Matrix44<T> &mat,
					     bool exc = true);

template <class T>  bool        extractAndRemoveScalingAndShear
                                            (Matrix44<T> &mat,
					     Vec3<T>     &scl,
					     Vec3<T>     &shr,
					     bool exc = true);

template <class T>  void	extractEulerXYZ 
                                            (const Matrix44<T> &mat,
					     Vec3<T> &rot);

template <class T>  void	extractEulerZYX 
                                            (const Matrix44<T> &mat,
					     Vec3<T> &rot);

template <class T>  Quat<T>	extractQuat (const Matrix44<T> &mat);

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);

template <class T>  bool	extractSHRT 
                                    (const Matrix44<T> &mat,
				     Vec3<T> &s,
				     Vec3<T> &h,
				     Vec3<T> &r,
				     Vec3<T> &t,
				     bool exc = true);

template <class T>  bool	extractSHRT 
                                    (const Matrix44<T> &mat,
				     Vec3<T> &s,
				     Vec3<T> &h,
				     Euler<T> &r,
				     Vec3<T> &t,
				     bool exc = true);

//
// Internal utility function.
//

template <class T>  bool	checkForZeroScaleInRow
                                            (const T       &scl, 
					     const Vec3<T> &row,
					     bool exc = true);

//
// Returns a matrix that rotates "fromDirection" vector to "toDirection"
// vector.
//

template <class T> Matrix44<T>	rotationMatrix (const Vec3<T> &fromDirection,
						const Vec3<T> &toDirection);



//
// Returns a matrix that rotates the "fromDir" vector 
// so that it points towards "toDir".  You may also 
// specify that you want the up vector to be pointing 
// in a certain direction "upDir".
//

template <class T> Matrix44<T>	rotationMatrixWithUpDir 
                                            (const Vec3<T> &fromDir,
					     const Vec3<T> &toDir,
					     const Vec3<T> &upDir);


//
// Returns a 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
//

template <class T> Matrix44<T>	alignZAxisWithTargetDir 
                                            (Vec3<T> targetDir, 
					     Vec3<T> upDir);


//----------------------------------------------------------------------


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

 
template <class T>  bool        extractScaling 
                                            (const Matrix33<T> &mat,
					     Vec2<T> &scl,
					     bool exc = true);
  
template <class T>  Matrix33<T> sansScaling (const Matrix33<T> &mat, 
					     bool exc = true);

template <class T>  bool        removeScaling 
                                            (Matrix33<T> &mat, 
					     bool exc = true);

template <class T>  bool        extractScalingAndShear 
                                            (const Matrix33<T> &mat,
					     Vec2<T> &scl,
					     T &h,
					     bool exc = true);
  
template <class T>  Matrix33<T> sansScalingAndShear 
                                            (const Matrix33<T> &mat, 
					     bool exc = true);

template <class T>  bool        removeScalingAndShear 
                                            (Matrix33<T> &mat,
					     bool exc = true);

template <class T>  bool        extractAndRemoveScalingAndShear
                                            (Matrix33<T> &mat,
					     Vec2<T>     &scl,
					     T           &shr,
					     bool exc = true);

template <class T>  void	extractEuler
                                            (const Matrix33<T> &mat,
					     T       &rot);

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

template <class T>  bool	checkForZeroScaleInRow
                                            (const T       &scl, 
					     const Vec2<T> &row,
					     bool exc = true);




//-----------------------------------------------------------------------------
// 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>
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::abs (row[i][j]) > maxVal)
		maxVal = Imath::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 = Math<T>::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;