求解三维空间中，点到一个有限平面的最短距离（附matlab代码）

在做项目的过程中，我们有时会遇到求解三维空间中某一个点到一个有限平面的最短距离，很多人应该第一时间会想到拉格朗日乘数法或者KKT条件，然而这些方法在求解时，列式很简单，但是求解起来较为繁琐，因此本文相处了一个较为“奇特”的思路来解决，在实际应用计算中，可以较高的提升运算速度。

在初中大家就知道，如果是无限大的平面，那么点V(x,y,z)到平面的最短距离，就是点V(x,y,z)到其在平面Ax+By+Cz+D=0上面投影点的距离，那么我们只需要考虑现在投影点是否还在问题中的有限平面上，如果在，那么就直接使用上面的无限大平面的思路即可，那么我们的问题归结为一下两点：

1.如何判断点V(x,y,z)的投影点V~是否在有限平面上；
2.当点V(x,y,z)的投影点V~不在有限平面上时，如何求最短距离；

一、针对第一个问题，我们只需要做两次内积即可，如下图：

我们只需要求出向量V’p1和向量p3p1的内积aa1，以及向量V’p1和向量p2p1的内积aa2，然后判断aa1,aa2,aa1+aa2的大小即可，如果aa1,aa2,aa1+aa2都在0~1之间，那么说明投影点V’在有限平面中，如果其中某一个不满足，那么就肯定不在。

二、针对第一个问题，我们知道此时最短距离对应的点一定在边界，因此，我们先需要求出点V(x,y,z)的投影点V’和三角面片三条边和三个顶点之间的距离，但是如何判断点到直线的距离是可以取到的呢，我引入了两个辅助变量check1和check2,那么我们只需要检查最短距时，相应的点和端点的距离（不需要求出点的坐标），判断和线段长度的关系即可，如果有一个大于，那么说明此时取不到，只能取端点值，为了简化代码，我直接把相应的距离赋值为maxx=9999999。

clear;clc;
format long
xp1=input('请输入三角面片中向量op1的x坐标：');
yp1=input('请输入三角面片中向量op1的y坐标：');
zp1=input('请输入三角面片中向量op1的z坐标：');
xp2=input('请输入三角面片中向量op2的x坐标：');
yp2=input('请输入三角面片中向量op2的y坐标：');
zp2=input('请输入三角面片中向量op2的z坐标：');
xp3=input('请输入三角面片中向量op3的x坐标：');
yp3=input('请输入三角面片中向量op3的y坐标：');
zp3=input('请输入三角面片中向量op3的z坐标：');
x=input('请输入空间中某一点的x坐标：');
y=input('请输入空间中某一点的y坐标：');
z=input('请输入空间中某一点的z坐标：');
 
x_f=(yp2-yp1)*(zp3-zp1)-(zp2-zp1)*(yp3-yp1);%法线
y_f=(zp2-zp1)*(xp3-xp1)-(xp2-xp1)*(zp3-zp1);
z_f=(xp2-xp1)*(yp3-yp1)-(yp2-yp1)*(xp3-xp1);
 
A=x_f;  %平面方程
B=y_f;
C=z_f;
D=-x_f*xp1-y_f*yp1-z_f*zp1;
t=(A*x+B*y+C*z+D)/(A^2+B^2+C^2);
 
x_new=x-A*t;%投影点坐标
y_new=y-B*t;
z_new=z-C*t;
fprintf('投影点坐标：(%f,%f,%f)\n',x_new,y_new,z_new);
 
aa1=(x_new-xp1)*(xp2-xp1)+(y_new-yp1)*(yp2-yp1)+(z_new-zp1)*(zp2-zp1);%x_f p1和p2 p1的内积
aa2=(x_new-xp1)*(xp3-xp1)+(y_new-yp1)*(yp3-yp1)+(z_new-zp1)*(zp3-zp1);%x_f p1和p3 p1的内积
if (aa1>=0)&&(aa1<=1)&&(aa2>=0)&&(aa2<=1)&&((aa1+aa2)<=1)
    min_dis=sqrt((x_new-x)^2+(y_new-y)^2+(z_new-z)^2);
    fprintf('\n投影点在三角面片内，点(%f,%f,%f)到三角面片的最短距离为：%f\n',x,y,z,min_dis);
else
    a1=xp2-xp1;
    a2=yp2-yp1;
    a3=zp2-zp1;
    b1=xp2-x_new;
    b2=yp2-y_new;
    b3=zp2-z_new;
    s1_x=a2*b3-a3*b2;
    s1_y=a3*b1-a1*b3;
    s1_z=a1*b2-a2*b1;
    min_dis1(1)=sqrt(s1_x^2+s1_y^2+s1_z^2)/sqrt(a1^2+a2^2+a3^2);
    check1=sqrt(b1^2+b2^2+b3^2-min_dis1(1)^2);
    check2=sqrt((xp1-x_new)^2+(yp1-y_new)^2+(zp1-z_new)^2-min_dis1(1)^2);
    if(check1>sqrt(a1^2+a2^2+a3^2))||(check2>sqrt(a1^2+a2^2+a3^2))
        min_dis1(1)=999999;
    end
    a1=xp3-xp1;
    a2=yp3-yp1;
    a3=zp3-zp1;
    b1=xp3-x_new;
    b2=yp3-y_new;
    b3=zp3-z_new;
    s2_x=a2*b3-a3*b2;
    s2_y=a3*b1-a1*b3;
    s2_z=a1*b2-a2*b1;
    min_dis1(2)=sqrt(s2_x^2+s2_y^2+s2_z^2)/sqrt(a1^2+a2^2+a3^2);
    check1=sqrt(b1^2+b2^2+b3^2-min_dis1(2)^2);
    check2=sqrt((xp1-x_new)^2+(yp1-y_new)^2+(zp1-z_new)^2-min_dis1(2)^2);
    if(check1>sqrt(a1^2+a2^2+a3^2))||(check2>sqrt(a1^2+a2^2+a3^2))
        min_dis1(2)=999999;
    end
    a1=xp3-xp2;
    a2=yp3-yp2;
    a3=zp3-zp2;
    b1=xp3-x_new;
    b2=yp3-y_new;
    b3=zp3-z_new;
    s3_x=a2*b3-a3*b2;
    s3_y=a3*b1-a1*b3;
    s3_z=a1*b2-a2*b1;
    min_dis1(3)=sqrt(s3_x^2+s3_y^2+s3_z^2)/sqrt(a1^2+a2^2+a3^2);
    check1=sqrt(b1^2+b2^2+b3^2-min_dis1(3)^2);
    check2=sqrt((xp2-x_new)^2+(yp2-y_new)^2+(zp2-z_new)^2-min_dis1(3)^2);
    if(check1>sqrt(a1^2+a2^2+a3^2))||(check2>sqrt(a1^2+a2^2+a3^2))
        min_dis1(3)=999999;
    end
    min_dis1(4)=sqrt((x_new-xp1)^2+(y_new-yp1)^2+(z_new-zp1)^2);
    min_dis1(5)=sqrt((x_new-xp2)^2+(y_new-yp2)^2+(z_new-zp2)^2);
    min_dis1(6)=sqrt((x_new-xp3)^2+(y_new-yp3)^2+(z_new-zp3)^2);
    min_dis=min_dis1(1);
    for i=2:6
        if min_dis1(i)<min_dis
            min_dis = min_dis1(i);
        end
    end
    min_dis=sqrt(min_dis^2+(x_new-x)^2+(y_new-y)^2+(z_new-z)^2);
    fprintf('\n投影点不在三角面片内，点(%f,%f,%f)到三角面片的最短距离为：%f\n',x,y,z,min_dis);
end