计算几何基础
计算几何学了一些，感觉学的太乱，整理一下。
文章目录
精度:
由于使用浮点数，所以会产生精度问题，所以在比较大小的时候不能直接采用等于。
const double eps=1e-8;
int sgn(double x)
{
if(fabs(x) < eps)return 0;
if(x < 0) return -1;
else return 1;
}
三角函数和角度转换：
c语言中自带三角函数
s
i
n
(
x
)
,
c
o
s
(
x
)
,
t
a
n
(
x
)
。
x
sin(x),cos(x),tan(x)。x
sin(x),cos(x),tan(x)。x是以弧度为单位。
假设输入得角度为
θ
°
\theta°
θ°，那么：
x
=
θ
π
180
x=\frac{\theta\pi}{180}
x=180θπ为了精确，采用：
π
=
a
c
o
s
(
−
1
)
\pi=acos(-1)
π=acos(−1)。
还有经常使用得一个函数atan2(double y,double x)。返回值为
θ
=
a
r
c
t
a
n
y
x
\begin{aligned}\theta =arctan\frac{y}{x}\end{aligned}
θ=arctanxy。
点、线、向量:
struct point
{
double x,y;
point(){}
point(double xx,double yy){x=xx;y=yy;}
point operator +(point b){return point(x+b.x,y+b.y);}
point operator -(point b){return point(x-b.x,y-b.y);}
double operator ^(point b){return x*b.y-y*b.x;}//叉乘
double operator *(point b){return x*b.x+y*b.y;}//点乘
point operator *(double b){return point(x*b,y*b);}//数乘
double gettan(){return atan2(y,x);}//角度
};
struct line
{
point s,e,l;
line(){}
line(point ss,point ee){s = ss;e = ee;l=e-s;}
double gettan(){return atan2(e.y-s.y,e.x-s.x);}
};
叉乘（叉积）：
a
⃗
×
b
⃗
=
∣
a
⃗
∣
∣
b
⃗
∣
s
i
n
θ
=
x
a
×
y
b
−
y
a
×
x
b
\vec a\times\vec b=|\vec a||\vec b|sin\theta=x_a\times y_b-y_a\times x_b
a×b=∣a∣∣b∣sinθ=xa×yb−ya×xb
如果
b
⃗
\vec b
b在
a
⃗
\vec a
a的逆时针方向，值为正；顺时针方向值为负；共线为0。
叉积也可以用来表示这两条向量构成的平行四边形的有向面积。
hdu多校第八场1003Clockwise or Counterclockwise就是用叉积判断 A B → \overrightarrow {AB} AB和 A C → \overrightarrow {AC} AC的位置关系。
向量旋转：
//a为弧度
point spin(point p,double a)
{
return point(p.x*cos(a)-p.y*sin(a),p.y*cos(a)+p.x*sin(a));
}
设
a
⃗
=
(
x
,
y
)
\vec a=(x,y)
a=(x,y),倾斜角为
θ
\theta
θ，则长度为
l
=
x
2
+
y
2
l=\sqrt{x^2+y^2}
l=x2+y2,所以
x
=
l
c
o
s
θ
,
y
=
l
s
i
n
θ
x=lcos\theta,y=lsin\theta
x=lcosθ,y=lsinθ。逆时针旋转
α
\alpha
α度角得到
b
⃗
\vec b
b。
b
⃗
=
(
l
c
o
s
(
θ
+
α
)
,
l
s
i
n
(
θ
+
α
)
)
=
(
l
(
c
o
s
θ
c
o
s
α
−
s
i
n
θ
s
i
n
α
)
,
l
(
s
i
n
θ
c
o
s
α
+
c
o
s
θ
s
i
n
α
)
)
=
(
l
c
o
s
θ
c
o
s
α
−
l
s
i
n
θ
s
i
n
α
,
l
s
i
n
θ
c
o
s
α
+
l
c
o
s
θ
s
i
n
α
)
=
(
x
c
o
s
α
−
y
s
i
n
α
,
y
c
o
s
α
+
x
s
i
n
α
)
\begin{aligned} \vec b&=(lcos(\theta+\alpha),lsin(\theta+\alpha))\\&=(l(cos\theta cos\alpha -sin\theta sin\alpha),l(sin\theta cos\alpha+cos\theta sin\alpha ))\\&=(lcos\theta cos\alpha -lsin\theta sin\alpha,lsin\theta cos\alpha+lcos\theta sin\alpha)\\&=(xcos\alpha -ysin\alpha,ycos\alpha+xsin\alpha)\end{aligned}
b=(lcos(θ+α),lsin(θ+α))=(l(cosθcosα−sinθsinα),l(sinθcosα+cosθsinα))=(lcosθcosα−lsinθsinα,lsinθcosα+lcosθsinα)=(xcosα−ysinα,ycosα+xsinα)
点到直线的投影点：
point projection(line l,point p)
{
return l.s+l.l*((l.l*(p-l.s))/((l.l)*(l.l)));
}
用点积来求。
设
C
C
C到直线
A
B
AB
AB的投影点为
C
1
C_1
C1。
∵
A
B
→
⋅
A
C
→
=
∣
A
B
→
∣
∣
A
C
→
∣
c
o
s
θ
=
∣
A
B
→
∣
∣
A
C
1
→
∣
\because \overrightarrow{AB}\cdot\overrightarrow{AC}=|\overrightarrow{AB}||\overrightarrow{AC}|cos\theta=|\overrightarrow{AB}||\overrightarrow{AC_1}|
∵AB⋅AC=∣AB∣∣AC∣cosθ=∣AB∣∣AC1∣
∴
∣
A
C
1
→
∣
=
A
B
→
⋅
A
C
→
∣
A
B
→
∣
\begin{aligned}\therefore |\overrightarrow{AC_1}|=\frac{\overrightarrow{AB}\cdot\overrightarrow{AC}}{|\overrightarrow{AB}|}\end{aligned}
∴∣AC1∣=∣AB∣AB⋅AC
∵
O
C
1
→
=
O
A
→
+
A
C
1
→
A
C
1
→
=
∣
A
C
1
→
∣
∣
A
B
→
∣
×
A
B
→
\begin{aligned}\because &\overrightarrow{OC_1}=\overrightarrow{OA}+\overrightarrow{AC_1}\\&\overrightarrow{AC_1}=\frac{|\overrightarrow{AC_1}|}{|\overrightarrow{AB}|}\times \overrightarrow{AB}\end{aligned}
∵OC1=OA+AC1AC1=∣AB∣∣AC1∣×AB
∴
O
C
1
→
=
O
A
→
+
∣
A
C
1
→
∣
∣
A
B
→
∣
×
A
B
→
=
O
A
→
+
A
B
→
⋅
A
C
→
∣
A
B
→
∣
∣
A
B
→
∣
×
A
B
→
=
O
A
→
+
A
B
→
⋅
A
C
→
∣
A
B
→
∣
2
×
A
B
→
\begin{aligned}\therefore \overrightarrow{OC_1}=\overrightarrow{OA}+\frac{|\overrightarrow{AC_1}|}{|\overrightarrow{AB}|}\times \overrightarrow{AB}=\overrightarrow{OA}+\frac{\frac{\overrightarrow{AB}\cdot\overrightarrow{AC}}{|\overrightarrow{AB}|}}{|\overrightarrow{AB}|}\times\overrightarrow{AB} =\overrightarrow{OA}+\frac{\overrightarrow{AB}\cdot\overrightarrow{AC}}{|\overrightarrow{AB}|^2}\times\overrightarrow{AB}\end{aligned}
∴OC1=OA+∣AB∣∣AC1∣×AB=OA+∣AB∣∣AB∣AB⋅AC×AB=OA+∣AB∣2AB⋅AC×AB
计算两点间距离：
double dis(point a,point b)
{
return sqrt((a-b)*(a-b));
}
判断点在线段上：
bool onseg(line l,point p)
{
bool x1=sgn((l.s-p)^(l.e-p))==0;
bool x2=sgn((l.s-p)*(l.e-p))<=0;
return x1&&x2;
}
因为点 P P P在线段 A B AB AB上,由叉积的性质可以得到 P A → × P B → = 0 \overrightarrow{PA}\times\overrightarrow{PB}=0 PA×PB=0。但是如果 P P P在 A B AB AB的延长线上，这也符合。所以要看 P A → ⋅ P B → < = 0 \overrightarrow{PA}\cdot\overrightarrow{PB}<=0 PA⋅PB<=0是否成立，因为此时 ∠ A P B = 180 ° , c o s 180 ° = − 1 ∠APB=180°,cos180°=-1 ∠APB=180°,cos180°=−1,所以点乘的值为负。
计算点到直线距离：
double dis_pl(point p,line l)
{
point v=p-l.s;
double ans=(v^l.l)/sqrt(l.l*l.l);
return fabs(ans);
}
以前都学过
(
x
0
,
y
0
)
(x_0,y_0)
(x0,y0)到
A
x
+
B
y
+
C
=
0
Ax+By+C=0
Ax+By+C=0的距离为
d
=
∣
A
x
0
+
B
y
0
+
C
∣
A
2
+
B
2
\begin{aligned}d=\frac{|Ax_0+By_0+C|}{\sqrt{A^2+B^2}}\end{aligned}
d=A2+B2∣Ax0+By0+C∣。但是这么计算比较复杂，我们还是通过叉乘和点乘来计算。
我们要求点
P
P
P到直线
A
B
AB
AB的距离。
∵
S
Δ
A
B
P
=
∣
A
P
→
×
A
B
→
2
∣
=
d
×
∣
A
B
→
∣
2
\begin{aligned}\because S_{\Delta ABP}=\left|\frac{\overrightarrow{AP}\times\overrightarrow{AB}}{2}\right|=\frac{d\times |\overrightarrow{AB}|}{2}\end{aligned}
∵SΔABP=∣∣∣∣∣2AP×AB∣∣∣∣∣=2d×∣AB∣
∴
d
=
∣
A
P
→
×
A
B
→
∣
A
B
→
∣
∣
\begin{aligned}\therefore d=\left|\frac{\overrightarrow{AP}\times\overrightarrow{AB}}{|\overrightarrow{AB}|}\right|\end{aligned}
∴d=∣∣∣∣∣∣AB∣AP×AB∣∣∣∣∣
如图，
2
S
Δ
A
B
P
=
∣
A
P
→
×
A
B
→
∣
=
∣
(
−
2
,
2
)
×
(
2
,
1
)
∣
=
∣
−
6
∣
=
6
\begin{aligned}2S_{\Delta ABP}=|\overrightarrow{AP}\times\overrightarrow{AB}|=|(-2,2)\times(2,1)|=|-6|=6\end{aligned}
2SΔABP=∣AP×AB∣=∣(−2,2)×(2,1)∣=∣−6∣=6
∣
A
B
→
∣
=
x
A
B
2
+
y
A
B
2
=
5
≈
2.23607
\begin{aligned}|\overrightarrow{AB}|=\sqrt{x_{AB}^2+y_{AB}^2}=\sqrt{5}\approx2.23607\end{aligned}
∣AB∣=xAB2+yAB2=5≈2.23607
d
=
∣
A
P
→
×
A
B
→
∣
A
B
→
∣
∣
=
6
5
≈
2.82843
\begin{aligned}d=\left|\frac{\overrightarrow{AP}\times\overrightarrow{AB}}{|\overrightarrow{AB}|}\right|=\frac{6}{\sqrt{5}}\approx2.82843\end{aligned}
d=∣∣∣∣∣∣AB∣AP×AB∣∣∣∣∣=56≈2.82843
判断线段相交：
bool judge(line l1,line l2)
{
point a=l1.s,b=l1.e,c=l2.s,d=l2.e;
double x1=((c-a)^l1.l)*((d-a)^l1.l);
double x2=((a-c)^l2.l)*((b-c)^l2.l);
x1&&x2&&
max(l1.s.x,l1.e.x) >= min(l2.s.x,l2.e.x) &&
max(l2.s.x,l2.e.x) >= min(l1.s.x,l1.e.x) &&
max(l1.s.y,l1.e.y) >= min(l2.s.y,l2.e.y) &&
max(l2.s.y,l2.e.y) >= min(l1.s.y,l1.e.y);
}
如果有两个线段AB和CD，判断相交就是判断一条线段的两个端点是否在另外一条线段的两端。也就是判断A，B是否在CD两端以及C，D是否在AB两端。
如果A，B在CD两端，那么
(
C
A
→
×
C
D
→
)
⋅
(
C
B
→
×
C
D
→
)
<
0
(\overrightarrow {CA}\times\overrightarrow{CD})\cdot(\overrightarrow{CB}\times\overrightarrow{CD})<0
(CA×CD)⋅(CB×CD)<0,因为如果A，B在CD两侧，那么CA，CB肯定分别在CD的顺时针和逆时针方向，这时候用叉积即可判断。由于平行或者点重合的情况存在，所以要特判一下叉积为0时，
C
,
D
C,D
C,D点是否在
A
B
AB
AB上，或者
A
,
B
A,B
A,B点是否在
C
D
CD
CD上。
题目：HDU1086
计算两直线交点:
point getcross(line l1,line l2)
{
double s1=((l2.s-l1.s)^(l2.e-l1.s));
double s2=((l2.e-l1.e)^(l2.s-l1.e));
if(s1+s2==0)//平行
return point{inf,inf};
return l1.s+l1.l*(s1/(s1+s2));
}
如图：
设
A
B
AB
AB与
C
D
CD
CD交点为
E
E
E。
用叉积计算出
S
Δ
A
C
D
S_{\Delta ACD}
SΔACD和
S
Δ
B
C
D
S_{\Delta BCD}
SΔBCD的有向面积的值
S
1
,
S
2
S_1,S_2
S1,S2。做垂线
A
H
1
,
B
H
2
AH_1,BH_2
AH1,BH2。
易得：
S
1
=
∣
C
D
∣
×
∣
A
H
1
∣
2
,
S
2
=
∣
C
D
∣
×
∣
B
H
2
∣
2
\begin{aligned}S_1=\frac{|CD|\times|AH_1|}{2},S_2=\frac{|CD|\times|BH_2|}{2}\end{aligned}
S1=2∣CD∣×∣AH1∣,S2=2∣CD∣×∣BH2∣
∴
S
1
S
2
=
A
H
1
B
H
2
\begin{aligned}\therefore \frac{S_1}{S_2}=\frac{AH_1}{BH_2}\end{aligned}
∴S2S1=BH2AH1
易证：
Δ
A
H
1
E
∼
Δ
B
H
2
E
\Delta AH_1E\sim\Delta BH_2E
ΔAH1E∼ΔBH2E
∴
S
1
S
1
+
S
2
=
A
H
1
A
H
1
+
B
H
2
=
A
E
A
E
+
B
E
=
A
E
A
B
\begin{aligned}\therefore\frac{S_1}{S_1+S_2}=\frac{AH_1}{AH_1+BH_2}=\frac{AE} {AE+BE}=\frac{AE}{AB}\end{aligned}
∴S1+S2S1=AH1+BH2AH1=AE+BEAE=ABAE
∴
O
E
→
=
O
A
→
+
S
1
S
1
+
S
2
A
B
→
\begin{aligned}\therefore\overrightarrow{OE}=\overrightarrow{OA}+\frac{S_1}{S_1+S_2}\overrightarrow{AB}\end{aligned}
∴OE=OA+S1+S2S1AB
这样就得到了
E
E
E的坐标了。
三角形面积：
- 两边叉积除以2取绝对值。
- 海伦公式：
S = p ( p − a ) ( p − b ) ( p − c ) S=\sqrt{p(p-a)(p-b)(p-c)} S=p(p−a)(p−b)(p−c)
p = a + b + c 2 \begin{aligned}p=\frac{a+b+c}{2}\end{aligned} p=2a+b+c
Pick定理：
如果一个多边形的顶点在格点上，那么他的面积可以表示为:
S
=
a
+
b
2
−
1
S=a+\frac{b}{2}-1
S=a+2b−1
a
a
a表示多边形内部的点数，
b
b
b表示多边形落在格点边界上的点数。
多边形面积：
point d[1000];
int n,i;
scanf("%d",&n);
for(i=1;i<=n;i++)
scanf("%lf%lf",&d[i].x,&d[i].y);
d[n+1]=d[1];
double area=0;
for(i=1;i<=n;i++)
area+=(d[i]^d[i+1]);
area/=2;
printf("%.2lf",fabs(area));
利用叉积（三角形剖分）进行计算。
如下图，
S
A
B
C
D
E
F
=
S
Δ
O
A
F
+
S
Δ
O
B
A
+
S
Δ
O
C
B
+
S
Δ
O
D
C
−
S
Δ
O
E
D
−
S
Δ
O
F
E
S_{ABCDEF}=S_{\Delta OAF}+S_{\Delta OBA}+S_{\Delta OCB}+S_{\Delta ODC}-S_{\Delta OED}-S_{\Delta OFE}
SABCDEF=SΔOAF+SΔOBA+SΔOCB+SΔODC−SΔOED−SΔOFE
由于叉积有方向，我们可以发现叉积之和除以2就是多边形的面积。这样我们只要逆时针求相邻点的叉积即可,多边形面积就可以表示为：
S
=
O
A
n
→
×
O
A
1
→
+
∑
i
=
1
n
−
1
O
A
i
→
×
O
A
i
+
1
→
2
\begin{aligned}S=\frac{\overrightarrow{OA_n}\times\overrightarrow{OA_1}+\sum_{i=1}^{n-1}\overrightarrow{OA_i}\times\overrightarrow{OA_{i+1}}}{2}\end{aligned}
S=2OAn×OA1+∑i=1n−1OAi×OAi+1
题目：HDU2036
凸包：
凸包可以理解为包含所有点得凸多边形，或者形象得想成能包围所有点得橡皮筋，也可以是说一个周长最小，能包含住所有点得多边形。
时间复杂度
O
(
n
l
o
g
n
)
O(nlogn)
O(nlogn)。
//graham扫描法 O(nlogn)
#include<bits/stdc++.h>
#define ll long long
#define cl(x,y) memset(x,y,sizeof(x))
const int N=1e6+10;
const double eps=1e-8;
const double pi=acos(-1);
using namespace std;
int n,top;
int sgn(double x)
{
if(fabs(x) < eps)return 0;
if(x < 0) return -1;
else return 1;
}
struct point
{
double x,y;
point(){}
point(double xx,double yy){x=xx;y=yy;}
point operator +(point b){return point(x+b.x,y+b.y);}
point operator -(point b){return point(x-b.x,y-b.y);}
double operator ^(point b){return x*b.y-y*b.x;}//叉乘
double operator *(point b){return x*b.x+y*b.y;}//点乘
point operator *(double b){return point(x*b,y*b);}//数乘
void read(){scanf("%lf%lf",&x,&y);}
}p[N],sta[N];
int cmp1(point a,point b)
{
if(a.y==b.y)
return a.x<b.x;
else
return a.y<b.y;
}
int cmp2(point a,point b)//极角排序
{
if(atan2(a.y-sta[1].y,a.x-sta[1].x)!=atan2(b.y-sta[1].y,b.x-sta[1].x))
return (atan2(a.y-sta[1].y,a.x-sta[1].x))<(atan2(b.y-sta[1].y,b.x-sta[1].x));
return a.x<b.x;
}
double dis(point a,point b)
{
return sqrt((a-b)*(a-b));
}
void graham()
{
int i;
sort(p+1,p+n+1,cmp1);
sta[1]=p[1];
sort(p+2,p+n+1,cmp2);
sta[2]=p[2];
top=2;
for(i=3;i<=n;i++)
{
while(top>1 && ((sta[top]-sta[top-1])^(p[i]-sta[top-1]))<0)
top--;
sta[++top]=p[i];
}
return;
}
int main()
{
int i;
scanf("%d",&n);
for(i=1;i<=n;i++)
p[i].read();
graham();
//计算凸包周长
double s=0;
for(i=2;i<=top;i++)
s+=dis(sta[i],sta[i-1]);
s+=dis(sta[top],sta[1]);
//计算凸包面积
double area=0;
for(i=3;i<=top;i++)
area+=((sta[i-1]-sta[1])^(sta[i]-sta[1]));
cout<<s<<endl<<area<<endl;
return 0;
}
旋转卡壳：
时间复杂度 O ( n ) O(n) O(n)。
求凸多边形直径（凸包上最远点对的距离）：
double rotating_calipers()
{
sta[top+1]=sta[1];
int now=3,i;
double ans=0;
for(i=1;i<=top;i++)
{
while(((sta[i+1]-sta[i])^(sta[now]-sta[i]))<((sta[i+1]-sta[i])^(sta[now+1]-sta[i])))
{
now++;
if(now==top+1)
now=1;
}
ans=max(ans,dis(sta[now],sta[i]));
}
return ans;
}
求凸包的宽：
double rc()
{
sta[top+1]=sta[1];
int now=3,i;
double ans=INF;
for(i=1;i<=top;i++)
{
while(((sta[i+1]-sta[i])^(sta[now]-sta[i]))<((sta[i+1]-sta[i])^(sta[now+1]-sta[i])))
{
now++;
if(now==top+1)
now=1;
}
ans=min(ans,((sta[i+1]-sta[i])^(sta[now]-sta[i]))/dis(sta[i],sta[i+1]));
}
return ans;
}
半平面交：
半平面是指一条直线会把平面分成两半，我们一般把直线左边的半平面视为需要的半平面。
求多边形的核（从这个区域可以看到多边形的任意位置）。
思路：
- 极角排序。
- 去掉极角相同时，取左边的边。因为显然右边线的半平面包含左边线的半平面，加上左边的线右边的线肯定要被舍弃。
- 维护双向队列，看加入一条线之后，头尾部线相交的交点是否在这条线的右侧，如果在，弹出线。
- 得到的线的集合就是要求的半平面交。
int cmp(line a,line b)
{
double t1=a.gettan(),t2=b.gettan();
if(sgn(t1-t2)==0)
return (a.l^(b.e-a.s))>0;
return t1<t2;
}
//判断bc交点是否在a的右边，如果在，去掉头部或者尾部的线
int onright(line a,line b,line c)
{
point p=getcross(b,c);
return (a.l^(p-a.s))<0;
}
//半平面交
void hpi(line li[],int n)
{
sort(li+1,li+n+1,cmp);
int l=1,r=1,cnt=0,i;//模拟deque
line que[N];
//去重，极角相同取左边，即保留最后直线
for(i=1;i<n;i++)
if(sgn(li[i].gettan()-li[i+1].gettan()))
li[++cnt]=li[i];
li[++cnt]=li[n];
//判断新加入线的影响
for(i=1;i<=cnt;i++)
{
while(r-l>1 && onright(li[i],que[r-1],que[r-2]))
r--;
while(r-l>1 && onright(li[i],que[l],que[l+1]))
l++;
que[r++]=li[i];
}
//判断最先和最后加入的直线的影响
while(r-l>1 && onright(que[l],que[r-1],que[r-2]))
r--;
while(r-l>1 && onright(que[r-1],que[l],que[l+1]))
l++;
}
2071
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
