计算几何初步

约定:以下内容中\(*\)表示普通乘法,\(\cdot\)为点乘,\(\times\)为叉乘
代码中的*表示叉积,&表示点积, ^表示数乘
(当然是从xzy学长那儿学来的习惯啦)
文末有惊喜

前置知识

向量

向量基本概念与运算详见《2020人教版高中数学·选修一》

点积

又叫内积、数量积, 计算出来的结果是数值
几何意义: 第一个向量在第二个向量上的投影再乘上第二个向量的模长

其运算法则为:(\(\theta\)表示\(a\),\(b\)的夹角,\(0^\circ\le\theta\le180^\circ\)

\[\vec{a}\cdot\vec{b}=\left|\vec{a}\right|\left|\vec{b}\right|\cos\theta \]

\(\vec{a}(x_1,y_1),\vec{b}(x_2,y_2)\)

\[\vec{a}\cdot\vec{b}=x_1x_2+y_1y_2 \]

根据三角函数关系,易知:

  • \(\theta=0^\circ\),那么它们的点积为他们的模长之积。
  • \(\theta<90^\circ\),那么它们的点积为
  • \(\theta=90^\circ\),那么它们的点积为0
  • \(\theta>90^\circ\),那么它们的点积为

叉积

又称外积、向量积,计算出来的结果是向量(垂直于两向量所在平面)
运算法则为:

\[\vec{a}\times\vec{b}=|\vec{a}||\vec{b}|\sin\theta, ( 0^\circ\le\theta\le180^\circ) \]

\(\vec{a}(x_1,y_1),\vec{b}(x_2,y_2)\)

\[\vec{a}\times\vec{b}=x_1y_2-x_2y_1 \]

几何意义是两向量围成的有向面积, 因此叉积计算是没有交换律的
其方向可以通过右手法则来判定: 让\(\vec{a}\)穿过手心,然后四指卷曲从\(\vec{a}\)弯向\(\vec{b}\),大拇指伸直所指的方向就是叉积方向

那么有:

  • \(\vec{a}\times\vec{b}=0\)\(\vec{a}\)\(\vec{b}\) 平行
  • \(\vec{a}\times\vec{b}>0\)\(\vec{b}\)\(\vec{a}\) 的逆时针方向
  • \(\vec{a}\times\vec{b}<0\)\(\vec{b}\)\(\vec{a}\) 的顺时针方向
    (大拇指朝纸里为负,朝纸外为正)

极坐标系


如上图
选取平面上一定点\(O\) 作为极点, 从极点引出一条射线(上图中的x轴) 作为极轴
那么对于一个点\(B\),以极轴逆时针方向旋转到\(B\)的角度称为极角记为 \(\theta\)\(B\)\(O\)的距离称为极径记为 \(\rho\)
那么此时平面上的任何一个点都可以通过((极径)\(\rho\),(极角)\(\theta\))的方式表示

极坐标系和平面直角坐标系(笛卡尔坐标系)的转化
若点\(A(x,y)\)的极坐标为\((\rho,\theta)\), 那么有:

\[\begin{cases} x=\rho \cos \theta\\ y=\rho \sin \theta \end{cases} \]

又可推出:

\[\rho^2=x^2+y^2\\ \tan \theta=\frac{y}{x}\ \ \ \ (x\not =0) \]

那么极角 \(\theta=\arctan \frac{y}{x}\),在代码中可直接用 atan2(y, x)

点、线、角的相关运算

点与线

点P到直线 l 的距离

在直线\(l\)上任取两点A,B, P到\(l\)的距离D为:

\[D=\frac{\vec{AB}\times\vec{AP}}{\vec{|AB|}} \]

点P到线段 AB 的距离

若点P向\(l_{AB}\)引一条垂线,垂足不在线段\(AB\)上,则距离为 \(min(|PA|,|PB|)\)

线与线

判断线段是否相交

可以用快速排斥实验跨立实验
其中快速排斥实验大意是:若两个以给定带测验的线段为对角线的矩形,矩形边框不相交,则两线段一定不相交。否则,不一定相交

主要来讲讲跨立实验
思路:两条线段相交,则对于两条线段各自满足,另一条线段的两端一定分别位于该点两侧(也就是“跨立”)

对于跨立的判定:

上图中,观察到,一条线段CD若满足跨立,则\(\vec{CA}\times\vec{CD}\)的方向应与\(\vec{CD}\times\vec{CB}\)方向相同(注意先后顺序)
对应表达式为:

\[0\le|\vec{CA}\times\vec{CD}|*|\vec{CD}\times\vec{CB}| \]

但是不难发现上图中仅有CD满足跨立,而AB不满足,即:

\[0>|\vec{AC}\times\vec{AB}|*|\vec{AB}\times\vec{AD}| \]

(搞不清楚的话,可以通过右手定则找一找方向)

求两线段交点

若存在线段交点,则线段交点即为直线交点,那么只需判断是否有交,再找到直线交点就好了

直线交点

(从xzy学长的博客里嫖来的图)

不妨设点\(P\),可通过 \(A1+\lambda\vec{a}\) 得到
根据相似三角形(面积比—>线段比),则

\[\because \frac{S_{B_1B_2A_1E}}{S_{B_1B_2DC}}=\frac{A_1F}{CG},\ \triangle B_1CG \sim\triangle PA_1F \]

\[\therefore \frac{A_1F}{CG}=\frac{A_1P}{CB_1}=\lambda \]

即:

\[\lambda=\frac{\vec{b}\times\vec{c}}{\vec{b}\times\vec{a}} \]

向量旋转


(上图中 \(l^2=x^2+y^2\))
根据三角关系得:

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

(不记得了,也可通过图中式子经三角恒等变换推出)

极角排序

  • 法1:(时间更快)
    根据前文提到的 atan2(y,x) 这个函数的返回值(极角)作为键值,利用sort排序

  • 法2:(精度更高)
    利用叉积的正负性作为键值,利用sort排序

多边形相关

叉积求面积

计算方法:
按顺序将相邻顶点与原点构成的向量的叉积除以2依次累加, 所得之和就是该封闭多边形的面积

具体过程如下图:

因为叉积是两向量所夹的平行四边形的有向面积,而对于一个封闭图形顶点A,一定会存在两个顶点B、C使得\(\vec{OA}\times\vec{OC}\)\(\vec{OB}\times\vec{OA}\)反向,每块多边形外的面积一定会被抵消

因此,对于一个多边形A,其顶点一次为{\(A_1,A_2...A_n\)},其面积公式为:

\[2S=\vec{OA_n}\times\vec{OA_1}+\sum_{i=1}^{n-1}\vec{OA_i}\times\vec{OA_{i+1}} \]

凸包

所谓凸包,就是在最小化周长的情况下,一个能将所有给定的点都覆盖的图形。

如下图:左图不是凸包,而改成右图后就成了凸包

凸包由逆时针方向遍历,每条边与下一条边的叉积一定为
(由此可以求出凸包)

求凸包

这里分享一下一次性求出整个凸包的方法:

  1. 以(y,x),找到最小的点
  2. 以该点建立极坐标系(即将该点当做(0,0))
  3. 将剩下的点按 (极角,极径)排序
  4. 利用相邻边叉积为正的性质,单调栈维护

(脑补一下画面,很好理解)

bool cmp(vec a,vec b){return a*b==0 ? dis(a) < dis(b) : a*b>0; }
void convexhull()
{
    for(int i=n; i>=1; --i) a[i]=a[i]-a[1];
    sort(a+2,a+n+1,cmp);
    for(int i=1;i<=n;++i) {
        st[++top]=a[i];
        while(top>2 && (st[top-1]-st[top-2])*(st[top]-st[top-1])<0) st[top-1]=st[top], top--;
    }
}

判断点是否在凸包内

将点A,加入上面的方法建立的极坐标系,lower_bound得到极角排序后在它之前与在它之后的,然后跨立实验即可

旋转卡壳

如何求平面上的最远点对?

思路:

不难发现,平面上的最远点对一定存在于凸包上(画图很好理解)。
凸包上,距离点A最远的点B,那么离A最远的,一定以B为端点。
换句话说,若能确定凸包上某条边与某个点相距最远,则最远的点对一定是该点与边的某个端点
感性理解一下: 若以该点向边作垂线,连端点,端点于该点连线一定是斜边,斜边大于直角边

如何找到离每条边最远的点(对踵点)呢?
用选定的边做平行直线,恰好接触到凸包时,平行线与凸包的交点就是了
而凸包上边的斜率变化是单调的, 因此不难发现对踵点的位置也是单调的
即:随着逆时针选取边,其对踵点也会逆时针旋转

Code

db dis_line(vec a,vec b,vec c) { return (b-a)*(c-a)/(b-a).dis(); }//点C到直线AB的距离
ll qiaqiao()
{
    if(top==1) return 0;
    if(top==2) return (st[1]-st[2]).dis();
    st[1+top]=st[1]; //便于找边
    int j=1; db ans=0;//j记录对踵点
    for(re int i=1;i<=top;++i) {
        while(dis_line(st[i],st[i+1],st[j])<dis_line(st[i],st[i+1],st[j+1])) j=j%top+1;
        ans=max(ans, max((st[i]-st[j]).dis(), (st[i+1]-st[j]).dis()));
    }
    return ans;
}

半平面交

常见应用

  1. 多边形面积交
  2. 多边形的核
    (能在该区域“看到”多边形内部的每个角落)
  3. 求解二元一次不等式组

算法本质

求解由二元一次不等式组的可行解
例题:射箭
如多条直线\(l_i: A_ix+B_iy+C_i\geq0\) 所包含的面积

因为,\(A_ix+B_iy+C_i\geq0\) 类似的不等式,其中合法的(x,y)只能是这条直线的一个半平面.

对于不等式到半平面交中的转化:
(以求直线左侧半平面交为例)

  1. \(l_i: A_ix+B_iy+C_i\geq0\), 则把直线的方向与x轴相同,若是垂直与x轴,方向竖直向上
  2. \(l_i: A_ix+B_iy+C_i\leq0\), 则把直线的方向与x轴相反,若是垂直与x轴,方向竖直向下

算法思路

以逆时针方向为正,则所有边的左边部分都是可能保留的半平面,多个图形的所有边的半平面的交集即为多边形面积交,单个图形的半平面交为多边形的核
具体求法如下:

  1. 将所有线段按方向向量的极角排序
  2. 将所有线段一次插入双端队列:
    1. 比较新加入的边是否遮盖队头,是则弹出队头
    2. 比较新加入的边是否遮盖队尾,是则弹出队尾
  3. 最后检查队列的头尾是否能遮盖

  1. 极角排序过程中,若有两线段平行(极角相同),保留更靠左的那条(靠左的遮盖了靠右的)
  2. 上述中的遮盖,指的是:若a能遮盖b,则双端队列中有了a,即使去除b后半平面交也不会有变化

关于遮盖与否的判断方式:

待加入的线段a与倒数第二个线段c的交点最后两个线段b,c的交点的左边,则代表线段a在半平面交中能完全替代线段b

Code

struct line{ 
    vec st, ed, c;//c是st->ed的向量
    bool operator <(const line l)const {
        db A=atan2(c.y,c.x), B=atan2(l.c.y,l.c.x);
        if(fabs(A-B)>eps) return A<B;
        return c*(l.ed-st)<0;
    }
}l[_*_], qwe[_*_];
vec intersection(line l1,line l2) //两直线求交
{
    vec st1=l1.st, st2=l2.st;
    vec a=l1.c, b=l2.c, c=st2-st1;
    if(fabs(b*a)<eps) return vec(-1e12,-1e12);
    return st1+(a^((b*c)/(b*a)));
}
int n, m, hd, tl;
line q[_*_];
bool check(line a,line b,line c) 
{//通过叉积来比较b,c交点与a向量的左右关系,进而推出线段是否遮盖
    vec p=intersection(b,c);
    return (a.c*(p-a.st))+eps<=0;
}
in void solve(int cnt)
{
    int tot=0;
    sort(l+1,l+cnt+1); qwe[++tot]=l[1];
    for(re int i=2;i<=cnt;++i)
    {
        if(fabs(atan2(l[i].c.y,l[i].c.x)-atan2(l[i-1].c.y, l[i-1].c.x))<eps) continue;
        qwe[++tot]=l[i];
    } //去重
    for(re int i=1;i<=tot;++i) l[i]=qwe[i];
    m=tot;
    for(re int i=1;i<=m;++i) //双端队列维护
    {
        while(tl>hd && check(l[i], q[tl], q[tl-1])) tl--;
        while(tl>hd && check(l[i], q[hd], q[hd+1])) hd++;
        q[++tl]=l[i];
    }
    while(tl>hd && check(q[tl], q[hd], q[hd+1])) hd++;
    while(tl>hd && check(q[hd], q[tl], q[tl-1])) tl--;
}

最小圆覆盖

算法流程

  1. random_shuffle
  2. 枚举第一个点i作为圆心
  3. \([1,i)\) 中枚举j, 并以 \(ij\) 作为直径
  4. \([1,j)\) 中枚举k,并找到 \(\triangle_{ijk}\) 的外接圆
    找外接圆,可以通过三角形两边的中垂线的交点

在上述234过程中不断维护最小圆的圆心和半径,若一个点已在当前圆内,则跳过它

Code(LG模板)

#include<bits/stdc++.h>
using namespace std;
#define re register
#define ll long long
#define get getchar()
#define in inline
#define db long double
const int _=2e5+23;
const db eps = 1e-15;
int n;
struct vec{
    db x,y;
    vec(){x=y=0;}
    vec(db a,db b){x=a, y=b;}
    vec operator+(const vec b)const{return vec(x+b.x,y+b.y);}
    vec operator-(const vec b)const{return vec(x-b.x,y-b.y);}
    vec operator/(const db b)const {return vec(x/b, y/b);}
    db operator*(const vec b)const{return x*b.y-y*b.x;}
    vec operator^(const db b)const{return vec(x*b, y*b);}
    db dis(){ return sqrt(x*x+y*y);}
    vec rot(){ return vec(-y,x);}
}a[_];
vec intersection(vec p,vec a,vec q,vec b) //p是点,a是从p出发的向量,q、b同理
{
    vec c=(q-p);
    return p+(a^((b*c)/(b*a)));
}
vec getcircle(vec A,vec B,vec C)//找到三角形abc的外接圆圆心
{ 
    return intersection((A+B)/2,(B-A).rot(),(C+A)/2,(C-A).rot());
}
int main()
{
    scanf("%d",&n);
    for(re int i=1;i<=n;++i) scanf("%Lf%Lf", &a[i].x, &a[i].y);
    if(n==1) { cout<<0<<endl<<a[1].x<<' '<<a[1].y<<endl; return 0; }
    vec o; db R=0;
    for(re int i=2;i<=n;++i)
    {
        if((a[i]-o).dis()<=R) continue;
        o=a[i];
        for(re int j=1;j<i;++j)
        {
            if((a[j]-o).dis()<=R) continue;
            o=(a[j]+a[i])/2; R=(o-a[i]).dis();
            for(re int k=1;k<j;++k)
            {
                if((a[k]-o).dis()<=R) continue;
                o=getcircle(a[i],a[j],a[k]); R=(o-a[i]).dis();
            }
        }
    }
    printf("%.10Lf\n%.10Lf %.10Lf\n",R,o.x,o.y);
    return 0;
}

附赠上述所有模板(除半平面交,圆覆盖)
#include<bits/stdc++.h>
using namespace std;
#define re register
#define ll long long
#define get getchar()
#define in inline
in int read()
{
    int t=0, x=1; char ch=get;
    while(ch!='-' && (ch<'0' || ch>'9')) ch=get;
    if(ch=='-') ch=get, x=-1;
    while(ch<='9' && ch>='0') t=t*10+ch-'0', ch=get;
    return t*x;
}
#define db double
const int _=1e5+23;
const db inf=1e9, pi = acos(-1), eps=1e-10;

struct vec{//存储向量或者点坐标
    db x,y;
    vec(){x=y=0;}
    vec(db a,db b){x=a,y=b;}
    vec operator+(const vec b) const{ return (vec){x+b.x, y+b.y}; }
    vec operator-(const vec b) const{ return (vec){x-b.x, y-b.y}; }
    db operator*(const vec b) const{ return x*b.y-y*b.x; } //叉积
    db operator&(const vec b) const{ return x*b.x+y*b.y; } //点积
    vec operator^(const db b) const{ return (vec){x*b,y*b}; } //数乘
    db dis(){ return sqrt(x*x+y*y);} //向量模长/距离
};

bool cmp(vec a,vec b){return a*b==0 ? a.dis()<b.dis() : a*b>0;}//极角排序
vec st[_]; int top;
void convexhull(vec *a,int n)//求凸包
{
    for(int i=n;i>=1;--i) a[i]=a[i]-a[1];
    sort(a+2,a+n+1,cmp); top=0;
    for(int i=1;i<=n;++i) {
        st[++top]=a[i];
        while(top>2 && (st[top-1]-st[top-2])*(st[top]-st[top-1])<0) st[top-1]=st[top], top--;
    }
}

in db dis_line(vec a,vec b,vec c) { return (b-a)*(c-a)/(b-a).dis(); }//点到直线的距离

in ll qiaqiao() //旋转卡壳(凸包直径)
{
    if(top==1) return 0;
    if(top==2) return (st[1]-st[2]).dis();
    st[1+top]=st[1]; int j=1; db ans=0;
    for(re int i=1;i<=top;++i)
    {
        while(dis_line(st[i],st[i+1],st[j])<dis_line(st[i],st[i+1],st[j+1])) j=j%top+1;
        ans=max(ans, max((st[i]-st[j]).dis(), (st[i+1]-st[j]).dis()));
    }
    return ans;
}
bool iscross(vec a,vec b,vec c,vec d) //ab,cd两线段是否相交
{
    if(((a-c)*(d-c))*((d-c)*(b-c))<0) return 0;
    if(((c-a)*(b-a))*((b-a)*(d-a))<0) return 0;
    return 1;
}

vec crossover(vec A1,vec A2,vec B1,vec B2) //A1A2,B1B2两直线交点
{
    vec a=(A2-A1), b=(B2-B1), c=(B1-A1);
    if(fabs(a*b)<eps) return (vec){-inf,-inf}; //两直线平行
    return A1+(a^((b*c)/(b*a)));
}

db getarea(vec *A,int n) //求以A点集为顶点的多边形面积
{
    db s=(A[n]*A[1])/2.000;
    for(re int i=1;i<n;++i) s+=(A[i]*A[i+1])/2.000;
    return fabs(s);
}

int main()
{
    return 0;
}

posted @ 2021-12-17 20:18  yzhx  阅读(272)  评论(4编辑  收藏  举报