1
2 /*
3 A* -------------------------------------------------------------------
4 B* This file contains source code for the PyMOL computer program
5 C* copyright 1998-2000 by Warren Lyford Delano of DeLano Scientific.
6 D* -------------------------------------------------------------------
7 E* It is unlawful to modify or remove this copyright notice.
8 F* -------------------------------------------------------------------
9 G* Please see the accompanying LICENSE file for further information.
10 H* -------------------------------------------------------------------
11 I* Additional authors of this source file include:
12 -* Larry Coopet (various optimizations)
13 -*
14 -*
15 Z* -------------------------------------------------------------------
16 */
17
18 #include"os_predef.h"
19 #include"os_std.h"
20
21 #include"MemoryDebug.h"
22 #include"Base.h"
23 #include"Basis.h"
24 #include"Err.h"
25 #include"Feedback.h"
26 #include"Util.h"
27 #include"MemoryCache.h"
28 #include"Character.h"
29
30 static const float kR_SMALL4 = 0.0001F;
31 static const float kR_SMALL5 = 0.0001F;
32 #define EPSILON 0.000001F
33
34 /*========================================================================*/
ZLineToSphere(float * base,float * point,float * dir,float radius,float maxial,float * sphere,float * asum,float * pre)35 static int ZLineToSphere(float *base, float *point, float *dir, float radius,
36 float maxial, float *sphere, float *asum, float *pre)
37 {
38 /* Strategy - find an imaginary sphere that lies at the correct point on
39 the line segment, then treat as a sphere reflection */
40
41 float intra_p0, intra_p1, intra_p2;
42 float radial, axial, axial_sum, dangle, ab_dangle, axial_perp;
43 float radialsq, tan_acos_dangle;
44 float vradial0, vradial1, vradial2;
45 const float _0 = 0.0f;
46 float point0 = point[0];
47 float point1 = point[1];
48 float point2 = point[2];
49 float perpAxis0 = pre[0]; /* was cross_product(MinusZ,dir,perpAxis),normalize */
50 float perpAxis1 = pre[1];
51 float intra0 = point0 - base[0];
52 float intra1 = point1 - base[1];
53 const float dir0 = dir[0];
54 const float dir1 = dir[1];
55 const float dir2 = dir[2];
56 float dot, perpDist;
57
58 /* the perpAxis defines a perp-plane which includes the cyl-axis */
59
60 /* get minimum distance between the lines */
61
62 perpDist = intra0 * perpAxis0 + intra1 * perpAxis1;
63 /* was dot_product3f(intra,perpAxis); */
64
65 if((perpDist < _0 ? -perpDist : perpDist) > radius)
66 return 0;
67
68 /*if(fabs(perpDist) > radius) return 0; */
69
70 dangle = -dir2; /* was dot(MinusZ,dir) */
71 ab_dangle = (float) fabs(dangle);
72 if(ab_dangle > (1.0f - kR_SMALL4)) {
73 if(dangle > _0) {
74 sphere[0] = point0;
75 sphere[1] = point1;
76 sphere[2] = point2;
77 } else {
78 sphere[0] = dir0 * maxial + point0;
79 sphere[1] = dir1 * maxial + point1;
80 sphere[2] = dir2 * maxial + point2;
81 }
82 return (1);
83 }
84
85 if(ab_dangle > kR_SMALL4)
86 tan_acos_dangle = (float) (sqrt1d(1.0 - dangle * dangle) / dangle);
87 else
88 tan_acos_dangle = FLT_MAX;
89
90 /* now we need to define the triangle in the perp-plane
91 to figure out where the projected line intersection point is */
92
93 intra_p2 = point2 - base[2];
94
95 /* first, compute radial distance in the perp-plane between the two starting points */
96
97 dot = intra0 * perpAxis0 + intra1 * perpAxis1;
98 intra_p0 = intra0 - perpAxis0 * dot;
99 intra_p1 = intra1 - perpAxis1 * dot;
100
101 dot = intra_p0 * dir0 + intra_p1 * dir1 + intra_p2 * dir2;
102 vradial0 = intra_p0 - dir0 * dot;
103 vradial1 = intra_p1 - dir1 * dot;
104 vradial2 = intra_p2 - dir2 * dot;
105
106 radialsq = ((vradial0 * vradial0) + (vradial1 * vradial1) + (vradial2 * vradial2));
107
108 /* now figure out the axial distance along the cyl-line that will give us
109 the point of closest approach */
110
111 if(ab_dangle < kR_SMALL4)
112 axial_perp = _0;
113 else
114 axial_perp = (float) (sqrt1f(radialsq) / tan_acos_dangle);
115
116 axial =
117 (float) sqrt1f(((intra_p0 * intra_p0) + (intra_p1 * intra_p1) + (intra_p2 * intra_p2))
118 - radialsq);
119
120 if((intra_p0 * dir0 + intra_p1 * dir1 + intra_p2 * dir2) >= _0)
121 axial = axial_perp - axial;
122 else
123 axial = axial_perp + axial;
124
125 /* now we have to think about where the vector will actually strike the cylinder */
126
127 /* by definition, it must be perpdist away from the perp-plane becuase the perp-plane
128 is parallel to the line, so we can compute the radial component to this point */
129
130 radial = radius * radius - perpDist * perpDist;
131 radial = (float) sqrt1f(radial);
132
133 /* now the trick is figuring out how to adjust the axial distance to get the actual
134 position along the cyl line which will give us a representative sphere */
135
136 if(ab_dangle > kR_SMALL4)
137 axial_sum = axial - radial / tan_acos_dangle;
138 else
139 axial_sum = axial;
140 /*
141 printf("radial2 %8.3f \n",radial); */
142
143 if(axial_sum < _0)
144 axial_sum = _0;
145 else if(axial_sum > maxial)
146 axial_sum = maxial;
147
148 sphere[0] = dir0 * axial_sum + point0;
149 sphere[1] = dir1 * axial_sum + point1;
150 sphere[2] = dir2 * axial_sum + point2;
151
152 *asum = axial_sum;
153 /* printf("==>%8.3f sphere %8.3f %8.3f %8.3f\n",base[1],sphere[1],axial_perp,axial); */
154 return (1);
155 }
156
LineToSphere(float * base,float * ray,float * point,float * dir,float radius,float maxial,float * sphere,float * asum)157 static int LineToSphere(float *base, float *ray, float *point, float *dir, float radius,
158 float maxial, float *sphere, float *asum)
159 {
160 /* Strategy - find an imaginary sphere that lies at the correct point on
161 the line segment, then treat as a sphere reflection */
162
163 float perpDist, radial, axial, axial_sum, dangle, ab_dangle, axial_perp;
164 float radialsq, tan_acos_dangle;
165 float perpAxis0, perpAxis1, perpAxis2;
166 float intra0, intra1, intra2;
167 float intra_p0, intra_p1, intra_p2;
168 float vradial0, vradial1, vradial2;
169 float dir0 = dir[0], dir1 = dir[1], dir2 = dir[2];
170 float ray0 = ray[0], ray1 = ray[1], ray2 = ray[2];
171 float point0 = point[0], point1 = point[1], point2 = point[2];
172 float dot;
173
174 const float _0 = 0.0F;
175 const float _1 = 1.0F;
176
177 /* cross_product3f(ray,dir,perpAxis); */
178
179 perpAxis0 = ray1 * dir2 - ray2 * dir1;
180 perpAxis1 = ray2 * dir0 - ray0 * dir2;
181 perpAxis2 = ray0 * dir1 - ray1 * dir0;
182
183 /* subtract3f(point,base,intra); */
184 intra0 = point0 - base[0];
185
186 /* normalize3f(perpAxis) */
187 {
188 float len =
189 (float) sqrt1d((perpAxis0 * perpAxis0) + (perpAxis1 * perpAxis1) +
190 (perpAxis2 * perpAxis2));
191 intra1 = point1 - base[1];
192 /* subtract3f(point,base,intra); */
193 if(len > R_SMALL8) {
194 len = _1 / len;
195 intra2 = point2 - base[2];
196 perpAxis0 *= len;
197 perpAxis1 *= len;
198 perpAxis2 *= len;
199 } else {
200 intra2 = point2 - base[2];
201 }
202 }
203 /* the perpAxis defines a perp-plane which includes the cyl-axis */
204
205 /* get minimum distance between the lines */
206
207 /* dot_product3f(intra, perpAxis); */
208 perpDist = intra0 * perpAxis0;
209 perpDist += intra1 * perpAxis1 + intra2 * perpAxis2;
210
211 /*if(fabs(perpDist) > radius) return 0; */
212
213 dangle = ray0 * dir0;
214
215 if((perpDist < _0 ? -perpDist : perpDist) > radius)
216 return 0;
217
218 /* dangle = dot_product3f(ray, dir); */
219 dangle += ray1 * dir1 + ray2 * dir2;
220
221 ab_dangle = (float) fabs(dangle);
222
223 if(ab_dangle > (1.0f - kR_SMALL4)) {
224 if(dangle > _0) {
225 sphere[0] = point0;
226 sphere[1] = point1;
227 sphere[2] = point2;
228 } else {
229 sphere[0] = dir0 * maxial + point0;
230 sphere[1] = dir1 * maxial + point1;
231 sphere[2] = dir2 * maxial + point2;
232 }
233 return (1);
234 }
235
236 if(ab_dangle > kR_SMALL4)
237 tan_acos_dangle = (float) (sqrt1d(1.0 - dangle * dangle) / dangle);
238 else
239 tan_acos_dangle = FLT_MAX;
240
241 /* now we need to define the triangle in the perp-plane
242 to figure out where the projected line intersection point is */
243
244 /* first, compute radial distance in the perp-plane between the two starting points */
245
246 /* dot = dot_product3f(intra,perpAxis); */
247 dot = intra0 * perpAxis0 + intra1 * perpAxis1 + intra2 * perpAxis2;
248
249 intra_p0 = intra0 - perpAxis0 * dot;
250 intra_p1 = intra1 - perpAxis1 * dot;
251 intra_p2 = intra2 - perpAxis2 * dot;
252
253 /* dot = dot_product3f(intra_p, dir); */
254 dot = intra_p0 * dir0 + intra_p1 * dir1 + intra_p2 * dir2;
255
256 vradial0 = intra_p0 - dir0 * dot;
257 vradial1 = intra_p1 - dir1 * dot;
258 vradial2 = intra_p2 - dir2 * dot;
259
260 radialsq = ((vradial0 * vradial0) + (vradial1 * vradial1) + (vradial2 * vradial2));
261
262 /* now figure out the axial distance along the cyl-line that will give us
263 the point of closest approach */
264
265 if(ab_dangle < kR_SMALL4)
266 axial_perp = _0;
267 else
268 axial_perp = (float) (sqrt1f(radialsq) / tan_acos_dangle);
269
270 axial = (float) sqrt1f(((intra_p0 * intra_p0) +
271 (intra_p1 * intra_p1) + (intra_p2 * intra_p2)) - radialsq);
272
273 if((intra_p0 * dir0 + intra_p1 * dir1 + intra_p2 * dir2) >= _0)
274 axial = axial_perp - axial;
275 else
276 axial = axial_perp + axial;
277
278 /* now we have to think about where the vector will actually strike the cylinder */
279
280 /* by definition, it must be perpdist away from the perp-plane becuase the perp-plane
281 is parallel to the line, so we can compute the radial component to this point */
282
283 radial = radius * radius - perpDist * perpDist;
284 radial = (float) sqrt1f(radial);
285
286 /* now the trick is figuring out how to adjust the axial distance to get the actual
287 position along the cyl line which will give us a representative sphere */
288
289 if(ab_dangle > kR_SMALL4)
290 axial_sum = axial - radial / tan_acos_dangle;
291 else
292 axial_sum = axial;
293 /*
294 printf("radial2 %8.3f \n",radial); */
295
296 if(axial_sum < _0)
297 axial_sum = _0;
298 else if(axial_sum > maxial)
299 axial_sum = maxial;
300
301 sphere[0] = dir0 * axial_sum + point0;
302 sphere[1] = dir1 * axial_sum + point1;
303 sphere[2] = dir2 * axial_sum + point2;
304
305 *asum = axial_sum;
306 /* printf("==>%8.3f sphere %8.3f %8.3f %8.3f\n",base[1],sphere[1],axial_perp,axial); */
307 return (1);
308 }
309
FrontToInteriorSphere(float * front,float * point,float * dir,float radius,float radius2,float maxial)310 static int FrontToInteriorSphere(float *front,
311 float *point,
312 float *dir, float radius, float radius2, float maxial)
313 {
314 float intra_p[3];
315 float axial;
316 float intra[3], axis[3];
317 float sphere[3];
318
319 subtract3f(point, front, intra);
320 remove_component3f(intra, dir, intra_p);
321 add3f(front, intra_p, intra_p);
322 subtract3f(point, intra_p, axis);
323 axial = -dot_product3f(axis, dir);
324
325 if(axial < 0.0F)
326 axial = 0.0F;
327 else if(axial > maxial)
328 axial = maxial;
329
330 sphere[0] = axial * dir[0] + point[0];
331 sphere[1] = axial * dir[1] + point[1];
332 sphere[2] = axial * dir[2] + point[2];
333
334 return (diffsq3f(sphere, front) <= radius2);
335 }
336
337
338 /*========================================================================*/
ZLineToSphereCapped(float * base,float * point,float * dir,float radius,float maxial,float * sphere,float * asum,int cap1,int cap2,float * pre)339 static int ZLineToSphereCapped(float *base, float *point,
340 float *dir, float radius, float maxial,
341 float *sphere, float *asum, int cap1, int cap2, float *pre)
342 {
343 /* Strategy - find an imaginary sphere that lies at the correct point on
344 the line segment, then treat as a sphere reflection */
345
346 float perpAxis[3], intra_p[3];
347 float perpDist, radial, axial, axial_sum, dangle, ab_dangle, axial_perp;
348 float radialsq, tan_acos_dangle;
349 float len_proj;
350 float intra[3], vradial[3];
351 float diff[3], fpoint[3];
352 float proj[3];
353
354 perpAxis[0] = pre[0]; /* was cross_product(MinusZ,dir,perpAxis),normalize */
355 perpAxis[1] = pre[1];
356
357 /* the perpAxis defines a perp-plane which includes the cyl-axis */
358
359 intra[0] = point[0] - base[0];
360 intra[1] = point[1] - base[1];
361
362 /* get minimum distance between the lines */
363
364 perpDist = intra[0] * perpAxis[0] + intra[1] * perpAxis[1];
365 /* was dot_product3f(intra,perpAxis); */
366
367 if(fabs(perpDist) > radius) {
368 return (0);
369 }
370
371 perpAxis[2] = 0.0;
372 intra[2] = point[2] - base[2];
373
374 dangle = -dir[2]; /* was dot(MinusZ,dir) */
375 ab_dangle = (float) fabs(dangle);
376 if(ab_dangle > (1 - kR_SMALL4)) { /* vector inline with light ray... */
377 vradial[0] = point[0] - base[0];
378 vradial[1] = point[1] - base[1];
379 vradial[2] = 0.0F;
380 radial = (float) length3f(vradial);
381 if(radial > radius)
382 return 0;
383 if(dangle > 0.0) {
384 switch (cap1) {
385 case cCylCapFlat:
386 sphere[0] = base[0];
387 sphere[1] = base[1];
388 sphere[2] = point[2] - radius;
389 break;
390 case cCylCapRound:
391 sphere[0] = point[0];
392 sphere[1] = point[1];
393 sphere[2] = point[2];
394 }
395 return (1);
396 } else {
397 switch (cap1) {
398 case cCylCapFlat:
399 sphere[0] = base[0];
400 sphere[1] = base[1];
401 sphere[2] = dir[2] * maxial + point[2] - radius;
402 break;
403 case cCylCapRound:
404 sphere[0] = dir[0] * maxial + point[0];
405 sphere[1] = dir[1] * maxial + point[1];
406 sphere[2] = dir[2] * maxial + point[2];
407 }
408 return (1);
409 }
410 }
411
412 /*tan_acos_dangle = tan(acos(dangle)); */
413 tan_acos_dangle = (float) sqrt1f(1 - dangle * dangle) / dangle;
414
415 /*
416 printf("perpDist %8.3f\n",perpDist);
417 printf("dir %8.3f %8.3f %8.3f\n",dir[0],dir[1],dir[2]);
418 printf("base %8.3f %8.3f %8.3f\n",base[0],base[1],base[2]);
419 printf("point %8.3f %8.3f %8.3f\n",point[0],point[1],point[2]);
420 printf("intra %8.3f %8.3f %8.3f\n",intra[0],intra[1],intra[2]);
421 printf("perpAxis %8.3f %8.3f %8.3f\n",perpAxis[0],perpAxis[1],perpAxis[2]);
422 */
423
424 /* now we need to define the triangle in the perp-plane
425 to figure out where the projected line intersection point is */
426
427 /* first, compute radial distance in the perp-plane between the two starting points */
428
429 remove_component3f(intra, perpAxis, intra_p);
430 remove_component3f(intra_p, dir, vradial);
431
432 radialsq = lengthsq3f(vradial);
433
434 /* now figure out the axial distance along the cyl-line that will give us
435 the point of closest approach */
436
437 if(ab_dangle < kR_SMALL4)
438 axial_perp = 0;
439 else
440 axial_perp = (float) sqrt1f(radialsq) / tan_acos_dangle;
441
442 axial = (float) lengthsq3f(intra_p) - radialsq;
443 axial = (float) sqrt1f(axial);
444
445 /*
446 printf("radial %8.3f\n",radial);
447 printf("vradial %8.3f %8.3f %8.3f\n",vradial[0],vradial[1],vradial[2]);
448 printf("radial %8.3f\n",radial);
449 printf("dangle %8.3f \n",dangle);
450 printf("axial_perp %8.3f \n",axial_perp);
451 printf("axial1 %8.3f \n",axial);
452 printf("%8.3f\n",dot_product3f(intra_p,dir));
453 */
454
455 if(dot_product3f(intra_p, dir) >= 0.0)
456 axial = axial_perp - axial;
457 else
458 axial = axial_perp + axial;
459
460 /*
461 printf("axial2 %8.3f\n",axial);
462 */
463
464 /* now we have to think about where the vector will actually strike the cylinder */
465
466 /* by definition, it must be perpdist away from the perp-plane becuase the perp-plane
467 is parallel to the line, so we can compute the radial component to this point */
468
469 radial = radius * radius - perpDist * perpDist;
470 radial = (float) sqrt1f(radial);
471
472 /* now the trick is figuring out how to adjust the axial distance to get the actual
473 position along the cyl line which will give us a representative sphere */
474
475 if(ab_dangle > kR_SMALL4)
476 axial_sum = axial - radial / tan_acos_dangle;
477 else
478 axial_sum = axial;
479
480 /* printf("ab_dangle %8.3f \n",ab_dangle);
481
482 printf("axial_sum %8.3f \n",axial_sum);
483 */
484 if(axial_sum < 0) {
485 switch (cap1) {
486 case cCylCapFlat:
487 subtract3f(point, base, diff);
488 project3f(diff, dir, proj);
489 len_proj = (float) length3f(proj);
490 dangle = -proj[2] / len_proj;
491 if(fabs(dangle) < kR_SMALL4)
492 return 0;
493 sphere[0] = base[0];
494 sphere[1] = base[1];
495 sphere[2] = base[2] - len_proj / dangle;
496 if(diff3f(sphere, point) > radius)
497 return 0;
498 sphere[0] += dir[0] * radius;
499 sphere[1] += dir[1] * radius;
500 sphere[2] += dir[2] * radius;
501 *asum = 0;
502 break;
503 case cCylCapRound:
504 axial_sum = 0;
505 sphere[0] = point[0];
506 sphere[1] = point[1];
507 sphere[2] = point[2];
508 *asum = axial_sum;
509 break;
510 case cCylCapNone:
511 default:
512 return 0;
513 break;
514 }
515 } else if(axial_sum > maxial) {
516 switch (cap2) {
517
518 case cCylCapFlat:
519 scale3f(dir, maxial, fpoint);
520 add3f(fpoint, point, fpoint);
521 subtract3f(fpoint, base, diff);
522 project3f(diff, dir, proj);
523 len_proj = (float) length3f(proj);
524 dangle = -proj[2] / len_proj;
525 if(fabs(dangle) < kR_SMALL4)
526 return 0;
527 sphere[0] = base[0];
528 sphere[1] = base[1];
529 sphere[2] = base[2] - len_proj / dangle;
530 if(diff3f(sphere, fpoint) > radius)
531 return 0;
532 sphere[0] -= dir[0] * radius;
533 sphere[1] -= dir[1] * radius;
534 sphere[2] -= dir[2] * radius;
535 *asum = maxial;
536 break;
537 case cCylCapRound:
538 axial_sum = maxial;
539 sphere[0] = dir[0] * axial_sum + point[0];
540 sphere[1] = dir[1] * axial_sum + point[1];
541 sphere[2] = dir[2] * axial_sum + point[2];
542 *asum = axial_sum;
543 break;
544 case cCylCapNone:
545 default:
546 return 0;
547 break;
548 }
549 } else {
550 sphere[0] = dir[0] * axial_sum + point[0];
551 sphere[1] = dir[1] * axial_sum + point[1];
552 sphere[2] = dir[2] * axial_sum + point[2];
553
554 *asum = axial_sum;
555
556 /* printf("==>%8.3f sphere %8.3f %8.3f %8.3f\n",base[1],sphere[1],axial_perp,axial); */
557 }
558 return (1);
559 }
560
LineToSphereCapped(float * base,float * ray,float * point,float * dir,float radius,float maxial,float * sphere,float * asum,int cap1,int cap2)561 static int LineToSphereCapped(float *base, float *ray,
562 float *point, float *dir, float radius, float maxial,
563 float *sphere, float *asum, int cap1, int cap2)
564 {
565 /* Strategy - find an imaginary sphere that lies at the correct point on
566 the line segment, then treat as a sphere reflection */
567
568 float perpAxis[3], intra_p[3];
569 float perpDist, radial, axial, axial_sum, dangle, ab_dangle, axial_perp;
570 float radialsq, tan_acos_dangle;
571 float len_proj;
572 float intra[3], vradial[3];
573 float diff[3], fpoint[3];
574 float proj[3];
575
576 subtract3f(point, base, intra);
577
578 cross_product3f(ray, dir, perpAxis);
579
580 normalize3f(perpAxis);
581
582 /* the perpAxis defines a perp-plane which includes the cyl-axis */
583
584 /* get minimum distance between the lines */
585
586 perpDist = dot_product3f(intra, perpAxis);
587
588 /* was dot_product3f(intra,perpAxis); */
589
590 if(fabs(perpDist) > radius) {
591 return (0);
592 }
593
594 dangle = dot_product3f(ray, dir);
595 ab_dangle = (float) fabs(dangle);
596
597 if(ab_dangle > (1 - kR_SMALL4)) { /* vector inline with light ray... */
598 vradial[0] = point[0] - base[0];
599 vradial[1] = point[1] - base[1];
600 vradial[2] = point[2] - base[2];
601 radial = (float) length3f(vradial);
602 if(radial > radius)
603 return 0;
604 if(dangle > 0.0) {
605 switch (cap1) {
606 case cCylCapFlat:
607 sphere[0] = dir[0] * radius + point[0];
608 sphere[1] = dir[1] * radius + point[1];
609 sphere[2] = dir[2] * radius + point[2];
610 break;
611 case cCylCapRound:
612 sphere[0] = point[0];
613 sphere[1] = point[1];
614 sphere[2] = point[2];
615 }
616 return (1);
617 } else {
618 switch (cap1) {
619 case cCylCapFlat:
620 maxial -= radius;
621 break;
622 }
623 sphere[0] = dir[0] * maxial + point[0];
624 sphere[1] = dir[1] * maxial + point[1];
625 sphere[2] = dir[2] * maxial + point[2];
626 return (1);
627 }
628 }
629
630 /*tan_acos_dangle = tan(acos(dangle)); */
631 tan_acos_dangle = (float) sqrt1f(1 - dangle * dangle) / dangle;
632
633 /*
634 printf("perpDist %8.3f\n",perpDist);
635 printf("dir %8.3f %8.3f %8.3f\n",dir[0],dir[1],dir[2]);
636 printf("base %8.3f %8.3f %8.3f\n",base[0],base[1],base[2]);
637 printf("point %8.3f %8.3f %8.3f\n",point[0],point[1],point[2]);
638 printf("intra %8.3f %8.3f %8.3f\n",intra[0],intra[1],intra[2]);
639 printf("perpAxis %8.3f %8.3f %8.3f\n",perpAxis[0],perpAxis[1],perpAxis[2]);
640 */
641
642 /* now we need to define the triangle in the perp-plane
643 to figure out where the projected line intersection point is */
644
645 /* first, compute radial distance in the perp-plane between the two starting points */
646
647 remove_component3f(intra, perpAxis, intra_p);
648 remove_component3f(intra_p, dir, vradial);
649
650 radialsq = lengthsq3f(vradial);
651
652 /* now figure out the axial distance along the cyl-line that will give us
653 the point of closest approach */
654
655 if(ab_dangle < kR_SMALL4)
656 axial_perp = 0;
657 else
658 axial_perp = (float) sqrt1f(radialsq) / tan_acos_dangle;
659
660 axial = (float) lengthsq3f(intra_p) - radialsq;
661 axial = (float) sqrt1f(axial);
662
663 /*
664 printf("radial %8.3f\n",radial);
665 printf("vradial %8.3f %8.3f %8.3f\n",vradial[0],vradial[1],vradial[2]);
666 printf("radial %8.3f\n",radial);
667 printf("dangle %8.3f \n",dangle);
668 printf("axial_perp %8.3f \n",axial_perp);
669 printf("axial1 %8.3f \n",axial);
670 printf("%8.3f\n",dot_product3f(intra_p,dir));
671 */
672
673 if(dot_product3f(intra_p, dir) >= 0.0)
674 axial = axial_perp - axial;
675 else
676 axial = axial_perp + axial;
677
678 /*
679 printf("axial2 %8.3f\n",axial);
680 */
681
682 /* now we have to think about where the vector will actually strike the cylinder */
683
684 /* by definition, it must be perpdist away from the perp-plane becuase the perp-plane
685 is parallel to the line, so we can compute the radial component to this point */
686
687 radial = radius * radius - perpDist * perpDist;
688 radial = (float) sqrt1f(radial);
689
690 /* now the trick is figuring out how to adjust the axial distance to get the actual
691 position along the cyl line which will give us a representative sphere */
692
693 if(ab_dangle > kR_SMALL4)
694 axial_sum = axial - radial / tan_acos_dangle;
695 else
696 axial_sum = axial;
697
698 /* printf("ab_dangle %8.3f \n",ab_dangle);
699
700 printf("axial_sum %8.3f \n",axial_sum);
701 */
702 if(axial_sum < 0) {
703 switch (cap1) {
704 case cCylCapFlat:
705 subtract3f(point, base, diff);
706 project3f(diff, dir, proj);
707 len_proj = (float) length3f(proj);
708 dangle = dot_product3f(proj, ray) / len_proj;
709 if(fabs(dangle) < kR_SMALL4)
710 return 0;
711 len_proj /= dangle;
712 sphere[0] = base[0] + ray[0] * len_proj;
713 sphere[1] = base[1] + ray[1] * len_proj;
714 sphere[2] = base[2] + ray[2] * len_proj;
715 if(diff3f(sphere, point) > radius)
716 return 0;
717 sphere[0] += dir[0] * radius;
718 sphere[1] += dir[1] * radius;
719 sphere[2] += dir[2] * radius;
720 *asum = 0;
721 break;
722 case cCylCapRound:
723 axial_sum = 0;
724 sphere[0] = point[0];
725 sphere[1] = point[1];
726 sphere[2] = point[2];
727 *asum = axial_sum;
728 break;
729 case cCylCapNone:
730 default:
731 return 0;
732 break;
733 }
734 } else if(axial_sum > maxial) {
735 switch (cap2) {
736
737 case cCylCapFlat:
738 scale3f(dir, maxial, fpoint);
739 add3f(fpoint, point, fpoint);
740 subtract3f(fpoint, base, diff);
741 project3f(diff, dir, proj);
742 len_proj = (float) length3f(proj);
743 dangle = dot_product3f(proj, ray) / len_proj;
744 if(fabs(dangle) < kR_SMALL4)
745 return 0;
746 len_proj /= dangle;
747 sphere[0] = base[0] + ray[0] * len_proj;
748 sphere[1] = base[1] + ray[1] * len_proj;
749 sphere[2] = base[2] + ray[2] * len_proj;
750 if(diff3f(sphere, fpoint) > radius)
751 return 0;
752 sphere[0] -= dir[0] * radius;
753 sphere[1] -= dir[1] * radius;
754 sphere[2] -= dir[2] * radius;
755 *asum = maxial;
756 break;
757 case cCylCapRound:
758 axial_sum = maxial;
759 sphere[0] = dir[0] * axial_sum + point[0];
760 sphere[1] = dir[1] * axial_sum + point[1];
761 sphere[2] = dir[2] * axial_sum + point[2];
762 *asum = axial_sum;
763 break;
764 case cCylCapNone:
765 default:
766 return 0;
767 break;
768 }
769 } else {
770 sphere[0] = dir[0] * axial_sum + point[0];
771 sphere[1] = dir[1] * axial_sum + point[1];
772 sphere[2] = dir[2] * axial_sum + point[2];
773
774 *asum = axial_sum;
775
776 /* printf("==>%8.3f sphere %8.3f %8.3f %8.3f\n",base[1],sphere[1],axial_perp,axial); */
777 }
778 return (1);
779 }
780
ConeLineToSphereCapped(float * base,float * ray,float * point,float * dir,float radius,float small_radius,float maxial,float * sphere,float * asum,float * sph_rad,float * sph_rad_sq,int cap1,int cap2)781 static int ConeLineToSphereCapped(float *base, float *ray,
782 float *point, float *dir, float radius,
783 float small_radius, float maxial, float *sphere,
784 float *asum, float *sph_rad, float *sph_rad_sq,
785 int cap1, int cap2)
786 {
787 /* Strategy - find an imaginary sphere that lies at the correct point on
788 the line segment, then treat as a sphere reflection */
789
790 float axial_sum, dangle, ab_dangle;
791 float len_proj;
792 float diff[3], fpoint[3];
793 float proj[3];
794
795 {
796 float perp_axis[3];
797 float intra[3];
798 float perp_dist;
799
800 subtract3f(point, base, intra);
801 cross_product3f(ray, dir, perp_axis);
802
803 normalize3f(perp_axis);
804
805 /* the perp_axis defines a perp-plane which includes the cyl-axis */
806
807 /* get minimum distance between the lines */
808
809 perp_dist = fabs(dot_product3f(intra, perp_axis));
810
811 if(perp_dist > radius) {
812 /* the infinite ray and the cone direction lines don't pass close
813 enough to intersect within the bounding cylinder */
814 return 0;
815 }
816 }
817
818 dangle = dot_product3f(ray, dir);
819 ab_dangle = (float) fabs(dangle);
820
821 /* set up the cone */
822
823 {
824 double spread = (radius - small_radius) / maxial;
825 float orig_axial_len = radius / spread;
826 float base2orig_radial[3];
827 float base2orig_normal[3];
828 float base2orig_radial_len, base2orig_radial_len_sq;
829 float base2orig_len_sq;
830 float base2orig_axial_len, base2orig_spread;
831 float orig[3];
832 float base2orig[3];
833 float near_p[3];
834
835 int base_inside_cone = false;
836 float shift1, shift2; /* this is what we are solving for --
837 the distance along the cone axis to the point of intersection */
838
839 scale3f(dir, orig_axial_len, orig);
840 add3f(point, orig, orig);
841 subtract3f(orig, base, base2orig);
842
843 remove_component3f(base2orig, dir, base2orig_radial);
844
845 base2orig_radial_len_sq = lengthsq3f(base2orig_radial);
846
847 base2orig_len_sq = lengthsq3f(base2orig);
848
849 base2orig_axial_len = sqrt1f(base2orig_len_sq - base2orig_radial_len_sq);
850
851 base2orig_radial_len = sqrt1f(base2orig_radial_len_sq);
852
853 base2orig_spread = base2orig_radial_len / base2orig_axial_len;
854
855 base_inside_cone = (base2orig_spread < spread);
856
857 normalize23f(base2orig, base2orig_normal);
858
859 if(ab_dangle > kR_SMALL4) {
860 float ray_extend = base2orig_axial_len / dangle;
861 if(dot_product3f(base2orig_normal, dir) < 0.0)
862 ray_extend = -ray_extend;
863
864 scale3f(ray, ray_extend, near_p);
865 add3f(base, near_p, near_p);
866
867 /* Now we punt entirely and throw the solution of this quadratic
868 relationship over to Mathematica. Surely this calculation
869 could be significantly optimized... */
870
871 {
872 double partA, partB, partC;
873
874 double dir0 = dir[0], dir1 = dir[1], dir2 = dir[2];
875 double ray0 = ray[0], ray1 = ray[1], ray2 = ray[2];
876 double dir0Sq = dir0 * dir0, dir1Sq = dir1 * dir1, dir2Sq = dir2 * dir2;
877 double ray0Sq = ray0 * ray0, ray1Sq = ray1 * ray1, ray2Sq = ray2 * ray2;
878
879 double cone0 = orig[0], cone1 = orig[1], cone2 = orig[2];
880 double near0 = near_p[0], near1 = near_p[1], near2 = near_p[2];
881 double cone0Sq = cone0 * cone0, cone1Sq = cone1 * cone1, cone2Sq = cone2 * cone2;
882 double near0Sq = near0 * near0, near1Sq = near1 * near1, near2Sq = near2 * near2;
883
884 double dAngle = ray0 * dir0 + ray1 * dir1 + ray2 * dir2;
885
886 double spreadSq = spread * spread;
887 double dAngleSq = dAngle * dAngle;
888
889 partB = dAngleSq * (4 * pow(cone0 * dAngle * dir0 + cone1 * dAngle * dir1 +
890 cone2 * dAngle * dir2 - dAngle * dir0 * near0 -
891 dAngle * dir1 * near1 - dAngle * dir2 * near2 -
892 cone0 * ray0 + near0 * ray0 - cone1 * ray1 +
893 near1 * ray1 - cone2 * ray2 + near2 * ray2,
894 2.0) - 4 * (cone0Sq + cone1Sq + cone2Sq -
895 2 * cone0 * near0 + near0Sq -
896 2 * cone1 * near1 + near1Sq -
897 2 * cone2 * near2 + near2Sq) * (ray0Sq +
898 ray1Sq -
899 2 *
900 dAngle *
901 (dir0 *
902 ray0 +
903 dir1 *
904 ray1 +
905 dir2 *
906 ray2) +
907 ray2Sq +
908 dAngleSq *
909 (dir0Sq +
910 dir1Sq +
911 dir2Sq -
912 spreadSq)));
913
914 if(partB < 0.0) { /* negative? then there are NO real solutions */
915 return 0;
916 } else {
917 double partBroot = sqrt(partB);
918
919 partA = -(cone0 * dAngleSq * dir0) - cone1 * dAngleSq * dir1 -
920 cone2 * dAngleSq * dir2 + dAngleSq * dir0 * near0 + dAngleSq * dir1 * near1 +
921 dAngleSq * dir2 * near2 + cone0 * dAngle * ray0 - dAngle * near0 * ray0 +
922 cone1 * dAngle * ray1 - dAngle * near1 * ray1 + cone2 * dAngle * ray2 -
923 dAngle * near2 * ray2;
924
925 partC = (ray0Sq + ray1Sq -
926 2 * dAngle * (dir0 * ray0 + dir1 * ray1 + dir2 * ray2) + ray2Sq +
927 dAngleSq * (dir0Sq + dir1Sq + dir2Sq - spreadSq));
928
929 shift1 = (float) ((partA + partBroot * 0.5) / partC);
930 shift2 = (float) ((partA - partBroot * 0.5) / partC);
931 }
932 }
933
934 {
935 float axial_sum1 = orig_axial_len + shift1;
936 float axial_sum2 = orig_axial_len + shift2;
937
938 if(dangle > 0.0F) { /* cone is narrowing in parallel with ray */
939 if(shift1 < shift2) {
940 axial_sum = axial_sum1;
941 } else {
942 axial_sum = axial_sum2;
943 }
944 if((axial_sum < 0.0F) || (base_inside_cone && (axial_sum < orig_axial_len))) {
945 switch (cap1) {
946 case cCylCapFlat:
947 subtract3f(point, base, diff);
948 project3f(diff, dir, proj);
949 len_proj = (float) length3f(proj);
950 dangle = dot_product3f(proj, ray) / len_proj;
951 if(fabs(dangle) < kR_SMALL5)
952 return 0;
953 len_proj /= dangle;
954 sphere[0] = base[0] + ray[0] * len_proj;
955 sphere[1] = base[1] + ray[1] * len_proj;
956 sphere[2] = base[2] + ray[2] * len_proj;
957 if(diff3f(sphere, point) > radius)
958 return 0;
959 sphere[0] += dir[0] * radius;
960 sphere[1] += dir[1] * radius;
961 sphere[2] += dir[2] * radius;
962 *sph_rad = radius;
963 *sph_rad_sq = radius * radius;
964 *asum = 0;
965 return 1;
966 break;
967 case cCylCapNone:
968 default:
969 return 0;
970 break;
971 }
972 } else if(axial_sum > maxial) {
973 return 0;
974 }
975 } else { /* cone is narrowing against ray */
976 if(shift1 < shift2) {
977 axial_sum = axial_sum2;
978 if(axial_sum > orig_axial_len) {
979 axial_sum = axial_sum1;
980 }
981 } else {
982 axial_sum = axial_sum1;
983 if(axial_sum > orig_axial_len) {
984 axial_sum = axial_sum2;
985 }
986 }
987 if(axial_sum < 0.0F) {
988 return 0;
989 } else if(axial_sum >= maxial) {
990 switch (cap2) {
991 case cCylCapFlat:
992 scale3f(dir, maxial, fpoint);
993 add3f(fpoint, point, fpoint);
994 subtract3f(fpoint, base, diff);
995 project3f(diff, dir, proj);
996 len_proj = (float) length3f(proj);
997 dangle = dot_product3f(proj, ray) / len_proj;
998 if(fabs(dangle) < kR_SMALL5)
999 return 0;
1000 len_proj /= dangle;
1001 sphere[0] = base[0] + ray[0] * len_proj;
1002 sphere[1] = base[1] + ray[1] * len_proj;
1003 sphere[2] = base[2] + ray[2] * len_proj;
1004 if(diff3f(sphere, fpoint) > small_radius)
1005 /* need to handle this case */
1006 return 0;
1007 sphere[0] -= dir[0] * small_radius;
1008 sphere[1] -= dir[1] * small_radius;
1009 sphere[2] -= dir[2] * small_radius;
1010 *sph_rad = small_radius;
1011 *sph_rad_sq = small_radius * small_radius;
1012 *asum = maxial;
1013 return 1;
1014 break;
1015 case cCylCapNone:
1016 default:
1017 return 0;
1018 break;
1019 }
1020 }
1021 }
1022 }
1023 } else {
1024 axial_sum = orig_axial_len - base2orig_axial_len;
1025 if((axial_sum < 0.0F) || (axial_sum > maxial))
1026 return 0;
1027 }
1028
1029 /* normal hit in mid-section of the cone */
1030
1031 {
1032 float radius_at_hit = radius - spread * axial_sum;
1033
1034 float adjustment = radius_at_hit * spread;
1035
1036 *asum = axial_sum; /* color blend based on actual hit location */
1037
1038 axial_sum -= adjustment;
1039
1040 /* printf(" %8.3f ",axial_sum); */
1041
1042 sphere[0] = dir[0] * axial_sum + point[0];
1043 sphere[1] = dir[1] * axial_sum + point[1];
1044 sphere[2] = dir[2] * axial_sum + point[2];
1045
1046 /* and provide virtual sphere radius info */
1047
1048 *sph_rad_sq = (float) ((radius_at_hit * radius_at_hit) + (adjustment * adjustment));
1049 *sph_rad = (float) sqrt(*sph_rad_sq);
1050 /* printf("%8.3f %8.3f ",*sph_rad, tan_acos_dangle); */
1051
1052 /* dump3f(sphere,"sph"); */
1053 return 1;
1054 }
1055 }
1056 return 0;
1057 }
1058
FrontToInteriorSphereCapped(float * front,float * point,float * dir,float radius,float radius2,float maxial,int cap1,int cap2)1059 static int FrontToInteriorSphereCapped(float *front,
1060 float *point,
1061 float *dir,
1062 float radius,
1063 float radius2, float maxial, int cap1, int cap2)
1064 {
1065 float intra_p[3];
1066 float axial;
1067 float intra[3], axis[3];
1068 float sphere[3];
1069
1070 subtract3f(point, front, intra);
1071 remove_component3f(intra, dir, intra_p);
1072 add3f(front, intra_p, intra_p);
1073 subtract3f(point, intra_p, axis);
1074 axial = -dot_product3f(axis, dir);
1075
1076 if(axial < 0.0F)
1077 return 0;
1078 else if(axial > maxial)
1079 return 0;
1080
1081 sphere[0] = axial * dir[0] + point[0];
1082 sphere[1] = axial * dir[1] + point[1];
1083 sphere[2] = axial * dir[2] + point[2];
1084
1085 return (diffsq3f(sphere, front) < radius2);
1086
1087 }
1088
1089
1090 /*========================================================================*/
ZLineClipPoint(float * base,float * point,float * alongNormalSq,float cutoff)1091 static float ZLineClipPoint(float *base, float *point, float *alongNormalSq, float cutoff)
1092 {
1093 float hyp0, hyp1, hyp2;
1094 float result = FLT_MAX;
1095
1096 /* this routine determines whether or not a vector starting at "base"
1097 heading in the direction "normal" intersects a sphere located at "point".
1098
1099 It returns how far along the vector the intersection with the plane is or
1100 MAXFLOAT if there isn't a relevant intersection
1101
1102 NOTE: this routine has been optimized for normals along Z
1103 Optimizes-out vectors that are more than "cutoff" from "point" in x,y plane
1104 */
1105
1106 hyp0 = point[0] - base[0];
1107 if(fabs(hyp0) > cutoff)
1108 return result;
1109
1110 hyp1 = point[1] - base[1];
1111 if(fabs(hyp1) > cutoff)
1112 return result;
1113
1114 hyp2 = point[2] - base[2];
1115
1116 if(hyp2 < 0.0) {
1117 (*alongNormalSq) = (hyp2 * hyp2);
1118 result = (hyp0 * hyp0) + (hyp1 * hyp1);
1119 }
1120 return result;
1121 }
1122
ZLineClipPointNoZCheck(float * base,float * point,float * alongNormalSq,float cutoff)1123 static float ZLineClipPointNoZCheck(float *base, float *point, float *alongNormalSq,
1124 float cutoff)
1125 {
1126 float hyp0, hyp1, hyp2;
1127 float result = FLT_MAX;
1128
1129 /* this routine determines whether or not a vector starting at "base"
1130 heading in the direction "normal" intersects a sphere located at "point".
1131
1132 It returns how far along the vector the intersection with the plane is or
1133 MAXFLOAT if there isn't a relevant intersection
1134
1135 NOTE: this routine has been optimized for normals along Z
1136 Optimizes-out vectors that are more than "cutoff" from "point" in x,y plane
1137 */
1138
1139 hyp0 = point[0] - base[0];
1140 if(fabs(hyp0) > cutoff)
1141 return result;
1142
1143 hyp1 = point[1] - base[1];
1144 if(fabs(hyp1) > cutoff)
1145 return result;
1146
1147 hyp2 = point[2] - base[2];
1148
1149 (*alongNormalSq) = (hyp2 * hyp2);
1150 result = (hyp0 * hyp0) + (hyp1 * hyp1);
1151
1152 return result;
1153 }
1154
LineClipPoint(float * base,float * ray,float * point,float * dist,float cutoff,float cutoff2)1155 static int LineClipPoint(float *base, float *ray,
1156 float *point, float *dist, float cutoff, float cutoff2)
1157 {
1158 float hyp0, hyp1, hyp2;
1159 float opp0, opp1, opp2;
1160 float adj0, adj1, adj2;
1161 float ray0, ray1, ray2;
1162 float proj;
1163 double dcutoff = (double) cutoff;
1164 float opp_len_sq;
1165
1166 /* this routine determines whether or not a vector starting at "base"
1167 heading in the direction "ray" intersects a sphere located at "point".
1168
1169 It returns how far along the vector the intersection with the plane is or
1170 MAXFLOAT if there isn't a relevant intersection
1171
1172 NOTE: this routine has been optimized for normals along Z
1173 Optimizes-out vectors that are more than "cutoff" from "point"
1174 */
1175
1176 /* compute the hypo */
1177
1178 hyp2 = point[2] - base[2];
1179 hyp1 = point[1] - base[1];
1180 hyp0 = point[0] - base[0];
1181
1182 ray0 = ray[0];
1183 ray1 = ray[1];
1184 ray2 = ray[2];
1185
1186 /* compute the adjacent edge (dot-projection) */
1187
1188 proj = (ray0 * hyp0) + (ray1 * hyp1) + (ray2 * hyp2);
1189
1190 adj0 = ray0 * proj;
1191 adj1 = ray1 * proj;
1192 opp0 = hyp0 - adj0;
1193 adj2 = ray2 * proj;
1194 if(fabs(opp0) > dcutoff)
1195 return 0;
1196 opp1 = hyp1 - adj1;
1197 opp2 = hyp2 - adj2;
1198 if(fabs(opp1) > dcutoff)
1199 return 0;
1200 if(fabs(opp2) > dcutoff)
1201 return 0;
1202
1203 opp_len_sq = (opp0 * opp0) + (opp1 * opp1) + (opp2 * opp2);
1204
1205 if(opp_len_sq <= cutoff2) {
1206 *dist = proj - (float) sqrt1f(cutoff2 - opp_len_sq);
1207 return 1;
1208 }
1209
1210 return 0;
1211 }
1212
LineClipEllipsoidPoint(float * base,float * ray,float * point,float * dist,float cutoff,float cutoff2,float * scale,float * n1,float * n2,float * n3)1213 static int LineClipEllipsoidPoint(float *base, float *ray,
1214 float *point, float *dist,
1215 float cutoff, float cutoff2,
1216 float *scale, float *n1, float *n2, float *n3)
1217 {
1218 float point_to_base[3];
1219 float scaled_base[3];
1220 float scaled_ray[3];
1221 float d1, d2, d3, s1, s2, s3;
1222 float comp1[3], comp2[3], comp3[3];
1223
1224 subtract3f(base, point, point_to_base);
1225
1226 /* project difference vector onto ellipsoid axes */
1227
1228 d1 = dot_product3f(point_to_base, n1);
1229 d2 = dot_product3f(point_to_base, n2);
1230 d3 = dot_product3f(point_to_base, n3);
1231
1232 s1 = d1 / scale[0];
1233 s2 = d2 / scale[1];
1234 s3 = d3 / scale[2];
1235
1236 scale3f(n1, s1, comp1);
1237 scale3f(n2, s2, comp2);
1238 scale3f(n3, s3, comp3);
1239
1240 copy3f(point, scaled_base);
1241 add3f(comp1, scaled_base, scaled_base);
1242 add3f(comp2, scaled_base, scaled_base);
1243 add3f(comp3, scaled_base, scaled_base);
1244
1245 d1 = dot_product3f(ray, n1);
1246 d2 = dot_product3f(ray, n2);
1247 d3 = dot_product3f(ray, n3);
1248
1249 s1 = d1 / scale[0];
1250 s2 = d2 / scale[1];
1251 s3 = d3 / scale[2];
1252
1253 scale3f(n1, s1, comp1);
1254 scale3f(n2, s2, comp2);
1255 scale3f(n3, s3, comp3);
1256
1257 copy3f(comp1, scaled_ray);
1258 add3f(comp2, scaled_ray, scaled_ray);
1259 add3f(comp3, scaled_ray, scaled_ray);
1260
1261 d1 = length3f(scaled_ray);
1262
1263 normalize3f(scaled_ray);
1264
1265 /* dump3f(ray,"ray");
1266 dump3f(scaled_ray,"scaled_ray");
1267 */
1268
1269 {
1270 float hyp0, hyp1, hyp2;
1271 float opp0, opp1, opp2;
1272 float adj0, adj1, adj2;
1273 float ray0, ray1, ray2;
1274 float proj;
1275 double dcutoff = (double) cutoff;
1276 float opp_len_sq;
1277
1278 /* this routine determines whether or not a vector starting at "base"
1279 heading in the direction "ray" intersects a sphere located at "point".
1280
1281 It returns how far along the vector the intersection with the plane is or
1282 MAXFLOAT if there isn't a relevant intersection
1283
1284 */
1285
1286 /* compute the hypo */
1287
1288 hyp2 = point[2] - scaled_base[2];
1289 hyp1 = point[1] - scaled_base[1];
1290 hyp0 = point[0] - scaled_base[0];
1291
1292 ray0 = scaled_ray[0];
1293 ray1 = scaled_ray[1];
1294 ray2 = scaled_ray[2];
1295
1296 /* compute the adjacent edge (dot-projection) */
1297
1298 proj = (ray0 * hyp0) + (ray1 * hyp1) + (ray2 * hyp2);
1299
1300 adj0 = ray0 * proj;
1301 adj1 = ray1 * proj;
1302 opp0 = hyp0 - adj0;
1303 adj2 = ray2 * proj;
1304 if(fabs(opp0) > dcutoff)
1305 return 0;
1306 opp1 = hyp1 - adj1;
1307 opp2 = hyp2 - adj2;
1308 if(fabs(opp1) > dcutoff)
1309 return 0;
1310 if(fabs(opp2) > dcutoff)
1311 return 0;
1312
1313 opp_len_sq = (opp0 * opp0) + (opp1 * opp1) + (opp2 * opp2);
1314
1315 if(opp_len_sq <= cutoff2) { /* line hits the virtual sphere */
1316
1317 float scaled_dist = proj - (float) sqrt1f(cutoff2 - opp_len_sq);
1318
1319 *(dist) = scaled_dist / d1;
1320
1321 return 1;
1322 }
1323
1324 return 0;
1325 }
1326 }
1327
1328
1329 /*========================================================================*/
BasisSetupMatrix(CBasis * I)1330 void BasisSetupMatrix(CBasis * I)
1331 {
1332 float oldZ[3] = { 0.0, 0.0, 1.0 };
1333 float newY[3];
1334 float dotgle, angle;
1335
1336 cross_product3f(oldZ, I->LightNormal, newY);
1337
1338 dotgle = dot_product3f(oldZ, I->LightNormal);
1339
1340 if((1.0 - fabs(dotgle)) < kR_SMALL4) {
1341 dotgle = (float) (dotgle / fabs(dotgle));
1342 newY[0] = 0.0;
1343 newY[1] = 1.0;
1344 newY[2] = 0.0;
1345 }
1346
1347 normalize3f(newY);
1348
1349 angle = (float) (-acos(dotgle));
1350
1351 /* now all we gotta do is effect a rotation about the new Y axis to line up new Z with Z */
1352
1353 rotation_to_matrix33f(newY, angle, I->Matrix);
1354
1355 /*
1356 printf("%8.3f %8.3f %8.3f %8.3f\n",angle*180.0/cPI,newY[0],newY[1],newY[2]);
1357
1358 matrix_transform33f3f(I->Matrix,newY,test);
1359 printf(" %8.3f %8.3f %8.3f\n",test[0],test[1],test[2]);
1360
1361 printf(" %8.3f %8.3f %8.3f\n",I->LightNormal[0],I->LightNormal[1],I->LightNormal[2]);
1362 matrix_transform33f3f(I->Matrix,I->LightNormal,test);
1363 printf(" %8.3f %8.3f %8.3f\n",test[0],test[1],test[2]);
1364
1365 printf(">%8.3f %8.3f %8.3f\n",I->Matrix[0][0],I->Matrix[0][1],I->Matrix[0][2]);
1366 printf(">%8.3f %8.3f %8.3f\n",I->Matrix[1][0],I->Matrix[1][1],I->Matrix[1][2]);
1367 printf(">%8.3f %8.3f %8.3f\n",I->Matrix[2][0],I->Matrix[2][1],I->Matrix[2][2]);
1368 */
1369 }
1370
1371
1372 /*========================================================================*/
BasisGetTriangleFlatDotgle(CBasis * I,RayInfo * r,int i)1373 void BasisGetTriangleFlatDotgle(CBasis * I, RayInfo * r, int i)
1374 {
1375 float *n0 = I->Normal + (3 * I->Vert2Normal[i]);
1376 r->flat_dotgle = n0[2];
1377 }
1378
BasisGetTriangleFlatDotglePerspective(CBasis * I,RayInfo * r,int i)1379 void BasisGetTriangleFlatDotglePerspective(CBasis * I, RayInfo * r, int i)
1380 {
1381 float *n0 = I->Normal + (3 * I->Vert2Normal[i]);
1382 r->flat_dotgle = -dot_product3f(r->dir, n0);
1383 }
1384
1385
1386 /*========================================================================*/
BasisGetEllipsoidNormal(CBasis * I,RayInfo * r,int i,int perspective)1387 void BasisGetEllipsoidNormal(CBasis * I, RayInfo * r, int i, int perspective)
1388 {
1389 if(perspective) {
1390 r->impact[0] = r->base[0] + r->dir[0] * r->dist;
1391 r->impact[1] = r->base[1] + r->dir[1] * r->dist;
1392 r->impact[2] = r->base[2] + r->dir[2] * r->dist;
1393 } else {
1394 r->impact[0] = r->base[0];
1395 r->impact[1] = r->base[1];
1396 r->impact[2] = r->base[2] - r->dist;
1397 }
1398
1399 {
1400 float *n1 = I->Normal + (3 * I->Vert2Normal[i]);
1401 float *n2 = n1 + 3, *n3 = n1 + 6;
1402 float *scale = r->prim->n0;
1403 float d1, d2, d3, s1, s2, s3;
1404 float comp1[3], comp2[3], comp3[3];
1405 float direct[3], surfnormal[3];
1406
1407 direct[0] = r->impact[0] - r->sphere[0];
1408 direct[1] = r->impact[1] - r->sphere[1];
1409 direct[2] = r->impact[2] - r->sphere[2];
1410
1411 normalize3f(direct);
1412
1413 d1 = dot_product3f(direct, n1);
1414 d2 = dot_product3f(direct, n2);
1415 d3 = dot_product3f(direct, n3);
1416
1417 if(scale[0] > R_SMALL8) {
1418 s1 = d1 / (scale[0] * scale[0]);
1419 } else {
1420 s1 = 0.0F;
1421 }
1422 if(scale[1] > R_SMALL8) {
1423 s2 = d2 / (scale[1] * scale[1]);
1424 } else {
1425 s2 = 0.0F;
1426 }
1427 if(scale[2] > R_SMALL8) {
1428 s3 = d3 / (scale[2] * scale[2]);
1429 } else {
1430 s3 = 0.0F;
1431 }
1432
1433 scale3f(n1, s1, comp1);
1434 scale3f(n2, s2, comp2);
1435 scale3f(n3, s3, comp3);
1436
1437 copy3f(comp1, surfnormal);
1438 add3f(comp2, surfnormal, surfnormal);
1439 add3f(comp3, surfnormal, surfnormal);
1440
1441 normalize23f(surfnormal, r->surfnormal);
1442 }
1443 }
1444
BasisGetTriangleNormal(CBasis * I,RayInfo * r,int i,float * fc,int perspective)1445 void BasisGetTriangleNormal(CBasis * I, RayInfo * r, int i, float *fc, int perspective)
1446 {
1447 float *n0, w2, fc0, fc1, fc2;
1448 float vt1[3];
1449 CPrimitive *lprim = r->prim;
1450
1451 if(perspective) {
1452 r->impact[0] = r->base[0] + r->dir[0] * r->dist;
1453 r->impact[1] = r->base[1] + r->dir[1] * r->dist;
1454 r->impact[2] = r->base[2] + r->dir[2] * r->dist;
1455 } else {
1456 r->impact[0] = r->base[0];
1457 r->impact[1] = r->base[1];
1458 r->impact[2] = r->base[2] - r->dist;
1459 }
1460
1461 n0 = I->Normal + (3 * I->Vert2Normal[i]) + 3; /* skip triangle normal */
1462 w2 = 1.0F - (r->tri1 + r->tri2);
1463 /* printf("%8.3f %8.3f\n",r->tri[1],r->tri[2]); */
1464
1465 fc0 = (lprim->c2[0] * r->tri1) + (lprim->c3[0] * r->tri2) + (lprim->c1[0] * w2);
1466 fc1 = (lprim->c2[1] * r->tri1) + (lprim->c3[1] * r->tri2) + (lprim->c1[1] * w2);
1467 fc2 = (lprim->c2[2] * r->tri1) + (lprim->c3[2] * r->tri2) + (lprim->c1[2] * w2);
1468
1469 r->trans = (lprim->tr[1] * r->tri1) + (lprim->tr[2] * r->tri2) + (lprim->tr[0] * w2);
1470
1471 scale3f(n0 + 3, r->tri1, r->surfnormal);
1472 scale3f(n0 + 6, r->tri2, vt1);
1473 add3f(vt1, r->surfnormal, r->surfnormal);
1474
1475 scale3f(n0, w2, vt1);
1476 add3f(vt1, r->surfnormal, r->surfnormal);
1477
1478 normalize3f(r->surfnormal);
1479
1480 fc[0] = fc0;
1481 fc[1] = fc1;
1482 fc[2] = fc2;
1483 }
1484
1485
1486 /*========================================================================*/
1487
1488 #ifdef PROFILE_BASIS
1489 int n_cells = 0;
1490 int n_prims = 0;
1491 int n_triangles = 0;
1492 int n_spheres = 0;
1493 int n_cylinders = 0;
1494 int n_sausages = 0;
1495 int n_skipped = 0;
1496 #endif
1497
BasisHitPerspective(BasisCallRec * BC)1498 int BasisHitPerspective(BasisCallRec * BC)
1499 {
1500 CBasis *BI = BC->Basis;
1501 MapType *map = BI->Map;
1502 int iMin0 = map->iMin[0];
1503 int iMin1 = map->iMin[1];
1504 int iMin2 = map->iMin[2];
1505 int iMax0 = map->iMax[0];
1506 int iMax1 = map->iMax[1];
1507 int iMax2 = map->iMax[2];
1508 int a, b, c;
1509
1510 float iDiv = map->recipDiv;
1511 float base0, base1, base2;
1512
1513 float min0 = map->Min[0] * iDiv;
1514 float min1 = map->Min[1] * iDiv;
1515 float min2 = map->Min[2] * iDiv;
1516
1517 int new_ray = !BC->pass;
1518 RayInfo *r = BC->rr;
1519
1520 MapCache *cache = &BC->cache;
1521 int *cache_cache = cache->Cache;
1522 int *cache_CacheLink = cache->CacheLink;
1523
1524 CPrimitive *r_prim = NULL;
1525
1526 if(new_ray) { /* see if we can eliminate this ray right away using the mask */
1527
1528 base0 = (r->base[0] * iDiv) - min0;
1529 base1 = (r->base[1] * iDiv) - min1;
1530
1531 a = (int) base0;
1532 b = (int) base1;
1533 a += MapBorder;
1534 b += MapBorder;
1535 if(a < iMin0)
1536 a = iMin0;
1537 else if(a > iMax0)
1538 a = iMax0;
1539 if(b < iMin1)
1540 b = iMin1;
1541 else if(b > iMax1)
1542 b = iMax1;
1543
1544 if(!*(map->EMask + a * map->Dim[1] + b))
1545 return -1;
1546 }
1547
1548 {
1549 int last_a = -1, last_b = -1, last_c = -1;
1550 int allow_break;
1551 int minIndex = -1;
1552
1553 float step0, step1, step2;
1554 float back_dist = BC->back_dist;
1555
1556 const float _0 = 0.0F, _1 = 1.0F;
1557 float r_tri1 = _0, r_tri2 = _0, r_dist, dist; /* zero inits to suppress compiler warnings */
1558 float r_sphere0 = _0, r_sphere1 = _0, r_sphere2 = _0;
1559 int h, *ip;
1560 int excl_trans_flag;
1561 int *elist, local_iflag = false;
1562 int terminal = -1;
1563 int *ehead = map->EHead;
1564 int d1d2 = map->D1D2;
1565 int d2 = map->Dim[2];
1566 const int *vert2prim = BC->vert2prim;
1567 const float excl_trans = BC->excl_trans;
1568 const float BasisFudge0 = BC->fudge0;
1569 const float BasisFudge1 = BC->fudge1;
1570 int v2p;
1571 int i, ii;
1572 int n_vert = BI->NVertex, n_eElem = map->NEElem;
1573 int except1 = BC->except1;
1574 int except2 = BC->except2;
1575 int check_interior_flag = BC->check_interior && !BC->pass;
1576 float sph[3], vt[3], tri1 = _0, tri2;
1577 CPrimitive *BC_prim = BC->prim;
1578 int *BI_Vert2Normal = BI->Vert2Normal;
1579 float *BI_Vertex = BI->Vertex;
1580 float *BI_Precomp = BI->Precomp;
1581 float *BI_Normal = BI->Normal;
1582 float *BI_Radius = BI->Radius;
1583 float *BI_Radius2 = BI->Radius2;
1584 copy3f(r->base, vt);
1585
1586 elist = map->EList;
1587
1588 r_dist = FLT_MAX;
1589
1590 excl_trans_flag = (excl_trans != _0);
1591
1592 if(except1 >= 0)
1593 except1 = vert2prim[except1];
1594 if(except2 >= 0)
1595 except2 = vert2prim[except2];
1596
1597 MapCacheReset(cache);
1598
1599 { /* take steps with a Z-size equil to the grid spacing */
1600 float div = iDiv * (-MapGetDiv(BI->Map) / r->dir[2]);
1601 step0 = r->dir[0] * div;
1602 step1 = r->dir[1] * div;
1603 step2 = r->dir[2] * div;
1604 }
1605
1606 base0 = (r->skip[0] * iDiv) - min0;
1607 base1 = (r->skip[1] * iDiv) - min1;
1608 base2 = (r->skip[2] * iDiv) - min2;
1609
1610 allow_break = false;
1611 while(1) {
1612 int inside_code;
1613 int clamped;
1614
1615 a = ((int) base0);
1616 b = ((int) base1);
1617 c = ((int) base2);
1618
1619 inside_code = 1;
1620 clamped = false;
1621
1622 a += MapBorder;
1623 b += MapBorder;
1624 c += MapBorder;
1625 #define EDGE_ALLOWANCE 1
1626
1627 if(a < iMin0) {
1628 if(((iMin0 - a) > EDGE_ALLOWANCE) && allow_break)
1629 break;
1630 else {
1631 a = iMin0;
1632 clamped = true;
1633 }
1634 } else if(a > iMax0) {
1635 if(((a - iMax0) > EDGE_ALLOWANCE) && allow_break)
1636 break;
1637 else {
1638 a = iMax0;
1639 clamped = true;
1640 }
1641 }
1642 if(b < iMin1) {
1643 if(((iMin1 - b) > EDGE_ALLOWANCE) && allow_break)
1644 break;
1645 else {
1646 b = iMin1;
1647 clamped = true;
1648 }
1649 } else if(b > iMax1) {
1650 if(((b - iMax1) > EDGE_ALLOWANCE) && allow_break)
1651 break;
1652 else {
1653 b = iMax1;
1654 clamped = true;
1655 }
1656 }
1657 if(c < iMin2) {
1658 if((iMin2 - c) > EDGE_ALLOWANCE)
1659 break;
1660 else {
1661 c = iMin2;
1662 clamped = true;
1663 }
1664 } else if(c > iMax2) {
1665 if((c - iMax2) > EDGE_ALLOWANCE)
1666 inside_code = 0;
1667 else {
1668 c = iMax2;
1669 clamped = true;
1670 }
1671 }
1672 if(inside_code && (((a != last_a) || (b != last_b) || (c != last_c)))) {
1673 int new_min_index;
1674 h = *(ehead + (d1d2 * a) + (d2 * b) + c);
1675
1676 new_min_index = -1;
1677
1678 if(!clamped) /* don't discard a ray until it has hit the objective at least once */
1679 allow_break = true;
1680
1681 if((terminal > 0) && (last_c != c)) {
1682 if(!terminal--)
1683 break;
1684 }
1685 if((h > 0) && (h < n_eElem)) {
1686 int do_loop;
1687
1688 ip = elist + h;
1689 last_a = a;
1690 i = *(ip++);
1691 last_b = b;
1692 do_loop = ((i >= 0) && (i < n_vert));
1693 last_c = c;
1694
1695 while(do_loop) { /* n_vert checking is a bug workaround */
1696 CPrimitive *prm;
1697 v2p = vert2prim[i];
1698 ii = *(ip++);
1699 prm = BC_prim + v2p;
1700 do_loop = ((ii >= 0) && (ii < n_vert));
1701 /* if((v2p != except1) && (v2p != except2) && (!MapCached(cache, v2p))) { */
1702 if((v2p != except1) && (v2p != except2) && (!cache_cache[v2p])) {
1703 int prm_type = prm->type;
1704
1705 /*MapCache(cache,v2p); */
1706 cache_cache[v2p] = 1;
1707 cache_CacheLink[v2p] = cache->CacheStart;
1708 cache->CacheStart = v2p;
1709
1710 switch (prm_type) {
1711 case cPrimTriangle:
1712 case cPrimCharacter:
1713 {
1714 float *dir = r->dir;
1715 float *d10 = BI_Precomp + BI_Vert2Normal[i] * 3;
1716 float *d20 = d10 + 3;
1717 float *v0;
1718 float det, inv_det;
1719 float pvec0, pvec1, pvec2;
1720 float dir0 = dir[0], dir1 = dir[1], dir2 = dir[2];
1721 float d20_0 = d20[0], d20_1 = d20[1], d20_2 = d20[2];
1722 float d10_0 = d10[0], d10_1 = d10[1], d10_2 = d10[2];
1723
1724 /* cross_product3f(dir, d20, pvec); */
1725
1726 pvec0 = dir1 * d20_2 - dir2 * d20_1;
1727 pvec1 = dir2 * d20_0 - dir0 * d20_2;
1728 pvec2 = dir0 * d20_1 - dir1 * d20_0;
1729
1730 /* det = dot_product3f(pvec, d10); */
1731
1732 det = pvec0 * d10_0 + pvec1 * d10_1 + pvec2 * d10_2;
1733
1734 v0 = BI_Vertex + prm->vert * 3;
1735 if((det >= EPSILON) || (det <= -EPSILON)) {
1736 float tvec0, tvec1, tvec2;
1737 float qvec0, qvec1, qvec2;
1738
1739 inv_det = _1 / det;
1740
1741 /* subtract3f(vt,v0,tvec); */
1742
1743 tvec0 = vt[0] - v0[0];
1744 tvec1 = vt[1] - v0[1];
1745 tvec2 = vt[2] - v0[2];
1746
1747 /* dot_product3f(tvec,pvec) * inv_det; */
1748 tri1 = (tvec0 * pvec0 + tvec1 * pvec1 + tvec2 * pvec2) * inv_det;
1749
1750 /* cross_product3f(tvec,d10,qvec); */
1751
1752 qvec0 = tvec1 * d10_2 - tvec2 * d10_1;
1753 qvec1 = tvec2 * d10_0 - tvec0 * d10_2;
1754
1755 if((tri1 >= BasisFudge0) && (tri1 <= BasisFudge1)) {
1756 qvec2 = tvec0 * d10_1 - tvec1 * d10_0;
1757
1758 /* dot_product3f(dir, qvec) * inv_det; */
1759 tri2 = (dir0 * qvec0 + dir1 * qvec1 + dir2 * qvec2) * inv_det;
1760
1761 /* dot_product3f(d20, qvec) * inv_det; */
1762 dist = (d20_0 * qvec0 + d20_1 * qvec1 + d20_2 * qvec2) * inv_det;
1763
1764 if((tri2 >= BasisFudge0) && (tri2 <= BasisFudge1)
1765 && ((tri1 + tri2) <= BasisFudge1)) {
1766 if((dist < r_dist) && (dist >= _0) && (dist <= back_dist)
1767 && (prm->trans != _1)) {
1768 new_min_index = prm->vert;
1769 r_tri1 = tri1;
1770 r_tri2 = tri2;
1771 r_dist = dist;
1772 }
1773 }
1774 }
1775 }
1776 }
1777 break;
1778 case cPrimSphere:
1779 {
1780 if(LineClipPoint(r->base, r->dir,
1781 BI_Vertex + i * 3, &dist,
1782 BI_Radius[i], BI_Radius2[i])) {
1783 if((dist < r_dist) && (prm->trans != _1)) {
1784 if((dist >= _0) && (dist <= back_dist)) {
1785 new_min_index = prm->vert;
1786 r_dist = dist;
1787 } else if(check_interior_flag && (dist <= back_dist)) {
1788 if(diffsq3f(vt, BI_Vertex + i * 3) < BI_Radius2[i]) {
1789
1790 local_iflag = true;
1791 r_prim = prm;
1792 r_dist = _0;
1793 new_min_index = prm->vert;
1794 }
1795 }
1796 }
1797 }
1798 }
1799 break;
1800 case cPrimEllipsoid:
1801 {
1802 if(LineClipPoint(r->base, r->dir,
1803 BI_Vertex + i * 3, &dist,
1804 BI_Radius[i], BI_Radius2[i])) {
1805 if((dist < r_dist) && (prm->trans != _1)) {
1806 float *n1 = BI_Normal + BI_Vert2Normal[i] * 3;
1807 if(LineClipEllipsoidPoint(r->base, r->dir,
1808 BI_Vertex + i * 3, &dist,
1809 BI_Radius[i], BI_Radius2[i],
1810 prm->n0, n1, n1 + 3, n1 + 6)) {
1811 if(dist < r_dist) {
1812 if((dist >= _0) && (dist <= back_dist)) {
1813 new_min_index = prm->vert;
1814 r_dist = dist;
1815 }
1816 }
1817 }
1818 }
1819 }
1820 }
1821 break;
1822
1823 case cPrimCylinder:
1824 if(LineToSphereCapped(r->base, r->dir, BI_Vertex + i * 3,
1825 BI_Normal + BI_Vert2Normal[i] * 3,
1826 BI_Radius[i], prm->l1, sph, &tri1,
1827 prm->cap1, prm->cap2)) {
1828 if(LineClipPoint
1829 (r->base, r->dir, sph, &dist, BI_Radius[i], BI_Radius2[i])) {
1830 if((dist < r_dist) && (prm->trans != _1)) {
1831 if((dist >= _0) && (dist <= back_dist)) {
1832 if(prm->l1 > kR_SMALL4)
1833 r_tri1 = tri1 / prm->l1;
1834
1835 r_sphere0 = sph[0];
1836 r_sphere1 = sph[1];
1837 r_sphere2 = sph[2];
1838 new_min_index = prm->vert;
1839 r_dist = dist;
1840 } else if(check_interior_flag && (dist <= back_dist)) {
1841 if(FrontToInteriorSphereCapped(vt,
1842 BI_Vertex + i * 3,
1843 BI_Normal + BI_Vert2Normal[i] * 3,
1844 BI_Radius[i],
1845 BI_Radius2[i],
1846 prm->l1, prm->cap1, prm->cap2)) {
1847 local_iflag = true;
1848 r_prim = prm;
1849 r_dist = _0;
1850
1851 new_min_index = prm->vert;
1852 }
1853 }
1854 }
1855 }
1856 }
1857 break;
1858 case cPrimCone:
1859 {
1860 float sph_rad, sph_rad_sq;
1861 if(ConeLineToSphereCapped(r->base, r->dir, BI_Vertex + i * 3,
1862 BI_Normal + BI_Vert2Normal[i] * 3,
1863 BI_Radius[i], prm->r2, prm->l1, sph, &tri1,
1864 &sph_rad, &sph_rad_sq,
1865 prm->cap1, prm->cap2)) {
1866
1867 if(LineClipPoint(r->base, r->dir, sph, &dist, sph_rad, sph_rad_sq)) {
1868 if((dist < r_dist) && (prm->trans != _1)) {
1869 if((dist >= _0) && (dist <= back_dist)) {
1870 if(prm->l1 > kR_SMALL4)
1871 r_tri1 = tri1 / prm->l1; /* color blending */
1872 r_sphere0 = sph[0];
1873 r_sphere1 = sph[1];
1874 r_sphere2 = sph[2];
1875 new_min_index = prm->vert;
1876 r_dist = dist;
1877 } else if(check_interior_flag && (dist <= back_dist)) {
1878 if(FrontToInteriorSphereCapped(vt,
1879 BI_Vertex + i * 3,
1880 BI_Normal +
1881 BI_Vert2Normal[i] * 3,
1882 BI_Radius[i], BI_Radius2[i],
1883 prm->l1, prm->cap1, prm->cap2)) {
1884 local_iflag = true;
1885 r_prim = prm;
1886 r_dist = _0;
1887 new_min_index = prm->vert;
1888 }
1889 }
1890 }
1891 }
1892 }
1893 }
1894 break;
1895 case cPrimSausage:
1896 if(LineToSphere(r->base, r->dir,
1897 BI_Vertex + i * 3, BI_Normal + BI_Vert2Normal[i] * 3,
1898 BI_Radius[i], prm->l1, sph, &tri1)) {
1899
1900 if(LineClipPoint
1901 (r->base, r->dir, sph, &dist, BI_Radius[i], BI_Radius2[i])) {
1902
1903 int tmp_flag = false;
1904 if((dist < r_dist) && (prm->trans != _1)) {
1905 if((dist >= _0) && (dist <= back_dist)) {
1906 tmp_flag = true;
1907 if(excl_trans_flag) {
1908 if((prm->trans > _0) && (dist < excl_trans))
1909 tmp_flag = false;
1910 }
1911 if(tmp_flag) {
1912
1913 if(prm->l1 > kR_SMALL4)
1914 r_tri1 = tri1 / prm->l1;
1915
1916 r_sphere0 = sph[0];
1917 r_sphere1 = sph[1];
1918 r_sphere2 = sph[2];
1919 new_min_index = prm->vert;
1920 r_dist = dist;
1921
1922 }
1923 } else if(check_interior_flag && (dist <= back_dist)) {
1924 if(FrontToInteriorSphere(vt, BI_Vertex + i * 3,
1925 BI_Normal + BI_Vert2Normal[i] * 3,
1926 BI_Radius[i], BI_Radius2[i], prm->l1)) {
1927 local_iflag = true;
1928 r_prim = prm;
1929 r_dist = _0;
1930 new_min_index = prm->vert;
1931 }
1932 }
1933 }
1934 }
1935 }
1936 break;
1937 } /* end of switch */
1938 }
1939 /* end of if */
1940 i = ii;
1941
1942 } /* end of while */
1943
1944 if(local_iflag) {
1945 r->prim = r_prim;
1946 r->dist = r_dist;
1947
1948 break;
1949 }
1950
1951 if(new_min_index > -1) {
1952
1953 minIndex = new_min_index;
1954
1955 r_prim = BC_prim + vert2prim[minIndex];
1956
1957 if((r_prim->type == cPrimSphere) || (r_prim->type == cPrimEllipsoid)) {
1958 const float *vv = BI->Vertex + minIndex * 3;
1959 r_sphere0 = vv[0];
1960 r_sphere1 = vv[1];
1961 r_sphere2 = vv[2];
1962 }
1963
1964 BC->interior_flag = local_iflag;
1965 r->tri1 = r_tri1;
1966 r->tri2 = r_tri2;
1967 r->prim = r_prim;
1968 r->dist = r_dist;
1969 r->sphere[0] = r_sphere0;
1970 r->sphere[1] = r_sphere1;
1971 r->sphere[2] = r_sphere2;
1972 }
1973 } /* if -- h valid */
1974 }
1975 /* end of if */
1976 if(minIndex > -1) {
1977 if(terminal < 0)
1978 terminal = EDGE_ALLOWANCE + 1;
1979 }
1980
1981 base0 += step0;
1982 base1 += step1;
1983 base2 += step2;
1984 /* advance through the map one block at a time -- note that this is a crappy way to walk through the map... */
1985 }
1986
1987 BC->interior_flag = local_iflag;
1988 return (minIndex);
1989 }
1990 }
1991
BasisHitOrthoscopic(BasisCallRec * BC)1992 int BasisHitOrthoscopic(BasisCallRec * BC)
1993 {
1994 const float _0 = 0.0F, _1 = 1.0F;
1995 float oppSq, dist = _0, sph[3], vt[3], tri1, tri2;
1996 int a, b, c, h, *ip;
1997 int excl_trans_flag;
1998 int check_interior_flag;
1999 int *elist, local_iflag = false;
2000 float minusZ[3] = { 0.0F, 0.0F, -1.0F };
2001
2002 CBasis *BI = BC->Basis;
2003 RayInfo *r = BC->rr;
2004
2005 if(MapInsideXY(BI->Map, r->base, &a, &b, &c)) {
2006 int minIndex = -1;
2007 int v2p;
2008 int i, ii;
2009 int *xxtmp;
2010 int do_loop;
2011 int except1 = BC->except1;
2012 int except2 = BC->except2;
2013 int n_vert = BI->NVertex, n_eElem = BI->Map->NEElem;
2014 const int *vert2prim = BC->vert2prim;
2015 const float front = BC->front;
2016 const float back = BC->back;
2017 const float excl_trans = BC->excl_trans;
2018 const float BasisFudge0 = BC->fudge0;
2019 const float BasisFudge1 = BC->fudge1;
2020
2021 MapCache *cache = &BC->cache;
2022
2023 float r_tri1 = _0, r_tri2 = _0, r_dist = _0; /* zero inits to suppress compiler warnings */
2024 float r_sphere0 = _0, r_sphere1 = _0, r_sphere2 = _0;
2025 CPrimitive *r_prim = NULL;
2026
2027 check_interior_flag = BC->check_interior && (!BC->pass);
2028
2029 /* assumption: always heading in the negative Z direction with our vector... */
2030 vt[0] = r->base[0];
2031 vt[1] = r->base[1];
2032 vt[2] = r->base[2] - front;
2033
2034 if(except1 >= 0)
2035 except1 = vert2prim[except1];
2036 if(except2 >= 0)
2037 except2 = vert2prim[except2];
2038
2039 excl_trans_flag = (excl_trans != _0);
2040
2041 r_dist = FLT_MAX;
2042
2043 xxtmp = BI->Map->EHead + (a * BI->Map->D1D2) + (b * BI->Map->Dim[2]) + c;
2044
2045 MapCacheReset(cache);
2046
2047 elist = BI->Map->EList;
2048
2049 while(c >= MapBorder) {
2050 h = *xxtmp;
2051 if((h > 0) && (h < n_eElem)) {
2052 ip = elist + h;
2053 i = *(ip++);
2054 do_loop = ((i >= 0) && (i < n_vert));
2055 while(do_loop) {
2056 ii = *(ip++);
2057 v2p = vert2prim[i];
2058 do_loop = ((ii >= 0) && (ii < n_vert));
2059
2060 if((v2p != except1) && (v2p != except2) && (!MapCached(cache, v2p))) {
2061 CPrimitive *prm = BC->prim + v2p;
2062 MapCache(cache, v2p);
2063
2064 switch (prm->type) {
2065 case cPrimTriangle:
2066 case cPrimCharacter:
2067 if(!prm->cull) {
2068 float *pre = BI->Precomp + BI->Vert2Normal[i] * 3;
2069
2070 if(pre[6]) {
2071 float *vert0 = BI->Vertex + prm->vert * 3;
2072
2073 float tvec0 = vt[0] - vert0[0];
2074 float tvec1 = vt[1] - vert0[1];
2075
2076 tri1 = (tvec0 * pre[4] - tvec1 * pre[3]) * pre[7];
2077 tri2 = -(tvec0 * pre[1] - tvec1 * pre[0]) * pre[7];
2078
2079 if(!((tri1 < BasisFudge0) || (tri2 < BasisFudge0) ||
2080 (tri1 > BasisFudge1) || ((tri1 + tri2) > BasisFudge1))) {
2081 dist = (r->base[2] - (tri1 * pre[2]) - (tri2 * pre[5]) - vert0[2]);
2082
2083 if((dist < r_dist) && (dist >= front) &&
2084 (dist <= back) && (prm->trans != _1)) {
2085 minIndex = prm->vert;
2086 r_tri1 = tri1;
2087 r_tri2 = tri2;
2088 r_dist = dist;
2089 }
2090 }
2091 }
2092 }
2093 break;
2094
2095 case cPrimSphere:
2096 oppSq = ZLineClipPoint(r->base, BI->Vertex + i * 3, &dist, BI->Radius[i]);
2097 if(oppSq <= BI->Radius2[i]) {
2098 dist = (float) (sqrt1f(dist) - sqrt1f((BI->Radius2[i] - oppSq)));
2099
2100 if((dist < r_dist) && (prm->trans != _1)) {
2101 if((dist >= front) && (dist <= back)) {
2102 minIndex = prm->vert;
2103 r_dist = dist;
2104 } else if(check_interior_flag) {
2105 if(diffsq3f(vt, BI->Vertex + i * 3) < BI->Radius2[i]) {
2106 local_iflag = true;
2107 r_prim = prm;
2108 r_dist = front;
2109 minIndex = prm->vert;
2110 }
2111 }
2112 }
2113 }
2114 break;
2115 case cPrimEllipsoid:
2116 oppSq = ZLineClipPoint(r->base, BI->Vertex + i * 3, &dist, BI->Radius[i]);
2117 if(oppSq <= BI->Radius2[i]) {
2118
2119 dist = (float) (sqrt1f(dist) - sqrt1f((BI->Radius2[i] - oppSq)));
2120
2121 if((dist < r_dist) && (prm->trans != _1)) {
2122 float *n1 = BI->Normal + BI->Vert2Normal[i] * 3;
2123 if(LineClipEllipsoidPoint(r->base, minusZ,
2124 BI->Vertex + i * 3, &dist,
2125 BI->Radius[i], BI->Radius2[i],
2126 prm->n0, n1, n1 + 3, n1 + 6)) {
2127 if(dist < r_dist) {
2128 if((dist >= _0) && (dist <= back)) {
2129 minIndex = prm->vert;
2130 r_dist = dist;
2131 }
2132 }
2133 }
2134 }
2135 }
2136 break;
2137
2138 case cPrimCylinder:
2139 if(ZLineToSphereCapped(r->base, BI->Vertex + i * 3,
2140 BI->Normal + BI->Vert2Normal[i] * 3,
2141 BI->Radius[i], prm->l1, sph, &tri1, prm->cap1,
2142 prm->cap2, BI->Precomp + BI->Vert2Normal[i] * 3)) {
2143 oppSq = ZLineClipPoint(r->base, sph, &dist, BI->Radius[i]);
2144 if(oppSq <= BI->Radius2[i]) {
2145 dist = (float) (sqrt1f(dist) - sqrt1f((BI->Radius2[i] - oppSq)));
2146
2147 if((dist < r_dist) && (prm->trans != _1)) {
2148 if((dist >= front) && (dist <= back)) {
2149 if(prm->l1 > kR_SMALL4)
2150 r_tri1 = tri1 / prm->l1;
2151
2152 r_sphere0 = sph[0];
2153 r_sphere1 = sph[1];
2154 r_sphere2 = sph[2];
2155 minIndex = prm->vert;
2156 r_dist = dist;
2157 } else if(check_interior_flag) {
2158 if(FrontToInteriorSphereCapped(vt,
2159 BI->Vertex + i * 3,
2160 BI->Normal + BI->Vert2Normal[i] * 3,
2161 BI->Radius[i],
2162 BI->Radius2[i],
2163 prm->l1, prm->cap1, prm->cap2)) {
2164 local_iflag = true;
2165 r_prim = prm;
2166 r_dist = front;
2167 minIndex = prm->vert;
2168 }
2169 }
2170 }
2171 }
2172 }
2173 break;
2174 case cPrimCone:
2175 {
2176 float sph_rad, sph_rad_sq;
2177 if(ConeLineToSphereCapped(r->base, minusZ, BI->Vertex + i * 3,
2178 BI->Normal + BI->Vert2Normal[i] * 3,
2179 BI->Radius[i], prm->r2, prm->l1, sph, &tri1,
2180 &sph_rad, &sph_rad_sq, prm->cap1, prm->cap2)) {
2181
2182 oppSq = ZLineClipPoint(r->base, sph, &dist, sph_rad);
2183 if(oppSq <= sph_rad_sq) {
2184 dist = (float) (sqrt1f(dist) - sqrt1f((sph_rad_sq - oppSq)));
2185
2186 if((dist < r_dist) && (prm->trans != _1)) {
2187 if((dist >= front) && (dist <= back)) {
2188 if(prm->l1 > kR_SMALL4)
2189 r_tri1 = tri1 / prm->l1;
2190
2191 r_sphere0 = sph[0];
2192 r_sphere1 = sph[1];
2193 r_sphere2 = sph[2];
2194 minIndex = prm->vert;
2195 r_dist = dist;
2196 } else if(check_interior_flag) {
2197 if(FrontToInteriorSphereCapped(vt,
2198 BI->Vertex + i * 3,
2199 BI->Normal +
2200 BI->Vert2Normal[i] * 3, sph_rad,
2201 sph_rad_sq, prm->l1, prm->cap1,
2202 prm->cap2)) {
2203 local_iflag = true;
2204 r_prim = prm;
2205 r_dist = front;
2206 minIndex = prm->vert;
2207 }
2208 }
2209 }
2210 }
2211 }
2212 }
2213 break;
2214 case cPrimSausage:
2215 if(ZLineToSphere
2216 (r->base, BI->Vertex + i * 3, BI->Normal + BI->Vert2Normal[i] * 3,
2217 BI->Radius[i], prm->l1, sph, &tri1,
2218 BI->Precomp + BI->Vert2Normal[i] * 3)) {
2219 oppSq = ZLineClipPoint(r->base, sph, &dist, BI->Radius[i]);
2220 if(oppSq <= BI->Radius2[i]) {
2221 int tmp_flag = false;
2222
2223 dist = (float) (sqrt1f(dist) - sqrt1f((BI->Radius2[i] - oppSq)));
2224 if((dist < r_dist) && (prm->trans != _1)) {
2225 if((dist >= front) && (dist <= back)) {
2226 tmp_flag = true;
2227 if(excl_trans_flag) {
2228 if((prm->trans > _0) && (dist < excl_trans))
2229 tmp_flag = false;
2230 }
2231 if(tmp_flag) {
2232 if(prm->l1 > kR_SMALL4)
2233 r_tri1 = tri1 / prm->l1;
2234
2235 r_sphere0 = sph[0];
2236 r_sphere1 = sph[1];
2237 r_sphere2 = sph[2];
2238 minIndex = prm->vert;
2239 r_dist = dist;
2240 }
2241 } else if(check_interior_flag) {
2242 if(FrontToInteriorSphere
2243 (vt, BI->Vertex + i * 3, BI->Normal + BI->Vert2Normal[i] * 3,
2244 BI->Radius[i], BI->Radius2[i], prm->l1)) {
2245 local_iflag = true;
2246 r_prim = prm;
2247 r_dist = front;
2248 minIndex = prm->vert;
2249 }
2250 }
2251 }
2252 }
2253 }
2254 break;
2255 } /* end of switch */
2256 }
2257 /* end of if */
2258 i = ii;
2259
2260 } /* end of while */
2261 }
2262
2263 /* and of course stop when we hit the edge of the map */
2264
2265 if(local_iflag)
2266 break;
2267
2268 /* we've processed all primitives associated with this voxel,
2269 so if an intersection has been found which occurs in front of
2270 the next voxel, then we can stop */
2271
2272 if(minIndex > -1) {
2273 int aa, bb, cc;
2274
2275 vt[2] = r->base[2] - r_dist;
2276 MapLocus(BI->Map, vt, &aa, &bb, &cc);
2277 if(cc > c)
2278 break;
2279 else
2280 vt[2] = r->base[2] - front;
2281 }
2282
2283 c--;
2284 xxtmp--;
2285
2286 } /* end of while */
2287
2288 if(minIndex > -1) {
2289 r_prim = BC->prim + vert2prim[minIndex];
2290
2291 if((r_prim->type == cPrimSphere) || (r_prim->type == cPrimEllipsoid)) {
2292 const float *vv = BI->Vertex + minIndex * 3;
2293 r_sphere0 = vv[0];
2294 r_sphere1 = vv[1];
2295 r_sphere2 = vv[2];
2296 }
2297 }
2298
2299 BC->interior_flag = local_iflag;
2300 r->tri1 = r_tri1;
2301 r->tri2 = r_tri2;
2302 r->prim = r_prim;
2303 r->dist = r_dist;
2304 r->sphere[0] = r_sphere0;
2305 r->sphere[1] = r_sphere1;
2306 r->sphere[2] = r_sphere2;
2307 return (minIndex);
2308 } /* end of if */
2309 BC->interior_flag = local_iflag;
2310 return (-1);
2311 }
2312
BasisHitShadow(BasisCallRec * BC)2313 int BasisHitShadow(BasisCallRec * BC)
2314 {
2315 const float _0 = 0.0F;
2316 const float _1 = 1.0F;
2317 float oppSq, dist = _0, tri1, tri2;
2318 float sph[3], vt[3];
2319 int h, *ip;
2320 int a, b, c;
2321 int *elist, local_iflag = false;
2322 float minusZ[3] = { 0.0F, 0.0F, -1.0F };
2323 /* local copies (eliminate these extra copies later on) */
2324
2325 CBasis *BI = BC->Basis;
2326 RayInfo *r = BC->rr;
2327
2328 if(MapInsideXY(BI->Map, r->base, &a, &b, &c)) {
2329 int minIndex = -1;
2330 int v2p;
2331 int i, ii;
2332 int *xxtmp;
2333
2334 int n_vert = BI->NVertex, n_eElem = BI->Map->NEElem;
2335 int except1 = BC->except1;
2336 int except2 = BC->except2;
2337 const int *vert2prim = BC->vert2prim;
2338 const int trans_shadows = BC->trans_shadows;
2339 const int nearest_shadow = BC->nearest_shadow;
2340 const float BasisFudge0 = BC->fudge0;
2341 const float BasisFudge1 = BC->fudge1;
2342 const int label_shadow_mode = BC->label_shadow_mode;
2343 MapCache *cache = &BC->cache;
2344 int *cache_cache = cache->Cache;
2345 int *cache_CacheLink = cache->CacheLink;
2346 CPrimitive *BC_prim = BC->prim;
2347
2348 float r_tri1 = _0, r_tri2 = _0, r_dist; /* zero inits to suppress compiler warnings */
2349 float r_sphere0 = _0, r_sphere1 = _0, r_sphere2 = _0;
2350 float r_trans = _0;
2351 CPrimitive *r_prim = NULL;
2352
2353 /* assumption: always heading in the negative Z direction with our vector... */
2354 vt[0] = r->base[0];
2355 vt[1] = r->base[1];
2356
2357 if(except1 >= 0)
2358 except1 = vert2prim[except1];
2359 if(except2 >= 0)
2360 except2 = vert2prim[except2];
2361
2362 r_trans = _1;
2363 r_dist = FLT_MAX;
2364
2365 xxtmp = BI->Map->EHead + (a * BI->Map->D1D2) + (b * BI->Map->Dim[2]) + c;
2366
2367 MapCacheReset(cache);
2368
2369 elist = BI->Map->EList;
2370
2371 while(c >= MapBorder) {
2372 h = *xxtmp;
2373 if((h > 0) && (h < n_eElem)) {
2374 int do_loop;
2375 ip = elist + h;
2376 i = *(ip++);
2377 do_loop = ((i >= 0) && (i < n_vert));
2378 while(do_loop) {
2379 ii = *(ip++);
2380 v2p = vert2prim[i];
2381 do_loop = ((ii >= 0) && (ii < n_vert));
2382 if((v2p != except1) && (v2p != except2) && !MapCached(cache, v2p)) {
2383 CPrimitive *prm = BC_prim + v2p;
2384 int prm_type;
2385
2386 /*MapCache(cache,v2p); */
2387 cache_cache[v2p] = 1;
2388 prm_type = prm->type;
2389 cache_CacheLink[v2p] = cache->CacheStart;
2390 cache->CacheStart = v2p;
2391
2392 switch (prm_type) {
2393 case cPrimCharacter: /* will need special handling for character shadows */
2394 if(label_shadow_mode & 0x2) { /* if labels case shadows... */
2395 float *pre = BI->Precomp + BI->Vert2Normal[i] * 3;
2396
2397 if(pre[6]) {
2398 float *vert0 = BI->Vertex + prm->vert * 3;
2399
2400 float tvec0 = vt[0] - vert0[0];
2401 float tvec1 = vt[1] - vert0[1];
2402
2403 tri1 = (tvec0 * pre[4] - tvec1 * pre[3]) * pre[7];
2404 tri2 = -(tvec0 * pre[1] - tvec1 * pre[0]) * pre[7];
2405
2406 if(!((tri1 < BasisFudge0) ||
2407 (tri2 < BasisFudge0) ||
2408 (tri1 > BasisFudge1) || ((tri1 + tri2) > BasisFudge1))) {
2409 dist = (r->base[2] - (tri1 * pre[2]) - (tri2 * pre[5]) - vert0[2]);
2410
2411 {
2412 float fc[3];
2413 float trans;
2414
2415 r->tri1 = tri1;
2416 r->tri2 = tri2;
2417 r->dist = dist;
2418 r->prim = prm;
2419
2420 {
2421 float w2;
2422 w2 = _1 - (r->tri1 + r->tri2);
2423
2424 fc[0] =
2425 (prm->c2[0] * r->tri1) + (prm->c3[0] * r->tri2) +
2426 (prm->c1[0] * w2);
2427 fc[1] =
2428 (prm->c2[1] * r->tri1) + (prm->c3[1] * r->tri2) +
2429 (prm->c1[1] * w2);
2430 fc[2] =
2431 (prm->c2[2] * r->tri1) + (prm->c3[2] * r->tri2) +
2432 (prm->c1[2] * w2);
2433 }
2434
2435 trans = CharacterInterpolate(BI->G, prm->char_id, fc);
2436
2437 if(trans == _0) { /* opaque? return immed. */
2438 if(dist > -kR_SMALL4) {
2439 if(nearest_shadow) {
2440 if(dist < r_dist) {
2441 minIndex = prm->vert;
2442 r_tri1 = tri1;
2443 r_tri2 = tri2;
2444 r_dist = dist;
2445 r_trans = (r->trans = trans);
2446 }
2447 } else {
2448 r->prim = prm;
2449 r->trans = _0;
2450 r->dist = dist;
2451 return (1);
2452 }
2453 }
2454 } else if(trans_shadows) {
2455 if((dist > -kR_SMALL4) &&
2456 ((r_trans > trans) ||
2457 (nearest_shadow && (dist < r_dist) && (r_trans >= trans)))) {
2458 minIndex = prm->vert;
2459 r_tri1 = tri1;
2460 r_tri2 = tri2;
2461 r_dist = dist;
2462 r_trans = (r->trans = trans);
2463 }
2464 }
2465 }
2466 }
2467 }
2468 }
2469 break;
2470
2471 case cPrimTriangle:
2472 {
2473 float *pre = BI->Precomp + BI->Vert2Normal[i] * 3;
2474
2475 if(pre[6]) {
2476 float *vert0 = BI->Vertex + prm->vert * 3;
2477
2478 float tvec0 = vt[0] - vert0[0];
2479 float tvec1 = vt[1] - vert0[1];
2480
2481 tri1 = (tvec0 * pre[4] - tvec1 * pre[3]) * pre[7];
2482 tri2 = -(tvec0 * pre[1] - tvec1 * pre[0]) * pre[7];
2483 if(!((tri1 < BasisFudge0) ||
2484 (tri2 < BasisFudge0) ||
2485 (tri1 > BasisFudge1) || ((tri1 + tri2) > BasisFudge1))) {
2486 float *tr = prm->tr;
2487 float trans = _0;
2488
2489 dist = (r->base[2] - (tri1 * pre[2]) - (tri2 * pre[5]) - vert0[2]);
2490
2491 if(prm->trans != _0) {
2492 trans =
2493 (tr[1] * tri1) + (tr[2] * tri2) + (tr[0] * (_1 - (tri1 + tri2)));
2494 }
2495
2496 if(trans == _0) {
2497 if(dist > -kR_SMALL4) {
2498 if(nearest_shadow) { /* do we need the nearest shadow? */
2499 if(dist < r_dist) {
2500 minIndex = prm->vert;
2501 r_tri1 = tri1;
2502 r_tri2 = tri2;
2503 r_dist = dist;
2504 r_trans = (r->trans = trans);
2505 }
2506 } else {
2507 r->prim = prm;
2508 r->trans = _0;
2509 r->dist = dist;
2510 return (1);
2511 }
2512 }
2513 } else if(trans_shadows) {
2514 if((dist > -kR_SMALL4) &&
2515 ((r_trans > trans) ||
2516 (nearest_shadow && (dist < r_dist) && (r_trans >= trans)))) {
2517 minIndex = prm->vert;
2518 r_tri1 = tri1;
2519 r_tri2 = tri2;
2520 r_dist = dist;
2521 r_trans = (r->trans = trans);
2522 }
2523 }
2524 }
2525 }
2526 }
2527 break;
2528
2529 case cPrimSphere:
2530
2531 oppSq = ZLineClipPoint(r->base, BI->Vertex + i * 3, &dist, BI->Radius[i]);
2532 if(oppSq <= BI->Radius2[i]) {
2533 dist = (float) (sqrt1f(dist) - sqrt1f((BI->Radius2[i] - oppSq)));
2534
2535 if(prm->trans == _0) {
2536 if(dist > -kR_SMALL4) {
2537 if(nearest_shadow) {
2538 if(dist < r_dist) {
2539 minIndex = prm->vert;
2540 r_dist = dist;
2541 r_trans = (r->trans = prm->trans);
2542 }
2543 } else {
2544 r->prim = prm;
2545 r->trans = prm->trans;
2546 r->dist = dist;
2547 return (1);
2548 }
2549 }
2550 } else if(trans_shadows) {
2551 if((dist > -kR_SMALL4) &&
2552 ((r_trans > prm->trans) ||
2553 (nearest_shadow && (dist < r_dist) && (r_trans >= prm->trans)))) {
2554 minIndex = prm->vert;
2555 r_dist = dist;
2556 r_trans = (r->trans = prm->trans);
2557 }
2558 }
2559 }
2560 break;
2561
2562 case cPrimEllipsoid:
2563
2564 oppSq =
2565 ZLineClipPointNoZCheck(r->base, BI->Vertex + i * 3, &dist, BI->Radius[i]);
2566 if(oppSq <= BI->Radius2[i]) {
2567 dist = (float) (sqrt1f(dist) - sqrt1f((BI->Radius2[i] - oppSq)));
2568
2569 if((dist < r_dist) || (trans_shadows && (r_trans != _0))) {
2570 float *n1 = BI->Normal + BI->Vert2Normal[i] * 3;
2571 if(LineClipEllipsoidPoint(r->base, minusZ,
2572 BI->Vertex + i * 3, &dist,
2573 BI->Radius[i], BI->Radius2[i],
2574 prm->n0, n1, n1 + 3, n1 + 6)) {
2575
2576 if(prm->trans == _0) {
2577 if(dist > -kR_SMALL4) {
2578 if(nearest_shadow) {
2579 if(dist < r_dist) {
2580 minIndex = prm->vert;
2581 r_dist = dist;
2582 r_trans = (r->trans = prm->trans);
2583 }
2584 } else {
2585 r->prim = prm;
2586 r->trans = prm->trans;
2587 r->dist = dist;
2588 return (1);
2589 }
2590 }
2591 } else if(trans_shadows) {
2592 if((dist > -kR_SMALL4) &&
2593 ((r_trans > prm->trans) ||
2594 (nearest_shadow && (dist < r_dist)
2595 && (r_trans >= prm->trans)))) {
2596 minIndex = prm->vert;
2597 r_dist = dist;
2598 r_trans = (r->trans = prm->trans);
2599 }
2600 }
2601 }
2602 }
2603 }
2604 break;
2605 case cPrimCone:
2606 {
2607 float sph_rad, sph_rad_sq;
2608 if(ConeLineToSphereCapped(r->base, minusZ, BI->Vertex + i * 3,
2609 BI->Normal + BI->Vert2Normal[i] * 3,
2610 BI->Radius[i], prm->r2, prm->l1, sph, &tri1,
2611 &sph_rad, &sph_rad_sq, 1, 1)) {
2612
2613 oppSq = ZLineClipPoint(r->base, sph, &dist, sph_rad);
2614 if(oppSq <= sph_rad_sq) {
2615 dist = (float) (sqrt1f(dist) - sqrt1f((sph_rad_sq - oppSq)));
2616
2617 if(prm->trans == _0) {
2618 if(dist > -kR_SMALL4) {
2619 if(nearest_shadow) {
2620 if(dist < r_dist) {
2621 if(prm->l1 > kR_SMALL4)
2622 r_tri1 = tri1 / prm->l1;
2623 r_sphere0 = sph[0];
2624 r_sphere1 = sph[1];
2625 r_sphere2 = sph[2];
2626 minIndex = prm->vert;
2627 r->trans = prm->trans;
2628 r_dist = dist;
2629 r_trans = (r->trans = prm->trans);
2630 }
2631 } else {
2632 r->prim = prm;
2633 r->trans = prm->trans;
2634 r->dist = dist;
2635 return (1);
2636 }
2637 }
2638 } else if(trans_shadows) {
2639 if((dist > -kR_SMALL4) &&
2640 ((r_trans > prm->trans) ||
2641 (nearest_shadow && (dist < r_dist)
2642 && (r_trans >= prm->trans)))) {
2643 if(prm->l1 > kR_SMALL4)
2644 r_tri1 = tri1 / prm->l1;
2645 r_sphere0 = sph[0];
2646 r_sphere1 = sph[1];
2647 r_sphere2 = sph[2];
2648 minIndex = prm->vert;
2649 r->trans = prm->trans;
2650 r_dist = dist;
2651 r_trans = (r->trans = prm->trans);
2652 }
2653 }
2654 }
2655 }
2656 }
2657 break;
2658 case cPrimCylinder:
2659 if(ZLineToSphereCapped(r->base, BI->Vertex + i * 3,
2660 BI->Normal + BI->Vert2Normal[i] * 3,
2661 BI->Radius[i], prm->l1, sph, &tri1, prm->cap1,
2662 prm->cap2, BI->Precomp + BI->Vert2Normal[i] * 3)) {
2663
2664 oppSq = ZLineClipPoint(r->base, sph, &dist, BI->Radius[i]);
2665 if(oppSq <= BI->Radius2[i]) {
2666 dist = (float) (sqrt1f(dist) - sqrt1f((BI->Radius2[i] - oppSq)));
2667
2668 if(prm->trans == _0) {
2669 if(dist > -kR_SMALL4) {
2670 if(nearest_shadow) {
2671 if(dist < r_dist) {
2672 if(prm->l1 > kR_SMALL4)
2673 r_tri1 = tri1 / prm->l1;
2674 r_sphere0 = sph[0];
2675 r_sphere1 = sph[1];
2676 r_sphere2 = sph[2];
2677 minIndex = prm->vert;
2678 r->trans = prm->trans;
2679 r_dist = dist;
2680 r_trans = (r->trans = prm->trans);
2681 }
2682 } else {
2683 r->prim = prm;
2684 r->trans = prm->trans;
2685 r->dist = dist;
2686 return (1);
2687 }
2688 }
2689 } else if(trans_shadows) {
2690 if((dist > -kR_SMALL4) &&
2691 ((r_trans > prm->trans) ||
2692 (nearest_shadow && (dist < r_dist) && (r_trans >= prm->trans)))) {
2693 if(prm->l1 > kR_SMALL4)
2694 r_tri1 = tri1 / prm->l1;
2695 r_sphere0 = sph[0];
2696 r_sphere1 = sph[1];
2697 r_sphere2 = sph[2];
2698 minIndex = prm->vert;
2699 r->trans = prm->trans;
2700 r_dist = dist;
2701 r_trans = (r->trans = prm->trans);
2702 }
2703 }
2704 }
2705 }
2706 break;
2707
2708 case cPrimSausage:
2709 if(ZLineToSphere
2710 (r->base, BI->Vertex + i * 3, BI->Normal + BI->Vert2Normal[i] * 3,
2711 BI->Radius[i], prm->l1, sph, &tri1,
2712 BI->Precomp + BI->Vert2Normal[i] * 3)) {
2713 oppSq = ZLineClipPoint(r->base, sph, &dist, BI->Radius[i]);
2714 if(oppSq <= BI->Radius2[i]) {
2715 dist = (float) (sqrt1f(dist) - sqrt1f((BI->Radius2[i] - oppSq)));
2716
2717 if(prm->trans == _0) {
2718 if(dist > -kR_SMALL4) {
2719 if(nearest_shadow) {
2720 if(dist < r_dist) {
2721 if(prm->l1 > kR_SMALL4)
2722 r_tri1 = tri1 / prm->l1;
2723 r_sphere0 = sph[0];
2724 r_sphere1 = sph[1];
2725 r_sphere2 = sph[2];
2726 minIndex = prm->vert;
2727 r_dist = dist;
2728 r_trans = (r->trans = prm->trans);
2729 }
2730 } else {
2731 r->prim = prm;
2732 r->trans = prm->trans;
2733 r->dist = dist;
2734 return (1);
2735 }
2736 }
2737 } else if(trans_shadows) {
2738 if((dist > -kR_SMALL4) &&
2739 ((r_trans > prm->trans) ||
2740 (nearest_shadow && (dist < r_dist) && (r_trans >= prm->trans)))) {
2741 if(prm->l1 > kR_SMALL4)
2742 r_tri1 = tri1 / prm->l1;
2743
2744 r_sphere0 = sph[0];
2745 r_sphere1 = sph[1];
2746 r_sphere2 = sph[2];
2747 minIndex = prm->vert;
2748 r_dist = dist;
2749 r_trans = (r->trans = prm->trans);
2750 }
2751 }
2752 }
2753 }
2754 break;
2755 } /* end of switch */
2756 }
2757 /* end of if */
2758 i = ii;
2759 } /* end of while */
2760 }
2761
2762 /* and of course stop when we hit the edge of the map */
2763
2764 if(local_iflag)
2765 break;
2766
2767 /* we've processed all primitives associated with this voxel,
2768 so if an intersection has been found which occurs in front of
2769 the next voxel, then we can stop */
2770
2771 /* this optimization invalid for transparent surfaces
2772
2773 if( minIndex > -1 )
2774 {
2775 int aa,bb,cc;
2776
2777 vt[2] = r->base[2] - r_dist;
2778 MapLocus(BI->Map,vt,&aa,&bb,&cc);
2779 if(cc > c)
2780 break;
2781 }
2782 */
2783
2784 c--;
2785 xxtmp--;
2786
2787 } /* end of while */
2788
2789 if(minIndex > -1) {
2790 r_prim = BC->prim + vert2prim[minIndex];
2791
2792 if((r_prim->type == cPrimSphere) || (r_prim->type == cPrimEllipsoid)) {
2793 const float *vv = BI->Vertex + minIndex * 3;
2794 r_sphere0 = vv[0];
2795 r_sphere1 = vv[1];
2796 r_sphere2 = vv[2];
2797 }
2798 }
2799
2800 BC->interior_flag = local_iflag;
2801 r->tri1 = r_tri1;
2802 r->tri2 = r_tri2;
2803 r->prim = r_prim;
2804 r->dist = r_dist;
2805 r->trans = r_trans;
2806 r->sphere[0] = r_sphere0;
2807 r->sphere[1] = r_sphere1;
2808 r->sphere[2] = r_sphere2;
2809 return (minIndex);
2810 } /* end of if */
2811 BC->interior_flag = local_iflag;
2812 return (-1);
2813 }
2814
2815 /*========================================================================*/
BasisMakeMap(CBasis * I,int * vert2prim,CPrimitive * prim,int n_prim,float * volume,int group_id,int block_base,int perspective,float front,float size_hint)2816 int BasisMakeMap(CBasis * I, int *vert2prim, CPrimitive * prim, int n_prim,
2817 float *volume,
2818 int group_id, int block_base,
2819 int perspective, float front, float size_hint)
2820 {
2821 float *v;
2822 float ll;
2823 CPrimitive *prm;
2824 int i;
2825 int *tempRef = NULL;
2826 int n = 0, h, q, x, y, z, j, k, l, e;
2827 int extra_vert = 0;
2828 float p[3], dd[3], *d1, *d2, vd[3], cx[3], cy[3];
2829 float *tempVertex = NULL;
2830 float xs, ys;
2831 int remapMode = true; /* remap mode means that some objects will span more
2832 * than one voxel, so we have to worry about populating
2833 * those voxels and also about eliminating duplicates
2834 * when traversing the neighbor lists */
2835 float min[3], max[3], extent[6];
2836 float sep;
2837 float diagonal[3];
2838 float l1, l2;
2839 float bh, ch;
2840 int n_voxel;
2841 int ok = true;
2842 const float _0 = 0.0;
2843 const float _p5 = 0.5;
2844
2845 PRINTFD(I->G, FB_Ray)
2846 " BasisMakeMap: I->NVertex %d [(%8.3f, %8.3f, %8.3f),...]\n", I->NVertex,
2847 I->Vertex[0], I->Vertex[1], I->Vertex[2]
2848 ENDFD;
2849
2850 sep = I->MinVoxel;
2851 if(sep == _0) {
2852 remapMode = false;
2853 sep = I->MaxRadius; /* also will imply no remapping of vertices */
2854 }
2855 /* we need to get a sense of the actual size in order to avoid sep being too small */
2856
2857 v = I->Vertex;
2858
2859 min[0] = max[0] = v[0];
2860 min[1] = max[1] = v[1];
2861 min[2] = max[2] = v[2];
2862
2863 v += 3;
2864
2865 {
2866 int a;
2867 for(a = 1; a < I->NVertex; a++) {
2868 if(min[0] > v[0])
2869 min[0] = v[0];
2870 if(max[0] < v[0])
2871 max[0] = v[0];
2872
2873 if(min[1] > v[1])
2874 min[1] = v[1];
2875 if(max[1] < v[1])
2876 max[1] = v[1];
2877
2878 if(min[2] > v[2])
2879 min[2] = v[2];
2880 if(max[2] < v[2])
2881 max[2] = v[2];
2882
2883 v += 3;
2884 }
2885 }
2886
2887 if(volume) {
2888
2889 /*
2890 if(min[0] > volume[0]) min[0] = volume[0];
2891 if(max[0] < volume[1]) max[0] = volume[1];
2892 if(min[1] > volume[2]) min[1] = volume[2];
2893 if(max[1] < volume[3]) max[1] = volume[3];
2894 if(min[2] > (-volume[5])) min[2] = (-volume[5]);
2895 if(max[2] < (-volume[4])) max[2] = (-volume[4]);
2896 */
2897
2898 if(min[0] < volume[0])
2899 min[0] = volume[0];
2900 if(max[0] > volume[1])
2901 max[0] = volume[1];
2902 if(min[1] < volume[2])
2903 min[1] = volume[2];
2904 if(max[1] > volume[3])
2905 max[1] = volume[3];
2906 if(min[2] > (-volume[5]))
2907 min[2] = (-volume[5]);
2908 if(max[2] < (-volume[4]))
2909 max[2] = (-volume[4]);
2910
2911 PRINTFB(I->G, FB_Ray, FB_Debugging)
2912 " BasisMakeMap: (%8.3f,%8.3f),(%8.3f,%8.3f),(%8.3f,%8.3f)\n",
2913 volume[0], volume[1], volume[2], volume[3], volume[4], volume[5]
2914 ENDFB(I->G);
2915 }
2916
2917 /* don't break up space unnecessarily if we only have a few vertices... */
2918
2919 if(I->NVertex) {
2920 l1 = (float) fabs(max[0] - min[0]);
2921 l2 = (float) fabs(max[1] - min[1]);
2922
2923 if(l1 < l2)
2924 l1 = l2;
2925
2926 l2 = (float) fabs(max[2] - min[2]);
2927
2928 if(l1 < l2)
2929 l1 = l2;
2930
2931 if(l1 < kR_SMALL4)
2932 l1 = 100.0;
2933
2934 if(I->NVertex < (l1 / sep))
2935 sep = (l1 / I->NVertex);
2936 }
2937
2938 if(Feedback(I->G, FB_Ray, FB_Debugging)) {
2939 dump3f(min, " BasisMakeMap: min");
2940 dump3f(max, " BasisMakeMap: max");
2941 dump3f(I->Vertex, " BasisMakeMap: I->Vertex");
2942 fflush(stdout);
2943 }
2944
2945 sep = MapGetSeparation(I->G, sep, max, min, diagonal); /* this needs to be a minimum
2946 * estimate of the actual value */
2947 /* here we have to carry out a complicated work-around in order to
2948 * efficiently encode our primitives into the map in a way that doesn't
2949 * require expanding the map cutoff to the size of the largest object*/
2950 if(remapMode) {
2951 int a, b, c;
2952
2953 if(sep < size_hint) /* this keeps us from wasting time & memory on unnecessary subdivision */
2954 sep = size_hint;
2955
2956 {
2957 int *vert2prim_a = vert2prim;
2958 for(a = 0; a < I->NVertex; a++) {
2959 prm = prim + *(vert2prim_a++);
2960
2961 switch (prm->type) {
2962 case cPrimTriangle:
2963 case cPrimCharacter:
2964 if(a == prm->vert) { /* only do this calculation for one of the three vertices */
2965 l1 = (float) length3f(I->Precomp + I->Vert2Normal[a] * 3);
2966 l2 = (float) length3f(I->Precomp + I->Vert2Normal[a] * 3 + 3);
2967 if((l1 >= sep) || (l2 >= sep)) {
2968 b = (int) ceil(l1 / sep) + 1;
2969 c = (int) ceil(l2 / sep) + 1;
2970 extra_vert += 4 * b * c;
2971 }
2972 }
2973 break;
2974
2975 case cPrimCone:
2976 case cPrimCylinder:
2977 case cPrimSausage:
2978 if((prm->l1 + 2 * prm->r1) >= sep) {
2979 q = ((int) (2 * (floor(prm->r1 / sep) + 1))) + 1;
2980 q = q * q * ((int) ceil((prm->l1 + 2 * prm->r1) / sep) + 2);
2981 extra_vert += q;
2982 }
2983 break;
2984 case cPrimEllipsoid:
2985 case cPrimSphere:
2986 if(prm->r1 >= sep) {
2987 b = (int) (2 * floor(prm->r1 / sep) + 1);
2988 extra_vert += (b * b * b);
2989 }
2990 break;
2991 }
2992 } /* for */
2993 }
2994
2995 extra_vert += I->NVertex;
2996 if (ok)
2997 tempVertex =
2998 CacheAlloc(I->G, float, extra_vert * 3, group_id, cCache_basis_tempVertex);
2999 CHECKOK(ok, tempVertex);
3000 if (ok)
3001 tempRef = CacheAlloc(I->G, int, extra_vert, group_id, cCache_basis_tempRef);
3002 CHECKOK(ok, tempRef);
3003
3004 ErrChkPtr(I->G, tempVertex); /* can happen if extra vert is unreasonable */
3005 ErrChkPtr(I->G, tempRef);
3006
3007 /* lower indexes->flags, top is ref->lower index */
3008 ok &= !I->G->Interrupt;
3009
3010 if (ok){
3011 float *vv, *d;
3012 int *vert2prim_a = vert2prim;
3013
3014 n = I->NVertex;
3015
3016 v = tempVertex;
3017 vv = I->Vertex;
3018
3019 memcpy(v, vv, n * sizeof(float) * 3);
3020 vv += n * 3;
3021 v += n * 3;
3022
3023 for(a = 0; ok && a < I->NVertex; a++) {
3024
3025 prm = prim + *(vert2prim_a++);
3026
3027 switch (prm->type) {
3028 case cPrimTriangle:
3029 case cPrimCharacter:
3030 if(a == prm->vert) {
3031 /* only do this calculation for one of the three vertices */
3032 d1 = I->Precomp + I->Vert2Normal[a] * 3;
3033 d2 = I->Precomp + I->Vert2Normal[a] * 3 + 3;
3034 vv = I->Vertex + a * 3;
3035 l1 = (float) length3f(d1);
3036 l2 = (float) length3f(d2);
3037 if((l1 >= sep) || (l2 >= sep)) {
3038 b = (int) floor(l1 / sep) + 1;
3039 c = (int) floor(l2 / sep) + 1;
3040 extra_vert += b * c;
3041 bh = (float) (b / 2) + 1;
3042 ch = (float) (c / 2) + 1;
3043
3044 for(x = 0; x < bh; x++) {
3045 const float xb = (float) x / b;
3046
3047 for(y = 0; y < ch; y++) {
3048 const float yc = (float) y / c;
3049
3050 *(v++) = vv[0] + (d1[0] * xb) + (d2[0] * yc);
3051 *(v++) = vv[1] + (d1[1] * xb) + (d2[1] * yc);
3052 *(v++) = vv[2] + (d1[2] * xb) + (d2[2] * yc);
3053
3054 tempRef[n++] = a;
3055 }
3056 }
3057
3058 for(x = 0; x < bh; x++) {
3059 const float xb = (float) x / b;
3060
3061 for(y = 0; y < ch; y++) {
3062 const float yc = (float) y / c;
3063
3064 if((xb + yc) < _p5) {
3065 *(v++) = vv[0] + d1[0] * (_p5 + xb) + (d2[0] * yc);
3066 *(v++) = vv[1] + d1[1] * (_p5 + xb) + (d2[1] * yc);
3067 *(v++) = vv[2] + d1[2] * (_p5 + xb) + (d2[2] * yc);
3068
3069 tempRef[n++] = a;
3070
3071 *(v++) = vv[0] + (d1[0] * xb) + d2[0] * (_p5 + yc);
3072 *(v++) = vv[1] + (d1[1] * xb) + d2[1] * (_p5 + yc);
3073 *(v++) = vv[2] + (d1[2] * xb) + d2[2] * (_p5 + yc);
3074
3075 tempRef[n++] = a;
3076 }
3077 }
3078 }
3079 }
3080 } /* if */
3081 break;
3082 case cPrimCone:
3083 case cPrimCylinder:
3084 case cPrimSausage:
3085
3086 ll = prm->l1 + 2 * prm->r1;
3087
3088 if(ll >= sep) {
3089 d = I->Normal + 3 * I->Vert2Normal[a];
3090 vv = I->Vertex + a * 3;
3091 get_system1f3f(d, cx, cy); /* creates an orthogonal system about d */
3092
3093 p[0] = vv[0] - d[0] * prm->r1;
3094 p[1] = vv[1] - d[1] * prm->r1;
3095 p[2] = vv[2] - d[2] * prm->r1;
3096 dd[0] = d[0] * sep;
3097 dd[1] = d[1] * sep;
3098 dd[2] = d[2] * sep;
3099
3100 q = (int) floor(prm->r1 / sep) + 1;
3101
3102 while(1) {
3103 vd[0] = (p[0] += dd[0]);
3104 vd[1] = (p[1] += dd[1]);
3105 vd[2] = (p[2] += dd[2]);
3106
3107 for(x = -q; x <= q; x++) {
3108 for(y = -q; y <= q; y++) {
3109 xs = x * sep;
3110 ys = y * sep;
3111 *(v++) = vd[0] + xs * cx[0] + ys * cy[0];
3112 *(v++) = vd[1] + xs * cx[1] + ys * cy[1];
3113 *(v++) = vd[2] + xs * cx[2] + ys * cy[2];
3114
3115 tempRef[n++] = a;
3116 }
3117 }
3118
3119 if(ll <= _0)
3120 break;
3121 ll -= sep;
3122 }
3123 }
3124 break;
3125 case cPrimEllipsoid:
3126 case cPrimSphere:
3127 if(prm->r1 >= sep) {
3128 q = (int) floor(prm->r1 / sep);
3129 vv = I->Vertex + a * 3;
3130
3131 for(x = -q; x <= q; x++) {
3132 for(y = -q; y <= q; y++) {
3133 for(z = -q; z <= q; z++) {
3134 *(v++) = vv[0] + x * sep;
3135 *(v++) = vv[1] + y * sep;
3136 *(v++) = vv[2] + z * sep;
3137
3138 tempRef[n++] = a;
3139 }
3140 }
3141 }
3142 }
3143 break;
3144 } /* end of switch */
3145 ok &= !I->G->Interrupt;
3146 }
3147 }
3148 if(n > extra_vert) {
3149 printf("BasisMakeMap: %d>%d\n", n, extra_vert);
3150 ErrFatal(I->G, "BasisMakeMap",
3151 "used too many extra vertices (this indicates a bug)...\n");
3152 }
3153 PRINTFB(I->G, FB_Ray, FB_Blather)
3154 " BasisMakeMap: %d total vertices\n", n ENDFB(I->G);
3155
3156 if (ok){
3157 if(volume) {
3158 v = tempVertex;
3159
3160 min[0] = max[0] = v[0];
3161 min[1] = max[1] = v[1];
3162 min[2] = max[2] = v[2];
3163
3164 v += 3;
3165
3166 if(Feedback(I->G, FB_Ray, FB_Debugging)) {
3167 dump3f(min, " BasisMakeMap: remapped min");
3168 dump3f(max, " BasisMakeMap: remapped max");
3169 fflush(stdout);
3170 }
3171
3172 for(a = 1; a < n; a++) {
3173 if(min[0] > v[0])
3174 min[0] = v[0];
3175 if(max[0] < v[0])
3176 max[0] = v[0];
3177
3178 if(min[1] > v[1])
3179 min[1] = v[1];
3180 if(max[1] < v[1])
3181 max[1] = v[1];
3182
3183 if(min[2] > v[2])
3184 min[2] = v[2];
3185 if(max[2] < v[2])
3186 max[2] = v[2];
3187 v += 3;
3188 }
3189
3190 if(Feedback(I->G, FB_Ray, FB_Debugging)) {
3191 dump3f(min, " BasisMakeMap: remapped min");
3192 dump3f(max, " BasisMakeMap: remapped max");
3193 fflush(stdout);
3194 }
3195 if(min[0] < volume[0])
3196 min[0] = volume[0];
3197 if(max[0] > volume[1])
3198 max[0] = volume[1];
3199 if(min[1] < volume[2])
3200 min[1] = volume[2];
3201 if(max[1] > volume[3])
3202 max[1] = volume[3];
3203
3204 if(min[2] < (-volume[5]))
3205 min[2] = (-volume[5]);
3206 if(max[2] > (-volume[4]))
3207 max[2] = (-volume[4]);
3208
3209 extent[0] = min[0];
3210 extent[1] = max[0];
3211 extent[2] = min[1];
3212 extent[3] = max[1];
3213 extent[4] = min[2];
3214 extent[5] = max[2];
3215 PRINTFB(I->G, FB_Ray, FB_Blather)
3216 " BasisMakeMap: Extent [%8.2f %8.2f] [%8.2f %8.2f] [%8.2f %8.2f]\n",
3217 extent[0], extent[1], extent[2], extent[3], extent[4], extent[5]
3218 ENDFB(I->G);
3219 I->Map = MapNewCached(I->G, -sep, tempVertex, n, extent, group_id, block_base);
3220 CHECKOK(ok, I->Map);
3221 } else {
3222 I->Map = MapNewCached(I->G, sep, tempVertex, n, NULL, group_id, block_base);
3223 CHECKOK(ok, I->Map);
3224 }
3225 }
3226 if (ok)
3227 n_voxel = I->Map->Dim[0] * I->Map->Dim[1] * I->Map->Dim[2];
3228
3229 if (!ok){
3230 } else if(perspective) {
3231
3232 /* this is a new optimization which prevents primitives
3233 contained entirely within a single voxel from spilling over
3234 into neighboring voxels for performance of intersection
3235 checks (NOTE: this currently only helps with triangles, and
3236 is only useful in optimizing between Z layers). The main
3237 impact of this optimization is to reduce the size of the
3238 express list (I->Map->Elist) array by about 3-fold */
3239
3240 int *prm_spanner = NULL;
3241 int *spanner = NULL;
3242 float *v = tempVertex;
3243 int j, k, l, jj, kk, ll;
3244 int nVertex = I->NVertex;
3245 int prm_index;
3246 MapType *map = I->Map;
3247
3248 prm_spanner = pymol::calloc<int>(n_prim);
3249 CHECKOK(ok, prm_spanner);
3250 if (ok)
3251 spanner = pymol::calloc<int>(n);
3252 CHECKOK(ok, spanner);
3253
3254 /* figure out which primitives span more than one voxel */
3255 for(a = 0; ok && a < n; a++) {
3256 if(a < nVertex)
3257 prm_index = vert2prim[a];
3258 else
3259 prm_index = vert2prim[tempRef[a]];
3260 if(!prm_spanner[prm_index]) {
3261 prm = prim + prm_index;
3262 {
3263 float *vv = prm->vert * 3 + I->Vertex;
3264 switch (prm->type) {
3265 case cPrimTriangle:
3266 case cPrimCharacter:
3267 MapLocus(map, v, &j, &k, &l);
3268 MapLocus(map, vv, &jj, &kk, &ll);
3269 if((j != jj) || (k != kk) || (l != ll)) {
3270 prm_spanner[prm_index] = 1;
3271 }
3272 break;
3273 default: /* currently we aren't optimizing other primitives */
3274 prm_spanner[prm_index] = 1;
3275 break;
3276 }
3277 }
3278 }
3279 v += 3;
3280 ok &= !I->G->Interrupt;
3281 }
3282 /* now flag associated vertices */
3283 for(a = 0; ok && a < n; a++) {
3284 if(a < nVertex)
3285 prm_index = vert2prim[a];
3286 else
3287 prm_index = vert2prim[tempRef[a]];
3288 spanner[a] = prm_spanner[prm_index];
3289 ok &= !I->G->Interrupt;
3290 }
3291 /* and do the optimized expansion */
3292 ok &= MapSetupExpressPerp(I->Map, tempVertex, front, n, true, spanner);
3293 FreeP(spanner);
3294 FreeP(prm_spanner);
3295 } else if(n_voxel < (3 * n)) {
3296 ok &= MapSetupExpressXY(I->Map, n, true);
3297 } else {
3298 ok &= MapSetupExpressXYVert(I->Map, tempVertex, n, true);
3299 }
3300
3301 if (ok) {
3302 MapType *map = I->Map;
3303 int *sp, *ip, *ip0, ii;
3304 int *elist = map->EList, *ehead = map->EHead;
3305 int *elist_new = elist, *ehead_new = ehead;
3306 int newelem = 0, neelem = -map->NEElem;
3307 int i_nVertex = I->NVertex;
3308 const int iMin0 = map->iMin[0];
3309 const int iMin1 = map->iMin[1];
3310 const int iMin2 = map->iMin[2];
3311 const int iMax0 = map->iMax[0];
3312 const int iMax1 = map->iMax[1];
3313 const int iMax2 = map->iMax[2];
3314
3315 /* now do a filter-reassignment pass to remap fake vertices
3316 to the original line vertex while deleting duplicate entries */
3317
3318 memset(tempRef, 0, sizeof(int) * i_nVertex);
3319
3320 if(n_voxel < (3 * n)) { /* faster to traverse the entire map */
3321 int *start;
3322 for(a = iMin0; a <= iMax0; a++) {
3323 for(b = iMin1; b <= iMax1; b++) {
3324 for(c = iMin2; c <= iMax2; c++) {
3325 start = MapEStart(map, a, b, c);
3326 h = *start;
3327 if((h < 0) && (h >= neelem)) {
3328 sp = elist - h;
3329 *(start) = -h; /* flip sign */
3330 i = *(sp++);
3331 if(ehead_new != ehead) {
3332 ehead_new[(start - ehead) - 1] = newelem;
3333 ip = ip0 = (elist_new + newelem);
3334 } else {
3335 ip = ip0 = (sp - 1);
3336 }
3337
3338 while(i >= 0) {
3339 int ii = *(sp++);
3340 if(i >= i_nVertex) /* reference -- remap */
3341 i = tempRef[i];
3342 if(!tempRef[i]) { /*eliminate duplicates */
3343 tempRef[i] = 1;
3344 *(ip++) = i;
3345 }
3346 i = ii;
3347 }
3348
3349 *(ip++) = -1; /* terminate list */
3350 newelem += (ip - ip0);
3351
3352 /* now reset flags efficiently */
3353 i = *(ip0++);
3354 while(i >= 0) { /* unroll */
3355 ii = *(ip0++);
3356 tempRef[i] = 0;
3357 if((i = ii) >= 0) {
3358 ii = *(ip0++);
3359 tempRef[i] = 0;
3360 if((i = ii) >= 0) {
3361 ii = *(ip0++);
3362 tempRef[i] = 0;
3363 i = ii;
3364 }
3365 }
3366 }
3367 }
3368 } /* for c */
3369 } /* for b */
3370 } /* for a */
3371 } else {
3372 int *start;
3373 int jm1, jp1, km1, kp1, lm1, lp1;
3374 MapType *map = I->Map;
3375
3376 v = tempVertex;
3377
3378 for(e = 0; e < n; e++) {
3379
3380 MapLocus(map, v, &j, &k, &l);
3381 jm1 = j - 1;
3382 jp1 = j + 1;
3383 km1 = k - 1;
3384 kp1 = k + 1;
3385 lm1 = l - 1;
3386 lp1 = l + 1;
3387
3388 if(perspective) { /* for normal maps */
3389
3390 int *iPtr1 =
3391 map->EHead + ((j - 1) * map->D1D2) + ((k - 1) * map->Dim[2]) + (l - 1);
3392 for(a = jm1; a <= jp1; a++) {
3393 int *iPtr2 = iPtr1;
3394 if((a >= iMin0) && (a <= iMax0)) {
3395 for(b = km1; b <= kp1; b++) {
3396 int *iPtr3 = iPtr2;
3397 if((b >= iMin1) && (b <= iMax1)) {
3398 for(c = lm1; c <= lp1; c++) {
3399 if((c >= iMin2) && (c <= iMax2)) {
3400
3401 start = iPtr3;
3402 /*start = MapEStart(map,a,b,c); */
3403 h = *start;
3404 if((h < 0) && (h >= neelem)) {
3405 sp = elist - h;
3406 (*start) = -h; /* no repeat visits */
3407 i = *(sp++);
3408 if(ehead_new != ehead) {
3409 ehead_new[(start - ehead)] = newelem;
3410 ip = ip0 = (elist_new + newelem);
3411 } else {
3412 ip = ip0 = (sp - 1);
3413 }
3414
3415 while(i >= 0) {
3416 int ii = *(sp++);
3417 if(i >= i_nVertex)
3418 i = tempRef[i];
3419 if(!tempRef[i]) { /*eliminate duplicates */
3420 tempRef[i] = 1;
3421 *(ip++) = i;
3422 }
3423 i = ii;
3424 }
3425
3426 *(ip++) = -1; /* terminate list */
3427 newelem += (ip - ip0);
3428
3429 /* now reset flags efficiently */
3430 i = *(ip0++);
3431
3432 while(i >= 0) { /* unroll */
3433 ii = *(ip0++);
3434 tempRef[i] = 0;
3435 if((i = ii) >= 0) {
3436 ii = *(ip0++);
3437 tempRef[i] = 0;
3438 if((i = ii) >= 0) {
3439 ii = *(ip0++);
3440 tempRef[i] = 0;
3441 i = ii;
3442 }
3443 }
3444 }
3445 } /* h > 0 */
3446 }
3447 iPtr3++;
3448 }
3449 }
3450 iPtr2 += map->Dim[2];
3451 }
3452 }
3453 iPtr1 += map->D1D2;
3454 }
3455 } else { /* for XY maps... */
3456
3457 c = l;
3458 {
3459 int *iPtr1 =
3460 ehead + ((j - 1) * map->D1D2) + ((k - 1) * map->Dim[2]) + c;
3461 if((c >= iMin2) && (c <= iMax2)) {
3462
3463 for(a = jm1; a <= jp1; a++) {
3464 int *iPtr2 = iPtr1;
3465 if((a >= iMin0) && (a <= iMax0)) {
3466
3467 for(b = km1; b <= kp1; b++) {
3468 if((b >= iMin1) && (b <= iMax1)) {
3469 start = iPtr2;
3470 /*start = MapEStart(map,a,b,c); */
3471 h = *start;
3472 if((h < 0) && (h >= neelem)) {
3473
3474 sp = elist - h;
3475 (*start) = -h; /* no repeat visits */
3476 i = *(sp++);
3477 if(ehead_new != ehead) {
3478 ehead_new[(start - ehead)] = newelem;
3479 ip = ip0 = (elist_new + newelem);
3480 } else {
3481 ip = ip0 = sp - 1;
3482 }
3483
3484 while(i >= 0) {
3485 int ii = *(sp++);
3486
3487 if(i >= i_nVertex)
3488 i = tempRef[i];
3489
3490 if(!tempRef[i]) { /*eliminate duplicates */
3491 tempRef[i] = 1;
3492 *(ip++) = i;
3493 }
3494 i = ii;
3495 }
3496
3497 *(ip++) = -1; /* terminate list */
3498 newelem += (ip - ip0);
3499
3500 /* now reset flags efficiently */
3501 i = *(ip0++);
3502 while(i >= 0) { /* unroll */
3503 ii = *(ip0++);
3504 tempRef[i] = 0;
3505 if((i = ii) >= 0) {
3506 ii = *(ip0++);
3507 tempRef[i] = 0;
3508 if((i = ii) >= 0) {
3509 ii = *(ip0++);
3510 tempRef[i] = 0;
3511 i = ii;
3512 }
3513 }
3514 }
3515 } /* h > 0 */
3516 }
3517 iPtr2 += map->Dim[2];
3518 } /* for b */
3519 }
3520 iPtr1 += map->D1D2;
3521 } /* for a */
3522 }
3523 }
3524 }
3525 v += 3; /* happens for EVERY e! */
3526 } /* for e */
3527 }
3528 if(ehead != ehead_new) {
3529 CacheFreeP(I->G, map->EHead, group_id, block_base + cCache_map_ehead_offset,
3530 false);
3531 VLACacheFreeP(I->G, map->EList, group_id, block_base + cCache_map_elist_offset,
3532 false);
3533 map->EList = elist_new;
3534 map->EHead = ehead_new;
3535 map->NEElem = newelem;
3536 MemoryCacheReplaceBlock(I->G, group_id, block_base + cCache_map_ehead_new_offset,
3537 block_base + cCache_map_ehead_offset);
3538 MemoryCacheReplaceBlock(I->G, group_id, block_base + cCache_map_elist_new_offset,
3539 block_base + cCache_map_elist_offset);
3540 VLACacheSize(G, map->EList, int, map->NEElem, group_id,
3541 block_base + cCache_map_elist_offset);
3542 }
3543 }
3544
3545 CacheFreeP(I->G, tempVertex, group_id, cCache_basis_tempVertex, false);
3546 CacheFreeP(I->G, tempRef, group_id, cCache_basis_tempRef, false);
3547
3548 } else {
3549 /* simple sphere mode */
3550 I->Map = MapNewCached(I->G, -sep, I->Vertex, I->NVertex, NULL, group_id, block_base);
3551 CHECKOK(ok, I->Map);
3552 if (ok){
3553 if(perspective) {
3554 ok &= MapSetupExpressPerp(I->Map, I->Vertex, front, I->NVertex, false, NULL);
3555 } else {
3556 ok &= MapSetupExpressXYVert(I->Map, I->Vertex, I->NVertex, false);
3557 }
3558 }
3559 }
3560 return ok;
3561 }
3562
3563
3564 /*========================================================================*/
BasisInit(PyMOLGlobals * G,CBasis * I,int group_id)3565 int BasisInit(PyMOLGlobals * G, CBasis * I, int group_id)
3566 {
3567 int ok = true;
3568 I->G = G;
3569 I->Radius = NULL;
3570 I->Radius2 = NULL;
3571 I->Normal = NULL;
3572 I->Vert2Normal = NULL;
3573 I->Precomp = NULL;
3574 I->Vertex = VLACacheAlloc(I->G, float, 1, group_id, cCache_basis_vertex);
3575 CHECKOK(ok, I->Vertex);
3576 if (ok)
3577 I->Radius = VLACacheAlloc(I->G, float, 1, group_id, cCache_basis_radius);
3578 CHECKOK(ok, I->Radius);
3579 if (ok)
3580 I->Radius2 = VLACacheAlloc(I->G, float, 1, group_id, cCache_basis_radius2);
3581 CHECKOK(ok, I->Radius2);
3582 if (ok)
3583 I->Normal = VLACacheAlloc(I->G, float, 1, group_id, cCache_basis_normal);
3584 CHECKOK(ok, I->Normal);
3585 if (ok)
3586 I->Vert2Normal = VLACacheAlloc(I->G, int, 1, group_id, cCache_basis_vert2normal);
3587 CHECKOK(ok, I->Vert2Normal);
3588 if (ok)
3589 I->Precomp = VLACacheAlloc(I->G, float, 1, group_id, cCache_basis_precomp);
3590 CHECKOK(ok, I->Precomp);
3591 I->Map = NULL;
3592 I->NVertex = 0;
3593 I->NNormal = 0;
3594 return ok;
3595 }
3596
3597
3598 /*========================================================================*/
BasisFinish(CBasis * I,int group_id)3599 void BasisFinish(CBasis * I, int group_id)
3600 {
3601 if(I->Map) {
3602 MapFree(I->Map);
3603 I->Map = NULL;
3604 }
3605 VLACacheFreeP(I->G, I->Radius2, group_id, cCache_basis_radius2, false);
3606 VLACacheFreeP(I->G, I->Radius, group_id, cCache_basis_radius, false);
3607 VLACacheFreeP(I->G, I->Vertex, group_id, cCache_basis_vertex, false);
3608 VLACacheFreeP(I->G, I->Vert2Normal, group_id, cCache_basis_vert2normal, false);
3609 VLACacheFreeP(I->G, I->Normal, group_id, cCache_basis_normal, false);
3610 VLACacheFreeP(I->G, I->Precomp, group_id, cCache_basis_precomp, false);
3611 I->Vertex = NULL;
3612 }
3613
3614
3615 /*========================================================================*/
3616
BasisTrianglePrecompute(float * v0,float * v1,float * v2,float * pre)3617 void BasisTrianglePrecompute(float *v0, float *v1, float *v2, float *pre)
3618 {
3619 float det;
3620
3621 subtract3f(v1, v0, pre);
3622 subtract3f(v2, v0, pre + 3);
3623
3624 det = pre[0] * pre[4] - pre[1] * pre[3];
3625
3626 if(fabs(det) < EPSILON)
3627 *(pre + 6) = 0.0F;
3628 else {
3629 *(pre + 6) = 1.0F;
3630 *(pre + 7) = 1.0F / det;
3631 }
3632 }
3633
BasisTrianglePrecomputePerspective(float * v0,float * v1,float * v2,float * pre)3634 void BasisTrianglePrecomputePerspective(float *v0, float *v1, float *v2, float *pre)
3635 {
3636 subtract3f(v1, v0, pre);
3637 subtract3f(v2, v0, pre + 3);
3638 }
3639
BasisCylinderSausagePrecompute(float * dir,float * pre)3640 void BasisCylinderSausagePrecompute(float *dir, float *pre)
3641 {
3642 float ln = (float) (1.0f / sqrt1f(dir[1] * dir[1] + dir[0] * dir[0]));
3643 pre[0] = dir[1] * ln;
3644 pre[1] = -dir[0] * ln;
3645 }
3646