|
1 /************************************************************************* |
|
2 * * |
|
3 * Open Dynamics Engine, Copyright (C) 2001,2002 Russell L. Smith. * |
|
4 * All rights reserved. Email: russ@q12.org Web: www.q12.org * |
|
5 * * |
|
6 * This library is free software; you can redistribute it and/or * |
|
7 * modify it under the terms of EITHER: * |
|
8 * (1) The GNU Lesser General Public License as published by the Free * |
|
9 * Software Foundation; either version 2.1 of the License, or (at * |
|
10 * your option) any later version. The text of the GNU Lesser * |
|
11 * General Public License is included with this library in the * |
|
12 * file LICENSE.TXT. * |
|
13 * (2) The BSD-style license that is included with this library in * |
|
14 * the file LICENSE-BSD.TXT. * |
|
15 * * |
|
16 * This library is distributed in the hope that it will be useful, * |
|
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of * |
|
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the files * |
|
19 * LICENSE.TXT and LICENSE-BSD.TXT for more details. * |
|
20 * * |
|
21 *************************************************************************/ |
|
22 |
|
23 /* |
|
24 |
|
25 some useful collision utility stuff. this includes some API utility |
|
26 functions that are defined in the public header files. |
|
27 |
|
28 */ |
|
29 |
|
30 #include <ode/common.h> |
|
31 #include <ode/collision.h> |
|
32 #include <ode/odemath.h> |
|
33 #include "collision_util.h" |
|
34 |
|
35 //**************************************************************************** |
|
36 |
|
37 int dCollideSpheres (dVector3 p1, dReal r1, |
|
38 dVector3 p2, dReal r2, dContactGeom *c) |
|
39 { |
|
40 // printf ("d=%.2f (%.2f %.2f %.2f) (%.2f %.2f %.2f) r1=%.2f r2=%.2f\n", |
|
41 // d,p1[0],p1[1],p1[2],p2[0],p2[1],p2[2],r1,r2); |
|
42 |
|
43 dReal d = dDISTANCE (p1,p2); |
|
44 if (d > (r1 + r2)) return 0; |
|
45 if (d <= 0) { |
|
46 c->pos[0] = p1[0]; |
|
47 c->pos[1] = p1[1]; |
|
48 c->pos[2] = p1[2]; |
|
49 c->normal[0] = 1; |
|
50 c->normal[1] = 0; |
|
51 c->normal[2] = 0; |
|
52 c->depth = r1 + r2; |
|
53 } |
|
54 else { |
|
55 dReal d1 = dRecip (d); |
|
56 c->normal[0] = dMUL((p1[0]-p2[0]),d1); |
|
57 c->normal[1] = dMUL((p1[1]-p2[1]),d1); |
|
58 c->normal[2] = dMUL((p1[2]-p2[2]),d1); |
|
59 dReal k = dMUL(REAL(0.5),(r2 - r1 - d)); |
|
60 c->pos[0] = p1[0] + dMUL(c->normal[0],k); |
|
61 c->pos[1] = p1[1] + dMUL(c->normal[1],k); |
|
62 c->pos[2] = p1[2] + dMUL(c->normal[2],k); |
|
63 c->depth = r1 + r2 - d; |
|
64 } |
|
65 return 1; |
|
66 } |
|
67 |
|
68 |
|
69 void dLineClosestApproach (const dVector3 pa, const dVector3 ua, |
|
70 const dVector3 pb, const dVector3 ub, |
|
71 dReal *alpha, dReal *beta) |
|
72 { |
|
73 dVector3 p; |
|
74 p[0] = pb[0] - pa[0]; |
|
75 p[1] = pb[1] - pa[1]; |
|
76 p[2] = pb[2] - pa[2]; |
|
77 dReal uaub = dDOT(ua,ub); |
|
78 dReal q1 = dDOT(ua,p); |
|
79 dReal q2 = -dDOT(ub,p); |
|
80 dReal d = 1-dMUL(uaub,uaub); |
|
81 if (d <= REAL(0.0001)) { |
|
82 // @@@ this needs to be made more robust |
|
83 *alpha = 0; |
|
84 *beta = 0; |
|
85 } |
|
86 else { |
|
87 d = dRecip(d); |
|
88 *alpha = dMUL((q1 + dMUL(uaub,q2)),d); |
|
89 *beta = dMUL((dMUL(uaub,q1) + q2),d); |
|
90 } |
|
91 } |
|
92 |
|
93 |
|
94 // given two line segments A and B with endpoints a1-a2 and b1-b2, return the |
|
95 // points on A and B that are closest to each other (in cp1 and cp2). |
|
96 // in the case of parallel lines where there are multiple solutions, a |
|
97 // solution involving the endpoint of at least one line will be returned. |
|
98 // this will work correctly for zero length lines, e.g. if a1==a2 and/or |
|
99 // b1==b2. |
|
100 // |
|
101 // the algorithm works by applying the voronoi clipping rule to the features |
|
102 // of the line segments. the three features of each line segment are the two |
|
103 // endpoints and the line between them. the voronoi clipping rule states that, |
|
104 // for feature X on line A and feature Y on line B, the closest points PA and |
|
105 // PB between X and Y are globally the closest points if PA is in V(Y) and |
|
106 // PB is in V(X), where V(X) is the voronoi region of X. |
|
107 |
|
108 EXPORT_C void dClosestLineSegmentPoints (const dVector3 a1, const dVector3 a2, |
|
109 const dVector3 b1, const dVector3 b2, |
|
110 dVector3 cp1, dVector3 cp2) |
|
111 { |
|
112 dVector3 a1a2,b1b2,a1b1,a1b2,a2b1,a2b2,n; |
|
113 dReal la,lb,k,da1,da2,da3,da4,db1,db2,db3,db4,det; |
|
114 |
|
115 #define SET2(a,b) a[0]=b[0]; a[1]=b[1]; a[2]=b[2]; |
|
116 #define SET3(a,b,op,c) a[0]=b[0] op c[0]; a[1]=b[1] op c[1]; a[2]=b[2] op c[2]; |
|
117 |
|
118 // check vertex-vertex features |
|
119 |
|
120 SET3 (a1a2,a2,-,a1); |
|
121 SET3 (b1b2,b2,-,b1); |
|
122 SET3 (a1b1,b1,-,a1); |
|
123 da1 = dDOT(a1a2,a1b1); |
|
124 db1 = dDOT(b1b2,a1b1); |
|
125 if (da1 <= 0 && db1 >= 0) { |
|
126 SET2 (cp1,a1); |
|
127 SET2 (cp2,b1); |
|
128 return; |
|
129 } |
|
130 |
|
131 SET3 (a1b2,b2,-,a1); |
|
132 da2 = dDOT(a1a2,a1b2); |
|
133 db2 = dDOT(b1b2,a1b2); |
|
134 if (da2 <= 0 && db2 <= 0) { |
|
135 SET2 (cp1,a1); |
|
136 SET2 (cp2,b2); |
|
137 return; |
|
138 } |
|
139 |
|
140 SET3 (a2b1,b1,-,a2); |
|
141 da3 = dDOT(a1a2,a2b1); |
|
142 db3 = dDOT(b1b2,a2b1); |
|
143 if (da3 >= 0 && db3 >= 0) { |
|
144 SET2 (cp1,a2); |
|
145 SET2 (cp2,b1); |
|
146 return; |
|
147 } |
|
148 |
|
149 SET3 (a2b2,b2,-,a2); |
|
150 da4 = dDOT(a1a2,a2b2); |
|
151 db4 = dDOT(b1b2,a2b2); |
|
152 if (da4 >= 0 && db4 <= 0) { |
|
153 SET2 (cp1,a2); |
|
154 SET2 (cp2,b2); |
|
155 return; |
|
156 } |
|
157 |
|
158 // check edge-vertex features. |
|
159 // if one or both of the lines has zero length, we will never get to here, |
|
160 // so we do not have to worry about the following divisions by zero. |
|
161 |
|
162 la = dDOT(a1a2,a1a2); |
|
163 if (da1 >= 0 && da3 <= 0) { |
|
164 k = dDIV(da1,la); |
|
165 a1a2[0] = dMUL(k,a1a2[0]); |
|
166 a1a2[1] = dMUL(k,a1a2[1]); |
|
167 a1a2[2] = dMUL(k,a1a2[2]); |
|
168 SET3 (n,a1b1,-,a1a2); |
|
169 if (dDOT(b1b2,n) >= 0) { |
|
170 a1a2[0] = dMUL(k,a1a2[0]); |
|
171 a1a2[1] = dMUL(k,a1a2[1]); |
|
172 a1a2[2] = dMUL(k,a1a2[2]); |
|
173 SET3 (cp1,a1,+,a1a2); |
|
174 SET2 (cp2,b1); |
|
175 return; |
|
176 } |
|
177 } |
|
178 |
|
179 if (da2 >= 0 && da4 <= 0) { |
|
180 k = dDIV(da2,la); |
|
181 a1a2[0] = dMUL(k,a1a2[0]); |
|
182 a1a2[1] = dMUL(k,a1a2[1]); |
|
183 a1a2[2] = dMUL(k,a1a2[2]); |
|
184 SET3 (n,a1b2,-,a1a2); |
|
185 if (dDOT(b1b2,n) <= 0) { |
|
186 a1a2[0] = dMUL(k,a1a2[0]); |
|
187 a1a2[1] = dMUL(k,a1a2[1]); |
|
188 a1a2[2] = dMUL(k,a1a2[2]); |
|
189 SET3 (cp1,a1,+,a1a2); |
|
190 SET2 (cp2,b2); |
|
191 return; |
|
192 } |
|
193 } |
|
194 |
|
195 lb = dDOT(b1b2,b1b2); |
|
196 if (db1 <= 0 && db2 >= 0) { |
|
197 k = -dDIV(db1,lb); |
|
198 b1b2[0] = dMUL(k,b1b2[0]); |
|
199 b1b2[1] = dMUL(k,b1b2[1]); |
|
200 b1b2[2] = dMUL(k,b1b2[2]); |
|
201 SET3 (n,-a1b1,-,b1b2); |
|
202 if (dDOT(a1a2,n) >= 0) { |
|
203 b1b2[0] = dMUL(k,b1b2[0]); |
|
204 b1b2[1] = dMUL(k,b1b2[1]); |
|
205 b1b2[2] = dMUL(k,b1b2[2]); |
|
206 SET2 (cp1,a1); |
|
207 SET3 (cp2,b1,+,b1b2); |
|
208 return; |
|
209 } |
|
210 } |
|
211 |
|
212 if (db3 <= 0 && db4 >= 0) { |
|
213 k = -dDIV(db3,lb); |
|
214 b1b2[0] = dMUL(k,b1b2[0]); |
|
215 b1b2[1] = dMUL(k,b1b2[1]); |
|
216 b1b2[2] = dMUL(k,b1b2[2]); |
|
217 SET3 (n,-a2b1,-,b1b2); |
|
218 if (dDOT(a1a2,n) <= 0) { |
|
219 b1b2[0] = dMUL(k,b1b2[0]); |
|
220 b1b2[1] = dMUL(k,b1b2[1]); |
|
221 b1b2[2] = dMUL(k,b1b2[2]); |
|
222 SET2 (cp1,a2); |
|
223 SET3 (cp2,b1,+,b1b2); |
|
224 return; |
|
225 } |
|
226 } |
|
227 |
|
228 // it must be edge-edge |
|
229 |
|
230 k = dDOT(a1a2,b1b2); |
|
231 det = dMUL(la,lb) - dMUL(k,k); |
|
232 if (det <= 0) { |
|
233 // this should never happen, but just in case... |
|
234 SET2(cp1,a1); |
|
235 SET2(cp2,b1); |
|
236 return; |
|
237 } |
|
238 det = dRecip (det); |
|
239 dReal alpha = dMUL((dMUL(lb,da1) - dMUL(k,db1)),det); |
|
240 a1a2[0] = dMUL(alpha,a1a2[0]); |
|
241 a1a2[1] = dMUL(alpha,a1a2[1]); |
|
242 a1a2[2] = dMUL(alpha,a1a2[2]); |
|
243 dReal beta = dMUL(( dMUL(k,da1) - dMUL(la,db1)),det); |
|
244 b1b2[0] = dMUL(beta,b1b2[0]); |
|
245 b1b2[1] = dMUL(beta,b1b2[1]); |
|
246 b1b2[2] = dMUL(beta,b1b2[2]); |
|
247 SET3 (cp1,a1,+,a1a2); |
|
248 SET3 (cp2,b1,+,b1b2); |
|
249 |
|
250 # undef SET2 |
|
251 # undef SET3 |
|
252 } |
|
253 |
|
254 |
|
255 // a simple root finding algorithm is used to find the value of 't' that |
|
256 // satisfies: |
|
257 // d|D(t)|^2/dt = 0 |
|
258 // where: |
|
259 // |D(t)| = |p(t)-b(t)| |
|
260 // where p(t) is a point on the line parameterized by t: |
|
261 // p(t) = p1 + t*(p2-p1) |
|
262 // and b(t) is that same point clipped to the boundary of the box. in box- |
|
263 // relative coordinates d|D(t)|^2/dt is the sum of three x,y,z components |
|
264 // each of which looks like this: |
|
265 // |
|
266 // t_lo / |
|
267 // ______/ -->t |
|
268 // / t_hi |
|
269 // / |
|
270 // |
|
271 // t_lo and t_hi are the t values where the line passes through the planes |
|
272 // corresponding to the sides of the box. the algorithm computes d|D(t)|^2/dt |
|
273 // in a piecewise fashion from t=0 to t=1, stopping at the point where |
|
274 // d|D(t)|^2/dt crosses from negative to positive. |
|
275 |
|
276 void dClosestLineBoxPoints (const dVector3 p1, const dVector3 p2, |
|
277 const dVector3 c, const dMatrix3 R, |
|
278 const dVector3 side, |
|
279 dVector3 lret, dVector3 bret) |
|
280 { |
|
281 int i; |
|
282 |
|
283 // compute the start and delta of the line p1-p2 relative to the box. |
|
284 // we will do all subsequent computations in this box-relative coordinate |
|
285 // system. we have to do a translation and rotation for each point. |
|
286 dVector3 tmp,s,v; |
|
287 tmp[0] = p1[0] - c[0]; |
|
288 tmp[1] = p1[1] - c[1]; |
|
289 tmp[2] = p1[2] - c[2]; |
|
290 dMULTIPLY1_331 (s,R,tmp); |
|
291 tmp[0] = p2[0] - p1[0]; |
|
292 tmp[1] = p2[1] - p1[1]; |
|
293 tmp[2] = p2[2] - p1[2]; |
|
294 dMULTIPLY1_331 (v,R,tmp); |
|
295 |
|
296 // mirror the line so that v has all components >= 0 |
|
297 dVector3 sign; |
|
298 for (i=0; i<3; i++) { |
|
299 if (v[i] < 0) { |
|
300 s[i] = -s[i]; |
|
301 v[i] = -v[i]; |
|
302 sign[i] = REAL(-1.0); |
|
303 } |
|
304 else sign[i] = REAL(1.0); |
|
305 } |
|
306 |
|
307 // compute v^2 |
|
308 dVector3 v2; |
|
309 v2[0] = dMUL(v[0],v[0]); |
|
310 v2[1] = dMUL(v[1],v[1]); |
|
311 v2[2] = dMUL(v[2],v[2]); |
|
312 |
|
313 // compute the half-sides of the box |
|
314 dReal h[3]; |
|
315 h[0] = dMUL(REAL(0.5),side[0]); |
|
316 h[1] = dMUL(REAL(0.5),side[1]); |
|
317 h[2] = dMUL(REAL(0.5),side[2]); |
|
318 |
|
319 // region is -1,0,+1 depending on which side of the box planes each |
|
320 // coordinate is on. tanchor is the next t value at which there is a |
|
321 // transition, or the last one if there are no more. |
|
322 int region[3]; |
|
323 dReal tanchor[3]; |
|
324 |
|
325 // Denormals are a problem, because we divide by v[i], and then |
|
326 // multiply that by 0. Alas, infinity times 0 is infinity (!) |
|
327 // We also use v2[i], which is v[i] squared. Here's how the epsilons |
|
328 // are chosen: |
|
329 // float epsilon = 1.175494e-038 (smallest non-denormal number) |
|
330 // double epsilon = 2.225074e-308 (smallest non-denormal number) |
|
331 // For single precision, choose an epsilon such that v[i] squared is |
|
332 // not a denormal; this is for performance. |
|
333 // For double precision, choose an epsilon such that v[i] is not a |
|
334 // denormal; this is for correctness. (Jon Watte on mailinglist) |
|
335 |
|
336 const dReal tanchor_eps = REAL(2e-5f); |
|
337 |
|
338 // find the region and tanchor values for p1 |
|
339 for (i=0; i<3; i++) { |
|
340 if (v[i] > tanchor_eps) { |
|
341 if (s[i] < -h[i]) { |
|
342 region[i] = -1; |
|
343 tanchor[i] = dDIV((-h[i]-s[i]),v[i]); |
|
344 } |
|
345 else { |
|
346 region[i] = (s[i] > h[i]); |
|
347 tanchor[i] = dDIV((h[i]-s[i]),v[i]); |
|
348 } |
|
349 } |
|
350 else { |
|
351 region[i] = 0; |
|
352 tanchor[i] = REAL(2.0); // this will never be a valid tanchor |
|
353 } |
|
354 } |
|
355 |
|
356 // compute d|d|^2/dt for t=0. if it's >= 0 then p1 is the closest point |
|
357 dReal t = 0; |
|
358 dReal dd2dt = 0; |
|
359 for (i=0; i<3; i++) dd2dt -= dMUL((region[i] ? v2[i] : 0),tanchor[i]); |
|
360 if (dd2dt >= 0) goto got_answer; |
|
361 |
|
362 do { |
|
363 // find the point on the line that is at the next clip plane boundary |
|
364 dReal next_t = REAL(1.0); |
|
365 for (i=0; i<3; i++) { |
|
366 if (tanchor[i] > t && tanchor[i] < 1 && tanchor[i] < next_t) |
|
367 next_t = tanchor[i]; |
|
368 } |
|
369 |
|
370 // compute d|d|^2/dt for the next t |
|
371 dReal next_dd2dt = REAL(0.0); |
|
372 for (i=0; i<3; i++) { |
|
373 next_dd2dt += dMUL((region[i] ? v2[i] : 0),(next_t - tanchor[i])); |
|
374 } |
|
375 |
|
376 // if the sign of d|d|^2/dt has changed, solution = the crossover point |
|
377 if (next_dd2dt >= 0) { |
|
378 dReal m = dDIV((next_dd2dt-dd2dt),(next_t - t)); |
|
379 t -= dDIV(dd2dt,m); |
|
380 goto got_answer; |
|
381 } |
|
382 |
|
383 // advance to the next anchor point / region |
|
384 for (i=0; i<3; i++) { |
|
385 if (tanchor[i] == next_t) { |
|
386 tanchor[i] = dDIV((h[i]-s[i]),v[i]); |
|
387 region[i]++; |
|
388 } |
|
389 } |
|
390 t = next_t; |
|
391 dd2dt = next_dd2dt; |
|
392 } |
|
393 while (t < REAL(1.0)); |
|
394 t = REAL(1.0); |
|
395 |
|
396 got_answer: |
|
397 |
|
398 // compute closest point on the line |
|
399 for (i=0; i<3; i++) lret[i] = p1[i] + dMUL(t,tmp[i]); // note: tmp=p2-p1 |
|
400 |
|
401 // compute closest point on the box |
|
402 for (i=0; i<3; i++) { |
|
403 tmp[i] = dMUL(sign[i],(s[i] + dMUL(t,v[i]))); |
|
404 if (tmp[i] < -h[i]) tmp[i] = -h[i]; |
|
405 else if (tmp[i] > h[i]) tmp[i] = h[i]; |
|
406 } |
|
407 dMULTIPLY0_331 (s,R,tmp); |
|
408 for (i=0; i<3; i++) bret[i] = s[i] + c[i]; |
|
409 } |
|
410 |
|
411 |
|
412 // given boxes (p1,R1,side1) and (p1,R1,side1), return 1 if they intersect |
|
413 // or 0 if not. |
|
414 |
|
415 EXPORT_C int dBoxTouchesBox (const dVector3 p1, const dMatrix3 R1, |
|
416 const dVector3 side1, const dVector3 p2, |
|
417 const dMatrix3 R2, const dVector3 side2) |
|
418 { |
|
419 // two boxes are disjoint if (and only if) there is a separating axis |
|
420 // perpendicular to a face from one box or perpendicular to an edge from |
|
421 // either box. the following tests are derived from: |
|
422 // "OBB Tree: A Hierarchical Structure for Rapid Interference Detection", |
|
423 // S.Gottschalk, M.C.Lin, D.Manocha., Proc of ACM Siggraph 1996. |
|
424 |
|
425 // Rij is R1'*R2, i.e. the relative rotation between R1 and R2. |
|
426 // Qij is abs(Rij) |
|
427 dVector3 p,pp; |
|
428 dReal A1,A2,A3,B1,B2,B3,R11,R12,R13,R21,R22,R23,R31,R32,R33, |
|
429 Q11,Q12,Q13,Q21,Q22,Q23,Q31,Q32,Q33; |
|
430 |
|
431 // get vector from centers of box 1 to box 2, relative to box 1 |
|
432 p[0] = p2[0] - p1[0]; |
|
433 p[1] = p2[1] - p1[1]; |
|
434 p[2] = p2[2] - p1[2]; |
|
435 dMULTIPLY1_331 (pp,R1,p); // get pp = p relative to body 1 |
|
436 |
|
437 // get side lengths / 2 |
|
438 A1 = dMUL(side1[0],REAL(0.5)); A2 = dMUL(side1[1],REAL(0.5)); A3 = dMUL(side1[2],REAL(0.5)); |
|
439 B1 = dMUL(side2[0],REAL(0.5)); B2 = dMUL(side2[1],REAL(0.5)); B3 = dMUL(side2[2],REAL(0.5)); |
|
440 |
|
441 // for the following tests, excluding computation of Rij, in the worst case, |
|
442 // 15 compares, 60 adds, 81 multiplies, and 24 absolutes. |
|
443 // notation: R1=[u1 u2 u3], R2=[v1 v2 v3] |
|
444 |
|
445 // separating axis = u1,u2,u3 |
|
446 R11 = dDOT44(R1+0,R2+0); R12 = dDOT44(R1+0,R2+1); R13 = dDOT44(R1+0,R2+2); |
|
447 Q11 = dFabs(R11); Q12 = dFabs(R12); Q13 = dFabs(R13); |
|
448 if (dFabs(pp[0]) > (A1 + dMUL(B1,Q11) + dMUL(B2,Q12) + dMUL(B3,Q13))) return 0; |
|
449 R21 = dDOT44(R1+1,R2+0); R22 = dDOT44(R1+1,R2+1); R23 = dDOT44(R1+1,R2+2); |
|
450 Q21 = dFabs(R21); Q22 = dFabs(R22); Q23 = dFabs(R23); |
|
451 if (dFabs(pp[1]) > (A2 + dMUL(B1,Q21) + dMUL(B2,Q22) + dMUL(B3,Q23))) return 0; |
|
452 R31 = dDOT44(R1+2,R2+0); R32 = dDOT44(R1+2,R2+1); R33 = dDOT44(R1+2,R2+2); |
|
453 Q31 = dFabs(R31); Q32 = dFabs(R32); Q33 = dFabs(R33); |
|
454 if (dFabs(pp[2]) > (A3 + dMUL(B1,Q31) + dMUL(B2,Q32) + dMUL(B3,Q33))) return 0; |
|
455 |
|
456 // separating axis = v1,v2,v3 |
|
457 if (dFabs(dDOT41(R2+0,p)) > (dMUL(A1,Q11) + dMUL(A2,Q21) + dMUL(A3,Q31) + B1)) return 0; |
|
458 if (dFabs(dDOT41(R2+1,p)) > (dMUL(A1,Q12) + dMUL(A2,Q22) + dMUL(A3,Q32) + B2)) return 0; |
|
459 if (dFabs(dDOT41(R2+2,p)) > (dMUL(A1,Q13) + dMUL(A2,Q23) + dMUL(A3,Q33) + B3)) return 0; |
|
460 |
|
461 // separating axis = u1 x (v1,v2,v3) |
|
462 if (dFabs(dMUL(pp[2],R21)-dMUL(pp[1],R31)) > dMUL(A2,Q31) + dMUL(A3,Q21) + dMUL(B2,Q13) + dMUL(B3,Q12)) return 0; |
|
463 if (dFabs(dMUL(pp[2],R22)-dMUL(pp[1],R32)) > dMUL(A2,Q32) + dMUL(A3,Q22) + dMUL(B1,Q13) + dMUL(B3,Q11)) return 0; |
|
464 if (dFabs(dMUL(pp[2],R23)-dMUL(pp[1],R33)) > dMUL(A2,Q33) + dMUL(A3,Q23) + dMUL(B1,Q12) + dMUL(B2,Q11)) return 0; |
|
465 |
|
466 // separating axis = u2 x (v1,v2,v3) |
|
467 if (dFabs(dMUL(pp[0],R31)-dMUL(pp[2],R11)) > dMUL(A1,Q31) + dMUL(A3,Q11) + dMUL(B2,Q23) + dMUL(B3,Q22)) return 0; |
|
468 if (dFabs(dMUL(pp[0],R32)-dMUL(pp[2],R12)) > dMUL(A1,Q32) + dMUL(A3,Q12) + dMUL(B1,Q23) + dMUL(B3,Q21)) return 0; |
|
469 if (dFabs(pp[0]*R33-pp[2]*R13) > A1*Q33 + A3*Q13 + B1*Q22 + B2*Q21) return 0; |
|
470 |
|
471 // separating axis = u3 x (v1,v2,v3) |
|
472 if (dFabs(dMUL(pp[1],R11)-dMUL(pp[0],R21)) > dMUL(A1,Q21) + dMUL(A2,Q11) + dMUL(B2,Q33) + dMUL(B3,Q32)) return 0; |
|
473 if (dFabs(dMUL(pp[1],R12)-dMUL(pp[0],R22)) > dMUL(A1,Q22) + dMUL(A2,Q12) + dMUL(B1,Q33) + dMUL(B3,Q31)) return 0; |
|
474 if (dFabs(dMUL(pp[1],R13)-dMUL(pp[0],R23)) > dMUL(A1,Q23) + dMUL(A2,Q13) + dMUL(B1,Q32) + dMUL(B2,Q31)) return 0; |
|
475 |
|
476 return 1; |
|
477 } |
|
478 |
|
479 //**************************************************************************** |
|
480 // other utility functions |
|
481 |
|
482 EXPORT_C void dInfiniteAABB (dxGeom */*geom*/, dReal aabb[6]) |
|
483 { |
|
484 aabb[0] = -dInfinity; |
|
485 aabb[1] = dInfinity; |
|
486 aabb[2] = -dInfinity; |
|
487 aabb[3] = dInfinity; |
|
488 aabb[4] = -dInfinity; |
|
489 aabb[5] = dInfinity; |
|
490 } |
|
491 |
|
492 |
|
493 //**************************************************************************** |
|
494 // Helpers for Croteam's collider - by Nguyen Binh |
|
495 |
|
496 int dClipEdgeToPlane( dVector3 &vEpnt0, dVector3 &vEpnt1, const dVector4& plPlane) |
|
497 { |
|
498 // calculate distance of edge points to plane |
|
499 dReal fDistance0 = dPointPlaneDistance( vEpnt0 ,plPlane ); |
|
500 dReal fDistance1 = dPointPlaneDistance( vEpnt1 ,plPlane ); |
|
501 |
|
502 // if both points are behind the plane |
|
503 if ( fDistance0 < 0 && fDistance1 < 0 ) |
|
504 { |
|
505 // do nothing |
|
506 return 0; |
|
507 // if both points in front of the plane |
|
508 } |
|
509 else if ( fDistance0 > 0 && fDistance1 > 0 ) |
|
510 { |
|
511 // accept them |
|
512 return 1; |
|
513 // if we have edge/plane intersection |
|
514 } else if ((fDistance0 > 0 && fDistance1 < 0) || ( fDistance0 < 0 && fDistance1 > 0)) |
|
515 { |
|
516 |
|
517 // find intersection point of edge and plane |
|
518 dVector3 vIntersectionPoint; |
|
519 vIntersectionPoint[0]= vEpnt0[0]-dMUL((vEpnt0[0]-vEpnt1[0]),dDIV(fDistance0,(fDistance0-fDistance1))); |
|
520 vIntersectionPoint[1]= vEpnt0[1]-dMUL((vEpnt0[1]-vEpnt1[1]),dDIV(fDistance0,(fDistance0-fDistance1))); |
|
521 vIntersectionPoint[2]= vEpnt0[2]-dMUL((vEpnt0[2]-vEpnt1[2]),dDIV(fDistance0,(fDistance0-fDistance1))); |
|
522 |
|
523 // clamp correct edge to intersection point |
|
524 if ( fDistance0 < REAL(0.0) ) |
|
525 { |
|
526 dVector3Copy(vIntersectionPoint,vEpnt0); |
|
527 } else |
|
528 { |
|
529 dVector3Copy(vIntersectionPoint,vEpnt1); |
|
530 } |
|
531 return 1; |
|
532 } |
|
533 return 1; |
|
534 } |
|
535 |
|
536 // clip polygon with plane and generate new polygon points |
|
537 void dClipPolyToPlane( const dVector3 avArrayIn[], const int ctIn, |
|
538 dVector3 avArrayOut[], int &ctOut, |
|
539 const dVector4 &plPlane ) |
|
540 { |
|
541 // start with no output points |
|
542 ctOut = 0; |
|
543 |
|
544 int i0 = ctIn-1; |
|
545 |
|
546 // for each edge in input polygon |
|
547 for (int i1=0; i1<ctIn; i0=i1, i1++) { |
|
548 |
|
549 |
|
550 // calculate distance of edge points to plane |
|
551 dReal fDistance0 = dPointPlaneDistance( avArrayIn[i0],plPlane ); |
|
552 dReal fDistance1 = dPointPlaneDistance( avArrayIn[i1],plPlane ); |
|
553 |
|
554 // if first point is in front of plane |
|
555 if( fDistance0 >= 0 ) { |
|
556 // emit point |
|
557 avArrayOut[ctOut][0] = avArrayIn[i0][0]; |
|
558 avArrayOut[ctOut][1] = avArrayIn[i0][1]; |
|
559 avArrayOut[ctOut][2] = avArrayIn[i0][2]; |
|
560 ctOut++; |
|
561 } |
|
562 |
|
563 // if points are on different sides |
|
564 if( (fDistance0 > 0 && fDistance1 < 0) || ( fDistance0 < 0 && fDistance1 > 0) ) { |
|
565 |
|
566 // find intersection point of edge and plane |
|
567 dVector3 vIntersectionPoint; |
|
568 vIntersectionPoint[0]= avArrayIn[i0][0] - |
|
569 dMUL((avArrayIn[i0][0]-avArrayIn[i1][0]),dDIV(fDistance0,(fDistance0-fDistance1))); |
|
570 vIntersectionPoint[1]= avArrayIn[i0][1] - |
|
571 dMUL((avArrayIn[i0][1]-avArrayIn[i1][1]),dDIV(fDistance0,(fDistance0-fDistance1))); |
|
572 vIntersectionPoint[2]= avArrayIn[i0][2] - |
|
573 dMUL((avArrayIn[i0][2]-avArrayIn[i1][2]),dDIV(fDistance0,(fDistance0-fDistance1))); |
|
574 |
|
575 // emit intersection point |
|
576 avArrayOut[ctOut][0] = vIntersectionPoint[0]; |
|
577 avArrayOut[ctOut][1] = vIntersectionPoint[1]; |
|
578 avArrayOut[ctOut][2] = vIntersectionPoint[2]; |
|
579 ctOut++; |
|
580 } |
|
581 } |
|
582 |
|
583 } |
|
584 |
|
585 void dClipPolyToCircle(const dVector3 avArrayIn[], const int ctIn, |
|
586 dVector3 avArrayOut[], int &ctOut, |
|
587 const dVector4 &plPlane ,dReal fRadius) |
|
588 { |
|
589 // start with no output points |
|
590 ctOut = 0; |
|
591 |
|
592 int i0 = ctIn-1; |
|
593 |
|
594 // for each edge in input polygon |
|
595 for (int i1=0; i1<ctIn; i0=i1, i1++) |
|
596 { |
|
597 // calculate distance of edge points to plane |
|
598 dReal fDistance0 = dPointPlaneDistance( avArrayIn[i0],plPlane ); |
|
599 dReal fDistance1 = dPointPlaneDistance( avArrayIn[i1],plPlane ); |
|
600 |
|
601 // if first point is in front of plane |
|
602 if( fDistance0 >= 0 ) |
|
603 { |
|
604 // emit point |
|
605 if (dVector3Length2(avArrayIn[i0]) <= dMUL(fRadius,fRadius)) |
|
606 { |
|
607 avArrayOut[ctOut][0] = avArrayIn[i0][0]; |
|
608 avArrayOut[ctOut][1] = avArrayIn[i0][1]; |
|
609 avArrayOut[ctOut][2] = avArrayIn[i0][2]; |
|
610 ctOut++; |
|
611 } |
|
612 } |
|
613 |
|
614 // if points are on different sides |
|
615 if( (fDistance0 > 0 && fDistance1 < 0) || ( fDistance0 < 0 && fDistance1 > 0) ) |
|
616 { |
|
617 |
|
618 // find intersection point of edge and plane |
|
619 dVector3 vIntersectionPoint; |
|
620 vIntersectionPoint[0]= avArrayIn[i0][0] - |
|
621 dMUL((avArrayIn[i0][0]-avArrayIn[i1][0]),dDIV(fDistance0,(fDistance0-fDistance1))); |
|
622 vIntersectionPoint[1]= avArrayIn[i0][1] - |
|
623 dMUL((avArrayIn[i0][1]-avArrayIn[i1][1]),dDIV(fDistance0,(fDistance0-fDistance1))); |
|
624 vIntersectionPoint[2]= avArrayIn[i0][2] - |
|
625 dMUL((avArrayIn[i0][2]-avArrayIn[i1][2]),dDIV(fDistance0,(fDistance0-fDistance1))); |
|
626 |
|
627 // emit intersection point |
|
628 if (dVector3Length2(avArrayIn[i0]) <= dMUL(fRadius,fRadius)) |
|
629 { |
|
630 avArrayOut[ctOut][0] = vIntersectionPoint[0]; |
|
631 avArrayOut[ctOut][1] = vIntersectionPoint[1]; |
|
632 avArrayOut[ctOut][2] = vIntersectionPoint[2]; |
|
633 ctOut++; |
|
634 } |
|
635 } |
|
636 } |
|
637 } |
|
638 |