计算几何浅谈

计算几何浅谈

注：此浅谈中运用到部分参考资料以及博客中的术语。

一、前置知识

　　计算几何的学习需要用到高中数学向量的知识。在高中数学中已经涉及到向量的点积，在这里就介绍一下向量的叉积。

　　我们定义两个平面向量：$\vec{a}=（x1,y1）、\vec{b}=（x2,y2）$。对于他们叉积运算可以运用以下公式来求解：$\vec{a}\times\vec{b}=x1\times y2-x2\times y1$。运用上面的式子，可以求出这两个向量的叉积，通过上面的式子，可得出来这个叉积是可正可负的。因为叉积有正有负，所以就可以通过叉积的正负来判断这两个向量的位置关系。

　　在这里我们引入一个叫右手螺旋的判定方法。对于这个判定方法，需要建立一个立体的直角坐标系（如下图），我们将要计算叉积的两个向量的起点平移到原点，并保证这两个向量是处于平面XOY上，伸出右手，四指方向为叉积运算顺序方向，拇指方向为叉积的正负，即拇指指上，叉积为正，拇指指下，叉积为负。对于我们得出的值可以成为有向面积。举两个例子：对于下图中的三个向量——$\vec{AB}、\vec{AC}、\vec{AD}$，我们进行叉积操作。对于$\vec{AB}\times\vec{AD}$，四指方向为逆时针方向，拇指向上，即叉积为正，对于$\vec{AB}\times\vec{AB}$，四指方向为顺时针方向，拇指向下，即叉积为负。

　　了解了右手螺旋，不难发现，当两个向量进行叉积运算的时候，运算顺序为顺时针时，叉积为负，运算顺序为逆时针时，叉积为正。

二、线段与线段的位置关系

　　既然只想知道这两个线的位置关系，即是不是相交，那么能否绕过求交点的过程，直接判断这两个线是否相交呢？我们先观察三张图片。

　　观察上图，不难发现不论两个线段处于怎样的位置，只要满足线段$l_1$的两个端点分列线段$l_2 $两侧，且线段$l_2$的两个端点分列线段$l_1$两侧，这两个线段一定相交（不规范相交，即这两条线段的交点在这两条线段的端点上不算相交，且这两条线段有重叠不算相交）。通过上述的判断方法，似乎可以运用向量的知识，对于线段$AB、CD$，若$\vec{AB}\times\vec{AC}$与$\vec{AB}\times\vec{AD}$的符号相反，说明$C、D$两点一定在线段$AB$的两侧，同理若$\vec{CD}\times\vec{CA}$与$\vec{CD}\times\vec{CB}$的符号相反，说明$A、B$两点一定在线段$CD$的两侧。根据这个性质，我们可以写出以下的判断方法：

int dblcmp(double x)
	{if(fabs(x)<=eps) return 0;return x>0?1:-1;} 
double det(double x1,double y1,double x2,double y2)
	{return x1*y2-x1*y1;}
double cross(Point a,Point b,Point c)
	{return det(b.x-a.x,b.y-a.y,c.x-a.x,c.y-a.y);}
bool check(Point A,Point B,Point C,Point D)
{
    return (dblcmp(cross(c,a,d))^dblcmp(cross(c,b,d)))==-2&&
        (dblcmp(cross(a,c,b))^dblcmp(cross(a,d,b)))==-2;
}

　　解析一下上述代码：函数$dblcmp$是我们自己定义的一种函数，用来判断数字的正负。函数$det$和$cross$就是求出两个向量的叉积。对于$check$似乎很好理解，就是两两求出叉积，并进行正负判断（显然一正一负就是$-1\oplus1=-2​$）。

　　考虑完不规范相交，现在思考规范相交。我们规范相交中的两种情况，即这两条线段的交点在这两条线段的端点上，且这两条线段有重叠。对于这两种情况来说，上述的四个叉积中会有等于$0$的情况（当且仅当两个向量方向相同时，叉积为零）。对于有叉积等于零时，只需要判断是否有交点在两条线段内就可以了。根据这个思路，我们进行编程。

int dblcmp(double x)
	{if(fabs(x)<eps) return 0; return x>0?1:-1;}
double det(double x1,double y1,double x2,double y2)
	{return x1*y2-x2*y1;}
double cross(Point a,Point b,Point c)
	{return det(b.x-a.x,b.y-a.y,c.x-a.x,c.y-a.y);}
int xycmp(double p,double mn,double mx)
	{return dblcmp(p-mn)*dblcmp(p-mx);}
int between(Point a,Point b,Point c)
{
    if(fabs(b.x-c.x)>fabs(b.y-c.y))
        return xycmp(a.x,min(b.x,c.x),max(b.x,c.x));
    return xycmp(a.y,min(b.y,c.y),max(b.y,c.y));
}
bool seg_cross(Point a,Point b,Point c,Point d)
{
    double s1,s2,s3,s4; int d1,d2,d3,d4;
    d1=dblcmp(s1=cross(a,b,c));
    d2=dblcmp(s2=cross(a,b,d));
    d3=dblcmp(s3=cross(c,d,a));
    d4=dblcmp(s4=cross(c,d,b));
    if((d1^d2)==-2&&(d3^d4)==-2) return true;
    if(d1==0&&between(c,a,b)<=0) return true;
    if(d2==0&&between(d,a,b)<=0) return true;
    if(d3==0&&between(a,c,d)<=0) return true;
    if(d4==0&&between(b,c,d)<=0) return true;
    return false;
}

　　对于上方的代码，进行一下解释说明：在函数$deg\_cross$中，第一个$if$判断是判断不规范相交的，后四个判断是判断规范相交中出去不规范相交的其余四种情况。上方的$between$函数，是判断点$a$是否在点$b、c$之间的。

　　以上就是判断两个线段相交的讲解。

三、制作凸包

　　什么是凸包？现对凸包有一个感性的认识：对于条件“在平面上给出$n$个点”想象成为“在一个木板上钉上$n$颗钉子”，对于这$n$个点的凸包想象成为，用橡皮筋加将这$n$颗钉子围起来的图形。如下图，六边形$IJDCLK$为这十二个点的凸包。

　　对于上面的图形，我们可以用一种叫卷包袱的做法，我们发现对于每一个点，与它连接的点一定是最外围的点。这样我们以每一个点为原点，进行极角排序，选出最外侧的点，进行连接，由于我们可以从一个点开始，逆时针转动一圈得到整个凸包，所以我们没有必要每一次连接两侧最外侧的点，只需要链接一侧。（如下图）

　　分析一下上文的时间复杂度：每一次排序为$O(nlog_2^n)$，一共$n$次，总时间复杂度$O(n^2log_2^n)$。但是每一次排序我们只要最外侧的点，所以我们可以$O(n)$扫一遍所有的点，得到此点，这样就可以将时间复杂度优化为$O(n^2) $。

　　说实话$O(n^2)$的做法不怎么样，我们思考有没有跟好的做法。

　　我们找到所有$n$个点中中坐标最低的点，不难想到，这个点一定是凸包上的点，我们以这个点为原点，将剩余的$n-1$个点进行极角排序。观察上图，我们发现，在排序后，排名靠后的点在凸包的逆时针序列（以最低点为开始，逆时针的点号序列）上一定排在后面。

　　运用上一段叙述的性质，我们就能运用单调栈来求出凸包，具体过程如下：

　　1.找到y坐标最小的点(如果有几个y坐标最小，就选其中x坐标最小的)，作为基点（H）。

　　2.把其他每个点到基点的线段与x轴的夹角，然后按这个夹角从小到大排序，然后逆时针扫一遍。

　　3.排序后依次是K,C,D,L,F,G,E,I,B,A,J。

　　4.首先是线段HK，假设他就在凸包上。然后我们扫到KC，因为HK到KC的旋转方向和我们开始扫过来的逆时针方向是相同的，所以我们假设KC也在凸包上。然后到了CD。现在问题来了，KC到CD的旋转方向不一致，所以我们确定C点不在凸包上。以此类推，等到点集中所有点遍历完成，我们就得到HKDGBJ几个凸包点。

　　5.总之确定某一个点是否在凸包上就依靠一个条件：上一个点连到该点的线段和该点连到下一个点的线段的旋转方向是否和开始扫描的方向相同（赤裸裸的凸包性质）

　　根据上述过程我们能写出下方的代码

double dis(Point a) {return sqrt(a.x*a.x+a.y*a.y);}
bool cmp(const Point &a,const Point &b)
	{return (atan2(a.y,a.x)<atan2(b.y,b.x))||
	(atan2(a.y,a.x)<atan2(b.y,b.x)&&dis(a)<dis(b));}
double det(double x1,double y1,double x2,double y2)
	{return x1*y2-x2*y1;}
double cross(Point a,Point b,Point c)
	{return det(b.x-a.x,b.y-a.y,c.x-a.x,c.y-a.y);}
void make()
{
	for(int i=1;i<=n;i++) if((point[i].y<point[1].y)||
		(point[i].y==point[1].y&&point[i].x<point[1].x))
			swap(point[i],point[1]);
	for(int i=2;i<=n;i++)
        point[i].x-=point[1].x,point[i].y-=point[1].y;
	sort(point+2,point+n+1,cmp);
	for(int i=2;i<=n;i++)
        point[i].x+=point[1].x,point[i].y+=point[1].y;
	sta[top=1]=point[1];
	for(int i=2;i<=n;i++)
	{
		while(top>1&&cross(point[i],sta[top-1],sta[top])>0) 
            top--;
		sta[++top]=point[i];
	}
}

　　上面就是构建凸包的讲解。

四、判断点在凸包内

方法一：

　　此方法是先取到一个十分诡异的值（下面举例子取值为$\pi$），将我们要验证的点和$（inf，inf\times\pi）$相连，看连接出的线段与凸包有几个交点，若焦点数为偶数个，则要验证的点在凸包外，反之在凸包内。

int dblcmp(double x)
	{if(fabs(x)<eps) return 0; return x>0?1:-1;}
double det(double x1,double y1,double x2,double y2)
	{return x1*y2-x2*y1;}
double cross(Point a,Point b,Point c)
	{return det(b.x-a.x,b.y-a.y,c.x-a.x,c.y-a.y);}
int xycmp(double p,double mn,double mx)
	{return dblcmp(p-mn)*dblcmp(p-mx);}
int between(Point a,Point b,Point c)
{
    if(fabs(b.x-c.x)>fabs(b.y-c.y))
        return xycmp(a.x,min(b.x,c.x),max(b.x,c.x));
    return xycmp(a.y,min(b.y,c.y),max(b.y,c.y));
}
bool seg_cross(Point a,Point b,Point c,Point d)
{
    double s1,s2,s3,s4; int d1,d2,d3,d4;
    d1=dblcmp(s1=cross(a,b,c));
    d2=dblcmp(s2=cross(a,b,d));
    d3=dblcmp(s3=cross(c,d,a));
    d4=dblcmp(s4=cross(c,d,b));
    if((d1^d2)==-2&&(d3^d4)==-2) return true;
    if(d1==0&&between(c,a,b)<=0) return true;
    if(d2==0&&between(d,a,b)<=0) return true;
    if(d3==0&&between(a,c,d)<=0) return true;
    if(d4==0&&between(b,c,d)<=0) return true;
    return false;
}
bool check(Point a)
{
    sta[top+1]=sta[1]; int times=0;
    Point inf_pla=Point(inf,inf*acos(-1));
    for(int i=1;i<=top;i++)
        if(seg_cross(a,inf_pla,sta[i],sta[i+1]))
	    	times++;
    return times%2;
}

　　这个验证方法是$O(n)$的，在一些题目中可以接受，但是当验证的点变得很多的时候就不好用了。

方法二：

　　我们将一个含有$n$个点的凸包划分为$n-2$个三角区域。这样我们就能将问题转化为要查询的点是否在这$n-2$个三角形区域内。

　　我们可以先二分出来，这个点属于哪一个大区间，即只查找两条射线所夹的区域，查找过后我们在判断是否属于三角形内。

　　看上图：假设我们查询绿色的点是否在凸包内，我们首先二分得到了它所在的区间，然后判断它和绿色的向量的关系，蓝色和紫色的点类似，蓝色的点在边界上，紫色的点在边界右边。

int dblcmp(double x)
	{if(fabs(x)<eps) return 0; return x>0?1:-1;}
double det(double x1,double y1,double x2,double y2)
	{return x1*y2-x2*y1;}
double cross(Point a,Point b,Point c)
	{return det(b.x-a.x,b.y-a.y,c.x-a.x,c.y-a.y);}
bool check(double x,double y)
{
    if(dblcmp(cross(point[1],Point(x,y),point[n]))<0)
        return false;
    if(dblcmp(cross(point[1],point[2],Point(x,y)))<0)
        return false;
    int l=2,r=n;
    while(l<r)
    {
        int mid=(l+r)>>1;
        if(dblcmp(cross(point[1],point[mid],Point(x,y)))>=0) 
            l=mid+1;
        else r=mid;
    } l--;
    return dblcmp(cross(point[l],Point(x,y),point[l+1]))<0;
}