LA 3616 How I Wonder What You Are! 三维几何 角度计算

题目地址：LA3616

直接O(nm)暴力枚举..  简单题

代码：

#include<iostream>
#include<cmath>
#include<cstdio>
#include<iomanip>
#include<ctime>
#include<cstdlib>
#include<vector>

using namespace std;
const double eps=1e-8;

// 1e-8?

const double PI=acos(-1.0);

int dcmp(double x)
{
    return (x>eps)-(x<-eps);
}

struct Point3
{
    double x,y,z;
    Point3(double x=0,double y=0,double z=0):x(x),y(y),z(z) {}
    
};

typedef Point3 Vector3;

Vector3 operator+(Point3 A,Point3 B)
{
    return Point3(A.x+B.x,A.y+B.y,A.z+B.z);
    
}

Vector3 operator-(Point3 A,Point3 B)
{
    return Point3(A.x-B.x,A.y-B.y,A.z-B.z);
    
}

Vector3 operator*(Point3 A,double p)
{
    return Point3(A.x*p,A.y*p,A.z*p);
    
}

Vector3 operator/(Point3 A,double p)
{
    return Point3(A.x/p,A.y/p,A.z/p);
    
}

bool operator<(Point3 &A,Point3 &B)
{
    return A.x<B.x||(A.x==B.x&&A.y<B.y)||(A.x==B.x&&A.y==B.y&&A.z<B.z);
    
}
bool operator==(Point3 &A,Point3 &B)
{
    return dcmp(A.x-B.x)==0&&dcmp(A.y-B.y)==0&&dcmp(A.z-B.z)==0;
}
ostream & operator<<(ostream & out,Vector3 A)
{
    out<<A.x<<' '<<A.y<<' '<<A.z<<endl;
    
    return out;
}

Point3 read_point3()
{
    Point3 p;
    scanf("%lf%lf%lf",&p.x,&p.y,&p.z);
    return p;
    
}


double Dot(Vector3 A,Vector3 B)
{
    return  A.x*B.x+A.y*B.y+A.z*B.z;
}

double Length(Vector3 A)
{
    return sqrt(Dot(A,A));
}

double Angle(Vector3 A,Vector3 B)
{
    return acos(Dot(A,B)/Length(A)/Length(B));
}


Vector3 Cross(Vector3 A,Vector3 B)   // 真正的叉积
{
    return  Vector3(A.y*B.z-A.z*B.y,A.z*B.x-A.x*B.z,A.x*B.y-A.y*B.x);
}


// 面积的两倍
double Area2(Point3 A,Point3 B,Point3 C)
{
    return Length(Cross(B-A,C-A));
}


double DistanceToPlane(const Point3 &p,const Point3 & p0,const Vector3 &n)
{
    return fabs(Dot(p-p0,n))/Length(n);
    
}

//  n不用管方向

Point3  GetPlaneProjection(Point3 p,Point3 p0,Vector3 n)
{
    return p-n*Dot(p-p0,n)/Dot(n,n);
}

//   三点式
Point3  GetPlaneProjection(Point3 p,Point3 p0,Point3 p1,Point3 p2)
{
       Vector3 n=Cross(p1-p0,p2-p0);
    
       return p-n*Dot(p-p0,n)/Dot(n,n);
}
Point3 LinePlaneIntersection(Point3 p1,Point3 p2,Point3 p0,Vector3 n)
{
    double t=Dot(n,p0-p1)/Dot(n,p2-p1);
    return p1+(p2-p1)*t;
    
}

//  点在三角形内部
bool   PointInTri(Point3 p,Point3 p0,Point3 p1,Point3 p2)
{
    double area1=Area2(p, p0, p1);
    double area2=Area2(p, p1, p2);
    double area3=Area2(p, p2, p0);
    
    return dcmp(area1+area2+area3-Area2(p0,p1,p2))==0;
}

//  线段与三角形相交 不考虑共面情形
bool  TriSegIntersection(Point3 p0,Point3 p1,Point3 p2,Point3 A,Point3 B,Point3 & p)
{
    Vector3 n=Cross(p1-p0, p2-p0);
    if(dcmp(Dot(n,B-A))==0)  return false;
    else{
        
        double t=Dot(n,p0-A)/Dot(n,B-A);
        if(dcmp(t)<0||dcmp(t-1)>0)  return false;
        else
        {
            p=A+(B-A)*t;
            return PointInTri(p,p0,p1,p2);
        }
        
    }
}

double DistanceToLine(Point3 P,Point3 A,Point3 B)
{
    Vector3 v1=B-A;
    Vector3 v2=P-A;
    return  Length(Cross(v1, v2))/Length(v1);
    
}

double  DistanceToSegment(Point3 P,Point3 A,Point3 B)
{
    if(A==B)  return Length(P-A);
    Vector3 v1=B-A,v2=P-A,v3=P-B;
    if(dcmp(Dot(v1,v2))<0) return Length(v2);
    else if(dcmp(Dot(v1,v3))>0) return Length(v3);
    else  return Length(Cross(v1,v2))/Length(v1);
    
}

double volume6(Point3 A,Point3 B,Point3 C,Point3 D)
{
    return  Dot(D-A,Cross(B-A,C-A));
}

// 0-1随机数

double rand01() { return rand()/(double)RAND_MAX;}
double randeps() {
    //srand((unsigned) time(NULL));
    return  (rand01()-0.5)*eps;}

Point3 add_noise(Point3 P)
{
    return Point3(P.x+randeps(),P.y+randeps(),P.z+randeps());
    
}


// 三维凸包  增量法

struct Face {
    int v[3];
    Face(int a, int b, int c) { v[0] = a; v[1] = b; v[2] = c; }
    Vector3 Normal(const vector<Point3>& P) const {
        return Cross(P[v[1]]-P[v[0]], P[v[2]]-P[v[0]]);
    }
    // f是否能看见P[i]
    int CanSee(const vector<Point3>& P, int i) const {
        return Dot(P[i]-P[v[0]], Normal(P)) > 0;
    }
    
};

// 增量法求三维凸包
// 注意：没有考虑各种特殊情况（如四点共面）。实践中，请在调用前对输入点进行微小扰动
vector<Face> CH3D(const vector<Point3>& P) {
    int n = P.size();
    vector<vector<int> > vis(n);
    for(int i = 0; i < n; i++) vis[i].resize(n);
    
    vector<Face> cur;
    cur.push_back(Face(0, 1, 2)); // 由于已经进行扰动，前三个点不共线
    cur.push_back(Face(2, 1, 0));
    for(int i = 3; i < n; i++) {
        vector<Face> next;
        // 计算每条边的“左面”的可见性
        for(int j = 0; j < cur.size(); j++) {
            Face& f = cur[j];
            int res = f.CanSee(P, i);
            if(!res) next.push_back(f);
            for(int k = 0; k < 3; k++) vis[f.v[k]][f.v[(k+1)%3]] = res;
        }
        for(int j = 0; j < cur.size(); j++)
            for(int k = 0; k < 3; k++) {
                int a = cur[j].v[k], b = cur[j].v[(k+1)%3];
                if(vis[a][b] != vis[b][a] && vis[a][b]) // (a,b)是分界线，左边对P[i]可见
                    next.push_back(Face(a, b, i));
            }
        cur = next;
    }
    return cur;
    
}

double DistanceToPlane( const Point3 &p,const Point3 & p0,const  Point3 & p1,const  Point3 & p2)
{
    return fabs(volume6(p,p0,p1,p2)/Area2(p0, p1, p2));
    
}

//  四面体质心  就是四个点的平均
Point3 Centroid(const Point3 & A,const Point3 & B,const Point3 & C,const Point3 & D)
{
    return (A+B+C+D)/4.0;
}


struct  ConvexPolyhedron
{
    int n;  // 总顶点数
    vector<Point3>  P,P2;  //  P是真实顶点， P1是扰乱后的顶点
    vector<Face>  faces;
    
    void init(int n)
    {
        this->n=n;
        P.resize(n);
        P2.resize(n);
        
        for(int i=0;i<n;i++)
        {
            P[i]=read_point3();
            P2[i]=add_noise(P[i]);
        }
        
        
        faces=CH3D(P2);
        
    }
    
    // 本题需要,写在类里面
    
    
    double mindist(Point3 C)
    {
        double ans=1e30;
        for(int i=0;i<faces.size();i++)
        {
            ans=min(ans,DistanceToPlane(C,P[faces[i].v[0]], P[faces[i].v[1]], P[faces[i].v[2]]));  // 都是P这个数组
            
        }
        return ans;
    }
    
    
    Point3 centroid()
    {
        Point3 C=Point3(0,0,0);
        
        Point3 tot=Point3(0,0,0);
        
        double totv=0;
        
        for(int i=0;i<faces.size();i++)
        {
            Point3  p0=P[faces[i].v[0]];
            Point3  p1=P[faces[i].v[1]];
            Point3  p2=P[faces[i].v[2]];
            
            double v=volume6(p0, p1, p2, C);
            tot=tot+Centroid(p0,p1,p2,C)*v;
            totv+=v;
            
        }
        
        return tot/totv;
        
    }
    
};



//  self
// p0,p1在面A,B,C 的同侧
bool SameSide(Point3 p0,Point3 p1,Point3 A,Point3 B,Point3 C)
{
    Vector3 n=Cross(B-A, C-A);
    Vector3 v1=p0-A;
    Vector3 v2=p1-A;
    
    return Dot(v1, n)*Dot(v2,n)>=0;
}

// 改写读点函数

bool  read_point3(Point3 & P)
{
   if(scanf("%lf%lf%lf",&P.x,&P.y,&P.z)==3)
   {
       return 1;
   }
   else return 0;
    
}


//   求质心
Point3 Centroid(Point3 &A,Point3 & B,Point3 &C,Point3 &D)
{
    return (A+B+C+D)/4.0;
}

Point3 Centroid(Point3 &A,Point3 & B,Point3 &C,Point3 &D,Point3 &E)  //  注意A,B,C 为四面体公共平面
{
    double  v1=volume6(A, B, C, D);
    double  v2=volume6(A, B, C, E);
    v1=fabs(v1);
    v2=fabs(v2);
    
    Point3 C1=Centroid(A, B, C, D);
    Point3 C2=Centroid(A, B, C, E);
  
    
    return (C1*v1+C2*v2)/(v1+v2);

}

bool  OnSamePlane(Point3 &A,Point3 & B,Point3 &C,Point3 &D)
{

    return dcmp(volume6(A, B, C, D))==0;
}

int n,m;
Point3 P[505];
Point3 T[55];
double T_ang[55];

bool  cansee(int x)
{
    for(int i=0;i<m;i++)
    {
        if(dcmp(Angle(P[x], T[i])-T_ang[i])<0)
        {
            return 1;
        }
    }
    return  0;
}



int main()
{
    while(cin>>n)
    {
        if(!n)  break;
        for(int i=0;i<n;i++)
        {
            P[i]=read_point3();
            
        }
        cin>>m;
        for(int i=0;i<m;i++)
        {
            T[i]=read_point3();
            scanf("%lf",&T_ang[i]);
        }
        
        int  cnt=0;
        for(int i=0;i<n;i++)
        {
           if(cansee(i))
           {
               cnt++;
           }
        }
        cout<<cnt<<endl;
        
    }
}