需要提前知道的知识

点积

其中 $\theta$ 是两个向量之间的夹角。

叉积

其中 $\theta$ 是两个向量之间的夹角，$\vec{n}$ 是 $\vec{a}$ 和 $\vec{b}$ 所在平面的法向量。

那么可以用行列式表达为：

这个行列式可以展开为：

直线的方程

其中 $\vec{r_0}$ 是直线上的一个点，$\vec{v}$ 是直线的方向向量，$t$ 是参数。

如果 $\vec{v}$ 的坐标知道，那么我们可以写成参数方程的形式：

而其写成 $t$ 的表达式为：

平面方程

其中 $\vec{n}$ 是平面的法向量，$\vec{r_0}$ 是平面上的一个点。

写成方程的形式为：

其中 $D = -\vec{n} \cdot \vec{r_0}$ ，即把 $\vec{r_0}$ 代入平面方程中即可，或者把平面中的一点坐标带入。

点到平面的距离

其中 $n_x,n_y,n_z$ 是平面的法向量，$x_0,y_0,z_0$ 是平面上的一个点，$D$ 是平面方程中的常数项。

题目

描述

已知三角形平面 $ABC$ 的三个顶点坐标和另外一点 $P$ 的坐标，求 $P$ 到三角形平面$ABC$ 的距离。

这里三角形平面 $ABC$ 是一个薄片

题目分析

要求 $P$ 到三角形平面 $ABC$ 的距离，我们可以先想想什么是点 $P$ 到平面 $ABC$ 的距离？

  1. 如果 $P$ 到三角形平面的投影点 $D$ 在 $ABC$ 内，那么 $PD$ 就是 $P$ 到三角形平面的距离（如图一）。
    图一
  2. 如果 $P$ 到三角形平面的投影点 $D$ 在 $ABC$ 外，那么点 $P$ 到 $AB$ 或 $BC$ 或 $CA$ 的距离的最小值就是 $P$ 到三角形平面的距离。但注意到 $AB$ 这些为线段，过 $P$ 的垂线与线段相交的点不一定在 $AB$ 上，所以又有以下两种情况：
    1. $P$ 到 $AB$ 的垂线与 $AB$ 相交于点 $E$，那么 $PE$ 就是 $P$ 到三角形平面的距离（如图二）。
      图二
    2. $P$ 到 $AB$ 的垂线与 $AB$ 相交于点 $E$，但 $E$ 不在 $AB$ 上，那么 $P$ 到这条线段的距离就是 $P$ 到 $A$ 或 $B$ 的距离的最小值。
      图三

求解过程

好了，现在我们知道大概的过程，接下来就是数学推导了。

1. 求 $P$ 到三角形平面 $ABC$ 的垂足

设 $P$ 到三角形平面 $ABC$ 的垂足为 $D$，那么 $PD$ 就是 $P$ 到三角形平面的距离。

根据平面方程，我们可以得到，三角形平面 $ABC$ 的方程为：

其中

作直线 $PD$ ，其中 $D$ 为 $P$ 到三角形平面 $ABC$ 的垂足，其方程为：

带入平面方程解得，那么直线 $PD$ 的参数方程为：

那么 $P$ 到三角形平面 $ABC$ 的垂足 $D$ 的坐标为：

我们这里把 $D=-(n_xx_A+n_yy_A+n_zz_A)$ 带入上式有:

2. 判断 $D$ 是否在 $ABC$ 内

判断 $D$ 是否在 $ABC$ 内，我们可以用向量叉积来判断。

我们知道，如果 $D$ 在 $ABC$ 内，那么三角形 $ABC$ 的面积等于三角形 $ABD$ ，$ACD$ ，$BCD$ 的面积之和。又$S_{ABC}=|AB|\cdot|AC|\cdot \sin(\theta)=|\vec{AB} \times \vec{AC}|$，那么如果$|\vec{DA} \times \vec{DB}|+|\vec{DB} \times \vec{DC}|+|\vec{DC} \times \vec{DA}|=|\vec{AB} \times \vec{AC}|$，那么$D$就在$ABC$内。

3. 求距离

1. $D$ 在 $ABC$ 内

图一

如果 $D$ 在 $ABC$ 内，那么 $PD$ 就是 $P$ 到三角形平面的距离，即：

这里由 $D=-(n_xx_A+n_yy_A+n_zz_A)$ 带入有：

2. $D$ 在 $ABC$ 外

我们先求 $P$ 到 $AB$ 的距离，再求 $P$ 到 $BC$ 的距离，再求 $P$ 到 $CA$ 的距离，取最小值即可。

$P$ 到 $AB$ 的距离

图三

我们先求 $P$ 到 $AB$ 的垂足 $E$ 的坐标。

设过 $P$ 点且垂直于 $AB$ 的平面为 $\pi$ ，那么 $\pi$ 的方程为：

那么我们联立直线 $AB$ 的方程：$\frac{x-x_A}{x_B-x_A} = \frac{y-y_A}{y_B-y_A} = \frac{z-z_A}{z_B-z_A}=t$ 和平面 $\pi$ 的方程，解得 $t$ 为：

其中我们限制 $t \in [0,1] $ ，如果其在外面，将其限制在靠近 $[0,1]$ 的边界。
这样我们就得到了 $P$ 到 $AB$ 的垂足 $E$ 的坐标为：

那么 $P$ 到 $AB$ 的距离为：

剩下的同理，这里不过多说明，值得注意的是，我们可以构造这样的函数来方便求三次距离。

代码实现

方法一

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
#include<bits/stdc++.h>
using namespace std;
#define dl double

dl dot(dl x1,dl y1,dl z1,dl x2,dl y2,dl z2)
{
    return x1*x2+y1*y2+z1*z2;
}

dl abscross(dl x1,dl y1,dl z1,dl x2,dl y2,dl z2)
{
    dl x=y1*z2-z1*y2,y=z1*x2-x1*z2,z=x1*y2-y1*x2;
    return sqrt(x*x+y*y+z*z);
}

dl dist(dl px,dl py,dl pz,dl x1,dl y1,dl z1,dl x2,dl y2,dl z2)
{
    dl dx=x2-x1,dy=y2-y1,dz=z2-z1;
    dl t=(dot(px-x1,py-y1,pz-z1,dx,dy,dz))/(dx*dx+dy*dy+dz*dz);
    t=max((dl)0.0,min((dl)1.0,t));
    dl xp=x1+t*dx,yp=y1+t*dy,zp=z1+t*dz;
    return sqrt((px-xp)*(px-xp)+(py-yp)*(py-yp)+(pz-zp)*(pz-zp));
}

int main()
{
    int t;cin>>t;
    while(t--)
    {
        dl px,py,pz,ax,ay,az,bx,by,bz,cx,cy,cz;
        cin>>px>>py>>pz>>ax>>ay>>az>>bx>>by>>bz>>cx>>cy>>cz;
        dl abx=bx-ax,aby=by-ay,abz=bz-az;
        dl acx=cx-ax,acy=cy-ay,acz=cz-az;
        dl nx=aby*acz-abz*acy,ny=abz*acx-abx*acz,nz=abx*acy-aby*acx;
        dl norm=sqrt(nx*nx+ny*ny+nz*nz);
        dl D=dot(nx,ny,nz,px-ax,py-ay,pz-az);
        dl h=fabs(D)/norm;
        dl tt=-D/(nx*nx+ny*ny+nz*nz);
        dl ox=px+nx*tt,oy=py+ny*tt,oz=pz+nz*tt;
        dl oax=ax-ox,oay=ay-oy,oaz=az-oz;
        dl obx=bx-ox,oby=by-oy,obz=bz-oz;
        dl ocx=cx-ox,ocy=cy-oy,ocz=cz-oz;
        dl tot=abscross(abx,aby,abz,acx,acy,acz);
        dl a1=abscross(oax,oay,oaz,obx,oby,obz);
        dl a2=abscross(obx,oby,obz,ocx,ocy,ocz);
        dl a3=abscross(ocx,ocy,ocz,oax,oay,oaz);
        if(fabs(a1+a2+a3-tot)<1e-9)printf("%.3lf\n",h);
        else
        {
            dl d1=dist(px,py,pz,ax,ay,az,bx,by,bz);
            dl d2=dist(px,py,pz,bx,by,bz,cx,cy,cz);
            dl d3=dist(px,py,pz,cx,cy,cz,ax,ay,az);
            printf("%.3lf\n",min({d1,d2,d3}));
        }
    }
}

这里说明提示一下： 1. 由于 $\frac{n_x(x_p-x_A)+n_y(y_p-y_A)+n_z(z_p-z_A)}{n_x^2 + n_y^2 + n_z^2}$ 在此后的运算中至少要用到3次，那么不妨用这个作为 $D$ 的值。 2. 在计算点到线段的距离时，由于我们要限制 $t \in [0,1]$ ，那么我们可以用C++中的max和min函数来限制 $t$ 的范围，先用min限制$t \leq 1$，然后再用max限制 $t \geq 0$ 。

方法二

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
#include<bits/stdc++.h>
using namespace std;
#define dl double

dl dot(dl x1,dl y1,dl z1,dl x2,dl y2,dl z2)
{
    return x1*x2+y1*y2+z1*z2;
}

dl abscross(dl x1,dl y1,dl z1,dl x2,dl y2,dl z2)
{
    dl x=y1*z2-z1*y2,y=z1*x2-x1*z2,z=x1*y2-y1*x2;
    return sqrt(x*x+y*y+z*z);
}

dl dist(dl px,dl py,dl pz,dl x1,dl y1,dl z1,dl x2,dl y2,dl z2)
{
    dl dx=x2-x1,dy=y2-y1,dz=z2-z1;
    dl t=dot(px-x1,py-y1,pz-z1,dx,dy,dz)/(dx*dx+dy*dy+dz*dz);
    t=max((dl)0.0,min((dl)1.0,t));
    dl xp=x1+t*dx,yp=y1+t*dy,zp=z1+t*dz;
    return sqrt((px-xp)*(px-xp)+(py-yp)*(py-yp)+(pz-zp)*(pz-zp));
}

int main()
{
    int t;cin>>t;
    while(t--)
    {
        dl px,py,pz,ax,ay,az,bx,by,bz,cx,cy,cz;
        cin>>px>>py>>pz>>ax>>ay>>az>>bx>>by>>bz>>cx>>cy>>cz;
        
        // 平移坐标系使P为原点
        ax-=px; ay-=py; az-=pz;
        bx-=px; by-=py; bz-=pz;
        cx-=px; cy-=py; cz-=pz;
        px=py=pz=0;  // P点现在坐标为(0,0,0)

        dl abx=bx-ax,aby=by-ay,abz=bz-az;
        dl acx=cx-ax,acy=cy-ay,acz=cz-az;
        dl nx=aby*acz-abz*acy;
        dl ny=abz*acx-abx*acz;
        dl nz=abx*acy-aby*acx;

        // 计算投影点O
        dl D=dot(nx,ny,nz,ax,ay,az); 
        dl tt=-D/(nx*nx+ny*ny+nz*nz);
        dl ox=nx*tt, oy=ny*tt, oz=nz*tt;  

        // 计算相对向量
        dl oax=ax-ox, oay=ay-oy, oaz=az-oz;
        dl obx=bx-ox, oby=by-oy, obz=bz-oz;
        dl ocx=cx-ox, ocy=cy-oy, ocz=cz-oz;

        // 面积判断
        dl tot=abscross(abx,aby,abz,acx,acy,acz);
        dl a1=abscross(oax,oay,oaz,obx,oby,obz);
        dl a2=abscross(obx,oby,obz,ocx,ocy,ocz);
        dl a3=abscross(ocx,ocy,ocz,oax,oay,oaz);
        
        dl h=fabs(D)/sqrt(nx*nx+ny*ny+nz*nz);
        if(fabs(a1+a2+a3-tot)<1e-9)
            printf("%.3lf\n",h);
        else 
        {
            dl d1=dist(0,0,0,ax,ay,az,bx,by,bz);
            dl d2=dist(0,0,0,bx,by,bz,cx,cy,cz);
            dl d3=dist(0,0,0,cx,cy,cz,ax,ay,az);
            printf("%.3lf\n",min({d1,d2,d3}));
        }
    }
}

方法二和方法一没有本质的区别，只是把坐标系平移了一下，这样在计算投影点时，可以少计算一些量。