imathmatrix.h

image converter source code

    x[2][3] /= a;
    x[3][0] /= a;
    x[3][1] /= a;
    x[3][2] /= a;
    x[3][3] /= a;

    return *this;
}

template <class T>
Matrix44<T>
Matrix44<T>::operator / (T a) const
{
    return Matrix44 (x[0][0] / a,
		     x[0][1] / a,
		     x[0][2] / a,
		     x[0][3] / a,
		     x[1][0] / a,
		     x[1][1] / a,
		     x[1][2] / a,
		     x[1][3] / a,
		     x[2][0] / a,
		     x[2][1] / a,
		     x[2][2] / a,
		     x[2][3] / a,
		     x[3][0] / a,
		     x[3][1] / a,
		     x[3][2] / a,
		     x[3][3] / a);
}

template <class T>
const Matrix44<T> &
Matrix44<T>::transpose ()
{
    Matrix44 tmp (x[0][0],
		  x[1][0],
		  x[2][0],
		  x[3][0],
		  x[0][1],
		  x[1][1],
		  x[2][1],
		  x[3][1],
		  x[0][2],
		  x[1][2],
		  x[2][2],
		  x[3][2],
		  x[0][3],
		  x[1][3],
		  x[2][3],
		  x[3][3]);
    *this = tmp;
    return *this;
}

template <class T>
Matrix44<T>
Matrix44<T>::transposed () const
{
    return Matrix44 (x[0][0],
		     x[1][0],
		     x[2][0],
		     x[3][0],
		     x[0][1],
		     x[1][1],
		     x[2][1],
		     x[3][1],
		     x[0][2],
		     x[1][2],
		     x[2][2],
		     x[3][2],
		     x[0][3],
		     x[1][3],
		     x[2][3],
		     x[3][3]);
}

template <class T>
const Matrix44<T> &
Matrix44<T>::gjInvert (bool singExc) throw (Iex::MathExc)
{
    *this = gjInverse (singExc);
    return *this;
}

template <class T>
Matrix44<T>
Matrix44<T>::gjInverse (bool singExc) const throw (Iex::MathExc)
{
    int i, j, k;
    Matrix44 s;
    Matrix44 t (*this);

    // Forward elimination

    for (i = 0; i < 3 ; i++)
    {
	int pivot = i;

	T pivotsize = t[i][i];

	if (pivotsize < 0)
	    pivotsize = -pivotsize;

	for (j = i + 1; j < 4; j++)
	{
	    T tmp = t[j][i];

	    if (tmp < 0)
		tmp = -tmp;

	    if (tmp > pivotsize)
	    {
		pivot = j;
		pivotsize = tmp;
	    }
	}

	if (pivotsize == 0)
	{
	    if (singExc)
		throw ::Imath::SingMatrixExc ("Cannot invert singular matrix.");

	    return Matrix44();
	}

	if (pivot != i)
	{
	    for (j = 0; j < 4; j++)
	    {
		T tmp;

		tmp = t[i][j];
		t[i][j] = t[pivot][j];
		t[pivot][j] = tmp;

		tmp = s[i][j];
		s[i][j] = s[pivot][j];
		s[pivot][j] = tmp;
	    }
	}

	for (j = i + 1; j < 4; j++)
	{
	    T f = t[j][i] / t[i][i];

	    for (k = 0; k < 4; k++)
	    {
		t[j][k] -= f * t[i][k];
		s[j][k] -= f * s[i][k];
	    }
	}
    }

    // Backward substitution

    for (i = 3; i >= 0; --i)
    {
	T f;

	if ((f = t[i][i]) == 0)
	{
	    if (singExc)
		throw ::Imath::SingMatrixExc ("Cannot invert singular matrix.");

	    return Matrix44();
	}

	for (j = 0; j < 4; j++)
	{
	    t[i][j] /= f;
	    s[i][j] /= f;
	}

	for (j = 0; j < i; j++)
	{
	    f = t[j][i];

	    for (k = 0; k < 4; k++)
	    {
		t[j][k] -= f * t[i][k];
		s[j][k] -= f * s[i][k];
	    }
	}
    }

    return s;
}

template <class T>
const Matrix44<T> &
Matrix44<T>::invert (bool singExc) throw (Iex::MathExc)
{
    *this = inverse (singExc);
    return *this;
}

template <class T>
Matrix44<T>
Matrix44<T>::inverse (bool singExc) const throw (Iex::MathExc)
{
    if (x[0][3] != 0 || x[1][3] != 0 || x[2][3] != 0 || x[3][3] != 1)
	return gjInverse(singExc);

    Matrix44 s (x[1][1] * x[2][2] - x[2][1] * x[1][2],
		x[2][1] * x[0][2] - x[0][1] * x[2][2],
		x[0][1] * x[1][2] - x[1][1] * x[0][2],
		0,

		x[2][0] * x[1][2] - x[1][0] * x[2][2],
		x[0][0] * x[2][2] - x[2][0] * x[0][2],
		x[1][0] * x[0][2] - x[0][0] * x[1][2],
		0,

		x[1][0] * x[2][1] - x[2][0] * x[1][1],
		x[2][0] * x[0][1] - x[0][0] * x[2][1],
		x[0][0] * x[1][1] - x[1][0] * x[0][1],
		0,

		0,
		0,
		0,
		1);

    T r = x[0][0] * s[0][0] + x[0][1] * s[1][0] + x[0][2] * s[2][0];

    if (Imath::abs (r) >= 1)
    {
	for (int i = 0; i < 3; ++i)
	{
	    for (int j = 0; j < 3; ++j)
	    {
		s[i][j] /= r;
	    }
	}
    }
    else
    {
	T mr = Imath::abs (r) / limits<T>::smallest();

	for (int i = 0; i < 3; ++i)
	{
	    for (int j = 0; j < 3; ++j)
	    {
		if (mr > Imath::abs (s[i][j]))
		{
		    s[i][j] /= r;
		}
		else
		{
		    if (singExc)
			throw SingMatrixExc ("Cannot invert singular matrix.");

		    return Matrix44();
		}
	    }
	}
    }

    s[3][0] = -x[3][0] * s[0][0] - x[3][1] * s[1][0] - x[3][2] * s[2][0];
    s[3][1] = -x[3][0] * s[0][1] - x[3][1] * s[1][1] - x[3][2] * s[2][1];
    s[3][2] = -x[3][0] * s[0][2] - x[3][1] * s[1][2] - x[3][2] * s[2][2];

    return s;
}

template <class T>
template <class S>
const Matrix44<T> &
Matrix44<T>::setEulerAngles (const Vec3<S>& r)
{
    S cos_rz, sin_rz, cos_ry, sin_ry, cos_rx, sin_rx;
    
    cos_rz = Math<T>::cos (r[2]);
    cos_ry = Math<T>::cos (r[1]);
    cos_rx = Math<T>::cos (r[0]);
    
    sin_rz = Math<T>::sin (r[2]);
    sin_ry = Math<T>::sin (r[1]);
    sin_rx = Math<T>::sin (r[0]);
    
    x[0][0] =  cos_rz * cos_ry;
    x[0][1] =  sin_rz * cos_ry;
    x[0][2] = -sin_ry;
    x[0][3] =  0;
    
    x[1][0] = -sin_rz * cos_rx + cos_rz * sin_ry * sin_rx;
    x[1][1] =  cos_rz * cos_rx + sin_rz * sin_ry * sin_rx;
    x[1][2] =  cos_ry * sin_rx;
    x[1][3] =  0;
    
    x[2][0] =  sin_rz * sin_rx + cos_rz * sin_ry * cos_rx;
    x[2][1] = -cos_rz * sin_rx + sin_rz * sin_ry * cos_rx;
    x[2][2] =  cos_ry * cos_rx;
    x[2][3] =  0;

    x[3][0] =  0;
    x[3][1] =  0;
    x[3][2] =  0;
    x[3][3] =  1;

    return *this;
}

template <class T>
template <class S>
const Matrix44<T> &
Matrix44<T>::setAxisAngle (const Vec3<S>& axis, S angle)
{
    Vec3<S> unit (axis.normalized());
    S sine   = Math<T>::sin (angle);
    S cosine = Math<T>::cos (angle);

    x[0][0] = unit[0] * unit[0] * (1 - cosine) + cosine;
    x[0][1] = unit[0] * unit[1] * (1 - cosine) + unit[2] * sine;
    x[0][2] = unit[0] * unit[2] * (1 - cosine) - unit[1] * sine;
    x[0][3] = 0;

    x[1][0] = unit[0] * unit[1] * (1 - cosine) - unit[2] * sine;
    x[1][1] = unit[1] * unit[1] * (1 - cosine) + cosine;
    x[1][2] = unit[1] * unit[2] * (1 - cosine) + unit[0] * sine;
    x[1][3] = 0;

    x[2][0] = unit[0] * unit[2] * (1 - cosine) + unit[1] * sine;
    x[2][1] = unit[1] * unit[2] * (1 - cosine) - unit[0] * sine;
    x[2][2] = unit[2] * unit[2] * (1 - cosine) + cosine;
    x[2][3] = 0;

    x[3][0] = 0;
    x[3][1] = 0;
    x[3][2] = 0;
    x[3][3] = 1;

    return *this;
}

template <class T>
template <class S>
const Matrix44<T> &
Matrix44<T>::rotate (const Vec3<S> &r)
{
    S cos_rz, sin_rz, cos_ry, sin_ry, cos_rx, sin_rx;
    S m00, m01, m02;
    S m10, m11, m12;
    S m20, m21, m22;

    cos_rz = Math<S>::cos (r[2]);
    cos_ry = Math<S>::cos (r[1]);
    cos_rx = Math<S>::cos (r[0]);
    
    sin_rz = Math<S>::sin (r[2]);
    sin_ry = Math<S>::sin (r[1]);
    sin_rx = Math<S>::sin (r[0]);

    m00 =  cos_rz *  cos_ry;
    m01 =  sin_rz *  cos_ry;
    m02 = -sin_ry;
    m10 = -sin_rz *  cos_rx + cos_rz * sin_ry * sin_rx;
    m11 =  cos_rz *  cos_rx + sin_rz * sin_ry * sin_rx;
    m12 =  cos_ry *  sin_rx;
    m20 = -sin_rz * -sin_rx + cos_rz * sin_ry * cos_rx;
    m21 =  cos_rz * -sin_rx + sin_rz * sin_ry * cos_rx;
    m22 =  cos_ry *  cos_rx;

    Matrix44<T> P (*this);

    x[0][0] = P[0][0] * m00 + P[1][0] * m01 + P[2][0] * m02;
    x[0][1] = P[0][1] * m00 + P[1][1] * m01 + P[2][1] * m02;
    x[0][2] = P[0][2] * m00 + P[1][2] * m01 + P[2][2] * m02;
    x[0][3] = P[0][3] * m00 + P[1][3] * m01 + P[2][3] * m02;

    x[1][0] = P[0][0] * m10 + P[1][0] * m11 + P[2][0] * m12;
    x[1][1] = P[0][1] * m10 + P[1][1] * m11 + P[2][1] * m12;
    x[1][2] = P[0][2] * m10 + P[1][2] * m11 + P[2][2] * m12;
    x[1][3] = P[0][3] * m10 + P[1][3] * m11 + P[2][3] * m12;

    x[2][0] = P[0][0] * m20 + P[1][0] * m21 + P[2][0] * m22;
    x[2][1] = P[0][1] * m20 + P[1][1] * m21 + P[2][1] * m22;
    x[2][2] = P[0][2] * m20 + P[1][2] * m21 + P[2][2] * m22;
    x[2][3] = P[0][3] * m20 + P[1][3] * m21 + P[2][3] * m22;

    return *this;
}

template <class T>
const Matrix44<T> &
Matrix44<T>::setScale (T s)
{
    x[0][0] = s;
    x[0][1] = 0;
    x[0][2] = 0;
    x[0][3] = 0;

    x[1][0] = 0;
    x[1][1] = s;
    x[1][2] = 0;
    x[1][3] = 0;

    x[2][0] = 0;
    x[2][1] = 0;
    x[2][2] = s;
    x[2][3] = 0;

    x[3][0] = 0;
    x[3][1] = 0;
    x[3][2] = 0;
    x[3][3] = 1;

    return *this;
}

template <class T>
template <class S>
const Matrix44<T> &
Matrix44<T>::setScale (const Vec3<S> &s)
{
    x[0][0] = s[0];
    x[0][1] = 0;
    x[0][2] = 0;
    x[0][3] = 0;

    x[1][0] = 0;
    x[1][1] = s[1];
    x[1][2] = 0;
    x[1][3] = 0;

    x[2][0] = 0;
    x[2][1] = 0;
    x[2][2] = s[2];
    x[2][3] = 0;

    x[3][0] = 0;
    x[3][1] = 0;
    x[3][2] = 0;
    x[3][3] = 1;

    return *this;
}

template <class T>
template <class S>
const Matrix44<T> &
Matrix44<T>::scale (const Vec3<S> &s)
{
    x[0][0] *= s[0];
    x[0][1] *= s[0];
    x[0][2] *= s[0];
    x[0][3] *= s[0];

    x[1][0] *= s[1];
    x[1][1] *= s[1];
    x[1][2] *= s[1];
    x[1][3] *= s[1];

    x[2][0] *= s[2];
    x[2][1] *= s[2];
    x[2][2] *= s[2];
    x[2][3] *= s[2];

    return *this;
}

template <class T>
template <class S>
const Matrix44<T> &
Matrix44<T>::setTranslation (const Vec3<S> &t)
{
    x[0][0] = 1;
    x[0][1] = 0;
    x[0][2] = 0;
    x[0][3] = 0;

    x[1][0] = 0;
    x[1][1] = 1;
    x[1][2] = 0;
    x[1][3] = 0;

    x[2][0] = 0;
    x[2][1] = 0;
    x[2][2] = 1;
    x[2][3] = 0;

    x[3][0] = t[0];
    x[3][1] = t[1];
    x[3][2] = t[2];
    x[3][3] = 1;

    return *this;
}

template <class T>
inline const Vec3<T>
Matrix44<T>::translation () const
{
    return Vec3<T> (x[3][0], x[3][1], x[3][2]);
}

template <class T>
template <class S>
const Matrix44<T> &
Matrix44<T>::translate (const Vec3<S> &t)
{
    x[3][0] += t[0] * x[0][0] + t[1] * x[1][0] + t[2] * x[2][0];
    x[3][1] += t[0] * x[0][1] + t[1] * x[1][1] + t[2] * x[2][1];
    x[3][2] += t[0] * x[0][2] + t[1] * x[1][2] + t[2] * x[2][2];
    x[3][3] += t[0] * x[0][3] + t[1] * x[1][3] + t[2] * x[2][3];

    return *this;
}

template <class T>
template <class S>
const Matrix44<T> &
Matrix44<T>::setShear (const Vec3<S> &h)
{
    x[0][0] = 1;
    x[0][1] = 0;
    x[0][2] = 0;
    x[0][3] = 0;

    x[1][0] = h[0];
    x[1][1] = 1;
    x[1][2] = 0;
    x[1][3] = 0;

    x[2][0] = h[1];
    x[2][1] = h[2];
    x[2][2] = 1;
    x[2][3] = 0;

    x[3][0] = 0;
    x[3][1] = 0;
    x[3][2] = 0;
    x[3][3] = 1;

    return *this;
}

template <class T>
template <class S>
const Matrix44<T> &
Matrix44<T>::setShear (const Shear6<S> &h)
{
    x[0][0] = 1;
    x[0][1] = h.yx;
    x[0][2] = h.zx;
    x[0][3] = 0;

    x[1][0] = h.xy;
    x[1][1] = 1;
    x[1][2] = h.zy;
    x[1][3] = 0;

    x[2][0] = h.xz;
    x[2][1] = h.yz;
    x[2][2] = 1;
    x[2][3] = 0;

    x[3][0] = 0;
    x[3][1] = 0;
    x[3][2] = 0;
    x[3][3] = 1;

    return *this;
}

template <class T>
template <class S>
const Matrix44<T> &
Matrix44<T>::shear (const Vec3<S> &h)
{
    //
    // In this case, we don't need a temp. copy of the matrix 
    // because we never use a value on the RHS after we've 
    // changed it on the LHS.
    // 

    for (int i=0; i < 4; i++)
    {
	x[2][i] += h[1] * x[0][i] + h[2] * x[1][i];
	x[1][i] += h[0] * x[0][i];
    }

    return *this;
}

template <class T>
template <class S>
const Matrix44<T> &
Matrix44<T>::shear (const Shear6<S> &h)
{
    Matrix44<T> P (*this);

    for (int i=0; i < 4; i++)
    {
	x[0][i] = P[0][i] + h.yx * P[1][i] + h.zx * P[2][i];
	x[1][i] = h.xy * P[0][i] + P[1][i] + h.zy * P[2][i];
	x[2][i] = h.xz * P[0][i] + h.yz * P[1][i] + P[2][i];
    }

    return *this;
}


//--------------------------------
// Implementation of stream output
//--------------------------------

template <class T>
std::ostream &
operator << (std::ostream &s, const Matrix33<T> &m)
{
    std::ios_base::fmtflags oldFlags = s.flags();
    int width;