/************************************************************************ **** * blob.cpp * * This module implements functions that manipulate blobs. * * The original file was written by Alexander Enzmann. * He wrote the code for blobs and generously provided us these enhancements. * * Modifications and enhancements by Dieter Bayer [DB]. * * from Persistence of Vision(tm) Ray Tracer * Copyright 1996-2002 Persistence of Vision Team *------------------------------------------------------------------------ --- * NOTICE: This source code file is provided so that users may experiment * with enhancements to POV-Ray and to port the software to platforms other * than those supported by the POV-Ray Team. There are strict rules under * which you are permitted to use this file. The rules are in the file * named POVLEGAL.DOC which should be distributed with this file. * If POVLEGAL.DOC is not available it may be found online at - * * http://www.povray.org/povlegal.html. * * This program is based on the popular DKB raytracer version 2.12. * DKBTrace was originally written by David K. Buck. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins. * * $File: //depot/povray/3.5/source/blob.cpp $ * $Revision: #12 $ * $Change: 1817 $ * $DateTime: 2002/07/27 10:45:37 $ * $Author: chrisc $ * $Log$ * ************************************************************************* ****/ /************************************************************************ **** * * Explanation: * * - * * Syntax: * * blob * { * threshold THRESHOLD_VALUE * * component STRENGTH, RADIUS, * * sphere { , RADIUS, [strength] STRENGTH * [ translate ] * [ rotate ] * [ scale ] * [ finish { ... } ] * [ pigment { ... } ] * [ tnormal { ... } ] * [ texture { ... } ] * } * * cylinder { , , RADIUS, [strength] STRENGTH * [ translate ] * [ rotate ] * [ scale ] * [ finish { ... } ] * [ pigment { ... } ] * [ tnormal { ... } ] * [ texture { ... } ] * } * * [ sturm ] * [ hierarchy FLAG ] * } * * --- * * Jul 1994 : Most functions rewritten, bounding hierarchy added. [DB] * * Aug 1994 : Cylindrical blobs added. [DB] * * Sep 1994 : Multi-texturing added (each component can have its own texture). * Translation, rotation and scaling of each component added. [DB] * * Oct 1994 : Adopted the method for the bounding slab creation to build the * bounding sphere hierarchy of the blob to get a much better * hierarchy. Improved bounding sphere calculation for tighter * bounds. [DB] * * Dec 1994 : Added code for dynamic blob queue allocation. [DB] * * Feb 1995 : Moved bounding sphere stuff into a seperate file. [DB] * ************************************************************************* ****/ #include "frame.h" #include "povray.h" #include "vector.h" #include "povproto.h" #include "blob.h" #include "bbox.h" #include "bsphere.h" #include "lighting.h" #include "matrices.h" #include "objects.h" #include "polysolv.h" #include "texture.h" #include using namespace std; /************************************************************************ ***** * Local preprocessor defines ************************************************************************* *****/ /* Minimal intersection depth for a valid intersection. */ #define DEPTH_TOLERANCE 1.0e-2 /* Tolerance for inside test. */ #define INSIDE_TOLERANCE 1.0e-6 /* Ray enters/exits a component. */ #define ENTERING 0 #define EXITING 1 /************************************************************************ ***** * Local typedefs ************************************************************************* *****/ /************************************************************************ ***** * Static functions ************************************************************************* *****/ static void element_normal (VECTOR Result, VECTOR P, BLOB_ELEMENT *Element); static int intersect_element (VECTOR P, VECTOR D, BLOB_ELEMENT *Element, DBL mindist, DBL *t0, DBL *t1); static void insert_hit (BLOB_ELEMENT *Element, DBL t0, DBL t1, BLOB_INTERVAL *intervals, int *cnt); static int determine_influences (VECTOR P, VECTOR D, BLOB *Blob, DBL mindist, BLOB_INTERVAL *intervals); static DBL calculate_field_value (BLOB *Blob, VECTOR P); static DBL calculate_element_field (BLOB_ELEMENT *Element, VECTOR P); static int intersect_cylinder (BLOB_ELEMENT *Element, VECTOR P, VECTOR D, DBL mindist, DBL *tmin, DBL *tmax); static int intersect_hemisphere (BLOB_ELEMENT *Element, VECTOR P, VECTOR D, DBL mindist, DBL *tmin, DBL *tmax); static int intersect_sphere (BLOB_ELEMENT *Element, VECTOR P, VECTOR D, DBL mindist, DBL *tmin, DBL *tmax); static int intersect_ellipsoid (BLOB_ELEMENT *Element, VECTOR P, VECTOR D, DBL mindist, DBL *tmin, DBL *tmax); static void get_element_bounding_sphere (BLOB_ELEMENT *Element, VECTOR Center, DBL *Radius2); static void build_bounding_hierarchy (BLOB *Blob); static void init_blob_element (BLOB_ELEMENT *Element); static void determine_element_texture (BLOB *Blob, BLOB_ELEMENT *Element, TEXTURE *Texture, VECTOR P, int *Count, TEXTURE **Textures, DBL *Weights); static void insert_node (BSPHERE_TREE *Node, int *size); static int All_Blob_Intersections (OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack); static int Inside_Blob (VECTOR point, OBJECT *Object); static void Blob_Normal (VECTOR Result, OBJECT *Object, INTERSECTION *Inter); static BLOB *Copy_Blob (OBJECT *Object); static void Translate_Blob (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans); static void Rotate_Blob (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans); static void Scale_Blob (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans); static void Invert_Blob (OBJECT *Object); static void Transform_Blob (OBJECT *Object, TRANSFORM *Trans); static void Destroy_Blob (OBJECT *Object); static void Compute_Blob_BBox (BLOB *Blob); /************************************************************************ ***** * Local variables ************************************************************************* *****/ METHODS Blob_Methods = { All_Blob_Intersections, Inside_Blob, Blob_Normal, Default_UVCoord, (COPY_METHOD)Copy_Blob, Translate_Blob, Rotate_Blob, Scale_Blob, Transform_Blob, Invert_Blob, Destroy_Blob }; static BSPHERE_TREE **Queue; /* Maximum number of entries in queue during hierarchy traversal. */ static unsigned Max_Queue_Size = 1024; /************************************************************************ ***** * * FUNCTION * * All_Blob_Intersections * * INPUT * * Object - Object * Ray - Ray * * OUTPUT * * Depth_Stack - Intersection stack * * RETURNS * * int - true, if a intersection was found * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Generate intervals of influence for each component. After these * are made, determine their aggregate effect on the ray. As the * individual intervals are checked, a quartic is generated * that represents the density at a particular point on the ray. * * Explanation for spherical components: * * After making the substitutions in MakeBlob, there is a formula * for each component that has the form: * * c0 * r^4 + c1 * r^2 + c2. * * In order to determine the influence on the ray of all of the * individual components, we start by determining the distance * from any point on the ray to the specified point. This can * be found using the pythagorean theorem, using C as the center * of this component, P as the start of the ray, and D as the * direction of travel of the ray: * * r^2 = (t * D + P - C) . (t * D + P - C) * * we insert this equation for each appearance of r^2 in the * components' formula, giving: * * r^2 = D.D t^2 + 2 t D . (P - C) + (P - C) . (P - C) * * Since the direction vector has been normalized, D.D = 1. * Using the substitutions: * * t0 = (P - C) . (P - C), * t1 = D . (P - C) * * We can write the formula as: * * r^2 = t0 + 2 t t1 + t^2 * * Taking r^2 and substituting into the formula for this component * of the Blob we get the formula: * * density = c0 * (r^2)^2 + c1 * r^2 + c2, * * or: * * density = c0 * (t0 + 2 t t1 + t^2)^2 + * c1 * (t0 + 2 t t1 + t^2) + * c2 * * Expanding terms and collecting with respect to "t" gives: * * t^4 * c0 + * t^3 * 4 c0 t1 + * t^2 * (c1 + 2 * c0 t0 + 4 c0 t1^2) * t * 2 (c1 t1 + 2 c0 t0 t1) + * c2 + c1*t0 + c0*t0^2 * * This formula can now be solved for "t" by any of the quartic * root solvers that are available. * * CHANGES * * Jul 1994 : Added code for cylindrical and ellipsoidical blobs. [DB] * * Oct 1994 : Added code to convert polynomial into a bezier curve for * a quick test if an intersection exists in an interval. [DB] * * Sep 1995 : Added code to avoid numerical problems with distant blobs. [DB] * ************************************************************************* *****/ static int All_Blob_Intersections(OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack) { int i, j, cnt; int root_count, in_flag; int Intersection_Found = false; DBL t0, t1, t2, c0, c1, c2, dist, len, start_dist; DBL *fcoeffs; DBL coeffs[5]; DBL roots[4]; VECTOR P, D, V1, PP, DD; VECTOR IPoint; BLOB_INTERVAL *intervals; BLOB_ELEMENT *Element; BLOB *Blob = (BLOB *)Object; DBL newcoeffs[5], dk[5]; DBL l, w; Increase_Counter(stats[Ray_Blob_Tests]); /* Transform the ray into blob space. */ if (Blob->Trans != NULL) { MInvTransPoint(P, Ray->Initial, Blob->Trans); MInvTransDirection(D, Ray->Direction, Blob->Trans); VLength(len, D); VInverseScaleEq(D, len); } else { Assign_Vector(P, Ray->Initial); Assign_Vector(D, Ray->Direction); len = 1.0; } /* Get the intervals along the ray where each component has an effect. */ intervals = Blob->Data->Intervals; if ((cnt = determine_influences(P, D, Blob, DEPTH_TOLERANCE, intervals)) == 0) { /* Ray doesn't hit any of the component elements. */ return (false); } /* To avoid numerical problems we start at the first interval. */ if ((start_dist = intervals[0].bound) < Small_Tolerance) { start_dist = 0.0; } DBL Max_Bound = intervals[0].bound; for (i = 0; i < cnt; i++) { if (intervals[i].bound > Max_Bound) Max_Bound = intervals[i].bound; intervals[i].bound -= start_dist; } /* Get the new starting point. */ VAddScaledEq(P, start_dist, D); if (Max_Bound != 0) { VScaleEq(D, Max_Bound); for (i = 0; i < cnt; i++) { intervals[i].bound /= Max_Bound; } } else { Max_Bound = 1; } /* Clear out the coefficients. */ coeffs[0] = coeffs[1] = coeffs[2] = coeffs[3] = 0.0; coeffs[4] = - Blob->Data->Threshold; /* * Step through the list of intersection points, adding the * influence of each component as it appears. */ fcoeffs = NULL; for (i = in_flag = 0; i < cnt; i++) { if ((intervals[i].type & 1) == ENTERING) { /* * Something is just starting to influence the ray, so calculate * its coefficients and add them into the pot. */ in_flag++; Element = intervals[i].Element; switch (Element->Type) { case BLOB_SPHERE: VSub(V1, P, Element->O); VDot(t0, V1, V1); VDot(t1, V1, D); t2 = Max_Bound * Max_Bound; c0 = Element->c[0]; c1 = Element->c[1]; c2 = Element->c[2]; fcoeffs = &(Element->f[0]); fcoeffs[0] = c0 * t2 * t2; fcoeffs[1] = 4.0 * c0 * t1 * t2; fcoeffs[2] = 2.0 * c0 * (2.0 * t1 * t1 + t0 * t2) + c1 * t2; fcoeffs[3] = 2.0 * t1 * (2.0 * c0 * t0 + c1); fcoeffs[4] = t0 * (c0 * t0 + c1) + c2; break; case BLOB_ELLIPSOID: MInvTransPoint(PP, P, Element->Trans); MInvTransDirection(DD, D, Element->Trans); VSub(V1, PP, Element->O); VDot(t0, V1, V1); VDot(t1, V1, DD); VDot(t2, DD, DD); c0 = Element->c[0]; c1 = Element->c[1]; c2 = Element->c[2]; fcoeffs = &(Element->f[0]); fcoeffs[0] = c0 * t2 * t2; fcoeffs[1] = 4.0 * c0 * t1 * t2; fcoeffs[2] = 2.0 * c0 * (2.0 * t1 * t1 + t0 * t2) + c1 * t2; fcoeffs[3] = 2.0 * t1 * (2.0 * c0 * t0 + c1); fcoeffs[4] = t0 * (c0 * t0 + c1) + c2; break; case BLOB_BASE_HEMISPHERE: case BLOB_APEX_HEMISPHERE: MInvTransPoint(PP, P, Element->Trans); MInvTransDirection(DD, D, Element->Trans); if (Element->Type == BLOB_APEX_HEMISPHERE) { PP[Z] -= Element->len; } VDot(t0, PP, PP); VDot(t1, PP, DD); VDot(t2, DD, DD); c0 = Element->c[0]; c1 = Element->c[1]; c2 = Element->c[2]; fcoeffs = &(Element->f[0]); fcoeffs[0] = c0 * t2 * t2; fcoeffs[1] = 4.0 * c0 * t1 * t2; fcoeffs[2] = 2.0 * c0 * (2.0 * t1 * t1 + t0 * t2) + c1 * t2; fcoeffs[3] = 2.0 * t1 * (2.0 * c0 * t0 + c1); fcoeffs[4] = t0 * (c0 * t0 + c1) + c2; break; case BLOB_CYLINDER: /* Transform ray into cylinder space. */ MInvTransPoint(PP, P, Element->Trans); MInvTransDirection(DD, D, Element->Trans); t0 = PP[X] * PP[X] + PP[Y] * PP[Y]; t1 = PP[X] * DD[X] + PP[Y] * DD[Y]; t2 = DD[X] * DD[X] + DD[Y] * DD[Y]; c0 = Element->c[0]; c1 = Element->c[1]; c2 = Element->c[2]; fcoeffs = &(Element->f[0]); fcoeffs[0] = c0 * t2 * t2; fcoeffs[1] = 4.0 * c0 * t1 * t2; fcoeffs[2] = 2.0 * c0 * (2.0 * t1 * t1 + t0 * t2) + c1 * t2; fcoeffs[3] = 2.0 * t1 * (2.0 * c0 * t0 + c1); fcoeffs[4] = t0 * (c0 * t0 + c1) + c2; break; default: Error("Unknown blob component in All_Blob_Intersections()."); } for (j = 0; j < 5; j++) { coeffs[j] += fcoeffs[j]; } } else { /* * We are losing the influence of a component --> * subtract off its coefficients. */ fcoeffs = intervals[i].Element->f; for (j = 0; j < 5; j++) { coeffs[j] -= fcoeffs[j]; } /* If no components are currently affecting the ray ---> skip ahead. */ if (--in_flag == 0) { continue; } } /* * If the following intersection lies close to the current intersection * then first add/subtract next region before testing. [DB 7/94] */ if ((i + 1 < cnt) && (fabs(intervals[i].bound - intervals[i + 1].bound) < EPSILON)) { continue; } /* * Transform polynomial in a way that the interval boundaries are moved * to 0 and 1, i. e. the roots of interest are between 0 and 1. [DB 10/94] */ l = intervals[i].bound; w = intervals[i+1].bound - l; newcoeffs[0] = coeffs[0] * w * w * w * w; newcoeffs[1] = (coeffs[1] + 4.0 * coeffs[0] * l) * w * w * w; newcoeffs[2] = (3.0 * l * (2.0 * coeffs[0] * l + coeffs[1]) + coeffs[2]) * w * w; newcoeffs[3] = (2.0 * l * (2.0 * l * (coeffs[0] * l + 0.75 * coeffs[1]) + coeffs[2]) + coeffs[3]) * w; newcoeffs[4] = l * (l * (l * (coeffs[0] * l + coeffs[1]) + coeffs[2]) + coeffs[3]) + coeffs[4]; /* Calculate coefficients of corresponding bezier curve. [DB 10/94] */ dk[0] = newcoeffs[4]; dk[1] = newcoeffs[4] + 0.25 * newcoeffs[3]; dk[2] = newcoeffs[4] + 0.50 * (newcoeffs[3] + newcoeffs[2] / 12.0); dk[3] = newcoeffs[4] + 0.50 * (0.375 * newcoeffs[3] + newcoeffs[2] + 0.125 * newcoeffs[1]); dk[4] = newcoeffs[4] + newcoeffs[3] + newcoeffs[2] + newcoeffs[1] + newcoeffs[0]; /* * Skip this interval if the ray doesn't intersect the convex hull of the * bezier curve, because no valid intersection will be found. [DB 10/94] */ if (((dk[0] >= 0.0) && (dk[1] >= 0.0) && (dk[2] >= 0.0) && (dk[3] >= 0.0) && (dk[4] >= 0.0)) || ((dk[0] <= 0.0) && (dk[1] <= 0.0) && (dk[2] <= 0.0) && (dk[3] <= 0.0) && (dk[4] <= 0.0))) { continue; } /* * Now we could do bezier clipping to find the roots * but I have no idea how this works. [DB 2/95] */ /* Solve polynomial. */ root_count = Solve_Polynomial(4, coeffs, roots, Test_Flag(Blob, STURM_FLAG), 1.0e-11); /* See if any of the roots are valid. */ for (j = 0; j < root_count; j++) { dist = roots[j]; /* * First see if the root is in the interval of influence of * the currently active components. */ if ((dist >= intervals[i].bound) && (dist <= intervals[i+1].bound)) { /* Correct distance. */ dist = (dist * Max_Bound + start_dist) / len; if ((dist > DEPTH_TOLERANCE) && (dist < Max_Distance)) { VEvaluateRay(IPoint, Ray->Initial, dist, Ray->Direction); if (Point_In_Clip(IPoint, Object->Clip)) { push_entry(dist, IPoint, Object, Depth_Stack); Intersection_Found = true; } } } } /* * If the blob isn't used inside a CSG and we have found at least * one intersection then we are ready, because all possible intersections * will be further away (we have a sorted list!). [DB 7/94] */ if (!(Blob->Type & IS_CHILD_OBJECT) && (Intersection_Found)) { break; } } if (Intersection_Found) { Increase_Counter(stats[Ray_Blob_Tests_Succeeded]); } return (Intersection_Found); } /************************************************************************ ***** * * FUNCTION * * insert_hit * * INPUT * * Blob - Pointer to blob structure * Element - Element to insert * t0, t1 - Intersection depths * * OUTPUT * * intervals - Pointer to sorted list of hits * cnt - Number of hits in intervals * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Store the points of intersection. Keep track of: whether this is * the start or end point of the hit, which component was pierced * by the ray, and the point along the ray that the hit occured at. * * CHANGES * * Oct 1994 : Modified to use memmove instead of loops for copying. [DB] * Sep 1995 : Changed to allow use of memcpy if memmove isn't available. [AED] * Jul 1996 : Changed to use POV_MEMMOVE, which can be memmove or pov_memmove. * Oct 1996 : Changed to avoid unnecessary compares. [DB] * ************************************************************************* *****/ static void insert_hit(BLOB_ELEMENT *Element, DBL t0, DBL t1, BLOB_INTERVAL *intervals, int *cnt) { int k; /* We are entering the component. */ intervals[*cnt].type = Element->Type | ENTERING; intervals[*cnt].bound = t0; intervals[*cnt].Element = Element; for (k = 0; t0 > intervals[k].bound; k++); if (k < *cnt) { /* * This hit point is smaller than one that already exists --> * bump the rest and insert it here. */ POV_MEMMOVE(&intervals[k+1], &intervals[k], (*cnt-k)*sizeof(BLOB_INTERVAL)); /* We are entering the component. */ intervals[k].type = Element->Type | ENTERING; intervals[k].bound = t0; intervals[k].Element = Element; (*cnt)++; /* We are exiting the component. */ intervals[*cnt].type = Element->Type | EXITING; intervals[*cnt].bound = t1; intervals[*cnt].Element = Element; for (k = k + 1; t1 > intervals[k].bound; k++); if (k < *cnt) { POV_MEMMOVE(&intervals[k+1], &intervals[k], (*cnt-k)*sizeof(BLOB_INTERVAL)); /* We are exiting the component. */ intervals[k].type = Element->Type | EXITING; intervals[k].bound = t1; intervals[k].Element = Element; } (*cnt)++; } else { /* Just plop the start and end points at the end of the list */ (*cnt)++; /* We are exiting the component. */ intervals[*cnt].type = Element->Type | EXITING; intervals[*cnt].bound = t1; intervals[*cnt].Element = Element; (*cnt)++; } } /************************************************************************ ***** * * FUNCTION * * intersect_cylinder * * INPUT * * Element - Pointer to element structure * P, D - Ray = P + t * D * mindist - Min. valid distance * * OUTPUT * * tmin, tmax - Intersection depths found * * RETURNS * * AUTHOR * * Dieter Bayer (with help from Alexander Enzmann) * * DESCRIPTION * * - * * CHANGES * * Jul 1994 : Creation. * ************************************************************************* *****/ static int intersect_cylinder(BLOB_ELEMENT *Element, VECTOR P, VECTOR D, DBL mindist, DBL *tmin, DBL *tmax) { DBL a, b, c, d, t, u, v, w, len; VECTOR PP, DD; /* Transform ray into cylinder space. */ MInvTransPoint(PP, P, Element->Trans); MInvTransDirection(DD, D, Element->Trans); VLength(len, DD); VInverseScaleEq(DD, len); /* Intersect ray with cylinder. */ a = DD[X] * DD[X] + DD[Y] * DD[Y]; if (a > EPSILON) { b = PP[X] * DD[X] + PP[Y] * DD[Y]; c = PP[X] * PP[X] + PP[Y] * PP[Y] - Element->rad2; d = b * b - a * c; if (d > EPSILON) { d = sqrt(d); t = ( - b + d) / a; w = PP[Z] + t * DD[Z]; if ((w >= 0.0) && (w <= Element->len)) { if (t < *tmin) { *tmin = t; } if (t > *tmax) { *tmax = t; } } t = ( - b - d) / a; w = PP[Z] + t * DD[Z]; if ((w >= 0.0) && (w <= Element->len)) { if (t < *tmin) { *tmin = t; } if (t > *tmax) { *tmax = t; } } } } /* Intersect base/cap plane. */ if (fabs(DD[Z]) > EPSILON) { /* Intersect base plane. */ t = - PP[Z] / DD[Z]; u = PP[X] + t * DD[X]; v = PP[Y] + t * DD[Y]; if ((u * u + v * v) <= Element->rad2) { if (t < *tmin) { *tmin = t; } if (t > *tmax) { *tmax = t; } } /* Intersect cap plane. */ t = (Element->len - PP[Z]) / DD[Z]; u = PP[X] + t * DD[X]; v = PP[Y] + t * DD[Y]; if ((u * u + v * v) <= Element->rad2) { if (t < *tmin) { *tmin = t; } if (t > *tmax) { *tmax = t; } } } /* Check if the intersections are valid. */ *tmin /= len; *tmax /= len; if (*tmin < mindist) { *tmin = 0.0; } if (*tmax < mindist) { *tmax = 0.0; } if (*tmin >= *tmax) { return (false); } return (true); } /************************************************************************ ***** * * FUNCTION * * intersect_ellipsoid * * INPUT * * Element - Pointer to element structure * P, D - Ray = P + t * D * mindist - Min. valid distance * * OUTPUT * * tmin, tmax - Intersection depths found * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * - * * CHANGES * * Sep 1994 : Creation. * ************************************************************************* *****/ static int intersect_ellipsoid(BLOB_ELEMENT *Element, VECTOR P, VECTOR D, DBL mindist, DBL *tmin, DBL *tmax) { DBL b, d, t, len; VECTOR V1, PP, DD; MInvTransPoint(PP, P, Element->Trans); MInvTransDirection(DD, D, Element->Trans); VLength(len, DD); VInverseScaleEq(DD, len); VSub(V1, PP, Element->O); VDot(b, V1, DD); VDot(t, V1, V1); d = b * b - t + Element->rad2; if (d < EPSILON) { return (false); } d = sqrt(d); *tmax = ( - b + d) / len; if (*tmax < mindist) { *tmax = 0.0; } *tmin = ( - b - d) / len; if (*tmin < mindist) { *tmin = 0.0; } if (*tmax == *tmin) { return (false); } else { if (*tmax < *tmin) { d = *tmin; *tmin = *tmax; *tmax = d; } } return (true); } /************************************************************************ ***** * * FUNCTION * * intersect_hemisphere * * INPUT * * Element - Pointer to element structure * P, D - Ray = P + t * D * mindist - Min. valid distance * * OUTPUT * * tmin, tmax - Intersection depths found * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * - * * CHANGES * * Jul 1994 : Creation (with help from Alexander Enzmann). * ************************************************************************* *****/ static int intersect_hemisphere(BLOB_ELEMENT *Element, VECTOR P, VECTOR D, DBL mindist, DBL *tmin, DBL *tmax) { DBL b, d, t, z1, z2, len; VECTOR PP, DD; /* Transform ray into hemisphere space. */ MInvTransPoint(PP, P, Element->Trans); MInvTransDirection(DD, D, Element->Trans); VLength(len, DD); VInverseScaleEq(DD, len); if (Element->Type == BLOB_BASE_HEMISPHERE) { VDot(b, PP, DD); VDot(t, PP, PP); d = b * b - t + Element->rad2; if (d < EPSILON) { return (false); } d = sqrt(d); *tmax = - b + d; *tmin = - b - d; if (*tmax < *tmin) { d = *tmin; *tmin = *tmax; *tmax = d; } /* Cut intersection at the plane. */ z1 = PP[Z] + *tmin * DD[Z]; z2 = PP[Z] + *tmax * DD[Z]; /* If both points are inside --> no intersection */ if ((z1 >= 0.0) && (z2 >= 0.0)) { return (false); } /* If both points are outside --> intersections found */ if ((z1 < 0.0) && (z2 < 0.0)) { *tmin /= len; *tmax /= len; return (true); } /* Determine intersection with plane. */ t = - PP[Z] / DD[Z]; if (z1 >= 0.0) { /* Ray is crossing the plane from inside to outside. */ *tmin = (t < mindist) ? 0.0 : t; } else { /* Ray is crossing the plane from outside to inside. */ *tmax = (t < mindist) ? 0.0 : t; } *tmin /= len; *tmax /= len; return (true); } else { PP[Z] -= Element->len; VDot(b, PP, DD); VDot(t, PP, PP); d = b * b - t + Element->rad2; if (d < EPSILON) { return (false); } d = sqrt(d); *tmax = - b + d; *tmin = - b - d; if (*tmax < *tmin) { d = *tmin; *tmin = *tmax; *tmax = d; } /* Cut intersection at the plane. */ z1 = PP[Z] + *tmin * DD[Z]; z2 = PP[Z] + *tmax * DD[Z]; /* If both points are inside --> no intersection */ if ((z1 <= 0.0) && (z2 <= 0.0)) { return (false); } /* If both points are outside --> intersections found */ if ((z1 > 0.0) && (z2 > 0.0)) { *tmin /= len; *tmax /= len; return (true); } /* Determine intersection with plane. */ t = - PP[Z] / DD[Z]; if (z1 <= 0.0) { /* Ray is crossing the plane from inside to outside. */ *tmin = (t < mindist) ? 0.0 : t; } else { /* Ray is crossing the plane from outside to inside. */ *tmax = (t < mindist) ? 0.0 : t; } *tmin /= len; *tmax /= len; return (true); } } /************************************************************************ ***** * * FUNCTION * * intersect_sphere * * INPUT * * Element - Pointer to element structure * P, D - Ray = P + t * D * mindist - Min. valid distance * * OUTPUT * * tmin, tmax - Intersection depths found * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * - * * CHANGES * * Jul 1994 : Creation (with help from Alexander Enzmann). * ************************************************************************* *****/ static int intersect_sphere(BLOB_ELEMENT *Element, VECTOR P, VECTOR D, DBL mindist, DBL *tmin, DBL *tmax) { DBL b, d, t; VECTOR V1; VSub(V1, P, Element->O); VDot(b, V1, D); VDot(t, V1, V1); d = b * b - t + Element->rad2; if (d < EPSILON) { return (false); } d = sqrt(d); *tmax = - b + d; if (*tmax < mindist) { *tmax = 0.0; } *tmin = - b - d; if (*tmin < mindist) { *tmin = 0.0; } if (*tmax == *tmin) { return (false); } else { if (*tmax < *tmin) { d = *tmin; *tmin = *tmax; *tmax = d; } } return (true); } /************************************************************************ ***** * * FUNCTION * * intersect_element * * INPUT * * P, D - Ray = P + t * D * Element - Pointer to element structure * mindist - Min. valid distance * * OUTPUT * * tmin, tmax - Intersection depths found * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * - * * CHANGES * * Jul 1994 : Creation. * ************************************************************************* *****/ static int intersect_element(VECTOR P, VECTOR D, BLOB_ELEMENT *Element, DBL mindist, DBL *tmin, DBL *tmax) { #ifdef BLOB_EXTRA_STATS Increase_Counter(stats[Blob_Element_Tests]); #endif *tmin = BOUND_HUGE; *tmax = - BOUND_HUGE; switch (Element->Type) { case BLOB_SPHERE: if (!intersect_sphere(Element, P, D, mindist, tmin, tmax)) { return (false); } break; case BLOB_ELLIPSOID: if (!intersect_ellipsoid(Element, P, D, mindist, tmin, tmax)) { return (false); } break; case BLOB_BASE_HEMISPHERE: case BLOB_APEX_HEMISPHERE: if (!intersect_hemisphere(Element, P, D, mindist, tmin, tmax)) { return (false); } break; case BLOB_CYLINDER: if (!intersect_cylinder(Element, P, D, mindist, tmin, tmax)) { return (false); } break; } #ifdef BLOB_EXTRA_STATS Increase_Counter(stats[Blob_Element_Tests_Succeeded]); #endif return (true); } /************************************************************************ ***** * * FUNCTION * * determine_influences * * INPUT * * P, D - Ray = P + t * D * Blob - Pointer to blob structure * mindist - Min. valid distance * * OUTPUT * * intervals - Sorted list of intersections found * * RETURNS * * int - Number of intersection found * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * Make a sorted list of points along the ray at which the various blob * components start and stop adding their influence. * * CHANGES * * Jul 1994 : Added code for bounding hierarchy traversal. [DB] * ************************************************************************* *****/ static int determine_influences(VECTOR P, VECTOR D, BLOB *Blob, DBL mindist, BLOB_INTERVAL *intervals) { int i; int cnt, size; DBL b, t, t0, t1; VECTOR V1; BSPHERE_TREE *Tree; cnt = 0; if (Blob->Data->Tree == NULL) { /* There's no bounding hierarchy so just step through all elements. */ for (i = 0; i < Blob->Data->Number_Of_Components; i++) { if (intersect_element(P, D, &Blob->Data->Entry[i], mindist, &t0, &t1)) { insert_hit(&Blob->Data->Entry[i], t0, t1, intervals, &cnt); } } } else { /* Use blob's bounding hierarchy. */ size = 0; Queue[size++] = Blob->Data->Tree; while (size > 0) { Tree = Queue[--size]; /* Test if current node is a leaf. */ if (Tree->Entries <= 0) { /* Test element. */ if (intersect_element(P, D, (BLOB_ELEMENT *)Tree->Node, mindist, &t0, &t1)) { insert_hit((BLOB_ELEMENT *)Tree->Node, t0, t1, intervals, &cnt); } } else { /* Test all sub-nodes. */ for (i = 0; i < (int)Tree->Entries; i++) { #ifdef BLOB_EXTRA_STATS Increase_Counter(stats[Blob_Bound_Tests]); #endif VSub(V1, Tree->Node[i]->C, P); VDot(b, V1, D); VDot(t, V1, V1); if ((t - Sqr(b)) <= Tree->Node[i]->r2) { #ifdef BLOB_EXTRA_STATS Increase_Counter(stats[Blob_Bound_Tests_Succeeded]); #endif insert_node(Tree->Node[i], &size); } } } } } return (cnt); } /************************************************************************ ***** * * FUNCTION * * calculate_element_field * * INPUT * * Element - Pointer to element structure * P - Point whos field value is calculated * * OUTPUT * * RETURNS * * DBL - Field value * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Calculate the field value of a single element in a given point P * (which must already have been transformed into blob space). * * CHANGES * * Jul 1994 : Added code for cylindrical and ellipsoidical blobs. [DB] * ************************************************************************* *****/ static DBL calculate_element_field(BLOB_ELEMENT *Element, VECTOR P) { DBL rad2, density; VECTOR V1, PP; density = 0.0; switch (Element->Type) { case BLOB_SPHERE: VSub(V1, P, Element->O); VDot(rad2, V1, V1); if (rad2 < Element->rad2) { density = rad2 * (rad2 * Element->c[0] + Element->c[1]) + Element->c[2]; } break; case BLOB_ELLIPSOID: MInvTransPoint(PP, P, Element->Trans); VSub(V1, PP, Element->O); VDot(rad2, V1, V1); if (rad2 < Element->rad2) { density = rad2 * (rad2 * Element->c[0] + Element->c[1]) + Element->c[2]; } break; case BLOB_BASE_HEMISPHERE: MInvTransPoint(PP, P, Element->Trans); if (PP[Z] <= 0.0) { VDot(rad2, PP, PP); if (rad2 <= Element->rad2) { density = rad2 * (rad2 * Element->c[0] + Element->c[1]) + Element->c[2]; } } break; case BLOB_APEX_HEMISPHERE: MInvTransPoint(PP, P, Element->Trans); PP[Z] -= Element->len; if (PP[Z] >= 0.0) { VDot(rad2, PP, PP); if (rad2 <= Element->rad2) { density = rad2 * (rad2 * Element->c[0] + Element->c[1]) + Element->c[2]; } } break; case BLOB_CYLINDER: MInvTransPoint(PP, P, Element->Trans); if ((PP[Z] >= 0.0) && (PP[Z] <= Element->len)) { if ((rad2 = Sqr(PP[X]) + Sqr(PP[Y])) <= Element->rad2) { density = rad2 * (rad2 * Element->c[0] + Element->c[1]) + Element->c[2]; } } break; } return (density); } /************************************************************************ ***** * * FUNCTION * * calculate_field_value * * INPUT * * Blob - Pointer to blob structure * P - Point whos field value is calculated * * OUTPUT * * RETURNS * * DBL - Field value * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * Calculate the field value of a blob in a given point P * (which must already have been transformed into blob space). * * CHANGES * * Jul 1994 : Added code for bounding hierarchy traversal. [DB] * ************************************************************************* *****/ static DBL calculate_field_value(BLOB *Blob, VECTOR P) { int i; int size; DBL density, rad2; VECTOR V1; BSPHERE_TREE *Tree; density = 0.0; if (Blob->Data->Tree == NULL) { /* There's no tree --> step through all elements. */ for (i = 0; i < Blob->Data->Number_Of_Components; i++) { density += calculate_element_field(&Blob->Data->Entry[i], P); } } else { /* A tree exists --> step through the tree. */ size = 0; Queue[size++] = Blob->Data->Tree; while (size > 0) { Tree = Queue[--size]; /* Test if current node is a leaf. */ if (Tree->Entries <= 0) { density += calculate_element_field((BLOB_ELEMENT *)Tree->Node, P); } else { /* Test all sub-nodes. */ for (i = 0; i < (int)Tree->Entries; i++) { /* Insert sub-node if we are inside. */ VSub(V1, P, Tree->Node[i]->C); VDot(rad2, V1, V1); if (rad2 <= Tree->Node[i]->r2) { insert_node(Tree->Node[i], &size); } } } } } return (density); } /************************************************************************ ***** * * FUNCTION * * Inside_Blob * * INPUT * * Test_Point - Point to test * Object - Pointer to blob structure * * OUTPUT * * RETURNS * * int - true if Test_Point is inside * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Calculate the density at the given point and then compare to * the threshold to see if we are in or out of the blob. * * CHANGES * * - * ************************************************************************* *****/ static int Inside_Blob(VECTOR Test_Point, OBJECT *Object) { VECTOR New_Point; BLOB *Blob = (BLOB *) Object; /* Transform the point into blob space. */ if (Blob->Trans != NULL) { MInvTransPoint(New_Point, Test_Point, Blob->Trans); } else { Assign_Vector(New_Point, Test_Point); } if (calculate_field_value(Blob, New_Point) > Blob->Data->Threshold - INSIDE_TOLERANCE) { /* We are inside. */ return (!Test_Flag(Blob, INVERTED_FLAG)); } else { /* We are outside. */ return (Test_Flag(Blob, INVERTED_FLAG)); } } /************************************************************************ ***** * * FUNCTION * * element_normal * * INPUT * * P - Surface point * Element - Pointer to element structure * * OUTPUT * * Result - Element's normal * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * Calculate the normal of a single element in the point P. * * CHANGES * * Jul 1994 : Creation (with help from Alexander Enzmann). * ************************************************************************* *****/ static void element_normal(VECTOR Result, VECTOR P, BLOB_ELEMENT *Element) { DBL val, dist; VECTOR V1, PP; switch (Element->Type) { case BLOB_SPHERE: VSub(V1, P, Element->O); VDot(dist, V1, V1); if (dist <= Element->rad2) { val = -2.0 * Element->c[0] * dist - Element->c[1]; VAddScaledEq(Result, val, V1); } break; case BLOB_ELLIPSOID: MInvTransPoint(PP, P, Element->Trans); VSub(V1, PP, Element->O); VDot(dist, V1, V1); if (dist <= Element->rad2) { val = -2.0 * Element->c[0] * dist - Element->c[1]; MTransNormal(V1, V1, Element->Trans); VAddScaledEq(Result, val, V1); } break; case BLOB_BASE_HEMISPHERE: MInvTransPoint(PP, P, Element->Trans); if (PP[Z] <= 0.0) { VDot(dist, PP, PP); if (dist <= Element->rad2) { val = -2.0 * Element->c[0] * dist - Element->c[1]; MTransNormal(PP, PP, Element->Trans); VAddScaledEq(Result, val, PP); } } break; case BLOB_APEX_HEMISPHERE: MInvTransPoint(PP, P, Element->Trans); PP[Z] -= Element->len; if (PP[Z] >= 0.0) { VDot(dist, PP, PP); if (dist <= Element->rad2) { val = -2.0 * Element->c[0] * dist - Element->c[1]; MTransNormal(PP, PP, Element->Trans); VAddScaledEq(Result, val, PP); } } break; case BLOB_CYLINDER: MInvTransPoint(PP, P, Element->Trans); if ((PP[Z] >= 0.0) && (PP[Z] <= Element->len)) { if ((dist = Sqr(PP[X]) + Sqr(PP[Y])) <= Element->rad2) { val = -2.0 * Element->c[0] * dist - Element->c[1]; PP[Z] = 0.0; MTransNormal(PP, PP, Element->Trans); VAddScaledEq(Result, val, PP); } } break; } } /************************************************************************ ***** * * FUNCTION * * Blob_Normal * * INPUT * * Object - Pointer to blob structure * Inter - Pointer to intersection * * OUTPUT * * Result - Blob's normal * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Calculate the blob's surface normal in the intersection point. * * CHANGES * * Jul 1994 : Added code for bounding hierarchy traversal. [DB] * ************************************************************************* *****/ static void Blob_Normal(VECTOR Result, OBJECT *Object, INTERSECTION *Inter) { int i; int size; DBL dist, val; VECTOR New_Point, V1; BLOB *Blob = (BLOB *) Object; BSPHERE_TREE *Tree; /* Transform the point into the blob space. */ if (Blob->Trans != NULL) { MInvTransPoint(New_Point, Inter->IPoint, Blob->Trans); } else { Assign_Vector(New_Point, Inter->IPoint); } Make_Vector(Result, 0.0, 0.0, 0.0); /* For each component that contributes to this point, add its bit to the normal */ if (Blob->Data->Tree == NULL) { /* There's no tree --> step through all elements. */ for (i = 0; i < Blob->Data->Number_Of_Components; i++) { element_normal(Result, New_Point, &(Blob->Data->Entry[i])); } } else { /* A tree exists --> step through the tree. */ size = 0; Queue[size++] = Blob->Data->Tree; while (size > 0) { Tree = Queue[--size]; /* Test if current node is a leaf. */ if (Tree->Entries <= 0) { element_normal(Result, New_Point, (BLOB_ELEMENT *)Tree->Node); } else { /* Test all sub-nodes. */ for (i = 0; i < (int)Tree->Entries; i++) { /* Insert sub-node if we are inside. */ VSub(V1, New_Point, Tree->Node[i]->C); VDot(dist, V1, V1); if (dist <= Tree->Node[i]->r2) { insert_node(Tree->Node[i], &size); } } } } } VDot(val, Result, Result); if (val == 0.0) { Make_Vector(Result, 1.0, 0.0, 0.0); } else { /* Normalize normal vector. */ val = 1.0 / sqrt(val); VScaleEq(Result, val); } /* Transform back to world space. */ if (Blob->Trans != NULL) { MTransNormal(Result, Result, Blob->Trans); VNormalize(Result, Result); } } /************************************************************************ ***** * * FUNCTION * * Translate_Blob * * INPUT * * Vector - Translation vector * * OUTPUT * * Object - Pointer to blob structure * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Translate a blob. * * CHANGES * * - * ************************************************************************* *****/ static void Translate_Blob(OBJECT *Object, VECTOR, TRANSFORM *Trans) { Transform_Blob(Object, Trans); } /************************************************************************ ***** * * FUNCTION * * Rotate_Blob * * INPUT * * Vector - Rotation vector * * OUTPUT * * Object - Pointer to blob structure * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Rotate a blob. * * CHANGES * * - * ************************************************************************* *****/ static void Rotate_Blob(OBJECT *Object, VECTOR, TRANSFORM *Trans) { Transform_Blob(Object, Trans); } /************************************************************************ ***** * * FUNCTION * * Scale_Blob * * INPUT * * Vector - Scaling vector * * OUTPUT * * Object - Pointer to blob structure * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Scale a blob. * * CHANGES * * - * ************************************************************************* *****/ static void Scale_Blob(OBJECT *Object, VECTOR, TRANSFORM *Trans) { Transform_Blob(Object, Trans); } /************************************************************************ ***** * * FUNCTION * * Transform_Blob * * INPUT * * Trans - Pointer to transformation * * OUTPUT * * Object - Pointer to blob structure * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Transform a blob. * * CHANGES * * - * ************************************************************************* *****/ static void Transform_Blob(OBJECT *Object, TRANSFORM *Trans) { int i; BLOB *Blob = (BLOB *)Object; if (Blob->Trans == NULL) { Blob->Trans = Create_Transform(); } Recompute_BBox(&Object->BBox, Trans); Compose_Transforms(Blob->Trans, Trans); for (i = 0; i < Blob->Data->Number_Of_Components; i++) { Transform_Textures(Blob->Element_Texture[i], Trans); } } /************************************************************************ ***** * * FUNCTION * * Invert_Blob * * INPUT * * Object - Pointer to blob structure * * OUTPUT * * Object * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Invert a blob. * * CHANGES * * - * ************************************************************************* *****/ static void Invert_Blob(OBJECT *Object) { Invert_Flag(Object, INVERTED_FLAG); } /************************************************************************ ***** * * FUNCTION * * Create_Blob * * INPUT * * Object - Pointer to blob structure * * OUTPUT * * Object * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Create a new blob. * * CHANGES * * - * ************************************************************************* *****/ BLOB *Create_Blob() { BLOB *New; New = (BLOB *)POV_MALLOC(sizeof (BLOB), "blob"); INIT_OBJECT_FIELDS(New, BLOB_OBJECT, &Blob_Methods) Set_Flag(New, HIERARCHY_FLAG); New->Trans = NULL; New->Data = NULL; New->Element_Texture = NULL; return (New); } /************************************************************************ ***** * * FUNCTION * * Copy_Blob * * INPUT * * Object - Pointer to blob structure * * OUTPUT * * Object * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Copy a blob. * * NOTE: The components are not copied, only the number of references is * counted, so that Destroy_Blob() knows if they can be destroyed. * * CHANGES * * Jul 1994 : Added code for blob data reference counting. [DB] * ************************************************************************* *****/ static BLOB *Copy_Blob(OBJECT *Object) { int i; BLOB *New, *Old = (BLOB *)Object; New = Create_Blob(); /* Copy blob. */ *New = *Old; New->Trans = Copy_Transform(New->Trans); New->Data->References++; New->Element_Texture = (TEXTURE **)POV_MALLOC(New->Data->Number_Of_Components*sizeof(TEXTURE *), "blob texture list"); for (i = 0; i < New->Data->Number_Of_Components; i++) { New->Element_Texture[i] = Copy_Textures(Old->Element_Texture[i]); } return (New); } /************************************************************************ ***** * * FUNCTION * * Create_Blob_List_Element * * INPUT * * OUTPUT * * RETURNS * * BLOB_LIST * - Pointer to blob element * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * Create a new blob element in the component list used during parsing. * * CHANGES * * Sep 1994 : Creation. * ************************************************************************* *****/ BLOB_LIST *Create_Blob_List_Element() { BLOB_LIST *New; New = (BLOB_LIST *)POV_MALLOC(sizeof(BLOB_LIST), "blob component"); init_blob_element(&New->elem); return (New); } /************************************************************************ ***** * * FUNCTION * * Destroy_Blob * * INPUT * * Object - Pointer to blob structure * * OUTPUT * * Object * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Destroy a blob. * * NOTE: The blob data is destroyed if they are no longer used by any copy. * * CHANGES * * Jul 1994 : Added code for blob data reference counting. [DB] * * Dec 1994 : Fixed memory leakage. [DB] * * Aug 1995 : Fixed freeing of already freed memory. [DB] * ************************************************************************* *****/ static void Destroy_Blob(OBJECT *Object) { int i; BLOB *Blob = (BLOB *)Object; Destroy_Transform(Blob->Trans); for (i = 0; i < Blob->Data->Number_Of_Components; i++) { Destroy_Textures(Blob->Element_Texture[i]); } POV_FREE(Blob->Element_Texture); if (--(Blob->Data->References) == 0) { Destroy_Bounding_Sphere_Hierarchy(Blob->Data->Tree); for (i = 0; i < Blob->Data->Number_Of_Components; i++) { /* * Make sure to destroy multiple references of a texture * and/or transformation only once. Multiple references * are only used with cylindrical blobs. Thus it's * enough to ignore all cylinder caps. */ if ((Blob->Data->Entry[i].Type == BLOB_SPHERE) || (Blob->Data->Entry[i].Type == BLOB_ELLIPSOID) || (Blob->Data->Entry[i].Type == BLOB_CYLINDER)) { Destroy_Transform(Blob->Data->Entry[i].Trans); Blob->Data->Entry[i].Trans = NULL; Blob->Data->Entry[i].Texture = NULL; } } POV_FREE(Blob->Data->Entry); POV_FREE(Blob->Data->Intervals); POV_FREE(Blob->Data); } POV_FREE(Object); } /************************************************************************ ***** * * FUNCTION * * Compute_Blob_BBox * * INPUT * * Blob - Blob * * OUTPUT * * Blob * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * Calculate the bounding box of a blob. * * CHANGES * * Aug 1994 : Creation. * ************************************************************************* *****/ static void Compute_Blob_BBox(BLOB *Blob) { int i; DBL radius, radius2; VECTOR Center, Min, Max; Make_Vector(Min, BOUND_HUGE, BOUND_HUGE, BOUND_HUGE); Make_Vector(Max, - BOUND_HUGE, - BOUND_HUGE, - BOUND_HUGE); for (i = 0; i < Blob->Data->Number_Of_Components; i++) { if (Blob->Data->Entry[i].c[2] > 0.0) { get_element_bounding_sphere(&Blob->Data->Entry[i], Center, &radius2); radius = sqrt(radius2); Min[X] = min(Min[X], Center[X] - radius); Min[Y] = min(Min[Y], Center[Y] - radius); Min[Z] = min(Min[Z], Center[Z] - radius); Max[X] = max(Max[X], Center[X] + radius); Max[Y] = max(Max[Y], Center[Y] + radius); Max[Z] = max(Max[Z], Center[Z] + radius); } } Make_BBox_from_min_max(Blob->BBox, Min, Max); if (Blob->Trans != NULL) { Recompute_BBox(&Blob->BBox, Blob->Trans); } } /************************************************************************ ***** * * FUNCTION * * get_element_bounding_sphere * * INPUT * * Element - Pointer to element * Center - Bounding sphere's center * Radius2 - Bounding sphere's squared radius * * OUTPUT * * Center, Radius2 * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * Calculate the bounding sphere of a blob element. * * CHANGES * * Sep 1994 : Creation. * ************************************************************************* *****/ static void get_element_bounding_sphere(BLOB_ELEMENT *Element, VECTOR Center, DBL *Radius2) { DBL r, r2 = 0.0; VECTOR C; BBOX local_BBox; switch (Element->Type) { case BLOB_SPHERE: case BLOB_ELLIPSOID: r2 = Element->rad2; Assign_Vector(C, Element->O); break; case BLOB_BASE_HEMISPHERE: r2 = Element->rad2; Make_Vector(C, 0.0, 0.0, 0.0); break; case BLOB_APEX_HEMISPHERE: r2 = Element->rad2; Make_Vector(C, 0.0, 0.0, Element->len); break; case BLOB_CYLINDER : Make_Vector(C, 0.0, 0.0, 0.5 * Element->len); r2 = Element->rad2 + Sqr(0.5 * Element->len); break; } /* Transform bounding sphere if necessary. */ if (Element->Trans != NULL) { r = sqrt(r2); MTransPoint(C, C, Element->Trans); Make_BBox(local_BBox, 0, 0, 0, r, r, r); Recompute_BBox(&local_BBox, Element->Trans); r = max(max(fabs(local_BBox.Lengths[X]), fabs(local_BBox.Lengths[Y])), fabs(local_BBox.Lengths[Z])); r2 = Sqr(r) + EPSILON; } Assign_Vector(Center, C); *Radius2 = r2; } /************************************************************************ ***** * * FUNCTION * * init_blob_element * * INPUT * * Element - Pointer to blob element * * OUTPUT * * Element * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * Init blob element. * * CHANGES * * Sep 1994 : Creation. * ************************************************************************* *****/ static void init_blob_element(BLOB_ELEMENT *Element) { Element->Type = 0; Element->index = 0; Element->len = Element->rad2 = 0.0; Element->c[0] = Element->c[1] = Element->c[2] = Element->f[0] = Element->f[1] = Element->f[2] = Element->f[3] = Element->f[4] = 0.0; Make_Vector(Element->O, 0.0, 0.0, 0.0); Element->Texture = NULL; Element->Trans = NULL; } /************************************************************************ ***** * * FUNCTION * * Make_Blob * * INPUT * * Blob - Pointer to blob structure * threshold - Blob's threshold * BlobList - Pointer to elements * npoints - Number of elements * * OUTPUT * * Blob * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Create a blob after it was read from the scene file. * * Starting with the density function: (1-r^2)^2, we have a field * that varies in strength from 1 at r = 0 to 0 at r = 1. By * substituting r/rad for r, we can adjust the range of influence * of a particular component. By multiplication by coeff, we can * adjust the amount of total contribution, giving the formula: * * coeff * (1 - (r/rad)^2)^2 * * This varies in strength from coeff at r = 0, to 0 at r = rad. * * CHANGES * * Jul 1994 : Added code for cylindrical and ellipsoidical blobs. [DB] * ************************************************************************* *****/ void Make_Blob(BLOB *Blob, DBL threshold, BLOB_LIST *BlobList, int npoints) { int i, count; DBL rad2, coeff; BLOB_LIST *temp; BLOB_ELEMENT *Entry; if (npoints < 1) { Error("Need at least one component in a blob."); } /* Figure out how many components there will be. */ temp = BlobList; for (i = count = 0; i < npoints; i++) { if (temp->elem.Type & BLOB_CYLINDER) { count += 3; } else { count++; } temp = temp->next; /* Test for too many components. [DB 12/94] */ if (count >= MAX_BLOB_COMPONENTS) { Error("There are more than %d components in a blob.", MAX_BLOB_COMPONENTS); } } /* Initialize the blob data. */ Blob->Data->Threshold = threshold; Entry = Blob->Data->Entry; for (i = 0; i < npoints; i++) { temp = BlobList; if ((fabs(temp->elem.c[2]) < EPSILON) || (temp->elem.rad2 < EPSILON)) { Warning(0, "Degenerate Blob element"); } /* Initialize component. */ *Entry = temp->elem; /* We have a multi-texture blob. */ if (Entry->Texture != NULL) { Set_Flag(Blob, MULTITEXTURE_FLAG); } /* Store blob specific information. */ switch (temp->elem.Type) { case BLOB_ELLIPSOID : case BLOB_SPHERE : rad2 = temp->elem.rad2; coeff = temp->elem.c[2]; Entry->c[0] = coeff / (rad2 * rad2); Entry->c[1] = -(2.0 * coeff) / rad2; Entry->c[2] = coeff; Entry++; break; case BLOB_CYLINDER : rad2 = temp->elem.rad2; coeff = temp->elem.c[2]; /* Create cylindrical component. */ Entry->c[0] = coeff / (rad2 * rad2); Entry->c[1] = -(2.0 * coeff) / rad2; Entry->c[2] = coeff; Entry++; /* Create hemispherical component at the base. */ *Entry = temp->elem; Entry->Type = BLOB_BASE_HEMISPHERE; Entry->c[0] = coeff / (rad2 * rad2); Entry->c[1] = -(2.0 * coeff) / rad2; Entry->c[2] = coeff; Entry++; /* Create hemispherical component at the apex. */ *Entry = temp->elem; Entry->Type = BLOB_APEX_HEMISPHERE; Entry->c[0] = coeff / (rad2 * rad2); Entry->c[1] = -(2.0 * coeff) / rad2; Entry->c[2] = coeff; Entry++; break; default : Error("Unknown blob component."); } /* Get rid of texture non longer needed. */ BlobList = BlobList->next; Destroy_Textures(temp->elem.Texture); POV_FREE(temp); } for (i = 0; i < count; i++) { Blob->Data->Entry[i].index = i; } /* Compute bounding box. */ Compute_Blob_BBox(Blob); /* Create bounding sphere hierarchy. */ if (Test_Flag(Blob, HIERARCHY_FLAG)) { build_bounding_hierarchy(Blob); } } /************************************************************************ ***** * * FUNCTION * * Test_Blob_Opacity * * INPUT * * Blob - Pointer to blob structure * * OUTPUT * * Blob * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * Set the opacity flag of the blob according to the opacity * of the blob's texture(s). * * CHANGES * * Apr 1996 : Creation. * ************************************************************************* *****/ void Test_Blob_Opacity(BLOB *Blob) { int i; /* Initialize opacity flag to the opacity of the object's texture. */ if ((Blob->Texture == NULL) || (Test_Opacity(Blob->Texture))) { Set_Flag(Blob, OPAQUE_FLAG); } if (Test_Flag(Blob, MULTITEXTURE_FLAG)) { for (i = 0; i < Blob->Data->Number_Of_Components; i++) { if (Blob->Element_Texture[i] != NULL) { /* If component's texture isn't opaque the blob is neither. */ if (!Test_Opacity(Blob->Element_Texture[i])) { Clear_Flag(Blob, OPAQUE_FLAG); } } } } } /************************************************************************ ***** * * FUNCTION * * build_bounding_hierarchy * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * Create the bounding sphere hierarchy. * * CHANGES * * Oct 1994 : Creation. (Derived from the bounding slab creation code) * ************************************************************************* *****/ static void build_bounding_hierarchy(BLOB *Blob) { int i, nElem, maxelements; BSPHERE_TREE **Elements; nElem = (int)Blob->Data->Number_Of_Components; maxelements = 2 * nElem; /* * Now allocate an array to hold references to these elements. */ Elements = (BSPHERE_TREE **)POV_MALLOC(maxelements*sizeof(BSPHERE_TREE *), "blob bounding hierarchy"); /* Init list with blob elements. */ for (i = 0; i < nElem; i++) { Elements[i] = (BSPHERE_TREE *)POV_MALLOC(sizeof(BSPHERE_TREE), "blob bounding hierarchy"); Elements[i]->Entries = 0; Elements[i]->Node = (BSPHERE_TREE **)&Blob->Data->Entry[i]; get_element_bounding_sphere(&Blob->Data->Entry[i], Elements[i]->C, &Elements[i]->r2); } Build_Bounding_Sphere_Hierarchy(&Blob->Data->Tree, nElem, Elements); /* Get rid of the Elements array. */ POV_FREE(Elements); } /************************************************************************ ***** * * FUNCTION * * Determine_Blob_Textures * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * Determine the textures and weights of all components affecting * the given intersection point. The weights are calculated from * the field values and sum to 1. * * CHANGES * * Sep 1994 : Creation. * * Mar 1996 : Make the call to resize the textures/weights list just once * at the beginning instead of doing it for every element. [DB] * ************************************************************************* *****/ void Determine_Blob_Textures(BLOB *Blob, VECTOR IPoint, int *Count, TEXTURE **Textures, DBL *Weights) { int i; int size; DBL rad2, sum; VECTOR V1, P; BLOB_ELEMENT *Element; BSPHERE_TREE *Tree; /* Make sure we have enough room in the textures/weights list. */ Reinitialize_Lighting_Code(Blob->Data->Number_Of_Components, &Textures, &Weights); /* Transform the point into the blob space. */ if (Blob->Trans != NULL) { MInvTransPoint(P, IPoint, Blob->Trans); } else { Assign_Vector(P, IPoint); } *Count = 0; if (Blob->Data->Tree == NULL) { /* There's no tree --> step through all elements. */ for (i = 0; i < Blob->Data->Number_Of_Components; i++) { Element = &Blob->Data->Entry[i]; determine_element_texture(Blob, Element, Blob->Element_Texture[i], P, Count, Textures, Weights); } } else { /* A tree exists --> step through the tree. */ size = 0; Queue[size++] = Blob->Data->Tree; while (size > 0) { Tree = Queue[--size]; /* Test if current node is a leaf. */ if (Tree->Entries <= 0) { determine_element_texture(Blob, (BLOB_ELEMENT *)Tree->Node, Blob->Element_Texture[((BLOB_ELEMENT *)Tree->Node)->index], P, Count, Textures, Weights); } else { /* Test all sub-nodes. */ for (i = 0; i < (int)Tree->Entries; i++) { /* Insert sub-node if we are inside. */ VSub(V1, P, Tree->Node[i]->C); VDot(rad2, V1, V1); if (rad2 <= Tree->Node[i]->r2) { insert_node(Tree->Node[i], &size); } } } } } /* Normalize weights so that their sum is 1. */ if (*Count > 0) { sum = 0.0; for (i = 0; i < *Count; i++) { sum += Weights[i]; } if (sum > 0.0) { for (i = 0; i < *Count; i++) { Weights[i] /= sum; } } } } /************************************************************************ ***** * * FUNCTION * * determine_element_texture * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * If the intersection point is inside the component calculate * the field density and store the element's texture and the field * value in the texture/weight list. * * CHANGES * * Sep 1994 : Creation. * ************************************************************************* *****/ static void determine_element_texture(BLOB *Blob, BLOB_ELEMENT *Element, TEXTURE *Texture, VECTOR P, int *Count, TEXTURE **Textures, DBL *Weights) { int i; DBL density; density = fabs(calculate_element_field(Element, P)); if (density > 0.0) { if (Texture == NULL) { Textures[*Count] = Blob->Texture; } else { Textures[*Count] = Texture; } /* Test if this texture is already used. */ for (i = 0; i < *Count; i++) { if (Textures[i] == Textures[*Count]) { /* Add current weight to already existing texture weight. */ Weights[i] += density; /* Any texture can only be in the list once --> exit. */ return; } } Weights[(*Count)++] = density; } } /************************************************************************ ***** * * FUNCTION * * Translate_Blob_Element * * INPUT * * Element - Pointer to blob element * Vector - Translation vector * * OUTPUT * * Object * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * Translate a blob element. * * CHANGES * * Sep 1994 : Creation. * ************************************************************************* *****/ void Translate_Blob_Element(BLOB_ELEMENT *Element, VECTOR Vector) { TRANSFORM Trans; Compute_Translation_Transform(&Trans, Vector); if (Element->Trans == NULL) { /* This is a sphere component. */ VAddEq(Element->O, Vector); Transform_Textures(Element->Texture, &Trans); } else { /* This is one of the other components. */ Transform_Blob_Element(Element, &Trans); } } /************************************************************************ ***** * * FUNCTION * * Rotate_Blob_Element * * INPUT * * Element - Pointer to blob element * Vector - Translation vector * * OUTPUT * * Object * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * Rotate a blob element. * * CHANGES * * Sep 1994 : Creation. * ************************************************************************* *****/ void Rotate_Blob_Element(BLOB_ELEMENT *Element, VECTOR Vector) { TRANSFORM Trans; Compute_Rotation_Transform(&Trans, Vector); if (Element->Trans == NULL) { /* This is a sphere component. */ MTransPoint(Element->O, Element->O, &Trans); Transform_Textures(Element->Texture, &Trans); } else { /* This is one of the other components. */ Transform_Blob_Element(Element, &Trans); } } /************************************************************************ ***** * * FUNCTION * * Scale_Blob_Element * * INPUT * * Element - Pointer to blob element * Vector - Translation vector * * OUTPUT * * Object * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * Scale a blob element. * * CHANGES * * Sep 1994 : Creation. * ************************************************************************* *****/ void Scale_Blob_Element(BLOB_ELEMENT *Element, VECTOR Vector) { TRANSFORM Trans; if ((Vector[X] != Vector[Y]) || (Vector[X] != Vector[Z])) { if (Element->Trans == NULL) { /* This is a sphere component --> change to ellipsoid component. */ Element->Type = BLOB_ELLIPSOID; Element->Trans = Create_Transform(); } } Compute_Scaling_Transform(&Trans, Vector); if (Element->Trans == NULL) { /* This is a sphere component. */ VScaleEq(Element->O, Vector[X]); Element->rad2 *= Sqr(Vector[X]); Transform_Textures(Element->Texture, &Trans); } else { /* This is one of the other components. */ Transform_Blob_Element(Element, &Trans); } } /************************************************************************ ***** * * FUNCTION * * Transform_Blob_Element * * INPUT * * Element - Pointer to blob element * Trans - Transformation * * OUTPUT * * Object * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * Transform a blob element. * * CHANGES * * Sep 1994 : Creation. * ************************************************************************* *****/ void Transform_Blob_Element(BLOB_ELEMENT *Element, TRANSFORM *Trans) { if (Element->Trans == NULL) { /* This is a sphere component --> change to ellipsoid component. */ Element->Type = BLOB_ELLIPSOID; Element->Trans = Create_Transform(); } Compose_Transforms(Element->Trans, Trans); Transform_Textures(Element->Texture, Trans); } /************************************************************************ ***** * * FUNCTION * * Invert_Blob_Element * * INPUT * * Element - Pointer to blob element * * OUTPUT * * Object * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * Invert blob element by negating its strength. * * CHANGES * * Sep 1994 : Creation. * ************************************************************************* *****/ void Invert_Blob_Element(BLOB_ELEMENT *Element) { Element->c[2] *= -1.0; } /************************************************************************ ***** * * FUNCTION * * Init_Blob_Queue * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * Init queues for blob intersections. * * CHANGES * * Dec 1994 : Creation. * ************************************************************************* *****/ void Init_Blob_Queue() { Queue = (BSPHERE_TREE **)POV_MALLOC(Max_Queue_Size*sizeof(BSPHERE_TREE *), "blob queue"); } /************************************************************************ ***** * * FUNCTION * * Destroy_Blob_Queue * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * Destroy queues for blob intersections. * * CHANGES * * Dec 1994 : Creation. * ************************************************************************* *****/ void Destroy_Blob_Queue() { if (Queue != NULL) { POV_FREE(Queue); } Queue = NULL; } /************************************************************************ ***** * * FUNCTION * * insert_node * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * Insert a node into the node queue. * * CHANGES * * Feb 1995 : Creation. * ************************************************************************* *****/ static void insert_node(BSPHERE_TREE *Node, int *size) { /* Resize queue if necessary. */ if (*size >= (int)Max_Queue_Size) { if (Max_Queue_Size >= INT_MAX/2) { Error("Blob queue overflow!"); } Max_Queue_Size *= 2; Queue = (BSPHERE_TREE **)POV_REALLOC(Queue, Max_Queue_Size*sizeof(BSPHERE_TREE *), "blob queue"); } Queue[(*size)++] = Node; } /************************************************************************ ***** * * FUNCTION * * Create_Blob_Element_Texture_List * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * Create a list of all textures in the blob. * * The list actually contains copies of the textures not * just references to them. * * CHANGES * * Mar 1996 : Created. * ************************************************************************* *****/ void Create_Blob_Element_Texture_List(BLOB *Blob, BLOB_LIST *BlobList, int npoints) { int i, element_count, count; BLOB_LIST *temp; if (npoints < 1) { Error("Need at least one component in a blob."); } /* Figure out how many components there will be. */ temp = BlobList; for (i = count = 0; i < npoints; i++) { if (temp->elem.Type & BLOB_CYLINDER) { count += 3; } else { count++; } temp = temp->next; /* Test for too many components. [DB 12/94] */ if (count >= MAX_BLOB_COMPONENTS) { Error("There are more than %d components in a blob.", MAX_BLOB_COMPONENTS); } } /* Allocate memory for components. */ Blob->Data = (BLOB_DATA *)POV_MALLOC(sizeof(BLOB_DATA), "blob data"); Blob->Data->References = 1; Blob->Data->Number_Of_Components = count; Blob->Data->Entry = (BLOB_ELEMENT *)POV_MALLOC(count*sizeof(BLOB_ELEMENT), "blob data"); Blob->Data->Intervals = (BLOB_INTERVAL *)POV_MALLOC(2*Blob->Data->Number_Of_Components*sizeof(BLOB_INTERVAL), "blob intervals"); Blob->Data->Tree = NULL; /* Init components. */ for (i = 0; i < count; i++) { init_blob_element(&Blob->Data->Entry[i]); } /* Allocate memory for list. */ Blob->Element_Texture = (TEXTURE **)POV_MALLOC(count*sizeof(TEXTURE *), "blob texture list"); for (i = 0; i < count; i++) { Blob->Element_Texture[i] = NULL; } for (i = element_count = 0; i < npoints; i++) { temp = BlobList; /* Copy textures. */ switch (temp->elem.Type) { case BLOB_ELLIPSOID : case BLOB_SPHERE : /* * Copy texture into element texture list. This is neccessary * because individual textures have to be transformed too if * copies of the blob are transformed. */ Blob->Element_Texture[element_count++] = Copy_Textures(temp->elem.Texture); break; case BLOB_CYLINDER : /* * Copy texture into element texture list. This is neccessary * because individual textures have to be transformed too if * copies of the blob are transformed. */ Blob->Element_Texture[element_count++] = Copy_Textures(temp->elem.Texture); Blob->Element_Texture[element_count++] = Copy_Textures(temp->elem.Texture); Blob->Element_Texture[element_count++] = Copy_Textures(temp->elem.Texture); break; } BlobList = BlobList->next; } }