Poj 2177 Ghost Busters（球交点 三维计算几何模板）
来源：互联网 发布：脸部表情软件 编辑：程序博客网 时间：2024/06/04 19:34
2015-4-27更新：由于百度空间即将关闭，故将提到的参考文章复制到了最下面。
题意：给出三维空间第一卦限内若干个个球。从原点引出一条射线，问最多能够穿过多少个球？
思路参考：http://hi.baidu.com/lccycc_acm/item/7ab19fd59d39c61d21e250d2
从原点看过去，球就是一个个圆。任取两个球，将其放到同一个距离上。显然如果两圆不相交，那么射线只可能是射球心，如果两圆相交，则通过两个圆的交点。
设两球心为p1,p2 先算出两个交点q1,q2的连线和连心线的交点（这两条线是垂直的），设为m。于是Om和p1p2,q1q2三条线就两两垂直了。用Om 叉乘p1p2,得到q1q2的向量，而|q1m|=|q2m|也可以算，于是q1,q2就算出来了。
#pragma warning(disable:4786)#include <set>#include <map>#include <cmath>#include <queue>#include <stack>#include <string>#include <vector>#include <cstdio>#include <cstdlib>#include <cstring>#include <sstream>#include <iostream>#include <algorithm>using namespace std;#define lowbit(x) ((x)&(-(x)))#define max(a,b) ((a)>(b)?(a):(b))#define max3(a,b,c) (max(a,max(b,c)))#define min(a,b) ((a)<(b)?(a):(b))const int N=105;const double dinf=1e15;const double EPS=1e-7;struct point3{ double x,y,z; point3(){} point3(double _x,double _y,double _z) {x=_x;y=_y;z=_z;} void Get () {scanf("%lf%lf%lf",&x,&y,&z);} point3 operator+(point3 a) { return point3(x+a.x,y+a.y,z+a.z); } point3 operator-(point3 a) { return point3(x-a.x,y-a.y,z-a.z); } point3 operator*(point3 a) //叉积 { return point3(y*a.z-z*a.y,z*a.x-x*a.z,x*a.y-y*a.x); } point3 operator*(double t) { return point3(x*t,y*t,z*t); } double operator^(point3 a) //点积 { return x*a.x+y*a.y+z*a.z; } point3 operator/(double t) { return point3(x/t,y/t,z/t); } double len () //长度 { return sqrt(x*x+y*y+z*z); } point3 adjust (double L) //调整为L长 { double t=len(); L/=t; return point3(x*L,y*L,z*L); } void print () { printf("%.10lf %.10lf %.10lf\n",x+EPS,y+EPS,z+EPS); }};int sgn (double x){ if (fabs(x) < EPS)return 0; return x>0?1:-1;}double len(point3 a){ return a.len();}double getArea(point3 a,point3 b,point3 c) //三角形面积{ double x=len((b-a)*(c-a)); return x/2;}double getVolume(point3 a,point3 b,point3 c,point3 d) //四面体体积{ double x=(b-a)*(c-a)^(d-a); return x/6;}point3 pShadowOnPlane(point3 p,point3 a,point3 b,point3 c) //点在平面的投影{ point3 v=(b-a)*(c-a); if(sgn(v^(a-p))<0) v=v*-1; v=v.adjust(1); double d=fabs(v^(a-p)); return p+v*d;}double lineToLine (point3 a,point3 b,point3 p,point3 q) //异面直线间距离{//参数：直线上两点 point3 v=(b-a)*(q-p); return fabs((a-p)^v)/v.len();}//两异面直线距离：两直线上的点的连线在其法向量上的投影double LineToLine (point3 p1,point3 k1,point3 p2,point3 k2){//参数：直线上一点及其法向量 point3 nV=k1*k2; return fabs(nV^(p1-p2))/nV.len();}int pInPlane(point3 p,point3 a,point3 b,point3 c) //点是否在面上{ double S=getArea(a,b,c); double S1=getArea(a,b,p); double S2=getArea(a,c,p); double S3=getArea(b,c,p); return sgn(S-S1-S2-S3)==0;}int opposite(point3 p,point3 q,point3 a,point3 b,point3 c) //p q在平面两侧{ point3 v=(b-a)*(c-a); double x=v^(p-a); double y=v^(q-a); return sgn(x*y)<0; //投影异号}int segCrossTri(point3 p,point3 q,point3 a,point3 b,point3 c) //pq是与三角形abc相交{ return opposite(p,q,a,b,c)&&opposite(a,b,p,q,c)&& opposite(a,c,p,q,b)&& opposite(b,c,p,q,a);}double pToPlane(point3 p,point3 a,point3 b,point3 c) //p到面abc的距离{ double v=((b-a)*(c-a)^(p-a))/6; double s=len((b-a)*(c-a))/2; return fabs(3*v/s);}double pToLine(point3 p,point3 a,point3 b) //点到直线的距离{ double S=len((a-p)*(b-p)); return S/len(a-b);}double pToSeg(point3 p,point3 a,point3 b) //点到线段的距离{ if(sgn((p-a)^(b-a))<=0) return len(a-p); if(sgn((p-b)^(a-b))<=0) return len(b-p); return pToLine(p,a,b);}double pToPlane1(point3 p,point3 a,point3 b,point3 c) //点到三角形的距离{ point3 k=pShadowOnPlane(p,a,b,c); if(pInPlane(k,a,b,c)) return pToPlane(p,a,b,c); double x=pToSeg(p,a,b); double y=pToSeg(p,a,c); double z=pToSeg(p,b,c); return min(x,min(y,z));}double getAng(point3 a,point3 b) //求两向量的夹角{ double x=(a^b)/len(a)/len(b); return acos(x);}double segToSeg(point3 a,point3 b,point3 p,point3 q) //线段到线段的距离{ point3 v=(b-a)*(q-p); double A,B,A1,B1; A=((b-a)*v)^(p-a); B=((b-a)*v)^(q-a); A1=((p-q)*v)^(a-q); B1=((p-q)*v)^(b-q); if(sgn(A*B)<=0&&sgn(A1*B1)<=0) { return lineToLine(a,b,p,q); } double x=min(pToSeg(a,p,q),pToSeg(b,p,q)); double y=min(pToSeg(p,a,b),pToSeg(q,a,b)); return min(x,y);}struct face{ int a,b,c,ok; face(){} face(int _a,int _b,int _c,int _ok) { a=_a; b=_b; c=_c; ok=_ok; }};struct _3DCH //三维凸包{ face F[N<<2]; int b[N][N],cnt,n; point3 p[N]; int getDir(point3 t,face F) { double x=(p[F.b]-p[F.a])*(p[F.c]-p[F.a])^(t-p[F.a]); return sgn(x); } void deal(int i,int x,int y) { int f=b[x][y]; if(!F[f].ok) return; if(getDir(p[i],F[f])==1) DFS(i,f); else { b[y][x]=b[x][i]=b[i][y]=cnt; F[cnt++]=face(y,x,i,1); } } void DFS(int i,int j) { F[j].ok=0; deal(i,F[j].b,F[j].a); deal(i,F[j].c,F[j].b); deal(i,F[j].a,F[j].c); } void construct() { int i,j,k=0; for(i=1;i<n;i++) if(sgn(len(p[i]-p[0]))) //找出两个点不共点 { swap(p[i],p[1]); k++; break; } if(k!=1) return; for(i=2;i<n;i++) if(sgn(getArea(p[0],p[1],p[i]))) //找出三点不共线 { swap(p[i],p[2]); k++; break; } if(k!=2) return; for(i=3;i<n;i++) if(sgn(getVolume(p[0],p[1],p[2],p[i]))) //找出四点不共面 { swap(p[i],p[3]); k++; break; } if(k!=3) return; cnt=0;for (i=0;i<4;i++) { face k=face((i+1)%4,(i+2)%4,(i+3)%4,1); if(getDir(p[i],k)==1) swap(k.b,k.c); b[k.a][k.b]=b[k.b][k.c]=b[k.c][k.a]=cnt; F[cnt++]=k; } for (i=4;i<n;i++) for (j=0;j<cnt;j++) { if(F[j].ok&&getDir(p[i],F[j])==1) { DFS(i,j); break; } } j=0;for (i=0;i<cnt;i++)if(F[i].ok) F[j++]=F[i]; cnt=j; } point3 getCenter() //求三维凸包的重心 { point3 ans=point3(0,0,0),o=point3(0,0,0); double s=0,temp; int i; for (i=0;i<cnt;i++) { face k=F[i]; temp=getVolume(o,p[k.a],p[k.b],p[k.c]); ans=ans+(o+p[k.a]+p[k.b]+p[k.c])/4*temp; s+=temp; } ans=ans/s; return ans; } double getMinDis(point3 a) //求凸包内一个点到面的最短距离 { double ans=dinf; int i; for (i=0;i<cnt;i++) { face k=F[i]; ans=min(ans,pToPlane(a,p[k.a],p[k.b],p[k.c])); } return ans; }};///////////////////////////////////////////int n;point3 data[N];double R[N];vector<point3> V; //储存直线有可能通过的点void Deal (point3 a,double ra,point3 b,double rb){//参数：球心，半径 //球变圆,折算到同一距离 double x=len(a),y=len(b); //球心到原点距离 b=b/y*x; rb=rb/y*x; double d=len(a-b); //连心线长度 if (sgn(d-ra-rb)==0) //相切 { V.push_back(a+(b-a)/d*ra); return; } if(sgn(d-ra-rb)>0) return; //相离 if(sgn(d-fabs(ra-rb))<=0) return; //内含 x=(ra*ra+d*d-rb*rb)/(2*d);//余弦定理 y=sqrt(ra*ra-x*x); point3 M=a+(b-a)/d*x;//两交点中点所在位置 point3 v=a*M; //叉积求得两交点所在直线的向量 v=v.adjust(y)+M; V.push_back(v); V.push_back(M*2-v);//平行四边形法则}int Judge (point3 a,point3 p,double r){//球：球心a，半径r与直线op的关系 double x=pToLine(p,point3(0,0,0),a); return sgn(x-r)<=0;}int main (){#ifdef ONLINE_JUDGE#elsefreopen("read.txt","r",stdin);#endif scanf("%d",&n); int i,j; for (i=1;i<=n;i++) { data[i].Get(); scanf("%lf",&R[i]); V.push_back(data[i]); } for (i=1;i<=n;i++) for (j=i+1;j<=n;j++) Deal(data[i],R[i],data[j],R[j]); int ans=0,size=V.size(),temp,k;//k记录选中的点 for (i=0;i<size;i++) { temp=0; for (j=1;j<=n;j++) temp+=Judge(V[i],data[j],R[j]); if (temp>ans) ans=temp,k=i; } printf("%d\n",ans); //枚举所有球，判断其与选定直线是否相交 for (i=1;i<=n && ans;i++) if (Judge (V[k],data[i],R[i])) { ans--; printf(ans!=0?"%d ":"%d\n",i); } return 0;}空间中给一堆球 求一条射线能射过的最多球数
题目比较人性化 球的三个坐标都是正的，而且也只要考虑正象限空间的那部分方向就可以了
从原点看过去，球就是一个个圆。任取两个球，将其放到同一个距离上，很容易就想到射线只可能是射球心，或者那两个圆的交点
求圆的交点和平面的求法一样。设两球心为p1,p2 先算出两个交点q1,q2的连线和连心线的交点（这两条线是垂直的），设为m。于是Om和p1p2,q1q2三条线就两两垂直了。用Om 叉乘p1p2,得到q1q2的向量，而|q1m|=|q2m|也可以算，于是q1,q2就算出来了。
wa了两次 数据中会有输入球数是0的要注意。还有就是由于两两圆的交点可以有n^2个，数组没开够大导致越界。
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <cmath>
using namespace std;
const double eps=1e-8;
int sgn(double a){return a>eps?1:(a<-eps?-1:0);}
struct Pt3
{
double x,y,z;
Pt3(){x=y=z=0;}
void init(){scanf("%lf%lf%lf",&x,&y,&z);}
Pt3 operator +(Pt3 a){a.x+=x;a.y+=y;a.z+=z;return a;}
Pt3 operator -(Pt3 a){a.x=x-a.x;a.y=y-a.y;a.z=z-a.z;return a;}
};
Pt3 operator *(Pt3 a,double d){a.x*=d;a.y*=d;a.z*=d;return a;}
Pt3 operator /(Pt3 a,double d){a.x/=d;a.y/=d;a.z/=d;return a;}
Pt3 mult(Pt3 a,Pt3 b){Pt3 r;r.x=a.y*b.z-a.z*b.y;r.y=a.z*b.x-a.x*b.z;r.z=a.x*b.y-a.y*b.x;return r;}
double dj(Pt3 a,Pt3 b){return a.x*b.x+a.y*b.y+a.z*b.z;}
double leng(Pt3 a){return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);}
//////////////////////////////////////////
const int maxn=110;
Pt3 p[maxn];
int n;
double r[maxn];
Pt3 a[maxn*maxn*4];
int aa;
void ct(Pt3 p1,double r1,Pt3 p2,double r2)
{
double dis1=leng(p1);
double dis2=leng(p2);
p2=p2/dis2*dis1;
r2=r2/dis2*dis1;
double d=leng(p1-p2);
if (sgn(d-fabs(r1+r2))==0)
{
Pt3 h;
h=(p1-p2)/d*r2+p2;
a[aa++]=h;
return ;
}
if (sgn(d-fabs(r1+r2))>0)
return ;
if (sgn(d-fabs(r1-r2))<=0)
return ;
Pt3 m;
double x= (r2*r2-r1*r1+d*d)/(2*d);
m=(p1-p2)/d*x+p2;
double h=sqrt(r2*r2-x*x);
Pt3 tw=mult(m,p2);
Pt3 q1=tw/leng(tw)*h+m;
Pt3 q2=m*2-q1;
a[aa++]=q1;
a[aa++]=q2;
}
int shoot(Pt3 tw,Pt3 h,double r)
{
tw=tw/leng(tw)*leng(h)/2;
double dtw=leng(tw);
double dh=leng(h);
double d=leng(tw-h);
double y=(dh*dh-d*d-dtw*dtw)/(2*dtw);
double x=sqrt(d*d-y*y);
return sgn(r-x)>=0;
}
int main()
{
int i,j,k;
scanf("%d",&n);
for (i=0;i<n;i++)
{
p[i].init();
scanf("%lf",r+i);
}
aa=0;
for (i=0;i<n;i++)
for (j=i+1;j<n;j++)
ct(p[i],r[i],p[j],r[j]);
for (i=0;i<n;i++)
a[aa++]=p[i];
int maxs=0,maxa=-1;
for (i=0;i<aa;i++)
{
int tot=0;
for (j=0;j<n;j++)
tot+=shoot(a[i],p[j],r[j]);
if (tot>maxs)
{
maxs=tot;
maxa=i;
}
}
printf("%d\n",maxs);
for (j=0;j<n;j++) if (maxs)
if (shoot(a[maxa],p[j],r[j]))
{
printf("%d",j+1);
maxs--;
if (maxs>0)
printf(" ");
}
puts("");
return 0;
}
- Poj 2177 Ghost Busters（球交点 三维计算几何模板）
- POJ 1264 SCUD Busters 计算几何
- 计算几何三维模板
- 三维计算几何模板
- 三维计算几何模板
- 三维计算几何模板整理
- 三维计算几何模板整理
- POJ 1039-Pipe（计算几何-线段相交、求交点）
- 计算几何模板 (三维几何函数库)
- 计算几何较全模板（有包含三维的）
- HDU 1086 计算几何 求线段交点（吉大模板）
- POJ 2826 计算几何模板
- [模板]三维几何模板
- 计算几何 POJ 2826 An Easy Problem?! （线段位置判断并且求交点）
- poj 3340 Segments（计算几何，直线跟线段的交点）
- POJ 2653 Pick-up sticks（计算几何 求线段交点）
- poj 1408 Fishnet（计算几何 叉积求面积 求两直线交点 暴力）
- POJ 1408-Fishnet（计算几何-根据交点求多边形面积）
- 百度质量部实习生面试
- 增强现实技术原理
- Terracotta设计原理分析--（部分内容来自官方描述）
- 面向对象的三个基本特征
- unresolved external symbol “symbol”（不确定的外部“符号”）。
- Poj 2177 Ghost Busters（球交点 三维计算几何模板）
- 数据结构1 顺序列表
- HDU 3790 最短路径问题
- char* 和 char[]的区别
- Dijkstra(迪杰斯特拉)算法
- 就算所有人都放弃我，我也不会放弃我自己，坚持梦想，拒绝国企，拒绝公务员
- 15个实用的Linux find命令示例
- JAVA反射
- 市场格言
