快速检测空间三角形相交算法的代码实现(Devillers & Guigue算法)

  Devillers & Guigue算法( 简称Devillers 算法) 通过三角形各顶点构成的行列式正负的几何意义来判断三角形中点、线、面之间的相对位置关系,从而判断两三角形是否相交。其基本原理如下:给定空间四个点：a(ax,ay,az),b=(bx,by,bz),c=(cx,cy,cz),d=(dx,dy,dz),定义行列式如下：

    [ a, b, c, d] 采用右手螺旋法则定义了四个空间点的位置关系。[ a, b, c, d] > 0 表示 d 在 a、b、c 按逆时针顺序所组成的三角形的正法线方向( 即上方) ; [ a, b, c, d] < 0 表示 d 在 △abc的下方; [ a, b, c, d] = 0 表示四点共面。

    设两个三角形T1和T2，顶点分别为：V10，V11，V12和V20，V21，V22，三角形所在的平面分别为π1和π2，其法向量分别为N1和N2.算法先判别三角形和另一个三角形所在的平面的相互位置关系,提前排除不相交的情况。通过计算[V20,V21,V22,V1i].(i=0,1,2)来判断T1和π2的关系：如果所有的行列式的值都不为零且同号，则T1和T2不相交；否则T1和π2相交。相交又分为如下几种情况：

     a)如果所有的行列式的值为零，则T1和T2共面，转化为共面的线段相交问题。

     b)如果其中一个行列式的值为零，而其他两个行列式同号，则只有一个点在平面内，测试顶点是否则T2内部，是则相交，否则不相交；

     c)否则T1的顶点位于平面π2两侧(包含T1的一条边在平面π2中的情况)。

     再按照类似的方法对 T 2 和 π 1 作进一步的测试。如果通过测试, 则每个三角形必有确定的一点位于另一个三角形所在平面的一侧, 而另外两点位于其另一侧。算法分别循环置换每个三角形的顶点, 以使V10(V20)位于π2(π1)的一侧，另两个点位于其另一侧；同时对顶点V21，V22(V11,V12)进行交换操作，以确保V10(V20)位于π2(π1)的上方，即正法线方向。经过以上的预排除和置换操作，V10的邻边V10V11，V10V12和V20的邻边V20V21和V20V22与两平面的交线L相交于固定形式的点上，分别记为i，j，k，l(i<j,k<l),如图：

    这些点在L上形成的封闭区间为i1=[i,j],i2=[k,l].至此，两个三角形的相交测试问题转换为封闭区间i1，i2的重叠问题。若重叠则相交，否则不相交。由于交点形式固定，只需满足条件k<=j且i<=l即表明区间重叠，条件还可进一步缩减为判别式(1)是否成立：

    [V10,V11,V20,V21]<=0 && [V10,V12,V22,V20]<=0                                                    (1)

以上资料来自《计算机应用研究》第25卷第10期，邹益胜等的《空间三角形快速相交检测算法》下面是具体代码实现，代码测试通过并且已经用于公司产品。

 

[cpp]  view plain copy print ?

1. typedef float float3[3];  
2.   
3. enum TopologicalStructure   
4. {   
5.     INTERSECT, NONINTERSECT   
6. };  
7.   
8. struct Triangle  
9. {  
10.     float3 Normal_0;   
11.     float3 Vertex_1, Vertex_2, Vertex_3;  
12. };  
13.   
14. struct point  
15. {  
16.     float x, y;  
17. };  
18.   
19.   
20. //四点行列式  
21. inline float get_vector4_det( float3 v1, float3 v2, float3 v3, float3 v4 )  
22. {  
23.     float a[3][3];  
24.     for ( int i = 0; i != 3; ++i )  
25.     {  
26.         a[0][i] = v1[i] - v4[i];  
27.         a[1][i] = v2[i] - v4[i];  
28.         a[2][i] = v3[i] - v4[i];  
29.     }  
30.   
31.     return a[0][0] * a[1][1] * a[2][2]   
32.         + a[0][1] * a[1][2] * a[2][0]   
33.         + a[0][2] * a[1][0] * a[2][1]   
34.         - a[0][2] * a[1][1] * a[2][0]   
35.         - a[0][1] * a[1][0] * a[2][2]   
36.         - a[0][0] * a[1][2] * a[2][1];  
37. }  
38.   
39.   
40. //利用叉积计算点p相对线段p1p2的方位  
41. inline double direction( point p1, point p2, point p ){  
42.     return ( p.x - p1.x ) * ( p2.y - p1.y ) - ( p2.x - p1.x ) * ( p.y - p1.y )   ;  
43. }  
44.   
45.   
46. //确定与线段p1p2共线的点p是否在线段p1p2上  
47. inline int on_segment( point p1, point p2, point p ){  
48.     double max = p1.x > p2.x ? p1.x : p2.x ;  
49.     double min = p1.x < p2.x ? p1.x : p2.x ;  
50.     double max1 = p1.y > p2.y ? p1.y : p2.y ;  
51.     double min1 = p1.y < p2.y ? p1.y : p2.y ;  
52.     if ( p.x >= min && p.x <= max && p.y >= min1 && p.y <= max1 )  
53.     {  
54.         return 1;  
55.     }  
56.     else  
57.     {  
58.         return 0;  
59.     }  
60. }  
61.   
62.   
63. //判断线段p1p2与线段p3p4是否相交的主函数  
64. inline int segments_intersert( point p1, point p2, point p3, point p4 ){  
65.     double d1, d2, d3, d4;  
66.     d1 = direction( p3, p4, p1 );  
67.     d2 = direction( p3, p4, p2 );  
68.     d3 = direction( p1, p2, p3 );  
69.     d4 = direction( p1, p2, p4 );  
70.     if ( d1 * d2 < 0 && d3 * d4 < 0 )  
71.     {  
72.         return 1;  
73.     }  
74.     else if ( d1 == 0 && on_segment( p3, p4, p1 ) == 1 )  
75.     {  
76.         return 1;  
77.     }  
78.     else if ( d2 == 0 && on_segment( p3, p4, p2 ) == 1 )  
79.     {  
80.         return 1;  
81.     }  
82.     else if ( d3 == 0 && on_segment( p1, p2, p3 ) == 1 )  
83.     {  
84.         return 1;  
85.     }  
86.     else if ( d4 == 0 && on_segment( p1, p2, p4 ) == 1 )  
87.     {  
88.         return 1;  
89.     }  
90.     return 0;  
91. }  
92.   
93.   
94. //判断同一平面的直线和三角形是否相交  
95. inline bool line_triangle_intersert_inSamePlane( Triangle* tri, float3 f1, float3 f2 )  
96. {  
97.     point p1, p2, p3, p4;  
98.   
99.     copy_point( p1, f1 );  
100.   
101.     copy_point( p2, f2 );  
102.   
103.     copy_point( p3, tri->Vertex_1 );  
104.   
105.     copy_point( p4, tri->Vertex_2 );  
106.   
107.     if ( segments_intersert( p1, p2, p3, p4 ) )  
108.     {  
109.         return true;  
110.     }  
111.   
112.     copy_point( p3, tri->Vertex_2 );  
113.   
114.     copy_point( p4, tri->Vertex_3 );  
115.   
116.     if ( segments_intersert( p1, p2, p3, p4 ) )  
117.     {  
118.         return true;  
119.     }  
120.   
121.     copy_point( p3, tri->Vertex_1 );  
122.   
123.     copy_point( p4, tri->Vertex_3 );  
124.   
125.     if ( segments_intersert( p1, p2, p3, p4 ) )  
126.     {  
127.         return true;  
128.     }  
129.   
130.     return false;  
131. }   
132.   
133.   
134. //判断同一平面内的三角形是否相交  
135. inline bool triangle_intersert_inSamePlane( Triangle* tri1, Triangle* tri2 )  
136. {  
137.     if ( line_triangle_intersert_inSamePlane( tri2, tri1->Vertex_1, tri1->Vertex_2 ) )  
138.     {  
139.         return true;  
140.     }  
141.     else if ( line_triangle_intersert_inSamePlane( tri2, tri1->Vertex_2, tri1->Vertex_3 ) )  
142.     {  
143.         return true;  
144.     }  
145.     else if ( line_triangle_intersert_inSamePlane( tri2, tri1->Vertex_1, tri1->Vertex_3 ) )  
146.     {  
147.         return true;  
148.     }  
149.     else  
150.     {  
151.         return false;  
152.     }  
153. }   
154.   
155. //向量之差  
156. inline void get_vector_diff( float3& aimV, const float3 a, const float3 b )  
157. {  
158.     aimV[0] = b[0] - a[0];  
159.   
160.     aimV[1] = b[1] - a[1];  
161.   
162.     aimV[2] = b[2] - a[2];  
163. }   
164.   
165. //向量内积  
166. inline float Dot( const float3& v1, const float3& v2 )  
167. {      
168.     return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2] ;  
169. }   
170.   
171. //重心法判断点是否在三角形内部  
172. inline bool is_point_within_triangle( Triangle* tri, float3 point )  
173. {  
174.     float3 v0;     
175.     get_vector_diff( v0, tri->Vertex_1, tri->Vertex_3 );  
176.     float3 v1;     
177.     get_vector_diff( v1, tri->Vertex_1, tri->Vertex_2 );  
178.     float3 v2;   
179.     get_vector_diff( v2, tri->Vertex_1, point );  
180.     float dot00 = Dot( v0, v0 ) ;     
181.     float dot01 = Dot( v0, v1 ) ;      
182.     float dot02 = Dot( v0, v2 ) ;      
183.     float dot11 = Dot( v1, v1 ) ;     
184.     float dot12 = Dot( v1, v2 ) ;      
185.     float inverDeno = 1 / ( dot00* dot11 - dot01* dot01 ) ;      
186.     float u = ( dot11* dot02 - dot01* dot12 ) * inverDeno ;      
187.     if ( u < 0 || u > 1 ) // if u out of range, return directly  
188.     {  
189.         return false ;  
190.     }      
191.     float v = ( dot00* dot12 - dot01* dot02 ) * inverDeno ;      
192.     if ( v < 0 || v > 1 ) // if v out of range, return directly  
193.     {  
194.         return false ;  
195.     }      
196.     return u + v <= 1 ;  
197. }    
198.   
199. //Devillers算法主函数  
200. inline TopologicalStructure judge_triangle_topologicalStructure( Triangle* tri1, Triangle* tri2 )  
201. {  
202.     //设tri1所在的平面为p1,tri2所在的平面为p2  
203.     float p1_tri2_vertex1 = get_vector4_det( tri1->Vertex_1, tri1->Vertex_2, tri1->Vertex_3, tri2->Vertex_1 );  
204.   
205.     float p1_tri2_vertex2 = get_vector4_det( tri1->Vertex_1, tri1->Vertex_2, tri1->Vertex_3, tri2->Vertex_2 );  
206.   
207.     float p1_tri2_vertex3 = get_vector4_det( tri1->Vertex_1, tri1->Vertex_2, tri1->Vertex_3, tri2->Vertex_3 );  
208.   
209.   
210.     if ( p1_tri2_vertex1 > 0 && p1_tri2_vertex2 > 0 && p1_tri2_vertex3 > 0 )  
211.     {  
212.         return NONINTERSECT;  
213.     }  
214.   
215.     if ( p1_tri2_vertex1 < 0 && p1_tri2_vertex2 < 0 && p1_tri2_vertex3 < 0 )  
216.     {  
217.         return NONINTERSECT;  
218.     }  
219.   
220.   
221.     if ( p1_tri2_vertex1 == 0 && p1_tri2_vertex2 == 0 && p1_tri2_vertex3 == 0 )  
222.     {  
223.         if ( triangle_intersert_inSamePlane( tri1, tri2 ) )  
224.         {  
225.             return INTERSECT;  
226.         }  
227.         else  
228.         {  
229.             return NONINTERSECT;  
230.         }  
231.     }  
232.   
233.   
234.     if ( p1_tri2_vertex1 == 0 && p1_tri2_vertex2 * p1_tri2_vertex3 > 0 )  
235.     {  
236.         if ( is_point_within_triangle( tri1, tri2->Vertex_1 ) )  
237.         {  
238.             return INTERSECT;  
239.         }  
240.         else  
241.         {  
242.             return NONINTERSECT;  
243.         }  
244.     }  
245.     else if ( p1_tri2_vertex2 == 0 && p1_tri2_vertex1 * p1_tri2_vertex3 > 0 )  
246.     {  
247.         if ( is_point_within_triangle( tri1, tri2->Vertex_2 ) )  
248.         {  
249.             return INTERSECT;  
250.         }  
251.         else  
252.         {  
253.             return NONINTERSECT;  
254.         }  
255.     }  
256.     else if ( p1_tri2_vertex3 == 0 && p1_tri2_vertex1 * p1_tri2_vertex2 > 0 )  
257.     {  
258.         if ( is_point_within_triangle( tri1, tri2->Vertex_3 ) )  
259.         {  
260.             return INTERSECT;  
261.         }  
262.         else  
263.         {  
264.             return NONINTERSECT;  
265.         }  
266.     }  
267.   
268.   
269.   
270.     float p2_tri1_vertex1 = get_vector4_det( tri2->Vertex_1, tri2->Vertex_2, tri2->Vertex_3, tri1->Vertex_1 );  
271.   
272.     float p2_tri1_vertex2 = get_vector4_det( tri2->Vertex_1, tri2->Vertex_2, tri2->Vertex_3, tri1->Vertex_2 );  
273.   
274.     float p2_tri1_vertex3 = get_vector4_det( tri2->Vertex_1, tri2->Vertex_2, tri2->Vertex_3, tri1->Vertex_3 );  
275.   
276.   
277.     if ( p2_tri1_vertex1 > 0 && p2_tri1_vertex2 > 0 && p2_tri1_vertex3 > 0 )  
278.     {  
279.         return NONINTERSECT;  
280.     }  
281.   
282.     if ( p2_tri1_vertex1 < 0 && p2_tri1_vertex2 < 0 && p2_tri1_vertex3 < 0 )  
283.     {  
284.         return NONINTERSECT;  
285.     }  
286.   
287.   
288.     if ( p2_tri1_vertex1 == 0 && p2_tri1_vertex2 * p2_tri1_vertex3 > 0 )  
289.     {  
290.         if ( is_point_within_triangle( tri2, tri1->Vertex_1 ) )  
291.         {  
292.             return INTERSECT;  
293.         }  
294.         else  
295.         {  
296.             return NONINTERSECT;  
297.         }  
298.     }  
299.   
300.     if ( p2_tri1_vertex2 == 0 && p2_tri1_vertex1 * p2_tri1_vertex3 > 0 )  
301.     {  
302.         if ( is_point_within_triangle( tri2, tri1->Vertex_2 ) )  
303.         {  
304.             return INTERSECT;  
305.         }  
306.         else  
307.         {  
308.             return NONINTERSECT;  
309.         }  
310.     }  
311.   
312.     if ( p2_tri1_vertex3 == 0 && p2_tri1_vertex1 * p2_tri1_vertex2 > 0 )  
313.     {  
314.         if ( is_point_within_triangle( tri2, tri1->Vertex_3 ) )  
315.         {  
316.             return INTERSECT;  
317.         }  
318.         else  
319.         {  
320.             return NONINTERSECT;  
321.         }  
322.     }  
323.   
324.   
325.   
326.     float* tri1_a = tri1->Vertex_1,* tri1_b = tri1->Vertex_2,* tri1_c = tri1->Vertex_3  
327.            ,* tri2_a = tri2->Vertex_1,* tri2_b = tri2->Vertex_2,* tri2_c = tri2->Vertex_3;  
328.   
329.     float* m;  
330.   
331.     float im;  
332.   
333.     if ( p2_tri1_vertex2 * p2_tri1_vertex3 >= 0 && p2_tri1_vertex1 != 0 )  
334.     {  
335.         if ( p2_tri1_vertex1 < 0 )  
336.         {  
337.             m = tri2_b;  
338.             tri2_b = tri2_c;  
339.             tri2_c = m;  
340.   
341.             im = p1_tri2_vertex2;  
342.             p1_tri2_vertex2 = p1_tri2_vertex3;  
343.             p1_tri2_vertex3 = im;  
344.         }  
345.     }  
346.     else if ( p2_tri1_vertex1 * p2_tri1_vertex3 >= 0 && p2_tri1_vertex2 != 0 )  
347.     {  
348.         m = tri1_a;  
349.         tri1_a = tri1_b;  
350.         tri1_b = tri1_c;  
351.         tri1_c = m;  
352.   
353.         if ( p2_tri1_vertex2 < 0 )  
354.         {  
355.             m = tri2_b;  
356.             tri2_b = tri2_c;  
357.             tri2_c = m;  
358.   
359.             im = p1_tri2_vertex2;  
360.             p1_tri2_vertex2 = p1_tri2_vertex3;  
361.             p1_tri2_vertex3 = im;  
362.         }  
363.     }  
364.     else if ( p2_tri1_vertex1 * p2_tri1_vertex2 >= 0 && p2_tri1_vertex3 != 0 )  
365.     {  
366.         m = tri1_a;  
367.   
368.         tri1_a = tri1_c;  
369.   
370.         tri1_c = tri1_b;  
371.   
372.         tri1_b = m;  
373.   
374.         if ( p2_tri1_vertex3 < 0 )  
375.         {  
376.             m = tri2_b;  
377.             tri2_b = tri2_c;  
378.             tri2_c = m;  
379.   
380.             im = p1_tri2_vertex2;  
381.             p1_tri2_vertex2 = p1_tri2_vertex3;  
382.             p1_tri2_vertex3 = im;  
383.         }  
384.     }  
385.   
386.     if ( p1_tri2_vertex2 * p1_tri2_vertex3 >= 0 && p1_tri2_vertex1 != 0 )  
387.     {  
388.         if ( p1_tri2_vertex1 < 0 )  
389.         {  
390.             m = tri1_b;  
391.             tri1_b = tri1_c;  
392.             tri1_c = m;  
393.         }  
394.     }  
395.     else if ( p1_tri2_vertex1 * p1_tri2_vertex3 >= 0 && p1_tri2_vertex2 != 0 )  
396.     {  
397.         m = tri2_a;  
398.   
399.         tri2_a = tri2_b;  
400.   
401.         tri2_b = tri2_c;  
402.   
403.         tri2_c = m;  
404.   
405.         if ( p1_tri2_vertex2 < 0 )  
406.         {  
407.             m = tri1_b;  
408.             tri1_b = tri1_c;  
409.             tri1_c = m;  
410.         }  
411.     }  
412.     else if ( p1_tri2_vertex1 * p1_tri2_vertex2 >= 0 && p1_tri2_vertex3 != 0 )  
413.     {  
414.         m = tri2_a;  
415.   
416.         tri2_a = tri2_c;  
417.   
418.         tri2_c = tri2_b;  
419.   
420.         tri2_b = m;  
421.   
422.         if ( p1_tri2_vertex3 < 0 )  
423.         {  
424.             m = tri1_b;  
425.             tri1_b = tri1_c;  
426.             tri1_c = m;  
427.         }  
428.     }  
429.   
430.     if ( get_vector4_det( tri1_a, tri1_b, tri2_a, tri2_b ) <= 0 && get_vector4_det( tri1_a, tri1_c, tri2_c, tri2_a ) <= 0 )  
431.     {  
432.         return INTERSECT;  
433.     }  
434.     else  
435.     {  
436.         return NONINTERSECT;  
437.     }  
438. }  

参考资料：

[1]     空间三角形快速相交检测算法     邹益胜, 丁国富, 何邕, 许明恒( 西南交通大学 机械工程学院)

[2]     http://www.cnblogs.com/graphics/archive/2010/08/05/1793393.html

更多 0