计算几何

计算几何

Posted by JU on May 21, 2019

计算几何

注意

一定要好好注意两向量的旋转方向的判断!!

比如判断点与向量,向量与向量的位置/旋转关系要特别注意叉积是大于还是小于 $0$ .

向量

数量积

也叫点积/内积.

几何意义为一个向量在另一个向量上的投影再乘上第二个向量的模长.

(通过几何意义可以比较方便地求一些东西,比如最小矩形覆盖的左右边界啥的)

$\vec{a}\cdot\vec{b}=|\vec{a}||\vec{b}|\cos\theta=(x_1*x_2+y_1*y_2)$

$1.共线同向,\vec{a}\cdot\vec{b}=|\vec{a}||\vec{b}|$

$2.共线反向,\vec{a}\cdot\vec{b}=-|\vec{a}||\vec{b}|$

$3.\theta<90^\circ,\vec{a}\cdot\vec{b}> 0$

$4.\theta=90^\circ,\vec{a}\cdot\vec{b}=0$

$3.\theta>90^\circ,\vec{a}\cdot\vec{b}< 0$

外积

也叫叉积.

$\vec{a}\times\vec{b}=|\vec{a}||\vec{b}|\sin\theta=x_1*y_2-y_1*x_2$

$\theta$ 是将 $\vec{a}$ 逆时针转向 $\vec{b}$ 的角度大小.

可以求向量的旋转关系和点到直线距离.

两个共起点向量的叉积(的绝对值)在数值上等于它们构成的平行四边形的面积.

线段交点

对于两条分别由两个向量共四个向量构成的(有向)线段,它们的交点为:

\[S_{\triangle ACD}/(S_{\triangle ACD}+S_{\triangle BCD})*(\vec{b}-\vec{a})+\vec{a}\]

直线交点

对于两条分别用两个向量表示的直线,它们的交点为:

po crosspoint (line a,line b)
{
    double k1,k2,t;
    k1=(b.b-a.a)*(a.b-a.a);
    k2=(a.b-a.a)*(b.a-a.a);
    t=k1/(k1+k2);
    return b.b+(b.a-b.b)*t;
}

这东西强到可以代替上面那个(大概).

(暂时没搞懂原理)

向量旋转

逆时针.

\[\vec{b}=(x\cos\theta-y\sin\theta,y\cos\theta+x\sin\theta)\]

多边形面积

\[S=\vec{a_n}\times\vec{a_1}+\sum_{i=1}^{n-1}\vec{a_i}\times\vec{a_{i+1}}\]

模板

struct po
{
	db x,y;
	po (db xx=0,db yy=0) { x=xx,y=yy; }
};
po operator + (po a,po b) { return po(a.x+b.x,a.y+b.y); }
po operator - (po a,po b) { return po(a.x-b.x,a.y-b.y); }
po operator * (po a,db b) { return po(a.x*b,a.y*b); }
po operator * (db b,po a) { return po(a.x*b,a.y*b); }
db operator * (po a,po b) { return a.x*b.y-a.y*b.x; }
db operator / (po a,po b) { return a.x*b.x+a.y*b.y; }

凸包

性质

逆时针标记凸包上的点, 并设 $\vec{b_i}=\vec{a_{i+1}}-\vec{a_i}$ ,可以发现 $\vec{b_i}$ 是不断向左旋转的.

凸包求法

于是先将凸包按横坐标排序,升序枚举求下凸壳,再降序求上凸壳,再合并(比较傻).

模板

sort (a+1,a+n+1,ly);
for (int i=1; i<=n; i++)
{
	while (z>=2&&(a[st[z]]-a[st[z-1]])*(a[i]-a[st[z]])<=0) --z;
	st[++z]=i;
}
for (int i=1; i<=z; i++)
	t[++tot]=a[st[i]];
z=0;
for (int i=n; i>=1; i--)
{
	while (z>=2&&(a[st[z]]-a[st[z-1]])*(a[i]-a[st[z]])<=0) --z;
	st[++z]=i;
}
for (int i=2; i< z; i++)
	t[++tot]=a[st[i]];

注意一个点(除了边界)可能重复出现在凸包上(比如所有点在一条直线上).

旋转卡壳

对踵点

在凸多边形上两点作平行线,若有一种作法使得整个多边形都在两条线之间,则这两点是一对对踵点.

凸多边形的直径一定在对踵点中产生.

旋转卡壳过程

对于某条边,这条边两端的点具有同一个对踵点(也可能还有其他的).

一条边的对踵点,实际上就是离该边最远的点.

枚举每条边,记录当前对踵点位置,动态改变.

模板

t[tot+1]=t[1];
for (int i=1; i<=tot; i++)
{
	po B=t[i+1]-t[i];
	while (B*(t[p+1]-t[i])-B*(t[p]-t[i])> -eps) p=p==tot?1:p+1;
	//此处的-eps有时可改为0,否则可能会死循环
	ans=max (ans,max (getlen (t[p]-t[i]),getlen (t[p]-t[+1])));
}

在卡壳的时候还可以维护其他东西,比如最小覆盖矩形的左右边界啥的.

若只求单条边的对踵点,由于点到该边的距离是单峰的,所以可以三分求.

最小圆覆盖

为了优雅的做出 $GDSOI2019D2T2$ …

二元一次方程解法

\[\begin{cases}a_1x + b_1y=c_1\\a_2x + b_2y=c_2\end{cases}\] \[\begin{cases}x=\frac{b_1c_2-b_2c_1}{b_1a_2-b_2a_1}\\y=\frac{a_1c_2-a_2c_1}{a_1b_2-a_2b_1}\end{cases}\]

三点定圆

已知 $\vec{a}=(x_1,y_1),\vec{b}=(x_2,y_2),\vec{c}=(x_3,y_3)$ .
\(\begin{cases}(x_1-x)^2+(y_1-y)^2=r^2\\(x_2-x)^2+(y_2-y)^2=r^2\\(x_3-x)^2+(y_3-y)^2=r^2\\\end{cases}\)

$②-①,③-①$ ,并化简得:

\[\begin{cases}(x_2-x_1)x+(y_2-y_1)y=\frac{x_2^2-x_1^2+y_2^2+y_1^2}{2}\\(x_3-x_1)x+(y_3-y_1)y=\frac{x_3^2-x_1^2+y_3^2+y_1^2}{2}\end{cases}\]

对应到上面二元一次方程中:

db S(db a) { return a*a; }
po circle (po a,po b,po c)
{
	db a1=b.x-a.x,a2=c.x-a.x,b1=b.y-a.y,b2=c.y-a.y;
	db c1=(S(b.x)-S(a.x)+S(b.y)-S(a.y))/2;
	db c2=(S(c.x)-S(a.x)+S(c.y)-S(a.y))/2;
	return po((b1*c2-b2*c1)/(b1*a2-b2*a1),(a1*c2-a2*c1)/(a1*b2-a2*b1));
}

随机增量法

过程

若点 $p$ 不在点集 $S$ 的最小覆盖圆内,则一定在 ${S|p}$ 的最小覆盖圆上.

于是可以如此确定最小覆盖圆:

$1.$ 令前 $i-1$ 个点的最小覆盖圆为 $C$ .

$2.$ 若点 $i$ 在 $C$ 内(或 $C$ 上),则忽略该点.

$3.$ 点 $i$ 在最小覆盖圆上,另设一个圆 $G$ ,开始时仅有点 $i$ .

$4.$ 枚举点 $j=1\sim i-1$ ,若在 $G$ 外,则用点 $i,j$ 构成一个新圆.

$5.$ 枚举点 $k=1\sim j-1$ ,若在 $G$ 外,则用点 $i,j,k$ 构成一个新圆.

时间复杂度

这个算法看起来是 $O(n^3)$ 的,实际也是 $O(n^3)$ …

然鹅只需要在算法开始前随机打乱点的顺序,它的期望时间就变成了 $O(n)$ .

模板

random_shuffle (a+1,a+n+1);
po cen=0.5*(a[1]+a[2]); db r=getlen (cen-a[1]);
for (int i=3; i<=n; i++)
{
	if (getlen (cen-a[i])<=r) continue;
	cen=a[i],r=0;
	for (int j=1; j< i; j++) if (getlen (cen-a[j])> r)
	{
		cen=0.5*(a[i]+a[j]),r=getlen (cen-a[i]);
		for (int k=1; k< j; k++) if (getlen (cen-a[k])> r)
			cen=circle (a[i],a[j],a[k]),r=getlen (cen-a[i]);
	}
}

半平面交

半平面就是直线某一侧的部分,半平面交就是一堆这些部分的交集(雾).

对于用向量表示的直线,半平面一般指向量左侧的部分.

下面说”直线的某侧”实际是指向量的某侧.

极角排序

要先对直线极角排序.

先对每条直线都求一下 $atan2$ ,再按其为关键字排序,相同的只保留最左边.

bool ly (line a,line b)
{
	if (a.slp!=b.slp) return a.slp< b.slp;
	return (a.b-a.a)*(b.b-a.a)< 0;
}
int main()
{
    for (int i=1; i<=n; i++)
        l[i].slp=atan2 (l[i].b.y-l[i].a.y,l[i].b.x-l[i].a.x);
    sort (l+1,l+n+1,ly);
    tot=0;
    for (int i=1; i<=n; i++)
    if (l[i].slp!=l[i-1].slp) l[++tot]=l[i];
}

排序是按照某种方便处理的顺序,查查 $atan2$ 是啥就知道了.

求半平面交

需要维护一个双端队列.

当新加入一条直线时:

$1.$ 对于队尾的两条直线,若它们的交点在当前直线的右侧,则该点没有用,弹出队尾(重复).

$2.$ 对于队首的两条直线,若它们的交点在当前直线的右侧,则该点没有用,弹出队首(重复).

(好像一毛一样)

最后再用队首约束队尾,队尾约束队首,使得队首队尾也合法.

bool check (line a,line b,line g)//a和b的交点是否在c的右侧
{
	po p=crosspoint (a,b);
	return (g.b-g.a)*(p-g.a)< 0;
}
int main()
{
    st=1,ed=0;
    for (int i=1; i<=n; i++)
    {
        while (st< ed&&check (q[ed-1],q[ed],l[i])) --ed;
        while (st< ed&&check (q[st],q[st+1],l[i])) ++st;
        q[++ed]=l[i];
    }
    while (st< ed&&check (q[ed-1],q[ed],q[st])) --ed;
    while (st< ed&&check (q[st],q[st+1],q[ed])) ++st;
}

好像由于某些玄学原因,需要先用队尾约束队首.

半平面交是一个凸包,于是可以算出凸包上的点,方便进行其他操作.

q[ed+1]=q[st];
tot=0;
for (int i=st; i<=ed; i++)
	t[++tot]=crosspoint (q[i],q[i+1]);

如果要求多边形面积,要特判 $n<=2$ ,否则可能爆炸.

CC 原创文章采用CC BY-NC-SA 4.0协议进行许可，转载请注明：“转载自：计算几何”
CATALOG
    FRIENDS