• 计算几何二维几何前置基础知识


    内容参考书籍——《算法竞赛入门经典训练指南》

      在平面坐标系下,向量和点一样用两个数,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 }
  • 相关阅读:
    关于 log4j.additivity
    JDK8新特性:使用Optional:解决NPE问题的更干净的写法
    异常处理和日志输出使用小结
    搭建DNS服务器
    git 使用技巧
    mysql
    linux学习记录
    nginx解析
    node npm pm2命令简析
    jenkins使用简析
  • 原文地址:https://www.cnblogs.com/125418a/p/11575865.html
Copyright © 2020-2023  润新知