home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Professional / OS2PRO194.ISO / os2 / graphic / csg_rt / shape.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-02-11  |  28.8 KB  |  1,245 lines
  1. /*
  2.  
  3. SHAPE.C  General shapes
  4.  
  5. */
  6.  
  7. /*...sincludes:0:*/
  8. #include <stdio.h>
  9. #include <stdlib.h>
  10. #include <stddef.h>
  11. #include <malloc.h>
  12. #include <memory.h>
  13. #include <math.h>
  14. #include "standard.h"
  15. #include "rt.h"
  16. #include "fio.h"
  17. #include "tex.h"
  18. #include "vector.h"
  19. #include "rgbvec.h"
  20. #include "col.h"
  21. #include "surf.h"
  22. #include "sil.h"
  23. #include "plane.h"
  24. #include "sphere.h"
  25. #include "quad.h"
  26. #define    _SHAPE_
  27. #include "shape.h"
  28.  
  29. /*...vrt\46\h:0:*/
  30. /*...vfio\46\h:0:*/
  31. /*...vtex\46\h:0:*/
  32. /*...vvector\46\h:0:*/
  33. /*...vrgbvec\46\h:0:*/
  34. /*...vcol\46\h:0:*/
  35. /*...vsurf\46\h:0:*/
  36. /*...vsil\46\h:0:*/
  37. /*...vplane\46\h:0:*/
  38. /*...vsphere\46\h:0:*/
  39. /*...vquad\46\h:0:*/
  40. /*...vshape\46\h:0:*/
  41. /*...e*/
  42.  
  43. /*...screate_plane_shape  \45\ create a shape tree consisting of a single plane:0:*/
  44. SHAPE    *create_plane_shape(PLANE *plane, SURF *surf)
  45.     {
  46.     SHAPE    *shape;
  47.  
  48.     if ( (shape = malloc(sizeof(SHAPE))) == NULL )
  49.         return ( NULL );
  50.  
  51.     shape -> stype   = STYPE_PLANE;
  52.     shape -> u.plane = plane;
  53.     shape -> surf    = surf;
  54.     return ( shape );
  55.     }
  56. /*...e*/
  57. /*...screate_sphere_shape \45\ create a shape tree consisting of a single sphere:0:*/
  58. SHAPE    *create_sphere_shape(SPHERE *sphere, SURF *surf)
  59.     {
  60.     SHAPE    *shape;
  61.  
  62.     if ( (shape = malloc(sizeof(SHAPE))) == NULL )
  63.         return ( NULL );
  64.  
  65.     shape -> stype    = STYPE_SPHERE;
  66.     shape -> u.sphere = sphere;
  67.     shape -> surf     = surf;
  68.     return ( shape );
  69.     }
  70. /*...e*/
  71. /*...screate_quad_shape   \45\ create a shape tree consisting of a single quadratic:0:*/
  72. SHAPE    *create_quad_shape(QUAD *quad, SURF *surf)
  73.     {
  74.     SHAPE    *shape;
  75.  
  76.     if ( (shape = malloc(sizeof(SHAPE))) == NULL )
  77.         return ( NULL );
  78.  
  79.     shape -> stype  = STYPE_QUAD;
  80.     shape -> u.quad = quad;
  81.     shape -> surf   = surf;
  82.     return ( shape );
  83.     }
  84. /*...e*/
  85. /*...screate_bin_shape    \45\ create a binary shape tree made of 2 shapes:0:*/
  86. SHAPE    *create_bin_shape(STYPE stype, SHAPE *shape0, SHAPE *shape1)
  87.     {
  88.     SHAPE    *shape;
  89.  
  90.     if ( (shape = malloc(sizeof(SHAPE))) == NULL )
  91.         return ( NULL );
  92.  
  93.     shape -> stype        = stype;
  94.     shape -> u.shapes [0] = shape0;
  95.     shape -> u.shapes [1] = shape1;
  96.     return ( shape );
  97.     }
  98. /*...e*/
  99. /*...scopy_shape          \45\ make a complete copy of a shape tree:0:*/
  100. SHAPE    *copy_shape(SHAPE *shape)
  101.     {
  102.     SHAPE    *copy;
  103.  
  104.     if ( (copy = malloc(sizeof(SHAPE))) == NULL )
  105.         return ( NULL );
  106.  
  107.     copy -> stype = shape -> stype;
  108.  
  109.     if ( shape -> stype <= STYPE_QUAD )
  110.         /* Leaf node */
  111.         {
  112.         if ( (copy -> surf = copy_surf(shape -> surf)) == NULL )
  113.             {
  114.             free(copy);
  115.             return ( NULL );
  116.             }
  117.         switch ( shape -> stype )
  118.             {
  119. /*...sSTYPE_PLANE:24:*/
  120. case STYPE_PLANE:
  121.     if ( (copy -> u.plane = copy_plane(shape -> u.plane)) == NULL )
  122.         {
  123.         destroy_surf(copy -> surf);
  124.         free(copy);
  125.         return ( NULL );
  126.         }
  127.     break;
  128. /*...e*/
  129. /*...sSTYPE_SPHERE:24:*/
  130. case STYPE_SPHERE:
  131.     if ( (copy -> u.sphere = copy_sphere(shape -> u.sphere)) == NULL )
  132.         {
  133.         destroy_surf(copy -> surf);
  134.         free(copy);
  135.         return ( NULL );
  136.         }
  137.     break;
  138. /*...e*/
  139. /*...sSTYPE_QUAD:24:*/
  140. case STYPE_QUAD:
  141.     if ( (copy -> u.quad = copy_quad(shape -> u.quad)) == NULL )
  142.         {
  143.         destroy_surf(copy -> surf);
  144.         free(copy);
  145.         return ( NULL );
  146.         }
  147.     break;
  148. /*...e*/
  149.             }
  150.         }
  151.     else
  152. /*...sboolean combinations:16:*/
  153. {
  154. if ( (copy -> u.shapes [0] = copy_shape(shape -> u.shapes [0])) == NULL )
  155.     {
  156.     free(copy);
  157.     return ( NULL );
  158.     }
  159. if ( (copy -> u.shapes [1] = copy_shape(shape -> u.shapes [1])) == NULL )
  160.     {
  161.     free(copy -> u.shapes [0]);
  162.     free(copy);
  163.     return ( NULL );
  164.     }
  165. }
  166. /*...e*/
  167.  
  168.     return ( copy );
  169.     }
  170. /*...e*/
  171. /*...sdestroy_shape       \45\ delete a shape tree:0:*/
  172. void    destroy_shape(SHAPE *shape)
  173.     {
  174.     if ( shape -> stype <= STYPE_QUAD )
  175.         /* Leaf node */
  176.         {
  177.         destroy_surf(shape -> surf);
  178.         switch ( shape -> stype )
  179.             {
  180. /*...sSTYPE_PLANE:24:*/
  181. case STYPE_PLANE:
  182.     destroy_plane(shape -> u.plane);
  183.     break;
  184. /*...e*/
  185. /*...sSTYPE_SPHERE:24:*/
  186. case STYPE_SPHERE:
  187.     destroy_sphere(shape -> u.sphere);
  188.     break;
  189. /*...e*/
  190. /*...sSTYPE_QUAD:24:*/
  191. case STYPE_QUAD:
  192.     destroy_quad(shape -> u.quad);
  193.     break;
  194. /*...e*/
  195.             }
  196.         }
  197.     else
  198.         {
  199.         destroy_shape(shape -> u.shapes [0]);
  200.         destroy_shape(shape -> u.shapes [1]);
  201.         }
  202.  
  203.     free(shape);
  204.     }
  205. /*...e*/
  206. /*...ssphere_to_quad      \45\ convert sphere to quad:0:*/
  207. /*
  208. Spheres cannot be rescaled in each direction and remain a sphere.
  209. Hence this code will convert a sphere into the equivelent quadratic.
  210.  
  211.                  2      2      2
  212.     Sphere:        (x-a) +(y-b) +(z-c) -d <=0
  213.  
  214.              2      2  2      2  2      2
  215.             x -2ax+a +y -2by+b +z -2cz+c -d <= 0
  216.  
  217.               2   2   2
  218.     Quadratic:    ax +by +cz +dxy+eyz+fzx+gx+hy+iz+j<=0
  219.  
  220. */
  221.  
  222. static BOOLEAN sphere_to_quad(SHAPE *shape)
  223.     {
  224.     SPHERE *sphere = shape -> u.sphere;
  225.     QUAD *quad;
  226.  
  227.     if ( (quad = create_quad(1.0, 1.0, 1.0, 0.0, 0.0, 0.0,
  228.          -2.0 * sphere -> a, -2.0 * sphere -> b, -2.0 * sphere -> c,
  229.          - sphere -> d)) == NULL )
  230.         return ( FALSE );
  231.  
  232.     destroy_sphere(shape -> u.sphere);
  233.     shape -> stype = STYPE_QUAD;
  234.     shape -> u.quad = quad;
  235.  
  236.     return ( TRUE );
  237.     }
  238. /*...e*/
  239.  
  240. /*...strans_x             \45\ translate by amount in x direction:0:*/
  241. void    trans_x(SHAPE *shape, double t)
  242.     {
  243.     if ( shape -> stype <= STYPE_QUAD )
  244.         {
  245.         switch ( shape -> stype )
  246.             {
  247.             case STYPE_PLANE:
  248.                 trans_x_plane(shape -> u.plane, t);
  249.                 break;
  250.             case STYPE_SPHERE:
  251.                 trans_x_sphere(shape -> u.sphere, t);
  252.                 break;
  253.             case STYPE_QUAD:
  254.                 trans_x_quad(shape -> u.quad, t);
  255.                 break;
  256.             }
  257.         trans_x_surf(shape -> surf, t);
  258.         }
  259.     else
  260.         {
  261.         trans_x(shape -> u.shapes [0], t);
  262.         trans_x(shape -> u.shapes [1], t);
  263.         }
  264.     }
  265. /*...e*/
  266. /*...strans_y             \45\ translate by amount in y direction:0:*/
  267. void    trans_y(SHAPE *shape, double t)
  268.     {
  269.     if ( shape -> stype <= STYPE_QUAD )
  270.         {
  271.         switch ( shape -> stype )
  272.             {
  273.             case STYPE_PLANE:
  274.                 trans_y_plane(shape -> u.plane, t);
  275.                 break;
  276.             case STYPE_SPHERE:
  277.                 trans_y_sphere(shape -> u.sphere, t);
  278.                 break;
  279.             case STYPE_QUAD:
  280.                 trans_y_quad(shape -> u.quad, t);
  281.                 break;
  282.             }
  283.         trans_y_surf(shape -> surf, t);
  284.         }
  285.     else
  286.         {
  287.         trans_y(shape -> u.shapes [0], t);
  288.         trans_y(shape -> u.shapes [1], t);
  289.         }
  290.     }
  291. /*...e*/
  292. /*...strans_z             \45\ translate by amount in z direction:0:*/
  293. void    trans_z(SHAPE *shape, double t)
  294.     {
  295.     if ( shape -> stype <= STYPE_QUAD )
  296.         {
  297.         switch ( shape -> stype )
  298.             {
  299.             case STYPE_PLANE:
  300.                 trans_z_plane(shape -> u.plane, t);
  301.                 break;
  302.             case STYPE_SPHERE:
  303.                 trans_z_sphere(shape -> u.sphere, t);
  304.                 break;
  305.             case STYPE_QUAD:
  306.                 trans_z_quad(shape -> u.quad, t);
  307.                 break;
  308.             }
  309.         trans_z_surf(shape -> surf, t);
  310.         }
  311.     else
  312.         {
  313.         trans_z(shape -> u.shapes [0], t);
  314.         trans_z(shape -> u.shapes [1], t);
  315.         }
  316.     }
  317. /*...e*/
  318. /*...strans               \45\ translate by vector:0:*/
  319. void    trans(SHAPE *shape, VECTOR v)
  320.     {
  321.     trans_x(shape, v.x);
  322.     trans_y(shape, v.y);
  323.     trans_z(shape, v.z);
  324.     }
  325. /*...e*/
  326. /*...sscale_x             \45\ scale by factor in x direction:0:*/
  327. BOOLEAN    scale_x(SHAPE *shape, double factor)
  328.     {
  329.     if ( factor == 1.0 )
  330.         return ( TRUE );
  331.  
  332.     if ( shape -> stype <= STYPE_QUAD )
  333.         {
  334.         switch ( shape -> stype )
  335.             {
  336.             case STYPE_PLANE:
  337.                 scale_x_plane(shape -> u.plane, factor);
  338.                 break;
  339.             case STYPE_SPHERE:
  340.                 if ( !sphere_to_quad(shape) )
  341.                     return ( FALSE );
  342.                 scale_x_quad(shape -> u.quad, factor);
  343.                 break;
  344.             case STYPE_QUAD:
  345.                 scale_x_quad(shape -> u.quad, factor);
  346.                 break;
  347.             }
  348.         scale_x_surf(shape -> surf, factor);
  349.         }
  350.     else
  351.         {
  352.         if ( !scale_x(shape -> u.shapes [0], factor) ||
  353.              !scale_x(shape -> u.shapes [1], factor) )
  354.             return ( FALSE );
  355.         }
  356.     return ( TRUE );
  357.     }
  358. /*...e*/
  359. /*...sscale_y             \45\ scale by factor in y direction:0:*/
  360. BOOLEAN    scale_y(SHAPE *shape, double factor)
  361.     {
  362.     if ( factor == 1.0 )
  363.         return ( TRUE );
  364.  
  365.     if ( shape -> stype <= STYPE_QUAD )
  366.         {
  367.         switch ( shape -> stype )
  368.             {
  369.             case STYPE_PLANE:
  370.                 scale_y_plane(shape -> u.plane, factor);
  371.                 break;
  372.             case STYPE_SPHERE:
  373.                 if ( !sphere_to_quad(shape) )
  374.                     return ( FALSE );
  375.                 scale_y_quad(shape -> u.quad, factor);
  376.                 break;
  377.             case STYPE_QUAD:
  378.                 scale_y_quad(shape -> u.quad, factor);
  379.                 break;
  380.             }
  381.         scale_y_surf(shape -> surf, factor);
  382.         }
  383.     else
  384.         {
  385.         if ( !scale_y(shape -> u.shapes [0], factor) ||
  386.              !scale_y(shape -> u.shapes [1], factor) )
  387.             return ( FALSE );
  388.         }
  389.     return ( TRUE );
  390.     }
  391. /*...e*/
  392. /*...sscale_z             \45\ scale by factor in z direction:0:*/
  393. BOOLEAN    scale_z(SHAPE *shape, double factor)
  394.     {
  395.     if ( factor == 1.0 )
  396.         return ( TRUE );
  397.  
  398.     if ( shape -> stype <= STYPE_QUAD )
  399.         {
  400.         switch ( shape -> stype )
  401.             {
  402.             case STYPE_PLANE:
  403.                 scale_z_plane(shape -> u.plane, factor);
  404.                 break;
  405.             case STYPE_SPHERE:
  406.                 if ( !sphere_to_quad(shape) )
  407.                     return ( FALSE );
  408.                 scale_z_quad(shape -> u.quad, factor);
  409.                 break;
  410.             case STYPE_QUAD:
  411.                 scale_z_quad(shape -> u.quad, factor);
  412.                 break;
  413.             }
  414.         scale_z_surf(shape -> surf, factor);
  415.         }
  416.     else
  417.         {
  418.         if ( !scale_z(shape -> u.shapes [0], factor) ||
  419.              !scale_z(shape -> u.shapes [1], factor) )
  420.             return ( FALSE );
  421.         }
  422.     return ( TRUE );
  423.     }
  424. /*...e*/
  425. /*...sscale               \45\ scale by vector factor:0:*/
  426. BOOLEAN    scale(SHAPE *shape, VECTOR factor)
  427.     {
  428.     if ( shape -> stype == STYPE_SPHERE &&
  429.          factor.x == factor.y && factor.y == factor.z )
  430.         {
  431.         scale_sphere(shape -> u.sphere, factor.x);
  432.         scale_x_surf(shape -> surf, factor.x);
  433.         scale_y_surf(shape -> surf, factor.y);
  434.         scale_z_surf(shape -> surf, factor.z);
  435.         }
  436.     else
  437.         {
  438.         if ( !scale_x(shape, factor.x) ||
  439.              !scale_y(shape, factor.y) ||
  440.              !scale_z(shape, factor.z) )
  441.             return ( FALSE );
  442.         }
  443.     return ( TRUE );
  444.     }
  445. /*...e*/
  446. /*...srot_x               \45\ rotate about x axis by given angle:0:*/
  447. void    rot_x(SHAPE *shape, double angle)
  448.     {
  449.     if ( shape -> stype <= STYPE_QUAD )
  450.         {
  451.         switch ( shape -> stype )
  452.             {
  453.             case STYPE_PLANE:
  454.                 rot_x_plane(shape -> u.plane, angle);
  455.                 break;
  456.             case STYPE_SPHERE:
  457.                 rot_x_sphere(shape -> u.sphere, angle);
  458.                 break;
  459.             case STYPE_QUAD:
  460.                 rot_x_quad(shape -> u.quad, angle);
  461.                 break;
  462.             }
  463.         rot_x_surf(shape -> surf, angle);
  464.         }
  465.     else
  466.         {
  467.         rot_x(shape -> u.shapes [0], angle);
  468.         rot_x(shape -> u.shapes [1], angle);
  469.         }
  470.     }
  471. /*...e*/
  472. /*...srot_y               \45\ rotate about y axis by given angle:0:*/
  473. void    rot_y(SHAPE *shape, double angle)
  474.     {
  475.     if ( shape -> stype <= STYPE_QUAD )
  476.         {
  477.         switch ( shape -> stype )
  478.             {
  479.             case STYPE_PLANE:
  480.                 rot_y_plane(shape -> u.plane, angle);
  481.                 break;
  482.             case STYPE_SPHERE:
  483.                 rot_y_sphere(shape -> u.sphere, angle);
  484.                 break;
  485.             case STYPE_QUAD:
  486.                 rot_y_quad(shape -> u.quad, angle);
  487.                 break;
  488.             }
  489.         rot_y_surf(shape -> surf, angle);
  490.         }
  491.     else
  492.         {
  493.         rot_y(shape -> u.shapes [0], angle);
  494.         rot_y(shape -> u.shapes [1], angle);
  495.         }
  496.     }
  497. /*...e*/
  498. /*...srot_z               \45\ rotate about z axis by given angle:0:*/
  499. void    rot_z(SHAPE *shape, double angle)
  500.     {
  501.     if ( shape -> stype <= STYPE_QUAD )
  502.         {
  503.         switch ( shape -> stype )
  504.             {
  505.             case STYPE_PLANE:
  506.                 rot_z_plane(shape -> u.plane, angle);
  507.                 break;
  508.             case STYPE_SPHERE:
  509.                 rot_z_sphere(shape -> u.sphere, angle);
  510.                 break;
  511.             case STYPE_QUAD:
  512.                 rot_z_quad(shape -> u.quad, angle);
  513.                 break;
  514.             }
  515.         rot_z_surf(shape -> surf, angle);
  516.         }
  517.     else
  518.         {
  519.         rot_z(shape -> u.shapes [0], angle);
  520.         rot_z(shape -> u.shapes [1], angle);
  521.         }
  522.     }
  523. /*...e*/
  524.  
  525. /*...sresurf              \45\ change the surface of a shape:0:*/
  526. BOOLEAN resurf(SHAPE *shape, SURF *surf)
  527.     {
  528.     if ( shape -> stype <= STYPE_QUAD )
  529.         {
  530.         SURF *surf_copy;
  531.  
  532.         if ( (surf_copy = copy_surf(surf)) == NULL )
  533.             return ( FALSE );
  534.  
  535.         destroy_surf(shape -> surf);
  536.         shape -> surf = surf_copy;
  537.         }
  538.     else
  539.         {
  540.         if ( !resurf(shape -> u.shapes [0], surf) )
  541.             return ( FALSE );
  542.         if ( !resurf(shape -> u.shapes [1], surf) )
  543.             return ( FALSE );
  544.         }
  545.  
  546.     return ( TRUE );
  547.     }
  548. /*...e*/
  549.  
  550. /*...sis_empty_isectl     \45\ is intersection list empty:0:*/
  551. BOOLEAN    is_empty_isectl(ISECTL *il)
  552.     {
  553.     return ( il -> n_isects == 0 );
  554.     }
  555. /*...e*/
  556. /*...sis_solid_isectl     \45\ is intersection list solid:0:*/
  557. /*
  558. We will allow t from -INFINITE onwards or t from -INFINITE to INFINITE.
  559. */
  560.  
  561. BOOLEAN    is_solid_isectl(ISECTL *il)
  562.     {
  563.     if ( il -> n_isects > 2 )
  564.         return ( FALSE );
  565.  
  566.     if ( il -> isects [0].t != -INFINITE )
  567.         return ( FALSE );
  568.  
  569.     return ( il -> n_isects == 1 || il -> isects [1].t == INFINITE );
  570.     }
  571. /*...e*/
  572. /*...st_after_isectl      \45\ eliminate all intersections before a t value:0:*/
  573. void    t_after_isectl(ISECTL *il, double t)
  574.     {
  575.     int    i;
  576.  
  577.     for ( i = 0; i < il -> n_isects; i++ )
  578.         if ( il -> isects [i].t >= t )
  579.             break;
  580.  
  581.     if ( i == il -> n_isects )
  582.         /* All behind t */
  583.         il -> n_isects = 0;
  584.     else if ( i != 0 )
  585.         /* Some behind and some after t */
  586.         {
  587.         int    j = 0;
  588.  
  589.         if ( il -> isects [i].entering == FALSE )
  590.             /* Had better make an isect case first */
  591.             {
  592.             il -> isects [j  ]   = il -> isects [i - 1];
  593.             il -> isects [j++].t = t;
  594.             }
  595.  
  596.         while ( i < il -> n_isects )
  597.             il -> isects [j++] = il -> isects [i++];
  598.         il -> n_isects = j;
  599.         }
  600.     }
  601. /*...e*/
  602.  
  603. /*...screate_isectl       \45\ create an isectl:0:*/
  604. ISECTL    *create_isectl(int n_isects)
  605.     {
  606.     return ( malloc(sizeof(ISECTL) + (n_isects - 1) * sizeof(ISECT)) );
  607.     }
  608. /*...e*/
  609. /*...sdestroy_isectl      \45\ destroy an isectl:0:*/
  610. void    destroy_isectl(ISECTL *il)
  611.     {
  612.     free(il);
  613.     }
  614. /*...e*/
  615.  
  616. /*...spreprocess_shape    \45\ preprocess shape tree to save work when tracing:0:*/
  617. /*...sextent_shape:0:*/
  618. /*
  619. If we can find shapes that do not overlap, we can avoid tracing time later.
  620. This type of function is only made possible because the code that implements
  621. spheres etc. export their internal representations of the shapes.
  622. */
  623.  
  624. /*...sspheres_overlap:0:*/
  625. /*
  626. 2 spheres overlap if the sum of the radii > distance between the centres.
  627. We square both sides of this equation.
  628. */
  629.  
  630. static BOOLEAN spheres_overlap(SPHERE *sphere_a, SPHERE *sphere_b)
  631.     {
  632.     double    dx = sphere_a -> a - sphere_b -> a;
  633.     double    dy = sphere_a -> b - sphere_b -> b;
  634.     double    dz = sphere_a -> c - sphere_b -> c;
  635.     double    rr = sphere_a -> d + sphere_b -> d;
  636.  
  637.     return ( rr*rr >= dx*dx + dy*dy + dz*dz );
  638.     }
  639. /*...e*/
  640.  
  641. typedef struct { VECTOR min, max; } EXTENT;
  642.  
  643. /*...sextent_infinite:0:*/
  644. static EXTENT extent_infinite(void)
  645.     {
  646.     EXTENT    extent;
  647.  
  648.     extent.min.x = extent.min.y = extent.min.z = -INFINITE;
  649.     extent.max.x = extent.max.y = extent.max.z =  INFINITE;
  650.     return ( extent );
  651.     }
  652. /*...e*/
  653. /*...sextent_of_sphere:0:*/
  654. static EXTENT extent_of_sphere(SPHERE *sphere)
  655.     {
  656.     double    a = sphere -> a;
  657.     double    b = sphere -> b;
  658.     double    c = sphere -> c;
  659.     double    d = sphere -> d;
  660.     EXTENT    extent;
  661.  
  662.     extent.min.x = a - d; extent.max.x = a + d;
  663.     extent.min.y = b - d; extent.max.y = b + d;
  664.     extent.min.z = c - d; extent.max.z = c + d;
  665.  
  666.     return ( extent );
  667.     }
  668. /*...e*/
  669. /*...sextents_overlap:0:*/
  670. static BOOLEAN extents_overlap(EXTENT *extent_a, EXTENT *extent_b)
  671.     {
  672.     return ( extent_a -> max.x > extent_b -> min.x &&
  673.          extent_a -> max.y > extent_b -> min.y &&
  674.          extent_a -> max.z > extent_b -> min.z &&
  675.          extent_b -> max.x > extent_a -> min.x &&
  676.          extent_b -> max.y > extent_a -> min.y &&
  677.          extent_b -> max.z > extent_a -> min.z );
  678.     }
  679. /*...e*/
  680. /*...sextents_union:0:*/
  681. static EXTENT extents_union(EXTENT *extent_a, EXTENT *extent_b)
  682.     {
  683.     EXTENT extent;
  684.  
  685.     extent.min.x = min(extent_a -> min.x, extent_b -> min.x);
  686.     extent.min.y = min(extent_a -> min.y, extent_b -> min.y);
  687.     extent.min.z = min(extent_a -> min.z, extent_b -> min.z);
  688.     extent.max.x = max(extent_a -> max.x, extent_b -> max.x);
  689.     extent.max.y = max(extent_a -> max.y, extent_b -> max.y);
  690.     extent.max.z = max(extent_a -> max.z, extent_b -> max.z);
  691.  
  692.     return ( extent );
  693.     }
  694. /*...e*/
  695. /*...sextents_isect:0:*/
  696. static EXTENT extents_isect(EXTENT *extent_a, EXTENT *extent_b)
  697.     {
  698.     EXTENT    extent;
  699.  
  700.     extent.min.x = max(extent_a -> min.x, extent_b -> min.x);
  701.     extent.min.y = max(extent_a -> min.y, extent_b -> min.y);
  702.     extent.min.z = max(extent_a -> min.z, extent_b -> min.z);
  703.     extent.max.x = min(extent_a -> max.x, extent_b -> max.x);
  704.     extent.max.y = min(extent_a -> max.y, extent_b -> max.y);
  705.     extent.max.z = min(extent_a -> max.z, extent_b -> max.z);
  706.  
  707.     return ( extent );
  708.     }
  709. /*...e*/
  710.  
  711. static EXTENT extent_shape(SHAPE *shape)
  712.     {
  713.     static EXTENT dummy_extent = { { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0 } };
  714.  
  715.     if ( shape -> stype <= STYPE_QUAD )
  716.         {
  717.         switch ( shape -> stype )
  718.             {
  719. /*...sSTYPE_PLANE\44\ STYPE_QUAD \45\ too difficult:24:*/
  720. case STYPE_PLANE:
  721. case STYPE_QUAD:
  722.     return ( extent_infinite() );
  723. /*...e*/
  724. /*...sSTYPE_SPHERE:24:*/
  725. case STYPE_SPHERE:
  726.     return ( extent_of_sphere(shape -> u.sphere) );
  727. /*...e*/
  728.             }
  729.         }
  730.     else if ( shape -> stype == STYPE_EXTENT )
  731.         {
  732.         EXTENT extent;
  733.  
  734.         extent = extent_shape(shape -> u.shapes [0]);
  735.         /* Performed for side effect of extent_shape() call */
  736.  
  737.         return ( extent_shape(shape -> u.shapes [1]) );
  738.         }
  739.     else
  740.         {
  741.         SHAPE    *shape_a = shape -> u.shapes [0];
  742.         SHAPE    *shape_b = shape -> u.shapes [1];
  743.         EXTENT    extent_a;
  744.         EXTENT    extent_b;
  745.         EXTENT    extent;
  746.  
  747.         extent_a = extent_shape(shape_a);
  748.         extent_b = extent_shape(shape_b);
  749.  
  750.         switch ( shape -> stype )
  751.             {
  752. /*...sSTYPE_UNION\44\ STYPE_SDIFF:24:*/
  753. case STYPE_UNION:
  754. case STYPE_SDIFF:
  755.     extent = extents_union(&extent_a, &extent_b);
  756.     break;
  757. /*...e*/
  758. /*...sSTYPE_ISECT:24:*/
  759. case STYPE_ISECT:
  760.     extent = extents_isect(&extent_a, &extent_b);
  761.     break;
  762. /*...e*/
  763. /*...sSTYPE_DIFF:24:*/
  764. case STYPE_DIFF:
  765.     extent = extent_a;
  766.     break;
  767. /*...e*/
  768.             }
  769.  
  770.         if ( shape_a -> stype == STYPE_SPHERE &&
  771.              shape_b -> stype == STYPE_SPHERE )
  772.             shape -> overlap = spheres_overlap(
  773.                 shape_a -> u.sphere, shape_b -> u.sphere);
  774.         else
  775.             shape -> overlap = extents_overlap(&extent_a, &extent_b);
  776.  
  777.         return ( extent );
  778.         }
  779.  
  780.     return ( dummy_extent ); /* Keep fussy C compiler happy */
  781.     }
  782. /*...e*/
  783. /*...sn_isectls_reqd_shape:0:*/
  784. /*
  785. We use this function so that we can work out how many intersection lists
  786. we will need to pre-allocate for tracing the shape.
  787. */
  788.  
  789. int    n_isectls_reqd_shape(SHAPE *shape)
  790.     {
  791.     int nest0, nest1;
  792.  
  793.     if ( shape -> stype <= STYPE_QUAD )
  794.         return ( 1 );
  795.  
  796.     nest0 = n_isectls_reqd_shape(shape -> u.shapes [0]);
  797.     nest1 = n_isectls_reqd_shape(shape -> u.shapes [1]);
  798.  
  799.     return ( 2 + max(nest0, nest1) );
  800.     }
  801. /*...e*/
  802. /*...sisects_reqd_shape:0:*/
  803. /*
  804. Determine largest needed intersection list.
  805. */
  806.  
  807. int    isects_reqd_shape(SHAPE *shape)
  808.     {
  809.     if ( shape -> stype <= STYPE_QUAD )
  810.         switch ( shape -> stype )
  811.             {
  812. /*...sSTYPE_PLANE:24:*/
  813. case STYPE_PLANE:
  814.     return ( isects_reqd_plane(shape -> u.plane) );
  815. /*...e*/
  816. /*...sSTYPE_SPHERE:24:*/
  817. case STYPE_SPHERE:
  818.     return ( isects_reqd_sphere(shape -> u.sphere) );
  819. /*...e*/
  820. /*...sSTYPE_QUAD:24:*/
  821. case STYPE_QUAD:
  822.     return ( isects_reqd_quad(shape -> u.quad) );
  823. /*...e*/
  824.             }
  825.     else
  826.         {
  827.         int    reqd0 = isects_reqd_shape(shape -> u.shapes [0]);
  828.         int    reqd1 = isects_reqd_shape(shape -> u.shapes [1]);
  829.  
  830.         switch ( shape -> stype )
  831.             {
  832. /*...sSTYPE_UNION:24:*/
  833. case STYPE_UNION:
  834.     return ( reqd0 + reqd1 );
  835. /*...e*/
  836. /*...sSTYPE_ISECT:24:*/
  837. case STYPE_ISECT:
  838.     return ( ( shape -> overlap ) ? reqd0 + reqd1 : 0 );
  839. /*...e*/
  840. /*...sSTYPE_DIFF:24:*/
  841. case STYPE_DIFF:
  842.     return ( ( shape -> overlap ) ? reqd0 + reqd1 : reqd0 );
  843. /*...e*/
  844. /*...sSTYPE_SDIFF:24:*/
  845. case STYPE_SDIFF:
  846.     return ( reqd0 + reqd1 );
  847. /*...e*/
  848. /*...sSTYPE_EXTENT:24:*/
  849. case STYPE_EXTENT:
  850.     return ( max(reqd0, reqd1) );
  851. /*...e*/
  852.             }
  853.         }
  854.     return ( 0 ); /* Keep fussy C compiler happy */
  855.     }
  856. /*...e*/
  857.  
  858. void    preprocess_shape(SHAPE *shape, int *n_isectls, int *n_isects)
  859.     {
  860.     EXTENT    extent;
  861.  
  862.     extent = extent_shape(shape);
  863.         /* The side effect is to label the shape tree */
  864.  
  865.     *n_isectls = n_isectls_reqd_shape(shape);
  866.     *n_isects  = isects_reqd_shape(shape);
  867.     }
  868. /*...e*/
  869. /*...sintersect_shape     \45\ intersect with a shape to give a ISECTL:0:*/
  870. /*
  871. Given a shape tree, find the intersection list of a vector with it.
  872. If the shape tree is not a leaf node, then combine the intersection lists of
  873. the subtrees. We can make several optimisations in this area.
  874. */
  875.  
  876. /*...smake_empty_isectl:0:*/
  877. static void make_empty_isectl(ISECTL *il)
  878.     {
  879.     il -> n_isects = 0;
  880.     }
  881. /*...e*/
  882. /*...scopy_isectl:0:*/
  883. static void copy_isectl(ISECTL *il_dst, ISECTL *il_src)
  884.     {
  885.     int    i;
  886.  
  887.     il_dst -> n_isects = il_src -> n_isects;
  888.     for ( i = 0; i < il_dst -> n_isects; i++ )
  889.         il_dst -> isects [i] = il_src -> isects [i];
  890.     }
  891. /*...e*/
  892. /*...scombine_isectl:0:*/
  893. /*
  894. Given 2 intersection lists, produce a new one that is a combination of them.
  895.  
  896. The in_old flag is used to determine if the point of intersection changes
  897. meaning from in->out to out->in, or vice-versa. If this happens, the
  898. sense of the normal vector must change :-
  899.  
  900.          ..... .....                .....
  901.         .     .     .              .     .
  902.        .  A  . .  B  .            . A-B .
  903.       .     .   .     .          .     .
  904.     <-.   <-.   .->   .->      <-.     .-> This vector changes direction!
  905.       .     .   .     .          .     .
  906.        .     . .     .            .     .
  907.         .     .     .              .     .
  908.          ..... .....                .....
  909.  
  910. */
  911.  
  912. static void combine_isectl(
  913.     BOOLEAN (*combine)(BOOLEAN in_a, BOOLEAN in_b),
  914.     ISECTL *il_a, ISECTL *il_b, ISECTL *il)
  915.     {
  916.     BOOLEAN    in_a = FALSE;    /* When t = -INFINITE, in_a = FALSE */
  917.     BOOLEAN    in_b = FALSE;    /* When t = -INFINITE, in_b = FALSE */
  918.     BOOLEAN    in = FALSE;    /* Therefore combination starts as FALSE */
  919.     int    ptr_a = 0;
  920.     int    ptr_b = 0;
  921.     BOOLEAN    in_new, in_old;
  922.  
  923.     il -> n_isects = 0;
  924.  
  925.     /* Work through both a and b, looking at nearest ones first */
  926.  
  927.     while ( ptr_a < il_a -> n_isects && ptr_b < il_b -> n_isects )
  928.         {
  929.         ISECT    *isect;
  930.         double    t_a = il_a -> isects [ptr_a].t;
  931.         double    t_b = il_b -> isects [ptr_b].t;
  932.  
  933.         if ( t_a < t_b )
  934.             {
  935.             in_old = in_a = !in_a;
  936.             isect = &(il_a -> isects [ptr_a++]);
  937.             }
  938.         else if ( t_a > t_b )
  939.             {
  940.             in_old = in_b = !in_b;
  941.             isect = &(il_b -> isects [ptr_b++]);
  942.             }
  943.         else
  944.             /* Two surfaces at exactly the same place */
  945.             /* Not a very frequent event, but problematical */
  946.             /* Just label intersection arbitrarily as with B */
  947.             {
  948.             in_a = !in_a; ptr_a++;
  949.             in_old = in_b = !in_b;
  950.             isect = &(il_b -> isects [ptr_b++]);
  951.             }
  952.  
  953.         if ( (in_new = (*combine)(in_a, in_b)) != in )
  954.             /* Need to keep a record of this transition */
  955.             {
  956.             il -> isects [il -> n_isects] = *isect;
  957.             il -> isects [il -> n_isects].entering = in = in_new;
  958.             if ( in_new != in_old )
  959.                 il -> isects [il -> n_isects].negate_normal ^= TRUE;
  960.             il -> n_isects++;
  961.             }
  962.         }
  963.  
  964.     /* Either a or b is exhausted, so one of a or b may be left */
  965.  
  966.     while ( ptr_a < il_a -> n_isects )
  967.         {
  968.         in_old = in_a = !in_a;
  969.         if ( (in_new = (*combine)(in_a, in_b)) != in )
  970.             /* Need to keep a record of this transition */
  971.             {
  972.             il -> isects [il -> n_isects] = il_a -> isects [ptr_a];
  973.             il -> isects [il -> n_isects].entering = in = in_new;
  974.             if ( in_new != in_old )
  975.                 il -> isects [il -> n_isects].negate_normal ^= TRUE;
  976.             il -> n_isects++;
  977.             }
  978.         ptr_a++;
  979.         }
  980.  
  981.     while ( ptr_b < il_b -> n_isects )
  982.         {
  983.         in_old = in_b = !in_b;
  984.         if ( (in_new = (*combine)(in_a, in_b)) != in )
  985.             /* Need to keep a record of this transition */
  986.             {
  987.             il -> isects [il -> n_isects] = il_b -> isects [ptr_b];
  988.             il -> isects [il -> n_isects].entering = in = in_new;
  989.             if ( in_new != in_old )
  990.                 il -> isects [il -> n_isects].negate_normal ^= TRUE;
  991.             il -> n_isects++;
  992.             }
  993.         ptr_b++;
  994.         }
  995.     }
  996. /*...e*/
  997. /*...sunion_isectl:0:*/
  998. /*...scombine_union:0:*/
  999. static BOOLEAN combine_union(BOOLEAN in_a, BOOLEAN in_b)
  1000.     {
  1001.     return ( in_a || in_b );
  1002.     }
  1003. /*...e*/
  1004.  
  1005. static void union_isectl(ISECTL *il_a, ISECTL *il_b, ISECTL *il)
  1006.     {
  1007.     combine_isectl(combine_union, il_a, il_b, il);
  1008.     }
  1009. /*...e*/
  1010. /*...sisect_isectl:0:*/
  1011. /*...scombine_isect:0:*/
  1012. static BOOLEAN combine_isect(BOOLEAN in_a, BOOLEAN in_b)
  1013.     {
  1014.     return ( in_a && in_b );
  1015.     }
  1016. /*...e*/
  1017.  
  1018. static void isect_isectl(ISECTL *il_a, ISECTL *il_b, ISECTL *il)
  1019.     {
  1020.     combine_isectl(combine_isect, il_a, il_b, il);
  1021.     }
  1022. /*...e*/
  1023. /*...sdiff_isectl:0:*/
  1024. /*...scombine_diff:0:*/
  1025. static BOOLEAN combine_diff(BOOLEAN in_a, BOOLEAN in_b)
  1026.     {
  1027.     return ( in_a && !in_b );
  1028.     }
  1029. /*...e*/
  1030.  
  1031. static void diff_isectl(ISECTL *il_a, ISECTL *il_b, ISECTL *il)
  1032.     {
  1033.     combine_isectl(combine_diff, il_a, il_b, il);
  1034.     }
  1035. /*...e*/
  1036. /*...ssdiff_isectl:0:*/
  1037. /*...scombine_sdiff:0:*/
  1038. static BOOLEAN combine_sdiff(BOOLEAN in_a, BOOLEAN in_b)
  1039.     {
  1040.     return ( in_a ^ in_b );
  1041.     }
  1042. /*...e*/
  1043.  
  1044. static void sdiff_isectl(ISECTL *il_a, ISECTL *il_b, ISECTL *il)
  1045.     {
  1046.     combine_isectl(combine_sdiff, il_a, il_b, il);
  1047.     }
  1048. /*...e*/
  1049. /*...sconcat_isectl:0:*/
  1050. /*
  1051. This is a quick case of unioning two intersection lists for use when it is
  1052. known that they do not overlap.
  1053. */
  1054.  
  1055. static void concat_isectl(ISECTL *il_a, ISECTL *il_b, ISECTL *il)
  1056.     {
  1057.     ISECTL    *il_rest;
  1058.     int    i, j;
  1059.  
  1060.     if ( il_a -> n_isects == 0 )
  1061.         {
  1062.         copy_isectl(il, il_b);
  1063.         return;
  1064.         }
  1065.  
  1066.     if ( il_b -> n_isects == 0 )
  1067.         {
  1068.         copy_isectl(il, il_a);
  1069.         return;
  1070.         }
  1071.  
  1072.     if ( il_a -> isects [0].t < il_b -> isects [0].t )
  1073.         {
  1074.         copy_isectl(il, il_a);
  1075.         il_rest = il_b;
  1076.         }
  1077.     else
  1078.         {
  1079.         copy_isectl(il, il_b);
  1080.         il_rest = il_a;
  1081.         }
  1082.  
  1083.     for ( i = 0, j = il -> n_isects; i < il_rest -> n_isects; i++, j++ )
  1084.         il -> isects [j] = il_rest -> isects [i];
  1085.  
  1086.     il -> n_isects = j;
  1087.     }
  1088. /*...e*/
  1089.  
  1090. void    intersect_shape(SHAPE *shape, VECTOR p, VECTOR q, ISECTL *ils [])
  1091.     {
  1092.     if ( shape -> stype <= STYPE_QUAD )
  1093.         /* Is a leaf node */
  1094.         {
  1095.         SIL    sil;
  1096.         int    i;
  1097.     
  1098.         switch ( shape -> stype )
  1099.             {
  1100. /*...sSTYPE_PLANE:24:*/
  1101. case STYPE_PLANE:
  1102.     intersect_plane(shape -> u.plane, p, q, &sil);
  1103.     break;
  1104. /*...e*/
  1105. /*...sSTYPE_SPHERE:24:*/
  1106. case STYPE_SPHERE:
  1107.     intersect_sphere(shape -> u.sphere, p, q, &sil);
  1108.     break;
  1109. /*...e*/
  1110. /*...sSTYPE_QUAD:24:*/
  1111. case STYPE_QUAD:
  1112.     intersect_quad(shape -> u.quad, p, q, &sil);
  1113.     break;
  1114. /*...e*/
  1115.             }
  1116.         ils [0] -> n_isects = sil.n_sis;
  1117.         for ( i = 0; i < sil.n_sis; i++ )
  1118.             {
  1119.             ils [0] -> isects [i].t             = sil.sis [i].t;
  1120.             ils [0] -> isects [i].entering      = sil.sis [i].entering;
  1121.             ils [0] -> isects [i].shape         = shape;
  1122.             ils [0] -> isects [i].negate_normal = FALSE;
  1123.             }
  1124.         }
  1125.     else
  1126.         /* Binary combination of two subtrees */
  1127.         {
  1128.         SHAPE    *shape_a = shape -> u.shapes [0];
  1129.         SHAPE    *shape_b = shape -> u.shapes [1];
  1130.  
  1131.         switch ( shape -> stype )
  1132.             {
  1133. /*...sSTYPE_UNION:24:*/
  1134. case STYPE_UNION:
  1135.     if ( shape -> overlap )
  1136.         {
  1137.         intersect_shape(shape_b, p, q, ils + 1);
  1138.         if ( is_empty_isectl(ils [1]) )
  1139.             intersect_shape(shape_a, p, q, ils);
  1140.         else if ( is_solid_isectl(ils [1]) )
  1141.             copy_isectl(ils [0], ils [1]);
  1142.         else
  1143.             {
  1144.             intersect_shape(shape_a, p, q, ils + 2);
  1145.             union_isectl(ils [2], ils [1], ils [0]);
  1146.             }
  1147.         }
  1148.     else
  1149.         /* No overlap, treat like concatentation */
  1150.         {
  1151.         intersect_shape(shape_a, p, q, ils + 1);
  1152.         intersect_shape(shape_b, p, q, ils + 2);
  1153.         concat_isectl(ils [1], ils [2], ils [0]);
  1154.         }
  1155.     break;
  1156. /*...e*/
  1157. /*...sSTYPE_ISECT:24:*/
  1158. case STYPE_ISECT:
  1159.     if ( shape -> overlap )
  1160.         {
  1161.         intersect_shape(shape_b, p, q, ils + 1);
  1162.         if ( is_empty_isectl(ils [1]) )
  1163.             make_empty_isectl(ils [0]);
  1164.         else if ( is_solid_isectl(ils [1]) )
  1165.             intersect_shape(shape_a, p, q, ils);
  1166.         else
  1167.             {
  1168.             intersect_shape(shape_a, p, q, ils + 2);
  1169.             isect_isectl(ils [2], ils [1], ils [0]);
  1170.             }
  1171.         }
  1172.     else
  1173.         make_empty_isectl(ils [0]);
  1174.     break;
  1175. /*...e*/
  1176. /*...sSTYPE_DIFF:24:*/
  1177. case STYPE_DIFF:
  1178.     if ( shape -> overlap )
  1179.         {
  1180.         intersect_shape(shape_b, p, q, ils + 1);
  1181.         if ( is_empty_isectl(ils [1]) )
  1182.             intersect_shape(shape_a, p, q, ils);
  1183.         else if ( is_solid_isectl(ils [1]) )
  1184.             make_empty_isectl(ils [0]);
  1185.         else
  1186.             {
  1187.             intersect_shape(shape_a, p, q, ils + 2);
  1188.             diff_isectl(ils [2], ils [1], ils [0]);
  1189.             }
  1190.         }
  1191.     else
  1192.         intersect_shape(shape_a, p, q, ils);
  1193.     break;
  1194. /*...e*/
  1195. /*...sSTYPE_SDIFF:24:*/
  1196. case STYPE_SDIFF:
  1197.     if ( shape -> overlap )
  1198.         {
  1199.         intersect_shape(shape_b, p, q, ils + 1);
  1200.         if ( is_empty_isectl(ils [1]) )
  1201.             intersect_shape(shape_a, p, q, ils);
  1202.         else
  1203.             {
  1204.             intersect_shape(shape_a, p, q, ils + 2);
  1205.             sdiff_isectl(ils [2], ils [1], ils [0]);
  1206.             }
  1207.         }
  1208.     else
  1209.         /* No overlap, treat like concatentation */
  1210.         {
  1211.         intersect_shape(shape_a, p, q, ils + 1);
  1212.         intersect_shape(shape_b, p, q, ils + 2);
  1213.         concat_isectl(ils [1], ils [2], ils [0]);
  1214.         }
  1215.     break;
  1216. /*...e*/
  1217. /*...sSTYPE_EXTENT:24:*/
  1218. case STYPE_EXTENT:
  1219.     intersect_shape(shape_b, p, q, ils);
  1220.     if ( !is_empty_isectl(ils [0]) )
  1221.         intersect_shape(shape_a, p, q, ils);
  1222.     break;
  1223. /*...e*/
  1224.             }
  1225.         }
  1226.     }
  1227. /*...e*/
  1228. /*...snormal_to_shape     \45\ find normal at point on edge of shape:0:*/
  1229. VECTOR    normal_to_shape(SHAPE *shape, VECTOR p)
  1230.     {
  1231.     static VECTOR dummy_vector = { 0.0, 0.0, 0.0 };
  1232.  
  1233.     switch ( shape -> stype )
  1234.         {
  1235.         case STYPE_PLANE:
  1236.             return ( normal_to_plane(shape -> u.plane) );
  1237.         case STYPE_SPHERE:
  1238.             return ( normal_to_sphere(shape -> u.sphere, p) );
  1239.         case STYPE_QUAD:
  1240.             return ( normal_to_quad(shape -> u.quad, p) );
  1241.         }
  1242.     return ( dummy_vector ); /* Keep fussy C compiler happy */
  1243.     }
  1244. /*...e*/
  1245.