# 计算几何学习笔记

### 基础定义

#### 1.精度控制

const double eps=1e-7;
int dcmp(double x)
{
if (fabs(x)<eps) return 0;
else
if (x>0)
return 1;
else
return -1;
}

#### 2.计算$$π$$

const double pi=acos(-1.0);

#### 3.定义向量

struct Vector
{
double x,y;
Vector (double a=0,double b=0)
{
x=a,y=b;
}
};

#### 4.定义点

typedef Vector Point;

#### 5.定义直线

struct Line
{
Point p;
Vector v;
Line(Point P=Point(0,0),Vector V=Vector(0,0))
{
p=P,v=V;
}
};

### 向量

#### 1.四则运算

• 加法：横纵坐标分别相加
• 减法：横纵坐标分别相减
• 数乘：横纵坐标分别乘以乘数
• 除法：横纵坐标分别除以除数
• 相等：横纵左边分别相等
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 b) {return Vector(a.x*b  ,a.y*b  );}
Vector operator / (Vector a,double b) {return Vector(a.x/b  ,a.y/b  );}
bool   operator ==(Vector a,Vector b) {return !dcmp(a.x-b.x)&&!dcmp(a.y-b.y);}

#### 2.旋转

$sin(a+b)=sina*cosb+sinb*cosa$

$cos(a+b)=cosa*cosb-sina*sinb$

Vector Rotate(Vector a,double rad)
{
}

#### 3.点积

$$a \cdot b$$的集合意义为a在b上的投影长度乘以b的长度

$$a \cdot b=|a||b|cos<a,b>$$

$a=(x_1,y_1),b=(x_2,y_2),a \cdot b=(x_1 * x_2,y_1 * y_2)$

$$V_1 \cdot V_2=h \times |V1|$$

$$\because$$    $$V_1=\frac{y_1}{\cos \alpha}\,V_2=\frac{y_2}{\cos \beta}$$

$$\therefore$$    $$V_1 \cdot V_2=\frac{y_2}{\cos \beta } \times \cos (\beta – \alpha ) \times \frac{y_1}{\cos \alpha }$$

$$V_1 \cdot V_2=\frac{y_1y_2}{\cos \beta \cos \alpha } \times \cos (\beta – \alpha )$$

$$\because$$    $$\cos (\beta – \alpha )= \cos \alpha \cos \beta + \sin \alpha \sin \beta$$

$$\therefore$$    $$V_1 \cdot V_2=y_1y_2 + \sin \alpha \sin \beta \frac{y_1y_2}{\cos \beta \cos \alpha }$$

$$\because$$    $$\frac{\sin \alpha }{\cos \alpha }=\tan \alpha$$

$$\therefore$$    $$V_1 \cdot V_2=y_1y_2 + \tan \alpha \tan \beta y_1y_2$$

$$\therefore$$    $$V_1 \cdot V_2=y_1y_2 + x_1x_2$$

$$证毕$$

double Dot(Vector a,Vector b)
{
return a.x*b.x+a.y*b.y;
}

• 点积满足交换律
• 两个向量夹角是一个锐角时，点积为正数；为直角时点积为零（因此可用点积判垂直）；为钝角时点积为负数

double Len(Vector a){return sqrt(Dot(a,a));}

double Angle(Vector a,Vector b){return acos(Dot(a,b)/Len(a)/Len(b));}

#### 4.叉积

$$a \times b$$的几何意义就是ab所形成的平行四边形的有向面积。

$a=(x_1,y_1),b=(x_2,y_2),a \times b=x_1 * y_2 – x_2 * y_1$

double Cross(Vector a,Vector b)
{
return a.x*b.y-a.y*b.x;
}

poj 2318
poj 2398

#### 5.单位向量与单位法向量

$$\frac{a}{|a|}$$为a的单位向量，方向相同，长度为1

Vector Normal(Vector a)
{
double L=Len(a);
return Vector(-a.y/L,a.x/L);
}

### 点与直线

poj 1269

#### 2.直线交点

$$t_1=\frac{cross(w,u)}{cross(v,w)},t_2=\frac{cross(v,u)}{cross(v,w)}$$

Point line_cross(Line a,Line b)
{
Vector u=a.p-b.p;
double t=cross(b.v,u)/cross(a.v,b.v);
return a.p+t*a.v;
}

#### 4.点到直线的距离

double Dis_To_Line(Point p,Line l)
{
return fabs(Cross(l.v-l.p,p-l.p)/Len(l.v-l.p));
}

#### 5.点到线段的距离

double Dis_To_Seg(Point p,Point a,Point b)
{
if (a==b) return Len(p-a);
Vector v=b-a,w=p-a,u=p-b;
if (dcmp(Dot(v,w)<0))	return Len(w);
else if (dcmp(Dot(v,u))>0)	return Len(v);
else return fabs(Cross(v,w)/Len(v));
}

#### 6.点到直线的投影

Point Line_Pro(Point p,Point a,Point b)
{
Vector v=b-a;
return a+v*(Dot(v,p-a)/Dot(v,v));
}

#### 7.判断直线和线段是否相交

bool Line_Seg_Cross(Line l,Point a,Point b)
{
Vector v=l.p-l.v,w=a-l.p,u=b-l.p;
return dcmp(Cross(v,w))!=dcmp(Cross(v,u));
}

poj 1269

#### 8.判断线段和线段是否相交

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


poj 2653
poj 2826
poj 1410

### 三角形

#### 1.三角形面积

double Area_T(point a, point b,point c)
{
return Cross(b-a,c-a)*0.5;
}

#### 2.求三角形的外心

Point TC(Point a,Point b,Point c)
{
Point p=(a+b)/2,q=(a+c)/2;
Vector v=Rotate(b-a,pi/2),w=Rotate(c-a,pi/2);
if (!dcmp(Cross(v,w)))
{
if (dcmp(Len(a-b)+Len(b-c)-Len(a-c))==0)
return (a+c)/2;
if (dcmp(Len(a-c)+Len(b-c)-Len(a-b))==0)
return (a+b)/2;
if (dcmp(Len(a-b)+Len(a-c)-Len(b-c))==0)
return (b+c)/2;
}
return Line_cross(Line(p,v),Line(q,w));
}

#### 3.求三角形的重心

Point TG(Point a,Point b,Point c)
{
return (a+b+c)/3;
}

### 多边形

#### 1.多边形的面积

double Area_P(Point ploy[],int num)
{
double area=0;
for (int i=2;i<num;i++)
area+=Cross(ploy[i]-ploy[1],ploy[i+1]-ploy[1]);
return area/2;
}

poj 1654

#### 2.点是否在多边形内

bool Point_In_Ploy(Point p,Point poly[],int num)
{
int ans=0,k,d1,d2;
for (int i=1;i<=num;i++)
{
if (!dcmp(Dis_To_Seg(p,poly[i],poly[i%num+1]))) return 0;
k=dcmp(Cross(poly[i%num+1]-poly[i],p-poly[i]));
d1=dcmp(poly[i].y-p.y);
d2=dcmp(poly[i%num+1].y-p.y);
if (k>0&&d1<=0&&d2>0) ++ans;
if (k<0&&d1>=0&&d2<0) --ans;
}
if (ans) return 1;
return 0;
}

poj 1584

### 凸包

#### 2.计算流程

1. 将点按x坐标第一关键字，y坐标第二关键字排序
2. 首先将p1,p2压入栈中
3. 判断下一个点和栈中最后一个点组成的向量是否在当前前进方向的左边（求下凸壳），如果是的话就入栈，否则弹栈
4. 这样求出了下凸壳，再倒着做一遍就求出来了正凸壳，合起来就是完整的凸包

COGS 896
pku 1113
pku 1228
pku 1873
pku 3348

#### 3.凸包应用：旋转卡壳

1.

2.

Compute the polygon’s extreme points in the y direction. Call them ymin and ymax. Construct two horizontal lines of support through ymin and ymax. Since this is already an anti-podal pair, compute the distance, and keep as maximum. Rotate the lines until one is flush with an edge of the polygon. A new anti-podal pair is determined. Compute the new distance, compare to old maximum, and update if necessary. Repeat steps 3 and 4 until the anti-podal pair considered is (ymin,ymax) again. Output the pair(s) determining the maximum as the diameter pair(s).

double rotating()
{
if (top==1) return 0;
if (top==2) return Len(stack[2]-stack[1]);
int now=1;
double ans=0;
for (int i=1;i<=top;++i)
{
while (dcmp(Dis_To_Line(stack[now],stack[i],stack[i%top+1])-Dis_To_Line(stack[now%top+1],stack[i],stack[i%top+1]))<=0)
now=now%top+1;
ans=max(ans,Len(stack[now]-stack[i]));
ans=max(ans,Len(stack[now]-stack[i%top+1]));
}
return ans;
}

poj 2187
poj 2079
bzoj 1069
bzoj 1185

### 半平面交

#### 5.如何求半平面交呢

Point poly[N],npoly[N];int cnt,ncnt;
void Crate_Square()
{
cnt=0;
poly[++cnt]=Point(INF,INF);
poly[++cnt]=Point(INF,-INF);
poly[++cnt]=Point(-INF,-INF);
poly[++cnt]=Point(-INF,INF);
}
void Cut_Ploy(Point a,Point b)
{
ncnt=0;
Point c,d;
for (int i=1;i<=cnt;i++)
{
c=poly[i%cnt+1];
d=poly[(i+1)%cnt+1];
if (dcmp(Cross(b-a,c-a))>=0)
npoly[++ncnt]=c;
if (Line_Seg_Cross(a,b,c,d))
npoly[++ncnt]=Line_cross(Line(a,b-a),Line(c,d-c));
}
cnt=ncnt;
for (int i=1;i<=ncnt;i++)
poly[i]=npoly[i];
}

pku 3130
pku 1474
pku 1279
pku 3525
pku 3384