【笔记篇】最良心的计算几何学习笔记（三）

广告:

1. 先是放一下本文的::github传送门:: (不知道为什么要放)
2. 今天发现了一个AMA(Ask me anything)的东西, 觉得非常好玩, 就fork了一个放到自己的github里面,
   估计没有人会来问, 所以就放到这里拉拢人气(虽然这里也拉拢不到) 欢迎大家来玩哦~
   地址请戳这个就是传送门啦~

今天继续计算几何(明明已经颓废了半下午了

计算多边形面积

我们先从最简单的多边形——三角形开始看.

如何计算\(\triangle ABC\)的面积? 这个问题数学课上老师应该说过..

• \(S=\frac{1}{2}ah\)
  这个是最最最常见的公式, 但是这里我们并不知道高, 要算起来就比较麻烦.

• \(S=\frac{1}{2}bcsinA\)(其他两角同理)
  这个看上去比较靠谱. 我们算一下\(\vec{AB}\times\vec{AC}\)的绝对值就好了...

• \(S=\sqrt{p(p-a)(p-b)(p-c)},p=\frac{a+b+c}{2}\)
  海伦-秦九韶公式啊OvO 这个是可以算的... 但是如果用在多边形就不是很好用了.

简单的三角形我们看完了, 我们来看看多边形..
有些多边形我们都已经熟知面积公式(比如长方形啊 平行四边形啊 梯形啊什么的)
就不再提了.
来看看凸多边形...
还记得上次可爱的凸多边形么_ (:з」∠)_?

我们要计算它的面积的时候, 只需要像图中一样划分成若干个三角形, 然后用公式
\[ S=\sum_{i=1}^{7}bcsinA=\sum_{i=1}^{7}|\frac{1}{2}\vec{E_i}\times\vec{E_{i+1}}| \]
计算总面积即可.

但是对于凹多边形呢? 我们还是先划一下三角形...

我们会发现如果再按照上面的方法计算的话黄色和紫色(似乎故意标淡了点)的面积会重复计算, 显然是大于多边形面积的. 但是数学老师教的面积公式毕竟还是和我们的叉积不一样的, 公式算的是绝对值, 而叉积是有正负的.
如果\(E_i\)在\(E_{i+1}\)的顺时针方向, 面积算出来的是正值; 否则算出来的是负值.
发现刚才重复的黄色和紫色部分如果用叉乘算一下刚好是一正一负, 多余的面积都不见了..
再再求一波总和就做完了, 轻松加愉快...
而且有了这种正负的定义之后, 我们又有了一种新操作:

以某个点为出发点向多边形做向量, 一路做叉积绕个圈求出来的和也等于面积~
为了方便起见, 完全可以让"某个点"取原点\(O\), 这样从数值上来说我们只需要把点的坐标做叉积即可了~

且慢! 不是还有一种自我重叠的多边形吗? 这个方法也适用吗?
这个我就不配图了(其实是嫌麻烦←_←) 完全可以自己画一下..
发现是完全适用的, 而且自我重叠的部分的面积会计算正确的次数哦~

然后就是最后的总面积有可能是个负值, 可以视情况取个绝对值什么的^_^

贴代码(仍然并没有找到板子题~)(似乎是因为代码太简单了?

//求任意多边形面积
double polyArea(point *pts,int pcnt,double s=0){
    pts[pcnt]=pts[0];
    for(int i=0;i<pcnt;++i)
        s=pts[i]*pts[i+1]+s;
    return 0.5*s;
}

这样就搞定了OvO

计算多边形重心

这个也分很多情况啊OvO
而且这个涉及到了高端的数学及物理知识(头疼ing...

质量集中在顶点上

那就是每个顶点的质量关于坐标的平均咯~

\[ X=\frac{\sum_{i=0}^{n-1}x_im_i}{\sum_{i=0}^{n-1}m_i}\Y=\frac{\sum_{i=0}^{n-1}y_im_i}{\sum_{i=0}^{n-1}m_i} \]

质量均匀分布

还是从简单开始, 三角形的重心.
懒得再推了, 数学老师说坐标应该是\((\frac{x_1+x_2+x_3}3,\frac{y_1+y_2+y_3}3)\)..
所以我们就同样可以把多边形三角剖分, 每个三角形的质量都等效到中心去.
然后就变成了质量集中在顶点上的情况, 质量就取三角形的面积(注意是有向面积)即可.
要注意的比如总面积是0的时候, 因为要做分母, 所以要特殊处理.
板子题hdu1115适合写一下.
不过又被-0.00卡翻.. 做了一波优化把\(\frac1 3\)提出来竟然就A掉了, 不是很懂OvO
代码:

#include <cmath> 
#include <cstdio>
const double eps=1e-9;
int dcmp(const double &a){
    if(fabs(a)<eps) return 0;
    return a<0?-1:1;
}
struct point{
    double x,y;
    point(double X=0,double Y=0):x(X),y(Y){}
}poly[1000010],s;
double operator *(const point &a,const point &b){
    return a.x*b.y-a.y*b.x;
}
point polyCenter(point *pts,int pcnt,double sx=0,double sy=0,double area=0){
    pts[pcnt]=pts[0]; double ar;
    for(int i=0;i<pcnt;++i){
        ar=pts[i]*pts[i+1];
        sx+=(pts[i].x+pts[i+1].x)*ar; //这里如果写sx+=(pts[i].x+pts[i+1].x)/3*ar;
        sy+=(pts[i].y+pts[i+1].y)*ar; //这个地方写sy+=(pts[i].y+pts[i+1].y)/3*ar;
        area+=ar;
    } area*=3; //而这个地方不写的话就会被卡精度:-(
    return point(sx/area,sy/area);
}
int main(){
    int T; scanf("%d",&T);
    while(T--){
        int n;scanf("%d",&n);
        for(int i=0;i<n;++i)
            scanf("%lf%lf",&poly[i].x,&poly[i].y);
        s=polyCenter(poly,n);
        printf("%.2lf %.2lf\n",s.x,s.y);
    }
}

质量不均匀分布

这个据说要用到积分?
反正我是不会的←_←
等见到再考虑学不学吧..
估计(希望)我是见不到了(flag

然后还有一些点或许因为太麻烦, 或许因为不常见还没有学到..
比如什么求多边形之内最大的圆之类的.
据说特别麻烦, 等到有空或者用得到的时候再研究吧.
下次就该学学"更计算几何"的一些知识了
比如凸包.

再随便多说几句:
遇到多边形的问题要先考虑(读题)看分不分凹凸, 是不是简单.
一般让多边形第n个点等于第0个点做起来会很舒服.
计算几何都是毒瘤题见到还是直接弃疗吧←_←

但是这篇文章的长度似乎不太够...
我们再加一丢丢内容吧...

平面最近点对

暴力枚举每一对显然就是\(O(n^2)\)的, 那是很显然过不了的.
似乎有一些玄学的做法比如随机转个角度防卡然后分块, 但是这种做法看着就不科学...
我们要思考科学的方法, 比如考虑分治解决问题.

先将所有点按横坐标排个序.
最近点对的这两个点的分布只可能有三种情况:
都在左边、都在右边、左右各一.
对于前两种情况递归下去即可.
我们主要来处理左右各一的情况.
我们假设左右两边递归后求出的值的较小者为\(d\).(图中的r)
那很显然我们只需要考虑[mid-d,mid]和[mid+d,mid]中的点.
如果还是暴力 比较坏的情况复杂度跟暴力并没有什么区别, 还是\(O(n^2)\)的.
但是因为要求的是最近点对, 所以我们可以限制一波.

对于左侧的P点来说, 假如说\(d\)是左右两边求的最小值, 显然我们要找的点和\(p\)之间的横坐标是要小于\(d\)的
而这个矩形中最多放多少个互相距离不超过\(d\)的点呢? 答案是6个.
为什么呢? 我们将宽平均分成2份, 高平均分成3份,(红色) 这样就形成了一个2*3的格子.
每个格子的宽就是\(\frac1 2d\), 高就是\(\frac2 3d\), 由勾股定理, 对角线(蓝色)长为\(\sqrt{(\frac1 2d)^2+(\frac2 3d)^2}=\frac5 6d<d\)
也就是说不可能存在一个格子中能存在两个距离大于\(d\)的点.
那么根据抽屉原理, 最多就只有6个点了.
所以我们只需要找这些点进行检索即可, 这样就保证复杂度不会太高了.

还有一点小细节就是我们按\(y\)找的时候可以采用归并的方式, 每次只排子序列, 这样可以把复杂度控制在\(O(nlogn)\)级别, 这样这个问题就得到了完美的解决.
代码:

#include <cmath>
#include <cstdio>
#include <algorithm>
using namespace std;
const double eps=1e-9;
int dcmp(const double &a){
    if(fabs(a)<eps) return 0;
    return a<0?-1:1; 
}
struct point{
    double x,y;
    point(double X=0,double Y=0){}
}p[200020];
int t[200020];
inline bool cmpx(const point &a,const point &b){
    if(a.x==b.x) return a.y<b.y;
    return a.x<b.x;
}
inline bool cmpy(const int &a,const int &b){
    return p[a].y<p[b].y;
}
inline double dist(const point &a,const point &b){
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
double solve(int l,int r){
    if(r==l) return 1e9;
    if(r-l==1) return dist(p[l],p[r]);
    int mid=(l+r)>>1;
    double dl=solve(l,mid);
    double dr=solve(mid+1,r);
    if(dr<dl) dl=dr;
    
    int tot=0; double dis=0;
    for(int i=l;i<=r;++i)
        if(dcmp(fabs(p[i].x-p[mid].x)-dl)<0)
            t[tot++]=i; //合法的点才加入数组
    sort(t,t+tot,cmpy);
    for(int i=0;i<tot;++i)
        for(int j=i+1;j<tot&&p[t[j]].y-p[t[i]].y<dl;++j){
            if((dis=dist(p[t[i]],p[t[j]]))<dl) dl=dis;
        } //左右两边都在搜所以只需要考虑下半个矩形
    return dl;
}
inline int gn(int a=0,char c=0){
    for(;c<'0'||c>'9';c=getchar());
    for(;c>47&&c<58;c=getchar())a=a*10+c-48;return a;
}
int main(){
    int n=gn();
    for(int i=1;i<=n;++i) p[i].x=gn(),p[i].y=gn();
    sort(p+1,p+n+1,cmpx);
    printf("%.4lf",solve(1,n));
} 

那么就这样咯~
这一篇我竟然拖了两天..