imathmatrix.h

对gif

		     x[1][1] * a,
		     x[1][2] * a,
		     x[2][0] * a,
		     x[2][1] * a,
		     x[2][2] * a);
}

template <class T>
inline Matrix33<T>
operator * (T a, const Matrix33<T> &v)
{
    return v * a;
}

template <class T>
const Matrix33<T> &
Matrix33<T>::operator *= (const Matrix33<T> &v)
{
    Matrix33 tmp (T (0));

    for (int i = 0; i < 3; i++)
	for (int j = 0; j < 3; j++)
	    for (int k = 0; k < 3; k++)
		tmp.x[i][j] += x[i][k] * v.x[k][j];

    *this = tmp;
    return *this;
}

template <class T>
Matrix33<T>
Matrix33<T>::operator * (const Matrix33<T> &v) const
{
    Matrix33 tmp (T (0));

    for (int i = 0; i < 3; i++)
	for (int j = 0; j < 3; j++)
	    for (int k = 0; k < 3; k++)
		tmp.x[i][j] += x[i][k] * v.x[k][j];

    return tmp;
}

template <class T>
template <class S>
void
Matrix33<T>::multVecMatrix(const Vec2<S> &src, Vec2<S> &dst) const
{
    S a, b, w;

    a = src[0] * x[0][0] + src[1] * x[1][0] + x[2][0];
    b = src[0] * x[0][1] + src[1] * x[1][1] + x[2][1];
    w = src[0] * x[0][2] + src[1] * x[1][2] + x[2][2];

    dst.x = a / w;
    dst.y = b / w;
}

template <class T>
template <class S>
void
Matrix33<T>::multDirMatrix(const Vec2<S> &src, Vec2<S> &dst) const
{
    S a, b;

    a = src[0] * x[0][0] + src[1] * x[1][0];
    b = src[0] * x[0][1] + src[1] * x[1][1];

    dst.x = a;
    dst.y = b;
}

template <class T>
const Matrix33<T> &
Matrix33<T>::operator /= (T a)
{
    x[0][0] /= a;
    x[0][1] /= a;
    x[0][2] /= a;
    x[1][0] /= a;
    x[1][1] /= a;
    x[1][2] /= a;
    x[2][0] /= a;
    x[2][1] /= a;
    x[2][2] /= a;
  
    return *this;
}

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

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

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

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

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

    // Forward elimination

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

	T pivotsize = t[i][i];

	if (pivotsize < 0)
	    pivotsize = -pivotsize;

	for (j = i + 1; j < 3; 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 Matrix33();
	}

	if (pivot != i)
	{
	    for (j = 0; j < 3; 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 < 3; j++)
	{
	    T f = t[j][i] / t[i][i];

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

    // Backward substitution

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

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

	    return Matrix33();
	}

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

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

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

    return s;
}

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

template <class T>
Matrix33<T>
Matrix33<T>::inverse (bool singExc) const throw (Iex::MathExc)
{
    if (x[0][2] != 0 || x[1][2] != 0 || x[2][2] != 1)
    {
	Matrix33 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],

		    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],

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

	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 Matrix33();
		    }
		}
	    }
	}

	return s;
    }
    else
    {
	Matrix33 s ( x[1][1],
		    -x[0][1],
		     0, 

		    -x[1][0],
		     x[0][0],
		     0,

		     0,
		     0,
		     1);

	T r = x[0][0] * x[1][1] - x[1][0] * x[0][1];

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

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

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

	return s;
    }
}

template <class T>
template <class S>
const Matrix33<T> &
Matrix33<T>::setRotation (S r)
{
    S cos_r, sin_r;

    cos_r = Math<T>::cos (r);
    sin_r = Math<T>::sin (r);

    x[0][0] =  cos_r;
    x[0][1] =  sin_r;
    x[0][2] =  0;

    x[1][0] =  -sin_r;
    x[1][1] =  cos_r;
    x[1][2] =  0;

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

    return *this;
}

template <class T>
template <class S>
const Matrix33<T> &
Matrix33<T>::rotate (S r)
{
    *this *= Matrix33<T>().setRotation (r);
    return *this;
}

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

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

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

    return *this;
}

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

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

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

    return *this;
}

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

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

    return *this;
}

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

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

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

    return *this;
}

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

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

    return *this;
}

template <class T>
template <class S>
const Matrix33<T> &
Matrix33<T>::setShear (const S &xy)
{
    x[0][0] = 1;
    x[0][1] = 0;
    x[0][2] = 0;

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

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

    return *this;
}

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

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

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

    return *this;
}

template <class T>
template <class S>
const Matrix33<T> &
Matrix33<T>::shear (const S &xy)
{
    //
    // 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.
    // 

    x[1][0] += xy * x[0][0];
    x[1][1] += xy * x[0][1];
    x[1][2] += xy * x[0][2];

    return *this;
}

template <class T>
template <class S>
const Matrix33<T> &
Matrix33<T>::shear (const Vec2<S> &h)
{
    Matrix33<T> P (*this);
    
    x[0][0] = P[0][0] + h[1] * P[1][0];
    x[0][1] = P[0][1] + h[1] * P[1][1];
    x[0][2] = P[0][2] + h[1] * P[1][2];
    
    x[1][0] = P[1][0] + h[0] * P[0][0];
    x[1][1] = P[1][1] + h[0] * P[0][1];
    x[1][2] = P[1][2] + h[0] * P[0][2];

    return *this;
}


//---------------------------
// Implementation of Matrix44
//---------------------------

template <class T>
inline T *
Matrix44<T>::operator [] (int i)
{
    return x[i];
}

template <class T>
inline const T *
Matrix44<T>::operator [] (int i) const