WebSVN – gelsvn – Diff – /branches/SPD/src/Geometry/Triangle.cpp


#include "CGLA/Vec4f.h"
#include "Triangle.h"
 
 
using namespace std;
using namespace CGLA;
 
namespace Geometry
{
 
const float EPSILON = 0.000001f;
 
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;
 
#ifdef COMPUTE_SIGN
	vert_norm[0] = _vn0;
	vert_norm[1] = _vn1;
	vert_norm[2] = _vn2;
 
	edge_norm[0] = _en0;
	edge_norm[1] = _en1;
	edge_norm[2] = _en2;
#endif
 
	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];
#ifdef COMPUTE_SIGN
							normal = vert_norm[idx_0];
#endif
							break;
						}
				}
			else if(t >= edge_len[idx_0])
				{
					++vertex_scores[idx_1];
					if(vertex_scores[idx_1] == 2)
						{
							closest_pnt = vert[idx_1];
#ifdef COMPUTE_SIGN
							normal = vert_norm[idx_1];
#endif
							break;
						}
				}
			else if(dot(tri_plane_edge_norm[idx_0], p-vert[idx_0]) <=0)
				{
					closest_pnt=vert[idx_0]+t*dir_3d;
#ifdef COMPUTE_SIGN
					normal = edge_norm[idx_0];
#endif
					break;
				}
		}
	if(idx_0 == 3)
		{
			closest_pnt = p - face_norm*(dot(p-vert[0],face_norm));
#ifdef COMPUTE_SIGN
			normal = face_norm;
#endif
		}
	sq_dist = sqr_length(p-closest_pnt);
	
#ifdef COMPUTE_SIGN
	// 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;
#else
	sgn = 1.0f;
#endif
 
	return true;
}
 
 
}
 
Generated by GNU Enscript 1.6.6.
 
 
 

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

Subversion Repositories gelsvn

(root)/branches/SPD/src/Geometry/Triangle.cpp – Rev 290 → 291