⭐ 欢迎来到虫虫下载站！ | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站
📤 上传资源

📄 xmath.h

📁 一个简单的数学库.但是很实用..大家可以下载过去研究一下.
💻 H
📖 第 1 页 / 共 3 页
字号:
🔍
123 
// 点积
//
template<class T>
T operator * (const CVector<T,2> &v1, const CVector<T,2> &v2)
{
	CVector<T,2> r=v1;
	return r.Dot(v2);
}

template<class T>
T operator * (const CVector<T,3> &v1, const CVector<T,3> &v2)
{
	CVector<T,3> r=v1;
	return r.Dot(v2);
}

template<class T>
T operator * (const CHomoVector4D<T> &v1, const CHomoVector4D<T> &v2)
{
	CHomoVector4D<T> r=v1;
	return r.Dot(v2);
}

//
// 数量积
//
template<class T>
CVector<T,2> operator * (const CVector<T,2> &v1, T val)
{
	CVector<T,2> r=v1;
	return r.Scale(val);
}

template<class T>
CVector<T,3> operator * (const CVector<T,3> &v1, T val)
{
	CVector<T,3> r=v1;
	return r.Scale(val);
}

template<class T>
CHomoVector4D<T> operator * (const CHomoVector4D<T> &v1, T val)
{
	CHomoVector4D<T> r=v1;
	return r.Scale(val);
}

////////////////////////////////////////////////////////////////////////////////////////

//
// 坐标转换
//

// 2D极坐标
template<class T>
struct CPolar
{
	T	r;			// 半径
	T	theta;		// 角度
};

// 3D柱面坐标
template<class T>
struct CCylindrical
{
	T	r;			// 半径
	T	theta;		// 角度
	T	z;			// Z坐标
};

// 3D球坐标
template<class T>
struct CSpherical
{
	T	r;			// 半径
	T	theta;		// 向量在X-Y面上的投影于正X轴间的夹角
	T	phi;		// 向量与正Z轴的夹角
};

//
// 2D笛卡尔坐标与极坐标之间的转换
//
template<class T>
inline void PolarToCart (const CPolar<T> &p, CVector<T,2> &v)
{
	v.x = p.r * cos(p.theta);
	v.y = p.r * sin(p.theta);
}

template<class T>
inline void CartToPolar (const CVector<T,2> &v, CPolar<T> &p)
{
	p.r = v.Length();
	p.theta = atan2(v.y, v.x);
}

//
// 3D（4D齐次）笛卡尔坐标与柱面坐标之间的转换
//
template<class T>
inline void CylindricalToCart (const CCylindrical<T> &c, CVector<T,3> &v)
{
	v.x = c.r * cos(c.theta);
	v.y = c.r * sin(c.theta);
	v.z = c.z;
}

template<class T>
inline void CylindricalToCart (const CCylindrical<T> &c, CHomoVector4D<T> &v)
{
	v.x = c.r * cos(c.theta);
	v.y = c.r * sin(c.theta);
	v.z = c.z;
	v.w = 1;
}

template<class T>
inline void CartToCylindrical (const CVector<T,3> &v, CCylindrical<T> &c)
{
	c.r = sqrt(v.x*v.x + v.y*v.y);
	c.theta = atan2(v.y,v.x);
	c.z = v.z;
}

template<class T>
inline void CartToCylindrical (const CHomoVector4D<T> &v, CCylindrical<T> &c)
{
	c.r = sqrt(v.x*v.x + v.y*v.y);
	c.theta = atan2(v.y,v.x);
	c.z = v.z;
}

//
// 3D（4D齐次）笛卡尔坐标与球坐标之间的转换
//
template<class T>
void SphericalToCart (const CSpherical<T> &s, CVector<T,3> &v)
{
	T tmp = s.r * sin(s.phi);
	v.x = tmp * cos(s.theta);
	v.y = tmp * sin(s.theta);
	v.z = s.r * cos(s.phi);
}

template<class T>
void SphericalToCart (const CSpherical<T> &s, CHomoVector4D<T> &v)
{
	T tmp = s.r * sin(s.phi);
	v.x = tmp * cos(s.theta);
	v.y = tmp * sin(s.theta);
	v.z = s.r * cos(s.phi);
	v.w = 1;
}

template<class T>
void CartToSpherical (const  CVector<T,3> &v, CSpherical<T> &s)
{
	s.r = v.Length();
	s.theta = atan2(v.y, v.x);
	s.phi = asin(v.x/(s.r*cos(s.theta)));
}

template<class T>
void CartToSpherical (const  CHomoVector4D<T> &v, CSpherical<T> &s)
{
	s.r = v.Length();
	s.theta = atan2(v.y, v.x);
	s.phi = asin(v.x/(s.r*cos(s.theta)));
}

///////////////////////////////////////////////////////////////////////////////

//
// 矩阵
//

//
// 通用矩阵模板
//
template<class T, int R, int C>
class CMatrix
{
//
// 数据
//
public:
	T	M[R][C];	// 按行优先的次序存储矩阵元素
//
// 构造/析构函数
//
public:
	CMatrix() {};

	~CMatrix() {};

	CMatrix(const CMatrix &m)
	{
		*this=m;
	}

	CMatrix& operator = (const CMatrix &m)
	{
		for(int i=0; i<R; ++i)
		{
			for(int j=0; j<C; ++j)
				M[i][j]=m.M[i][j];
		}
		return *this;
	}
//
// 外部接口
//
public:
	// 设置为单位矩阵
	void SetIdentity()
	{
		memset(&(M[0][0]),0,sizeof(T)*R*C);
		if(R!=C)
			return;
		for(int i=0; i<R; ++i)
			M[i][i]=1;
	}
	// 安全的存取操作
	inline void SetVal(int r, int c, T val)
	{
		if(r>=0 && r<R && c>=0 && c<C)
			M[r][c]=val;
	}
	inline void GetVal(int r, int c, T &val) const
	{
		if(r>=0 && r<R && c>=0 && c<C)
			val=M[r][c];
	}
	// 以下的运算符重载使得存取元素更方便
	// 由于不作任何检查，小心内存越界
	inline T operator () (int r, int c) const
	{
		return M[r][c];
	}
	inline T& operator () (int r, int c)
	{
		return M[r][c];
	}
	// 返回行数和列数
	inline int Row() const
	{
		return R;
	}
	inline int Col() const
	{
		return C;
	}
	// 检测是否方阵
	inline bool IsSquare() const
	{
		return (R==C);
	}
	//
	// 矩阵的基本运算
	//
	inline CMatrix& operator += (const CMatrix &m)
	{
		return this->Add(m);
	}
	inline CMatrix& operator -= (const CMatrix &m)
	{
		return this->Sub(m);
	}
	// 矩阵加
	inline CMatrix& operator *= (T val)
	{
		return this->Scale(val);
	}
	CMatrix& Add (const CMatrix &m)
	{
		for(int i=0; i<R; ++i)
		{
			for(int j=0; j<C; ++j)
				M[i][j]+=m.M[i][j];
		}
		return *this;
	}
	CMatrix& Add (const CMatrix &m1, const CMatrix &m2)
	{
		for(int i=0; i<R; ++i)
		{
			for(int j=0; j<C; ++j)
				M[i][j]=m1.M[i][j]+m2.M[i][j];
		}
		return *this;
	}
	// 矩阵减
	CMatrix& Sub (const CMatrix &m)
	{
		for(int i=0; i<R; ++i)
		{
			for(int j=0; j<C; ++j)
				M[i][j]-=m.M[i][j];
		}
		return *this;
	}
	CMatrix& Sub (const CMatrix &m1, const CMatrix &m2)
	{
		for(int i=0; i<R; ++i)
		{
			for(int j=0; j<C; ++j)
				M[i][j]=m1.M[i][j]-m2.M[i][j];
		}
		return *this;
	}
	// 标量乘法
	CMatrix& Scale (T val)
	{
		for(int i=0; i<R; ++i)
		{
			for(int j=0; j<C; ++j)
				M[i][j]*=val;
		}
		return *this;
	}
	// 矩阵乘法，this=m1*m2
	template<int I>
	CMatrix& Mul (const CMatrix<T,R,I> &m1, const CMatrix<T,I,C> &m2)
	{
		T sum;
		for(int i=0; i<R; ++i)
		{
			for(int j=0; j<C; ++j)
			{
				sum=0;
				for(int k=0; k<I; ++k)
				{
					sum+=m1(i,k)*m2(k,j);
				}
				M[i][j]=sum;
			}
		}
		return *this;
	}
	// 将this设为m的转置
	CMatrix& Transpose (const CMatrix<T,C,R> &m)
	{
		for(int i=0; i<R; ++i)
		{
			for(int j=0; j<C; ++j)
			{
				M[i][j]=m.M[j][i];
			}
		}
		return *this;
	}
	// 根据拉普拉斯展开式计算行列式
	// 注意，只有方阵才能调用该函数
	T Det () const
	{
		ASSERT(IsSquare());
		// 1/2/3阶方阵直接计算
		if(Row()==1)
		{
			return M[0][0];
		}
		else if(Row()==2)
		{
			return (M[0][0]*M[1][1]-M[0][1]*M[1][0]);
		}
		else if(Row()==3)
		{
			return (
				M[0][0]*M[1][1]*M[2][2]+
				M[0][1]*M[1][2]*M[2][0]+
				M[0][2]*M[1][0]*M[2][1]-
				M[2][0]*M[1][1]*M[0][2]-
				M[2][1]*M[1][2]*M[0][0]-
				M[2][2]*M[1][0]*M[0][1]
				);
		}
		// 4阶以上递归展开
		CMatrix<T,R-1,C-1> sub;
		T v, r=0;
		for(UINT i=0;i<C;++i)
		{
			// 填充子矩阵
			SubMatrix(sub,0,i);
			// 计算该子阵的代数余子式
			v=(i%2) ? (-M[0][i]) : (M[0][i]);
			v*=sub.Det();
			// 将结果累加
			r+=v;
		}
		return r;
	}
	// 基于Gauss-Jordan消元法计算m的逆阵并保存在this中
	// 如果成功返回true，否则返回false
	bool Inverse (const CMatrix &m)
	{
		if(!IsSquare())
			return false;
		*this=m;
		int i, j, k;
		int row[R], col[C];
		for(k=0; k<R; ++k)
		{
			T max=0;
			// 寻找主元
			for(i=k;i<R;++i)
			{
				for(j=k;j<C;++j)
				{
					if(abs(M[i][j])>max)
					{
						max=M[i][j];
						row[k]=i;
						col[k]=j;
					}
				}
			}
			if(max==0)	// 行列式为0，逆阵不存在
				return false;
			if(row[k]!=k)
			{
				// 行交换
				T tmp;
				for(i=row[k],j=0;j<C;++j)
					SWAP(M[i][j], M[k][j], tmp);
			}
			if(col[k]!=k)
			{
				// 列交换
				T tmp;
				for(j=col[k],i=0;i<R;++i)
					SWAP(M[i][j], M[i][k], tmp);
			}
			T mkk=M[k][k];
			// 计算逆阵元素
			mkk=1/mkk;
			M[k][k]=mkk;
			for(j=0;j<C;++j)
			{
				if(j!=k)
					M[k][j]*=mkk;
			}
			for(i=0;i<R;++i)
			{
				if(i==k)
					continue;
				for(j=0;j<C;++j)
				{
					if(j==k)
						continue;
					M[i][j]-=M[k][j]*M[i][k];
				}
			}
			for(i=0;i<R;++i)
			{
				if(i!=k)
					M[i][k]*=-mkk;
			}
		}
		// 如果之前执行了行/列交换，则需要执行逆交换
		// 交换次序相反，且行（列）交换用列（行）交换代替
		for(k=R-1;k>=0;--k)
		{
			T tmp;
			if(col[k]!=k)
			{
				for(i=col[k],j=0;j<C;++j)
					SWAP(M[i][j], M[k][j], tmp);
			}
			if(row[k]!=k)
			{
				for(j=row[k],i=0;i<R;++i)
					SWAP(M[i][j], M[i][k], tmp);
			}
		}
		return true;
	}
//
// 内部接口
//
protected:
	// 获取去除第rr行rc列后的子阵
	void SubMatrix (CMatrix<T,R-1,C-1> &m, int rr, int rc) const
	{
		UINT i,j,k,l;
		for(i=0,k=0;i<R;++i)
		{
			if(i==rr)
				continue;
			for(j=0,l=0;j<C;++j)
			{
				if(j==rc)
					continue;
				m.M[k][l]=M[i][j];
				++l;
			}
			++k;
		}
	}
};

//////////////////////////////////////////////////////////////////////////

//
// 2X2方阵
//
template<class T>
class CMatrix<T,2,2>
{
// 数据
public:
	// 按行优先的次序存储矩阵元素
	union
	{
		T	M[2][2];
		struct
		{
			T 
			M00, M01,
			M10, M11;
		};
	};

// 构造/析构函数
public:
	CMatrix() {};

	~CMatrix() {};

	CMatrix(const CMatrix &m)
	{
		*this=m;
	}

	CMatrix& operator = (const CMatrix &m)
	{
		M00=m.M00; M01=m.M01;
		M10=m.M10; M11=m.M11;
		return *this;
	}

// 外部接口
public:
	// 
	inline void SetIdentity ()
	{
		M[0][0]=M[1][1]=0;
		M[0][1]=M[1][0]=1;
	}
	// 安全的存取操作
	inline void SetVal(int r, int c, T val)
	{
		if(r>=0 && r<2 && c>=0 && c<2)
			M[r][c]=val;
	}
	inline void GetVal(int r, int c, T &val) const
	{
		if(r>=0 && r<2 && c>=0 && c<2)
			val=M[r][c];
	}
	// 以下的运算符重载使得存取元素更方便
	// 由于不作任何检查，小心内存越界
	inline T operator () (int r, int c) const
	{
		return M[r][c];
	}
	inline T& operator () (int r, int c)
	{
		return M[r][c];
	}
	//
	// 矩阵的基本运算
	//
	inline CMatrix& operator += (const CMatrix &m)
	{
		return this->Add(m);
	}
	inline CMatrix& operator -= (const CMatrix &m)
	{
		return this->Sub(m);
	}
	CMatrix& operator *= (T val)
	{
		return this->Scale(val);
	}
	CMatrix& Add (const CMatrix &m)
	{
		M00+=m.M00; M01+=m.M01; 
		M10+=m.M10; M11+=m.M11; 
		return *this;
	}
	CMatrix& Sub (const CMatrix &m)
	{
		M00-=m.M00; M01-=m.M01; 
		M10-=m.M10; M11-=m.M11; 
		return *this;
	}
	CMatrix& Add (const CMatrix &m1, const CMatrix &m2)
	{
		M00=m1.M00+m2.M00; M01=m1.M01+m2.M01; 
		M10=m1.M10+m2.M10; M11=m1.M11+m2.M11; 
		return *this;
	}
	CMatrix& Sub (const CMatrix &m1, const CMatrix &m2)
	{
		M00=m1.M00-m2.M00; M01=m1.M01-m2.M01; 
		M10=m1.M10-m2.M10; M11=m1.M11-m2.M11; 
		return *this;
	}
	CMatrix& Scale (T val)
	{
		M00*=val; M01*=val;
		M10*=val; M11*=val;
		return *this;
	}
	CMatrix& Mul (const CMatrix &m1, const CMatrix &m2)
	{
		M00=m1.M00*m2.M00+m1.M01*m2.M10;
		M01=m1.M00*m2.M01+m1.M01*m2.M11;
		M10=m1.M10*m2.M00+m1.M11*m2.M10;
		M11=m1.M10*m2.M01+m1.M11*m2.M11;
		return *this;
	}
	CMatrix& Transpose (const CMatrix &m)
	{
		M00=m.M00; M01=m.M10;
		M10=m.M01; M11=m.M11;
		return *this;
	}
	inline T Det () const
	{
		return (M00*M11-M01*M10); 
	}
	bool Inverse (const CMatrix &m)
	{
		T det=m.Det();
		if(det == 0)
			return false;		// 行列式为0，无逆阵
		M00=m.M11/det; M01=-m.M10/det;
		M10=-m.M01/det; M11=m.M00/det;
		return true;
	}
};

///////////////////////////////////////////////////////////////////////////////////

//
// 常用矩阵的类型定义
//
typedef CMatrix<int,2,2>	MATRIX2Di;
typedef CMatrix<int,3,3>	MATRIX3Di;
typedef CMatrix<int,4,4>	MATRIX4Di;

typedef CMatrix<float,2,2>	MATRIX2Df;
typedef CMatrix<float,3,3>	MATRIX3Df;
typedef CMatrix<float,4,4>	MATRIX4Df;

typedef CMatrix<double,2,2>	MATRIX2Dd;
typedef CMatrix<double,3,3>	MATRIX3Dd;
typedef CMatrix<double,4,4>	MATRIX4Dd;

/////////////////////////////////////////////////////////////////////////////////////

//
// 向量与矩阵之间的运算
//

// 向量与矩阵的乘法
template<class T, int R, int C>
void vec_mul_mat(const CVector<T,R> &v, const CMatrix<T,R,C> &m, CVector<T,C> &r)
{
	T sum;
	for(int j=0; j<C; ++j)
	{
		sum=0;
		for(int i=0; i<R; ++i)
			sum+=v(i)*m(i,j);
		r(j)=sum;
	}
}

// 2D向量和3D齐次矩阵相乘
template<class T>
void vec_mul_mat(const CVector<T,2> &v, const CMatrix<T,3,2> &m, CVector<T,2> &r)
{
	// 假设m的最后一列为[0,0,1]T，并使用3X2矩阵表示
	r.x=v.x*m(0,0)+v.y*m(1,0)+m(2,0);
	r.y=v.x*m(0,1)+v.y*m(1,1)+m(2,1);
}

template<class T>
void vec_mul_mat(const CVector<T,2> &v, const CMatrix<T,3,3> &m, CVector<T,2> &r)
123

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -