二维计算几何基础

⑴二维计算几何中常用的定义

struct Point
{
    double x,y;
    Point(double x=0,double y=0):x(x),y(y) {}//构造函数，方便函数编写
};

typedef Point Vector;//Vector只是Point的别名

Vector operator + (Vector A,Vector B)
{
    return Vector (A.x+B.x,A.y+B.y);
}

Vector operator - (Vector A,Vector B)
{
    return Vector (A.x-B.x,A.y-B.y);
}

Vector operator * (Vector A,double p)
{
    return Vector (A.x*p,A.y*p);
}

Vector operator / (Vector A,double p)
{
    return Vector (A.x/p,A.y/p);
}

Vector operator < (const Point &A,const Point &B)
{
    return A.x<B.x||(A.x==B.x&&A.y<B.y);
}

Vector operator > (const Point &A,const Point &B)
{
    return A.x>B.x||(A.x==B.x&&A.y>B.y);
}

const double eps=1e-10;
int dcmp(double x)//三态函数，减少了精度问题
{
    if (fabs(x)<eps)
        return 0;
    else
        return x<0?-1:1;
}
bool operator == (const Point &A,const Point &B)
{
    return dcmp(A.x-B.x)==0&&dcmp(A.y-B.y)==0;
}

⑵向量有一个所谓的“极角”，即从x轴正半轴旋转到该向量方向所需要的角度。C标准库里的atan2函数就是用来求极角的，如向量(x,y)的极角就是atan2(y,x)(单位：弧度)

⑶点积。两个向量v和w的点积等于二者长度的乘积再乘上它们的夹角的余弦。向量v和w的夹角是指从向量v到向量w逆时针旋转的角。下面是点积计算的方法，以及利用点积计算向量的长度和夹角的函数。

double Dot(Vector A,Vector B)
{
    return A.x*B.x+A.y*B.y;
}

double Length(Vector A)
{
    return sqrt(Dot(A,A));
}

double Angle(Vector A,Vector B)
{
    return acos(Dot(A,B)/Length(A)/Length(B));
}

⑷叉积。两个向量v和w的叉积等于v和w组成的三角形的有向面积的两倍。

double Cross(Vector A,Vector B)
{
    return A.x*B.y-A.y*B.x;
}

double Area2(Point A,Point B,Point C)
{
    return Cross(B-A,C-A);
}

⑸两个向量的位置关系。把叉积和点积组合到一起，就可以更细致地判断两个向量的位置关系。如图，括号里的第一个数是点积符号，第二个数是叉积的符号，第一个向量v总是水平向右，另一个向量w的各种情况都包含在了图中。

                                             

⑹向量旋转。向量旋转的公式为：

x'=x*cosa - y*sina

y'=x*sina + y*cosa

其中a为逆时针旋转的角度。推到方法，单位圆。

Vector Rotate(Vector A,double rad)
{
    return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));
}

Vector Normal (Vector A)//计算单位法向量
{
    double L=Length(A);
    return Vector(-A.y/L,A.x/L);
}

⑺基于复数的几何计算。

#include <complex>
using namespace std;
typedef complex<double> Point;
typedef Point Vector;
//这样定义之后，我们就自动拥有了构造函数、加减法和数量积。
//用real(p)和imag(p)访问实部和虚部，conj(p)返回共轭复数，即conj(a+bi)=a-bi。
double Dot(Vector A,Vector B)
{
    return real(conj(A)*B);
}

double Cross(Vector A,Vector B)
{
    return imag(conj(A)*B);
}
Vector Rotate(Vector A,double rad)
{
    return A*exp(Point(0,rad));
}

⑻直线的参数表示。直线可以用直线上一点P和方向向量v表示(向量的大小没有什么用处)。直线上所有点Q满足Q=P+tv，其中t称为参数。如果已知直线上的两个不同点A和B,则方向向量为B-A，所以参数方程为A+(B-A)t。射线的t>0，线段的t在0~1之间。

⑼直线交点。参数方程下的直线交点公式，设直线分别为P+tv和Q+tw，设向量u=P-Q，交点在第一条直线的参数为t1，第二条直线上的参数为t2，则x和y坐标可以列出一个方程，解得:

t1=cross(w,u)/cross(v,w);

t2=cross(v,u)/cross(v,w);

Point GetLineIntersection(Point P,Vector V,Point Q,Vector w)
{
    Vector u=P-Q;
    double t=Cross(w,u)/Cross(V,w);
    return P+V*t;
}

⑽点到直线的距离。点到直线的距离是一个常用的函数，可以用叉积计算，即用平行四边形的面积除以底。

double DistanceToLine(Pint P,Point A,Point B)
{
    Vector v1=B-A,v2=P-A;
    return fabs(Cross(v1,v2)/Length(v1));
}

⑾点到线段的距离。点到线段的距离有两种可能。如求点P到线段AB的距离，设投影点为Q，如果Q在线段AB上，则所求距离就是PQ；如果Q在射线BA上，则所求的距离为PA；如果Q在射线AB上，则所求的距离为PB。判断Q的位置可以用点积进行判断（第(5)点的图）。

double DistanceToSegment(Point P,Point A,Point B)
{
    if (A==B)
        return Length(P-A);
    Vector v1=B-A,v2=P-A,v3=P-B;
    if (dcmp(Dot(v1,v2))<0)
    {
        return Length(v2);
    }
    else if (dcmp(Dot(v1,v2))>0)
    {
        return Length(v3);
    }
    else
    {
        return fabs(Cross(v1,v2))/Length(v1);
    }
}

⑿点在直线上的投影。将直线AB写成参数式A+Tv(v为向量AB)，并且设Q的参数为t，那么Q=A+tv。根据PQ垂直于AB，DOt(v,P-(A+tv))，根据分配率，有Dot(v,P-A)-t*Dot(v,v)=0。

Point GetLineProjection(Point P,Point A,Point B)
{
    Vector v=B-A;
    return A+v*(Dot(v,P-A)/Dot(v,v));
}

⒀线段判交。首先定义“规范相交”为两线段恰好有一个公共点，且不再任何一条线段的端点。线段规范相交的充分必要条件是：每条线段的两个端点都在另一条线段的两侧(这里的两侧是指叉积的符号不同)

bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2)
{
    double c1=Cross(a2-a1,b1-a1),c2=Cross(a2-a1,b2-a1);
    double c3=Cross(b2-b1,a1-b1),c4=Cross(b2-b1,a2-b1);
    return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;
}

       如果允许在端点处相交，如果c1和c2都是0，表示两线段共线，这时可能会有部分重叠的情况；如果c1和c2不都是0，则只有一个相交的方法，即某个端点在另一条线段上。为了判断这种情况，还需要一个判断一个点是否在一条线段上(不包含端点)的代码：

bool OnSemgment(Point p,Point a1,Point a2)
{
    return dcmp(Cross(a1-p,a2-p))==0&&dcmp(Dot(a1-p,a2-p))<0;
}