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

最新推荐文章于 2023-02-05 19:28:34 发布
最新推荐文章于 2023-02-05 19:28:34 发布
阅读量9.5k 收藏 26
点赞数 8
文章标签： 空间三角形 快速 相交 c++
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/fourierfeng/article/details/11969915
版权

    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期，邹益胜等的《空间三角形快速相交检测算法》下面是具体代码实现。

 

typedef float float3[3];

enum TopologicalStructure 
{ 
	INTERSECT, NONINTERSECT 
};

struct Triangle
{
	float3 Normal_0; 
	float3 Vertex_1, Vertex_2, Vertex_3;
};

struct point
{
	float x, y;
};


//三维点拷贝为二维点
static void copy_point( point& p, float3 f )
{
	p.x = f[0];
	p.y = f[1];
}

//四点行列式
inline float get_vector4_det( float3 v1, float3 v2, float3 v3, float3 v4 )
{
	float a[3][3];
	for ( int i = 0; i != 3; ++i )
	{
		a[0][i] = v1[i] - v4[i];
		a[1][i] = v2[i] - v4[i];
		a[2][i] = v3[i] - v4[i];
	}

	return a[0][0] * a[1][1] * a[2][2] 
	+ a[0][1] * a[1][2] * a[2][0] 
	+ a[0][2] * a[1][0] * a[2][1] 
	- a[0][2] * a[1][1] * a[2][0] 
	- a[0][1] * a[1][0] * a[2][2] 
	- a[0][0] * a[1][2] * a[2][1];
}


//利用叉积计算点p相对线段p1p2的方位
inline double direction( point p1, point p2, point p ){
	return ( p.x - p1.x ) * ( p2.y - p1.y ) - ( p2.x - p1.x ) * ( p.y - p1.y )   ;
}


//确定与线段p1p2共线的点p是否在线段p1p2上
inline int on_segment( point p1, point p2, point p ){
	double max = p1.x > p2.x ? p1.x : p2.x ;
	double min = p1.x < p2.x ? p1.x : p2.x ;
	double max1 = p1.y > p2.y ? p1.y : p2.y ;
	double min1 = p1.y < p2.y ? p1.y : p2.y ;
	if ( p.x >= min && p.x <= max && p.y >= min1 && p.y <= max1 )
	{
		return 1;
	}
	else
	{
		return 0;
	}
}


//判断线段p1p2与线段p3p4是否相交的主函数
inline int segments_intersert( point p1, point p2, point p3, point p4 ){
	double d1, d2, d3, d4;
	d1 = direction( p3, p4, p1 );
	d2 = direction( p3, p4, p2 );
	d3 = direction( p1, p2, p3 );
	d4 = direction( p1, p2, p4 );
	if ( d1 * d2 < 0 && d3 * d4 < 0 )
	{
		return 1;
	}
	else if ( d1 == 0 && on_segment( p3, p4, p1 ) == 1 )
	{
		return 1;
	}
	else if ( d2 == 0 && on_segment( p3, p4, p2 ) == 1 )
	{
		return 1;
	}
	else if ( d3 == 0 && on_segment( p1, p2, p3 ) == 1 )
	{
		return 1;
	}
	else if ( d4 == 0 && on_segment( p1, p2, p4 ) == 1 )
	{
		return 1;
	}
	return 0;
}


//判断同一平面的直线和三角形是否相交
inline bool line_triangle_intersert_inSamePlane( Triangle* tri, float3 f1, float3 f2 )
{
	point p1, p2, p3, p4;

	copy_point( p1, f1 );

	copy_point( p2, f2 );

	copy_point( p3, tri->Vertex_1 );

	copy_point( p4, tri->Vertex_2 );

	if ( segments_intersert( p1, p2, p3, p4 ) )
	{
		return true;
	}

	copy_point( p3, tri->Vertex_2 );

	copy_point( p4, tri->Vertex_3 );

	if ( segments_intersert( p1, p2, p3, p4 ) )
	{
		return true;
	}

	copy_point( p3, tri->Vertex_1 );

	copy_point( p4, tri->Vertex_3 );

	if ( segments_intersert( p1, p2, p3, p4 ) )
	{
		return true;
	}

	return false;
} 


inline void get_central_point(float3 ¢ralPoint,Triangle* tri)
{
	centralPoint[0]=(tri->Vertex_1[0]+tri->Vertex_2[0]+tri->Vertex_3[0])/3;

	centralPoint[1]=(tri->Vertex_1[1]+tri->Vertex_2[1]+tri->Vertex_3[1])/3;

	centralPoint[2]=(tri->Vertex_1[2]+tri->Vertex_2[2]+tri->Vertex_3[2])/3;
}




//向量之差
inline void get_vector_diff( float3& aimV, const float3 a, const float3 b )
{
	aimV[0] = b[0] - a[0];

	aimV[1] = b[1] - a[1];

	aimV[2] = b[2] - a[2];
} 

//向量内积
inline float Dot( const float3& v1, const float3& v2 )
{    
	return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2] ;
} 

//重心法判断点是否在三角形内部
inline bool is_point_within_triangle( Triangle* tri, float3 point )
{
	float3 v0;   
	get_vector_diff( v0, tri->Vertex_1, tri->Vertex_3 );
	float3 v1;   
	get_vector_diff( v1, tri->Vertex_1, tri->Vertex_2 );
	float3 v2; 
	get_vector_diff( v2, tri->Vertex_1, point );
	float dot00 = Dot( v0, v0 ) ;   
	float dot01 = Dot( v0, v1 ) ;    
	float dot02 = Dot( v0, v2 ) ;    
	float dot11 = Dot( v1, v1 ) ;   
	float dot12 = Dot( v1, v2 ) ;    
	float inverDeno = 1 / ( dot00* dot11 - dot01* dot01 ) ;    
	float u = ( dot11* dot02 - dot01* dot12 ) * inverDeno ;    
	if ( u < 0 || u > 1 ) // if u out of range, return directly
	{
		return false ;
	}    
	float v = ( dot00* dot12 - dot01* dot02 ) * inverDeno ;    
	if ( v < 0 || v > 1 ) // if v out of range, return directly
	{
		return false ;
	}    
	return u + v <= 1 ;
} 


//判断同一平面内的三角形是否相交
inline bool triangle_intersert_inSamePlane( Triangle* tri1, Triangle* tri2 )
{
	if ( line_triangle_intersert_inSamePlane( tri2, tri1->Vertex_1, tri1->Vertex_2 ) )
	{
		return true;
	}
	else if ( line_triangle_intersert_inSamePlane( tri2, tri1->Vertex_2, tri1->Vertex_3 ) )
	{
		return true;
	}
	else if ( line_triangle_intersert_inSamePlane( tri2, tri1->Vertex_1, tri1->Vertex_3 ) )
	{
		return true;
	}
	else
	{
		float3 centralPoint1,centralPoint2;

		get_central_point(centralPoint1,tri1);

		get_central_point(centralPoint2,tri2);

		if(is_point_within_triangle( tri2, centralPoint1 ) || is_point_within_triangle( tri1, centralPoint2 ) )
		{
			return true;
		}

		return false;
	}
} 

//Devillers算法主函数
inline TopologicalStructure judge_triangle_topologicalStructure( Triangle* tri1, Triangle* tri2 )
{
	//设tri1所在的平面为p1,tri2所在的平面为p2
	float p1_tri2_vertex1 = get_vector4_det( tri1->Vertex_1, tri1->Vertex_2, tri1->Vertex_3, tri2->Vertex_1 );

	float p1_tri2_vertex2 = get_vector4_det( tri1->Vertex_1, tri1->Vertex_2, tri1->Vertex_3, tri2->Vertex_2 );

	float p1_tri2_vertex3 = get_vector4_det( tri1->Vertex_1, tri1->Vertex_2, tri1->Vertex_3, tri2->Vertex_3 );


	if ( p1_tri2_vertex1 > 0 && p1_tri2_vertex2 > 0 && p1_tri2_vertex3 > 0 )
	{
		return NONINTERSECT;
	}

	if ( p1_tri2_vertex1 < 0 && p1_tri2_vertex2 < 0 && p1_tri2_vertex3 < 0 )
	{
		return NONINTERSECT;
	}


	if ( p1_tri2_vertex1 == 0 && p1_tri2_vertex2 == 0 && p1_tri2_vertex3 == 0 )
	{
		if ( triangle_intersert_inSamePlane( tri1, tri2 ) )
		{
			return INTERSECT;
		}
		else
		{
			return NONINTERSECT;
		}
	}


	if ( p1_tri2_vertex1 == 0 && p1_tri2_vertex2 * p1_tri2_vertex3 > 0 )
	{
		if ( is_point_within_triangle( tri1, tri2->Vertex_1 ) )
		{
			return INTERSECT;
		}
		else
		{
			return NONINTERSECT;
		}
	}
	else if ( p1_tri2_vertex2 == 0 && p1_tri2_vertex1 * p1_tri2_vertex3 > 0 )
	{
		if ( is_point_within_triangle( tri1, tri2->Vertex_2 ) )
		{
			return INTERSECT;
		}
		else
		{
			return NONINTERSECT;
		}
	}
	else if ( p1_tri2_vertex3 == 0 && p1_tri2_vertex1 * p1_tri2_vertex2 > 0 )
	{
		if ( is_point_within_triangle( tri1, tri2->Vertex_3 ) )
		{
			return INTERSECT;
		}
		else
		{
			return NONINTERSECT;
		}
	}



	float p2_tri1_vertex1 = get_vector4_det( tri2->Vertex_1, tri2->Vertex_2, tri2->Vertex_3, tri1->Vertex_1 );

	float p2_tri1_vertex2 = get_vector4_det( tri2->Vertex_1, tri2->Vertex_2, tri2->Vertex_3, tri1->Vertex_2 );

	float p2_tri1_vertex3 = get_vector4_det( tri2->Vertex_1, tri2->Vertex_2, tri2->Vertex_3, tri1->Vertex_3 );


	if ( p2_tri1_vertex1 > 0 && p2_tri1_vertex2 > 0 && p2_tri1_vertex3 > 0 )
	{
		return NONINTERSECT;
	}

	if ( p2_tri1_vertex1 < 0 && p2_tri1_vertex2 < 0 && p2_tri1_vertex3 < 0 )
	{
		return NONINTERSECT;
	}


	if ( p2_tri1_vertex1 == 0 && p2_tri1_vertex2 * p2_tri1_vertex3 > 0 )
	{
		if ( is_point_within_triangle( tri2, tri1->Vertex_1 ) )
		{
			return INTERSECT;
		}
		else
		{
			return NONINTERSECT;
		}
	}

	if ( p2_tri1_vertex2 == 0 && p2_tri1_vertex1 * p2_tri1_vertex3 > 0 )
	{
		if ( is_point_within_triangle( tri2, tri1->Vertex_2 ) )
		{
			return INTERSECT;
		}
		else
		{
			return NONINTERSECT;
		}
	}

	if ( p2_tri1_vertex3 == 0 && p2_tri1_vertex1 * p2_tri1_vertex2 > 0 )
	{
		if ( is_point_within_triangle( tri2, tri1->Vertex_3 ) )
		{
			return INTERSECT;
		}
		else
		{
			return NONINTERSECT;
		}
	}



	float* tri1_a = tri1->Vertex_1,* tri1_b = tri1->Vertex_2,* tri1_c = tri1->Vertex_3
		,* tri2_a = tri2->Vertex_1,* tri2_b = tri2->Vertex_2,* tri2_c = tri2->Vertex_3;

	float* m;

	float im;

	if ( p2_tri1_vertex2 * p2_tri1_vertex3 >= 0 && p2_tri1_vertex1 != 0 )
	{
		if ( p2_tri1_vertex1 < 0 )
		{
			m = tri2_b;
			tri2_b = tri2_c;
			tri2_c = m;

			im = p1_tri2_vertex2;
			p1_tri2_vertex2 = p1_tri2_vertex3;
			p1_tri2_vertex3 = im;
		}
	}
	else if ( p2_tri1_vertex1 * p2_tri1_vertex3 >= 0 && p2_tri1_vertex2 != 0 )
	{
		m = tri1_a;
		tri1_a = tri1_b;
		tri1_b = tri1_c;
		tri1_c = m;

		if ( p2_tri1_vertex2 < 0 )
		{
			m = tri2_b;
			tri2_b = tri2_c;
			tri2_c = m;

			im = p1_tri2_vertex2;
			p1_tri2_vertex2 = p1_tri2_vertex3;
			p1_tri2_vertex3 = im;
		}
	}
	else if ( p2_tri1_vertex1 * p2_tri1_vertex2 >= 0 && p2_tri1_vertex3 != 0 )
	{
		m = tri1_a;

		tri1_a = tri1_c;

		tri1_c = tri1_b;

		tri1_b = m;

		if ( p2_tri1_vertex3 < 0 )
		{
			m = tri2_b;
			tri2_b = tri2_c;
			tri2_c = m;

			im = p1_tri2_vertex2;
			p1_tri2_vertex2 = p1_tri2_vertex3;
			p1_tri2_vertex3 = im;
		}
	}

	if ( p1_tri2_vertex2 * p1_tri2_vertex3 >= 0 && p1_tri2_vertex1 != 0 )
	{
		if ( p1_tri2_vertex1 < 0 )
		{
			m = tri1_b;
			tri1_b = tri1_c;
			tri1_c = m;
		}
	}
	else if ( p1_tri2_vertex1 * p1_tri2_vertex3 >= 0 && p1_tri2_vertex2 != 0 )
	{
		m = tri2_a;

		tri2_a = tri2_b;

		tri2_b = tri2_c;

		tri2_c = m;

		if ( p1_tri2_vertex2 < 0 )
		{
			m = tri1_b;
			tri1_b = tri1_c;
			tri1_c = m;
		}
	}
	else if ( p1_tri2_vertex1 * p1_tri2_vertex2 >= 0 && p1_tri2_vertex3 != 0 )
	{
		m = tri2_a;

		tri2_a = tri2_c;

		tri2_c = tri2_b;

		tri2_b = m;

		if ( p1_tri2_vertex3 < 0 )
		{
			m = tri1_b;
			tri1_b = tri1_c;
			tri1_c = m;
		}
	}

	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 )
	{
		return INTERSECT;
	}
	else
	{
		return NONINTERSECT;
	}
}


参考资料：

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

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

玉名卍
关注
directx 三角片相交快速检测算法综述
10-12
三角形是构成三维模型的基本单元，因此能够快速准确地检测两个三角形是否相交对于场景中的碰撞检测、光线追踪、物理模拟等领域具有非常高的实用价值。对于需要进行大量三角形相交检测的应用，如三维游戏引擎、虚拟...
devillers.wso2is::locked_with_key: 修改 WSO2 身份服务器
06-17
6. **自动部署配置**：这可能涉及配置 CI/CD 管道，如 Jenkins 或 GitLab CI，使得每次代码更改都能自动构建和部署 WSO2IS，以实现快速迭代和更新。 7. **将现有测试用户迁移到 WSO2IS MSSQL 格式**：表明有数据...
三维空间三角形与线段相交判断
love code, love life.
06-17 967
#include <iostream> #include <numeric> #include <math.h> #include <stdio.h> class vec3 { public: float x, y, z; float dot(const vec3& b) { return vec3::x * b.x + vec3::y * b.y + vec3::z * b.z; } vec3 .
[算法]检测空间三角形相交算法(Devillers & Guigue算法)
weixin_33946020的博客
10-19 960
#pragma once //GYDevillersTriangle.h /* 快速检测空间三角形相交算法的代码实现(Devillers &amp; Guigue算法) 博客原地址:http://blog.csdn.net/fourierfeng/article/details/11969915# Devillers &amp; Guigue算法(简称Devillers 算法) 通过...
【JS实现】检测空间三角形相交算法(Devillers & Guigue算法)
sunqian的博客
11-07 668
用javascript实现判断空间三角形是否相交(Devillers & Guigue算法)
UVA 11275 判断空间三角形是否相交
smilesundream的博客
05-20 921
题意：训练指南292页#include#includeusing namespace std;struct Point3 { double x, y, z; Point3(double x=0, double y=0, double z=0):x(x),y(y),z(z) { } void read() { scanf("%lf%lf%lf", &x, &y, &z); }};t
三角形快速相交算法.doc
05-06
三角形快速相交算法是计算机图形学中的一个重要概念，它主要应用于三维空间中两个三角形是否相交的判断。在游戏开发、物理模拟、碰撞检测等领域有着广泛的应用。以下是对几种常见三角形相交算法的详细说明： 1. **...
空间数据质量
05-18
《空间数据质量：理论、分析与应用》是空间信息技术领域内一部权威且深入的著作，由Rodolphe Devillers和Robert Jeansoulin编辑。本书不仅涵盖了空间数据质量的基础概念，还深入探讨了其在实际应用中的理论与方法，...
三角形相交检测（碰撞）
04-22
描述了两个三角形相交的判断过程。vc++6.0界面设计。
检测空间三角形相交 C#
09-29
检测空间三角形相交 C#,通过计算行列式的值来判断空间内三角形点的分布情况，再计算是否相交，快速排除不相交的情况。
射线和三角形的相交检测（ray triangle intersection test）
03-15
射线和三角形的相交检测（ray triangle intersection test）
超斐波那契螺旋：快速、低差异采样SO(3)
该算法简单且快速。与其他方法进行的综合评估表明，生成的方向集具有低差异性，在功率谱中几乎没有杂散分量，并且几乎具有相同的Voronoi体积。这使它们在各种应用中非常有用，特别是在蒙特卡罗采样中。01. 引言0...
Olivier Devillers的空间三角形相交判断算法
谢公屐的博客
04-06 777
参考 Olivier Devillers 代数定义 二维空间下关于行列式的应用 定义点a=(ax,ay),b=(bx,by),c=(cx,cy)ca⃗=(ax−cx,ay−cy),bc⃗=(bx−cx,by−cy) 定义点a=(a_x,a_y),b=(b_x,b_y),c=(c_x,c_y) \\ \vec{ca}=(a_x-c_x,a_y-c_y),\vec{bc}=(b_x-c_x,b_y-c_y) 定义点a=(ax​,ay​),b=(bx​,by​),c=(cx​,cy​)ca=(ax​−cx​,
C++三角形相交算法(Devillers & Guigue算法)
u011208918的专栏
07-28 694
C++ 三角形相交检测
【数理知识】直线，射线与平面，三角形相交算法
最新发布
赵继超的笔记
02-05 436
直线，射线与平面，三角形相交算法
Möller-Trumbore算法-射线三角形相交算法
zhanxi1992的专栏
11-21 8325
Möller-Trumbore算法一、概述二、准备知识三、Möller-Trumbore 算法推导推导过程 一、概述 Möller-Trumbore 射线三角相交算法是一种快速计算射线与三角形在三个维度上的交点的方法，通过向量与矩阵计算可以快速得出交点与重心坐标，而无需对包含三角形的平面方程进行预计算。另外，它还应用于计算机图形学中以实现涉及三角形网格的光线跟踪计算。算法名字是以发明者 TomasMöller 和 Ben Trumbore 的名字来命名的。 二、准备知识 1、三阶方阵的行列式等于三个列向量的
C++实现空间中两个三角形位置关系（相交、平行）的判断
weixin_39655480的博客
03-21 1813
两个空间三角形相交关系的判断: 分别判断三角形trangle_1 的三条边与三角形trangle_2 的位置关系，三角形trangle_2 的三条边与三角形trangle_1 的位置关系。 --求解三角形trangle_1 的其中一条边的直线方程 --求解三角形trangle_2 所在平面的方程 --根据以上方程求得x,y,z --判断点(x,y...
评论 12
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值