计算几何模板补充（三维空间体积、平面、直线、向量相关计算。附上hdu4741，求异面直线的最短距离与交点）

转自：http://www.cnblogs.com/wulangzhou/p/3326187.html

Link：http://acm.hdu.edu.cn/showproblem.php?pid=4741

Save Labman No.004

Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 1839    Accepted Submission(s): 603

Problem Description

Due to the preeminent research conducted by Dr. Kyouma, human beings have a breakthrough in the understanding of time and universe. According to the research, the universe in common sense is not the only one. Multi World Line is running simultaneously. In simplicity, let us use a straight line in three-dimensional coordinate system to indicate a single World Line.

During the research in World Line Alpha, the assistant of Dr. Kyouma, also the Labman No.004, Christina dies. Dr. Kyouma wants to save his assistant. Thus, he has to build a Time Tunnel to jump from World Line Alpha to World Line Beta in which Christina can be saved. More specifically, a Time Tunnel is a line connecting World Line Alpha and World Line Beta. In order to minimizing the risks, Dr. Kyouma wants you, Labman No.003 to build a Time Tunnel with shortest length.

 

Input

The first line contains an integer T, indicating the number of test cases. 

Each case contains only one line with 12 float numbers (x1, y1, z1), (x2, y2, z2), (x3, y3, z3), (x4, y4, z4), correspondingly indicating two points in World Line Alpha and World Line Beta. Note that a World Line is a three-dimensional line with infinite length. 

Data satisfy T <= 10000, |x, y, z| <= 10,000.

 

Output

For each test case, please print two lines.

The first line contains one float number, indicating the length of best Time Tunnel. 

The second line contains 6 float numbers (xa, ya, za), (xb, yb, zb), seperated by blank, correspondingly indicating the endpoints of the best Time Tunnel in World Line Alpha and World Line Beta. 

All the output float number should be round to 6 digits after decimal point. Test cases guarantee the uniqueness of the best Time Tunnel.

 

Sample Input

  
  

   
   1
1 0 1 0 1 1 0 0 0 1 1 1
  
  

 

Sample Output

  
  

   
   0.408248
0.500000 0.500000 1.000000 0.666667 0.666667 0.666667
  
  

 

Source

2013 ACM/ICPC Asia Regional Hangzhou Online

 

编程思想： 求解 空间两条异面直线的距离 问题。

利用空间向量投影的原理求解：（转自：http://blog.sina.com.cn/s/blog_648868460100nen1.html）

假设用全站仪采集空间四个点的坐标为:

    A（x1,y1,z1）,B（x2,y2,z2）,C（x3,y3,z3）,D（x4,y4,z4）

    过点A和点B的直线记为L1,过点C和点D的直线记为L2，L1和L2为空间异面直线，利用空间向量及投影原理求解L1和L2的距离d。

    首先求出直线L1和L2的方向向量：

    

    求L1和L2的公垂线方向向量：

    
   在直线L1和L2上各任取一点，这里就取点A和点C,则线段AC在公垂线方向向量上的投影就是这两个异面直线的距离。

   线段AC的方向向量为：

   

   则异面直线L1和L2的距离d为：

   

而题目不止要求距离，还有求交点，怎么求呢？通过构造平面，然后求平面与另外一条直线的交点即可得，解题思路详见如下代码。

AC code：

#include<iostream>
#include<stdio.h>
#include<cstring>
#include<algorithm>
#include<cmath>
#define eps 1e-10
#define inf 0x1f1f1f1f
using namespace std;
 
int dcmp( double a ){ if( abs(a) < eps ) return 0; if(a > 0) return 1;return -1;}
struct point{
   double x,y,z;
   point(){}
   point( double x,double y,double z ):x(x),y(y),z(z){}
};
typedef point Vector ;
point operator + ( point a,point b ){
    return point( a.x+b.x , a.y+b.y , a.z+b.z );
}
// 确定一个平面；一般式:ax+by+cz+d=0;
struct plane{
    double a,b,c,d;
    plane( double a = 0,double b = 0,double c = 0,double d = 0 ):a(a),b(b),c(c),d(d){}
};
point operator - ( point a,point b ){
    return point( a.x-b.x , a.y-b.y, a.z-b.z );
}
point operator * ( point a,double p ){
    return point( a.x*p , a.y*p , a.z*p );
}
point operator / ( point a,double p ){
    return point( a.x/p , a.y/p , a.z/p );
}
bool operator == ( point a,point b ){
    if( dcmp(a.x-b.x)  == 0&& dcmp(a.y-b.y) == 0 && dcmp(a.z-b.z) == 0 )return 1;
    return 0;
}
double Dot( point a,point b ){ return a.x*b.x + a.y*b.y + a.z*b.z; }
double Length( point a ){ return sqrt(Dot(a,a));}
double Angle(  point a,point b ){ return acos(Dot(a,b))/Length(a)/Length(b); }
//  点到  平面p0 - n 的距离；n 必须为单位向量 n 是平面法向量
double p_po_n( const point &p,const point &po,const point n ){
    return abs(Dot( p-po,n ));
}
//  点在 平面p0 - n 的投影 ； n 必须为单位向量 n 是平面法向量
point p_po_n_jec( const point &p,const point &p0,const point &n ){
    return p-n*Dot( p-p0,n );
}
//  直线p1 p2 在平面 p0 - n 的交点，假设交点存在；唯一
point p_p_p( point p1,point p2,point p0,point n ){
    point  v = p2 - p1;
    double t = ( Dot(n,p0-p1)/Dot(n,p2-p1) );
    return p1 + v*t;
}
point cross( point a,point b ){
    return point( 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 area( point a,point b,point c ){ return Length(cross(b-a,c-a)); }
double dis_to_line( point p, point a,point b ){
    point v1 = b-a, v2 = p-a;
    return Length( cross(v1,v2)/Length(v1) );
}
// p 到 ab 线段的距离
double dis_to_segm( point p,point a,point b ){
    if( a == b )return Length(p-a);
    point v1 = b-a, v2 = p-a, v3 = p-b;
        if( dcmp( Dot(v1,v2) ) < 0 ) return Length(v2);
   else if( dcmp( Dot(v1,v3) ) > 0 ) return Length(v3);
   else return Length(cross(v1,v2))/Length(v1);
}
// 四面体的体积
double dis_l_to_l( point a,point b,point lt,point rt )
{
    while( abs(rt.x - lt.x) > eps || abs(rt.y - lt.y) > eps || abs( rt.z - lt.z) > eps )
    {
        point temp; temp = lt + ( rt - lt )/3.0;
        point now;   now = lt + ( rt - lt )/3.0*2.0;
        if( dis_to_segm(temp,a,b) < dis_to_segm(now,a,b) )
             rt = now;
        else lt = temp;
    }
    return dis_to_segm( rt,a,b );
}
// 获取平面 a*x + b*y + c*z + d; 由两个向量和一个平面上一个点确定一个平面
plane get_plane( Vector a,Vector b,point c )
{
    Vector p = cross( a,b ); // 由平面两个向量a,b,  可以确定平面法向量；
    double d = ( p.x*c.x + p.y*c.y + p.z*c.z)*-1;
    return plane( p.x,p.y,p.z,d );
}
// 获取 直线cd与平面p的交点；
point get_point( plane p,point c,point d )
{
    Vector v = c - d;
    double t =  p.a*c.x+p.b*c.y + p.c*c.z + p.d ;
    double s =  p.a*(d.x - c.x) + p.b*(d.y-c.y) + p.c*(d.z - c.z);
    double k = t/s;
    return c + v*k;
}
 
point A[5];
int main( )
{
   int T; scanf("%d",&T);
   while( T-- )
   {
       for( int i = 1; i <= 4; i++ )
         scanf("%lf %lf %lf",&A[i].x,&A[i].y,&A[i].z);
       Vector v1 = A[2] - A[1];
       Vector v2 = A[4] - A[3];
       Vector v3 = cross( v1,v2 ); // 获取公垂线；
        plane a1 = get_plane( v1,v3,A[1] ); // 两个向量和一个点确定一个平面；
        plane a2 = get_plane( v2,v3,A[3] );
       point p1 = get_point( a1,A[4],A[3] ); // 获取交点；
       point p2 = get_point( a2,A[1],A[2] );
       printf("%.6lf\n",Length(p1-p2));
       printf("%.6lf %.6lf %.6lf %.6lf %.6lf %.6lf\n",p2.x,p2.y,p2.z,p1.x,p1.y,p1.z);
   }
   return 0;
}

附上向量叉积公式：

a(x1,y1,z1),b(x2,y2,z2);

            |  i    j   k   |

a*b=    |x1  y1  z1|  ={(y1*z2-z1*y2)i, (x2*z1-z2*x1)j, (x1*y2-y1*x2)k}

           |x2  y2  z2|