1 //******************************************************************************
2 ///
3 /// @file core/shape/blob.cpp
4 ///
5 /// Implementation of the blob geometric primitive.
6 ///
7 /// @author Alexander Enzmann (original file)
8 /// @author Dieter Bayer (modifications and enhancements)
9 ///
10 /// @copyright
11 /// @parblock
12 ///
13 /// Persistence of Vision Ray Tracer ('POV-Ray') version 3.8.
14 /// Copyright 1991-2018 Persistence of Vision Raytracer Pty. Ltd.
15 ///
16 /// POV-Ray is free software: you can redistribute it and/or modify
17 /// it under the terms of the GNU Affero General Public License as
18 /// published by the Free Software Foundation, either version 3 of the
19 /// License, or (at your option) any later version.
20 ///
21 /// POV-Ray is distributed in the hope that it will be useful,
22 /// but WITHOUT ANY WARRANTY; without even the implied warranty of
23 /// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 /// GNU Affero General Public License for more details.
25 ///
26 /// You should have received a copy of the GNU Affero General Public License
27 /// along with this program. If not, see <http://www.gnu.org/licenses/>.
28 ///
29 /// ----------------------------------------------------------------------------
30 ///
31 /// POV-Ray is based on the popular DKB raytracer version 2.12.
32 /// DKBTrace was originally written by David K. Buck.
33 /// DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
34 ///
35 /// @endparblock
36 ///
37 //******************************************************************************
38
39 /****************************************************************************
40 *
41 * Explanation:
42 *
43 * -
44 *
45 * Syntax:
46 *
47 * blob
48 * {
49 * threshold THRESHOLD_VALUE
50 *
51 * component STRENGTH, RADIUS, <CENTER>
52 *
53 * sphere { <CENTER>, RADIUS, [strength] STRENGTH
54 * [ translate <VECTOR> ]
55 * [ rotate <VECTOR> ]
56 * [ scale <VECTOR> ]
57 * [ finish { ... } ]
58 * [ pigment { ... } ]
59 * [ tnormal { ... } ]
60 * [ texture { ... } ]
61 * }
62 *
63 * cylinder { <END1>, <END2>, RADIUS, [strength] STRENGTH
64 * [ translate <VECTOR> ]
65 * [ rotate <VECTOR> ]
66 * [ scale <VECTOR> ]
67 * [ finish { ... } ]
68 * [ pigment { ... } ]
69 * [ tnormal { ... } ]
70 * [ texture { ... } ]
71 * }
72 *
73 * [ sturm ]
74 * [ hierarchy FLAG ]
75 * }
76 *
77 * ---
78 *
79 * Jul 1994 : Most functions rewritten, bounding hierarchy added. [DB]
80 *
81 * Aug 1994 : Cylindrical blobs added. [DB]
82 *
83 * Sep 1994 : Multi-texturing added (each component can have its own texture).
84 * Translation, rotation and scaling of each component added. [DB]
85 *
86 * Oct 1994 : Adopted the method for the bounding slab creation to build the
87 * bounding sphere hierarchy of the blob to get a much better
88 * hierarchy. Improved bounding sphere calculation for tighter
89 * bounds. [DB]
90 *
91 * Dec 1994 : Added code for dynamic blob queue allocation. [DB]
92 *
93 * Feb 1995 : Moved bounding sphere stuff into a seperate file. [DB]
94 *
95 *****************************************************************************/
96
97 // Unit header file must be the first file included within POV-Ray *.cpp files (pulls in config)
98 #include "core/shape/blob.h"
99
100 #include <algorithm>
101
102 #include "base/pov_err.h"
103
104 #include "core/bounding/boundingbox.h"
105 #include "core/bounding/boundingsphere.h"
106 #include "core/material/texture.h"
107 #include "core/math/matrix.h"
108 #include "core/math/polynomialsolver.h"
109 #include "core/render/ray.h"
110 #include "core/scene/tracethreaddata.h"
111
112 // this must be the last file included
113 #include "base/povdebug.h"
114
115 namespace pov
116 {
117
118 /*****************************************************************************
119 * Local preprocessor defines
120 ******************************************************************************/
121
122 /* Minimal intersection depth for a valid intersection. */
123 const DBL DEPTH_TOLERANCE = 1.0e-2;
124
125 /* Tolerance for inside test. */
126 const DBL INSIDE_TOLERANCE = 1.0e-6;
127
128 /* Ray enters/exits a component. */
129 const int ENTERING = 0;
130 const int EXITING = BLOB_ENTER_EXIT_FLAG;
131
132
133 /*****************************************************************************
134 *
135 * FUNCTION
136 *
137 * All_Blob_Intersections
138 *
139 * INPUT
140 *
141 * Object - Object
142 * Ray - Ray
143 *
144 * OUTPUT
145 *
146 * Depth_Stack - Intersection stack
147 *
148 * RETURNS
149 *
150 * int - true, if a intersection was found
151 *
152 * AUTHOR
153 *
154 * Alexander Enzmann
155 *
156 * DESCRIPTION
157 *
158 * Generate intervals of influence for each component. After these
159 * are made, determine their aggregate effect on the ray. As the
160 * individual intervals are checked, a quartic is generated
161 * that represents the density at a particular point on the ray.
162 *
163 * Explanation for spherical components:
164 *
165 * After making the substitutions in MakeBlob, there is a formula
166 * for each component that has the form:
167 *
168 * c0 * r^4 + c1 * r^2 + c2.
169 *
170 * In order to determine the influence on the ray of all of the
171 * individual components, we start by determining the distance
172 * from any point on the ray to the specified point. This can
173 * be found using the pythagorean theorem, using C as the center
174 * of this component, P as the start of the ray, and D as the
175 * direction of travel of the ray:
176 *
177 * r^2 = (t * D + P - C) . (t * D + P - C)
178 *
179 * we insert this equation for each appearance of r^2 in the
180 * components' formula, giving:
181 *
182 * r^2 = D.D t^2 + 2 t D . (P - C) + (P - C) . (P - C)
183 *
184 * Since the direction vector has been normalized, D.D = 1.
185 * Using the substitutions:
186 *
187 * t0 = (P - C) . (P - C),
188 * t1 = D . (P - C)
189 *
190 * We can write the formula as:
191 *
192 * r^2 = t0 + 2 t t1 + t^2
193 *
194 * Taking r^2 and substituting into the formula for this component
195 * of the Blob we get the formula:
196 *
197 * density = c0 * (r^2)^2 + c1 * r^2 + c2,
198 *
199 * or:
200 *
201 * density = c0 * (t0 + 2 t t1 + t^2)^2 +
202 * c1 * (t0 + 2 t t1 + t^2) +
203 * c2
204 *
205 * Expanding terms and collecting with respect to "t" gives:
206 *
207 * t^4 * c0 +
208 * t^3 * 4 c0 t1 +
209 * t^2 * (c1 + 2 * c0 t0 + 4 c0 t1^2)
210 * t * 2 (c1 t1 + 2 c0 t0 t1) +
211 * c2 + c1*t0 + c0*t0^2
212 *
213 * This formula can now be solved for "t" by any of the quartic
214 * root solvers that are available.
215 *
216 * CHANGES
217 *
218 * Jul 1994 : Added code for cylindrical and ellipsoidical blobs. [DB]
219 *
220 * Oct 1994 : Added code to convert polynomial into a bezier curve for
221 * a quick test if an intersection exists in an interval. [DB]
222 *
223 * Sep 1995 : Added code to avoid numerical problems with distant blobs. [DB]
224 *
225 ******************************************************************************/
226
All_Intersections(const Ray & ray,IStack & Depth_Stack,TraceThreadData * Thread)227 bool Blob::All_Intersections(const Ray& ray, IStack& Depth_Stack, TraceThreadData *Thread)
228 {
229 int i, j, cnt;
230 int root_count, in_flag;
231 int Intersection_Found = false;
232 DBL t0, t1, t2, c0, c1, c2, dist, len, start_dist;
233 DBL *fcoeffs;
234 DBL coeffs[5];
235 DBL roots[4];
236 Vector3d P, D, V1, PP, DD;
237 Vector3d IPoint;
238 const Blob_Element *Element;
239 DBL newcoeffs[5], dk[5];
240 DBL max_bound;
241 DBL l, w;
242 DBL depthTolerance = (ray.IsSubsurfaceRay()? 0 : DEPTH_TOLERANCE);
243
244 Thread->Stats()[Ray_Blob_Tests]++;
245
246 /* Transform the ray into blob space. */
247
248 if (Trans != nullptr)
249 {
250 MInvTransPoint(P, ray.Origin, Trans);
251 MInvTransDirection(D, ray.Direction, Trans);
252
253 len = D.length();
254 D /= len;
255 }
256 else
257 {
258 P = ray.Origin;
259 D = ray.Direction;
260
261 len = 1.0;
262 }
263
264 /* Get the intervals along the ray where each component has an effect. */
265 Blob_Interval_Struct* intervals = Thread->Blob_Intervals;
266 if ((cnt = determine_influences(P, D, depthTolerance, intervals, Thread)) == 0)
267 {
268 /* Ray doesn't hit any of the component elements. */
269 return (false);
270 }
271
272 /* To avoid numerical problems we start at the first interval. */
273
274 if ((start_dist = intervals[0].bound) < SMALL_TOLERANCE)
275 {
276 start_dist = 0.0;
277 }
278
279 for (i = 0; i < cnt; i++)
280 {
281 intervals[i].bound -= start_dist;
282 }
283
284 P += start_dist * D;
285
286 max_bound = intervals[0].bound;
287 for (i = 0; i < cnt; i++)
288 {
289 if (intervals[i].bound > max_bound)
290 max_bound = intervals[i].bound;
291 }
292
293 /* Get the new starting point. */
294
295 if (max_bound != 0)
296 {
297 D *= max_bound;
298
299 for (i = 0; i < cnt; i++)
300 intervals[i].bound /= max_bound;
301 }
302 else
303 max_bound = 1;
304
305 /* Clear out the coefficients. */
306
307 coeffs[0] =
308 coeffs[1] =
309 coeffs[2] =
310 coeffs[3] = 0.0;
311 coeffs[4] = -Data->Threshold;
312
313 /*
314 * Step through the list of intersection points, adding the
315 * influence of each component as it appears.
316 */
317
318 fcoeffs = nullptr;
319
320 for (i = in_flag = 0; i < cnt; i++)
321 {
322 if ((intervals[i].type & 1) == ENTERING)
323 {
324 /*
325 * Something is just starting to influence the ray, so calculate
326 * its coefficients and add them into the pot.
327 */
328
329 in_flag++;
330
331 Element = intervals[i].Element;
332 fcoeffs = &Thread->Blob_Coefficients[Element->index * 5];
333
334 switch (Element->Type)
335 {
336 case BLOB_SPHERE:
337
338 V1 = P - Element->O;
339
340 t0 = V1.lengthSqr();
341 t1 = dot(V1, D);
342 t2 = max_bound * max_bound;
343
344 c0 = Element->c[0];
345 c1 = Element->c[1];
346 c2 = Element->c[2];
347
348 fcoeffs[0] = c0 * t2 * t2;
349 fcoeffs[1] = 4.0 * c0 * t1 * t2;
350 fcoeffs[2] = 2.0 * c0 * (2.0 * t1 * t1 + t0 * t2) + c1 * t2;
351 fcoeffs[3] = 2.0 * t1 * (2.0 * c0 * t0 + c1);
352 fcoeffs[4] = t0 * (c0 * t0 + c1) + c2;
353
354 break;
355
356 case BLOB_ELLIPSOID:
357
358 MInvTransPoint(PP, P, Element->Trans);
359 MInvTransDirection(DD, D, Element->Trans);
360
361 V1 = PP - Element->O;
362
363 t0 = V1.lengthSqr();
364 t1 = dot(V1, DD);
365 t2 = DD.lengthSqr();
366
367 c0 = Element->c[0];
368 c1 = Element->c[1];
369 c2 = Element->c[2];
370
371 fcoeffs[0] = c0 * t2 * t2;
372 fcoeffs[1] = 4.0 * c0 * t1 * t2;
373 fcoeffs[2] = 2.0 * c0 * (2.0 * t1 * t1 + t0 * t2) + c1 * t2;
374 fcoeffs[3] = 2.0 * t1 * (2.0 * c0 * t0 + c1);
375 fcoeffs[4] = t0 * (c0 * t0 + c1) + c2;
376
377 break;
378
379 case BLOB_BASE_HEMISPHERE:
380 case BLOB_APEX_HEMISPHERE:
381
382 MInvTransPoint(PP, P, Element->Trans);
383 MInvTransDirection(DD, D, Element->Trans);
384
385 if (Element->Type == BLOB_APEX_HEMISPHERE)
386 {
387 PP[Z] -= Element->len;
388 }
389
390 t0 = PP.lengthSqr();
391 t1 = dot(PP, DD);
392 t2 = DD.lengthSqr();
393
394 c0 = Element->c[0];
395 c1 = Element->c[1];
396 c2 = Element->c[2];
397
398 fcoeffs[0] = c0 * t2 * t2;
399 fcoeffs[1] = 4.0 * c0 * t1 * t2;
400 fcoeffs[2] = 2.0 * c0 * (2.0 * t1 * t1 + t0 * t2) + c1 * t2;
401 fcoeffs[3] = 2.0 * t1 * (2.0 * c0 * t0 + c1);
402 fcoeffs[4] = t0 * (c0 * t0 + c1) + c2;
403
404 break;
405
406 case BLOB_CYLINDER:
407
408 /* Transform ray into cylinder space. */
409
410 MInvTransPoint(PP, P, Element->Trans);
411 MInvTransDirection(DD, D, Element->Trans);
412
413 t0 = PP[X] * PP[X] + PP[Y] * PP[Y];
414 t1 = PP[X] * DD[X] + PP[Y] * DD[Y];
415 t2 = DD[X] * DD[X] + DD[Y] * DD[Y];
416
417 c0 = Element->c[0];
418 c1 = Element->c[1];
419 c2 = Element->c[2];
420
421 fcoeffs[0] = c0 * t2 * t2;
422 fcoeffs[1] = 4.0 * c0 * t1 * t2;
423 fcoeffs[2] = 2.0 * c0 * (2.0 * t1 * t1 + t0 * t2) + c1 * t2;
424 fcoeffs[3] = 2.0 * t1 * (2.0 * c0 * t0 + c1);
425 fcoeffs[4] = t0 * (c0 * t0 + c1) + c2;
426
427 break;
428
429 default:
430
431 throw POV_EXCEPTION_STRING("Unknown blob component in All_Blob_Intersections().");
432 }
433
434 for (j = 0; j < 5; j++)
435 {
436 coeffs[j] += fcoeffs[j];
437 }
438 }
439 else
440 {
441 /*
442 * We are losing the influence of a component -->
443 * subtract off its coefficients.
444 */
445 fcoeffs = &Thread->Blob_Coefficients[intervals[i].Element->index * 5];
446
447 for (j = 0; j < 5; j++)
448 {
449 coeffs[j] -= fcoeffs[j];
450 }
451
452 /* If no components are currently affecting the ray ---> skip ahead. */
453
454 if (--in_flag == 0)
455 {
456 continue;
457 }
458 }
459
460 /*
461 * If the following intersection lies close to the current intersection
462 * then first add/subtract next region before testing. [DB 7/94]
463 */
464
465 if ((i + 1 < cnt) && (fabs(intervals[i].bound - intervals[i + 1].bound) < EPSILON))
466 {
467 continue;
468 }
469
470 /*
471 * Transform polynomial in a way that the interval boundaries are moved
472 * to 0 and 1, i. e. the roots of interest are between 0 and 1. [DB 10/94]
473 */
474
475 l = intervals[i].bound;
476 w = intervals[i+1].bound - l;
477
478 newcoeffs[0] = coeffs[0] * w * w * w * w;
479 newcoeffs[1] = (coeffs[1] + 4.0 * coeffs[0] * l) * w * w * w;
480 newcoeffs[2] = (3.0 * l * (2.0 * coeffs[0] * l + coeffs[1]) + coeffs[2]) * w * w;
481 newcoeffs[3] = (2.0 * l * (2.0 * l * (coeffs[0] * l + 0.75 * coeffs[1]) + coeffs[2]) + coeffs[3]) * w;
482 newcoeffs[4] = l * (l * (l * (coeffs[0] * l + coeffs[1]) + coeffs[2]) + coeffs[3]) + coeffs[4];
483
484 /* Calculate coefficients of corresponding bezier curve. [DB 10/94] */
485
486 dk[0] = newcoeffs[4];
487 dk[1] = newcoeffs[4] + 0.25 * newcoeffs[3];
488 dk[2] = newcoeffs[4] + 0.50 * (newcoeffs[3] + newcoeffs[2] / 3.0);
489 dk[3] = newcoeffs[4] + 0.50 * (1.5 * newcoeffs[3] + newcoeffs[2] + 0.5 * newcoeffs[1]);
490 dk[4] = newcoeffs[4] + newcoeffs[3] + newcoeffs[2] + newcoeffs[1] + newcoeffs[0];
491
492 /*
493 * Skip this interval if the ray doesn't intersect the convex hull of the
494 * bezier curve, because no valid intersection will be found. [DB 10/94]
495 */
496
497 if (((dk[0] >= 0.0) && (dk[1] >= 0.0) && (dk[2] >= 0.0) && (dk[3] >= 0.0) && (dk[4] >= 0.0)) ||
498 ((dk[0] <= 0.0) && (dk[1] <= 0.0) && (dk[2] <= 0.0) && (dk[3] <= 0.0) && (dk[4] <= 0.0)))
499 {
500 continue;
501 }
502
503 /*
504 * Now we could do bezier clipping to find the roots
505 * but I have no idea how this works. [DB 2/95]
506 */
507
508
509 /* Solve polynomial. */
510
511 root_count = Solve_Polynomial(4, coeffs, roots, Test_Flag(this, STURM_FLAG), 1.0e-11, Thread->Stats());
512
513 /* See if any of the roots are valid. */
514
515 for (j = 0; j < root_count; j++)
516 {
517 dist = roots[j];
518
519 /*
520 * First see if the root is in the interval of influence of
521 * the currently active components.
522 */
523
524 if ((dist >= intervals[i].bound) &&
525 (dist <= intervals[i+1].bound))
526 {
527 /* Correct distance. */
528
529 dist = (dist * max_bound + start_dist) / len;
530
531 if ((dist > depthTolerance) && (dist < MAX_DISTANCE))
532 {
533 IPoint = ray.Evaluate(dist);
534
535 if (Clip.empty() || Point_In_Clip(IPoint, Clip, Thread))
536 {
537 Depth_Stack->push(Intersection(dist, IPoint, this));
538
539 Intersection_Found = true;
540 }
541 }
542 }
543 }
544
545 /*
546 * If the blob isn't used inside a CSG and we have found at least
547 * one intersection then we are ready, because all possible intersections
548 * will be further away (we have a sorted list!). [DB 7/94]
549 */
550
551 if (!(Type & IS_CHILD_OBJECT) && (Intersection_Found))
552 {
553 break;
554 }
555 }
556
557 if (Intersection_Found)
558 Thread->Stats()[Ray_Blob_Tests_Succeeded]++;
559
560 return (Intersection_Found);
561 }
562
563
564
565 /*****************************************************************************
566 *
567 * FUNCTION
568 *
569 * insert_hit
570 *
571 * INPUT
572 *
573 * Blob - Pointer to blob structure
574 * Element - Element to insert
575 * t0, t1 - Intersection depths
576 *
577 * OUTPUT
578 *
579 * intervals - Pointer to sorted list of hits
580 * cnt - Number of hits in intervals
581 *
582 * RETURNS
583 *
584 * AUTHOR
585 *
586 * Alexander Enzmann
587 *
588 * DESCRIPTION
589 *
590 * Store the points of intersection. Keep track of: whether this is
591 * the start or end point of the hit, which component was pierced
592 * by the ray, and the point along the ray that the hit occured at.
593 *
594 * CHANGES
595 *
596 * Oct 1994 : Modified to use memmove instead of loops for copying. [DB]
597 * Sep 1995 : Changed to allow use of memcpy if memmove isn't available. [AED]
598 * Jul 1996 : Changed to use POV_MEMMOVE, which can be memmove or pov_memmove.
599 * Oct 1996 : Changed to avoid unnecessary compares. [DB]
600 *
601 ******************************************************************************/
602
insert_hit(const Blob_Element * Element,DBL t0,DBL t1,Blob_Interval_Struct * intervals,unsigned int * cnt)603 void Blob::insert_hit(const Blob_Element *Element, DBL t0, DBL t1, Blob_Interval_Struct *intervals, unsigned int *cnt)
604 {
605 unsigned int k;
606
607 /* We are entering the component. */
608
609 intervals[*cnt].type = Element->Type | ENTERING;
610 intervals[*cnt].bound = t0;
611 intervals[*cnt].Element = Element;
612
613 for (k = 0; t0 > intervals[k].bound; k++);
614
615 if (k < *cnt)
616 {
617 /*
618 * This hit point is smaller than one that already exists -->
619 * bump the rest and insert it here.
620 */
621
622 POV_MEMMOVE(&intervals[k+1], &intervals[k], (*cnt-k)*sizeof(Blob_Interval_Struct));
623
624 /* We are entering the component. */
625
626 intervals[k].type = Element->Type | ENTERING;
627 intervals[k].bound = t0;
628 intervals[k].Element = Element;
629
630 (*cnt)++;
631
632 /* We are exiting the component. */
633
634 intervals[*cnt].type = Element->Type | EXITING;
635 intervals[*cnt].bound = t1;
636 intervals[*cnt].Element = Element;
637
638 for (k = k + 1; t1 > intervals[k].bound; k++);
639
640 if (k < *cnt)
641 {
642 POV_MEMMOVE(&intervals[k+1], &intervals[k], (*cnt-k)*sizeof(Blob_Interval_Struct));
643
644 /* We are exiting the component. */
645
646 intervals[k].type = Element->Type | EXITING;
647 intervals[k].bound = t1;
648 intervals[k].Element = Element;
649 }
650
651 (*cnt)++;
652 }
653 else
654 {
655 /* Just plop the start and end points at the end of the list */
656
657 (*cnt)++;
658
659 /* We are exiting the component. */
660
661 intervals[*cnt].type = Element->Type | EXITING;
662 intervals[*cnt].bound = t1;
663 intervals[*cnt].Element = Element;
664
665 (*cnt)++;
666 }
667 }
668
669
670
671 /*****************************************************************************
672 *
673 * FUNCTION
674 *
675 * intersect_cylinder
676 *
677 * INPUT
678 *
679 * Element - Pointer to element structure
680 * P, D - Ray = P + t * D
681 * mindist - Min. valid distance
682 *
683 * OUTPUT
684 *
685 * tmin, tmax - Intersection depths found
686 *
687 * RETURNS
688 *
689 * AUTHOR
690 *
691 * Dieter Bayer (with help from Alexander Enzmann)
692 *
693 * DESCRIPTION
694 *
695 * -
696 *
697 * CHANGES
698 *
699 * Jul 1994 : Creation.
700 *
701 ******************************************************************************/
702
intersect_cylinder(const Blob_Element * Element,const Vector3d & P,const Vector3d & D,DBL mindist,DBL * tmin,DBL * tmax)703 int Blob::intersect_cylinder(const Blob_Element *Element, const Vector3d& P, const Vector3d& D, DBL mindist, DBL *tmin, DBL *tmax)
704 {
705 DBL a, b, c, d, t, u, v, w, len;
706 Vector3d PP, DD;
707
708 /* Transform ray into cylinder space. */
709
710 MInvTransPoint(PP, P, Element->Trans);
711 MInvTransDirection(DD, D, Element->Trans);
712
713 len = DD.length();
714 DD /= len;
715
716 /* Intersect ray with cylinder. */
717
718 a = DD[X] * DD[X] + DD[Y] * DD[Y];
719
720 if (a > EPSILON)
721 {
722 b = PP[X] * DD[X] + PP[Y] * DD[Y];
723 c = PP[X] * PP[X] + PP[Y] * PP[Y] - Element->rad2;
724
725 d = b * b - a * c;
726
727 if (d > EPSILON)
728 {
729 d = sqrt(d);
730
731 t = ( - b + d) / a;
732
733 w = PP[Z] + t * DD[Z];
734
735 if ((w >= 0.0) && (w <= Element->len))
736 {
737 if (t < *tmin) { *tmin = t; }
738 if (t > *tmax) { *tmax = t; }
739 }
740
741 t = ( - b - d) / a;
742
743 w = PP[Z] + t * DD[Z];
744
745 if ((w >= 0.0) && (w <= Element->len))
746 {
747 if (t < *tmin) { *tmin = t; }
748 if (t > *tmax) { *tmax = t; }
749 }
750 }
751 }
752
753 /* Intersect base/cap plane. */
754
755 if (fabs(DD[Z]) > EPSILON)
756 {
757 /* Intersect base plane. */
758
759 t = - PP[Z] / DD[Z];
760
761 u = PP[X] + t * DD[X];
762 v = PP[Y] + t * DD[Y];
763
764 if ((u * u + v * v) <= Element->rad2)
765 {
766 if (t < *tmin) { *tmin = t; }
767 if (t > *tmax) { *tmax = t; }
768 }
769
770 /* Intersect cap plane. */
771
772 t = (Element->len - PP[Z]) / DD[Z];
773
774 u = PP[X] + t * DD[X];
775 v = PP[Y] + t * DD[Y];
776
777 if ((u * u + v * v) <= Element->rad2)
778 {
779 if (t < *tmin) { *tmin = t; }
780 if (t > *tmax) { *tmax = t; }
781 }
782 }
783
784 /* Check if the intersections are valid. */
785
786 *tmin /= len;
787 *tmax /= len;
788
789 if (*tmin < mindist) { *tmin = 0.0; }
790 if (*tmax < mindist) { *tmax = 0.0; }
791
792 if (*tmin >= *tmax)
793 {
794 return (false);
795 }
796
797 return (true);
798 }
799
800
801
802 /*****************************************************************************
803 *
804 * FUNCTION
805 *
806 * intersect_ellipsoid
807 *
808 * INPUT
809 *
810 * Element - Pointer to element structure
811 * P, D - Ray = P + t * D
812 * mindist - Min. valid distance
813 *
814 * OUTPUT
815 *
816 * tmin, tmax - Intersection depths found
817 *
818 * RETURNS
819 *
820 * AUTHOR
821 *
822 * Dieter Bayer
823 *
824 * DESCRIPTION
825 *
826 * -
827 *
828 * CHANGES
829 *
830 * Sep 1994 : Creation.
831 *
832 ******************************************************************************/
833
intersect_ellipsoid(const Blob_Element * Element,const Vector3d & P,const Vector3d & D,DBL mindist,DBL * tmin,DBL * tmax)834 int Blob::intersect_ellipsoid(const Blob_Element *Element, const Vector3d& P, const Vector3d& D, DBL mindist, DBL *tmin, DBL *tmax)
835 {
836 DBL b, d, t, len;
837 Vector3d V1, PP, DD;
838
839 MInvTransPoint(PP, P, Element->Trans);
840 MInvTransDirection(DD, D, Element->Trans);
841
842 len = DD.length();
843 DD /= len;
844
845 V1 = PP - Element->O;
846 b = dot(V1, DD);
847 t = V1.lengthSqr();
848
849 d = b * b - t + Element->rad2;
850
851 if (d < EPSILON)
852 {
853 return (false);
854 }
855
856 d = sqrt(d);
857
858 *tmax = ( - b + d) / len; if (*tmax < mindist) { *tmax = 0.0; }
859 *tmin = ( - b - d) / len; if (*tmin < mindist) { *tmin = 0.0; }
860
861 if (*tmax == *tmin)
862 {
863 return (false);
864 }
865 else
866 {
867 if (*tmax < *tmin)
868 {
869 d = *tmin; *tmin = *tmax; *tmax = d;
870 }
871 }
872
873 return (true);
874 }
875
876
877
878 /*****************************************************************************
879 *
880 * FUNCTION
881 *
882 * intersect_hemisphere
883 *
884 * INPUT
885 *
886 * Element - Pointer to element structure
887 * P, D - Ray = P + t * D
888 * mindist - Min. valid distance
889 *
890 * OUTPUT
891 *
892 * tmin, tmax - Intersection depths found
893 *
894 * RETURNS
895 *
896 * AUTHOR
897 *
898 * Dieter Bayer
899 *
900 * DESCRIPTION
901 *
902 * -
903 *
904 * CHANGES
905 *
906 * Jul 1994 : Creation (with help from Alexander Enzmann).
907 *
908 ******************************************************************************/
909
intersect_hemisphere(const Blob_Element * Element,const Vector3d & P,const Vector3d & D,DBL mindist,DBL * tmin,DBL * tmax)910 int Blob::intersect_hemisphere(const Blob_Element *Element, const Vector3d& P, const Vector3d& D, DBL mindist, DBL *tmin, DBL *tmax)
911 {
912 DBL b, d, t, z1, z2, len;
913 Vector3d PP, DD;
914
915 /* Transform ray into hemisphere space. */
916
917 MInvTransPoint(PP, P, Element->Trans);
918 MInvTransDirection(DD, D, Element->Trans);
919
920 len = DD.length();
921 DD /= len;
922
923 if (Element->Type == BLOB_BASE_HEMISPHERE)
924 {
925 b = dot(PP, DD);
926 t = PP.lengthSqr();
927
928 d = b * b - t + Element->rad2;
929
930 if (d < EPSILON)
931 {
932 return (false);
933 }
934
935 d = sqrt(d);
936
937 *tmax = - b + d;
938 *tmin = - b - d;
939
940 if (*tmax < *tmin)
941 {
942 d = *tmin; *tmin = *tmax; *tmax = d;
943 }
944
945 /* Cut intersection at the plane. */
946
947 z1 = PP[Z] + *tmin * DD[Z];
948 z2 = PP[Z] + *tmax * DD[Z];
949
950 /* If both points are inside --> no intersection */
951
952 if ((z1 >= 0.0) && (z2 >= 0.0))
953 {
954 return (false);
955 }
956
957 /* If both points are outside --> intersections found */
958
959 if ((z1 < 0.0) && (z2 < 0.0))
960 {
961 *tmin /= len;
962 *tmax /= len;
963
964 return (true);
965 }
966
967 /* Determine intersection with plane. */
968
969 t = - PP[Z] / DD[Z];
970
971 if (z1 >= 0.0)
972 {
973 /* Ray is crossing the plane from inside to outside. */
974
975 *tmin = (t < mindist) ? 0.0 : t;
976 }
977 else
978 {
979 /* Ray is crossing the plane from outside to inside. */
980
981 *tmax = (t < mindist) ? 0.0 : t;
982 }
983
984 *tmin /= len;
985 *tmax /= len;
986
987 return (true);
988 }
989 else
990 {
991 PP[Z] -= Element->len;
992
993 b = dot(PP, DD);
994 t = PP.lengthSqr();
995
996 d = b * b - t + Element->rad2;
997
998 if (d < EPSILON)
999 {
1000 return (false);
1001 }
1002
1003 d = sqrt(d);
1004
1005 *tmax = - b + d;
1006 *tmin = - b - d;
1007
1008 if (*tmax < *tmin)
1009 {
1010 d = *tmin; *tmin = *tmax; *tmax = d;
1011 }
1012
1013 /* Cut intersection at the plane. */
1014
1015 z1 = PP[Z] + *tmin * DD[Z];
1016 z2 = PP[Z] + *tmax * DD[Z];
1017
1018 /* If both points are inside --> no intersection */
1019
1020 if ((z1 <= 0.0) && (z2 <= 0.0))
1021 {
1022 return (false);
1023 }
1024
1025 /* If both points are outside --> intersections found */
1026
1027 if ((z1 > 0.0) && (z2 > 0.0))
1028 {
1029 *tmin /= len;
1030 *tmax /= len;
1031
1032 return (true);
1033 }
1034
1035 /* Determine intersection with plane. */
1036
1037 t = - PP[Z] / DD[Z];
1038
1039 if (z1 <= 0.0)
1040 {
1041 /* Ray is crossing the plane from inside to outside. */
1042
1043 *tmin = (t < mindist) ? 0.0 : t;
1044 }
1045 else
1046 {
1047 /* Ray is crossing the plane from outside to inside. */
1048
1049 *tmax = (t < mindist) ? 0.0 : t;
1050 }
1051
1052 *tmin /= len;
1053 *tmax /= len;
1054
1055 return (true);
1056 }
1057 }
1058
1059
1060
1061 /*****************************************************************************
1062 *
1063 * FUNCTION
1064 *
1065 * intersect_sphere
1066 *
1067 * INPUT
1068 *
1069 * Element - Pointer to element structure
1070 * P, D - Ray = P + t * D
1071 * mindist - Min. valid distance
1072 *
1073 * OUTPUT
1074 *
1075 * tmin, tmax - Intersection depths found
1076 *
1077 * RETURNS
1078 *
1079 * AUTHOR
1080 *
1081 * Dieter Bayer
1082 *
1083 * DESCRIPTION
1084 *
1085 * -
1086 *
1087 * CHANGES
1088 *
1089 * Jul 1994 : Creation (with help from Alexander Enzmann).
1090 *
1091 ******************************************************************************/
1092
intersect_sphere(const Blob_Element * Element,const Vector3d & P,const Vector3d & D,DBL mindist,DBL * tmin,DBL * tmax)1093 int Blob::intersect_sphere(const Blob_Element *Element, const Vector3d& P, const Vector3d& D, DBL mindist, DBL *tmin, DBL *tmax)
1094 {
1095 DBL b, d, t;
1096 Vector3d V1;
1097
1098 V1 = P - Element->O;
1099 b = dot(V1, D);
1100 t = V1.lengthSqr();
1101
1102 d = b * b - t + Element->rad2;
1103
1104 if (d < EPSILON)
1105 {
1106 return (false);
1107 }
1108
1109 d = sqrt(d);
1110
1111 *tmax = - b + d; if (*tmax < mindist) { *tmax = 0.0; }
1112 *tmin = - b - d; if (*tmin < mindist) { *tmin = 0.0; }
1113
1114 if (*tmax == *tmin)
1115 {
1116 return (false);
1117 }
1118 else
1119 {
1120 if (*tmax < *tmin)
1121 {
1122 d = *tmin; *tmin = *tmax; *tmax = d;
1123 }
1124 }
1125
1126 return (true);
1127 }
1128
1129
1130
1131 /*****************************************************************************
1132 *
1133 * FUNCTION
1134 *
1135 * intersect_element
1136 *
1137 * INPUT
1138 *
1139 * P, D - Ray = P + t * D
1140 * Element - Pointer to element structure
1141 * mindist - Min. valid distance
1142 *
1143 * OUTPUT
1144 *
1145 * tmin, tmax - Intersection depths found
1146 *
1147 * RETURNS
1148 *
1149 * AUTHOR
1150 *
1151 * Dieter Bayer
1152 *
1153 * DESCRIPTION
1154 *
1155 * -
1156 *
1157 * CHANGES
1158 *
1159 * Jul 1994 : Creation.
1160 *
1161 ******************************************************************************/
1162
intersect_element(const Vector3d & P,const Vector3d & D,const Blob_Element * Element,DBL mindist,DBL * tmin,DBL * tmax,TraceThreadData * Thread)1163 int Blob::intersect_element(const Vector3d& P, const Vector3d& D, const Blob_Element *Element, DBL mindist, DBL *tmin, DBL *tmax, TraceThreadData *Thread)
1164 {
1165 #ifdef BLOB_EXTRA_STATS
1166 Thread->Stats()[Blob_Element_Tests]++;
1167 #endif
1168
1169 *tmin = BOUND_HUGE;
1170 *tmax = - BOUND_HUGE;
1171
1172 switch (Element->Type)
1173 {
1174 case BLOB_SPHERE:
1175
1176 if (!intersect_sphere(Element, P, D, mindist, tmin, tmax))
1177 {
1178 return (false);
1179 }
1180
1181 break;
1182
1183 case BLOB_ELLIPSOID:
1184
1185 if (!intersect_ellipsoid(Element, P, D, mindist, tmin, tmax))
1186 {
1187 return (false);
1188 }
1189
1190 break;
1191
1192 case BLOB_BASE_HEMISPHERE:
1193 case BLOB_APEX_HEMISPHERE:
1194
1195 if (!intersect_hemisphere(Element, P, D, mindist, tmin, tmax))
1196 {
1197 return (false);
1198 }
1199
1200 break;
1201
1202 case BLOB_CYLINDER:
1203
1204 if (!intersect_cylinder(Element, P, D, mindist, tmin, tmax))
1205 {
1206 return (false);
1207 }
1208
1209 break;
1210 }
1211
1212 #ifdef BLOB_EXTRA_STATS
1213 Thread->Stats()[Blob_Element_Tests_Succeeded]++;
1214 #endif
1215
1216 return (true);
1217 }
1218
1219
1220
1221 /*****************************************************************************
1222 *
1223 * FUNCTION
1224 *
1225 * determine_influences
1226 *
1227 * INPUT
1228 *
1229 * P, D - Ray = P + t * D
1230 * Blob - Pointer to blob structure
1231 * mindist - Min. valid distance
1232 *
1233 * OUTPUT
1234 *
1235 * intervals - Sorted list of intersections found
1236 *
1237 * RETURNS
1238 *
1239 * int - Number of intersection found
1240 *
1241 * AUTHOR
1242 *
1243 * Dieter Bayer
1244 *
1245 * DESCRIPTION
1246 *
1247 * Make a sorted list of points along the ray at which the various blob
1248 * components start and stop adding their influence.
1249 *
1250 * CHANGES
1251 *
1252 * Jul 1994 : Added code for bounding hierarchy traversal. [DB]
1253 *
1254 ******************************************************************************/
1255
determine_influences(const Vector3d & P,const Vector3d & D,DBL mindist,Blob_Interval_Struct * intervals,TraceThreadData * Thread) const1256 int Blob::determine_influences(const Vector3d& P, const Vector3d& D, DBL mindist, Blob_Interval_Struct *intervals, TraceThreadData *Thread) const
1257 {
1258 unsigned int cnt, size;
1259 DBL b, t, t0, t1;
1260 Vector3d V1;
1261 BSPHERE_TREE *Tree;
1262 BSPHERE_TREE **Queue = reinterpret_cast<BSPHERE_TREE **>(Thread->Blob_Queue);
1263
1264 cnt = 0;
1265
1266 if (Data->Tree == nullptr)
1267 {
1268 /* There's no bounding hierarchy so just step through all elements. */
1269
1270 for (vector<Blob_Element>::iterator i = Data->Entry.begin(); i != Data->Entry.end(); ++i)
1271 {
1272 if (intersect_element(P, D, &(*i), mindist, &t0, &t1, Thread))
1273 {
1274 insert_hit(&(*i), t0, t1, intervals, &cnt);
1275 }
1276 }
1277 }
1278 else
1279 {
1280 /* Use blob's bounding hierarchy. */
1281
1282 size = 0;
1283
1284 Queue[size++] = Data->Tree;
1285
1286 while (size > 0)
1287 {
1288 Tree = Queue[--size];
1289
1290 /* Test if current node is a leaf. */
1291
1292 if (Tree->Entries <= 0)
1293 {
1294 /* Test element. */
1295
1296 if (intersect_element(P, D, reinterpret_cast<Blob_Element *>(Tree->Node), mindist, &t0, &t1, Thread))
1297 {
1298 insert_hit(reinterpret_cast<Blob_Element *>(Tree->Node), t0, t1, intervals, &cnt);
1299 }
1300 }
1301 else
1302 {
1303 /* Test all sub-nodes. */
1304
1305 for (int i = 0; i < (int)Tree->Entries; i++)
1306 {
1307 #ifdef BLOB_EXTRA_STATS
1308 Thread->Stats()[Blob_Bound_Tests]++;
1309 #endif
1310
1311 V1 = Tree->Node[i]->C - P;
1312 b = dot(V1, D);
1313 t = V1.lengthSqr();
1314
1315 if ((t - Sqr(b)) <= Tree->Node[i]->r2)
1316 {
1317 #ifdef BLOB_EXTRA_STATS
1318 Thread->Stats()[Blob_Bound_Tests_Succeeded]++;
1319 #endif
1320
1321 if (insert_node(Tree->Node[i], &size, Thread))
1322 Queue = reinterpret_cast<BSPHERE_TREE **>(Thread->Blob_Queue);
1323 }
1324 }
1325 }
1326 }
1327 }
1328
1329 return (cnt);
1330 }
1331
1332
1333
1334 /*****************************************************************************
1335 *
1336 * FUNCTION
1337 *
1338 * calculate_element_field
1339 *
1340 * INPUT
1341 *
1342 * Element - Pointer to element structure
1343 * P - Point whos field value is calculated
1344 *
1345 * OUTPUT
1346 *
1347 * RETURNS
1348 *
1349 * DBL - Field value
1350 *
1351 * AUTHOR
1352 *
1353 * Alexander Enzmann
1354 *
1355 * DESCRIPTION
1356 *
1357 * Calculate the field value of a single element in a given point P
1358 * (which must already have been transformed into blob space).
1359 *
1360 * CHANGES
1361 *
1362 * Jul 1994 : Added code for cylindrical and ellipsoidical blobs. [DB]
1363 *
1364 ******************************************************************************/
1365
calculate_element_field(const Blob_Element * Element,const Vector3d & P)1366 DBL Blob::calculate_element_field(const Blob_Element *Element, const Vector3d& P)
1367 {
1368 DBL rad2, density;
1369 Vector3d V1, PP;
1370
1371 density = 0.0;
1372
1373 switch (Element->Type)
1374 {
1375 case BLOB_SPHERE:
1376
1377 V1 = P - Element->O;
1378
1379 rad2 = V1.lengthSqr();
1380
1381 if (rad2 < Element->rad2)
1382 {
1383 density = rad2 * (rad2 * Element->c[0] + Element->c[1]) + Element->c[2];
1384 }
1385
1386 break;
1387
1388 case BLOB_ELLIPSOID:
1389
1390 MInvTransPoint(PP, P, Element->Trans);
1391
1392 V1 = PP - Element->O;
1393
1394 rad2 = V1.lengthSqr();
1395
1396 if (rad2 < Element->rad2)
1397 {
1398 density = rad2 * (rad2 * Element->c[0] + Element->c[1]) + Element->c[2];
1399 }
1400
1401 break;
1402
1403 case BLOB_BASE_HEMISPHERE:
1404
1405 MInvTransPoint(PP, P, Element->Trans);
1406
1407 if (PP[Z] <= 0.0)
1408 {
1409 rad2 = PP.lengthSqr();
1410
1411 if (rad2 <= Element->rad2)
1412 {
1413 density = rad2 * (rad2 * Element->c[0] + Element->c[1]) + Element->c[2];
1414 }
1415 }
1416
1417 break;
1418
1419 case BLOB_APEX_HEMISPHERE:
1420
1421 MInvTransPoint(PP, P, Element->Trans);
1422
1423 PP[Z] -= Element->len;
1424
1425 if (PP[Z] >= 0.0)
1426 {
1427 rad2 = PP.lengthSqr();
1428
1429 if (rad2 <= Element->rad2)
1430 {
1431 density = rad2 * (rad2 * Element->c[0] + Element->c[1]) + Element->c[2];
1432 }
1433 }
1434
1435 break;
1436
1437 case BLOB_CYLINDER:
1438
1439 MInvTransPoint(PP, P, Element->Trans);
1440
1441 if ((PP[Z] >= 0.0) && (PP[Z] <= Element->len))
1442 {
1443 if ((rad2 = Sqr(PP[X]) + Sqr(PP[Y])) <= Element->rad2)
1444 {
1445 density = rad2 * (rad2 * Element->c[0] + Element->c[1]) + Element->c[2];
1446 }
1447 }
1448
1449 break;
1450 }
1451
1452 return (density);
1453 }
1454
1455
1456
1457 /*****************************************************************************
1458 *
1459 * FUNCTION
1460 *
1461 * calculate_field_value
1462 *
1463 * INPUT
1464 *
1465 * Blob - Pointer to blob structure
1466 * P - Point whos field value is calculated
1467 *
1468 * OUTPUT
1469 *
1470 * RETURNS
1471 *
1472 * DBL - Field value
1473 *
1474 * AUTHOR
1475 *
1476 * Dieter Bayer
1477 *
1478 * DESCRIPTION
1479 *
1480 * Calculate the field value of a blob in a given point P
1481 * (which must already have been transformed into blob space).
1482 *
1483 * CHANGES
1484 *
1485 * Jul 1994 : Added code for bounding hierarchy traversal. [DB]
1486 *
1487 ******************************************************************************/
1488
calculate_field_value(const Vector3d & P,TraceThreadData * Thread) const1489 DBL Blob::calculate_field_value(const Vector3d& P, TraceThreadData *Thread) const
1490 {
1491 unsigned int size;
1492 DBL density, rad2;
1493 Vector3d V1;
1494 BSPHERE_TREE *Tree;
1495 BSPHERE_TREE **Queue = reinterpret_cast<BSPHERE_TREE **>(Thread->Blob_Queue);
1496
1497 density = 0.0;
1498
1499 if (Data->Tree == nullptr)
1500 {
1501 /* There's no tree --> step through all elements. */
1502
1503 for (vector<Blob_Element>::iterator i = Data->Entry.begin(); i != Data->Entry.end(); ++i)
1504 {
1505 density += calculate_element_field(&(*i), P);
1506 }
1507 }
1508 else
1509 {
1510 /* A tree exists --> step through the tree. */
1511
1512 size = 0;
1513
1514 Queue[size++] = Data->Tree;
1515
1516 while (size > 0)
1517 {
1518 Tree = Queue[--size];
1519
1520 /* Test if current node is a leaf. */
1521
1522 if (Tree->Entries <= 0)
1523 {
1524 density += calculate_element_field(reinterpret_cast<Blob_Element *>(Tree->Node), P);
1525 }
1526 else
1527 {
1528 /* Test all sub-nodes. */
1529
1530 for (int i = 0; i < (int)Tree->Entries; i++)
1531 {
1532 /* Insert sub-node if we are inside. */
1533
1534 V1 = P - Tree->Node[i]->C;
1535
1536 rad2 = V1.lengthSqr();
1537
1538 if (rad2 <= Tree->Node[i]->r2)
1539 if (insert_node(Tree->Node[i], &size, Thread))
1540 Queue = reinterpret_cast<BSPHERE_TREE **>(Thread->Blob_Queue);
1541 }
1542 }
1543 }
1544 }
1545
1546 return (density);
1547 }
1548
1549
1550
1551 /*****************************************************************************
1552 *
1553 * FUNCTION
1554 *
1555 * Inside_Blob
1556 *
1557 * INPUT
1558 *
1559 * Test_Point - Point to test
1560 * Object - Pointer to blob structure
1561 *
1562 * OUTPUT
1563 *
1564 * RETURNS
1565 *
1566 * int - true if Test_Point is inside
1567 *
1568 * AUTHOR
1569 *
1570 * Alexander Enzmann
1571 *
1572 * DESCRIPTION
1573 *
1574 * Calculate the density at the given point and then compare to
1575 * the threshold to see if we are in or out of the blob.
1576 *
1577 * CHANGES
1578 *
1579 * -
1580 *
1581 ******************************************************************************/
1582
Inside(const Vector3d & Test_Point,TraceThreadData * Thread) const1583 bool Blob::Inside(const Vector3d& Test_Point, TraceThreadData *Thread) const
1584 {
1585 Vector3d New_Point;
1586
1587 /* Transform the point into blob space. */
1588
1589 if (Trans != nullptr)
1590 {
1591 MInvTransPoint(New_Point, Test_Point, Trans);
1592 }
1593 else
1594 {
1595 New_Point = Test_Point;
1596 }
1597
1598 if (calculate_field_value(New_Point, Thread) > Data->Threshold - INSIDE_TOLERANCE)
1599 {
1600 /* We are inside. */
1601
1602 return (!Test_Flag(this, INVERTED_FLAG));
1603 }
1604 else
1605 {
1606 /* We are outside. */
1607
1608 return (Test_Flag(this, INVERTED_FLAG));
1609 }
1610 }
1611
1612
GetPotential(const Vector3d & globalPoint,bool subtractThreshold,TraceThreadData * pThread) const1613 double Blob::GetPotential (const Vector3d& globalPoint, bool subtractThreshold, TraceThreadData *pThread) const
1614 {
1615 Vector3d localPoint;
1616
1617 if (Trans != nullptr)
1618 MInvTransPoint (localPoint, globalPoint, Trans);
1619 else
1620 localPoint = globalPoint;
1621
1622 double potential = calculate_field_value (localPoint, pThread);
1623 if (subtractThreshold)
1624 potential -= Data->Threshold;
1625
1626 if (Test_Flag (this, INVERTED_FLAG))
1627 return -potential;
1628 else
1629 return potential;
1630 }
1631
1632
1633 /*****************************************************************************
1634 *
1635 * FUNCTION
1636 *
1637 * element_normal
1638 *
1639 * INPUT
1640 *
1641 * P - Surface point
1642 * Element - Pointer to element structure
1643 *
1644 * OUTPUT
1645 *
1646 * Result - Element's normal
1647 *
1648 * RETURNS
1649 *
1650 * AUTHOR
1651 *
1652 * Dieter Bayer
1653 *
1654 * DESCRIPTION
1655 *
1656 * Calculate the normal of a single element in the point P.
1657 *
1658 * CHANGES
1659 *
1660 * Jul 1994 : Creation (with help from Alexander Enzmann).
1661 *
1662 ******************************************************************************/
1663
element_normal(Vector3d & Result,const Vector3d & P,const Blob_Element * Element)1664 void Blob::element_normal(Vector3d& Result, const Vector3d& P, const Blob_Element *Element)
1665 {
1666 DBL val, dist;
1667 Vector3d V1, PP;
1668
1669 switch (Element->Type)
1670 {
1671 case BLOB_SPHERE:
1672
1673 V1 = P - Element->O;
1674
1675 dist = V1.lengthSqr();
1676
1677 if (dist <= Element->rad2)
1678 {
1679 val = -2.0 * Element->c[0] * dist - Element->c[1];
1680
1681 Result += val * V1;
1682 }
1683
1684 break;
1685
1686 case BLOB_ELLIPSOID:
1687
1688 MInvTransPoint(PP, P, Element->Trans);
1689
1690 V1 = PP - Element->O;
1691
1692 dist = V1.lengthSqr();
1693
1694 if (dist <= Element->rad2)
1695 {
1696 val = -2.0 * Element->c[0] * dist - Element->c[1];
1697
1698 MTransNormal(V1, V1, Element->Trans);
1699
1700 Result += val * V1;
1701 }
1702
1703 break;
1704
1705 case BLOB_BASE_HEMISPHERE:
1706
1707 MInvTransPoint(PP, P, Element->Trans);
1708
1709 if (PP[Z] <= 0.0)
1710 {
1711 dist = PP.lengthSqr();
1712
1713 if (dist <= Element->rad2)
1714 {
1715 val = -2.0 * Element->c[0] * dist - Element->c[1];
1716
1717 MTransNormal(PP, PP, Element->Trans);
1718
1719 Result += val * PP;
1720 }
1721 }
1722
1723 break;
1724
1725 case BLOB_APEX_HEMISPHERE:
1726
1727 MInvTransPoint(PP, P, Element->Trans);
1728
1729 PP[Z] -= Element->len;
1730
1731 if (PP[Z] >= 0.0)
1732 {
1733 dist = PP.lengthSqr();
1734
1735 if (dist <= Element->rad2)
1736 {
1737 val = -2.0 * Element->c[0] * dist - Element->c[1];
1738
1739 MTransNormal(PP, PP, Element->Trans);
1740
1741 Result += val * PP;
1742 }
1743 }
1744
1745 break;
1746
1747 case BLOB_CYLINDER:
1748
1749 MInvTransPoint(PP, P, Element->Trans);
1750
1751 if ((PP[Z] >= 0.0) && (PP[Z] <= Element->len))
1752 {
1753 if ((dist = Sqr(PP[X]) + Sqr(PP[Y])) <= Element->rad2)
1754 {
1755 val = -2.0 * Element->c[0] * dist - Element->c[1];
1756
1757 PP[Z] = 0.0;
1758
1759 MTransNormal(PP, PP, Element->Trans);
1760
1761 Result += val * PP;
1762 }
1763 }
1764
1765 break;
1766 }
1767 }
1768
1769
1770
1771 /*****************************************************************************
1772 *
1773 * FUNCTION
1774 *
1775 * Blob_Normal
1776 *
1777 * INPUT
1778 *
1779 * Object - Pointer to blob structure
1780 * Inter - Pointer to intersection
1781 *
1782 * OUTPUT
1783 *
1784 * Result - Blob's normal
1785 *
1786 * RETURNS
1787 *
1788 * AUTHOR
1789 *
1790 * Alexander Enzmann
1791 *
1792 * DESCRIPTION
1793 *
1794 * Calculate the blob's surface normal in the intersection point.
1795 *
1796 * CHANGES
1797 *
1798 * Jul 1994 : Added code for bounding hierarchy traversal. [DB]
1799 *
1800 ******************************************************************************/
1801
Normal(Vector3d & Result,Intersection * Inter,TraceThreadData * Thread) const1802 void Blob::Normal(Vector3d& Result, Intersection *Inter, TraceThreadData *Thread) const
1803 {
1804 int i;
1805 unsigned int size;
1806 DBL dist, val;
1807 Vector3d New_Point, V1;
1808 BSPHERE_TREE *Tree;
1809 BSPHERE_TREE **Queue = reinterpret_cast<BSPHERE_TREE **>(Thread->Blob_Queue);
1810
1811 /* Transform the point into the blob space. */
1812 getLocalIPoint(New_Point, Inter);
1813
1814 Result = Vector3d(0.0, 0.0, 0.0);
1815
1816 /* For each component that contributes to this point, add its bit to the normal */
1817
1818 if (Data->Tree == nullptr)
1819 {
1820 /* There's no tree --> step through all elements. */
1821
1822 for (vector<Blob_Element>::iterator i = Data->Entry.begin(); i != Data->Entry.end(); ++i)
1823 {
1824 element_normal(Result, New_Point, &(*i));
1825 }
1826 }
1827 else
1828 {
1829 /* A tree exists --> step through the tree. */
1830
1831 size = 0;
1832
1833 Queue[size++] = Data->Tree;
1834
1835 while (size > 0)
1836 {
1837 Tree = Queue[--size];
1838
1839 /* Test if current node is a leaf. */
1840
1841 if (Tree->Entries <= 0)
1842 {
1843 element_normal(Result, New_Point, reinterpret_cast<Blob_Element *>(Tree->Node));
1844 }
1845 else
1846 {
1847 /* Test all sub-nodes. */
1848
1849 for (i = 0; i < (int)Tree->Entries; i++)
1850 {
1851 /* Insert sub-node if we are inside. */
1852
1853 V1 = New_Point - Tree->Node[i]->C;
1854
1855 dist = V1.lengthSqr();
1856
1857 if (dist <= Tree->Node[i]->r2)
1858 if (insert_node(Tree->Node[i], &size, Thread))
1859 Queue = reinterpret_cast<BSPHERE_TREE **>(Thread->Blob_Queue);
1860 }
1861 }
1862 }
1863 }
1864
1865 val = Result.lengthSqr();
1866
1867 if (val == 0.0)
1868 {
1869 Result = Vector3d(1.0, 0.0, 0.0);
1870 }
1871 else
1872 {
1873 /* Normalize normal vector. */
1874
1875 val = 1.0 / sqrt(val);
1876
1877 Result *= val;
1878 }
1879
1880 /* Transform back to world space. */
1881
1882 if (Trans != nullptr)
1883 {
1884 MTransNormal(Result, Result, Trans);
1885
1886 Result.normalize();
1887 }
1888 }
1889
1890
1891
1892 /*****************************************************************************
1893 *
1894 * FUNCTION
1895 *
1896 * Translate_Blob
1897 *
1898 * INPUT
1899 *
1900 * Vector - Translation vector
1901 *
1902 * OUTPUT
1903 *
1904 * Object - Pointer to blob structure
1905 *
1906 * RETURNS
1907 *
1908 * AUTHOR
1909 *
1910 * Alexander Enzmann
1911 *
1912 * DESCRIPTION
1913 *
1914 * Translate a blob.
1915 *
1916 * CHANGES
1917 *
1918 * -
1919 *
1920 ******************************************************************************/
1921
Translate(const Vector3d &,const TRANSFORM * tr)1922 void Blob::Translate(const Vector3d&, const TRANSFORM *tr)
1923 {
1924 Transform(tr);
1925 }
1926
1927
1928
1929 /*****************************************************************************
1930 *
1931 * FUNCTION
1932 *
1933 * Rotate_Blob
1934 *
1935 * INPUT
1936 *
1937 * Vector - Rotation vector
1938 *
1939 * OUTPUT
1940 *
1941 * Object - Pointer to blob structure
1942 *
1943 * RETURNS
1944 *
1945 * AUTHOR
1946 *
1947 * Alexander Enzmann
1948 *
1949 * DESCRIPTION
1950 *
1951 * Rotate a blob.
1952 *
1953 * CHANGES
1954 *
1955 * -
1956 *
1957 ******************************************************************************/
1958
Rotate(const Vector3d &,const TRANSFORM * tr)1959 void Blob::Rotate(const Vector3d&, const TRANSFORM *tr)
1960 {
1961 Transform(tr);
1962 }
1963
1964
1965
1966 /*****************************************************************************
1967 *
1968 * FUNCTION
1969 *
1970 * Scale_Blob
1971 *
1972 * INPUT
1973 *
1974 * Vector - Scaling vector
1975 *
1976 * OUTPUT
1977 *
1978 * Object - Pointer to blob structure
1979 *
1980 * RETURNS
1981 *
1982 * AUTHOR
1983 *
1984 * Alexander Enzmann
1985 *
1986 * DESCRIPTION
1987 *
1988 * Scale a blob.
1989 *
1990 * CHANGES
1991 *
1992 * -
1993 *
1994 ******************************************************************************/
1995
Scale(const Vector3d &,const TRANSFORM * tr)1996 void Blob::Scale(const Vector3d&, const TRANSFORM *tr)
1997 {
1998 Transform(tr);
1999 }
2000
2001
2002
2003 /*****************************************************************************
2004 *
2005 * FUNCTION
2006 *
2007 * Transform_Blob
2008 *
2009 * INPUT
2010 *
2011 * Trans - Pointer to transformation
2012 *
2013 * OUTPUT
2014 *
2015 * Object - Pointer to blob structure
2016 *
2017 * RETURNS
2018 *
2019 * AUTHOR
2020 *
2021 * Alexander Enzmann
2022 *
2023 * DESCRIPTION
2024 *
2025 * Transform a blob.
2026 *
2027 * CHANGES
2028 *
2029 * -
2030 *
2031 ******************************************************************************/
2032
Transform(const TRANSFORM * tr)2033 void Blob::Transform(const TRANSFORM *tr)
2034 {
2035 if (Trans == nullptr)
2036 Trans = Create_Transform();
2037
2038 Recompute_BBox(&BBox, tr);
2039
2040 Compose_Transforms(Trans, tr);
2041
2042 for(vector<TEXTURE*>::iterator i = Element_Texture.begin(); i != Element_Texture.end(); ++i)
2043 Transform_Textures(*i, tr);
2044 }
2045
2046
2047
2048 /*****************************************************************************
2049 *
2050 * FUNCTION
2051 *
2052 * Create_Blob
2053 *
2054 * INPUT
2055 *
2056 * Object - Pointer to blob structure
2057 *
2058 * OUTPUT
2059 *
2060 * Object
2061 *
2062 * RETURNS
2063 *
2064 * AUTHOR
2065 *
2066 * Alexander Enzmann
2067 *
2068 * DESCRIPTION
2069 *
2070 * Create a new blob.
2071 *
2072 * CHANGES
2073 *
2074 * -
2075 *
2076 ******************************************************************************/
2077
Blob()2078 Blob::Blob() : ObjectBase(BLOB_OBJECT)
2079 {
2080 Set_Flag(this, HIERARCHY_FLAG);
2081 Trans = nullptr;
2082 Data = nullptr;
2083 }
2084
2085 /*****************************************************************************
2086 *
2087 * FUNCTION
2088 *
2089 * Copy_Blob
2090 *
2091 * INPUT
2092 *
2093 * Object - Pointer to blob structure
2094 *
2095 * OUTPUT
2096 *
2097 * Object
2098 *
2099 * RETURNS
2100 *
2101 * AUTHOR
2102 *
2103 * Alexander Enzmann
2104 *
2105 * DESCRIPTION
2106 *
2107 * Copy a blob.
2108 *
2109 * NOTE: The components are not copied, only the number of references is
2110 * counted, so that Destroy_Blob() knows if they can be destroyed.
2111 *
2112 * CHANGES
2113 *
2114 * Jul 1994 : Added code for blob data reference counting. [DB]
2115 *
2116 ******************************************************************************/
2117
Copy()2118 ObjectPtr Blob::Copy()
2119 {
2120 Blob *New = new Blob();
2121
2122 /* Copy blob. */
2123
2124 Destroy_Transform(New->Trans);
2125 New->Trans = Copy_Transform(Trans);
2126 New->Data = Data->AcquireReference();
2127 New->Element_Texture.reserve(Element_Texture.size());
2128 for (vector<TEXTURE*>::iterator i = Element_Texture.begin(); i != Element_Texture.end(); ++i)
2129 New->Element_Texture.push_back(Copy_Textures(*i));
2130 return (New);
2131 }
2132
2133
2134
2135 /*****************************************************************************
2136 *
2137 * FUNCTION
2138 *
2139 * Create_Blob_List_Element
2140 *
2141 * INPUT
2142 *
2143 * OUTPUT
2144 *
2145 * RETURNS
2146 *
2147 * Blob_List_Struct * - Pointer to blob element
2148 *
2149 * AUTHOR
2150 *
2151 * Dieter Bayer
2152 *
2153 * DESCRIPTION
2154 *
2155 * Create a new blob element in the component list used during parsing.
2156 *
2157 * CHANGES
2158 *
2159 * Sep 1994 : Creation.
2160 *
2161 ******************************************************************************/
2162
Create_Blob_List_Element()2163 Blob_List_Struct *Blob::Create_Blob_List_Element()
2164 {
2165 return (new Blob_List_Struct);
2166 }
2167
2168
2169
2170 /*****************************************************************************
2171 *
2172 * FUNCTION
2173 *
2174 * Destroy_Blob
2175 *
2176 * INPUT
2177 *
2178 * Object - Pointer to blob structure
2179 *
2180 * OUTPUT
2181 *
2182 * Object
2183 *
2184 * RETURNS
2185 *
2186 * AUTHOR
2187 *
2188 * Alexander Enzmann
2189 *
2190 * DESCRIPTION
2191 *
2192 * Destroy a blob.
2193 *
2194 * NOTE: The blob data is destroyed if they are no longer used by any copy.
2195 *
2196 * CHANGES
2197 *
2198 * Jul 1994 : Added code for blob data reference counting. [DB]
2199 *
2200 * Dec 1994 : Fixed memory leakage. [DB]
2201 *
2202 * Aug 1995 : Fixed freeing of already freed memory. [DB]
2203 *
2204 ******************************************************************************/
2205
~Blob()2206 Blob::~Blob()
2207 {
2208 for (vector<TEXTURE*>::iterator i = Element_Texture.begin(); i != Element_Texture.end(); ++i)
2209 Destroy_Textures(*i);
2210 if (Data != nullptr)
2211 Data->ReleaseReference () ;
2212 }
2213
2214 /*****************************************************************************
2215 *
2216 * FUNCTION
2217 *
2218 * Compute_Blob_BBox
2219 *
2220 * INPUT
2221 *
2222 * Blob - Blob
2223 *
2224 * OUTPUT
2225 *
2226 * Blob
2227 *
2228 * RETURNS
2229 *
2230 * AUTHOR
2231 *
2232 * Dieter Bayer
2233 *
2234 * DESCRIPTION
2235 *
2236 * Calculate the bounding box of a blob.
2237 *
2238 * CHANGES
2239 *
2240 * Aug 1994 : Creation.
2241 *
2242 ******************************************************************************/
2243
Compute_BBox()2244 void Blob::Compute_BBox()
2245 {
2246 DBL radius, radius2;
2247 Vector3d Center, Min, Max;
2248
2249 Min = Vector3d(BOUND_HUGE);
2250 Max = Vector3d(-BOUND_HUGE);
2251
2252 for (vector<Blob_Element>::iterator i = Data->Entry.begin(); i != Data->Entry.end(); ++i)
2253 {
2254 if (i->c[2] > 0.0)
2255 {
2256 get_element_bounding_sphere(&(*i), Center, &radius2);
2257
2258 radius = sqrt(radius2);
2259
2260 Min[X] = min(Min[X], Center[X] - radius);
2261 Min[Y] = min(Min[Y], Center[Y] - radius);
2262 Min[Z] = min(Min[Z], Center[Z] - radius);
2263 Max[X] = max(Max[X], Center[X] + radius);
2264 Max[Y] = max(Max[Y], Center[Y] + radius);
2265 Max[Z] = max(Max[Z], Center[Z] + radius);
2266 }
2267 }
2268
2269 Make_BBox_from_min_max(BBox, Min, Max);
2270
2271 if (Trans != nullptr)
2272 {
2273 Recompute_BBox(&BBox, Trans);
2274 }
2275 }
2276
2277
2278
2279 /*****************************************************************************
2280 *
2281 * FUNCTION
2282 *
2283 * get_element_bounding_sphere
2284 *
2285 * INPUT
2286 *
2287 * Element - Pointer to element
2288 * Center - Bounding sphere's center
2289 * Radius2 - Bounding sphere's squared radius
2290 *
2291 * OUTPUT
2292 *
2293 * Center, Radius2
2294 *
2295 * RETURNS
2296 *
2297 * AUTHOR
2298 *
2299 * Dieter Bayer
2300 *
2301 * DESCRIPTION
2302 *
2303 * Calculate the bounding sphere of a blob element.
2304 *
2305 * CHANGES
2306 *
2307 * Sep 1994 : Creation.
2308 *
2309 ******************************************************************************/
2310
get_element_bounding_sphere(const Blob_Element * Element,Vector3d & Center,DBL * Radius2)2311 void Blob::get_element_bounding_sphere(const Blob_Element *Element, Vector3d& Center, DBL *Radius2)
2312 {
2313 DBL r, r2 = 0.0;
2314 Vector3d C;
2315 BoundingBox local_BBox;
2316
2317 switch (Element->Type)
2318 {
2319 case BLOB_SPHERE:
2320 case BLOB_ELLIPSOID:
2321
2322 r2 = Element->rad2;
2323
2324 C = Element->O;
2325
2326 break;
2327
2328 case BLOB_BASE_HEMISPHERE:
2329
2330 r2 = Element->rad2;
2331
2332 C = Vector3d(0.0, 0.0, 0.0);
2333
2334 break;
2335
2336 case BLOB_APEX_HEMISPHERE:
2337
2338 r2 = Element->rad2;
2339
2340 C = Vector3d(0.0, 0.0, Element->len);
2341
2342 break;
2343
2344 case BLOB_CYLINDER :
2345
2346 C = Vector3d(0.0, 0.0, 0.5 * Element->len);
2347
2348 r2 = Element->rad2 + Sqr(0.5 * Element->len);
2349
2350 break;
2351 }
2352
2353 /* Transform bounding sphere if necessary. */
2354 if (Element->Trans != nullptr)
2355 {
2356 r = sqrt(r2);
2357
2358 MTransPoint(C, C, Element->Trans);
2359
2360 Make_BBox(local_BBox, 0, 0, 0, r, r, r);
2361 Recompute_BBox(&local_BBox, Element->Trans);
2362 r = max(max(fabs(local_BBox.size[X]), fabs(local_BBox.size[Y])),
2363 fabs(local_BBox.size[Z]));
2364
2365 r2 = Sqr(r) + EPSILON;
2366 }
2367
2368 Center = C;
2369
2370 *Radius2 = r2;
2371 }
2372
2373
2374
2375 /*****************************************************************************
2376 *
2377 * FUNCTION
2378 *
2379 *
2380 * INPUT
2381 *
2382 *
2383 * OUTPUT
2384 *
2385 * RETURNS
2386 *
2387 * AUTHOR
2388 *
2389 * Dieter Bayer + others
2390 *
2391 * DESCRIPTION
2392 *
2393 *
2394 * CHANGES
2395 *
2396 * Sep 1994 : Creation.
2397 *
2398 ******************************************************************************/
2399
Blob_Element(void)2400 Blob_Element::Blob_Element (void)
2401 {
2402 Type = 0;
2403 index = 0;
2404 len = 0.0;
2405 rad2 = 0.0;
2406 c[0] = 0.0;
2407 c[1] = 0.0;
2408 c[2] = 0.0;
2409 Texture = nullptr;
2410 Trans = nullptr;
2411 }
2412
~Blob_Element()2413 Blob_Element::~Blob_Element ()
2414 {
2415 }
2416
Blob_Data(int Count)2417 Blob_Data::Blob_Data (int Count)
2418 {
2419 References = 1;
2420 Tree = nullptr;
2421 Entry.resize(Count);
2422 }
2423
AcquireReference(void)2424 Blob_Data *Blob_Data::AcquireReference (void)
2425 {
2426 References++;
2427 return (this) ;
2428 }
2429
ReleaseReference(void)2430 void Blob_Data::ReleaseReference (void)
2431 {
2432 if (--References == 0)
2433 {
2434 Destroy_Bounding_Sphere_Hierarchy(Tree);
2435
2436 /*
2437 * Make sure to destroy multiple references of a texture
2438 * and/or transformation only once. Multiple references
2439 * are only used with cylindrical blobs. Thus it's
2440 * enough to ignore all cylinder caps.
2441 */
2442 for (vector<Blob_Element>::iterator i = Entry.begin(); i != Entry.end(); ++i)
2443 if ((i->Type == BLOB_SPHERE) || (i->Type == BLOB_ELLIPSOID) || (i->Type == BLOB_CYLINDER))
2444 Destroy_Transform(i->Trans);
2445
2446 delete this ;
2447 }
2448 }
2449
~Blob_Data()2450 Blob_Data::~Blob_Data ()
2451 {
2452 // NOTE: ReleaseReference will call "delete this" if References == 1, so
2453 // we need to ensure we don't recurse when that happens. Only call
2454 // ReleaseReference() if References > 0.
2455 POV_SHAPE_ASSERT(References <= 1) ;
2456 if (References > 0)
2457 ReleaseReference () ;
2458 }
2459
2460 /*****************************************************************************
2461 *
2462 * FUNCTION
2463 *
2464 * Make_Blob
2465 *
2466 * INPUT
2467 *
2468 * Blob - Pointer to blob structure
2469 * threshold - Blob's threshold
2470 * BlobList - Pointer to elements
2471 * npoints - Number of elements
2472 *
2473 * OUTPUT
2474 *
2475 * Blob
2476 *
2477 * RETURNS
2478 *
2479 * AUTHOR
2480 *
2481 * Alexander Enzmann
2482 *
2483 * DESCRIPTION
2484 *
2485 * Create a blob after it was read from the scene file.
2486 *
2487 * Starting with the density function: (1-r^2)^2, we have a field
2488 * that varies in strength from 1 at r = 0 to 0 at r = 1. By
2489 * substituting r/rad for r, we can adjust the range of influence
2490 * of a particular component. By multiplication by coeff, we can
2491 * adjust the amount of total contribution, giving the formula:
2492 *
2493 * coeff * (1 - (r/rad)^2)^2
2494 *
2495 * This varies in strength from coeff at r = 0, to 0 at r = rad.
2496 *
2497 * CHANGES
2498 *
2499 * Jul 1994 : Added code for cylindrical and ellipsoidical blobs. [DB]
2500 *
2501 ******************************************************************************/
2502
Make_Blob(DBL threshold,Blob_List_Struct * BlobList,int npoints,TraceThreadData * Thread)2503 int Blob::Make_Blob(DBL threshold, Blob_List_Struct *BlobList, int npoints, TraceThreadData *Thread)
2504 {
2505 int i, count;
2506 DBL rad2, coeff;
2507 Blob_List_Struct *temp;
2508
2509 if (npoints < 1)
2510 throw POV_EXCEPTION_STRING("Need at least one component in a blob.");
2511
2512 /* Figure out how many components there will be. */
2513
2514 for (i = 0, count = npoints, temp = BlobList; i < npoints; i++, temp = temp->next)
2515 if (temp->elem.Type & BLOB_CYLINDER)
2516 count += 2;
2517
2518 #ifdef MAX_BLOB_COMPONENTS
2519 /* Test for too many components. [DB 12/94] */
2520 if (count >= MAX_BLOB_COMPONENTS)
2521 throw POV_EXCEPTION_STRING("There are more than the maximum supported components in a blob.");
2522 #endif
2523
2524 /* Initialize the blob data. */
2525
2526 Data->Threshold = threshold;
2527 Data->Entry.resize(count) ;
2528
2529 vector<Blob_Element>::iterator it = Data->Entry.begin() ;
2530 for (i = 0; i < npoints; i++)
2531 {
2532 temp = BlobList;
2533 POV_SHAPE_ASSERT(it != Data->Entry.end()) ;
2534 Blob_Element* Entry = &*it++ ;
2535
2536 if ((fabs(temp->elem.c[2]) < EPSILON) || (temp->elem.rad2 < EPSILON))
2537 {
2538 ;// TODO MESSAGE Warning("Degenerate Blob element");
2539 }
2540
2541 /* Initialize component. */
2542 *Entry = temp->elem;
2543
2544 /* We have a multi-texture blob. */
2545 if (Entry->Texture != nullptr)
2546 Set_Flag(this, MULTITEXTURE_FLAG);
2547
2548 /* Store blob specific information. */
2549 rad2 = temp->elem.rad2;
2550 coeff = temp->elem.c[2];
2551 Entry->c[0] = coeff / (rad2 * rad2);
2552 Entry->c[1] = -(2.0 * coeff) / rad2;
2553 Entry->c[2] = coeff;
2554
2555 if (temp->elem.Type == BLOB_CYLINDER)
2556 {
2557 /* Create hemispherical component at the base. */
2558 Entry = &*it++ ;
2559 *Entry = temp->elem;
2560 Entry->Type = BLOB_BASE_HEMISPHERE;
2561 Entry->c[0] = coeff / (rad2 * rad2);
2562 Entry->c[1] = -(2.0 * coeff) / rad2;
2563 Entry->c[2] = coeff;
2564
2565 /* Create hemispherical component at the apex. */
2566 Entry = &*it++ ;
2567 *Entry = temp->elem;
2568 Entry->Type = BLOB_APEX_HEMISPHERE;
2569 Entry->c[0] = coeff / (rad2 * rad2);
2570 Entry->c[1] = -(2.0 * coeff) / rad2;
2571 Entry->c[2] = coeff;
2572 }
2573
2574 /* Get rid of texture non longer needed. */
2575
2576 BlobList = BlobList->next;
2577 Destroy_Textures(temp->elem.Texture);
2578 delete temp ;
2579 }
2580
2581 for (i = 0; i < count; i++)
2582 Data->Entry[i].index = i;
2583
2584 /* Compute bounding box. */
2585
2586 Compute_BBox();
2587
2588 /* Create bounding sphere hierarchy. */
2589
2590 if (Test_Flag(this, HIERARCHY_FLAG))
2591 build_bounding_hierarchy();
2592
2593 if (count * 5 >= Thread->Blob_Coefficient_Count)
2594 {
2595 POV_FREE(Thread->Blob_Coefficients);
2596 Thread->Blob_Coefficient_Count = count * 7;
2597 Thread->Blob_Coefficients = reinterpret_cast<DBL *>(POV_MALLOC(sizeof(DBL) * Thread->Blob_Coefficient_Count, "Blob Coefficients"));
2598 }
2599
2600 if (Data->Entry.size() * 2 >= Thread->Blob_Interval_Count)
2601 {
2602 delete[] Thread->Blob_Intervals;
2603 Thread->Blob_Interval_Count = Data->Entry.size() * 5 / 2;
2604 Thread->Blob_Intervals = new Blob_Interval_Struct [Thread->Blob_Interval_Count];
2605 }
2606
2607 return (count) ;
2608 }
2609
2610 /*****************************************************************************
2611 *
2612 * FUNCTION
2613 *
2614 * Test_Blob_Opacity
2615 *
2616 * INPUT
2617 *
2618 * Blob - Pointer to blob structure
2619 *
2620 * OUTPUT
2621 *
2622 * Blob
2623 *
2624 * RETURNS
2625 *
2626 * AUTHOR
2627 *
2628 * Dieter Bayer
2629 *
2630 * DESCRIPTION
2631 *
2632 * Set the opacity flag of the blob according to the opacity
2633 * of the blob's texture(s).
2634 *
2635 * CHANGES
2636 *
2637 * Apr 1996 : Creation.
2638 *
2639 ******************************************************************************/
2640
IsOpaque() const2641 bool Blob::IsOpaque() const
2642 {
2643 if (Test_Flag(this, MULTITEXTURE_FLAG))
2644 {
2645 for (auto&& elementTexture : Element_Texture)
2646 {
2647 // If component's texture isn't opaque the blob is neither.
2648 if ((elementTexture != nullptr) && !Test_Opacity(elementTexture))
2649 return false;
2650 }
2651 }
2652
2653 // Otherwise it's a question of whether the common texture is opaque or not.
2654 // TODO FIXME - other objects report as non-opaque if Texture == nullptr.
2655 // TODO FIXME - other objects report as non-opaque if Interior_Texture present and non-opaque.
2656 // What we probably really want here is `return ObjectBase::IsOpaque()`.
2657 return (Texture == nullptr) || Test_Opacity(Texture);
2658 }
2659
2660
2661
2662 /*****************************************************************************
2663 *
2664 * FUNCTION
2665 *
2666 * build_bounding_hierarchy
2667 *
2668 * INPUT
2669 *
2670 * OUTPUT
2671 *
2672 * RETURNS
2673 *
2674 * AUTHOR
2675 *
2676 * Dieter Bayer
2677 *
2678 * DESCRIPTION
2679 *
2680 * Create the bounding sphere hierarchy.
2681 *
2682 * CHANGES
2683 *
2684 * Oct 1994 : Creation. (Derived from the bounding slab creation code)
2685 *
2686 ******************************************************************************/
2687
build_bounding_hierarchy()2688 void Blob::build_bounding_hierarchy()
2689 {
2690 int i, nElem, maxelements;
2691 BSPHERE_TREE **Elements;
2692
2693 nElem = (int)Data->Entry.size();
2694
2695 maxelements = 2 * nElem;
2696
2697 /*
2698 * Now allocate an array to hold references to these elements.
2699 */
2700
2701 Elements = reinterpret_cast<BSPHERE_TREE **>(POV_MALLOC(maxelements*sizeof(BSPHERE_TREE *), "blob bounding hierarchy"));
2702
2703 /* Init list with blob elements. */
2704
2705 for (i = 0; i < nElem; i++)
2706 {
2707 Elements[i] = reinterpret_cast<BSPHERE_TREE *>(POV_MALLOC(sizeof(BSPHERE_TREE), "blob bounding hierarchy"));
2708
2709 Elements[i]->Entries = 0;
2710 Elements[i]->Node = reinterpret_cast<BSPHERE_TREE **>(&Data->Entry[i]);
2711
2712 get_element_bounding_sphere(&Data->Entry[i], Elements[i]->C, &Elements[i]->r2);
2713 }
2714
2715 Build_Bounding_Sphere_Hierarchy(&Data->Tree, nElem, &Elements);
2716
2717 /* Get rid of the Elements array. */
2718
2719 POV_FREE(Elements);
2720 }
2721
2722
2723
2724 /*****************************************************************************
2725 *
2726 * FUNCTION
2727 *
2728 * Determine_Blob_Textures
2729 *
2730 * INPUT
2731 *
2732 * OUTPUT
2733 *
2734 * RETURNS
2735 *
2736 * AUTHOR
2737 *
2738 * Dieter Bayer
2739 *
2740 * DESCRIPTION
2741 *
2742 * Determine the textures and weights of all components affecting
2743 * the given intersection point. The weights are calculated from
2744 * the field values and sum to 1.
2745 *
2746 * CHANGES
2747 *
2748 * Sep 1994 : Creation.
2749 *
2750 * Mar 1996 : Make the call to resize the textures/weights list just once
2751 * at the beginning instead of doing it for every element. [DB]
2752 *
2753 ******************************************************************************/
2754
Determine_Textures(Intersection * isect,bool hitinside,WeightedTextureVector & textures,TraceThreadData * Thread)2755 void Blob::Determine_Textures(Intersection *isect, bool hitinside, WeightedTextureVector& textures, TraceThreadData *Thread)
2756 {
2757 unsigned int size;
2758 DBL rad2;
2759 Vector3d V1, P;
2760 Blob_Element *Element;
2761 BSPHERE_TREE *Tree;
2762 size_t firstinserted = textures.size();
2763 BSPHERE_TREE **Queue = reinterpret_cast<BSPHERE_TREE **>(Thread->Blob_Queue);
2764
2765 /* Transform the point into the blob space. */
2766 getLocalIPoint(P, isect);
2767
2768 if (Data->Tree == nullptr)
2769 {
2770 /* There's no tree --> step through all elements. */
2771
2772 for (int i = 0; i < Data->Entry.size(); i++)
2773 {
2774 Element = &Data->Entry[i];
2775 determine_element_texture(Element, Element_Texture[i], P, textures);
2776 }
2777 }
2778 else
2779 {
2780 /* A tree exists --> step through the tree. */
2781
2782 size = 0;
2783
2784 Queue[size++] = Data->Tree;
2785
2786 while (size > 0)
2787 {
2788 Tree = Queue[--size];
2789
2790 /* Test if current node is a leaf. */
2791
2792 if (Tree->Entries <= 0)
2793 {
2794 determine_element_texture(reinterpret_cast<Blob_Element *>(Tree->Node), Element_Texture[((Blob_Element *)Tree->Node)->index], P, textures);
2795 }
2796 else
2797 {
2798 /* Test all sub-nodes. */
2799
2800 for (int i = 0; i < (int)Tree->Entries; i++)
2801 {
2802 /* Insert sub-node if we are inside. */
2803
2804 V1 = P - Tree->Node[i]->C;
2805
2806 rad2 = V1.lengthSqr();
2807
2808 if (rad2 <= Tree->Node[i]->r2)
2809 if (insert_node(Tree->Node[i], &size, Thread))
2810 Queue = reinterpret_cast<BSPHERE_TREE **>(Thread->Blob_Queue);
2811 }
2812 }
2813 }
2814 }
2815
2816 /* Normalize weights so that their sum is 1. */
2817
2818 if((textures.size() - firstinserted) > 0)
2819 {
2820 COLC sum = 0.0;
2821
2822 for(size_t i = firstinserted; i < textures.size(); i++)
2823 sum += textures[i].weight;
2824
2825 sum = 1.0 / sum;
2826
2827 for(size_t i = firstinserted; i < textures.size(); i++)
2828 textures[i].weight *= sum;
2829 }
2830 }
2831
2832
2833 /*****************************************************************************
2834 *
2835 * FUNCTION
2836 *
2837 * determine_element_texture
2838 *
2839 * INPUT
2840 *
2841 * OUTPUT
2842 *
2843 * RETURNS
2844 *
2845 * AUTHOR
2846 *
2847 * Dieter Bayer
2848 *
2849 * DESCRIPTION
2850 *
2851 * If the intersection point is inside the component calculate
2852 * the field density and store the element's texture and the field
2853 * value in the texture/weight list.
2854 *
2855 * CHANGES
2856 *
2857 * Sep 1994 : Creation.
2858 *
2859 ******************************************************************************/
2860
determine_element_texture(const Blob_Element * Element,TEXTURE * ElementTex,const Vector3d & P,WeightedTextureVector & textures)2861 void Blob::determine_element_texture(const Blob_Element *Element, TEXTURE *ElementTex, const Vector3d& P, WeightedTextureVector& textures)
2862 {
2863 DBL density = fabs(calculate_element_field(Element, P));
2864
2865 if(density > 0.0)
2866 textures.push_back(WeightedTexture(density, ElementTex != nullptr ? ElementTex : Texture));
2867 }
2868
2869
2870
2871 /*****************************************************************************
2872 *
2873 * FUNCTION
2874 *
2875 * Translate_Blob_Element
2876 *
2877 * INPUT
2878 *
2879 * Element - Pointer to blob element
2880 * Vector - Translation vector
2881 *
2882 * OUTPUT
2883 *
2884 * Object
2885 *
2886 * RETURNS
2887 *
2888 * AUTHOR
2889 *
2890 * Dieter Bayer
2891 *
2892 * DESCRIPTION
2893 *
2894 * Translate a blob element.
2895 *
2896 * CHANGES
2897 *
2898 * Sep 1994 : Creation.
2899 *
2900 ******************************************************************************/
2901
Translate_Blob_Element(Blob_Element * Element,const Vector3d & Vector)2902 void Blob::Translate_Blob_Element(Blob_Element *Element, const Vector3d& Vector)
2903 {
2904 TRANSFORM Trans;
2905
2906 Compute_Translation_Transform(&Trans, Vector);
2907
2908 if (Element->Trans == nullptr)
2909 {
2910 /* This is a sphere component. */
2911
2912 Element->O += Vector;
2913 Transform_Textures(Element->Texture, &Trans);
2914 }
2915 else
2916 {
2917 /* This is one of the other components. */
2918
2919 Transform_Blob_Element(Element, &Trans);
2920 }
2921 }
2922
2923
2924
2925 /*****************************************************************************
2926 *
2927 * FUNCTION
2928 *
2929 * Rotate_Blob_Element
2930 *
2931 * INPUT
2932 *
2933 * Element - Pointer to blob element
2934 * Vector - Translation vector
2935 *
2936 * OUTPUT
2937 *
2938 * Object
2939 *
2940 * RETURNS
2941 *
2942 * AUTHOR
2943 *
2944 * Dieter Bayer
2945 *
2946 * DESCRIPTION
2947 *
2948 * Rotate a blob element.
2949 *
2950 * CHANGES
2951 *
2952 * Sep 1994 : Creation.
2953 *
2954 ******************************************************************************/
2955
Rotate_Blob_Element(Blob_Element * Element,const Vector3d & Vector)2956 void Blob::Rotate_Blob_Element(Blob_Element *Element, const Vector3d& Vector)
2957 {
2958 TRANSFORM Trans;
2959
2960 Compute_Rotation_Transform(&Trans, Vector);
2961
2962 if (Element->Trans == nullptr)
2963 {
2964 /* This is a sphere component. */
2965
2966 MTransPoint(Element->O, Element->O, &Trans);
2967 Transform_Textures(Element->Texture, &Trans);
2968 }
2969 else
2970 {
2971 /* This is one of the other components. */
2972
2973 Transform_Blob_Element(Element, &Trans);
2974 }
2975 }
2976
2977
2978
2979 /*****************************************************************************
2980 *
2981 * FUNCTION
2982 *
2983 * Scale_Blob_Element
2984 *
2985 * INPUT
2986 *
2987 * Element - Pointer to blob element
2988 * Vector - Translation vector
2989 *
2990 * OUTPUT
2991 *
2992 * Object
2993 *
2994 * RETURNS
2995 *
2996 * AUTHOR
2997 *
2998 * Dieter Bayer
2999 *
3000 * DESCRIPTION
3001 *
3002 * Scale a blob element.
3003 *
3004 * CHANGES
3005 *
3006 * Sep 1994 : Creation.
3007 *
3008 ******************************************************************************/
3009
Scale_Blob_Element(Blob_Element * Element,const Vector3d & Vector)3010 void Blob::Scale_Blob_Element(Blob_Element *Element, const Vector3d& Vector)
3011 {
3012 TRANSFORM Trans;
3013
3014 if ((Vector[X] != Vector[Y]) || (Vector[X] != Vector[Z]))
3015 {
3016 if (Element->Trans == nullptr)
3017 {
3018 /* This is a sphere component --> change to ellipsoid component. */
3019
3020 Element->Type = BLOB_ELLIPSOID;
3021
3022 Element->Trans = Create_Transform();
3023 }
3024 }
3025
3026 Compute_Scaling_Transform(&Trans, Vector);
3027
3028 if (Element->Trans == nullptr)
3029 {
3030 /* This is a sphere component. */
3031
3032 Element->O *= Vector[X];
3033
3034 Element->rad2 *= Sqr(Vector[X]);
3035
3036 Transform_Textures(Element->Texture, &Trans);
3037 }
3038 else
3039 {
3040 /* This is one of the other components. */
3041
3042 Transform_Blob_Element(Element, &Trans);
3043 }
3044 }
3045
3046
3047
3048 /*****************************************************************************
3049 *
3050 * FUNCTION
3051 *
3052 * Transform_Blob_Element
3053 *
3054 * INPUT
3055 *
3056 * Element - Pointer to blob element
3057 * Trans - Transformation
3058 *
3059 * OUTPUT
3060 *
3061 * Object
3062 *
3063 * RETURNS
3064 *
3065 * AUTHOR
3066 *
3067 * Dieter Bayer
3068 *
3069 * DESCRIPTION
3070 *
3071 * Transform a blob element.
3072 *
3073 * CHANGES
3074 *
3075 * Sep 1994 : Creation.
3076 *
3077 ******************************************************************************/
3078
Transform_Blob_Element(Blob_Element * Element,const TRANSFORM * Trans)3079 void Blob::Transform_Blob_Element(Blob_Element *Element, const TRANSFORM *Trans)
3080 {
3081 if (Element->Trans == nullptr)
3082 {
3083 /* This is a sphere component --> change to ellipsoid component. */
3084
3085 Element->Type = BLOB_ELLIPSOID;
3086
3087 Element->Trans = Create_Transform();
3088 }
3089
3090 Compose_Transforms(Element->Trans, Trans);
3091
3092 Transform_Textures(Element->Texture, Trans);
3093 }
3094
3095
3096
3097 /*****************************************************************************
3098 *
3099 * FUNCTION
3100 *
3101 * Invert_Blob_Element
3102 *
3103 * INPUT
3104 *
3105 * Element - Pointer to blob element
3106 *
3107 * OUTPUT
3108 *
3109 * Object
3110 *
3111 * RETURNS
3112 *
3113 * AUTHOR
3114 *
3115 * Dieter Bayer
3116 *
3117 * DESCRIPTION
3118 *
3119 * Invert blob element by negating its strength.
3120 *
3121 * CHANGES
3122 *
3123 * Sep 1994 : Creation.
3124 *
3125 ******************************************************************************/
3126
Invert_Blob_Element(Blob_Element * Element)3127 void Blob::Invert_Blob_Element(Blob_Element *Element)
3128 {
3129 Element->c[2] *= -1.0;
3130 }
3131
3132
3133 /*****************************************************************************
3134 *
3135 * FUNCTION
3136 *
3137 * insert_node
3138 *
3139 * INPUT
3140 *
3141 * OUTPUT
3142 *
3143 * RETURNS
3144 *
3145 * AUTHOR
3146 *
3147 * Dieter Bayer
3148 *
3149 * DESCRIPTION
3150 *
3151 * Insert a node into the node queue.
3152 *
3153 * CHANGES
3154 *
3155 * Feb 1995 : Creation.
3156 *
3157 ******************************************************************************/
3158
insert_node(BSPHERE_TREE * Node,unsigned int * size,TraceThreadData * Thread)3159 bool Blob::insert_node(BSPHERE_TREE *Node, unsigned int *size, TraceThreadData *Thread)
3160 {
3161 /* Resize queue if necessary. */
3162 bool rval = false ;
3163 BSPHERE_TREE **Queue = reinterpret_cast<BSPHERE_TREE **>(Thread->Blob_Queue);
3164
3165 if (*size >= Thread->Max_Blob_Queue_Size)
3166 {
3167 Thread->Max_Blob_Queue_Size = (*size + 1) * 3 / 2;
3168 Queue = reinterpret_cast<BSPHERE_TREE **>(POV_REALLOC(Queue, Thread->Max_Blob_Queue_Size*sizeof(BSPHERE_TREE *), "blob queue"));
3169 Thread->Blob_Queue = reinterpret_cast<void **>(Queue);
3170 rval = true ;
3171 }
3172
3173 Queue[(*size)++] = Node;
3174 return (rval) ;
3175 }
3176
3177
3178
3179 /*****************************************************************************
3180 *
3181 * FUNCTION
3182 *
3183 * Create_Blob_Element_Texture_List
3184 *
3185 * INPUT
3186 *
3187 * OUTPUT
3188 *
3189 * RETURNS
3190 *
3191 * AUTHOR
3192 *
3193 * Dieter Bayer
3194 *
3195 * DESCRIPTION
3196 *
3197 * Create a list of all textures in the blob.
3198 *
3199 * The list actually contains copies of the textures not
3200 * just references to them.
3201 *
3202 * CHANGES
3203 *
3204 * Mar 1996 : Created.
3205 *
3206 ******************************************************************************/
3207
Create_Blob_Element_Texture_List(Blob_List_Struct * BlobList,int npoints)3208 void Blob::Create_Blob_Element_Texture_List(Blob_List_Struct *BlobList, int npoints)
3209 {
3210 int i, count;
3211 Blob_List_Struct *bl;
3212
3213 if (npoints < 1)
3214 throw POV_EXCEPTION_STRING("Need at least one component in a blob.");
3215
3216 /* Figure out how many components there will be. */
3217 for (i = 0, count = npoints, bl = BlobList ; i < npoints; i++, bl = bl->next)
3218 if (bl->elem.Type & BLOB_CYLINDER)
3219 count += 2;
3220
3221 #ifdef MAX_BLOB_COMPONENTS
3222 /* Test for too many components. [DB 12/94] */
3223 if (count >= MAX_BLOB_COMPONENTS)
3224 throw POV_EXCEPTION_STRING("There are more than the maximum supported components in a blob.");
3225 #endif
3226
3227 Data = new Blob_Data (count) ;
3228
3229 /* Allocate memory for list. */
3230 Element_Texture.reserve(count);
3231 for (i = 0, bl = BlobList; i < npoints; i++, bl = bl->next)
3232 {
3233 /*
3234 * Copy texture into element texture list. This is neccessary
3235 * because individual textures have to be transformed too if
3236 * copies of the blob are transformed.
3237 */
3238 Element_Texture.push_back(Copy_Textures(bl->elem.Texture));
3239 if (bl->elem.Type == BLOB_CYLINDER)
3240 {
3241 Element_Texture.push_back(Copy_Textures(bl->elem.Texture));
3242 Element_Texture.push_back(Copy_Textures(bl->elem.Texture));
3243 }
3244 }
3245 }
3246
3247 /*****************************************************************************/
3248
getLocalIPoint(Vector3d & lip,Intersection * isect) const3249 void Blob::getLocalIPoint(Vector3d& lip, Intersection *isect) const
3250 {
3251 if(isect->haveLocalIPoint == false)
3252 {
3253 if (Trans != nullptr)
3254 MInvTransPoint(isect->LocalIPoint, isect->IPoint, Trans);
3255 else
3256 isect->LocalIPoint = isect->IPoint;
3257
3258 isect->haveLocalIPoint = true;
3259 }
3260
3261 lip = isect->LocalIPoint;
3262 }
3263
3264 }
3265