home *** CD-ROM | disk | FTP | other *** search
/ Otherware / Otherware_1_SB_Development.iso / mac / developm / source / macraysh.sit / Code / Source / blob.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-03-25  |  17.8 KB  |  714 lines
  1. /*
  2.  * blob.c
  3.  *
  4.  * Copyright (C) 1990, 1991, Mark Podlipec, Craig E. Kolb
  5.  * All rights reserved.
  6.  *
  7.  * This software may be freely copied, modified, and redistributed
  8.  * provided that this copyright notice is preserved on all copies.
  9.  *
  10.  * You may not distribute this software, in whole or in part, as part of
  11.  * any commercial product without the express consent of the authors.
  12.  *
  13.  * There is no warranty or other guarantee of fitness of this software
  14.  * for any purpose.  It is provided solely "as is".
  15.  *
  16.  * $Id: blob.c,v 4.0 91/07/17 14:36:07 kolb Exp Locker: kolb $
  17.  *
  18.  * $Log:    blob.c,v $
  19.  * Revision 4.0  91/07/17  14:36:07  kolb
  20.  * Initial version.
  21.  * 
  22.  */
  23. #include "geom.h"
  24. #include "blob.h"
  25.  
  26. Methods *iBlobMethods = NULL;
  27. static char blobName[] = "blob";
  28.  
  29. unsigned long BlobTests, BlobHits;
  30.  
  31. /*
  32.  * Blob/Metaball Description
  33.  *
  34.  * In this implementation a Blob is made up of a threshold T, and a 
  35.  * group of 1 or more Metaballs.  Each metaball 'i'  is defined by 
  36.  * its position (xi,yi,zi), its radius ri, and its strength si.
  37.  * Around each Metaball is a density function di(Ri) defined by:
  38.  *
  39.  *         di(Ri) =  c4i * Ri^4  +  c2i * Ri^2  + c0i     0 <= Ri <= ri
  40.  *         di(Ri) =  0                                    ri < Ri 
  41.  *
  42.  * where
  43.  *       c4i  = si / (ri^4)
  44.  *       c2i  = -(2 * si) / (ri^2)
  45.  *       c0i  = si
  46.  *       Ri^2 is the distance from a point (x,y,z) to the center of
  47.  *            the Metaball - (xi,yi,zi).
  48.  *
  49.  * The density function looks sort of like a W (Float-u) with the 
  50.  * center hump crossing the d-axis at  si  and the two humps on the
  51.  * side being tangent to the R-axis at  -ri  and  +ri. Only the 
  52.  * range  [0,ri] is being used. 
  53.  * I chose this function so that derivative = 0  at the points  0 and +ri. 
  54.  * This keeps the surface smooth when the densities are added. I also 
  55.  * wanted no  Ri^3  and  Ri  terms, since the equations are easier to 
  56.  * solve without them. This led to the symmetry about the d-axis and
  57.  * the derivative being equal to zero at -ri as well.
  58.  *
  59.  * The surface of the Blob is defined by:
  60.  *
  61.  *
  62.  *                  N
  63.  *                 ---  
  64.  *      F(x,y,z) = \    di(Ri)  = T
  65.  *                 /
  66.  *                 ---
  67.  *                 i=0
  68.  *
  69.  * where
  70.  *
  71.  *     di(Ri)    is   given above
  72.  *     Ri^2      =    (x-xi)^2  +  (y-yi)^2  + (z-zi)^2
  73.  *      N        is   number of Metaballs in Blob
  74.  *      T        is   the threshold
  75.  *    (xi,yi,zi) is   the center of Metaball i
  76.  *
  77.  */
  78.  
  79. /*****************************************************************************
  80.  * Create & return reference to a metaball.
  81.  */
  82. Blob *
  83. BlobCreate(T, mlist, npoints)
  84. Float T;
  85. MetaList *mlist;
  86. int npoints;
  87. {
  88.     Blob *blob;
  89.     int i;
  90.     MetaList *cur;
  91.  
  92. /* 
  93.  * There has to be at least one metaball in the blob.
  94.  * Note: if there's only one metaball, the blob 
  95.  * will be a sphere.
  96.  */
  97.     if (npoints < 1)
  98.     {
  99.         RLerror(RL_WARN, "blob field not correct.\n",0,0,0);
  100.         return (Blob *)NULL;
  101.     }
  102.  
  103. /*  
  104.  * Initialize primitive and Geom structures
  105.  */
  106.     blob = (Blob *)Malloc(sizeof(Blob));
  107.     blob->T = T;
  108.     blob->list=(MetaVector *) Malloc( (unsigned)(npoints*sizeof(MetaVector)) );
  109.     blob->num = npoints;
  110.  
  111. /*
  112.  * Initialize Metaball list
  113.  */
  114.     for(i=0;i<npoints;i++)
  115.     {
  116.         cur = mlist;
  117.         if ( (cur->mvec.c0 < EPSILON) || (cur->mvec.rs < EPSILON) )
  118.         {
  119.             RLerror(RL_WARN,"Degenerate metaball in blob (sr).\n",0,0,0);
  120.             return (Blob *)NULL;
  121.         }
  122.         /* store radius squared */
  123.         blob->list[i].rs = cur->mvec.rs * cur->mvec.rs;
  124.         /* Calculate and store coefficients for each metaball */
  125.         blob->list[i].c0 = cur->mvec.c0;
  126.         blob->list[i].c2 = -(2.0 * cur->mvec.c0) / blob->list[i].rs;
  127.         blob->list[i].c4 = cur->mvec.c0 /
  128.                 (blob->list[i].rs * blob->list[i].rs);
  129.         blob->list[i].x = cur->mvec.x;
  130.         blob->list[i].y = cur->mvec.y;
  131.         blob->list[i].z = cur->mvec.z;
  132.         mlist = mlist->next;
  133.         Free((voidstar)cur);
  134.     }
  135.     /* 
  136.      * Allocate room for Intersection Structures and
  137.      * Allocate room for an array of pointers to point to
  138.      * the Intersection Structures.
  139.      */
  140.     blob->ilist = (MetaInt *)Malloc( 2 * blob->num * sizeof(MetaInt));
  141.     blob->iarr = (MetaInt **)Malloc( 2 * blob->num * sizeof(MetaInt *));
  142.     return blob;
  143. }
  144.  
  145. Methods *
  146. BlobMethods()
  147. {
  148.     if (iBlobMethods == (Methods *)NULL) {
  149.         iBlobMethods = MethodsCreate();
  150.         iBlobMethods->create = (GeomCreateFunc *)BlobCreate;
  151.         iBlobMethods->methods = BlobMethods;
  152.         iBlobMethods->name = BlobName;
  153.         iBlobMethods->intersect = BlobIntersect;
  154.         iBlobMethods->normal = BlobNormal;
  155.         iBlobMethods->bounds = BlobBounds;
  156.         iBlobMethods->stats = BlobStats;
  157.         iBlobMethods->checkbounds = TRUE;
  158.         iBlobMethods->closed = TRUE;
  159.     }
  160.     return iBlobMethods;
  161. }
  162.  
  163. /*****************************************************************************
  164.  * Function used by qsort() when sorting the Ray/Blob intersection list
  165.  */
  166. int
  167. MetaCompare(A,B)
  168. char *A,*B;
  169. {
  170.     MetaInt **AA,**BB;
  171.  
  172.     AA = (MetaInt **) A;
  173.     BB = (MetaInt **) B;
  174.     if (AA[0]->bound == BB[0]->bound) return(0);
  175.     if (AA[0]->bound <  BB[0]->bound) return(-1);
  176.     return(1);  /* AA[0]->bound is > BB[0]->bound */
  177. }
  178.  
  179. /*****************************************************************************
  180.  * Ray/metaball intersection test.
  181.  */
  182. int
  183. BlobIntersect(blob, ray, mindist, maxdist)
  184. Blob *blob;
  185. Ray *ray;
  186. Float mindist, *maxdist;
  187. {
  188.     double c[5], s[4];
  189.     Float dist;
  190.     MetaInt *ilist,**iarr;
  191.     register int i,j,inum;
  192.     extern void qsort();
  193.  
  194.     BlobTests++;
  195.  
  196.     ilist = blob->ilist;
  197.     iarr = blob->iarr;
  198.  
  199. /*
  200.  * The first step in calculating the Ray/Blob intersection is to 
  201.  * divide the Ray into intervals such that only a fixed set of 
  202.  * Metaballs contribute to the density function along that interval. 
  203.  * This is done by finding the set of intersections between the Ray 
  204.  * and each Metaball's Sphere/Region of influence, which has a 
  205.  * radius  ri  and is centered at (xi,yi,zi).
  206.  * Intersection information is kept track of in the MetaInt 
  207.  * structure and consists of:
  208.  *
  209.  *   type    indicates whether this intersection is the start(R_START)
  210.  *           of a Region or the end(R_END) of one. 
  211.  *   pnt     the Metaball of this intersection
  212.  *   bound   the distance from Ray origin to this intersection
  213.  *
  214.  * This list is then sorted by  bound  and used later to find the Ray/Blob 
  215.  * intersection.
  216.  */
  217.  
  218.     inum = 0;
  219.     for(i=0; i < blob->num; i++)
  220.     {
  221.         register Float xadj, yadj, zadj;
  222.         register Float b, t, rs;
  223.         register Float dmin,dmax;
  224.  
  225.         rs   = blob->list[i].rs;
  226.         xadj = blob->list[i].x - ray->pos.x;
  227.         yadj = blob->list[i].y - ray->pos.y;
  228.         zadj = blob->list[i].z - ray->pos.z;
  229.  
  230.     /*
  231.      * Ray/Sphere of Influence intersection
  232.      */
  233.         b = xadj * ray->dir.x + yadj * ray->dir.y + zadj * ray->dir.z;
  234.         t = b * b - xadj * xadj - yadj * yadj - zadj * zadj + rs;
  235.  
  236.     /* 
  237.      * don't except imaginary or single roots. A single root is a ray
  238.      * tangent to the Metaball's Sphere/Region. The Metaball's
  239.      * contribution to the overall density function at this point is
  240.      * zero anyway.
  241.      */
  242.         if (t > 0.0) /* we have two roots */
  243.         {
  244.             t = sqrt(t);
  245.             dmin = b - t;
  246.     /* 
  247.      * only interested in stuff in front of ray origin 
  248.      */
  249.             if (dmin < mindist) dmin = mindist;
  250.             dmax = b + t;
  251.             if (dmax > dmin) /* we have a valid Region */
  252.             {
  253.         /*
  254.              * Initialize min/start and max/end Intersections Structures
  255.              * for this Metaball
  256.              */
  257.                 ilist[inum].type = R_START;
  258.                 ilist[inum].pnt = i;
  259.                 ilist[inum].bound = dmin;
  260.                 for(j=0;j<5;j++) ilist[inum].c[j] = 0.0;
  261.                 iarr[inum] = &(ilist[inum]);
  262.                 inum++;
  263.  
  264.                 ilist[inum].type = R_END;
  265.                 ilist[inum].pnt = i;
  266.                 ilist[inum].bound = dmax;
  267.                 for(j=0;j<5;j++) ilist[inum].c[j] = 0.0;
  268.                 iarr[inum] = &(ilist[inum]);
  269.                 inum++;
  270.             } /* End of valid Region */
  271.         } /* End of two roots */
  272.     } /* End of loop through metaballs */
  273.  
  274.     /*
  275.      * If there are no Ray/Metaball intersections there will 
  276.      * not be a Ray/Blob intersection. Exit now.
  277.      */
  278.     if (inum == 0)
  279.     {
  280.         return FALSE;
  281.     }
  282.  
  283.  
  284.     /* 
  285.      * Sort Intersection list. No sense using qsort if there's only
  286.      * two intersections.
  287.      *
  288.      * Note: we actually aren't sorting the Intersection structures, but
  289.      * an array of pointers to the Intersection structures. 
  290.      * This is faster than sorting the Intersection structures themselves.
  291.      */
  292.     if (inum==2)
  293.     {
  294.         MetaInt *t;
  295.         if (ilist[0].bound > ilist[1].bound)
  296.         {
  297.             t = iarr[0];
  298.             iarr[0] = iarr[1];
  299.             iarr[1] = t;
  300.         }
  301.     }
  302.     else
  303.         qsort((void *)iarr, (size_t)inum, (size_t)sizeof(MetaInt *),(__cmp_func)MetaCompare);
  304.  
  305.  
  306. /*
  307. * Finding the Ray/Blob Intersection
  308. * 
  309. * The non-zero part of the density function for each Metaball is
  310. * 
  311. *   di(Ri) = c4i * Ri^4  +  c2i * Ri^2  +  c0i 
  312. *
  313. * To find find the Ray/Blob intersection for one metaball
  314. * substitute for distance   
  315. *
  316. *     Ri^2 = (x-xi)^2 + (y-yi)^2 + (z-zi)^2
  317. * 
  318. * and then substitute the Ray equations:
  319. *   
  320. *     x  = x0 + x1 * t
  321. *     y  = y0 + y1 * t
  322. *     z  = z0 + z1 * t
  323. *
  324. * to get a big mess :^). Actually, it's a Quartic in t and it's fully 
  325. * listed farther down. Here's a short version:
  326. *
  327. *   c[4] * t^4  +  c[3] * t^3  +  c[2] * t^2  +  c[1] * t  +  c[0]  =  T
  328. *
  329. * Don't forget that the Blob is defined by the density being equal to 
  330. * the threshold T.
  331. * To calculate the intersection of a Ray and two or more Metaballs, 
  332. * the coefficients are calculated for each Metaball and then added 
  333. * together. We can do this since we're working with polynomials.
  334. * The points of intersection are the roots of the resultant equation.
  335. *
  336. * The algorithm loops through the intersection list, calculating
  337. * the coefficients if an intersection is the start of a Region and 
  338. * adding them to all intersections in that region.
  339. * When it detects a valid interval, it calculates the 
  340. * roots from the starting intersection's coefficients and if any of 
  341. * the roots are in the interval, the smallest one is returned.
  342. *
  343. */
  344.  
  345.     {
  346.         register Float *tmpc;
  347.         MetaInt *strt,*tmp;
  348.         register int istrt,itmp;
  349.         register int num,exitflag,inside;
  350.  
  351.     /*
  352.      * Start the algorithm outside the first interval
  353.      */
  354.         inside = 0;
  355.         istrt = 0; 
  356.         strt = iarr[istrt];
  357.         if (strt->type != R_START)
  358.             RLerror(RL_WARN,"MetaInt sanity check FAILED!\n",0,0,0);
  359.  
  360.         /*
  361.          * Loop through intersection. If a root is found the code
  362.          * will return at that point.
  363.          */
  364.         while( istrt < inum )
  365.         {
  366.             /* 
  367.              * Check for multiple intersections  at the same point.
  368.              * This is also where the coefficients are calculated
  369.              * and spread throughout that Metaball's sphere of
  370.              * influence.
  371.              */
  372.             do
  373.             {
  374.                 /* find out which metaball */
  375.                 i = strt->pnt;
  376.                 /* only at starting boundaries do this */
  377.                 if (strt->type == R_START)
  378.                 {
  379.                     register MetaVector *ml;
  380.                     register Float a1,a0;
  381.                     register Float xd,yd,zd;
  382.                     register Float c4,c2,c0;
  383.  
  384.                     /* we're inside */
  385.                     inside++;
  386.  
  387.     /*  As promised, the full equations
  388.      *
  389.      *   c[4] = c4*a2*a2; 
  390.      *   c[3] = 4.0*c4*a1*a2;
  391.      *   c[2] = 4.0*c4*a1*a1 + 2.0*c4*a2*a0 + c2*a2;
  392.      *   c[1] = 4.0*c4*a1*a0 + 2.0*c2*a1;
  393.      *   c[0] = c4*a0*a0 + c2*a0 + c0;
  394.      *
  395.      * where
  396.      *        a2 = (x1*x1 + y1*y1 + z1*z1) = 1.0 because the ray 
  397.      *                                           is normalized
  398.      *        a1 = (xd*x1 + yd*y1 + zd*z1)
  399.      *        a0 = (xd*xd + yd*yd + zd*zd)
  400.      *        xd = (x0 - xi)
  401.      *        yd = (y0 - yi)
  402.      *        zd = (z0 - zi)
  403.      *        (xi,yi,zi) is center of Metaball
  404.      *        (x0,y0,z0) is Ray origin
  405.      *        (x1,y1,z1) is normalized Ray direction
  406.      *        c4,c2,c0   are the coefficients for the
  407.      *                       Metaball's density function
  408.      *
  409.      */
  410.  
  411.                     /*
  412.                      * Some compilers would recalculate
  413.                      * this each time its used.
  414.                      */
  415.                     ml = &(blob->list[i]);
  416.  
  417.                     xd = ray->pos.x - ml->x;
  418.                     yd = ray->pos.y - ml->y;
  419.                     zd = ray->pos.z - ml->z;
  420.                     a1 = (xd * ray->dir.x + yd * ray->dir.y
  421.                         + zd * ray->dir.z);
  422.                     a0 = (xd * xd + yd * yd + zd * zd);
  423.  
  424.                     c0 = ml->c0;
  425.                     c2 = ml->c2;
  426.                     c4 = ml->c4;
  427.  
  428.                     c[4] = c4;
  429.                     c[3] = 4.0*c4*a1;
  430.                     c[2] = 2.0*c4*(2.0*a1*a1 + a0) + c2;
  431.                     c[1] = 2.0*a1*(2.0*c4*a0 + c2);
  432.                     c[0] = c4*a0*a0 + c2*a0 + c0;
  433.  
  434.         /* 
  435.          * Since we've gone through the trouble of calculating the 
  436.          * coefficients, put them where they'll be used.
  437.          * Starting at the current intersection(It's also the start of
  438.          * a region), add the coefficients to each intersection until
  439.          * we reach the end of this region.
  440.          */
  441.                     itmp = istrt; 
  442.                     tmp = strt;
  443.                     while( (tmp->pnt != i) ||
  444.                            (tmp->type != R_END) )
  445.                     {
  446.                         tmpc = tmp->c;
  447.                         for(j=0;j<5;j++) {
  448.                             *tmpc += c[j];
  449.                             tmpc++ ;
  450.                         }
  451.                         itmp++; 
  452.                         tmp = iarr[itmp];
  453.                     }
  454.  
  455.                 } /* End of start of a Region */ 
  456.                 /* 
  457.                  * Since the intersection wasn't the start
  458.                  * of a region, it must the end of one.
  459.                  */
  460.                 else inside--;
  461.  
  462.         /*
  463.          * If we are inside a region(or regions) and the next
  464.          * intersection is not at the same place, then we have
  465.          * the start of an interval. Set the exitflag.
  466.          */
  467.                 if ((inside > 0)
  468.                     && (strt->bound != iarr[istrt+1]->bound) )
  469.                     exitflag = 1;
  470.                 else 
  471.                 /* move to next intersection */
  472.                 {
  473.                     istrt++; 
  474.                     strt = iarr[istrt];
  475.                     exitflag = 0;
  476.                 }
  477.         /* 
  478.          * Check to see if we're at the end. If so then exit with
  479.          * no intersection found.
  480.          */
  481.                 if (istrt >= inum)
  482.                 {
  483.                     return FALSE;
  484.                 }
  485.             } while(!exitflag);
  486.                 /* End of Search for valid interval */
  487.  
  488.         /*
  489.          * Find Roots along this interval
  490.          */
  491.  
  492.             /* get coefficients from Intersection structure */
  493.             tmpc = strt->c;
  494.             for(j=0;j<5;j++) c[j] = *tmpc++;
  495.  
  496.             /* Don't forget to put in threshold */
  497.             c[0] -= blob->T;
  498.  
  499.             /* Use Graphic Gem's Quartic Root Finder */
  500.             num = SolveQuartic(c,s);
  501.  
  502.         /*
  503.          * If there are roots, search for roots within the interval.
  504.          */
  505.             dist = 0.0;
  506.             if (num>0)
  507.             {
  508.                 for(j=0;j<num;j++)
  509.                 {
  510.         /* 
  511.          * Not sure if EPSILON is truly needed, but it might cause
  512.          * small cracks between intervals in some cases. In any case 
  513.          * I don't believe it hurts.
  514.          */
  515.                     if ((s[j] >= (strt->bound - EPSILON))
  516.                         && (s[j] <= (iarr[istrt+1]->bound +
  517.                             EPSILON) ) )
  518.                     {
  519.                         if (dist == 0.0)
  520.                         /* first valid root */
  521.                             dist = s[j];
  522.                         else if (s[j] < dist)
  523.                          /* else only if smaller */
  524.                             dist = s[j];
  525.                     }
  526.                 }
  527.                 /*
  528.                  * Found a valid root 
  529.                  */
  530.                 if (dist > mindist && dist < *maxdist)
  531.                 {
  532.                     *maxdist = dist;
  533.                     BlobHits++;
  534.                     return TRUE;
  535.                     /* Yeah! Return valid root */
  536.                 }
  537.             }
  538.  
  539.         /* 
  540.          * Move to next intersection
  541.          */
  542.             istrt++;
  543.             strt = iarr[istrt];
  544.  
  545.         } /* End of loop through Intersection List */
  546.     } /* End of Solving for Ray/Blob Intersection */
  547.  
  548.     /* 
  549.      * return negative
  550.      */
  551.     return FALSE;
  552. }
  553.  
  554.  
  555. /***********************************************
  556.  * Find the Normal of a Blob at a given point
  557.  *
  558.  * The normal of a surface at a point (x0,y0,z0) is the gradient of that
  559.  * surface at that point. The gradient is the vector 
  560.  *
  561.  *            Fx(x0,y0,z0) , Fy(x0,y0,z0) , Fz(x0,y0,z0)
  562.  *
  563.  * where Fx(),Fy(),Fz() are the partial dervitives of F() with respect
  564.  * to x,y and z respectively. Since the surface of a Blob is made up
  565.  * of the sum of one or more polynomials, the normal of a Blob at a point
  566.  * is the sum of the gradients of the individual polynomials at that point.
  567.  * The individual polynomials in this case are di(Ri) i = 0,...,num .
  568.  *
  569.  * To find the gradient of the contributing polynomials
  570.  * take di(Ri) and substitute U = Ri^2;
  571.  *
  572.  *      di(U)    = c4i * U^2  +  c2i * U  + c0i
  573.  *
  574.  *      dix(U)   = 2 * c4i * Ux * U  +  c2i * Ux
  575.  *
  576.  *        U      = (x-xi)^2 + (y-yi)^2 + (z-zi)^2 
  577.  *
  578.  *        Ux     = 2 * (x-xi) 
  579.  *
  580.  * Rearranging we get
  581.  *  
  582.  *    dix(x0,y0,z0) = 4 * c4i * xdi * disti + 2 * c2i * xdi
  583.  *    diy(x0,y0,z0) = 4 * c4i * ydi * disti + 2 * c2i * ydi
  584.  *    diz(x0,y0,z0) = 4 * c4i * zdi * disti + 2 * c2i * zdi
  585.  * 
  586.  * where
  587.  *         xdi       =   x0-xi
  588.  *         ydi       =   y0-yi
  589.  *         zdi       =   y0-yi
  590.  *        disti      =   xdi*xdi + ydi*ydi + zdi*zdi   
  591.  *       (xi,yi,zi)  is  the center of Metaball i
  592.  *       c4i,c2i,c0i are the coefficients of Metaball i's density function
  593.  */
  594. int
  595. BlobNormal(blob, pos, nrm, gnrm)
  596. Blob *blob;
  597. Vector *pos, *nrm, *gnrm;
  598. {
  599.     register int i;
  600.  
  601.     /*  
  602.      * Initialize normals to zero 
  603.      */
  604.     nrm->x = nrm->y = nrm->z = 0.0;
  605.     /*
  606.      * Loop through Metaballs. If the point is within a Metaball's
  607.      * Sphere of influence, calculate the gradient and add it to the
  608.      * normals
  609.      */
  610.     for(i=0;i < blob->num; i++)
  611.     {
  612.         register MetaVector *sl;
  613.         register Float dist,xd,yd,zd;
  614.  
  615.         sl = &(blob->list[i]);
  616.         xd = pos->x - sl->x;
  617.         yd = pos->y - sl->y;
  618.         zd = pos->z - sl->z;
  619.  
  620.         dist  = xd*xd + yd*yd + zd*zd;
  621.         if (dist <= sl->rs )
  622.         {
  623.             register Float temp;
  624.  
  625.             /* temp is negative so normal points out of blob */
  626.             temp = -2.0 * (2.0 * sl->c4 * dist  +  sl->c2);
  627.             nrm->x += xd * temp;
  628.             nrm->y += yd * temp;
  629.             nrm->z += zd * temp;
  630.  
  631.         }
  632.     }
  633.     (void)VecNormalize(nrm);
  634.     *gnrm = *nrm;
  635.     return FALSE;
  636. }
  637.  
  638.  
  639. /*****************************************************************************
  640.  * Calculate the extent of the Blob
  641.  */
  642. void
  643. BlobBounds(blob, bounds)
  644. Blob *blob;
  645. Float bounds[2][3];
  646. {
  647.     Float r;
  648.     Float minx,miny,minz,maxx,maxy,maxz;
  649.     int i;
  650.  
  651.     /*
  652.      * Loop through Metaballs to find the minimun and maximum in each
  653.      * direction.
  654.      */
  655.     for( i=0;  i< blob->num;  i++)
  656.     {
  657.         register Float d;
  658.  
  659.         r = sqrt(blob->list[i].rs);
  660.         if (i==0)
  661.         {
  662.             minx = blob->list[i].x - r;
  663.             miny = blob->list[i].y - r;
  664.             minz = blob->list[i].z - r;
  665.             maxx = blob->list[i].x + r;
  666.             maxy = blob->list[i].y + r;
  667.             maxz = blob->list[i].z + r;
  668.         }
  669.         else
  670.         {
  671.             d = blob->list[i].x - r; 
  672.             if (d<minx) minx = d;
  673.             d = blob->list[i].y - r; 
  674.             if (d<miny) miny = d;
  675.             d = blob->list[i].z - r; 
  676.             if (d<minz) minz = d;
  677.             d = blob->list[i].x + r; 
  678.             if (d>maxx) maxx = d;
  679.             d = blob->list[i].y + r; 
  680.             if (d>maxy) maxy = d;
  681.             d = blob->list[i].z + r; 
  682.             if (d>maxz) maxz = d;
  683.         }
  684.     }
  685.     bounds[LOW][X]  = minx;
  686.     bounds[HIGH][X] = maxx;
  687.     bounds[LOW][Y]  = miny;
  688.     bounds[HIGH][Y] = maxy;
  689.     bounds[LOW][Z]  = minz;
  690.     bounds[HIGH][Z] = maxz;
  691. }
  692.  
  693. char *
  694. BlobName()
  695. {
  696.     return blobName;
  697. }
  698.  
  699. void
  700. BlobStats(tests, hits)
  701. unsigned long *tests, *hits;
  702. {
  703.     *tests = BlobTests;
  704.     *hits = BlobHits;
  705. }
  706.  
  707. void
  708. BlobMethodRegister(meth)
  709. UserMethodType meth;
  710. {
  711.     if (iBlobMethods)
  712.         iBlobMethods->user = meth;
  713. }
  714.