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