正月十一hu测 T2.阴影（计算几何【线面交点）

最新推荐文章于 2024-03-15 20:15:39 发布

Coco_T_

最新推荐文章于 2024-03-15 20:15:39 发布

阅读量274

点赞数

分类专栏：计算几何日常hu测文章标签：计算几何

本文链接：https://blog.csdn.net/wu_tongtong/article/details/79376694

版权

日常hu测同时被 2 个专栏收录

61 篇文章 0 订阅

订阅专栏

计算几何

42 篇文章 3 订阅

订阅专栏

这里写图片描述

分析：
30%还是比较好搞的，用相似三角形把点投影到xoy平面上，二维凸包即可

因为不会坐标系旋转，所以剩下的70%我用了一个三维凸包
先把点投射到地面上
用投影点和光源做一个三维凸包
$S={3*V \over h}$
（h就是点光源到地面的距离）

这里xue微说一下，三维投影时需要用到的：直线和平面的交点

首先我们需要把面的方向确定
（根据右手定则定义面的方向）
保证点光源和地面可以互相看见

//一个混合积就搞定了
if (Dot(o-face[0],Cross(face[1]-face[0],face[2]-face[0]))<eps) swap(face[1],face[2]);

我们想要求线面交点，首先需要用点法式表示平面
对于一个平面，我们可以讲其表示为 $(p_0,n)$
其中点 $p_0$ 是平面上一点， $n$ 是平面法向量
平面的法向量我们用Cross即可：

Node FA=Cross(face[1]-face[0],face[2]-face[0]);

设平面方程为 $Dot(n,p-p_0)=0$ ，过点 $p_1,p_2$ 的直线参数方程为 $p=p_1+t(p_2-p_1)$ ，则与平面方程联立得到：
$t=Dot(n,p_0-p_1)/Dot(n,p_2-p_1)$

Node jiao(Node p1,Node p2,Node p0,Node pp) {
    Node v=p2-p1;
    double t=Dot(pp,p0-p1)/Dot(pp,p2-p1);
    return p1+v*t;
}

#include<cstdio>
#include<cstring>
#include<iostream>
#include<cmath>
#include<algorithm>

using namespace std;

const int N=110;
const double eps=1e-8;
const double Pi=acos(-1.0);
struct node{
    double x,y;
    node (double xx=0,double yy=0) {
        x=xx;y=yy;
    }
};
node po[N];

node operator +(const node &a,const node &b) {return node(a.x+b.x,a.y+b.y);}
node operator -(const node &a,const node &b) {return node(a.x-b.x,a.y-b.y);}
node operator *(const node &a,const double &b) {return node(a.x*b,a.y*b);}
node operator /(const node &a,const double &b) {return node(a.x/b,a.y/b);}

double Dot(node a,node b) {return a.x*b.x+a.y*b.y;}
double Cross(node a,node b) {return a.x*b.y-a.y*b.x;}
double lenth(node a) {return sqrt(Dot(a,a));}

int dcmp(double x) {
    if (fabs(x)<eps) return 0;
    else if (x>0) return 1;
    else return -1;
}

struct Node{
    double x,y,z;
    Node (double xx=0,double yy=0,double zz=0)
    {
        x=xx;y=yy;z=zz;
    }
};
Node p[N],q[N];

Node operator +(const Node &A,const Node &B) {return Node(A.x+B.x,A.y+B.y,A.z+B.z);}
Node operator -(const Node &A,const Node &B) {return Node(A.x-B.x,A.y-B.y,A.z-B.z);}
Node operator *(const Node &A,const double &B) {return Node(A.x*B,A.y*B,A.z*B);}
Node operator /(const Node &A,const double &B) {return Node(A.x/B,A.y/B,A.z/B);}

double Dot(Node A,Node B) {return A.x*B.x+A.y*B.y+A.z*B.z;}
Node Cross(Node A,Node B) {return Node(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 lenth(Node a) {return sqrt(Dot(a,a));}

int n;
Node o,face[3],_;

struct CH3D{
    struct face{
        int a,b,c;   
        //凸包一个面上的三个点编号 
        bool ff;
        //该面是否属于最终凸包上的面 
    };

    int n;            //初始顶点数
    Node P[N]; 
    int num;          //凸包表面的三角形数
    face F[8*N];
    int vis[N][N];    

    //三角形面积*2 
    double S(Node a,Node b,Node c) {
        return lenth(Cross(b-a,c-a));    //注意要取lenth 
    }
    //四面体有向体积*6  混合积 
    double V(Node a,Node b,Node c,Node d) {
        return Dot(d-a,Cross(b-a,c-a));
    } 
    //面是否可见 
    int cansee(Node a,face f) {
        double x=Dot(a-P[f.a],Cross(P[f.b]-P[f.a],P[f.c]-P[f.a]));
        if (x>eps) return 1;
        else return 0;
    }

    //搜索与该边相邻的另一个平面
    void deal(int p,int a,int b) {
        int f=vis[a][b];
        face add;
        if (F[f].ff) {
            if (cansee(P[p],F[f]))
                dfs(p,f);
            else {
                add.a=b;
                add.b=a;
                add.c=p;
                add.ff=1;
                vis[p][b]=vis[a][p]=vis[b][a]=num;
                F[num++]=add;
            }
        }
    }
    //递归搜索所有应该从凸包内删除的面
    void dfs(int p,int now) {
        F[now].ff=0;
        deal(p,F[now].b,F[now].a);
        deal(p,F[now].c,F[now].b);
        deal(p,F[now].a,F[now].c);
    }
    //构建三维凸包
    void create() {
        face add;

        num=0;
        if (n<4) return;

        //保证前四个点不共面 
        bool flag=1;
        for (int i=1;i<n;i++) {
            if (lenth(P[0]-P[i])>eps) {
                swap(P[1],P[i]);
                flag=0; break;
            }
        }
        if (flag) return;
        flag=1;
        for (int i=2;i<n;i++) {
            if (lenth(Cross(P[1]-P[0],P[i]-P[0]))>eps) {
                swap(P[2],P[i]);
                flag=0; break;
            }
        }
        if (flag) return;
        flag=1;
        for (int i=3;i<n;i++) {
            if (fabs(Dot(P[i]-P[0],Cross(P[1]-P[0],P[2]-P[0])))>eps) {
                swap(P[3],P[i]);
                flag=0; break;
            }
        }
        if (flag) return;

        for (int i=0;i<4;i++)
        {
            add.a=(i+1)%4;
            add.b=(i+2)%4;
            add.c=(i+3)%4;
            add.ff=1;
            if (cansee(P[i],add)) swap(add.b,add.c);
            vis[add.a][add.b]=vis[add.b][add.c]=vis[add.c][add.a]=num;
            F[num++]=add;
        }
        for (int i=4;i<n;i++)
            for (int j=0;j<num;j++)
                if (F[j].ff&&cansee(P[i],F[j])){
                    dfs(i,j);
                    break;
                }
        int tmp=num;
        for (int i=num=0;i<tmp;i++)
            if (F[i].ff)
                F[num++]=F[i];
    }

    //体积
    double Vsum() {
        double ans=0;
        Node tmp(0,0,0);
        for (int i=0;i<num;i++)
            ans+=V(tmp,P[F[i].a],P[F[i].b],P[F[i].c]);
        return fabs(ans);
    }
    //点到面的距离
    double dis_face(Node p,Node a,Node b,Node c) {
        return fabs(V(a,b,c,p)/lenth(Cross(b-a,c-a)));
    }
};
CH3D hull;

int sta[N],top;

double SS()
{
    int i;
    double area=0;
    for (i=2;i<top;i++)  //这里要注意是i<top 
        area+=Cross(po[sta[i]]-po[sta[1]],po[sta[i+1]]-po[sta[1]]);
    return fabs(area)/2.0;
}

int cmp(const node &a,const node &b) {
    return a.x<b.x||(a.x==b.x&&a.y<b.y);
}

void TuB()
{
    sort(po+1,po+1+n,cmp);
    int i,j;
    top=0;
    for (i=1;i<=n;i++) {
        while (top>1&&dcmp(Cross(po[i]-po[sta[top-1]],po[sta[top]]-po[sta[top-1]]))<=0) top--;
        sta[++top]=i;
    }
    int k=top;
    for (i=n-1;i>=1;i--) {
        while (top>k&&dcmp(Cross(po[i]-po[sta[top-1]],po[sta[top]]-po[sta[top-1]]))<=0) top--;
        sta[++top]=i;
    }
    if (n>1) top--;
    if (top<3) printf("0");
    else printf("%0.2lf",SS());
    return;
}

void Shadow_2D() {    //po
    for (int i=1;i<=n;i++) {
        double za=fabs(o.z-p[i].z);
        double xa=o.x-p[i].x;
        double ya=o.y-p[i].y;
        po[i].x=p[i].x-xa*p[i].z/za;
        po[i].y=p[i].y-ya*p[i].z/za;
    }
}

Node jiao(Node p1,Node p2,Node p0,Node pp) {
    Node v=p2-p1;
    double t=Dot(pp,p0-p1)/Dot(pp,p2-p1);
    return p1+v*t;
}

void Shadow_3D() {    //hull.P
    hull.n=n;
    Node FA=Cross(face[1]-face[0],face[2]-face[0]);
    Node pp=face[0];   //面上一点     
    for (int i=1;i<=n;i++) 
        hull.P[i-1]=jiao(o,p[i],pp,FA);
    hull.P[n]=o;
    hull.n++;
}

void solve_2D() {
    Shadow_2D();
    TuB();
}

void solve_3D() {
    Shadow_3D();
    hull.create();
    double V=hull.Vsum();
    double x=hull.dis_face(o,face[0],face[1],face[2]);
    V=V/(2.0*x);
    printf("%0.2lf\n",V);
}

int main()
{
    freopen("shadow.in","r",stdin);
    freopen("shadow.out","w",stdout);

    for (int i=0;i<3;i++) scanf("%lf%lf%lf",&face[i].x,&face[i].y,&face[i].z);
    scanf("%lf%lf%lf",&o.x,&o.y,&o.z);
    if (Dot(o-face[0],Cross(face[1]-face[0],face[2]-face[0]))<eps) swap(face[1],face[2]);   //注意面的方向 

    scanf("%d",&n);
    for (int i=1;i<=n;i++) scanf("%lf%lf%lf",&p[i].x,&p[i].y,&p[i].z);

    if (face[0].z==0&&face[1].z==0&&face[2].z==0) solve_2D();
    else solve_3D();

    return 0;
}