home *** CD-ROM | disk | FTP | other *** search
/ Chestnut's Multimedia Mania / MM_MANIA.ISO / graphics / povsrc20 / bound.c < prev    next >
C/C++ Source or Header  |  1993-08-22  |  20KB  |  724 lines
  1. /****************************************************************************
  2. *                bound.c
  3. *
  4. *  This module implements the bounding slab calculations.
  5. *  This file was written by Alexander Enzmann.    He wrote the code for
  6. *  POV-Ray's bounding slabs and generously provided us these enhancements.
  7. *  The slab intersection code was further hacked by Eric Haines to speed it up.
  8. *
  9. *  Just so everyone knows where this came from, the code is VERY heavily
  10. *  based on the slab code from Mark VandeWettering's MTV raytracer.
  11. *  POV-Ray is just joining the crowd of admirers of Mark's contribution to
  12. *  the public domain. [ARE]
  13. *
  14. *  from Persistence of Vision Raytracer
  15. *  Copyright 1993 Persistence of Vision Team
  16. *---------------------------------------------------------------------------
  17. *  NOTICE: This source code file is provided so that users may experiment
  18. *  with enhancements to POV-Ray and to port the software to platforms other
  19. *  than those supported by the POV-Ray Team.  There are strict rules under
  20. *  which you are permitted to use this file.  The rules are in the file
  21. *  named POVLEGAL.DOC which should be distributed with this file. If
  22. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  23. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  24. *  Forum.  The latest version of POV-Ray may be found there as well.
  25. *
  26. * This program is based on the popular DKB raytracer version 2.12.
  27. * DKBTrace was originally written by David K. Buck.
  28. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  29. *
  30. *****************************************************************************/
  31.  
  32. #include "frame.h"
  33. #include "vector.h"
  34. #include "povproto.h"
  35.  
  36.   typedef struct {
  37.   int x,y,z ;
  38.   } 
  39. VECTORI, *pVECTORI ;
  40.  
  41.   typedef struct {
  42.   VECTOR slab_num ;
  43.   VECTOR slab_den ;
  44.   VECTORI nonzero ;
  45.   VECTORI positive ;
  46.   } 
  47. RAYINFO, *pRAYINFO ;
  48.  
  49.  
  50. extern FRAME Frame;
  51. extern long Bounds_Threshold;
  52. extern int Use_Slabs;
  53. static int Axis = 0;
  54. static unsigned long maxprimcount = 0;
  55.   typedef struct t_qelem {
  56.   DBL     q_key;
  57.   OBJECT *q_obj;
  58.   } 
  59. Qelem;
  60.  
  61. static int FindAxis PARAMS((OBJECT **Prims, unsigned long first,
  62. unsigned long last));
  63. static COMPOSITE *Create_Composite PARAMS((void));
  64. static int SortAndSplit PARAMS((OBJECT **Root, OBJECT **Prims,
  65. unsigned long *nPrims, unsigned long first, unsigned long last));
  66. static void PriorityQueueInsert PARAMS((Qelem *Queue, unsigned *Qsize,
  67. DBL key, OBJECT *obj));
  68. static void CheckAndEnqueue PARAMS((Qelem *Queue, unsigned *Qsize,
  69. OBJECT *obj, RAYINFO *rayinfo));
  70. static void PriorityQueueDelete PARAMS((Qelem *Queue, unsigned *Qsize,
  71. DBL *key, OBJECT **obj));
  72.  
  73. /* Should move these out of here... */
  74. unsigned long totalQueues = 0;
  75. unsigned long totalQueueResets = 0;
  76. unsigned long nChecked = 0;
  77. unsigned long nEnqueued = 0;
  78.  
  79. unsigned MAXQUEUE = 256;
  80.  
  81. METHODS Composite_Methods =
  82.   { 
  83.   NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, Destroy_Composite 
  84. };
  85.  
  86. void Destroy_Composite (Object)
  87. OBJECT *Object;
  88.   {
  89.   free (Object);
  90.   }
  91.  
  92. static COMPOSITE *Create_Composite ()
  93.   {
  94.   COMPOSITE *New;
  95.  
  96.   if ((New = (COMPOSITE *) malloc (sizeof (COMPOSITE))) == NULL)
  97.     MAError ("composite");
  98.  
  99.   INIT_OBJECT_FIELDS(New, COMPOSITE_OBJECT, &Composite_Methods)
  100.     return New;
  101.   }
  102.  
  103. void
  104. recompute_bbox(bbox, trans)
  105. BBOX *bbox;
  106. TRANSFORM *trans;
  107.   {
  108.   VECTOR lower_left, lengths, corner;
  109.   VECTOR mins, maxs;
  110.   int i;
  111.  
  112.   lower_left = bbox->Lower_Left;
  113.   lengths    = bbox->Lengths;
  114.   Make_Vector(&mins,  BOUND_HUGE,  BOUND_HUGE,  BOUND_HUGE);
  115.   Make_Vector(&maxs, -BOUND_HUGE, -BOUND_HUGE, -BOUND_HUGE);
  116.   for (i=1;i<=8;i++) 
  117.     {
  118.     corner = lower_left;
  119.     corner.x += ((i & 1) ? lengths.x : 0.0);
  120.     corner.y += ((i & 2) ? lengths.y : 0.0);
  121.     corner.z += ((i & 4) ? lengths.z : 0.0);
  122.     MTransPoint(&corner, &corner, trans);
  123.     if (corner.x < mins.x) mins.x = corner.x;
  124.     if (corner.x > maxs.x) maxs.x = corner.x;
  125.     if (corner.y < mins.y) mins.y = corner.y;
  126.     if (corner.y > maxs.y) maxs.y = corner.y;
  127.     if (corner.z < mins.z) mins.z = corner.z;
  128.     if (corner.z > maxs.z) maxs.z = corner.z;
  129.     }
  130.   bbox->Lower_Left = mins;
  131.   VSub(bbox->Lengths, maxs, mins);
  132.   }
  133.  
  134. void
  135. Recompute_Inverse_BBox(bbox, trans)
  136. BBOX *bbox;
  137. TRANSFORM *trans;
  138.   {
  139.   VECTOR lower_left, lengths, corner;
  140.   VECTOR mins, maxs;
  141.   int i;
  142.  
  143.   lower_left = bbox->Lower_Left;
  144.   lengths = bbox->Lengths;
  145.   Make_Vector(&mins,  BOUND_HUGE,  BOUND_HUGE,  BOUND_HUGE);
  146.   Make_Vector(&maxs, -BOUND_HUGE, -BOUND_HUGE, -BOUND_HUGE);
  147.   for (i=1;i<=8;i++) 
  148.     {
  149.     corner = lower_left;
  150.     corner.x += ((i & 1) ? lengths.x : 0.0);
  151.     corner.y += ((i & 2) ? lengths.y : 0.0);
  152.     corner.z += ((i & 4) ? lengths.z : 0.0);
  153.     MInvTransPoint(&corner, &corner, trans);
  154.     if (corner.x < mins.x) mins.x = corner.x;
  155.     if (corner.x > maxs.x) maxs.x = corner.x;
  156.     if (corner.y < mins.y) mins.y = corner.y;
  157.     if (corner.y > maxs.y) maxs.y = corner.y;
  158.     if (corner.z < mins.z) mins.z = corner.z;
  159.     if (corner.z > maxs.z) maxs.z = corner.z;
  160.     }
  161.   bbox->Lower_Left = mins;
  162.   VSub(bbox->Lengths, maxs, mins);
  163.   }
  164.  
  165. int CDECL compslabs(in_a, in_b)
  166. void *in_a;
  167. void *in_b;
  168.   {
  169.  
  170.   OBJECT **a, **b;
  171.   DBL am, bm;
  172.  
  173.   a = (OBJECT **)in_a;
  174.   b = (OBJECT **)in_b;
  175.  
  176.   switch (Axis) 
  177.   {
  178.   case 0:
  179.     am = 2.0 * (*a)->Bounds.Lower_Left.x + (*a)->Bounds.Lengths.x;
  180.     bm = 2.0 * (*b)->Bounds.Lower_Left.x + (*b)->Bounds.Lengths.x;
  181.     break;
  182.   case 1:
  183.     am = 2.0 * (*a)->Bounds.Lower_Left.y + (*a)->Bounds.Lengths.y;
  184.     bm = 2.0 * (*b)->Bounds.Lower_Left.y + (*b)->Bounds.Lengths.y;
  185.     break;
  186.   case 2:
  187.     am = 2.0 * (*a)->Bounds.Lower_Left.z + (*a)->Bounds.Lengths.z;
  188.     bm = 2.0 * (*b)->Bounds.Lower_Left.z + (*b)->Bounds.Lengths.z;
  189.     break;
  190.   default:
  191.     Error("Bad axis in compslabs\n");
  192.   }
  193.  
  194.   if (am < bm)
  195.     return -1;
  196.   else if (am == bm)
  197.     return 0;
  198.   else
  199.     return 1;
  200.   }
  201.  
  202. static int
  203. FindAxis(Prims, first, last)
  204. OBJECT **Prims;
  205. unsigned long first, last;
  206.   {
  207.   BBOX *bbox;
  208.   VECTOR mins, maxs;
  209.   unsigned long i;
  210.   int which;
  211.   DBL d = -BOUND_HUGE, e;
  212.  
  213.   Make_Vector(&mins, BOUND_HUGE, BOUND_HUGE, BOUND_HUGE);
  214.   Make_Vector(&maxs, -BOUND_HUGE, -BOUND_HUGE, -BOUND_HUGE);
  215.  
  216.   for (i=first;i<last;i++) 
  217.     {
  218.     bbox = &(Prims[i]->Bounds);
  219.     if (bbox->Lower_Left.x < mins.x)
  220.       mins.x = bbox->Lower_Left.x;
  221.     if (bbox->Lower_Left.x + bbox->Lengths.x > maxs.x)
  222.       maxs.x = bbox->Lower_Left.x;
  223.     if (bbox->Lower_Left.y < mins.y)
  224.       mins.y = bbox->Lower_Left.y;
  225.     if (bbox->Lower_Left.y + bbox->Lengths.y > maxs.y)
  226.       maxs.y = bbox->Lower_Left.y;
  227.     if (bbox->Lower_Left.z < mins.z)
  228.       mins.z = bbox->Lower_Left.z;
  229.     if (bbox->Lower_Left.z + bbox->Lengths.z > maxs.z)
  230.       maxs.z = bbox->Lower_Left.z;
  231.     }
  232.  
  233.   e = maxs.x - mins.x;
  234.   if (e > d) 
  235.     { 
  236.     d = e; which = 0; 
  237.   }
  238.   e = maxs.y - mins.y;
  239.   if (e > d) 
  240.     { 
  241.     d = e; which = 1; 
  242.   }
  243.   e = maxs.z - mins.z;
  244.   if (e > d) 
  245.     { 
  246.     d = e; which = 2; 
  247.   }
  248.  
  249.   return which;
  250.   }
  251.  
  252. static int
  253. SortAndSplit(Root, Prims, nPrims, first, last)
  254. OBJECT **Root;
  255. OBJECT **Prims;
  256. unsigned long *nPrims;
  257. unsigned long first;
  258. unsigned long last;
  259.   {
  260.   COMPOSITE *cd;
  261.   unsigned long size, i, j, m;
  262.   DBL dmin, dmax, tmin, tmax;
  263.  
  264.   Axis = FindAxis(Prims, first, last);
  265.   size = last - first;
  266.  
  267.   /* Actually, we could do this faster in several ways. we could use a
  268.       logn algorithm to find the median along the given axis, and then a
  269.       linear algorithm to partition along the axis. Oh well. */
  270.  
  271.   qsort((char *) (Prims + first), (int)size, sizeof(OBJECT *), compslabs);
  272.  
  273.   if (size <= BUNCHING_FACTOR) 
  274.     {
  275.     cd = Create_Composite();
  276.     cd->Entries = (unsigned short)size;
  277.  
  278.     for (i=0;i<size;i++) 
  279.       {
  280.       cd->Objects[i] = Prims[first+i];
  281.       /*
  282. printf("Extent of object %ld/%d: <%g, %g, %g> -> <%g, %g, %g>\n",
  283.        first+i, cd->Objects[i]->Type,
  284.        cd->Objects[i]->Bounds.Lower_Left.x,
  285.        cd->Objects[i]->Bounds.Lower_Left.y,
  286.        cd->Objects[i]->Bounds.Lower_Left.z,
  287.        cd->Objects[i]->Bounds.Lower_Left.x + cd->Objects[i]->Bounds.Lengths.x,
  288.        cd->Objects[i]->Bounds.Lower_Left.y + cd->Objects[i]->Bounds.Lengths.y,
  289.        cd->Objects[i]->Bounds.Lower_Left.z + cd->Objects[i]->Bounds.Lengths.z);
  290. */
  291.       }
  292.  
  293.     /* Check bounds in each direction */
  294.     /* First along the x axis */
  295.     dmin = BOUND_HUGE; dmax = -BOUND_HUGE;
  296.     for (j=0;j<size;j++) 
  297.       {
  298.       tmin = cd->Objects[j]->Bounds.Lower_Left.x;
  299.       tmax = tmin + cd->Objects[j]->Bounds.Lengths.x;
  300.       if (tmin < dmin) dmin = tmin;
  301.       if (tmax > dmax) dmax = tmax;
  302.       }
  303.     cd->Bounds.Lower_Left.x = dmin;
  304.     cd->Bounds.Lengths.x = dmax - dmin;
  305.  
  306.     /* Now along the y axis */
  307.     dmin = BOUND_HUGE; dmax = -BOUND_HUGE;
  308.     for (j=0;j<size;j++) 
  309.       {
  310.       tmin = cd->Objects[j]->Bounds.Lower_Left.y;
  311.       tmax = tmin + cd->Objects[j]->Bounds.Lengths.y;
  312.       if (tmin < dmin) dmin = tmin;
  313.       if (tmax > dmax) dmax = tmax;
  314.       }
  315.     cd->Bounds.Lower_Left.y = dmin;
  316.     cd->Bounds.Lengths.y = dmax - dmin;
  317.  
  318.     /* Lastly along the z axis */
  319.     dmin = BOUND_HUGE; dmax = -BOUND_HUGE;
  320.     for (j=0;j<size;j++) 
  321.       {
  322.       tmin = cd->Objects[j]->Bounds.Lower_Left.z;
  323.       tmax = tmin + cd->Objects[j]->Bounds.Lengths.z;
  324.       if (tmin < dmin) dmin = tmin;
  325.       if (tmax > dmax) dmax = tmax;
  326.       }
  327.     cd->Bounds.Lower_Left.z = dmin;
  328.     cd->Bounds.Lengths.z = dmax - dmin;
  329.  
  330.     *Root = (OBJECT *)cd;
  331.     if (*nPrims <= maxprimcount) 
  332.       {
  333.       Prims[*nPrims] = (OBJECT *)cd;
  334.       *nPrims += 1;
  335.       return 1;
  336.       }
  337.     else
  338.       Error("Too many primitives\n");
  339.     }
  340.   else 
  341.     {
  342.     m = (first + last) / 2;
  343.     SortAndSplit(Root, Prims, nPrims, first, m);
  344.     SortAndSplit(Root, Prims, nPrims, m , last);
  345.     return 0;
  346.     }
  347.   return -1;
  348.   }
  349.  
  350.   void
  351.   BuildBoundingSlabs(Root)
  352.     OBJECT **Root;
  353.   {
  354.   OBJECT **Prims, **prim, *head;
  355.   unsigned long nPrims;
  356.   unsigned long low, high;
  357.  
  358.   /* We have to start by counting how many frame level object there are */
  359.   head   = Frame.Objects;
  360.   nPrims = 0;
  361.   while (head != NULL) 
  362.     {
  363.     nPrims++;
  364.     head = head->Sibling;
  365.     }
  366.  
  367.   /* The total # of prims inflates around 150% when bounding objects are
  368.       generated.  If the 1.8 below proves to be too small, use 2.0.  The
  369.       inflation is never 200%. */
  370.   maxprimcount = (unsigned long)(1.8 * (nPrims + 1));
  371.  
  372.   /* Now allocate an array to hold references to these prims & any new
  373.      composite objects we may generate */
  374.   Prims = (OBJECT **)malloc((unsigned)maxprimcount * sizeof(OBJECT *));
  375.   if (Prims == NULL)
  376.     Error("Failed to allocate bounding slab reference information\n");
  377.  
  378.   /* Copy pointers to the objects into the array */
  379.   prim = Prims;
  380.   for (head=Frame.Objects;head!=NULL;head=head->Sibling) 
  381.     {
  382.     if (head->Type & LIGHT_SOURCE_OBJECT) 
  383.       {
  384.       /* Only bother with lights if they have an attached shape */
  385.       if (((LIGHT_SOURCE *)head)->Children != NULL)
  386.         *prim++ = ((LIGHT_SOURCE *)head)->Children;
  387.       else
  388.         nPrims--;
  389.       }
  390.     else
  391.       /* Normal sort of object - add it to the list */
  392.       *prim++ = head;
  393.     }
  394.  
  395.   /* Now do a sort on the objects, with the end result being a tree of
  396.       objects sorted along the x, y, and z axes */
  397.   low  = 0;
  398.   high = nPrims;
  399.   while (SortAndSplit(Root, Prims, &nPrims, low, high) == 0) 
  400.     {
  401.     low  = high;
  402.     high = nPrims;
  403.     }
  404.  
  405.   Use_Slabs = (nPrims >= Bounds_Threshold);
  406.  
  407.   /* Test */
  408.   /*
  409. printf("Extent of scene: <%g, %g, %g> -> <%g, %g, %g>\n",
  410.        (*Root)->Bounds.Lower_Left.x,
  411.        (*Root)->Bounds.Lower_Left.y,
  412.        (*Root)->Bounds.Lower_Left.z,
  413.        (*Root)->Bounds.Lower_Left.x + (*Root)->Bounds.Lengths.x,
  414.        (*Root)->Bounds.Lower_Left.y + (*Root)->Bounds.Lengths.y,
  415.        (*Root)->Bounds.Lower_Left.z + (*Root)->Bounds.Lengths.z);
  416. */
  417.  
  418.   /* Now we can get rid of the Prim array, and just use Root */
  419.   free(Prims);
  420.   }
  421.  
  422. static void
  423. PriorityQueueInsert(Queue, Qsize, key, obj)
  424. Qelem *Queue;
  425. unsigned *Qsize;
  426. DBL key;
  427. OBJECT *obj;
  428.   {
  429.   unsigned size;
  430.   int i;
  431.   Qelem tmp;
  432.  
  433.   totalQueues++;
  434.   (*Qsize)++;
  435.   size = *Qsize;
  436.   /* if (size > maxQueueSize) maxQueueSize = size; */
  437.   if (size >= MAXQUEUE)
  438.     Error("Priority queue overflow");
  439.   Queue[size].q_key = key;
  440.   Queue[size].q_obj = obj;
  441.  
  442.   i = size;
  443.   while (i > 1 && Queue[i].q_key < Queue[i/2].q_key) 
  444.     {
  445.     tmp = Queue[i];
  446.     Queue[i] = Queue[i/2];
  447.     Queue[i/2] = tmp;
  448.     i = i / 2;
  449.     }
  450.   }
  451.  
  452. static void
  453. PriorityQueueDelete(Queue, Qsize, key, obj)
  454. Qelem *Queue;
  455. unsigned *Qsize;
  456. DBL *key;
  457. OBJECT **obj;
  458.   {
  459.   Qelem tmp;
  460.   int i, j;
  461.   unsigned size;
  462.  
  463.   if (*Qsize == 0)
  464.     Error("priority queue is empty");
  465.  
  466.   *key = Queue[1].q_key;
  467.   *obj = Queue[1].q_obj;
  468.   Queue[1] = Queue[*Qsize];
  469.   (*Qsize)--;
  470.   size = *Qsize;
  471.  
  472.   i = 1 ;
  473.  
  474.   while (2 * i <= (int)size) 
  475.     {
  476.     if (2 * i == (int)size)
  477.       j = 2 * i;
  478.     else if (Queue[2*i].q_key < Queue[2*i+1].q_key)
  479.       j = 2 * i;
  480.     else
  481.       j = 2 * i + 1;
  482.  
  483.     if (Queue[i].q_key > Queue[j].q_key) 
  484.       {
  485.       tmp = Queue[i];
  486.       Queue[i] = Queue[j];
  487.       Queue[j] = tmp;
  488.       i = j;
  489.       }
  490.     else
  491.       break;
  492.     }
  493.   }
  494.  
  495. static void
  496. CheckAndEnqueue(Queue, Qsize, obj, rayinfo)
  497. Qelem *Queue;
  498. unsigned *Qsize;
  499. OBJECT *obj;
  500. RAYINFO *rayinfo;
  501.   {
  502.   DBL tmin, tmax;
  503.   DBL dmin, dmax ;
  504.  
  505.   nChecked++;
  506.  
  507.   if (rayinfo->nonzero.x ) 
  508.     {
  509.     if (rayinfo->positive.x ) 
  510.       {
  511.       dmin = (obj->Bounds.Lower_Left.x - rayinfo->slab_num.x) *
  512.       rayinfo->slab_den.x;
  513.       dmax = dmin + (obj->Bounds.Lengths.x  * rayinfo->slab_den.x);
  514.       if ( dmax < EPSILON ) return ;
  515.       } 
  516.     else 
  517.       {
  518.       dmax = (obj->Bounds.Lower_Left.x - rayinfo->slab_num.x) *
  519.         rayinfo->slab_den.x;
  520.       if ( dmax < EPSILON ) return ;
  521.       dmin = dmax + (obj->Bounds.Lengths.x  * rayinfo->slab_den.x);
  522.       }
  523.     if ( dmin > dmax ) return ;
  524.     }
  525.   else 
  526.     {
  527.     if ( ( rayinfo->slab_num.x < obj->Bounds.Lower_Left.x ) ||
  528.       ( rayinfo->slab_num.x >
  529.       obj->Bounds.Lengths.x + obj->Bounds.Lower_Left.x ) ) return ;
  530.     dmin = -BOUND_HUGE; dmax = BOUND_HUGE;
  531.     }
  532.  
  533.   if (rayinfo->nonzero.y ) 
  534.     {
  535.     if (rayinfo->positive.y ) 
  536.       {
  537.       tmin = (obj->Bounds.Lower_Left.y - rayinfo->slab_num.y) *
  538.       rayinfo->slab_den.y;
  539.       tmax = tmin + (obj->Bounds.Lengths.y  * rayinfo->slab_den.y);
  540.       } 
  541.     else 
  542.       {
  543.       tmax = (obj->Bounds.Lower_Left.y - rayinfo->slab_num.y) *
  544.         rayinfo->slab_den.y;
  545.       tmin = tmax + (obj->Bounds.Lengths.y  * rayinfo->slab_den.y);
  546.       }
  547.     /* unwrap the logic - do the dmin and dmax checks only when tmin and
  548.       tmax actually affect anything, also try to escape ASAP.  Better
  549.       yet, fold the logic below into the two branches above so as to
  550.       compute only what is needed.
  551.     */
  552.     /* you might even try tmax < EPSILON first (instead of second) for an
  553.           early quick out
  554.     */
  555.     if ( tmax < dmax ) 
  556.       {
  557.       if ( tmax < EPSILON ) return;
  558.       /* check bounds only if tmax changes dmax */
  559.       if ( tmin > dmin ) 
  560.         {
  561.         if ( tmin > tmax ) return ;
  562.         /* do this last in case it's not needed! */
  563.         dmin = tmin ;
  564.         } 
  565.       else 
  566.         {
  567.         if ( dmin > tmax ) return ;
  568.         }
  569.       /* do this last in case it's not needed! */
  570.       dmax = tmax ;
  571.       } 
  572.     else 
  573.       {
  574.       if ( tmin > dmin ) 
  575.         {
  576.         if ( tmin > dmax ) return ;
  577.         /* do this last in case it's not needed! */
  578.         dmin = tmin ;
  579.         } /* else nothing needs to happen, since dmin and dmax did not
  580.            change! */
  581.              
  582.       }
  583.     } 
  584.   else 
  585.     {
  586.     if (rayinfo->slab_num.y < obj->Bounds.Lower_Left.y ||
  587.       rayinfo->slab_num.y >
  588.       obj->Bounds.Lengths.y + obj->Bounds.Lower_Left.y ) return ;
  589.     }
  590.  
  591.   if (rayinfo->nonzero.z ) 
  592.     {
  593.     if (rayinfo->positive.z ) 
  594.       {
  595.       tmin = (obj->Bounds.Lower_Left.z - rayinfo->slab_num.z) *
  596.       rayinfo->slab_den.z;
  597.       tmax = tmin + (obj->Bounds.Lengths.z  * rayinfo->slab_den.z);
  598.       } 
  599.     else 
  600.       {
  601.       tmax = (obj->Bounds.Lower_Left.z - rayinfo->slab_num.z) *
  602.         rayinfo->slab_den.z;
  603.       tmin = tmax + (obj->Bounds.Lengths.z  * rayinfo->slab_den.z);
  604.       }
  605.     if ( tmax < dmax ) 
  606.       {
  607.       if ( tmax < EPSILON ) return;
  608.       /* check bounds only if tmax changes dmax */
  609.       if ( tmin > dmin ) 
  610.         {
  611.         if ( tmin > tmax ) return ;
  612.         /* do this last in case it's not needed! */
  613.         dmin = tmin ;
  614.         } 
  615.       else 
  616.         {
  617.         if ( dmin > tmax ) return ;
  618.         }
  619.       /* do this last in case it's not needed! */
  620.       dmax = tmax ;
  621.       } 
  622.     else 
  623.       {
  624.       if ( tmin > dmin ) 
  625.         {
  626.         if ( tmin > dmax ) return ;
  627.         /* do this last in case it's not needed! */
  628.         dmin = tmin ;
  629.         } /* else nothing needs to happen, since dmin and dmax did not
  630.            change! */
  631.              
  632.       }
  633.     } 
  634.   else 
  635.     {
  636.     if (rayinfo->slab_num.z < obj->Bounds.Lower_Left.z ||
  637.       rayinfo->slab_num.z >
  638.       obj->Bounds.Lengths.z + obj->Bounds.Lower_Left.z ) return ;
  639.     }
  640.  
  641.   PriorityQueueInsert(Queue, Qsize, dmin, obj);
  642.   nEnqueued++;
  643.   }
  644.  
  645. int
  646. Bounds_Intersect(Root, ray, Best_Intersection, Best_Object)
  647. OBJECT *Root;
  648. RAY *ray;
  649. INTERSECTION *Best_Intersection;
  650. OBJECT **Best_Object;
  651.   {
  652.   Qelem *Queue;
  653.   unsigned Qsize = 0;
  654.   int i;
  655.   OBJECT *cobj;
  656.   COMPOSITE *cdp;
  657.   RAYINFO rayinfo;
  658.   DBL t, key;
  659.   INTERSECTION New_Intersection;
  660.   int Intersection_Found = 0;
  661.  
  662.   Queue = (Qelem *)malloc(MAXQUEUE * sizeof(Qelem));
  663.   if (Queue == NULL)
  664.     Error("Failed to allocate priority queue\n");
  665.  
  666.   /* Create the direction vectors for this ray */
  667.   rayinfo.slab_num.x = ray->Initial.x;
  668.   rayinfo.slab_num.y = ray->Initial.y;
  669.   rayinfo.slab_num.z = ray->Initial.z;
  670.   if ( rayinfo.nonzero.x = ((t = ray->Direction.x) != 0.0) ) 
  671.     {
  672.     rayinfo.slab_den.x = 1.0 / t;
  673.     rayinfo.positive.x = ( ray->Direction.x > 0.0 ) ;
  674.     }
  675.   if ( rayinfo.nonzero.y = ((t = ray->Direction.y) != 0.0) ) 
  676.     {
  677.     rayinfo.slab_den.y = 1.0 / t;
  678.     rayinfo.positive.y = ( ray->Direction.y > 0.0 ) ;
  679.     }
  680.   if ( rayinfo.nonzero.z = ((t = ray->Direction.z) != 0.0) ) 
  681.     {
  682.     rayinfo.slab_den.z = 1.0 / t;
  683.     rayinfo.positive.z = ( ray->Direction.z > 0.0 ) ;
  684.     }
  685.  
  686.   /* start with an empty priority queue */
  687.   Qsize = 0;
  688.   totalQueueResets++;
  689.  
  690.   CheckAndEnqueue(Queue, &Qsize, Root, &rayinfo );
  691.  
  692.   for (;;)
  693.     {
  694.     if (Qsize == 0)
  695.       break;
  696.  
  697.     PriorityQueueDelete(Queue, &Qsize, &key, &cobj);
  698.  
  699.     if (key > Best_Intersection->Depth)
  700.       break;
  701.     else
  702.       if (cobj->Type & BOUNDING_OBJECT)
  703.         {
  704.         cdp = (COMPOSITE *)cobj;
  705.         for (i=0;(unsigned short)i < cdp->Entries;i++)
  706.           CheckAndEnqueue(Queue, &Qsize, cdp->Objects[i],
  707.             &rayinfo) ;
  708.         }
  709.       else
  710.         {
  711.         if (Intersection(&New_Intersection, cobj, ray))
  712.           if (New_Intersection.Depth < Best_Intersection->Depth)
  713.             {
  714.             *Best_Intersection = New_Intersection;
  715.             *Best_Object = cobj;
  716.             Intersection_Found = TRUE;
  717.             }
  718.         }
  719.     }
  720.  
  721.   free(Queue);
  722.   return Intersection_Found;
  723.   }
  724.