判断空间点是否在一个四面体(tetrahedron）内部

计算几何----判断空间点是否在一个四面体(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个小四面体的之和（判断点 与 四面体的四个面各自组成的 小四面体）

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

1. // copyright @ L.J.SHOU Dec.18, 2013  
2. // test whether a point is in a tet  
3. #include "include/cmatrix"  
4. #include "pt.h"  
5. #include <cassert>  
6. #include <vector>  
7. #include <iostream>  
8. using namespace std;  
9. typedef techsoft::matrix<double> Matrix;//class for matrix  
10. typedef cpt<double> CPt; //class for points   
11. enum SpaceRelation{ IN, OUT, ONSURFACE};  
12.   
13.   
14. /* 
15. * tell whether a point is in a tetrahedran or not 
16. * return IN, OUT, ONSURFACE 
17. */  
18. SpaceRelation TestPointInTet(vector<CPt>& tet, CPt& point)  
19. {  
20.   assert(tet.size() == 4);  
21.   Matrix mat[5];  
22.   for(int i=0; i<5; ++i)  
23.     mat[i].resize(4,4);  
24.   double det[5];  
25.   
26.   for(int i=0; i<4; ++i)  
27.   {  
28.     mat[0](i,0) = tet[i].x;  
29.     mat[0](i,1) = tet[i].y;  
30.     mat[0](i,2) = tet[i].z;  
31.     mat[0](i,3) = 1;  
32.   }  
33.   
34.   if(mat[0].det() < 0)  
35.   {  
36.     swap(tet[0].x, tet[1].x);  
37.     swap(tet[0].y, tet[1].y);  
38.     swap(tet[0].z, tet[1].z);  
39.   
40.     for(int i=0; i<4; ++i)  
41.     {  
42.       mat[0](i,0) = tet[i].x;  
43.       mat[0](i,1) = tet[i].y;  
44.       mat[0](i,2) = tet[i].z;  
45.       mat[0](i,3) = 1;  
46.     }  
47.   }  
48.   
49.   mat[1](0,0) = point.x;  
50.   mat[1](0,1) = point.y;  
51.   mat[1](0,2) = point.z;  
52.   mat[1](0,3) = 1;  
53.   for(int i=0; i<4; ++i)  
54.   {  
55.     if(i == 0) continue;  
56.     mat[1](i,0) = tet[i].x;  
57.     mat[1](i,1) = tet[i].y;  
58.     mat[1](i,2) = tet[i].z;  
59.     mat[1](i,3) = 1;  
60.   }  
61.   
62.   mat[2](1,0) = point.x;  
63.   mat[2](1,1) = point.y;  
64.   mat[2](1,2) = point.z;  
65.   mat[2](1,3) = 1;  
66.   for(int i=0; i<4; ++i)  
67.   {  
68.     if(i == 1) continue;  
69.     mat[2](i,0) = tet[i].x;  
70.     mat[2](i,1) = tet[i].y;  
71.     mat[2](i,2) = tet[i].z;  
72.     mat[2](i,3) = 1;  
73.   }  
74.     
75.   mat[3](2,0) = point.x;  
76.   mat[3](2,1) = point.y;  
77.   mat[3](2,2) = point.z;  
78.   mat[3](2,3) = 1;  
79.   for(int i=0; i<4; ++i)  
80.   {  
81.     if(i == 2) continue;  
82.     mat[3](i,0) = tet[i].x;  
83.     mat[3](i,1) = tet[i].y;  
84.     mat[3](i,2) = tet[i].z;  
85.     mat[3](i,3) = 1;  
86.   }  
87.   
88.   mat[4](3,0) = point.x;  
89.   mat[4](3,1) = point.y;  
90.   mat[4](3,2) = point.z;  
91.   mat[4](3,3) = 1;  
92.   for(int i=0; i<4; ++i)  
93.   {  
94.     if(i == 3) continue;  
95.     mat[4](i,0) = tet[i].x;  
96.     mat[4](i,1) = tet[i].y;  
97.     mat[4](i,2) = tet[i].z;  
98.     mat[4](i,3) = 1;  
99.   }  
100.   
101.   double volume = 0;  
102.   for(int i=0; i<5; ++i)  
103.   {  
104.     det[i] = mat[i].det();  
105.     //cout << det[i] << endl;  
106.   }  
107.   
108.   for(int i=1; i<=4; ++i)  
109.     volume += fabs(det[i]);  
110.   
111.   if(fabs(det[0]-volume) < 1e-15)  
112.   {  
113.     for(int i=1; i<=4; ++i)  
114.     {  
115.       if(fabs(det[i]) < 1e-15)  
116.         return ONSURFACE;  
117.     }  
118.     return IN;  
119.   }  
120.   else  
121.     return OUT;  
122. }  

转载：http://www.cnblogs.com/shoulinjun/p/3815627.html