计算几何进阶
The Skyline Problem, UVa105
description
给出\(n\) 个建筑物的位置，求出其轮廓线。每个建筑物均为长方形且共底。
solution
对于每个位置，不难发现它的轮廓线纵坐标就是包含它的所有建筑物中最高的纵坐标。因此注意下边界情况，直接模拟即可。
code
#include<bits/stdc++.h>
using namespace std;
const int N=10005;
int l,r,t,h[N],mxr;
int main()
{
while(scanf("%d%d%d",&l,&t,&r)==3)
{
for(int i=l;i<r;++i)h[i]=max(h[i],t);
mxr=max(mxr,r);
}
for(int i=1;i<mxr;++i)
if(h[i]!=h[i-1])printf("%d %d ",i,h[i]);
printf("%d %d\n",mxr,0);
return 0;
}
HDU 1542 Atlantis
description
求\(n\) 个矩形的面积并。（所有矩形的边都和\(x\) 轴垂直或平行）
solution
原题数据范围较小，可以直接暴力做，这里我们还是考虑时间复杂度较低的做法。
考虑有一个平行于\(y\) 轴的直线从最左侧扫到最右侧。对于每个矩形的两个宽，称横坐标较小的为入边，横坐标较大的为出边。扫描线扫到入边时将对应区间覆盖一次，扫到出边时将对应区间清除一次。我们称相邻两次修改（即某两个区间的覆盖或清除）间的时段为一个过程。显然的是对于每个过程，覆盖的总长度是不变的，而此过程也有自己的长度，二者的乘积就是此过程所对应的面积。我们只需不断累加，就可以求出面积并了。
具体实现时可以采用线段树。考虑将所有矩形的两个长的横坐标进行离散化，而后将线段树上的叶子节点对应于一个区间（先前离散化后相邻的两个横坐标形成一个区间）。记下此区间的长度。现在只用维护全局的覆盖长度。对于修改（覆盖或清除一次的操作），也可以通过在线段树上打标记来实现（注意此标记并不需要下传，类似标记永久化）。
time complexity
\(\mathcal O(n\log_2n)\)
code
#include<cstdio>
#include<algorithm>
using namespace std;
typedef double db;
const int N=100005;
int n,cnt,tot;
db cpy[N<<1],len[N<<1];
struct Seg{db l,r,p;int v;}e[N<<1];
inline bool operator<(const Seg&x,const Seg&y){return x.p<y.p;}
inline void lsh()
{
sort(cpy+1,cpy+cnt+1);
cnt=unique(cpy+1,cpy+cnt+1)-cpy-1;
for(int i=1;i<cnt;++i)
len[i]=cpy[i+1]-cpy[i];
}
namespace SGT
{
#define lc (rt<<1)
#define rc (rt<<1|1)
db s[N<<3],t[N<<3];int ct[N<<3];
inline void up(int rt,int l,int r)
{
if(l==r)s[rt]=ct[rt]?t[rt]:0;
else s[rt]=ct[rt]?t[rt]:s[lc]+s[rc];
}
void build(int rt,int l,int r)
{
s[rt]=0;ct[rt]=0;
if(l==r){t[rt]=len[l];return;}
int mid=l+r>>1;
build(lc,l,mid),build(rc,mid+1,r);
t[rt]=t[lc]+t[rc];
}
void update(int rt,int l,int r,int ql,int qr,int v)
{
if(ql<=l&&r<=qr){ct[rt]+=v;return up(rt,l,r);}
int mid=l+r>>1;
if(ql<=mid)update(lc,l,mid,ql,qr,v);
if(mid<qr)update(rc,mid+1,r,ql,qr,v);
up(rt,l,r);
}
#undef lc
#undef rc
}
using namespace SGT;
int main()
{
int kase=0;
while(1)
{
scanf("%d",&n);if(!n)break;
cnt=tot=0;
for(int i=1;i<=n;++i)
{
db a,b,c,d;scanf("%lf%lf%lf%lf",&a,&b,&c,&d);
cpy[++cnt]=b,cpy[++cnt]=d;
e[++tot].l=b,e[tot].r=d,e[tot].p=a,e[tot].v=1;
e[++tot].l=b,e[tot].r=d,e[tot].p=c,e[tot].v=-1;
}
lsh();build(1,1,cnt-1);
sort(e+1,e+tot+1);db ans=0,pos=0;
for(int i=1;i<=tot;++i)
{
int l=lower_bound(cpy+1,cpy+cnt+1,e[i].l)-cpy;
int r=lower_bound(cpy+1,cpy+cnt+1,e[i].r)-cpy-1;
ans+=(e[i].p-pos)*s[1];
update(1,1,cnt-1,l,r,e[i].v);pos=e[i].p;
}
printf("Test case #%d\nTotal explored area: %.2f\n\n",++kase,ans);
}
return 0;
}
POJ 1177, IOI 1998, Picture
description
求\(n\) 个矩形的周长并。（所有矩形的边都和\(x\) 轴垂直或平行）
solution
与上一题类似，仍然考虑扫描线。
考虑每一个过程对最终答案的贡献。对于平行于\(x\) 轴的边，它们的贡献就是\(2\times\) 该过程长度$\times $ 该过程中被覆盖的不相交的线段个数。对于平行于\(y\) 轴的边，它们的贡献就是过程末被覆盖的总长度减去过程初被覆盖的总长度的绝对值。（这些东西说起来有点抽象，画个图就很好理解）。这些东西用线段树维护仍然是十分方便的。
time complexity
\(\mathcal O(n\log_2n)\)
code
#include<bits/stdc++.h>
using namespace std;
const int N=1e5+5;
int len[N<<1];
namespace SGT
{
int t[N<<3],s[N<<3],ct[N<<3],num[N<<3];bool lp[N<<3],rp[N<<3];
#define lc (rt<<1)
#define rc (rt<<1|1)
inline void up(int rt,int l,int r)
{
if(ct[rt])
s[rt]=t[rt],num[rt]=1,lp[rt]=rp[rt]=1;
else
{
if(l==r)
{
s[rt]=num[rt]=0;
lp[rt]=rp[rt]=0;
return;
}
s[rt]=s[lc]+s[rc];
num[rt]=num[lc]+num[rc]-(lp[rc]&rp[lc]);
lp[rt]=lp[lc],rp[rt]=rp[rc];
}
}
void build(int rt,int l,int r)
{
if(l==r){t[rt]=len[l];return;}
int mid=l+r>>1;
build(lc,l,mid),build(rc,mid+1,r);
t[rt]=t[lc]+t[rc];
}
void update(int rt,int l,int r,int ql,int qr,int v)
{
if(ql<=l&&r<=qr){ct[rt]+=v;return up(rt,l,r);}
int mid=l+r>>1;
if(ql<=mid)update(lc,l,mid,ql,qr,v);
if(mid<qr)update(rc,mid+1,r,ql,qr,v);
up(rt,l,r);
}
#undef lc
#undef rc
}
using namespace SGT;
int n,cpy[N<<1],cnt,tot;
inline void lsh()
{
sort(cpy+1,cpy+cnt+1);
cnt=unique(cpy+1,cpy+cnt+1)-cpy-1;
for(int i=1;i<cnt;++i)
len[i]=cpy[i+1]-cpy[i];
}
struct Seg{int l,r,p,v;}e[N<<1];
inline bool operator<(const Seg&x,const Seg&y){return x.p==y.p?x.v>y.v:x.p<y.p;}
int main()
{
scanf("%d",&n);
for(int i=1,a,b,c,d;i<=n;++i)
{
scanf("%d%d%d%d",&a,&b,&c,&d);
cpy[++cnt]=b,cpy[++cnt]=d;
e[++tot].l=b,e[tot].r=d,e[tot].p=a,e[tot].v=1;
e[++tot].l=b,e[tot].r=d,e[tot].p=c,e[tot].v=-1;
}lsh();
build(1,1,cnt-1);
sort(e+1,e+tot+1);int ans=0,pc=0,pre=0;
for(int i=1;i<=tot;++i)
{
int l=lower_bound(cpy+1,cpy+cnt+1,e[i].l)-cpy;
int r=lower_bound(cpy+1,cpy+cnt+1,e[i].r)-cpy-1;
ans+=2*num[1]*(e[i].p-pre);
update(1,1,cnt-1,l,r,e[i].v);
ans+=abs(s[1]-pc);
pc=s[1],pre=e[i].p;
}
printf("%d\n",ans);
return 0;
}
正方形(Squares, Seoul 2009, UVA 1453)
description
平面上给出若干个点，要求找到两个欧式距离相距最大的点。输出距离的平方。
solution
求出凸包然后旋转卡壳即可。
time complexity
\(\mathcal O(n\log_2n)\)
code
#include<bits/stdc++.h>
using namespace std;
const int N=1e5+5;
struct pt{int x,y;}p[N<<2],q[N<<2];
inline pt operator-(const pt&x,const pt&y){return {x.x-y.x,x.y-y.y};}
inline int cross(const pt&x,const pt&y){return x.x*y.y-x.y*y.x;}
inline bool cmp(const pt&x,const pt&y){return x.x==y.x?x.y<y.y:x.x<y.x;}
inline int Andrew(int n)
{
sort(p+1,p+n+1,cmp);
int k=0;
for(int i=1;i<=n;++i)
{
while(k>1&&cross(q[k]-p[i],q[k-1]-p[i])>=0)--k;
q[++k]=p[i];
}
int m=k;
for(int i=n-1;i;--i)
{
while(k>m&&cross(q[k]-p[i],q[k-1]-p[i])>=0)--k;
q[++k]=p[i];
}
return k-=n>1;
}
inline int S(const pt&x,const pt&y,const pt&z)
{
return abs(cross(z-x,y-x));
}
inline int len(const pt&x){return x.x*x.x+x.y*x.y;}
inline int nxt(int x,int n){return x==n?1:x+1;}
inline int solve(int n)
{
int ans=0;
for(int i=1,j=2;i<=n;++i)
for(;;j=nxt(j,n))
{
int dif=S(q[nxt(j,n)],q[i],q[nxt(i,n)])-S(q[j],q[i],q[nxt(i,n)]);
if(dif<=0)
{
ans=max(ans,len(q[j]-q[i]));
if(dif==0)ans=max(ans,len(q[nxt(j,n)]-q[i]));
break;
}
}
return ans;
}
int n;
int main()
{
int T;scanf("%d",&T);
while(T--)
{
scanf("%d",&n);int cnt=0;
for(int i=1,x,y,d;i<=n;++i)
{
scanf("%d%d%d",&x,&y,&d);
p[++cnt]={x,y};
p[++cnt]={x,y+d};
p[++cnt]={x+d,y};
p[++cnt]={x+d,y+d};
}
printf("%d\n",solve(Andrew(cnt)));
}
return 0;
}
离海最远的点(Most Distant Point from the Sea, Tokyo 2007, UVA 1396)
description
给出一个凸多边形，在其内部找出一个点使得它到凸多边形边界的距离最大。
solution
二分答案\(x\) ，而后问题转化为一个判定性问题。
将凸多边形的所有边向其内部移动\(x\) 的距离，如果最后仍然有面积则可行，而这实际上就是一次求半平面交的过程。
time complexity
\(\mathcal O(n\log_2n)\)
code
#include<bits/stdc++.h>
using namespace std;
const int N=205;
const double eps=1e-10;
typedef double db;
inline int sig(db x){return x>eps?1:(x<-eps?-1:0);}
struct pt{db x,y;}t[N];
inline pt operator+(const pt&x,const pt&y){return {x.x+y.x,x.y+y.y};}
inline pt operator-(const pt&x,const pt&y){return {x.x-y.x,x.y-y.y};}
inline pt operator*(const pt&x,db t){return {x.x*t,x.y*t};}
inline pt operator/(const pt&x,db t){return {x.x/t,x.y/t};}
inline db ang(const pt&x){return atan2(x.y,x.x);}
inline db len(const pt&x){return sqrt(x.x*x.x+x.y*x.y);}
inline pt rot90(const pt&x){return {-x.y,x.x};}
inline db cross(const pt&x,const pt&y){return x.x*y.y-x.y*y.x;}
struct Line{pt p,f;}li[N],q[N],o[N];
inline pt inr(const Line&x,const Line&y)
{
pt u=x.p-y.p;
db t=cross(u,y.f)/cross(y.f,x.f);
return x.p+x.f*t;
}
int n;
inline bool cmp(const Line&x,const Line&y)
{
db ag1=ang(x.f),ag2=ang(y.f);
if(sig(ag1-ag2)!=0)return ag1<ag2;
return sig(cross(y.p-x.p,y.p+y.f-x.p))>0;
}
inline bool ck(db x)
{
for(int i=1;i<=n;++i)
{
o[i]=li[i];
pt u=rot90(li[i].f);u=u/len(u);
o[i].p=li[i].p+u*x;
}
sort(o+1,o+n+1,cmp);
int cnt=0;
for(int i=1;i<=n;++i)
{
if(i!=1&&sig(ang(o[i].f)-ang(o[i-1].f))==0)continue;
o[++cnt]=o[i];
}
int l=1,r=0;
for(int i=1;i<=n;++i)
{
while(l<r&&sig(cross(o[i].p+o[i].f-t[r],o[i].p-t[r]))>=0)--r;
while(l<r&&sig(cross(o[i].p+o[i].f-t[l+1],o[i].p-t[l+1]))>=0)++l;
q[++r]=o[i];
if(r>l)t[r]=inr(q[r],q[r-1]);
}
while(l<r&&sig(cross(q[l].p+q[l].f-t[r],q[l].p-t[r]))>=0)--r;
t[r+1]=inr(q[r],q[l]);++r;
return r-l>=3;
}
int main()
{
while(1)
{
scanf("%d",&n);if(!n)break;
pt pre,st,now;
for(int i=1;i<=n;++i)
{
scanf("%lf%lf",&now.x,&now.y);
if(i==1)st=now;
else li[i]={pre,now-pre};pre=now;
}
li[1]={now,st-now};
db l=0,r=1e9;ck(1);
for(int i=1;i<=70;++i)
{
db mid=(l+r)/2;
ck(mid)?l=mid:r=mid;
}
printf("%.6lf\n",l);
}
return 0;
}
NO PAIN NO GAIN
