计算几何 val.2
2024-08-29 03:44:43
计算几何 val.2
前置芝士：基础操作以及凸包
本文主要写旋转卡壳、半平面交、最小圆覆盖要注意的内容
几何单位结构体板子
不全（我知道
struct point{
double x,y;
point(double x=0,double y=0): x(x),y(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-(point b){
return point(x-b.x,y-b.y);
}
point operator+(point b){
return point(x+b.x,y+b.y);
}//向量运算
point operator*(double b){ //数乘
return point(x*b,y*b);
}
db dis(){ //模长
return sqrt(x*x+y*y);
}
}b[N];
int comp0(double x){
return fabs(x)<=eps?0:(x>0?1:-1);
}//判0，防止精度误差
struct line{
point p,v;//p起点，v终点 ，表示线段
double theta;
bool operator <(line y){
return comp0(theta-y.theta)==0?comp0((y.v-p)*(v-p))<0:comp0(theta-y.theta)<0;
}//排序，保证第一关键字是极角，第二关键字是与右边的距离（用叉积判相对关系，在左边的放后面）
};
point inter(line a,line b){//交点 intersection
point v1=a.v-a.p,v2=b.v-b.p;
return b.p+v2*(((b.p-a.p)*v1)/(v1*v2));
}//此处默认了没有平行的情况 , 判断线段有没有交点就是判断点是否在线段上，但是半平面交是直线（只是用线段表示）
旋转卡壳
\[\Large{\text{旋↗转↘↗卡↘↗壳(ko↗)}}
\]
\]
对，非常正确
此算法用来求凸多边形直径
基础概念
切线
过顶点的一条线，这个多边形都在这条线的一侧
对踵点
多边形上两个点作一对平行线，整个多边形都在这两条线之间
求法
考虑逆时针枚举每条边并找到距离这条边最远的点，那么这个点和这条边的一个顶点构成的直径是可能的答案
我们惊奇地发现，点出现的顺序也是逆时针的，那么就可以\(O(n)\)求出了
（当然求凸包还要\(n\log n\)
补充：如果要求单独一条边的最远点的话，可以三分求（单峰函数）
模板
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define db double
using namespace std;
struct point{
double x,y;
point(double x=0,double y=0): x(x),y(y){}
double operator*(point b){
return x*b.y-y*b.x;
}
point operator-(point b){
return point(x-b.x,y-b.y);
}
point operator+(point b){
return point(x+b.x,y+b.y);
}
db dis(){
return sqrt(x*x+y*y);
}
};
const int N = 50021;
point p[N],h[N];
int tp=0,stk[N],usd[N];
int cmp(point a,point b){
return a.x==b.x?a.y<b.y:a.x<b.x;
}
int n=0;
db Fabs(db a){
return a>0?a:-a;
}
db disl(point a,point b,point x){
return Fabs((a-x)*(b-x)/(a-b).dis());
}
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++){
scanf("%lf%lf",&p[i].x,&p[i].y);
}
sort(p+1,p+n+1,cmp);
stk[++tp]=1;
for(int i=2;i<=n;i++){
while(tp>1&&(p[stk[tp]]-p[stk[tp-1]])*(p[i]-p[stk[tp]])<=0) usd[stk[tp--]]=0;
usd[i]=1;stk[++tp]=i;
}
int ntp=tp;
for(int i=n-1;i>=1;i--){
if(!usd[i]){
while(tp>ntp&&(p[stk[tp]]-p[stk[tp-1]])*(p[i]-p[stk[tp]])<=0) usd[stk[tp--]]=0;
usd[i]=1;stk[++tp]=i;
}
}
for(int i=1;i<=tp;i++){
h[i]=p[stk[i]];
}
if(tp<=2){ //只有一个点
puts("0");
return 0;
}
if(tp==3){
printf("%.0lf\n",(h[1]-h[2]).dis()*(h[1]-h[2]).dis());
return 0;
}
db ans=0;
int t=1;
//注意最后一个点（就是1号点）要保留，因为后面有h[i+1],[t+1]，方便操作
for(int i=1;i<tp;i++){
while(disl(h[i],h[i+1],h[t])<disl(h[i],h[i+1],h[t+1])) t=t%(tp-1)+1; //逆时针枚举点
ans=max(ans,max((h[i]-h[t]).dis(),(h[i+1]-h[t]).dis())); //两端点到此点
}
printf("%.0f",ans*ans);
return 0;
}
半平面交
半平面此处我们用向量表示，以左侧或右侧为半平面
前置芝士：线段交
画图吧，用相似得出
\[h_1=\frac{fabs((v_2-p_1)*(p_2-p_1))}{dis(v_2-p_2)},h_2=\frac{fabs((v_2-v_1)*(p_2-v_1))}{dis(v_2-p_2)}
\]
\]
\[line(p_1,v_1) \cap line(p_2,v_2)=p_1+(v_1-p_1)*\frac{h_1}{h_1+h_2}
\]
\]
S&I算法
极角排序
令\(\theta=arctan\frac y x\)，以其从小到大排序，得到新的向量集
算法流程
以逆时针为正方向，建边(线段)，注意要搞一个巨大的框然后来切这个框
对线段根据极角排序
去除极角相同的情况下，位置在右边的边（如果求的是左半平面交（逆时针凸多边形）的话）
用双端队列储存线段集合 \(L\)，遍历所有线段
判断该线段加入后对半平面交的影响，(对双端队列的头部和尾部进行判断，因为线段加入是有序的，有影响是指交点在这条线的另一边（不需要的一边）)，即判断之前的交点是否在次线段右侧
上一步一定要先删尾部再删头部，具体见wjyyy大佬的博客，后文会附上
最后判断形成环的影响（尾部再更新头部一次，头部再更新尾部一次）
最后剩下的线段集合 \(L\)，即使最后要求的半平面交
全程都可以用叉积判方向
模板
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define db double
const double eps=1e-9;
using namespace std;
const int N = 10001;
struct point{
double x,y;
point(double x=0,double y=0): x(x),y(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-(point b){
return point(x-b.x,y-b.y);
}
point operator+(point b){
return point(x+b.x,y+b.y);
}
point operator*(double b){ //数乘
return point(x*b,y*b);
}
db dis(){ //模长
return sqrt(x*x+y*y);
}
}b[N];
int comp0(double x){
return fabs(x)<=eps?0:(x>0?1:-1);
}
struct line{
point p,v;//p起点，v终点 ，表示线段
double theta;
bool operator <(line y){
return comp0(theta-y.theta)==0?comp0((y.v-p)*(v-p))<0:comp0(theta-y.theta)<0;
}//排序，保证第一关键字是极角，第二关键字是与右边的距离（用叉积判相对关系，在左边的放后面）
};
point inter(line a,line b){//交点 intersection
point v1=a.v-a.p,v2=b.v-b.p;
return b.p+v2*(((b.p-a.p)*v1)/(v1*v2));
}//此处默认了没有平行的情况 , 判断线段有没有交点就是判断点是否在线段上，但是半平面交是直线（只是用线段表示）
int out(line a,line b,line k){
point ins=inter(a,b); //此处要判断的是ins是否在k的右边,画个图
return comp0((k.v-k.p)*(ins-k.p))<0;
}
int n;
line a[N],q[N];int top;
void work(){
sort(a+1,a+top+1);
int cnt=0;
for(int i=1;i<=top;i++){
if(comp0(a[i].theta-a[i-1].theta)!=0) cnt++;
a[cnt]=a[i];//极角相同时，右边的更优（排序保证了在左边）
}
int h=1,t=0;
q[++t]=a[1],q[++t]=a[2];//既然是交，至少包含两个元素
for(int i=3;i<=cnt;i++){
while(h<t&&out(q[t-1],q[t],a[i])) t--; //踢掉一些点 ，注意一定要先踢后面，不然会有些错
while(h<t&&out(q[h+1],q[h],a[i])) h++;
q[++t]=a[i];
}
while(h<t&&out(q[t-1],q[t],q[h])) t--; //由于是环形,判断影响
while(h<t&&out(q[h+1],q[h],q[t])) h++;
q[t+1]=q[h];
top=0;
for(int i=h;i<=t;i++){
b[++top]=inter(q[i],q[i+1]);
}
}
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++){
int k;scanf("%d",&k);
for(int j=1;j<=k;j++){
scanf("%lf%lf",&b[j].x,&b[j].y);
}
b[k+1]=b[1];
for(int j=1;j<=k;j++){
a[++top].p=b[j];a[top].v=b[j+1];//逆时针给出，如果不知道的话可以叉积判断
}
}
for(int i=1;i<=top;i++){
a[i].theta=atan2(a[i].v.y-a[i].p.y,a[i].v.x-a[i].p.x);
}
work();
double ans=0;
if(top<=2){
printf("%.3lf",0.0);return 0;//二 边 形
}
b[top+1]=b[1]; //同样，环形处理方法
for(int i=1;i<=top;i++){
ans+=(b[i]*b[i+1]);
}
ans/=2.0;
if(comp0(ans)==0) printf("%.3lf",0.0);
else{
printf("%.3lf",fabs(ans));
}
return 0;
}
最小圆覆盖
使用随机增量法
随机增量法
考虑钦定\(i,j\)一定在圆上，\(j < i\)
对于\(\forall k<j\)，考虑构造覆盖\(i,j,k\)的圆
如果当前\(k\)在圆内，计算\(k+1\)
否则更新为\(i,j,k\) 的外接圆（三点确定圆）
对于前面两层，如果某个点在圆内，直接跳过
否则枚举\(j\in [1,i)\) ，重新计算
时间复杂度
看起来是三层循环啊QwQ，在圆内的点再多也优化不了多少吧？
但是我们最开始 \(\text{random_suffle}\) 了一下，于是复杂度变成了期望意义下的
那么是多少呢？男默女泪的是，它达到了惊人的\(O(n)\)！
分析：
- 对于过\(P_i,P_j\)的圆，每个点至少循环了一遍，为\(O(j)\)
- 对于过\(P_i\)的圆，考虑覆盖前\(i\)个点的圆是由3个点确定的，于是在前\(i-1\)个点中，期望有2个点计算了下一层
- 对于所有点的最小圆覆盖，期望有三个点能做出贡献
- 总复杂度：\(O\left(\sum_{i=1}^{n} \frac{3}{i} \sum_{j=1}^{i} \frac{2}{j} \cdot j\right)=O(6\cdot n)=O(n)\)
模板
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define db double
#include<cstdlib>
#include<ctime>
const double eps=1e-17;
using namespace std;
const int N = 100001;
struct point{
double x,y;
point(double x=0,double y=0): x(x),y(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-(point b){
return point(x-b.x,y-b.y);
}
point operator+(point b){
return point(x+b.x,y+b.y);
}
point operator*(double b){ //数乘
return point(x*b,y*b);
}
db dis(){ //模长
return sqrt(x*x+y*y);
}
}p[N];
int comp0(double x){
return fabs(x)<=eps?0:(x>0?1:-1);
}
struct line{
point p,v;//p起点，v终点 ，表示线段
line(point a,point b): p(a),v(b){}
double theta;
bool operator <(line y){
return comp0(theta-y.theta)==0?comp0((y.v-p)*(v-p))<0:comp0(theta-y.theta)<0;
}//排序，保证第一关键字是极角，第二关键字是与右边的距离（用叉积判相对关系，在左边的放后面）
};
point inter(line a,line b){//交点 intersection
point v1=a.v-a.p,v2=b.v-b.p;
return b.p+v2*(((b.p-a.p)*v1)/(v1*v2));
}//此处默认了没有平行的情况 , 判断线段有没有交点就是判断点是否在线段上，但是半平面交是直线（只是用线段表示）
int out(line a,line b,line k){
point ins=inter(a,b); //此处要判断的是ins是否在k的右边,画个图
return comp0((k.v-k.p)*(ins-k.p))<0;
}
int n;
int in(point k,point c,double r){
return comp0(r-(k-c).dis())>=0;
}
point chrt(point a){
return point(-a.y,a.x);
}
void randomShuffle(){
srand(19260817+time(0)); //知道为什么不用STL的吗？
for (int i=1;i<=n;i++){
int j=(rand()%n+1926)%n+1;//这里随便rand()%n就行了，但我们要把玄学发扬光大【滑稽】
swap(p[i],p[j]);
}
}
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++){
scanf("%lf%lf",&p[i].x,&p[i].y);
}
randomShuffle(); //这句很重要
point c;
db r=0.0;
for(int i=1;i<=n;i++){
if(in(p[i],c,r)) continue;
c=p[i],r=0;//重新计算
for(int j=1;j<i;++j){
if(in(p[j],c,r)) continue;
r=(p[j]-p[i]).dis()/2.0;
c=p[i]*0.5+p[j]*0.5;
for(int k=1;k<j;k++){
if(in(p[k],c,r)) continue;
line a=line((p[i]+p[j])*0.5,chrt(p[i]-p[j])+(p[i]+p[j])*0.5);
line b=line((p[k]+p[j])*0.5,chrt(p[j]-p[k])+(p[j]+p[k])*0.5);
c=inter(a,b);r=(c-p[i]).dis();
}
}
}
printf("%.10f\n%.10f %.10f",r,c.x,c.y);
return 0;
}
后记
应该还有val.3 （学习nekopara
就写一写 定积分基础 / 闵可夫斯基和 / 辛普森积分 (自适应辛普森) 啥的
最新文章
- Angular2 组件通信
- JAVA程序中SQL语句无法传递中文参数
- 常见.NET功能代码汇总 （2）
- 总结-Intellij Idea (快捷键 配置修改)
- 【转】Hibernate 常见异常
- 网页中的CSS换行控制
- C语言学习笔记（一）_hello world
- ORACLE创建表空间、创建用户、更改用户默认表空间以及授权、查看权限
- VB动态添加WebBrowser控件，并拦截弹出窗口(不用引用任何组件)
- hbase日常运维管用命令，region管理
- POJ3278 Catch That Cow(BFS)
- javascript的本地存储 cookies、localStorage
- Asp.net MVC razor语法参考
- 在windows平台下忘记了root的密码如何解决?
- Java设计模式菜鸟系列(九)外观模式建模与实现
- 一个栗子上手CSS3动画
- 团队作业8——第二次项目冲刺（Beta阶段）--第三天
- Linux系列教程（二）——Linux系统安装（手把手学安装centos6.8）
- QT Designer基础——登录界面设计基础版2
- dpkg制作deb包