【三维几何】ICPC Dhaka Regional 2018 - I - Triangles

题目链接https://codeforces.com/gym/250668

---

题意

给出三维空间中两个三角形平面，问这两个三角形最近点对距离。

---

题解

枚举三角形每条线段到另一个三角形的距离。我是三分写法，三分出线段上的一个点，计算出点到三角形平面的距离。

点到三角形平面的距离分两种情况：

1. 点到平面的投影在三角形内，计算点到投影的距离
2. 点到平面的投影在三角形外，计算点到三条线段的最短距离

点到线段距离也分两种情况：

1. 点到线段的投影在线段内，计算点到线段的距离
2. 点到线段的投影在线段外，计算点到两个端点的最短距离

应该是有不三分的方法，待补。

---

#include<bits/stdc++.h>
using namespace std;
typedef long double db;
const int N=1e4+7;
const db eps=1e-8;
const db pi=acos(-1);
const db inf=1e18;
int sign(db k){if(k>eps) return 1;if(k<-eps) return -1;return 0;}
int dcmp(db k1,db k2){return sign(k1-k2);}
int inmid(db k1,db k2,db k3){return sign(k1-k2)*sign(k3-k2)<=0;}
struct Point{
    db x,y,z;
    Point operator - (const Point &k)const{
        return (Point){x-k.x,y-k.y,z-k.z};
    }
    Point operator + (const Point &k)const{
        return (Point){x+k.x,y+k.y,z+k.z};
    }
    Point operator * (const db k)const{
        return (Point){x*k,y*k,z*k};
    }
    Point operator / (const db k)const{
        return (Point){x/k,y/k,z/k};
    }
    db operator * (const Point &k)const{
        return x*k.x+y*k.y+z*k.z;
    }
    bool operator == (const Point &k)const{
        return (dcmp(x,k.x)==0)&&(dcmp(y,k.y)==0)&&(dcmp(z,k.z)==0);
    }
    Point operator ^ (const Point &k)const{
        return (Point){y*k.z-z*k.y,z*k.x-x*k.z,x*k.y-y*k.x};
    }
    db abs(){return sqrt(x*x+y*y+z*z);}
    db abs2(){return x*x+y*y+z*z;}
    Point unit(){return (*this)/abs();}
    db dis(Point k){return ((*this)-k).abs();}
    db angle(Point k){
        return acos((*this)*k/abs()/k.abs());
    }
    void input(){
        double xx,yy,zz;
        scanf("%lf%lf%lf",&xx,&yy,&zz);
        x=xx;y=yy;z=zz;
    }
    void output(){printf("(%f,%f,%f)\n",(double)x,(double)y,(double)z);}
};
bool inmid(Point k1,Point k2,Point k3){
    return inmid(k1.x,k2.x,k3.x)&&
        inmid(k1.y,k2.y,k3.y)&&
        inmid(k1.z,k2.z,k3.z);
}

struct Line{
    Point s,e;
    Point proj(Point q){
        return s+(((e-s)*((e-s)*(q-s)))/((e-s).abs2()));
    }
    db disSP(Point q){
        if(inmid(s,proj(q),e)) return proj(q).dis(q);
        return min(q.dis(s),q.dis(e));
    }
};
struct Plane{
    Point a,b,c,o;
    Point pvec(){
        return (b-a)^(c-a);
    }
    int crossLine(Line u,Point &p){
        db x=o*(u.e-a);
        db y=o*(u.s-a);
        db d=x-y;
        if(sign(d)==0) return 0;
        p=((u.s*x)-(u.e*y))/d;
        return 1;
    }
    Point pointToPlane(const Point p){
        Point tmp=p;
        Line u=(Line){tmp,tmp+o};
        crossLine(u,tmp);
        return tmp;
    }
    bool inPoly(const Point p){
        db alpha=0;
        if(p==a||p==b||p==c) return true;
        alpha+=((p-a)^(p-b)).abs();
        alpha+=((p-a)^(p-c)).abs();
        alpha+=((p-b)^(p-c)).abs();
        db area=0;
        area+=((a-b)^(a-c)).abs();
        return dcmp(alpha,area)==0;
    }
    db disSP(const Point k){
        Point p=pointToPlane(k);
        db res=inf;
        if(inPoly(p)){
            res=min(res,(p-k).abs());
        }
        else{
            res=min(res,(Line){a,b}.disSP(k));
            res=min(res,(Line){a,c}.disSP(k));
            res=min(res,(Line){b,c}.disSP(k));
        }
        return res;
    }
    db disSL(Line k){
        Point lo=k.s,hi=k.e;
        for(int i=1;i<=90;i++){
            Point lm=lo+(hi-lo)/3;
            Point rm=hi-(hi-lo)/3;
            if(dcmp(disSP(lm),disSP(rm))<=0) hi=rm;
            else lo=lm;
        }
        return disSP(lo);
    }
}s1,s2;

int main()
{
    //freopen("in.in","r",stdin);
    int T;
    scanf("%d",&T);
    while(T--){
        s2.a.input();s2.b.input();s2.c.input();s2.o=s2.pvec();
        s1.a.input();s1.b.input();s1.c.input();s1.o=s1.pvec();
        db ans=inf;
        ans=min(ans,s1.disSL((Line){s2.a,s2.b}));
        ans=min(ans,s1.disSL((Line){s2.a,s2.c}));
        ans=min(ans,s1.disSL((Line){s2.b,s2.c}));

        ans=min(ans,s2.disSL((Line){s1.a,s1.b}));
        ans=min(ans,s2.disSL((Line){s1.a,s1.c}));
        ans=min(ans,s2.disSL((Line){s1.b,s1.c}));
        cout<<fixed<<setprecision(10)<<ans<<endl;
        //printf("%.10f\n",ans);
    }
    return 0;
}