home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsv / ch3_5 / bsppartition.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-11-22  |  11.3 KB  |  347 lines
  1. /* partition.c: module to partition 3D convex face with an arbitrary plane.
  2.  * Copyright (c) Norman Chin
  3.  */ 
  4. #include "bsp.h"
  5.  
  6. #define PREPEND_FACE(f,fl) (f->fnext= fl, fl= f)
  7. #define ISVERTEX_EQ(v1,v2) \
  8.                           (IS_EQ((v1)->xx,(v2)->xx) && \
  9.                            IS_EQ((v1)->yy,(v2)->yy) && \
  10.                            IS_EQ((v1)->zz,(v2)->zz))
  11.  
  12. /* local functions */
  13. static VERTEX *findNextIntersection(VERTEX *vstart, const PLANE *plane,
  14.                     float *ixx, float *iyy, float *izz,
  15.                     SIGN *sign);
  16. static FACE *createOtherFace(FACE *face, VERTEX *v1, float ixx1, 
  17.                  float iyy1, float izz1, VERTEX *v2,
  18.                  float ixx2, float iyy2, float izz2
  19.                  );
  20. static SIGN whichSideIsFaceWRTplane(FACE *face, const PLANE *plane); 
  21.  
  22. /* Partitions a 3D convex polygon (face) with an arbitrary plane into its 
  23.  * negative and positive fragments, if any, w.r.t. the partitioning plane.
  24.  * Note that faceList is unusable afterwards since its vertex list has been
  25.  * parceled out to the other faces. It's set to null to avoid dangling
  26.  * pointer problem. Faces embedded in the plane are separated into two lists,
  27.  * one facing the same direction as the partitioning plane, faceSameDir, and 
  28.  * the other facing the opposite direction, faceOppDir.
  29.  */
  30. void BSPpartitionFaceListWithPlane(const PLANE *plane,FACE **faceList,
  31.                    FACE **faceNeg, FACE **facePos,
  32.                    FACE **faceSameDir, FACE **faceOppDir)
  33. {
  34. FACE *ftrav= *faceList;
  35.  
  36. *faceSameDir= *faceOppDir = *faceNeg= *facePos= NULL_FACE;
  37.  
  38. while (ftrav != NULL_FACE) {
  39.    VERTEX *v1, *v2; FACE *nextFtrav;
  40.    float ixx1,iyy1,izz1, ixx2,iyy2,izz2;
  41.    SIGN signV1, signV2;
  42.  
  43.    nextFtrav= ftrav->fnext;    /* unchain current face from list */
  44.    ftrav->fnext= NULL_FACE;
  45.  
  46.    /* find first intersection */
  47.    v1= findNextIntersection(ftrav->vhead,plane,&ixx1,&iyy1,&izz1,&signV1);
  48.    if (v1 != NULL_VERTEX) {
  49.       assert(signV1 != ZERO);
  50.  
  51.       /* first one found, find the 2nd one, if any */
  52.       v2= findNextIntersection(v1->vnext,plane,&ixx2,&iyy2,&izz2,&signV2);
  53.  
  54.       /* Due to numerical instabilities, both intersection points may
  55.        * have the same sign such as in the case when splitting very close
  56.        * to a vertex. This should not count as a split.
  57.        */
  58.       if (v2 != NULL_VERTEX && signV1 == signV2) v2= NULL_VERTEX; 
  59.  
  60.    }
  61.    else v2= NULL_VERTEX;    /* No first intersection found,
  62.                                  * therefore no second intersection.
  63.                  */
  64.                                 
  65.    /* an intersection? */
  66.    if (v1 != NULL_VERTEX && v2 != NULL_VERTEX) { /* yes, intersection */
  67.       FACE *newOtherFace;
  68.  
  69.       /* ftrav's vertex list will be modified */
  70.       newOtherFace= createOtherFace(ftrav,v1,ixx1,iyy1,izz1,
  71.                     v2,ixx2,iyy2,izz2
  72.                     );
  73.  
  74.       /* return split faces on appropriate lists */
  75.       if (signV1 == NEGATIVE) { 
  76.      PREPEND_FACE(ftrav,*faceNeg);
  77.      PREPEND_FACE(newOtherFace,*facePos);
  78.       }
  79.       else {
  80.      assert(signV1 == POSITIVE);
  81.      PREPEND_FACE(newOtherFace,*faceNeg);
  82.      PREPEND_FACE(ftrav,*facePos);
  83.       }
  84.                     
  85.    }
  86.    else {            /* no intersection  */
  87.       SIGN side;
  88.  
  89.       /* Face is embedded or wholly to one side of partitioning plane. */
  90.       side= whichSideIsFaceWRTplane(ftrav,plane);
  91.       if (side == NEGATIVE) 
  92.      PREPEND_FACE(ftrav,*faceNeg);
  93.       else if (side == POSITIVE) 
  94.      PREPEND_FACE(ftrav,*facePos);
  95.       else {
  96.      assert(side == ZERO);
  97.      if (IS_EQ(ftrav->plane.aa,plane->aa) && 
  98.          IS_EQ(ftrav->plane.bb,plane->bb) &&
  99.          IS_EQ(ftrav->plane.cc,plane->cc)) 
  100.             PREPEND_FACE(ftrav,*faceSameDir);
  101.      else PREPEND_FACE(ftrav,*faceOppDir);
  102.       }
  103.    }
  104.    ftrav= nextFtrav;        /* get next */
  105. } /* while ftrav != NULL_FACE */
  106.  
  107.    /* faceList's vertex list has been parceled out to other lists so
  108.     * set this to null.
  109.     */
  110.    *faceList= NULL_FACE;        
  111. #if 0 /* only true for BSPconstructTree() */
  112.    /* there's at least one face embedded in the partitioning plane */
  113.    assert(*faceSameDir != NULL_FACE); 
  114. #endif
  115. } /* BSPpartitionFaceListWithPlane() */
  116.                 
  117. /* Finds next intersection on or after vstart. 
  118.  * 
  119.  * If an intersection is found, 
  120.  *    a pointer to first vertex of the edge is returned, 
  121.  *    the intersection point (ixx,iyy,izz) and its sign is updated. 
  122.  * Otherwise a null pointer is returned.
  123.  */
  124. static VERTEX *findNextIntersection(VERTEX *vstart,const PLANE *plane,
  125.                     float *ixx,float *iyy,float *izz,
  126.                     SIGN *sign)
  127. {
  128.    VERTEX *vtrav;
  129.  
  130.    /* for all edges starting from vstart ... */
  131.    for (vtrav= vstart; vtrav->vnext != NULL_VERTEX; vtrav= vtrav->vnext) {
  132.       if ((*sign= anyEdgeIntersectWithPlane(vtrav->xx,vtrav->yy,vtrav->zz,
  133.                         vtrav->vnext->xx,vtrav->vnext->yy,
  134.                         vtrav->vnext->zz,plane,
  135.                         ixx,iyy,izz))) {
  136.      return(vtrav);
  137.       }
  138.    }
  139.  
  140.    return(NULL_VERTEX);
  141. } /* findNextIntersection() */
  142.  
  143. /* Memory allocated for split face's vertices and pointers tediously updated.
  144.  *
  145.  * face - face to be split
  146.  * v1   - 1st vertex of edge of where 1st intersection was found 
  147.  * (ixx1,iyy1,izz1) - 1st intersection
  148.  * v2   - 1st vertex of edge of where 2nd intersection was found 
  149.  * (ixx2,iyy2,izz2) - 2nd intersection
  150.  */
  151. static FACE *createOtherFace(FACE *face, 
  152.                  VERTEX *v1, float ixx1, float iyy1, float izz1,
  153.                  VERTEX *v2, float ixx2, float iyy2, float izz2
  154.                  )
  155. {
  156.    VERTEX *i1p1, *i2p1;        /* new vertices for original face  */
  157.    VERTEX *i1p2, *i2p2;        /* new vertices for new face */
  158.    VERTEX *p2end;        /* new vertex for end of new face */
  159.    VERTEX *vtemp;        /* place holders */
  160.    register VERTEX *beforeV2;    /* place holders */
  161.    FACE *newFace;        /* newly allocated face */
  162.  
  163.    /* new intersection vertices */
  164.    i1p1= allocVertex(ixx1,iyy1,izz1); 
  165.    i2p1= allocVertex(ixx2,iyy2,izz2);
  166.    i1p2= allocVertex(ixx1,iyy1,izz1);
  167.    i2p2= allocVertex(ixx2,iyy2,izz2); 
  168.  
  169.    /* duplicate 1st vertex of 2nd list to close it up */
  170.    p2end= allocVertex(v2->xx,v2->yy,v2->zz);
  171.  
  172.    vtemp= v1->vnext;
  173.  
  174.    /* merge both intersection vertices i1p1 & i2p1 into 1st list */
  175.    if (ISVERTEX_EQ(i2p1,v2->vnext)) { /* intersection vertex coincident? */
  176.       assert(i2p1->vnext == NULL_VERTEX);
  177.       freeVertexList(&i2p1);    /* yes, so free it */
  178.       i1p1->vnext= v2->vnext;
  179.    }
  180.    else {
  181.       i2p1->vnext= v2->vnext;    /* attach intersection list onto 1st list */
  182.       i1p1->vnext= i2p1;    /* attach both intersection vertices */
  183.    }
  184.    v1->vnext= i1p1; /* attach front of 1st list to intersection vertices */
  185.  
  186.    /* merge intersection vertices i1p2, i2p2 & p2end into second list */
  187.    i2p2->vnext= i1p2;        /* attach both intersection vertices */
  188.    v2->vnext= i2p2;        /* attach 2nd list to interection vertices */
  189.    if (vtemp == v2) {
  190.       i1p2->vnext= p2end;    /* close up 2nd list */
  191.    }
  192.    else {
  193.       if (ISVERTEX_EQ(i1p2,vtemp)) { /* intersection vertex coincident? */
  194.      assert(i1p2->vnext == NULL_VERTEX);
  195.      freeVertexList(&i1p2);    /* yes, so free it */
  196.      i2p2->vnext= vtemp;    /* attach intersection vertex to 2nd list */
  197.       }
  198.       else {
  199.      i1p2->vnext= vtemp;    /* attach intersection list to 2nd list */
  200.       }
  201.       /* find previous vertex to v2 */
  202.       for (beforeV2= vtemp; beforeV2->vnext != v2; beforeV2= beforeV2->vnext)
  203.      ;            /* lone semi-colon */
  204.       beforeV2->vnext= p2end;    /* and attach it to complete the 2nd list */
  205.    }
  206.  
  207.    /* copy original face info but with new vertex list */
  208.    newFace= allocFace(&face->color,v2,&face->plane);
  209.  
  210.    return(newFace);
  211. } /* createOtherFace() */
  212.  
  213. /* Determines which side a face is with respect to a plane.
  214.  *
  215.  * However, due to numerical problems, when a face is very close to the plane,
  216.  * some vertices may be misclassified. 
  217.  * There are several solutions, two of which are mentioned here:
  218.  *    1- classify the one vertex furthest away from the plane, (note that
  219.  *       one need not compute the actual distance) and use that side.
  220.  *    2- count how many vertices lie on either side and pick the side
  221.  *       with the maximum. (this is the one implemented).
  222.  */
  223. static SIGN whichSideIsFaceWRTplane(FACE *face, const PLANE *plane)
  224. {
  225.    register VERTEX *vtrav;
  226.    float value;
  227.    boolean isNeg, isPos;
  228.  
  229.    isNeg= isPos= FALSE;
  230.    
  231.    for (vtrav= face->vhead; vtrav->vnext != NULL_VERTEX; vtrav= vtrav->vnext){
  232.       value= (plane->aa*vtrav->xx) + (plane->bb*vtrav->yy) + 
  233.          (plane->cc*vtrav->zz) + plane->dd;
  234.       if (value < -TOLER) 
  235.      isNeg= TRUE;
  236.       else if (value > TOLER)
  237.      isPos= TRUE;
  238.       else assert(IS_EQ(value,0.0));
  239.    } /* for vtrav */ 
  240.  
  241.    /* in the very rare case that some vertices slipped thru to other side of
  242.     * plane due to round-off errors, execute the above again but count the 
  243.     * vertices on each side instead and pick the maximum.
  244.     */
  245.    if (isNeg && isPos) {    /* yes so handle this numerical problem */
  246.       int countNeg, countPos;
  247.       
  248.       /* count how many vertices are on either side */
  249.       countNeg= countPos= 0;
  250.       for (vtrav= face->vhead; vtrav->vnext != NULL_VERTEX; 
  251.        vtrav= vtrav->vnext) {
  252.      value= (plane->aa*vtrav->xx) + (plane->bb*vtrav->yy) + 
  253.             (plane->cc*vtrav->zz) + plane->dd;
  254.      if (value < -TOLER)
  255.         countNeg++;
  256.      else if (value > TOLER)
  257.         countPos++;
  258.      else assert(IS_EQ(value,0.0));
  259.       } /* for */
  260.  
  261.       /* return the side corresponding to the maximum */
  262.       if (countNeg > countPos) return(NEGATIVE);
  263.       else if (countPos > countNeg) return(POSITIVE);
  264.       else return(ZERO);
  265.    }
  266.    else {            /* this is the usual case */
  267.       if (isNeg) return(NEGATIVE);
  268.       else if (isPos) return(POSITIVE);
  269.       else return(ZERO);
  270.    }
  271. } /* whichSideIsFaceWRTplane() */
  272.  
  273. /* Determines if an edge bounded by (x1,y1,z1)->(x2,y2,z2) intersects
  274.  * the plane.
  275.  * 
  276.  * If there's an intersection, 
  277.  *    the sign of (x1,y1,z1), NEGATIVE or POSITIVE, w.r.t. the plane is
  278.  *    returned with the intersection (ixx,iyy,izz) updated.
  279.  * Otherwise ZERO is returned.    
  280.  */
  281. SIGN anyEdgeIntersectWithPlane(float x1,float y1,float z1,
  282.                    float x2,float y2,float z2,
  283.                    const PLANE *plane,
  284.                    float *ixx, float *iyy, float *izz)
  285. {
  286.    float temp1, temp2;
  287.    int sign1, sign2;        /* must be int since gonna do a bitwise ^ */
  288.    float aa= plane->aa; float bb= plane->bb; float cc= plane->cc;
  289.    float dd= plane->dd;
  290.  
  291.    /* get signs */
  292.    temp1= (aa*x1) + (bb*y1) + (cc*z1) + dd;
  293.    if (temp1 < -TOLER)
  294.       sign1= -1;
  295.    else if (temp1 > TOLER)
  296.       sign1= 1;
  297.    else {
  298.       /* edges beginning with a 0 sign are not considered valid intersections
  299.        * case 1 & 6 & 7, see Gems III.
  300.        */
  301.       assert(IS_EQ(temp1,0.0));
  302.       return(ZERO);
  303.    }
  304.  
  305.    temp2= (aa*x2) + (bb*y2) + (cc*z2) + dd;
  306.    if (temp2 < -TOLER)
  307.       sign2= -1;
  308.    else if (temp2 > TOLER)
  309.       sign2= 1;
  310.    else {            /* case 8 & 9, see Gems III */
  311.       assert(IS_EQ(temp2,0.0));
  312.       *ixx= x2;
  313.       *iyy= y2;
  314.       *izz= z2;
  315.  
  316.       return( (sign1 == -1) ? NEGATIVE : POSITIVE);
  317.    }
  318.  
  319.    /* signs different? 
  320.     * recall: -1^1 == 1^-1 ==> 1    case 4 & 5, see Gems III
  321.     *         -1^-1 == 1^1 ==> 0    case 2 & 3, see Gems III
  322.     */
  323.    if (sign1 ^ sign2) {
  324.       float dx,dy,dz;
  325.       float denom, tt;
  326.  
  327.       /* compute intersection point */
  328.       dx= x2-x1;
  329.       dy= y2-y1;
  330.       dz= z2-z1;
  331.  
  332.       denom= (aa*dx) + (bb*dy) + (cc*dz);
  333.       tt= - ((aa*x1) + (bb*y1) + (cc*z1) + dd) / denom;
  334.  
  335.       *ixx= x1 + (tt * dx);
  336.       *iyy= y1 + (tt * dy);
  337.       *izz= z1 + (tt * dz);
  338.  
  339.       assert(sign1 == 1 || sign1 == -1);
  340.  
  341.       return( (sign1 == -1) ? NEGATIVE : POSITIVE );
  342.    }
  343.    else return(ZERO);
  344.  
  345. } /* anyEdgeIntersectWithPlane() */
  346. /*** bspPartition.c ***/
  347.