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