也许是史上最不良心的低阶计算几何讲解和习题集？？

时间 2019-11-07

标签也许史上最不良心低阶计算几何讲解习题集繁體版

原文原文链接

-3.在此声明：

笔者极端厌恶计算几何，以致于直到今天以前都只打过几个计算几何的模板~~~~~算法

不过鉴于今年是18年，因此感受SD颇有可能考计算几何（18年是什么理由啊喂）数组

因而滚过来整理计算几何的资料了......ide

-2.可能会用到的资料：

《计算几何——算法与应用》学习

　　Mark·Overmars 著邓俊辉译优化

《计算几何导论》搜索引擎

　　[美] F·P·普霍帕拉塔 M·I·沙莫斯著spa

——好像是一些没什么用的东西呢3d

《算法导论——第33章》code

　　[美]Tomas H.Cormen Charles E.Leiserson Ronald L.Rivest Clifford Stein 著orm

——据说是能够提高文章逼格的东西呢

本文会尽可能（大概吧）把上述文献中对OI有用的部分扒拉出来

《蓝书》

　　刘汝佳著

《紫书》

　　刘汝佳著

某些大佬的blog——%%%%%%%%%%%%%%%%%%%%%%%%%%%

-1.可能会用到的OJ：

BZOJ——尽管他很爆炸

POJ——尽管全是英文

HDU——我好像没怎么用过

luogu——很良心的OJ吧，大概......

LOJ——几乎没用过，应该也用不到吧

UOJ——不用搜索引擎就能翻出须要的题——若是须要的题真的凑巧在这个OJ里的话

0.因而正文开始：

（废话真多）

计算几何最重要的是(拼命调)关注点和向量

1.几何的准备工做（能够直接略过而和后文没什么联系的工做）：

1.1.通常定义和记法以及一些无关紧要且并不全面的说明：

计算几何的研究对象每每是基于欧几里得空间的点集的，

由于，当咱们存在一个坐标系时

一个点能够有一个坐标，因而与原点构成了向量——实际上，在许多大佬的计算几何模板中都将点和向量宏定义为一个东西，如刘汝佳的《蓝书》

两个点能够构成一个线段或者直线，直线能够定义半平面交

多个点能够构成多边形，或者凸包——事实上凸包和半平面交存在着诸多类似性

多边形能够各类剖分

并且，相比于用其余元素描述上述内容而言，坐标系中的点是更加适用于计算机语言的，由于一个点表示成有序数对以后，与任何数对(pair)无异

固然以上胡扯也能够扯到扩展到更高的维度上去

因而约定有以下不经常使用的记法：

$E^d$表示d维欧几里得空间，（在下文中，d通常大于等于2且不超过3,——因为笔者水平而有的限制）

在d维欧几里得空间中，有以下基本对象：

1.1.1.点、向量——

一个d元组$(x_1,x_2,...,x_d)$，他能够表示一个点P，也能够表示一个起点为$E^d$的原点O，终点为点P的一个有d个份量的向量

1.1.2.线、平面——

给定$E^d$中两个不一样的点$P_1$,$P_1$的线性组合

$$aP_1+(1-a)P_2(a∈R)$$

能够表示$E^d$中的一条直线（若是把P看作向量的话，上式是高中数学中常见的向量基本定理推论）

给定$E^d$中三个不一样的点$P_1$,$P_2$,$P_3$的线性组合

$$aP_1+bP_2+(1-a-b)P_3$$

能够表示$E^d$中的一个平面

相似的上面关于线面的定义也能够推广的更高的维度而成为一种叫作线性簇的东西

（没有用）

线段——

给上文中的直线表达式中加一个对a的范围限制，能够获得一条线段

如，加入限制（0≤a≤1）能够获得一个以$P_1$,$P_2$位端点的线段，能够表示为$\overline{P_1P_2}$或$\overline{P_2P_1}$

1.1.3.凸集——

若$E^d$中的点集D中的任何两点$P_1,P_2$，$\overline{P_1P_2}$中的全部点属于D，则称D是$E^d$中的一个凸集

凸集的交是凸集

直线中的任意线段是凸集

平面内的任意凸多边形是凸集

空间中的任意凸多面体是凸集

关于更高维的凸集的相关性质，可移步《线性规划》——反正笔者看了一点就直接弃疗了的说

凸壳与凸包——

$E^d$中点集S的凸包是$E^d$中包含S的最小凸集，

凸壳是凸包的边界

1.1.4. 多边形——

在$E^2$中多边形定义为线段的一个有限集合使得每一个端点恰为两条边所共有，且没有边的子集具备该性质（多边形仅指这一个边界部分）

注意：

是多边形

而

不是

多边形的其余内容没有什么值得强调的

简单多边形——

若不相邻的边无公共点，则多边形是简单的

如上上图中的多边形就不是简单的

简单多边形把平面划成两个不相交的区域

称为为多边形的内部和外部

凸多边形与星形多边形——

若简单多边形P的内部是个凸集，则此简单多边形是凸多边形

若存在点z在简单多边形P的内部或边上，知足对于P的全部点p,线段$\overline{zp}$属于P的内部或属于P边，则说明多边形P是星形的

（所以每一个凸多边形都是星形的）

然而：

如上图就不是星形多边形

在星形多边形P中，全部z构成的集合为P的核（he?hu?）容易证实，核是一个凸集

（能够用凸集的定义来证）

2.二维几何——点与向量与线段的性质：

这一部分将在$E^2$内

以向量的计算为基础

讨论点与线段的关系，线段与线段的关系

将介绍如何定义一个点，如何定义一个向量

如何实现向量的基本计算

如何用点和向量表示直线、线段，进而用点和向量实现对直线线段的操做

这大概是计算几何中最简单的一部分，将会粘贴一些模板

2.1.二维向量的基本运算：

向量的基本运算做为处理计算几何点线的基础出现

2.1.1.因为一些众所周知的缘由，首先定义精度：

const double eps=1e-10;
int dcmp(double x){
    if(fabs(x)<eps)return 0;
    else return x<0?-1:1;
}//精度

eps&dcmp

2.1.2.定义点与向量：

struct Point{
    double x,y;
    Point(double x=0,double y=0):x(x),y(y){ };
};//定义点
typedef Point Vector;//定义向量

Point&Vector

2.1.3.定义向量加：

Vector operator + (Point A,Point B){
    return Vector(A.x+B.x,A.y+B.y);
}//定义向量加符合三角形定则

Vector +

2.1.4.定义向量减：

Vector operator - (Point A,Point B){
    return Vector(A.x-B.x,A.y-B.y);
}//定义向量减符合三角形定则

Vector -

2.1.5.定义向量数乘：

Vector operator * (Vector A,double p){
    return Vector(A.x*p,A.y*p);
}//定义向量数乘
Vector operator * (double p,Vector A){
    return Vector(A.x*p,A.y*p);
}//定义向量数乘

Vector *

（直接分别定义左右乘了，并不会其余更高端的方式）

2.1.6.定义向量数除：

Vector operator / (Vector A,double p){
    return Vector(A.x/p,A.y/p);
}//定义向量数除

Vector /

2.1.7.定义向量相等（用到了精度）：

bool operator == (const Point& a,const Point& b){
    return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;
}//定义向量相等为彻底相同

Vector ==

2.1.8.向量极角：

即从x轴正半轴旋转到该向量方向所需的弧度，（单位：弧度）

double Polar_Angle (Vector A){
    return atan2(A.y,A.x);
}//计算向量极角

Polar Angel

2.1.9.向量旋转：

将向量p(x,y)绕起点逆时针旋转弧度a，获得$p'(xcosa-ysina,xsina+ycosa)$

2.1.10.二维向量的点积：

若是你已经学习了《数学必修》的内容，则你应该对点积十分了解

double Dot(Vector A,Vector B){
    return A.x*B.x+A.y*B.y;
}//计算向量点积

Dot

坐标对应相乘再相加，点积的结果是一个实数，可认为是模长乘夹角cos

2.1.11.二维向量的叉积：

若是你有必定的物理或数学水平，则你应该明白二维向量的叉积结果应该是三维空间中的一个垂直于这个两个向量向量，方向按照右手法则断定

然而由于使用范围的限制，这里对叉积的定义十分浅显：

考虑向量$p_1(x_1,y_1),p_2(x_1,y_1)$的叉积定义为把两个向量起点移动到同一处后获得的三点围成的三角形的有向面积的两倍

或者定义为下述矩阵的行列式的值：

若其为正，则相对于原点来讲$p_1$位于$p_2$的顺时针方向，若其为0，则两向量共线

double Cross(Vector A,Vector B){
    return A.x*B.y-A.y*B.x;
}//计算向量叉积

Cross

交叉相乘再相减，可看作模长乘夹角sin

注意：叉积没有交换律

2.1.12.由点积叉积拓展来：

因而有了一个简单的基于叉积的计算三点围成三角形面积的算法：

double Area(Point A,Point B,Point C){
    return fabs(Cross(B-A,C-A)/2);
}//计算abc三点三角形的数值面积

Area

因而有了一个简单的基于点积的计算向量长的算法：

double Length(Vector A){
    return sqrt(Dot(A,A));
}//计算向量模长

Length

因而有了一个简单的基于点积的计算向量夹角的算法：

double Angle(Vector A,Vector B){
    return acos(Dot(A,B)/Lenth(A)/Lenth(B));
}//计算向量夹角（有序的）

Angle

2.2.经过向量来断定线段之间的关系：

线段取其两端点，定义为无序的一对点

而直线则取其中一点和直线的方向向量，定义为一个点和一个向量的二元组

（也能够看作两个向量，也能够看作两个点）

这一部分的存在应用价值（吗?）

2.2.1.肯定连续线段的偏转方向：

讨论在点$p_1$处，连续线段$\overline{p_0p_1}$,$\overline{p_1p_2}$的旋转方向，即断定给定角$∠p_0p_1p_2$的转向，而咱们能够经过计算$\overrightarrow{p_0p_1}$$\overrightarrow{p_0p_2}$叉积结果，来避免直接对角的计算，若其为正，则说明后一条线段向逆时针偏折；若其为负，则说明后一条线段向顺时针偏折；若其为零，则两条线段共线

int Direction(Vector p0,Vector p1,Vector p2){
    double num=Cross(Vector(p1-p0),Vector(p2-p0));
    return dcmp(num);
}//判断p0->p1,p1->p2的偏转方向，1为逆时针，-1为顺时针，0为共线

Direction

2.2.2.断定直线交点：

在以前，咱们定义的线性组合$a\vec{A}+(1-a)\vec{B}(a∈R)$为直线，

简单地，咱们定义直线的方向向量为$\vec{v}=\vec{A}-\vec{B}$，

化简，因而有了$B+a\vec{v}$

方向向量界定了直线的方向，向量B的终点B界定了直线在笛卡尔坐标系中的位置，使用规范的符号,因而有了参数方程：

$$P+t\vec v$$

其中t的不一样取值表明了直线上的不一样点

因而，设$l_1:P+t\vec v $，$l_2:Q+t\vec w $，联立两参数方程可解得交点在$l_1$的参数$t_1$，在$l_2$的参数$t_2$，若设$\vec u =P-Q$则有关于$t_1,t_2$的公式以下：

$$t_1={cross(\vec w,\vec u)\over cross(\vec v,\vec w)},t_2={cross(\vec v,\vec u)\over cross(\vec v,\vec 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;
}//用参数方程求取直线交点——若是他存在的话

Get Line Intersection

实际上，运用简单的高中解析几何知识就足以解决此问题了吧

2.3.经过向量处理点与直线、点与线段的关系：

这一部分有应用的价值（吗）

2.3.1.点到直线的距离：

直接叉积求解便可

double DistanceToLine(Point P,Point A,Point B){
    return abs(Cross(B-A,P-A)/Length(B-A));
}//P到直线AB的距离

Distance To Line

2.3.2.点到线段的距离：

须要使人十分不愉快的分类讨论：

设P在直线AB上的投影点为Q，

1：Q在线段外，据A较近，须要求解PA的长

2：Q在线段AB上，直接使用点到直线距离

3：Q在线段外，据B较近，须要求解PB的长

断定能够经过点积来实现

double DistanceToSigment(Point P,Point A,Point B){
    if(A==B)return Length(P-A);
    if(dcmp(Dot(P-A,B-A))<0)    return Length(P-A);
    if(dcmp(Dot(P-B,A-B))<0)    return Length(P-B);
    return DistanceToLine(P,A,B);
}//P到线段AB的距离

Distance To Sigment

习题：

POJ P2318

POJ P1269

UVa11796(UVA的小题目挺有意思的，适合心情好的时候)

3.二维几何——简单多边形与圆与二维几何对点集和平面的处理：

因而咱们很愉快地从对点线的讨论正式转入图形了

这一部分将讨论对多边形问题的解决，

进而研究凸包凸壳半平面交等基于点线对象

而后将给出圆的定义

3.1.简单多边形部分：

多边形是在二维几何中常常研究的图形，

多边形取其全部顶点和他们的逆时针链接顺序而定义为一个有序点集

具体的定义已经在本文的第一部分说起，因此不做赘述，

接下来将介绍简单多边形围成有向面积的求法

点与多边形关系的断定

3.1.1.简单多边形的有向面积：

考虑凸多边形的有向面积？

能够从一个点出发把凸多边形划成n-2个三角形，求面积再求和

这个方法对非凸的简单多边形一样适用——由于三角形面积的有向性使得在多边形外部的部分被抵消掉

具体的作法是，把多边形的n个顶点按照链接顺序从0到n-1编号，从$P_0$依次向$P_1$到$P_{n-1}$连向量$\vec{v}$，计算$\sum_{i=1}^{n-2}{cross(\vec{v_i},\vec{v_{i+1}})\over 2}$

（深色部分正负相抵）

double PolygonArea(Point *p,int n){
    double area=0;
    int i;
    for(i=1;i<n-1;i++)
        area+=Cross(p[i]-p[0],p[i+1]-p[0]);
    return area/2;
}//计算多边形的有向面积，也许在结果上加个||更有用些？？

Polygon Area

3.1.2.点在多边形内部的断定：

一下算法适用于全部点按链接顺序逆时针排列的简单多边形

如何断定p点是否在多边形Poly内部呢

粗略地想，从点p向右引一条平行于x轴的射线

因为射线无限长而多边形内部范围有限，

因此射线无限远处的尽头所在的区域必定是多边形外部

因为多边形每一条边两侧必定分属多边形的内外

因此从射线远离p的一侧向p走，

第一次通过一条边将到达多边形内

第二次通过一条边将到达多边形外

以此类推通过——

若该射线奇数次穿过边，则点p在多边形内，反之亦然

一些细节，

若是射线的某段与某条边重合，则这条边不该计入次数

若是射线穿过某端点，应该如何断定呢

比较好的方法是定义边的方向为从编号小的点到编号大的点

这样之后再定义

从下往上穿过射线的边包含起点不包含终点

从上往下穿过射线的边包含终点不包含起点

这样若穿过的端点所在的两边同向则只被计算一次

若穿过的端点所在的边反向则要么一次都不计算，要么直接算两次

int isPointInPolygon(Point p,Point poly[],int n){
    int wn=0,i,k,d1,d2;
    for(i=0;i<n;i++){
        if(!dcmp(DistanceToSigment(p,poly[i],poly[(i+1)%n])))    return -1;//在边界上 
        k=dcmp(Cross(poly[(i+1)%n]-poly[i],p-poly[i]));
        d1=dcmp(poly[i].y-p.y);
        d2=dcmp(poly[(i+1)%n].y-p.y);
        if(k>0&&d1<=0&&d2>0)wn++;//点在左，下上穿
        if(k<0&&d2<=0&&d1>0)wn++;//点在右，上下穿
    }
    return wn&1;//T:内部；F：外部 
}//断定点是否在多边形内部

is Point In Polygon

3.2.凸多边形、对踵点和旋转卡壳算法：

凸多边形的对踵点是凸多边形上存在的不仅一对的点对，

每对对踵点$(p_1,p_1)$知足都是凸多边形的顶点且存在过$p_1$的某条凸多边形的切线与某条过$p_2$的凸多边形切线平行

对踵点个数不超过(3n/2)

旋转卡壳算法是一种经过处理凸多边形的对踵点来解决一系列凸多边形问题的算法

笔者将旋转卡壳算法从多边形中剥离出来（主要是由于他内容比较多）

这里给出用旋转卡壳枚举全部对踵点的方法：

先找到y最小的点$p=p_{ymin}$，和y最大的点$q=p_{ymax}$，他们是一对对踵点，

从这对对踵点分别向逆时针方向引出平行于x轴的射线，

按相同的角速度逆时针转动这对射线，直到其中一条穿过多边形的下一端点$p_{next}$，{p_{next}}将代替他所在的射线的原端点成为射线的新端点，并和另外一射线的端点构成新的对踵点，

继续旋转新的这对射线以得到新的对踵点

当咱们枚举了全部角度的平行线以后，天然也就枚举出了全部对踵点

枚举角度的效率是INF的（十分高效）

然而并非全部角度的平行切线都切与不一样的对踵点，因此在转动的过程当中有许多角度能够直接略过

具体地说，同一角速度转动这一说法太过玄学，具体实现应该是，

把上述过程按照两个射线的顶点是否有其一被替代来划分为许多阶段，

能够经过叉积来断定$\overrightarrow {pp_{next}}$和$\overrightarrow{q_{next}q}$的极角大小，

进而肯定p与q被替代的前后顺序

而后直接在两个阶段之间跳跃，由于这期间没有新的对踵点产生

3.2.1.求凸多边形的直径：

凸多边形的直径是凸多边形边上全部点中距离最远的一对点的距离，

他显然是某对对踵点的距离，枚举全部对踵点便可

效率$O(n)$（对每一对对踵点O(1)）

double RotatingCaliper_diameter(Point poly[],int n){
    int p=0,q=n-1,fl,i;
    double ret=0;
    for(i=0;i<n;i++){
        if(poly[i].y>poly[q].y)    q=i;
        if(poly[i].y<poly[p].y)    p=i;
    }
    Point ymin=poly[p],ymax=poly[q];
    for(i=0;i<n*3;i++){
        if(ret<Length(poly[p%n]-poly[q%n]))    ret=Length(poly[p%n]-poly[q%n]);
        fl=dcmp(Cross(poly[(p+1)%n]-poly[p%n],poly[q%n]-poly[(q+1)%n]));
        if(!fl){
            if(ret<Length(poly[p%n]-poly[(q+1)%n]))    ret=Length(poly[p%n]-poly[(q+1)%n]);
            if(ret<Length(poly[q%n]-poly[(p+1)%n]))    ret=Length(poly[q%n]-poly[(p+1)%n]);
            p++,q++;
        }
        else{
            if(fl>0)    p++;
            else        q++;
        }
    }
    return ret;
}//用旋转卡壳求解凸多边形直径

R.C_diameter

习题：

POJ P2187——须要求凸包，而后求凸包的直径的平方输出

3.2.2.求凸多边形的宽度：

凸多边形的宽度是凸多边形的全部平行切线对的距离的最小值

对任意一对平行切线$l_1^a,l_2^a$，必有一对平行切线$l_1^b,l_2^b$，其中至少一条与某一整条边重合

知足$dis(l_1^b,l_2^b)<=dis(l_1^a,l_2^a)$

因而枚举此类平行切线便可

因为，旋转卡壳的过程正是由出现此类平行切线的时间点来划分的，因此能够枚举全部此类平行切线

效率$O(n)$（对每一对对踵点O(1)）

double RotatingCaliper_width(Point poly[],int n){
    int p=0,q=n-1,fl,i;
    double ret;
    for(i=0;i<n;i++){
        if(poly[i].y>poly[q].y)    q=i;
        if(poly[i].y<poly[p].y)    p=i;
    }
    ret=Length(poly[p]-poly[q]);
    for(i=0;i<n*3;i++){
        fl=dcmp(Cross(poly[(p+1)%n]-poly[p%n],poly[q%n]-poly[(q+1)%n]));
        if(!fl){
            if(ret>DistanceToLine(poly[p%n],poly[q%n],poly[(q+1)%n]))    ret=DistanceToLine(poly[p%n],poly[q%n],poly[(q+1)%n]);
            p++,q++;
        }
        else{
            if(fl>0){
                if(ret>DistanceToLine(poly[q%n],poly[p%n],poly[(p+1)%n]))    ret=DistanceToLine(poly[q%n],poly[p%n],poly[(p+1)%n]);
                p++;
            }
            else{
                if(ret>DistanceToLine(poly[p%n],poly[q%n],poly[(q+1)%n]))    ret=DistanceToLine(poly[p%n],poly[q%n],poly[(q+1)%n]);
                q++;
            }
        }
    }
    return ret;
}//用旋转卡壳求解凸多边形宽度

Rotating Caliper Width

（没找题，正好yhzq大神作的某次比赛有这个题，正好用那个题实验了板子的正确性2333）

3.2.3.求凸多边形的最小面积外接矩形：

这样的矩形必定与凸多边形至少一边发生重合，因而这一条重合的边以及这条边的对边能够经过旋转卡壳来枚举，

在枚举出这么一组有重合平行边以后，如何枚举另一对平行边呢，

在咱们枚举出第一对有重合平行边并找到其对应的另外一对边以后（这个另外一对边方向与有重合平行边垂直，因此这对边其实能够存成多边形上的一对点）

若是咱们旋转这组有重合平行边获得另外一组有重合平行边的话，

新的一组的对应另外一对边能够由旧的一组经过向相同的方向旋转得来，

断定是否完成旋转的方法能够是向量叉积

效率$O(n)$（对每一对对踵点O(1),另外一对平行线的旋转也是单向的因而也是O(n)）

其实求最小周长外接矩形也是同理

void RC(Poi poly[],int n){
    int p=0,q=n-1,l=0,r=n-1;
    int fl,i,j;
    Vec nm,dr;
    for(i=0;i<n;i++){
        if(poly[i].y<poly[p].y)    p=i;
        if(poly[i].y>poly[q].y)    q=i;
        if(poly[i].x<poly[l].x)    l=i;
        if(poly[i].x>poly[r].x)    r=i;
    }
    an[0]=intersect_p(poly[p],Vec(1,0),poly[l],Vec(0,1)),an[1]=intersect_p(poly[p],Vec(1,0),poly[r],Vec(0,1));
    an[2]=intersect_p(poly[r],Vec(0,1),poly[q],Vec(1,0)),an[3]=intersect_p(poly[l],Vec(0,1),poly[q],Vec(1,0));
    ans=fabs(Area(an[0],an[1],an[2]));
    for(i=0;i<n*3;i++){
        fl=dcmp(Cross(poly[(p+1)%n]-poly[p%n],poly[q%n]-poly[(q+1)%n]));
        if(!fl){
            dr=poly[(p+1)%n]-poly[p%n],dr=dr/Len(dr);
            nm=Vec(dr.y,-dr.x);
            while(dcmp(Cross(poly[(l+1)%n]-poly[l%n],nm))>0)    l++;
            nm=Vec(0,0)-nm;
            while(dcmp(Cross(poly[(r+1)%n]-poly[r%n],nm))>0)    r++;
            cha[0]=intersect_p(poly[p%n],dr,poly[l%n],nm),cha[1]=intersect_p(poly[p%n],dr,poly[r%n],nm);
            cha[2]=intersect_p(poly[r%n],nm,poly[q%n],dr),cha[3]=intersect_p(poly[l%n],nm,poly[q%n],dr);
            if(fabs(Area(cha[0],cha[1],cha[2]))<ans){
                for(j=0;j<4;j++)    an[j]=cha[j];
                ans=fabs(Area(an[0],an[1],an[2]));
            }
        }
        else{
            if(fl>0){
                dr=poly[(p+1)%n]-poly[p%n],dr=dr/Len(dr);
                nm=Vec(dr.y,-dr.x);
                while(dcmp(Cross(poly[(l+1)%n]-poly[l%n],nm))>0)    l++;
                nm=Vec(0,0)-nm;
                while(dcmp(Cross(poly[(r+1)%n]-poly[r%n],nm))>0)    r++;
                cha[0]=intersect_p(poly[p%n],dr,poly[l%n],nm),cha[1]=intersect_p(poly[p%n],dr,poly[r%n],nm);
                cha[2]=intersect_p(poly[r%n],nm,poly[q%n],dr),cha[3]=intersect_p(poly[l%n],nm,poly[q%n],dr);
                if(fabs(Area(cha[0],cha[1],cha[2]))<ans){
                    for(j=0;j<4;j++)    an[j]=cha[j];
                    ans=fabs(Area(an[0],an[1],an[2]));
                }
                p++;
            }
            else{
                dr=poly[(q+1)%n]-poly[q%n],dr=dr/Len(dr);
                nm=Vec(-dr.y,dr.x);
                while(dcmp(Cross(poly[(l+1)%n]-poly[l%n],nm))>0)    l++;
                nm=Vec(0,0)-nm;
                while(dcmp(Cross(poly[(r+1)%n]-poly[r%n],nm))>0)    r++;
                cha[0]=intersect_p(poly[p%n],dr,poly[l%n],nm),cha[1]=intersect_p(poly[p%n],dr,poly[r%n],nm);
                cha[2]=intersect_p(poly[r%n],nm,poly[q%n],dr),cha[3]=intersect_p(poly[l%n],nm,poly[q%n],dr);
                if(fabs(Area(cha[0],cha[1],cha[2]))<ans){
                    for(j=0;j<4;j++)    an[j]=cha[j];
                    ans=fabs(Area(an[0],an[1],an[2]));
                }
                q++;
            }
        }
    }
}

Minimum external rectangle area

习题：

bzoj P1185==luogu P3187 [HNOI2007]最小矩形覆盖

(luogu没有SPJ，卡精度；BZOJ输出九个nan不知还能不能过)

3.3.凸包部分：

凸包的定义前面已经给出了，并无什么值得强调的

值得一提的是某个问题——

定义一个物品集有n个物品，每一个物品彻底由一个大小为m的元素集中的一些元素按照某个比例融合而成

设在i物品中j元素占比为$x_{i,j}$

则有$\sum_{j=1}^mx_{i,j}=1(i=1...n)$

每次询问给出一种新的物品中各元素的占比，问他可否被原来物品集中的物品按必定比例融合来获得

这个问题在m为2时能够看作一个点是否在凸包内的断定问题

m为3时能够看作一个点是否在一个三维凸包内的断定问题

。。。

好吧，咱们仍是先看看二维凸包的求法吧

3.3.1.凸包、凸壳的求法——Andrew算法：

求凸包，实际是求在凸包边上的端点，或者说，是求凸壳的过程

其流程是：

先把全部点按x为第一关键字，y为第二关键字排序而后去重

而后将第一第二个点加入凸壳，

从第三个点开始，

当新加入多是凸壳第i个点的$p_i$时，断定$p_{i-1}$是否须要被剔出凸壳

若连续线段$\overline{p_{i-2}p_{i-1}}$和$\overline{p_{i-1}p_i}$是向顺时针偏折的，则把$p_{i-1}$剔出凸包，

而后继续断定$p_{i-2}$是否须要被剔出凸包......直到当咱们断定$p_j$时，发现他不须要被剔出凸包，

而后把$p_i$放入凸壳，发现他如今多是凸壳第j+1个点（固然咯，也可能什么都不是）

而后继续加入凸壳的第j+2个点，重复上述过程

当咱们结束对n的处理后，咱们获得了一个下凸壳——他的每条线都成下凸的趋势

这样再从n开始倒着来一遍就求出了上凸壳，合起来就是完整的凸壳啦

时间复杂度为排序的O(nlogn).

int ConvexHull(Point*inp,int n,Point*outp){
    sort(inp,inp+n);
    int m=-1,k,i,j;
    for(i=0;i<n;i++){
        while(m>0&&Direction(outp[m-1],outp[m],inp[i])<=0)m--;
        outp[++m]=inp[i];
    }
    k=m;
    for(i=n-2;i>=0;i--){
        while(m>k&&Direction(outp[m-1],outp[m],inp[i])<=0)m--;
        outp[++m]=inp[i];
    }
    if(m&&outp[m]==outp[0])m--;
    return m+1;
}//假设输入的inp数组已经去太重了，我还要写个离散吗（笑）

Convex Hull

习题：

HDU1348——读入n个点，和一个R，求n个点凸壳的周长+$2\pi R$做为答案

3.3.2.动态凸包：

通常须要排序的东西，搞成动态的话，好像都是套个堆或者平衡树的吧

用set|平衡树维护凸包能够支持动态加入新点

（若是只有删除操做的话，能够离线倒序变成加入）

不能同时支持删除......

具体的说，在set分别维护上下凸壳中的点按先x再y的顺序排完序后的结果，

当加入一个点后，先判断他可能加入上凸壳仍是下凸壳，

而后若是他须要被加入凸壳的话，找到他从按先x再y的顺序排完序后应该在位置，

判断他在set中的前驱和后继是否须要被删除，相似下图的状况是须要删除的：

删除以后判断新的前驱和后继，接着删除须要删除的，直到前驱后继合法为止

效率:$O(nlogn)$（每一个点只会被插入和删除最多一次，实际上用set，而后删除时用迭代器遍历可使得常数变得更小？？）

这里给出一个十分简单的例题：

bzoj P2300==luogu P2521 [HAOI2011]防线修建

这个题目只要求维护上凸壳以及上凸壳的总长，并且由于题目的限制使得这个题的细节比模板简单

用这个题的代码代替模板的缘由是[数据删除]

......好吧，是笔者set用得太烂了，给出的代码实在很差意思当模板

#include<set>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define db double
using namespace std;
const db eps=1e-10;
int dcmp(db x){
    if(fabs(x)<eps)return 0;
    return x>0?1:-1;
}
struct Poi{
    db x,y;
    Poi (db x=0,db y=0):x(x),y(y){    };
};
typedef Poi Vec;
Vec operator + (Vec a,Vec b){    return Vec(a.x+b.x,a.y+b.y);    }
Vec operator - (Vec a,Vec b){    return Vec(a.x-b.x,a.y-b.y);    }
Vec operator * (Vec a,db t){    return Vec(a.x*t,a.y*t);    }
Vec operator * (db t,Vec a){    return Vec(a.x*t,a.y*t);    }
Vec operator / (Vec a,db t){    return Vec(a.x/t,a.y/t);    }
bool operator < (Vec a,Vec b){    return dcmp(a.x-b.x)<0||(!dcmp(a.x-b.x)&&dcmp(a.y-b.y)<0);    };
bool operator > (Vec a,Vec b){    return dcmp(a.x-b.x)>0||(!dcmp(a.x-b.x)&&dcmp(a.y-b.y)>0);    };
bool operator == (Vec a,Vec b){    return !dcmp(a.x-b.x)&&!dcmp(a.y-b.y);    };
set <Poi > C_Hull_set;
set <int >aaa;
int n,x,y,m,q,top;
db an_stk[200010];
Vec city[100010],c0,c1,c2;
bool lim[100010];
int ask[200010][2];
db Dot(Vec ,Vec );
db Cross(Vec ,Vec );
db Len(Vec );
db Area(Poi ,Poi ,Poi );
int drc(Poi ,Poi ,Poi );
db C_ins(Poi ,db );
int main()
{
    int i,j,k;
    Poi now_l,now_r;
    scanf("%d%d%d",&n,&x,&y);
    c0=Poi(0,0),c1=Poi(x,y),c2=Poi(n,0);
    an_stk[0]=Len(c1-c0)+Len(c2-c1);
    C_Hull_set.insert(c0);
    C_Hull_set.insert(c1);
    C_Hull_set.insert(c2);
    scanf("%d",&m);
    for(i=1;i<=m;i++)
        scanf("%lf%lf",&city[i].x,&city[i].y);
    scanf("%d",&q);
    for(i=1;i<=q;i++){
        scanf("%d",&ask[i][0]);
        if(ask[i][0]==1)    scanf("%d",&ask[i][1]),lim[ask[i][1]]=true;
    }
    for(i=1;i<=m;i++)
        if(!lim[i])
            an_stk[0]=C_ins(city[i],an_stk[0]);
    an_stk[++top]=an_stk[0];
    for(i=q;i>0;i--){
        if(ask[i][0]==2)
            an_stk[++top]=an_stk[top-1];
        else
            an_stk[top]=C_ins(city[ask[i][1]],an_stk[top]);
    }
    for(i=top-1;i>0;i--)
        printf("%.2lf\n",an_stk[i]);
    return 0;
}
db Dot(Vec a,Vec b){
    return a.x*b.x+a.y*b.y;
}
db Cross(Vec a,Vec b){
    return a.x*b.y-b.x*a.y;
}
db Len(Vec v){
    return sqrt(Dot(v,v));
}
db Area(Poi a,Poi b,Poi c){
    return Cross(b-a,c-a);
}
int drc(Poi p0,Poi p1,Poi p2){
    return dcmp(Cross(p1-p0,p2-p0));
}
db C_ins(Poi p,db las_ans){
    set <Poi >:: iterator iter1;
    set <Poi >:: iterator iter2;
    set <Poi >:: iterator iter3;
    iter2=C_Hull_set.upper_bound(p);
    iter1=iter2,iter1--;
    if(Cross(p-*iter1,*iter2-*iter1)>=0)return las_ans;
    las_ans-=Len(*iter2-*iter1);
    iter2=iter1,iter2--;
    for(;!(*iter1==c0)&&drc(*iter2,*iter1,p)>=0;iter1--,iter2--,C_Hull_set.erase(iter3)){
        las_ans-=Len(*iter1-*iter2);
        iter3=iter1;
    }
    iter1=C_Hull_set.upper_bound(p);
    iter2=iter1,iter2++;
    for(;!(*iter1==c2)&&drc(p,*iter1,*iter2)>=0;iter1++,iter2++,C_Hull_set.erase(iter3)){
        las_ans-=Len(*iter1-*iter2);
        iter3=iter1;
    }
    iter2=C_Hull_set.upper_bound(p);
    iter1=iter2,iter1--;
    las_ans+=Len(*iter2-p)+Len(*iter1-p);
    C_Hull_set.insert(p);
    return las_ans;
}

[HAOI2011]set+C_Hull

3.4.有向直线、半平面、与半平面交：

给直线加一个方向，他便成为了有向直线，

有向直线能够区分左右——逆时针转动方向向量，率先扫过的一侧为左，另外一侧为右，

规定有向直线左侧平面为该有向直线肯定的半平面，

半平面的交集为半平面交，

3.4.1.有向直线的定义方式：

因为咱们对直线的定义方式是线上一点与方向向量，

因此有向直线的定义只要沿用直线的定义，并确保更规范地使用方向向量便可

struct Line{
    Point P;
    Vector v;
    Line (Point P,Vector v):P(P),v(v){    };
};//定义有向直线

定义有向直线

bool operator  == (Line l1,Line l2){
    return !dcmp(atan2(l1.v.y,l1.v.x)-atan2(l2.v.y,l2.v.x));
}//定义直线比较为极角比较 
bool operator < (Line l1,Line l2){
    return dcmp(atan2(l1.v.y,l1.v.x)-atan2(l2.v.y,l2.v.x))<0;
}//定义直线比较为极角比较 
bool operator > (Line l1,Line l2){
    return dcmp(atan2(l1.v.y,l1.v.x)-atan2(l2.v.y,l2.v.x))>0;
}//定义直线比较为极角比较

用极角比较直线

3.4.2.半平面交：

半平面是平面中的一个点集知足全部点都在某有向直线的左侧，显然是个凸集

因此半平面交也是个凸集

求半平面交的算法在咱们解决斜率优化时已经有所涉猎，

然而因为斜率优化中自有的限制因此这里给出的规范算法比以前的有更多的细节

这个算法与求凸包的方法相似：

排序，而后逐一扫描全部直线试图将其加入表明半平面交边界的队列，在加入以前先肯定队尾是否有直线会被彻底取代

方法是取出队尾直线$l_i$和队尾倒数第二条直线$l_{i-1}$判断$l_i$与新加入的直线L的交点若他在$l_{i-1}$的右侧则$l_i$能够被剔除

重复此操做

须要注意的是因为极角序是环形的，因此新加入的半平面可能绕了一圈把队首的半平面剔除，因此还要额外判断队首是否须要被剔除

习题：

bzoj P1007==luogu P3194 因为极角范围使得这个题比模板简单

3.4.对圆的处理：

圆还须要介绍吗？

3.4.1.圆的表示：

用圆心C和半径r能够肯定一个圆

引入参数方程能够用一个弧度来肯定圆上一个点

圆的参数方程：

$$P=C+(r~cos\theta,r~sin\theta)$$

struct Circle{
    Point c;
    double r;
    Circle (Point c,double r):c(c),r(r){    }
    Point poi(double cita){    return c+Point(cos(cita)*r,sin(cita)*r);    }
};//定义圆，引入圆的参数方程

Circle

（To be continue......）