矩阵类模板

        这是我在做游戏三维图形处理中编制的一个矩阵类模板，可以产生整型，浮点，和自定义类型的任意阶矩阵和行列式，并且定义了一些常用的三维向量运算函数模板，希望大家给出修改意见。建议大家只看函数生命，不要看具体实现，因为贴上来之后格式不对头，由于文件比较长，我将其分为两个帖子

                                                  矩阵类模板

//***************************************************************************/
//声明：                                                                    */
//                                                                          */
//这个版本纯粹是处于个人爱好而编写，其中可能有不少错误。主要是为了让各位爱好*/
//者相互交流共同提高编程水平。因此，你可以任意修改复制其中的代码，并可以应用*/
//于任何场合，但由此带来的问题本人不负任何责任。如果你有什么建议的话，欢迎与*/
//我联系，我们共同探讨。                                                    */
//我的邮箱是：[e-mail]                                          */
//***************************************************************************/
//                                                                          */
//Matrix.h                                                                  */
//本文件定义了基本的矩阵模板，可以由此复制出各种int,float,double等          */
//等各种类型的矩阵类型。重载了基本的矩阵运算符。为了满足矢量处理的          */
//要求，定义了一系列坐标变换函数                                            */
//                                                                          */
//                                                                          */
//***************************************************************************/
//版本升级说明：原来的版本中矩阵这一单一类型可以满足一般的需要，但是当在结构*/
//没有默认的构造函数，而且在编写坐标变换函数模板时非常不方便，且无法定义各种*/
//常用的特殊矢量。归根结蒂是因为原来的模板中只有一个矩阵内数据类型的参数。在*/
//这个版本中将模板参数扩展为3个，有两个是矩阵的行和列，取消了构造时的init参 */
//数，统一初始化为0。这样矩阵就有了默认的构造函数，且模板参数中有足够的信息 */
//有效的解决了上面提出的问题。                                              */
//***************************************************************************/
//编写人：王云                                                              */
//修改人：王云
//开始时间：2002.10.18                                                      */
//结束时间：2002.10.19                                                      */
//***************************************************************************/
#include<windows.h>
#include<iostream.h>
#include<math.h>
template<typename ElemType,int r,int c> class Matrix;
//***************************************************************************/
//二维向量和点                                                              */
//***************************************************************************/
typedef Matrix<int,1,2>      Vector2i;
typedef Matrix<int,1,2>      Vertex2i;
typedef Matrix<float,1,2>    Vector2f;
typedef Matrix<float,1,2>    Vertex2f;
typedef Matrix<double,1,2>   Vector2d;
typedef Matrix<double,1,2>   Vertex2d;
//***************************************************************************/
//三维向量和点                                                              */
//***************************************************************************/
typedef Matrix<int,1,3>      Vector3i;
typedef Matrix<int,1,3>      Vertex3i;
typedef Matrix<float,1,3>    Vector3f;
typedef Matrix<float,1,3>    Vertex3f;
typedef Matrix<double,1,3>   Vector3d;
typedef Matrix<double,1,3>   Vertex3d;
//***************************************************************************/
//四维齐次坐标和向量                                                        */
//***************************************************************************/
typedef Matrix<int,1,4>      Vector4i;
typedef Matrix<int,1,4>      Vertex4i;
typedef Matrix<float,1,4>    Vector4f;
typedef Matrix<float,1,4>    Vertex4f;
typedef Matrix<double,1,4>   Vector4d;
typedef Matrix<double,1,4>   Vertex4d;
//***************************************************************************/
//常用的向量运算函数模板定义                                                */
//***************************************************************************/
/*
template<typename ElemType,int r,int c>                    //向量在各个局部坐标轴上的分量
int VectorProjection(Matrix<ElemType,r,c>& vector,Matrix<ElemType,r,c>*pVectorX,Matrix<ElemType,r,c>*pVectorY,Matrix<ElemType,r,c>* pVectorZ);
*/
template<typename ElemType,int r,int c>                    //计算向量的模
ElemType VectorMo(Matrix<ElemType,r,c>& vector);
template<typename ElemType,int r,int c>                    //向量单位化
Matrix<ElemType,r,c> VectorUnit(Matrix<ElemType,r,c>& vector);
template<typename ElemType,int r,int c>                    //绕x轴旋转
Matrix<ElemType,r,c> RotateOX(Matrix<ElemType,r,c>& vector,float beta);
template<typename ElemType,int r,int c>                    //绕x轴旋转
Matrix<ElemType,r,c> RotateOY(Matrix<ElemType,r,c>& vector,float beta);
template<typename ElemType,int r,int c>                    //绕y轴旋转
Matrix<ElemType,r,c> RotateOZ(Matrix<ElemType,r,c>& vector,float beta);
template<typename ElemType,int r,int c>                    //绕z轴旋转
Matrix<ElemType,r,c>  RotateON(Matrix<ElemType,r,c>& reel,Matrix<ElemType,r,c>& vector,float beta);
template<typename ElemType,int r,int c>                    //绕任意轴旋转
Matrix<ElemType,r,c> VectorProduct(Matrix<ElemType,r,c>& va,Matrix<ElemType,r,c>& vb);             //向量积
template<typename ElemType,int r,int c>                    //坐标的变换
int CoordinateTransfer(Matrix<ElemType,r,c>& origin,Matrix<ElemType,r,c>& xDirect,Matrix<ElemType,r,c>& yDirect,Matrix<ElemType,r,c>& zDirect,
        Matrix<ElemType,r,c>*& pVertex,int nVertex);    
template<typename ElemType,int r,int c>                     //两个向量的夹角
int VectorAngle(Matrix<ElemType,r,c>& b,Matrix<ElemType,r,c>& a,float* sin,float* cos);
//***************************************************************************/
//矩阵类模板的定义                                                          */
//                                                                          */
//                                                                          */
//***************************************************************************/
template<typename ElemType,int r,int c> class Matrix
{
private:
 int row;                                               //行
 int col;                                               //列
 ElemType* pElemType;                                   //存储矩阵数据的类型数组

public:
 Matrix();
 Matrix(const Matrix<ElemType,r,c>&A);                  //深度堆拷贝函数必须在所有的运算函数之
                                                        //前定义，否则其它函数的返回给参数传递
                                                        //将采用默认的拷贝构造函数，从而得到错
                                                        //误的结果
 ~Matrix();

public:
 inline int GetRow() {return row;}
 inline int GetCol() {return col;}

    inline ElemType& operator ()(const int& i,const int& j) {return pElemType[(i-1)*row+j-1];}
 Matrix<ElemType,r,c>& operator = (const Matrix<ElemType,r,c>&A);

 friend Matrix<ElemType,r,c> operator + (const Matrix<ElemType,r,c>& A,const Matrix<ElemType,r,c>& B);
 friend Matrix<ElemType,r,c> operator - (const Matrix<ElemType,r,c>& A,const Matrix<ElemType,r,c>& B);
 friend Matrix<ElemType,r,c> operator * (const Matrix<ElemType,r,c>& A,const Matrix<ElemType,r,c>& B);

                                                        //下面这三个函数是针对常数和矩阵的乘法而
                                                        //定义的，显然，冗余非常大，未来的目标是
                                                        //在内部定一个函数模板，来自动匹配类型
    friend Matrix<ElemType,r,c> operator * (const double& scale,const Matrix<ElemType,r,c>& B);
 friend Matrix<ElemType,r,c> operator * (const float& scale,const Matrix<ElemType,r,c>& B);
 friend Matrix<ElemType,r,c> operator * (const int& scale,const Matrix<ElemType,r,c>& B);
   
                                                           //＋＝操作符之所以要返回对A的常量引用，
                                                        //是为了能够使用连等式，即C=A=D+E ,同
                                                        //时，返回值不能作为左值
 friend const Matrix<ElemType,r,c>& operator +=(Matrix<ElemType,r,c>& A,const Matrix<ElemType,r,c>& B);
 friend const Matrix<ElemType,r,c>& operator -=(Matrix<ElemType,r,c>& A,const Matrix<ElemType,r,c>& B);
 friend const Matrix<ElemType,r,c>& operator *=(Matrix<ElemType,r,c>& A,const Matrix<ElemType,r,c>& B);

                                                           //下面两个操作符的重载主要是为了字符界
                                                        //面下的使用和测试
 friend istream& operator >> (istream& in,Matrix<ElemType,r,c>& A);
 friend ostream& operator << (ostream& out,Matrix<ElemType,r,c>& A);

};

 

 

 

 

//***************************************************************************/
//常用的向量运算函数模板实现代码                                            */
//***************************************************************************/
template<typename ElemType,int r,int c>
ElemType VectorMo(Matrix<ElemType,r,c>& vector)
{
 ElemType a=vector(1,1),
       b=vector(1,2),
    c=vector(1,3);
 return sqrt(a*a+b*b+c*c);
}
template<typename ElemType,int r,int c>                    //向量单位化
Matrix<ElemType,r,c> VectorUnit(Matrix<ElemType,r,c>& vector)
{
 /*
 Matrix<ElemType,r,c> result;
 ElemType mo=1/VectorMo(vector);
 result=mo*vector;
    return result;
 */
 return 1/VectorMo(vector)*vector;
}
template<typename ElemType,int r,int c>  
Matrix<ElemType,r,c> RotateOX(Matrix<ElemType,r,c>& vector,float beta)

{
 Matrix<ElemType,r,c> result;
                                                           //为了减少下面调用函数的次数,增加空间暂
                                                        //时存放旋转轴的各个分量
 ElemType x=vector(1,1);
 ElemType y=vector(1,2);
 ElemType z=vector(1,3);
 
 result(1,1)=x;
 result(1,2)=y*cos(beta)-z*sin(beta);
 result(1,3)=y*sin(beta)+z*cos(beta);

    return result;
}
template<typename ElemType,int r,int c>  
Matrix<ElemType,r,c> RotateOY(Matrix<ElemType,r,c>& vector,float beta)
{
 Matrix<ElemType,r,c> result;
                                                           //为了减少下面调用函数的次数,增加空间暂
                                                        //时存放旋转轴的各个分量
 ElemType x=vector(1,1);
 ElemType y=vector(1,2);
 ElemType z=vector(1,3);

 result(1,1)=x*cos(beta)+z*sin(beta);
 result(1,2)=y;
 result(1,3)=-x*sin(beta)+z*cos(beta);

 return result;
}
template<typename ElemType,int r,int c>
Matrix<ElemType,r,c> RotateOZ(Matrix<ElemType,r,c>& vector,float beta)
{
 Matrix<ElemType,r,c> result;
                                                           //为了减少下面调用函数的次数,增加空间暂
                                                        //时存放旋转轴的各个分量
 ElemType x=vector(1,1);
 ElemType y=vector(1,2);
 ElemType z=vector(1,3);

 result(1,1)=x*cos(beta)-y*sin(beta);
 result(1,2)=x*sin(beta)+y*cos(beta);
 result(1,3)=z;

 return result;
}
template<typename ElemType,int r,int c>
Matrix<ElemType,r,c>  RotateON(Matrix<ElemType,r,c>& reel,Matrix<ElemType,r,c>& vector,float beta)
{
 Matrix<ElemType,r,c> result;
                                                           //为了减少下面调用函数的次数,增加空间暂
                                                        //时存放旋转轴的各个分量
 ElemType reelX=reel(1,1);
 ElemType reelY=reel(1,2);
 ElemType reelZ=reel(1,3);
                                                           //为了减少频繁的大量运算，增加空间暂时
                                                        //存放过原点的旋转轴与x,y,z轴的夹角的余
                                                        //弦值
 ElemType s=sqrt( reelX*reelX + reelY*reelY + reelZ*reelZ );
 ElemType n1=reelX/s;
 ElemType n2=reelY/s;
 ElemType n3=reelZ/s;
                                                           //为了减少计算次数，暂时存放常用的数据
 ElemType n1n1=n1*n1,
          n2n2=n2*n2,
    n3n3=n3*n3,
    n1n3=n1*n3,
    n1n2=n1*n2,
    n2n3=n2*n3;
                                                        //为了减少计算次数，暂时存放旋转角度的
                                                        //正弦值和余弦值
 ElemType cosBeta=cos(beta),
       sinBeta=sin(beta);
                                                           //核心，计算旋转后的各个分量
    result(1,1)=reelX*( n1n1+(1-n1n1)*cosBeta )+reelY*( n1n2*(1-cosBeta)-n3*sinBeta )+reelZ*( n1n3*(1-cosBeta)+n2*sinBeta );
 result(1,2)=reelX*( n1n2*(1-cosBeta)+n3*sinBeta )+reelY*( n2n2+(1-n2n2)*cosBeta )+reelZ*( n2n3*(1-cosBeta)-n1*sinBeta );
 result(1,3)=reelX*( n1n3*(1-cosBeta)-n2*sinBeta )+reelY*( n2n3*(1-cosBeta)+n1*sinBeta )+reelZ*( n3n3+(1-n3n3)*cosBeta );

 return result;
}
template<typename ElemType,int r,int c>
Matrix<ElemType,r,c> VectorProduct(Matrix<ElemType,r,c>& va,Matrix<ElemType,r,c>& vb)             //向量积
{
 Matrix<ElemType,r,c> result(1,3,0);
 result(1,1)=va(1,2)*vb(1,3)-va(1,3)*vb)*vb(1,2);
 result(1,2)=va(1,3)*vb(1,1)-va(1,1)*vb(1,3);
 result(1,3)=va(1,1)*vb(1,2)-va(1,2)*vb(1,1);
 return result;
}
template<typename ElemType,int r,int c>
int CoordinateTransfer(Matrix<ElemType,r,c>& origin,Matrix<ElemType,r,c>& xDirect,Matrix<ElemType,r,c>& yDirect,Matrix<ElemType,r,c>& zDirect,
        Matrix<ElemType,r,c>*& pVertex,int nVertex)
                                            //坐标变换
{
 return 0;
}
template<typename ElemType,int r,int c>
int VectorAngle(Matrix<ElemType,r,c>& b,Matrix<ElemType,r,c>& a,float* sin,float* cos)
{                                                          //求两个向量夹角的正弦和余弦
 ElemType ax=a(1,1),
       ay=a(1,2),
       az=a(1,3),
       bx=b(1,1),
       by=b(1,2),
       bz=b(1,3);
 *cos=(ax*bx+ay*by+az*bz)/sqrt(ax*ax+ay*ay+az*az)/sqrt(bx*bx+by*by+bz*bz);
 *sin=sqrt( 1-(*cos * *cos) );
 return 1;
}

 

//***************************************************************************/
//矩阵类模板的实现                                                          */
//***************************************************************************/
template<typename ElemType,int r,int c>
Matrix<ElemType,r,c>::Matrix()
{
 this->row=r;
 this->col=c;
 while(  !(this->pElemType=new ElemType[row*col])  );   //要构造一个矩阵，一定能够成功
 for(int i=0;i<row*col;i++)
  pElemType[i]=0;
}

template<typename ElemType,int r,int c>
Matrix<ElemType,r,c>::Matrix(const Matrix<ElemType,r,c>&A) //深度堆拷贝构造函数
{
 row=A.row;
 col=A.col;
 while(  !(pElemType=new ElemType[row*col])  );
 for(int i=0;i<row*col;i++)
  pElemType[i]=A.pElemType[i];
}

template<typename ElemType,int r,int c>
Matrix<ElemType,r,c>::~Matrix()
{
 if(pElemType)
  delete[]pElemType;
}

template<typename ElemType,int r,int c>
Matrix<ElemType,r,c>& Matrix<ElemType,r,c>::operator = (const Matrix<ElemType,r,c>&A)
{
                                                        //因该在这里插入矩阵行列检查的语句,下同
 row=A.row;
 col=A.col;
 //while(  !(pElemType=new ElemType[row*col])  );
                                                        //为提高效率，可以考虑采用内存区拷贝技术
 for(int i=0;i<row*col;i++)
  pElemType[i]=A.pElemType[i];
 cout<<"asdkfjkasdfjdklsfjasdklfj"<<endl;
 return *this;
}

template<typename ElemType,int r,int c>
Matrix<ElemType,r,c> operator + (const Matrix<ElemType,r,c>&A,const Matrix<ElemType,r,c>&B)
{
 Matrix<ElemType,r,c> result=A;                         //注意，这里使用const类型来构造对象的前
                                                        //提条件是堆拷贝函数的参数也必须是const
                                                        //类型，否则要发生错误
 for(int i=0;i<result.row*result.col;i++)
  //result.pElemType[i]=A.pElemType[i]+B.pElemType[i];
  result.pElemType[i]+=B.pElemType[i];
 return result;
}

template<typename ElemType,int r,int c>
Matrix<ElemType,r,c> operator - (const Matrix<ElemType,r,c>&A,const Matrix<ElemType,r,c>&B)
{
 Matrix<ElemType,r,c> result=A;
 for(int i=0;i<result.row*result.col;i++)
  result.pElemType[i]-=B.pElemType[i];
 return result;
}

template<typename ElemType,int r,int c>
Matrix<ElemType,r,c> operator * (const Matrix<ElemType,r,c>&A,const Matrix<ElemType,r,c>&B)
{
                                                        //对于矩阵的乘法，这里应该插入对两个操
                                                        //作数行列的检查
 Matrix<ElemType,r,c> result;
 for(int i=0;i<result.row;i++)
  for(int j=0;j<result.col;j++)
  {                                                  //这里定义了一个temp，避免了累加时，对
                                                  //矩阵数据数组下标的频繁引用，省去了频
                                                  //繁的查找，求出结果时，进行一次引用，
                                                  //较少了乘法的次数，提高了效率
   ElemType temp=0;
   for(int k=0;k<A.col;k++)
    temp+=A.pElemType[i*A.row+k]*B.pElemType[i*k+j];
   result.pElemType[i*result.row+j]=temp;
  }
 return result;
}

template<typename ElemType,int r,int c>
Matrix<ElemType,r,c> operator * (const double& scale,const Matrix<ElemType,r,c>& B)
{
 Matrix<ElemType,r,c> result=B;
 for(int i=0;i<B.row*B.col;i++)
  result.pElemType[i]*=(ElemType)scale;
 return result;
}

template<typename ElemType,int r,int c>
Matrix<ElemType,r,c> operator * (const float& scale,const Matrix<ElemType,r,c>& B)
{
 Matrix<ElemType,r,c> result=B;
 for(int i=0;i<B.row*B.col;i++)
  result.pElemType[i]*=(ElemType)scale;
 return result;
}

template<typename ElemType,int r,int c>
Matrix<ElemType,r,c> operator * (const int& scale,const Matrix<ElemType,r,c>& B)
{
 Matrix<ElemType,r,c> result=B;
 for(int i=0;i<B.row*B.col;i++)
  result.pElemType[i]*=(ElemType)scale;
 return result;
}

template<typename ElemType,int r,int c>
const Matrix<ElemType,r,c>& operator +=(Matrix<ElemType,r,c>& A,const Matrix<ElemType,r,c>& B)
{
 A=A+B;
 return A;
}

template<typename ElemType,int r,int c>
const Matrix<ElemType,r,c>& operator -=(Matrix<ElemType,r,c>& A,const Matrix<ElemType,r,c>& B)
{
 A=A-B;
 return A;
}

template<typename ElemType,int r,int c>
const Matrix<ElemType,r,c>& operator *=(Matrix<ElemType,r,c>& A,const Matrix<ElemType,r,c>& B)
{
 A=A*B;
 return A;
}
template<typename ElemType,int r,int c>
istream& operator >>(istream& in,Matrix<ElemType,r,c>& A)
{
 cout<<"please input the matrix:"<<endl;
 for(int i=0;i<A.row;i++)
  for(int j=0;j<A.col;j++)
   in>>A.pElemType[i*row+j]; 
 return in;
}
template<typename ElemType,int r,int c>
ostream& operator <<(ostream& out,Matrix<ElemType,r,c>& A)
{
 for(int i=0;i<A.row;i++)
 {
  for(int j=0;j<A.col;j++)
  {
   out<<A.pElemType[i*A.row+j];
   out<<"   ";
  }
  out<<endl;
 }
 return out;
}

 

 

//***************************************************************************/
//测试函数                                                                  */
//***************************************************************************/

void main()
{
 Vector3f v;
 v(1,1)=4;
 v(1,2)=2;
 Vector3f r;
 r(1,1)=1;
 cout<<RotateOX(v,3);
 cout<<RotateOY(v,3);

 cout<<v;
 Matrix<int,4,4> matrix;
 matrix(4,4)=235;
 matrix(2,1)=3;
 v=v+v;
 cout<<v;
 cout<<matrix+matrix;
 cout<<matrix*matrix;

 cout<<v;
 cout<<VectorMo(v);
 cout<<VectorUnit(v);
 

}