WebSVN – gelsvn – Diff – /trunk/src/Geometry/Triangle.cpp


/* ----------------------------------------------------------------------- *
 * This file is part of GEL, http://www.imm.dtu.dk/GEL
 * Copyright (C) the authors and DTU Informatics
 * For license and list of authors, see ../../doc/intro.pdf
 * ----------------------------------------------------------------------- */
 
#include "CGLA/Vec4f.h"
#include "Triangle.h"
 
 
using namespace std;
using namespace CGLA;
 
namespace Geometry
{
 
const float EPSILON = 1e-10f;
 
Triangle::Triangle(const CGLA::Vec3f& _v0, 
									 const CGLA::Vec3f& _v1, 
									 const CGLA::Vec3f& _v2,
									 
									 const CGLA::Vec3f& _vn0,
									 const CGLA::Vec3f& _vn1,
									 const CGLA::Vec3f& _vn2,
									 
									 const CGLA::Vec3f& _en0,
									 const CGLA::Vec3f& _en1,
									 const CGLA::Vec3f& _en2)
{
	vert[0]  =_v0;
	vert[1]  =_v1;
	vert[2]  =_v2;
 
 
	vert_norm[0] = _vn0;
	vert_norm[1] = _vn1;
	vert_norm[2] = _vn2;
 
	edge_norm[0] = _en0;
	edge_norm[1] = _en1;
	edge_norm[2] = _en2;
 
 
	face_norm = normalize(cross(vert[1]-vert[0], vert[2]-vert[0]));
	for(int i=0;i<3;++i)
		{
			int j= (i+1)%3;
			edge[i] = vert[j]-vert[i];
			tri_plane_edge_norm[i] = cross(face_norm, edge[i]);
			edge_len[i] = edge[i].length();
		}
}
 
 
// Moellers method
bool Triangle::intersect(const CGLA::Vec3f& orig,
												 const CGLA::Vec3f& dir, float&t) const
{
	Vec3f tvec, pvec, qvec;
	float det,inv_det;
 
   /* begin calculating determinant - also used to calculate U parameter */
   pvec = cross(dir, -edge[2]);
 
   /* if determinant is near zero, ray lies in plane of triangle */
   det = dot(edge[0], pvec);
 
   if (det > -EPSILON && det < EPSILON)
     return 0;
   inv_det = 1.0 / det;
 
   /* calculate distance from v0 to ray origin */
   tvec =  orig - vert[0];
 
   /* calculate U parameter and test bounds */
   float u = dot(tvec, pvec) * inv_det;
   if (u < 0.0 || u > 1.0)
     return false;
 
   /* prepare to test V parameter */
   qvec = cross(tvec, edge[0]);
 
   /* calculate V parameter and test bounds */
   float v = dot(dir, qvec) * inv_det;
   if (v < 0.0 || u + v > 1.0)
     return false;
 
   /* calculate t, ray intersects triangle */
   t = dot(-edge[2], qvec) * inv_det;
 
   return true;
}
 
 
 
 
bool Triangle::signed_distance(const Vec3f& p, 
															 float& sq_dist, float& sgn) const
{
	int vertex_scores[3] = {0,0,0};
	Vec3f closest_pnt, normal;
	int idx_0;
 
	// Loop over all three edges.
	for(idx_0=0; idx_0<3; ++idx_0)
		{
			const int idx_1 = (idx_0+1) % 3;
			const Vec3f dir_3d = edge[idx_0]/edge_len[idx_0];
			const float t = dot(p - vert[idx_0], dir_3d);
			if(t <= 0)
				{
					++vertex_scores[idx_0];
					if(vertex_scores[idx_0] == 2)
						{
							closest_pnt = vert[idx_0];
 
							normal = vert_norm[idx_0];
 
							break;
						}
				}
			else if(t >= edge_len[idx_0])
				{
					++vertex_scores[idx_1];
					if(vertex_scores[idx_1] == 2)
						{
							closest_pnt = vert[idx_1];
 
							normal = vert_norm[idx_1];
 
							break;
						}
				}
			else if(dot(tri_plane_edge_norm[idx_0], p-vert[idx_0]) <=0)
				{
					closest_pnt=vert[idx_0]+t*dir_3d;
 
					normal = edge_norm[idx_0];
 
					break;
				}
		}
	if(idx_0 == 3)
		{
			closest_pnt = p - face_norm*(dot(p-vert[0],face_norm));
 
			normal = face_norm;
 
		}
	sq_dist = sqr_length(p-closest_pnt);
	
 
	// Compute dot product with angle weighted normal, and
	// assign the sign based on the result.
	if(dot(normal, p-closest_pnt) >=0)
		sgn = 1.0f;
	else
		sgn = -1.0f;
 
 
 
 
	return true;
}
 
}
 

Rev 373	Rev 595
-		1	`/* ----------------------------------------------------------------------- *`
-		2	`* This file is part of GEL, http://www.imm.dtu.dk/GEL`
-		3	`* Copyright (C) the authors and DTU Informatics`
-		4	`* For license and list of authors, see ../../doc/intro.pdf`
-		5	`* ----------------------------------------------------------------------- */`
-		6
1	`#include "CGLA/Vec4f.h"`	7	`#include "CGLA/Vec4f.h"`
2	`#include "Triangle.h"`	8	`#include "Triangle.h"`
3		9
4		10
5	`using namespace std;`	11	`using namespace std;`
6	`using namespace CGLA;`	12	`using namespace CGLA;`
7		13
8	`namespace Geometry`	14	`namespace Geometry`
9	`{`	15	`{`
10		16
11	`const float EPSILON = 1e-10f;`	17	`const float EPSILON = 1e-10f;`
12		18
13	`Triangle::Triangle(const CGLA::Vec3f& _v0,`	19	`Triangle::Triangle(const CGLA::Vec3f& _v0,`
14	`const CGLA::Vec3f& _v1,`	20	`const CGLA::Vec3f& _v1,`
15	`const CGLA::Vec3f& _v2,`	21	`const CGLA::Vec3f& _v2,`
16		22
17	`const CGLA::Vec3f& _vn0,`	23	`const CGLA::Vec3f& _vn0,`
18	`const CGLA::Vec3f& _vn1,`	24	`const CGLA::Vec3f& _vn1,`
19	`const CGLA::Vec3f& _vn2,`	25	`const CGLA::Vec3f& _vn2,`
20		26
21	`const CGLA::Vec3f& _en0,`	27	`const CGLA::Vec3f& _en0,`
22	`const CGLA::Vec3f& _en1,`	28	`const CGLA::Vec3f& _en1,`
23	`const CGLA::Vec3f& _en2)`	29	`const CGLA::Vec3f& _en2)`
24	`{`	30	`{`
25	`vert[0] =_v0;`	31	`vert[0] =_v0;`
26	`vert[1] =_v1;`	32	`vert[1] =_v1;`
27	`vert[2] =_v2;`	33	`vert[2] =_v2;`
28		34
29	`#ifdef COMPUTE_SIGN`	-
30	`vert_norm[0] = _vn0;`	35	`vert_norm[0] = _vn0;`
31	`vert_norm[1] = _vn1;`	36	`vert_norm[1] = _vn1;`
32	`vert_norm[2] = _vn2;`	37	`vert_norm[2] = _vn2;`
33		38
34	`edge_norm[0] = _en0;`	39	`edge_norm[0] = _en0;`
35	`edge_norm[1] = _en1;`	40	`edge_norm[1] = _en1;`
36	`edge_norm[2] = _en2;`	41	`edge_norm[2] = _en2;`
37	`#endif`	-
38		42
39	`face_norm = normalize(cross(vert[1]-vert[0], vert[2]-vert[0]));`	43	`face_norm = normalize(cross(vert[1]-vert[0], vert[2]-vert[0]));`
40	`for(int i=0;i<3;++i)`	44	`for(int i=0;i<3;++i)`
41	`{`	45	`{`
42	`int j= (i+1)%3;`	46	`int j= (i+1)%3;`
43	`edge[i] = vert[j]-vert[i];`	47	`edge[i] = vert[j]-vert[i];`
44	`tri_plane_edge_norm[i] = cross(face_norm, edge[i]);`	48	`tri_plane_edge_norm[i] = cross(face_norm, edge[i]);`
45	`edge_len[i] = edge[i].length();`	49	`edge_len[i] = edge[i].length();`
46	`}`	50	`}`
47	`}`	51	`}`
48		52
49		53
50	`// Moellers method`	54	`// Moellers method`
51	`bool Triangle::intersect(const CGLA::Vec3f& orig,`	55	`bool Triangle::intersect(const CGLA::Vec3f& orig,`
52	`const CGLA::Vec3f& dir, float&t) const`	56	`const CGLA::Vec3f& dir, float&t) const`
53	`{`	57	`{`
54	`Vec3f tvec, pvec, qvec;`	58	`Vec3f tvec, pvec, qvec;`
55	`float det,inv_det;`	59	`float det,inv_det;`
56		60
57	`/* begin calculating determinant - also used to calculate U parameter */`	61	`/* begin calculating determinant - also used to calculate U parameter */`
58	`pvec = cross(dir, -edge[2]);`	62	`pvec = cross(dir, -edge[2]);`
59		63
60	`/* if determinant is near zero, ray lies in plane of triangle */`	64	`/* if determinant is near zero, ray lies in plane of triangle */`
61	`det = dot(edge[0], pvec);`	65	`det = dot(edge[0], pvec);`
62		66
63	`if (det > -EPSILON && det < EPSILON)`	67	`if (det > -EPSILON && det < EPSILON)`
64	`return 0;`	68	`return 0;`
65	`inv_det = 1.0 / det;`	69	`inv_det = 1.0 / det;`
66		70
67	`/* calculate distance from v0 to ray origin */`	71	`/* calculate distance from v0 to ray origin */`
68	`tvec = orig - vert[0];`	72	`tvec = orig - vert[0];`
69		73
70	`/* calculate U parameter and test bounds */`	74	`/* calculate U parameter and test bounds */`
71	`float u = dot(tvec, pvec) * inv_det;`	75	`float u = dot(tvec, pvec) * inv_det;`
72	`if (u < 0.0 \|\| u > 1.0)`	76	`if (u < 0.0 \|\| u > 1.0)`
73	`return false;`	77	`return false;`
74		78
75	`/* prepare to test V parameter */`	79	`/* prepare to test V parameter */`
76	`qvec = cross(tvec, edge[0]);`	80	`qvec = cross(tvec, edge[0]);`
77		81
78	`/* calculate V parameter and test bounds */`	82	`/* calculate V parameter and test bounds */`
79	`float v = dot(dir, qvec) * inv_det;`	83	`float v = dot(dir, qvec) * inv_det;`
80	`if (v < 0.0 \|\| u + v > 1.0)`	84	`if (v < 0.0 \|\| u + v > 1.0)`
81	`return false;`	85	`return false;`
82		86
83	`/* calculate t, ray intersects triangle */`	87	`/* calculate t, ray intersects triangle */`
84	`t = dot(-edge[2], qvec) * inv_det;`	88	`t = dot(-edge[2], qvec) * inv_det;`
85		89
86	`return true;`	90	`return true;`
87	`}`	91	`}`
88		92
89		93
90		94
91		95
92	`bool Triangle::signed_distance(const Vec3f& p,`	96	`bool Triangle::signed_distance(const Vec3f& p,`
93	`float& sq_dist, float& sgn) const`	97	`float& sq_dist, float& sgn) const`
94	`{`	98	`{`
95	`int vertex_scores[3] = {0,0,0};`	99	`int vertex_scores[3] = {0,0,0};`
96	`Vec3f closest_pnt, normal;`	100	`Vec3f closest_pnt, normal;`
97	`int idx_0;`	101	`int idx_0;`
98		102
99	`// Loop over all three edges.`	103	`// Loop over all three edges.`
100	`for(idx_0=0; idx_0<3; ++idx_0)`	104	`for(idx_0=0; idx_0<3; ++idx_0)`
101	`{`	105	`{`
102	`const int idx_1 = (idx_0+1) % 3;`	106	`const int idx_1 = (idx_0+1) % 3;`
103	`const Vec3f dir_3d = edge[idx_0]/edge_len[idx_0];`	107	`const Vec3f dir_3d = edge[idx_0]/edge_len[idx_0];`
104	`const float t = dot(p - vert[idx_0], dir_3d);`	108	`const float t = dot(p - vert[idx_0], dir_3d);`
105	`if(t <= 0)`	109	`if(t <= 0)`
106	`{`	110	`{`
107	`++vertex_scores[idx_0];`	111	`++vertex_scores[idx_0];`
108	`if(vertex_scores[idx_0] == 2)`	112	`if(vertex_scores[idx_0] == 2)`
109	`{`	113	`{`
110	`closest_pnt = vert[idx_0];`	114	`closest_pnt = vert[idx_0];`
111	`#ifdef COMPUTE_SIGN`	-
112	`normal = vert_norm[idx_0];`	115	`normal = vert_norm[idx_0];`
113	`#endif`	-
114	`break;`	116	`break;`
115	`}`	117	`}`
116	`}`	118	`}`
117	`else if(t >= edge_len[idx_0])`	119	`else if(t >= edge_len[idx_0])`
118	`{`	120	`{`
119	`++vertex_scores[idx_1];`	121	`++vertex_scores[idx_1];`
120	`if(vertex_scores[idx_1] == 2)`	122	`if(vertex_scores[idx_1] == 2)`
121	`{`	123	`{`
122	`closest_pnt = vert[idx_1];`	124	`closest_pnt = vert[idx_1];`
123	`#ifdef COMPUTE_SIGN`	-
124	`normal = vert_norm[idx_1];`	125	`normal = vert_norm[idx_1];`
125	`#endif`	-
126	`break;`	126	`break;`
127	`}`	127	`}`
128	`}`	128	`}`
129	`else if(dot(tri_plane_edge_norm[idx_0], p-vert[idx_0]) <=0)`	129	`else if(dot(tri_plane_edge_norm[idx_0], p-vert[idx_0]) <=0)`
130	`{`	130	`{`
131	`closest_pnt=vert[idx_0]+t*dir_3d;`	131	`closest_pnt=vert[idx_0]+t*dir_3d;`
132	`#ifdef COMPUTE_SIGN`	-
133	`normal = edge_norm[idx_0];`	132	`normal = edge_norm[idx_0];`
134	`#endif`	-
135	`break;`	133	`break;`
136	`}`	134	`}`
137	`}`	135	`}`
138	`if(idx_0 == 3)`	136	`if(idx_0 == 3)`
139	`{`	137	`{`
140	`closest_pnt = p - face_norm*(dot(p-vert[0],face_norm));`	138	`closest_pnt = p - face_norm*(dot(p-vert[0],face_norm));`
141	`#ifdef COMPUTE_SIGN`	-
142	`normal = face_norm;`	139	`normal = face_norm;`
143	`#endif`	-
144	`}`	140	`}`
145	`sq_dist = sqr_length(p-closest_pnt);`	141	`sq_dist = sqr_length(p-closest_pnt);`
146		142
147	`#ifdef COMPUTE_SIGN`	-
148	`// Compute dot product with angle weighted normal, and`	143	`// Compute dot product with angle weighted normal, and`
149	`// assign the sign based on the result.`	144	`// assign the sign based on the result.`
150	`if(dot(normal, p-closest_pnt) >=0)`	145	`if(dot(normal, p-closest_pnt) >=0)`
151	`sgn = 1.0f;`	146	`sgn = 1.0f;`
152	`else`	147	`else`
153	`sgn = -1.0f;`	148	`sgn = -1.0f;`
154	`#else`	-
155	`sgn = 1.0f;`	-
156	`#endif`	-
157		149
158	`return true;`	150	`return true;`
159	`}`	151	`}`
160		152
161	`}`	153	`}`
162		154

Subversion Repositories gelsvn

(root)/trunk/src/Geometry/Triangle.cpp – Rev 373 → 595