计算几何----判断空间点是否在一个四面体(tetrahedron）内部

DESCRIPTION：

判断空间点 P(x, y, z)是否在一个四面体的内部？

Let the tetrahedron have vertices

        V1 = (x1, y1, z1)
        V2 = (x2, y2, z2)
        V3 = (x3, y3, z3)
        V4 = (x4, y4, z4)

and your test point be

        P = (x, y, z).

Then the point P is in the tetrahedron if following fivedeterminants all have the same sign.

             |x1 y1 z1 1|
        D0 = |x2 y2 z2 1|
             |x3 y3 z3 1|
             |x4 y4 z4 1|

             |x  y  z  1|
        D1 = |x2 y2 z2 1|
             |x3 y3 z3 1|
             |x4 y4 z4 1|

             |x1 y1 z1 1|
        D2 = |x  y  z  1|
             |x3 y3 z3 1|
             |x4 y4 z4 1|

             |x1 y1 z1 1|
        D3 = |x2 y2 z2 1|
             |x  y  z  1|
             |x4 y4 z4 1|

             |x1 y1 z1 1|
        D4 = |x2 y2 z2 1|
             |x3 y3 z3 1|
             |x  y  z  1|

简单地对上面的算法进行分析:

其实上述算法的核心思想是 四面体的体积 = 4个小四面体的之和（判断点 与 四面体的四个面各自组成的 小四面体）

但是注意： 一个四面体的体积可有上述的行列式计算， 但是行列式的值可能是负的，只有保证点的顺序是左手法则是才能保证是正的。

// copyright @ L.J.SHOU Dec.18, 2013
// test whether a point is in a tet
#include "include/cmatrix"
#include "pt.h"
#include <cassert>
#include <vector>
#include <iostream>
using namespace std;
typedef techsoft::matrix<double> Matrix;//class for matrix
typedef cpt<double> CPt; //class for points 
enum SpaceRelation{ IN, OUT, ONSURFACE};


/*
* tell whether a point is in a tetrahedran or not
* return IN, OUT, ONSURFACE
*/
SpaceRelation TestPointInTet(vector<CPt>& tet, CPt& point)
{
  assert(tet.size() == 4);
  Matrix mat[5];
  for(int i=0; i<5; ++i)
    mat[i].resize(4,4);
  double det[5];

  for(int i=0; i<4; ++i)
  {
    mat[0](i,0) = tet[i].x;
    mat[0](i,1) = tet[i].y;
    mat[0](i,2) = tet[i].z;
    mat[0](i,3) = 1;
  }

  if(mat[0].det() < 0)
  {
    swap(tet[0].x, tet[1].x);
    swap(tet[0].y, tet[1].y);
    swap(tet[0].z, tet[1].z);

    for(int i=0; i<4; ++i)
    {
      mat[0](i,0) = tet[i].x;
      mat[0](i,1) = tet[i].y;
      mat[0](i,2) = tet[i].z;
      mat[0](i,3) = 1;
    }
  }

  mat[1](0,0) = point.x;
  mat[1](0,1) = point.y;
  mat[1](0,2) = point.z;
  mat[1](0,3) = 1;
  for(int i=0; i<4; ++i)
  {
    if(i == 0) continue;
    mat[1](i,0) = tet[i].x;
    mat[1](i,1) = tet[i].y;
    mat[1](i,2) = tet[i].z;
    mat[1](i,3) = 1;
  }

  mat[2](1,0) = point.x;
  mat[2](1,1) = point.y;
  mat[2](1,2) = point.z;
  mat[2](1,3) = 1;
  for(int i=0; i<4; ++i)
  {
    if(i == 1) continue;
    mat[2](i,0) = tet[i].x;
    mat[2](i,1) = tet[i].y;
    mat[2](i,2) = tet[i].z;
    mat[2](i,3) = 1;
  }
  
  mat[3](2,0) = point.x;
  mat[3](2,1) = point.y;
  mat[3](2,2) = point.z;
  mat[3](2,3) = 1;
  for(int i=0; i<4; ++i)
  {
    if(i == 2) continue;
    mat[3](i,0) = tet[i].x;
    mat[3](i,1) = tet[i].y;
    mat[3](i,2) = tet[i].z;
    mat[3](i,3) = 1;
  }

  mat[4](3,0) = point.x;
  mat[4](3,1) = point.y;
  mat[4](3,2) = point.z;
  mat[4](3,3) = 1;
  for(int i=0; i<4; ++i)
  {
    if(i == 3) continue;
    mat[4](i,0) = tet[i].x;
    mat[4](i,1) = tet[i].y;
    mat[4](i,2) = tet[i].z;
    mat[4](i,3) = 1;
  }

  double volume = 0;
  for(int i=0; i<5; ++i)
  {
	det[i] = mat[i].det();
	//cout << det[i] << endl;
  }

  for(int i=1; i<=4; ++i)
    volume += fabs(det[i]);

  if(fabs(det[0]-volume) < 1e-15)
  {
    for(int i=1; i<=4; ++i)
	{
	  if(fabs(det[i]) < 1e-15)
	    return ONSURFACE;
    }
	return IN;
  }
  else
    return OUT;
}