计算几何入门
向量(也称为欧几里得向量、几何向量、矢量),指具有大小(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]\)
注意事项:
- asin()和acos()的输入必须在\([-1,1]\)内!否则返回的是nan(非数),nan加减乘除任何数还是nan
- C++浮点有-0.0!而且默认输出不会把负号去掉!所以要注意asin(-0.0),atan(-0.0)等情况出现
- 比较浮点数一定要控制精度!
- 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);
}

浙公网安备 33010602011771号