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

广告:

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

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

计算多边形面积

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

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

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

• (S=frac{1}{2}bcsinA)(其他两角同理)
  这个看上去比较靠谱. 我们算一下(vec{AB} imesvec{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} imesvec{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
",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));
} 

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