Poj 2177 Ghost Busters（球交点三维计算几何模板）

待我称王封你为后i 2022-09-15 11:43 163阅读 0赞

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
    #else
    	freopen("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（球交点三维计算几何模板）

发表评论取消回复

还没有评论，来说两句吧...

相关阅读

相关三维计算几何之三维凸包增量法

相关 ACM计算几何模板（模板）

相关 Poj 2177 Ghost Busters（球交点三维计算几何模板）

相关 POJ 1408-Fishnet（计算几何-根据交点求多边形面积）

相关 POJ 1039-Pipe（计算几何-线段相交、求交点）

相关 poj-2504 (计算几何)

相关 POJ 2318(计算几何 )

相关 POJ 3528(计算几何+三维凸包表面积)

相关 POJ - 1113 Wall ————计算几何

相关 POJ - 2318 TOYS————计算几何