「2017 山东一轮集训 Day2 / SDWC2018 Day1」Shadow 题解

可恶的计算几何

思路

Subtask1

先考虑部分分，设 $A,B,C$ 是题目给定的确定平面的 3 个点，此时这三个点的 $z$ 相等。

那么此时可以先将所有点的 $z$ 坐标减去 $A$ 的 $z$ 坐标，此时这个平面即为 $z$ 坐标为 0 水平面。

此时投影即为光源 $L$ 与所有凸多面体的点 $p_i$ 的连线与水平面的交点，求交点很简单，由于水平面 $z$ 坐标为 0，求出直线解析式的常数项即可。此时阴影即为交点形成的凸包，求面积可以剖成三角形用二位叉积求面积，最后记得 $\times \frac{1}{2}$。

Subtask1 代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const double eps=1e-8;
int n;
struct N{
	double x,y,z;
}a[110]; 
struct P{
	double x,y;
}p[110]; 
double cj(P t,P a,P b){
	return (a.x-t.x)*(b.y-t.y)-(a.y-t.y)*(b.x-t.x);
}
bool cmp(P a,P b){
	if(fabs(a.x-b.x)>eps)return a.x<b.x;
	return a.y<b.y;
}
P stk[110];
int sti;
int main(){
	ios::sync_with_stdio(0);
	cin.tie(0);
	double Z;
	cin>>Z>>Z>>Z>>Z>>Z>>Z>>Z>>Z>>Z;
	cin>>a[0].x>>a[0].y>>a[0].z;
	a[0].z-=Z;
	cin>>n;
	for(int i=1;i<=n;i++){
		cin>>a[i].x>>a[i].y>>a[i].z;
		a[i].z-=Z; 
		p[i].x=a[i].x-(a[i].x-a[0].x)/(a[i].z-a[0].z)*a[i].z;
		p[i].y=a[i].y-(a[i].y-a[0].y)/(a[i].z-a[0].z)*a[i].z;
	}
	sort(p+1,p+1+n,cmp);
	for(int i=1;i<=n;i++){
		while(sti>1&&cj(stk[sti-1],stk[sti],p[i])<=eps)sti--;
		stk[++sti]=p[i];
	}
	vector<P> v;
	for(int i=1;i<=sti;i++)v.push_back(stk[i]);
	sti=0;
	for(int i=n;i;i--){
		while(sti>1&&cj(stk[sti-1],stk[sti],p[i])<=eps)sti--;
		stk[++sti]=p[i];
	}
	for(int i=2;i<sti;i++)v.push_back(stk[i]);
	double ans=0;
	for(int i=1;i<v.size()-1;i++)ans+=fabs(cj(v[0],v[i],v[i+1]));
	printf("%.2lf",ans/2);
	return 0;
}

正解

现在主要难点是求光源 $L$ 与顶点的连线与平面的交点。

我们首先要求出平面的正交基，设 $\vec{ZA}=\vec{CA},\vec{ZC}=\vec{CA}\times\vec{CB},\vec{ZB}=\vec{ZC}\times\vec{ZA}$，则 $\{\vec{ZA},\vec{ZB}\}$ 即为平面的正交基，$\vec{ZC}$ 为平面的法向量。

设 $Ld=\vec{AL}\cdot\vec{ZC}$，可以理解为 $\text{光源至平面的垂直距离}\times|\vec{ZC}|$。

接下来要求光源 $L$ 与每个顶点 $a$ 与平面交点的三维坐标。

首先由于交点一定在 $L,a$ 所在直线上，所以交点 $b$ 有 $b=L+t(a-L)$。

又由于 $b$ 在平面上，有 $\vec{Ab}\cdot\vec{ZC}=0$。

等价于

$$(b-A)\cdot\vec{ZC}=0$$

代入 $b=L+t(a-L)$

$$(L+t(a-L)-A)\cdot\vec{ZC}=0$$

拆开来

$$(L-A)\cdot\vec{ZC}+t(a-L)\cdot\vec{ZC}=0$$

注意到

$$(a-L)=(a-A)-(L-A)$$

所以有

$$Ld+t((a-A)\cdot\vec{ZC}-Ld)=0$$ $$t=\frac{Ld}{Ld-(a-A)\cdot\vec{ZC}}$$

代码：t=Ld/(Ld-dot(a-A,ZC));

那么我们就求出了 $b$ 的坐标，现在要将 $b$ 降至二维坐标的点 $P$。

由于 $\{\vec{ZA},\vec{ZB}\}$ 是平面的正交基，$C$ 为原点，则 $P$ 的坐标 $(x,y)$ 可以表示为 $\vec{CP}=x\frac{\vec{ZA}}{|\vec{ZA}|}+y\frac{\vec{ZB}}{|\vec{ZB}|}$（注意$\vec{ZA},\vec{ZB}$不是单位向量所以要除以向量的模长）。

首先对等式两边同乘 $\vec{ZA}$（内积）

$$\vec{CP}\cdot\vec{ZA}=x\frac{\vec{ZA}\cdot\vec{ZA}}{|\vec{ZA}|}+y\frac{\vec{ZB}\cdot\vec{ZA}}{|\vec{ZB}|}$$

有 $\vec{ZA}\cdot\vec{ZA}=|\vec{ZA}|^2$，由于 $\vec{ZA}$ 和 $\vec{ZB}$ 垂直，有 $\vec{ZA}\cdot\vec{ZB}=0$，则原式可以化为

$$\vec{CP}\cdot\vec{ZA}=x|\vec{ZA}|$$ $$x=\frac{\vec{CP}\cdot\vec{ZA}}{|\vec{ZA}|}$$

同理可得

$$y=\frac{\vec{CP}\cdot\vec{ZB}}{|\vec{ZB}|}$$

代码：p[i]={dot(b,ZA)/lenA,dot(b,ZB)/lenB};

那么我们就求出了 $P$ 的坐标，最后按照Subtask1的方法求凸包面积即可。

代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const double eps=1e-8;
int n;
struct N{double x,y,z;};//向量{x,y,z}（也可以表示点） 
N operator-(N a,N b){return {a.x-b.x,a.y-b.y,a.z-b.z};}
N cj3(N a,N b){return {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 dot(N a,N b){return a.x*b.x+a.y*b.y+a.z*b.z;}//向量的内积 
double dis(N a){return sqrt(dot(a,a));}//向量的模长 
struct P{double x,y;}p[110];//二维点坐标 
double cj(P t,P a,P b){return (a.x-t.x)*(b.y-t.y)-(a.y-t.y)*(b.x-t.x);}//二位叉积 
bool cmp(P a,P b){
	if(fabs(a.x-b.x)>eps)return a.x<b.x;
	return a.y<b.y;
}
P stk[110];
int sti;
int main(){
	ios::sync_with_stdio(0);
	cin.tie(0);
	N A,B,C,L;
	cin>>A.x>>A.y>>A.z>>B.x>>B.y>>B.z>>C.x>>C.y>>C.z>>L.x>>L.y>>L.z;
	N ZA=A-C,ZC=cj3(ZA,B-C),ZB=cj3(ZA,ZC);//正交基和法向量 
	double lenA=dis(ZA),lenB=dis(ZB),Ld=dot(L-A,ZC);
	//ZA和ZB的模长和AL·ZC 
	cin>>n;
	for(int i=1;i<=n;i++){
		N a;
		cin>>a.x>>a.y>>a.z;
		double t=Ld/(Ld-dot(a-A,ZC));
		N b=(N){L.x+t*(a.x-L.x),L.y+t*(a.y-L.y),L.z+t*(a.z-L.z)}-C;//求出光源与a连线与平面的交点 
		p[i]={dot(b,ZA)/lenA,dot(b,ZB)/lenB};//降至二维平面上的点 
	}
	//求凸包 
	sort(p+1,p+1+n,cmp);
	for(int i=1;i<=n;i++){
		while(sti>1&&cj(stk[sti-1],stk[sti],p[i])<=eps)sti--;
		stk[++sti]=p[i];
	}
	vector<P> v;
	for(int i=1;i<=sti;i++)v.push_back(stk[i]);
	sti=0;
	for(int i=n;i;i--){
		while(sti>1&&cj(stk[sti-1],stk[sti],p[i])<=eps)sti--;
		stk[++sti]=p[i];
	}
	for(int i=2;i<sti;i++)v.push_back(stk[i]);
	double ans=0;
	//用叉积求每个三角形的面积 
	for(int i=1;i<v.size()-1;i++)ans+=fabs(cj(v[0],v[i],v[i+1]));
	printf("%.2lf",ans/2);
	return 0;
}