Line data Source code
1 : /* Copyright (C) Dan Ginsburg, 2000.
2 : * All rights reserved worldwide.
3 : *
4 : * This software is provided "as is" without express or implied
5 : * warranties. You may freely copy and compile this source into
6 : * applications you distribute provided that the copyright text
7 : * below is included in the resulting source code, for example:
8 : * "Portions Copyright (C) Dan Ginsburg, 2000"
9 : */
10 : /*
11 : * This Code has been modified by Sebastian Reiter.
12 : * Changes were made to the underlying types and to the
13 : * return-value of TriangleBoxIntersection.
14 : * The original source-code was distributed with the book
15 : * 'Game Programming Gems'.
16 : */
17 :
18 : ////////////////////////////////////////////////////////////////////////////////////////
19 : //
20 : // TriBox.cpp
21 : //
22 : // Description:
23 : //
24 : // This Triangle-Box intersection code was adapated from the Graphics Gems III
25 : // source archive 'triangleCube.c' available at:
26 : //
27 : // http://www.acm.org/tog/GraphicsGems/
28 : //
29 : // The main modification is that the original code performed only intersection
30 : // with a voxel (that is, a unit cube centered at the origin). This code is
31 : // modified to take an arbitrary box and triangle and perform the scale/translation
32 : // necessary to get the triangle into voxel space.
33 : //
34 : //
35 : #include <math.h>
36 : #include "../ugmath.h"
37 :
38 : ///
39 : // Macros
40 : //
41 : #define INSIDE 0
42 : #define OUTSIDE 1
43 :
44 : #ifndef FALSE
45 : #define FALSE 0
46 : #endif
47 : #ifndef TRUE
48 : #define TRUE 1
49 : #endif
50 :
51 : #define EPS 10e-5
52 : #define SIGN3( A ) \
53 : (((A).x() < EPS ? 4 : 0) | ((A).x() > -EPS ? 32 : 0) | \
54 : ((A).y() < EPS ? 2 : 0) | ((A).y() > -EPS ? 16 : 0) | \
55 : ((A).z() < EPS ? 1 : 0) | ((A).z() > -EPS ? 8 : 0))
56 :
57 : #define CROSS( A, B, C ) { \
58 : (C).x() = (A).y() * (B).z() - (A).z() * (B).y(); \
59 : (C).y() = -(A).x() * (B).z() + (A).z() * (B).x(); \
60 : (C).z() = (A).x() * (B).y() - (A).y() * (B).x(); \
61 : }
62 : #define SUB( A, B, C ) { \
63 : (C).x() = (A).x() - (B).x(); \
64 : (C).y() = (A).y() - (B).y(); \
65 : (C).z() = (A).z() - (B).z(); \
66 : }
67 : #define LERP( A, B, C) ((B)+(A)*((C)-(B)))
68 : #define MIN3(a,b,c) ((((a)<(b))&&((a)<(c))) ? (a) : (((b)<(c)) ? (b) : (c)))
69 : #define MAX3(a,b,c) ((((a)>(b))&&((a)>(c))) ? (a) : (((b)>(c)) ? (b) : (c)))
70 :
71 :
72 : namespace ug
73 : {
74 :
75 : // a very simple triangle type
76 0 : struct TRI
77 : {
78 : vector3 m_P[3];
79 : };
80 :
81 : //////////////////////////////////////////////////////////////////////////////////////////
82 : //
83 : // Private Functions
84 : //
85 : //
86 :
87 : ///
88 : // FacePlane()
89 : //
90 : // Which of the six face-plane(s) is point P outside of?
91 : //
92 : static
93 0 : int FacePlane(const vector3& p)
94 : {
95 : int outcode;
96 :
97 : outcode = 0;
98 0 : if (p.x() > .5) outcode |= 0x01;
99 0 : if (p.x() < -.5) outcode |= 0x02;
100 0 : if (p.y() > .5) outcode |= 0x04;
101 0 : if (p.y() < -.5) outcode |= 0x08;
102 0 : if (p.z() > .5) outcode |= 0x10;
103 0 : if (p.z() < -.5) outcode |= 0x20;
104 :
105 0 : return(outcode);
106 : }
107 :
108 :
109 : ///
110 : // Bevel2d()
111 : //
112 : // Which of the twelve edge plane(s) is point P outside of?
113 : //
114 : static
115 0 : int Bevel2d(const vector3& p)
116 : {
117 : int outcode;
118 :
119 : outcode = 0;
120 0 : if ( p.x() + p.y() > 1.0) outcode |= 0x001;
121 0 : if ( p.x() - p.y() > 1.0) outcode |= 0x002;
122 0 : if (-p.x() + p.y() > 1.0) outcode |= 0x004;
123 0 : if (-p.x() - p.y() > 1.0) outcode |= 0x008;
124 0 : if ( p.x() + p.z() > 1.0) outcode |= 0x010;
125 0 : if ( p.x() - p.z() > 1.0) outcode |= 0x020;
126 0 : if (-p.x() + p.z() > 1.0) outcode |= 0x040;
127 0 : if (-p.x() - p.z() > 1.0) outcode |= 0x080;
128 0 : if ( p.y() + p.z() > 1.0) outcode |= 0x100;
129 0 : if ( p.y() - p.z() > 1.0) outcode |= 0x200;
130 0 : if (-p.y() + p.z() > 1.0) outcode |= 0x400;
131 0 : if (-p.y() - p.z() > 1.0) outcode |= 0x800;
132 0 : return(outcode);
133 : }
134 :
135 : ///
136 : // Bevel3d()
137 : //
138 : // Which of the eight corner plane(s) is point P outside of?
139 : //
140 : static
141 0 : int Bevel3d(const vector3& p)
142 : {
143 : int outcode;
144 :
145 : outcode = 0;
146 0 : if (( p.x() + p.y() + p.z()) > 1.5) outcode |= 0x01;
147 0 : if (( p.x() + p.y() - p.z()) > 1.5) outcode |= 0x02;
148 0 : if (( p.x() - p.y() + p.z()) > 1.5) outcode |= 0x04;
149 0 : if (( p.x() - p.y() - p.z()) > 1.5) outcode |= 0x08;
150 0 : if ((-p.x() + p.y() + p.z()) > 1.5) outcode |= 0x10;
151 0 : if ((-p.x() + p.y() - p.z()) > 1.5) outcode |= 0x20;
152 0 : if ((-p.x() - p.y() + p.z()) > 1.5) outcode |= 0x40;
153 0 : if ((-p.x() - p.y() - p.z()) > 1.5) outcode |= 0x80;
154 0 : return(outcode);
155 : }
156 :
157 : ///
158 : // CheckPoint()
159 : //
160 : // Test the point "alpha" of the way from P1 to P2
161 : // See if it is on a face of the cube
162 : // Consider only faces in "mask"
163 : static
164 0 : int CheckPoint(const vector3& p1, const vector3& p2, number alpha, long mask)
165 : {
166 : vector3 plane_point;
167 :
168 0 : plane_point.x() = LERP(alpha, p1.x(), p2.x());
169 0 : plane_point.y() = LERP(alpha, p1.y(), p2.y());
170 0 : plane_point.z() = LERP(alpha, p1.z(), p2.z());
171 0 : return(FacePlane(plane_point) & mask);
172 : }
173 :
174 : ///
175 : // CheckLine()
176 : //
177 : // Compute intersection of P1 --> P2 line segment with face planes
178 : // Then test intersection point to see if it is on cube face
179 : // Consider only face planes in "outcode_diff"
180 : // Note: Zero bits in "outcode_diff" means face line is outside of */
181 : //
182 : static
183 0 : int CheckLine(const vector3& p1, const vector3& p2, int outcode_diff)
184 : {
185 :
186 0 : if ((0x01 & outcode_diff) != 0)
187 0 : if (CheckPoint(p1,p2,( .5f-p1.x())/(p2.x()-p1.x()),0x3e) == INSIDE) return(INSIDE);
188 0 : if ((0x02 & outcode_diff) != 0)
189 0 : if (CheckPoint(p1,p2,(-.5f-p1.x())/(p2.x()-p1.x()),0x3d) == INSIDE) return(INSIDE);
190 0 : if ((0x04 & outcode_diff) != 0)
191 0 : if (CheckPoint(p1,p2,( .5f-p1.y())/(p2.y()-p1.y()),0x3b) == INSIDE) return(INSIDE);
192 0 : if ((0x08 & outcode_diff) != 0)
193 0 : if (CheckPoint(p1,p2,(-.5f-p1.y())/(p2.y()-p1.y()),0x37) == INSIDE) return(INSIDE);
194 0 : if ((0x10 & outcode_diff) != 0)
195 0 : if (CheckPoint(p1,p2,( .5f-p1.z())/(p2.z()-p1.z()),0x2f) == INSIDE) return(INSIDE);
196 0 : if ((0x20 & outcode_diff) != 0)
197 0 : if (CheckPoint(p1,p2,(-.5f-p1.z())/(p2.z()-p1.z()),0x1f) == INSIDE) return(INSIDE);
198 : return(OUTSIDE);
199 : }
200 :
201 : ///
202 : // PointTriangleIntersection()
203 : //
204 : // Test if 3D point is inside 3D triangle
205 : static
206 0 : int PointTriangleIntersection(const vector3& p, const TRI& t)
207 : {
208 : int sign12,sign23,sign31;
209 : vector3 vect12,vect23,vect31,vect1h,vect2h,vect3h;
210 : vector3 cross12_1p,cross23_2p,cross31_3p;
211 :
212 : ///
213 : // First, a quick bounding-box test:
214 : // If P is outside triangle bbox, there cannot be an intersection.
215 : //
216 0 : if (p.x() > MAX3(t.m_P[0].x(), t.m_P[1].x(), t.m_P[2].x())) return(OUTSIDE);
217 0 : if (p.y() > MAX3(t.m_P[0].y(), t.m_P[1].y(), t.m_P[2].y())) return(OUTSIDE);
218 0 : if (p.z() > MAX3(t.m_P[0].z(), t.m_P[1].z(), t.m_P[2].z())) return(OUTSIDE);
219 0 : if (p.x() < MIN3(t.m_P[0].x(), t.m_P[1].x(), t.m_P[2].x())) return(OUTSIDE);
220 0 : if (p.y() < MIN3(t.m_P[0].y(), t.m_P[1].y(), t.m_P[2].y())) return(OUTSIDE);
221 0 : if (p.z() < MIN3(t.m_P[0].z(), t.m_P[1].z(), t.m_P[2].z())) return(OUTSIDE);
222 :
223 : ///
224 : // For each triangle side, make a vector out of it by subtracting vertexes;
225 : // make another vector from one vertex to point P.
226 : // The crossproduct of these two vectors is orthogonal to both and the
227 : // signs of its X,Y,Z components indicate whether P was to the inside or
228 : // to the outside of this triangle side.
229 : //
230 0 : SUB(t.m_P[0], t.m_P[1], vect12);
231 0 : SUB(t.m_P[0], p, vect1h);
232 0 : CROSS(vect12, vect1h, cross12_1p)
233 0 : sign12 = SIGN3(cross12_1p); /* Extract X,Y,Z signs as 0..7 or 0...63 integer */
234 :
235 0 : SUB(t.m_P[1], t.m_P[2], vect23)
236 0 : SUB(t.m_P[1], p, vect2h);
237 0 : CROSS(vect23, vect2h, cross23_2p)
238 0 : sign23 = SIGN3(cross23_2p);
239 :
240 0 : SUB(t.m_P[2], t.m_P[0], vect31)
241 0 : SUB(t.m_P[2], p, vect3h);
242 0 : CROSS(vect31, vect3h, cross31_3p)
243 0 : sign31 = SIGN3(cross31_3p);
244 :
245 : ///
246 : // If all three crossproduct vectors agree in their component signs, /
247 : // then the point must be inside all three.
248 : // P cannot be OUTSIDE all three sides simultaneously.
249 : //
250 0 : return (((sign12 & sign23 & sign31) == 0) ? OUTSIDE : INSIDE);
251 : }
252 :
253 : ///
254 : // TriCubeIntersection()
255 : //
256 : // Triangle t is compared with a unit cube,
257 : // centered on the origin.
258 : // It returns INSIDE (0) or OUTSIDE(1) if t
259 : // Intersects or does not intersect the cube.
260 : static
261 0 : int TriCubeIntersection(const TRI& t)
262 : {
263 : int v1_test,v2_test,v3_test;
264 : number d;
265 : vector3 vect12,vect13,norm;
266 : vector3 hitpp,hitpn,hitnp,hitnn;
267 :
268 : ///
269 : // First compare all three vertexes with all six face-planes
270 : // If any vertex is inside the cube, return immediately!
271 : //
272 0 : if ((v1_test = FacePlane(t.m_P[0])) == INSIDE) return(INSIDE);
273 0 : if ((v2_test = FacePlane(t.m_P[1])) == INSIDE) return(INSIDE);
274 0 : if ((v3_test = FacePlane(t.m_P[2])) == INSIDE) return(INSIDE);
275 :
276 : ///
277 : // If all three vertexes were outside of one or more face-planes,
278 : // return immediately with a trivial rejection!
279 : //
280 0 : if ((v1_test & v2_test & v3_test) != 0) return(OUTSIDE);
281 :
282 : ///
283 : // Now do the same trivial rejection test for the 12 edge planes
284 : //
285 0 : v1_test |= Bevel2d(t.m_P[0]) << 8;
286 0 : v2_test |= Bevel2d(t.m_P[1]) << 8;
287 0 : v3_test |= Bevel2d(t.m_P[2]) << 8;
288 0 : if ((v1_test & v2_test & v3_test) != 0) return(OUTSIDE);
289 :
290 : ///
291 : // Now do the same trivial rejection test for the 8 corner planes
292 : //
293 0 : v1_test |= Bevel3d(t.m_P[0]) << 24;
294 0 : v2_test |= Bevel3d(t.m_P[1]) << 24;
295 0 : v3_test |= Bevel3d(t.m_P[2]) << 24;
296 0 : if ((v1_test & v2_test & v3_test) != 0) return(OUTSIDE);
297 :
298 : ///
299 : // If vertex 1 and 2, as a pair, cannot be trivially rejected
300 : // by the above tests, then see if the v1-->v2 triangle edge
301 : // intersects the cube. Do the same for v1-->v3 and v2-->v3.
302 : // Pass to the intersection algorithm the "OR" of the outcode
303 : // bits, so that only those cube faces which are spanned by
304 : // each triangle edge need be tested.
305 : //
306 0 : if ((v1_test & v2_test) == 0)
307 0 : if (CheckLine(t.m_P[0],t.m_P[1],v1_test|v2_test) == INSIDE) return(INSIDE);
308 0 : if ((v1_test & v3_test) == 0)
309 0 : if (CheckLine(t.m_P[0],t.m_P[2],v1_test|v3_test) == INSIDE) return(INSIDE);
310 0 : if ((v2_test & v3_test) == 0)
311 0 : if (CheckLine(t.m_P[1],t.m_P[2],v2_test|v3_test) == INSIDE) return(INSIDE);
312 :
313 : ///
314 : // By now, we know that the triangle is not off to any side,
315 : // and that its sides do not penetrate the cube. We must now
316 : // test for the cube intersecting the interior of the triangle.
317 : // We do this by looking for intersections between the cube
318 : // diagonals and the triangle...first finding the intersection
319 : // of the four diagonals with the plane of the triangle, and
320 : // then if that intersection is inside the cube, pursuing
321 : // whether the intersection point is inside the triangle itself.
322 :
323 : // To find plane of the triangle, first perform crossproduct on
324 : // two triangle side vectors to compute the normal vector.
325 0 : SUB(t.m_P[0],t.m_P[1],vect12);
326 0 : SUB(t.m_P[0],t.m_P[2],vect13);
327 0 : CROSS(vect12,vect13,norm)
328 :
329 : ///
330 : // The normal vector "norm" X,Y,Z components are the coefficients
331 : // of the triangles AX + BY + CZ + D = 0 plane equation. If we
332 : // solve the plane equation for X=Y=Z (a diagonal), we get
333 : // -D/(A+B+C) as a metric of the distance from cube center to the
334 : // diagonal/plane intersection. If this is between -0.5 and 0.5,
335 : // the intersection is inside the cube. If so, we continue by
336 : // doing a point/triangle intersection.
337 : // Do this for all four diagonals.
338 0 : d = norm.x() * t.m_P[0].x() + norm.y() * t.m_P[0].y() + norm.z() * t.m_P[0].z();
339 0 : hitpp.x() = hitpp.y() = hitpp.z() = d / (norm.x() + norm.y() + norm.z());
340 0 : if (fabs(hitpp.x()) <= 0.5)
341 0 : if (PointTriangleIntersection(hitpp,t) == INSIDE) return(INSIDE);
342 0 : hitpn.z() = -(hitpn.x() = hitpn.y() = d / (norm.x() + norm.y() - norm.z()));
343 0 : if (fabs(hitpn.x()) <= 0.5)
344 0 : if (PointTriangleIntersection(hitpn,t) == INSIDE) return(INSIDE);
345 0 : hitnp.y() = -(hitnp.x() = hitnp.z() = d / (norm.x() - norm.y() + norm.z()));
346 0 : if (fabs(hitnp.x()) <= 0.5)
347 0 : if (PointTriangleIntersection(hitnp,t) == INSIDE) return(INSIDE);
348 0 : hitnn.y() = hitnn.z() = -(hitnn.x() = d / (norm.x() - norm.y() - norm.z()));
349 0 : if (fabs(hitnn.x()) <= 0.5)
350 0 : if (PointTriangleIntersection(hitnn,t) == INSIDE) return(INSIDE);
351 :
352 : ///
353 : // No edge touched the cube; no cube diagonal touched the triangle.
354 : // We're done...there was no intersection.
355 : //
356 : return(OUTSIDE);
357 : }
358 :
359 : //////////////////////////////////////////////////////////////////////////////
360 : //
361 : // Public Functions
362 : //
363 : //
364 :
365 : ///
366 : // TriangleBoxIntersection()
367 : //
368 : // Determine if a bounding box and triangle intersect
369 : //
370 0 : bool TriangleBoxIntersection(const MathVector<3>& p0, const MathVector<3>& p1,
371 : const MathVector<3>& p2,
372 : const MathVector<3>& boxMin, const MathVector<3>& boxMax)
373 : {
374 : vector3 Trans;
375 : vector3 Scale(1.0, 1.0, 1.0);
376 : vector3 TransMax;
377 : TRI TestTri;
378 :
379 : ///
380 : // Compute the scale and transform required to make BBox
381 : // a voxel
382 : //
383 0 : Trans.x() = (boxMax.x() + boxMin.x()) / 2;
384 0 : Trans.y() = (boxMax.y() + boxMin.y()) / 2;
385 0 : Trans.z() = (boxMax.z() + boxMin.z()) / 2;
386 :
387 : VecSubtract(TransMax, boxMax, Trans);
388 :
389 0 : if(TransMax.x() != 0)
390 0 : Scale.x() = 0.5f / TransMax.x();
391 0 : if(TransMax.y() != 0)
392 0 : Scale.y() = 0.5f / TransMax.y();
393 0 : if(TransMax.z() != 0)
394 0 : Scale.z() = 0.5f / TransMax.z();
395 :
396 : ///
397 : // Put the triangle in voxel space
398 : //
399 0 : TestTri.m_P[0].x() = (p0.x() - Trans.x()) * Scale.x();
400 0 : TestTri.m_P[0].y() = (p0.y() - Trans.y()) * Scale.y();
401 0 : TestTri.m_P[0].z() = (p0.z() - Trans.z()) * Scale.z();
402 :
403 0 : TestTri.m_P[1].x() = (p1.x() - Trans.x()) * Scale.x();
404 0 : TestTri.m_P[1].y() = (p1.y() - Trans.y()) * Scale.y();
405 0 : TestTri.m_P[1].z() = (p1.z() - Trans.z()) * Scale.z();
406 :
407 0 : TestTri.m_P[2].x() = (p2.x() - Trans.x()) * Scale.x();
408 0 : TestTri.m_P[2].y() = (p2.y() - Trans.y()) * Scale.y();
409 0 : TestTri.m_P[2].z() = (p2.z() - Trans.z()) * Scale.z();
410 :
411 : ///
412 : // Test against the voxel
413 : //
414 0 : return(TriCubeIntersection(TestTri) == INSIDE);
415 : }
416 :
417 : }// end of namespace
|
