内容参考书籍——《算法竞赛入门经典训练指南》
在平面坐标系下,向量和点一样用两个数,x和y表示。判断相等使用“三态函数”dcmp,减少精度问题。另外,向量有一个所谓的“极角”,即从x轴正半轴旋转到该方向所需的角度。c标准库里的atan2函数就是用来求极角的,例如:向量(x,y)的极角就是atan2(y,x)(单位:弧度)。
点积,高中知识。
叉积。简单地说,两个向量的叉积等于两个向量组成的三角形的有向面积的两倍。(图中带箭头的向量为v,和v同一个起点的另外一条向量为w,忘记标v和w了尴尬)
cross(v,w)>0 cross(v,w)<0
顺着第一个向量v看,,如果w在左边,那么v和w的叉积大于0,否则小于0.如果两个向量共线(方向相同),那么叉积等于0(三角形退化成为线段)。不难发现,叉积不满足交换律。事实上,cross(v,w)=-cross(w,v)。在坐标系下,两个向量的叉积等于xAyB-xByA。
位置。当我们把叉积和点积组合在一起时,我们便可以细致地判断两个向量的位置关系。第一个向量v总是水平向右的,那么另外一个点就可以用这样的方式确定位置关系(点积符号,叉积符号)。例如,当w在第二象限时点积为负但叉积均为正,用(-,+)表示。
旋转。向量可以绕起点旋转,公式为x,= xcosa-ysina , y,= xsina+ycosa。其中a为逆时针旋转的角。
单位法线。旋转90°后把长度归一化即可。
以上内容采用如下代码实现:
1 struct Point 2 { 3 double x,y; 4 Point(double x=0, double y=0):x(x),y(y) {}//构造函数,方便编写代码 5 }; 6 typedef Point Vector;//别名 7 //向量+向量=向量,点+向量=点 8 Vector operator + (Vector A, Vector B){return Vector(A.x+B.x,A.y+B.y);} 9 //点-点=向量 10 Vector operator - (Vector A, Vector B){return Vector(A.x-B.x,A.y-B.y);} 11 //向量*数=向量 12 Vector operator * (Vector A, double p){return Vector(A.x*p,A.y*p);} 13 //向量/数=向量 14 Vector operator / (Vector A, double p){return Vector(A.x/p,A.y/p);} 15 16 bool operator < (const Point& a, const Point& b){return a.x<b.x||(a.x==b.x&&a.y<b.y);} 17 18 const double eps = 1e-10; 19 int dcmp(double x) 20 { 21 if (fabs(x) < eps) return 0; 22 else return x < 0 ? -1 : 1; 23 } 24 25 bool operator == (const Point& a, const Point &b){return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0;} 26 27 //计算点积、向量长度、夹角函数 28 double Dot(Vector A, Vector B){return A.x*B.x+A.y*B.y;} 29 double Length(Vector A){return sqrt(Dot(A,A));} 30 double Angle(Vector A, Vector B){return acos(Dot(A,B)/Length(A)/Length(B));} 31 32 //计算叉积 33 double Cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;} 34 double Area2(Point A, Point B, Point C){return Cross(B-A,C-A);} 35 36 //计算旋转和单位法向量 37 //rad为弧度 38 Vector Rotate(Vector A, double rad){return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));} 39 Vector Normal(Vector A) 40 { 41 double L =Length(A); 42 return Vector(-A.y/L,A.x/L); 43 }
下面是2019.9.24更新
点和直线。采用参数方程表示直线(参数方程在表示直线,线段,射线时其方程的区别仅仅在于参数,射线t>0,线段在0-1之间,直线无限制) 。直线用直线上一点P0和方向向量v表示,则直线上所有点P满足P=P0+tv,其中t为参数。如果已知直线上两不同点A和B,则方向向量为B-A,所以参数方程为A+(B-A)t。
直线交点。设直线为P+tv和Q+tw,设向量u=P-Q,交点在第一条直线的参数为t1,在第二条直线上的参数为t2,则交点的x和y可以列出一个方程,解得:
\[{t_1} = \frac{{cross\left( {w,u} \right)}}{{cross\left( {v,w} \right)}},\quad {t_2} = \frac{{cross\left( {v,u} \right)}}{{cross\left( {v,w} \right)}}\]
此处需要注意的是,从上述公式可以看到在精度要求极高的情况下需要自定义分数类。
点到直线的距离。利用叉积计算,即平行四边形的面积除以底。
点到线段的距离。分两种可能:该点的投影在线段上和不在线段上。若P点的投影在线段(AB)上(Q)则距离为PQ;若不在线段上,则所求距离为PA否则为PB。判断Q的位置可以用前文提到的(点积符号,叉积符号)进行判断。
点的投影。尽管上面的方法避免了求投影点Q,但是有些时候我们不得不求出投影点Q。为此,我们把直线AB写成参数式A+tv(v为向量AB),并且设Q的参数为t0,那么Q=A+t0v,接下来解出t0就能得到Q点,根据PQ垂直于AB,两个向量的点积应该为0,因此Dot(v,P-(A+t0v))=0。根据分配律,有:Dot(v,P-A)-t0*Dot(v,v)=0,这样便解出t0,从而得到Q点。
线段相交。判定两线段是否相交,我们定义“相交规范”为两线段恰好有一个公共点,且不在任何一条线段的端点。其充要条件为:每条线段的两个端点都在另一条线段的两侧(这里的“两侧”是指叉积的符号不同)。
如果允许在端点处相交呢?那么,情况就比较复杂了,c1和c2不都是0,表示两线段共线,这时可能还会出现部分重叠的情况哦~如果c1和c2不都是0,则只有一种相交方法,即某个端点在另外一条线段上。为了判断上述情况是否发生,还需要判断一个点是否在一条线段上(不包含端点),见代码。
多边形有向面积。如果多边形是凸的,可以从第一个顶点出发,把凸多边形分成n-2个三角形,然后把面积加起来。事实上这个方法对非凸多边形也适用,因为三角形的面积是有向的,在外面的部分可以正负抵消。
取P[0]点为划分顶点,一方面可以少算两个叉积(0和任意向量的叉积都等于0),另一方面也减少了乘法溢出的可能性,还不用特殊处理(i=n-1的时候,下一个顶点是P[0]而不是P[n],因为P[n]不存在)。
更新内容代码如下:
1 //交点,注意精度问题 2 //调用该函数前请确保两条直线P+tv和Q+tw有唯一交点。当且仅当Cross(v,w)非0 3 Point GetLineIntersection(Point P, Vector v, Point Q, Vector w) 4 { 5 Vector u = P-Q; 6 double t = Cross(w,u) / Cross(v,w); 7 return P+v*t; 8 } 9 //点到直线的距离。用叉积计算。平行四边形的面积除以底 10 double DistanceToLine(Point P, Point A, Point B) 11 { 12 Vector v1 = B-A,v2= P-A; 13 return fabs(Cross(v1,v2)) / Length(v1); 14 } 15 //点到线段的距离 16 double DistanceToSegment(Point P, Point A, Point B) 17 { 18 if (A==B) return Length(P-A); 19 Vector v1 = B-A,v2 = P-A,v3=P-B; 20 if (dcmp(Dot(v1,v2))<0) return Length(v2); 21 else if (dcmp(Dot(v1,v3))>0) return Length(v3); 22 else return fabs(Cross(v1,v2)) / Length(v1); 23 } 24 //计算投影点 25 Point GetLineProjection(Point P, Point A, Point B) 26 { 27 Vector v = B-A; 28 return A+v*(Dot(v,P-A) / Dot(v,v)); 29 } 30 //按“相交规范”判断相交 31 bool SegmentProperIntersection(Point a1, Point a2, Point b1, Point b2) 32 { 33 double c1 = Cross(a2-a1,b1-a1),c2 = Cross(a2-a1,b2-a1), 34 c3 = Cross(b2-b1,a1-b1),c4 = Cross(b2-b1,a2-b1); 35 return dcmp(c1)*dcmp(c2)<0 && dcmp(c3)*dcmp(c4)<0; 36 } 37 //允许在端点处相交,判断是否发生加入下面一段判断一个点是否在一条线段上的代码 38 bool OnSegment(Point P,Point a1, Point a2) 39 { 40 return dcmp(Cross(a1-P,a2-P)) == 0 && dcmp(Dot(a1-P,a2-P))<0; 41 } 42 //多边形的有向面积 43 double PolygonArea(Point* p, int n) 44 { 45 double area = 0; 46 for (int i = 1; i < n-1; ++i) 47 { 48 area +=Cross(p[i]-p[0],p[i+1]-p[0]); 49 } 50 return area/2; 51 }