计算几何入门

向量(也称为欧几里得向量、几何向量、矢量),指具有大小(magnitude)和方向的量,在程序中一般用两个浮点数来表示x方向和y方向的分量,本文默认的向量起点是坐标原点。

struct VT{
	double x,y;
	VT(){}
	VT(double a,double b){x=a;y=b;}
	VT(const VT &t){x=t.x;y=t.y;}
	VT operator-(){return VT(-x,-y);}
};

向量的长度称为模,为\(\sqrt{x^2+y^2}\)

double len()const{return sqrt(x*x+y*y);}
double operator--(int)const{return sqrt(x*x+y*y);}

单位向量即长为1的向量,一个向量的单位向量方向与该向量相同

VT unit()const{return *this/this->len();}

向量的加减即是让两分量分别加减即可

VT operator+(const VT &t)const{return VT(x+t.x,y+t.y);}
VT operator-(const VT &t)const{return VT(x-t.x,y-t.y);}

向量与数的乘除在几何上就是向量的缩放

VT operator*(double t)const{return VT(x*t,y*t);}
VT operator/(double t)const{return VT(x/t,y/t);}

点积,对于两个平面向量\(a,b\),其点积\(a\cdot b=|a|\cdot|b|\cdot\cos\theta=x_a x_b + y_ay_b\),其中\(\theta\)为两向量间的夹角。
这个公式可以用三角函数的角度加减来证明,其意义在于有时需要求两个向量间的夹角,然而用三角函数的运算的常数是非常大的(远比浮点数的加减乘除要大)。
点积还有一个用处是计算投影,如果\(b\)是单位向量,那么\(a\cdot b=a在b方向上的投影长度\)

double operator*(const VT &t)const{return x*t.x+y*t.y;}
inline double cos2(const VT &a,const VT &b){return (a*b)/a.len()/b.len();}

叉积,对于两个平面向量\(a,b\),其叉积\(a\times b=|a|\cdot|b|\cdot sin\theta =x_ay_b-x_by_a\)。这个公式也可以用三角函数角度加减来证明。
在二维空间中,两向量的叉积等于所构成的平行四边形的面积。而且叉积可以用来判断两向量夹角的正负。

double operator^(const VT &t)const{return x*t.y-t.x*y;}
inline double cross(const VT &a,const VT &b){return a.x*b.y-b.x*a.y;}

向量的旋转,将一个向量逆时针旋转\(\theta\),那么\((x,y)\)就会变成\((x\cos\theta-y\sin\theta,x\sin\theta+y\cos\theta)\)

VT operator+(double t)const{return VT(x*cos(t)-y*sin(t), x*sin(t)+y*cos(t));}
VT operator-(double t)const{return VT(x*cos(t)+y*sin(t),-x*sin(t)+y*cos(t));}

接下来是cmath库里的一些内容。
M_PI:double类型的\(\pi\)常量
asin(x):反正弦,值域为\([-\pi/2,\pi/2]\)
acos(x);反余弦,值域为\([0,\pi]\)
atan(x);反正切,值域为\([-\pi/2,\pi/2]\)
atan2(y,x);反正弦2,值域为\([-\pi,\pi]\)
注意事项:

  1. asin()和acos()的输入必须在\([-1,1]\)内!否则返回的是nan(非数),nan加减乘除任何数还是nan
  2. C++浮点有-0.0!而且默认输出不会把负号去掉!所以要注意asin(-0.0),atan(-0.0)等情况出现
  3. 比较浮点数一定要控制精度!
  4. atan2是先y后x!(可以理解为左边除右边的)
const double eps=1e-8;
inline int dcmp(double x){return fabs(x)<eps?0:(x>0?1:-1);}

有了atan2后就可以求两向量夹角了!\(\tan(\theta)=\frac{a\times b}{a\cdot b}\)
注意值域是\([-\pi,\pi]\),且返回的是a转向b的角度(逆时针)

inline double operator%(const VT &a,const VT &b){return atan2(a^b,a*b);}

接下来是有关线的,还是定义个结构体:

struct LN{
	VT a,b;//a,b为直线的两点/线段的两端点
	LN(){}
	LN(double x1,double y1,double x2,double y2){a.x=x1;a.y=y1;b.x=x2;b.y=y2;}
	LN(const VT &x,const VT &y){a=x;b=y;}
};

要求点到直线的面积,用点到一端点向量叉乘上线段的向量,再除线段长度再取绝对值即可。(这是用平行四边形的面积除以底,得到的就是高)

inline double cdis(const VT &p,const LN &l){
	VT v=p-l.a,u=l.b-l.a;
	return fabs(v^u)/u.len();
}

a,b两条线段相交,当且仅当a线段的两个端点在b线段所在直线的两边,且b线段的两个端点在a线段所在直线的两边。这个就可以用叉积来完成判断。

inline bool its(const LN &u,const LN &v){//intersect
	return dcmp(cross(u.b-u.a,v.a-u.a))*dcmp(cross(u.b-u.a,v.b-u.a))<=0&&
		dcmp(cross(v.b-v.a,u.a-v.a))*dcmp(cross(v.b-v.a,u.b-v.a))<=0;
}

计算两个直线的交点,以其中一条为基准,将其长度进行缩放到与另一直线刚好相交,还是用叉积来计算平行四边形面积:

inline VT cl(const LN &l1,const LN &l2){//cross line
	VT u=l1.b-l1.a,v=l2.b-l2.a;
	return l1.a+u*cross(v,l1.a-l2.a)/cross(u,v);
}

完整模板:

const double eps=1e-8,inf=1e20;
struct VT{
	double x,y;
	VT(){x=y=0.0;}
	VT(double a,double b){x=a;y=b;}
	VT(const VT &t){x=t.x;y=t.y;}
	VT operator-(){return VT(-x,-y);}
	double len()const{return sqrt(x*x+y*y);}
	double operator--(int)const{return sqrt(x*x+y*y);}
	VT unit()const{return *this/this->len();}
	VT operator+(const VT &t)const{return VT(x+t.x,y+t.y);}
	VT operator-(const VT &t)const{return VT(x-t.x,y-t.y);}
	VT operator*(double t)const{return VT(x*t,y*t);}
	VT operator/(double t)const{return VT(x/t,y/t);}
	double operator*(const VT &t)const{return x*t.x+y*t.y;}
	double operator^(const VT &t)const{return x*t.y-t.x*y;}
	VT operator+(double t)const{return VT(x*cos(t)-y*sin(t), x*sin(t)+y*cos(t));}
	VT operator-(double t)const{return VT(x*cos(t)+y*sin(t),-x*sin(t)+y*cos(t));}
};
inline double cos2(const VT &a,const VT &b){return (a*b)/a.len()/b.len();}
inline double operator%(const VT &a,const VT &b){return atan2(a^b,a*b);}
inline double cross(const VT &a,const VT &b){return a.x*b.y-b.x*a.y;}
inline int dcmp(double x){return fabs(x)<eps?0:(x>0?1:-1);}

struct LN{
	VT a,b;//a,b为直线的两点/线段的两端点
	LN(){}
	LN(double x1,double y1,double x2,double y2){a.x=x1;a.y=y1;b.x=x2;b.y=y2;}
	LN(const VT &x,const VT &y){a=x;b=y;}
};
inline double cdis(const VT &p,const LN &l){
	VT v=p-l.a,u=l.b-l.a;
	return fabs(v^u)/u.len();
}
inline bool its(const LN &l1,const LN &l2){//intersect
	return dcmp(cross(l1.b-l1.a,l2.a-l1.a))*dcmp(cross(l1.b-l1.a,l2.b-l1.a))<=0&&
		dcmp(cross(l2.b-l2.a,l1.a-l2.a))*dcmp(cross(l2.b-l2.a,l1.b-l2.a))<=0;
}
inline VT cl(const LN &l1,const LN &l2){//cross line
	VT u=l1.b-l1.a,v=l2.b-l2.a;
	return l1.a+u*cross(v,l1.a-l2.a)/cross(u,v);
}
posted @ 2021-01-16 20:47  Future_Tomorrow  阅读(219)  评论(0)    收藏  举报