Geant4 LXR

Cross-Referencing   Geant4
Geant4/externals/g4tools/include/tools/glutess/normal

   // see license file for original license.
   
   #ifndef tools_glutess_normal
   #define tools_glutess_normal
   
   #include "_tess"
   
   /* __gl_projectPolygon( tess ) determines the polygon normal
    * and project vertices onto the plane of the polygon.
   */
  //void __gl_projectPolygon( GLUtesselator *tess );
  
  ////////////////////////////////////////////////////////
  /// inlined C code : ///////////////////////////////////
  ////////////////////////////////////////////////////////
  #include <cmath>
  
  #define Dot(u,v)  (u[0]*v[0] + u[1]*v[1] + u[2]*v[2])
  
  inline/*static*/ int static_LongAxis( GLUdouble v[3] )
  {
    int i = 0;
  
    if( GLU_ABS(v[1]) > GLU_ABS(v[0]) ) { i = 1; }
    if( GLU_ABS(v[2]) > GLU_ABS(v[i]) ) { i = 2; }
    return i;
  }
  
  inline/*static*/ void static_ComputeNormal( GLUtesselator *tess, GLUdouble norm[3] )
  {
    GLUvertex *v, *v1, *v2;
    GLUdouble c, tLen2, maxLen2;
    GLUdouble maxVal[3], minVal[3], d1[3], d2[3], tNorm[3];
    GLUvertex *maxVert[3], *minVert[3];
    GLUvertex *vHead = &tess->mesh->vHead;
    int i;
  
    maxVal[0] = maxVal[1] = maxVal[2] = -2 * GLU_TESS_MAX_COORD;
    minVal[0] = minVal[1] = minVal[2] = 2 * GLU_TESS_MAX_COORD;
    
    minVert[0] = 0;minVert[1] = 0;minVert[2] = 0; //G.Barrand : to quiet Coverity.
    maxVert[0] = 0;maxVert[1] = 0;maxVert[2] = 0; //G.Barrand : to quiet Coverity.
  
    for( v = vHead->next; v != vHead; v = v->next ) {
      for( i = 0; i < 3; ++i ) {
        c = v->coords[i];
        if( c < minVal[i] ) { minVal[i] = c; minVert[i] = v; }
        if( c > maxVal[i] ) { maxVal[i] = c; maxVert[i] = v; }
      }
    }
  
    /* Find two vertices separated by at least 1/sqrt(3) of the maximum
     * distance between any two vertices
     */
    i = 0;
    if( maxVal[1] - minVal[1] > maxVal[0] - minVal[0] ) { i = 1; }
    if( maxVal[2] - minVal[2] > maxVal[i] - minVal[i] ) { i = 2; }
    if( minVal[i] >= maxVal[i] ) {
      /* All vertices are the same -- normal doesn't matter */
      norm[0] = 0; norm[1] = 0; norm[2] = 1;
      return;
    }
  
    /* Look for a third vertex which forms the triangle with maximum area
     * (Length of normal == twice the triangle area)
     */
    maxLen2 = 0;
    v1 = minVert[i];
    v2 = maxVert[i];
    if( !v1 || !v2 ) {norm[0] = 0; norm[1] = 0; norm[2] = 1;return;} //G.Barrand.
    d1[0] = v1->coords[0] - v2->coords[0];
    d1[1] = v1->coords[1] - v2->coords[1];
    d1[2] = v1->coords[2] - v2->coords[2];
    for( v = vHead->next; v != vHead; v = v->next ) {
      d2[0] = v->coords[0] - v2->coords[0];
      d2[1] = v->coords[1] - v2->coords[1];
      d2[2] = v->coords[2] - v2->coords[2];
      tNorm[0] = d1[1]*d2[2] - d1[2]*d2[1];
      tNorm[1] = d1[2]*d2[0] - d1[0]*d2[2];
      tNorm[2] = d1[0]*d2[1] - d1[1]*d2[0];
      tLen2 = tNorm[0]*tNorm[0] + tNorm[1]*tNorm[1] + tNorm[2]*tNorm[2];
      if( tLen2 > maxLen2 ) {
        maxLen2 = tLen2;
        norm[0] = tNorm[0];
        norm[1] = tNorm[1];
        norm[2] = tNorm[2];
      }
    }
  
    if( maxLen2 <= 0 ) {
      /* All points lie on a single line -- any decent normal will do */
      norm[0] = norm[1] = norm[2] = 0;
      norm[static_LongAxis(d1)] = 1;
    }
  }
  
  
  inline/*static*/ void static_CheckOrientation( GLUtesselator *tess )
  {
   GLUdouble area;
   GLUface *f, *fHead = &tess->mesh->fHead;
   GLUvertex *v, *vHead = &tess->mesh->vHead;
   GLUhalfEdge *e;
 
   /* When we compute the normal automatically, we choose the orientation
    * so that the sum of the signed areas of all contours is non-negative.
    */
   area = 0;
   for( f = fHead->next; f != fHead; f = f->next ) {
     e = f->anEdge;
     if( e->winding <= 0 ) continue;
     do {
       area += (e->Org->s - e->Dst->s) * (e->Org->t + e->Dst->t);
       e = e->Lnext;
     } while( e != f->anEdge );
   }
   if( area < 0 ) {
     /* Reverse the orientation by flipping all the t-coordinates */
     for( v = vHead->next; v != vHead; v = v->next ) {
       v->t = - v->t;
     }
     tess->tUnit[0] = - tess->tUnit[0];
     tess->tUnit[1] = - tess->tUnit[1];
     tess->tUnit[2] = - tess->tUnit[2];
   }
 }
 
 #if defined(SLANTED_SWEEP)
 /* The "feature merging" is not intended to be complete.  There are
  * special cases where edges are nearly parallel to the sweep line
  * which are not implemented.  The algorithm should still behave
  * robustly (ie. produce a reasonable tesselation) in the presence
  * of such edges, however it may miss features which could have been
  * merged.  We could minimize this effect by choosing the sweep line
  * direction to be something unusual (ie. not parallel to one of the
  * coordinate axes).
  */
 #define S_UNIT_X  0.50941539564955385 /* Pre-normalized */
 #define S_UNIT_Y  0.86052074622010633
 #else
 #define S_UNIT_X  1.0
 #define S_UNIT_Y  0.0
 #endif
 
 /* Determine the polygon normal and project vertices onto the plane
  * of the polygon.
  */
 inline void __gl_projectPolygon( GLUtesselator *tess )
 {
   GLUvertex *v, *vHead = &tess->mesh->vHead;
   GLUdouble norm[3];
   GLUdouble *sUnit, *tUnit;
   int i, computedNormal = TOOLS_GLU_FALSE;
 
   norm[0] = tess->normal[0];
   norm[1] = tess->normal[1];
   norm[2] = tess->normal[2];
   if( norm[0] == 0 && norm[1] == 0 && norm[2] == 0 ) {
     static_ComputeNormal( tess, norm );
     computedNormal = TOOLS_GLU_TRUE;
   }
   sUnit = tess->sUnit;
   tUnit = tess->tUnit;
   i = static_LongAxis( norm );
 
 #if defined(FOR_TRITE_TEST_PROGRAM) || defined(TRUE_PROJECT)
   /* Choose the initial sUnit vector to be approximately perpendicular
    * to the normal.
    */
   Normalize( norm );
 
   sUnit[i] = 0;
   sUnit[(i+1)%3] = S_UNIT_X;
   sUnit[(i+2)%3] = S_UNIT_Y;
 
   /* Now make it exactly perpendicular */
   w = Dot( sUnit, norm );
   sUnit[0] -= w * norm[0];
   sUnit[1] -= w * norm[1];
   sUnit[2] -= w * norm[2];
   Normalize( sUnit );
 
   /* Choose tUnit so that (sUnit,tUnit,norm) form a right-handed frame */
   tUnit[0] = norm[1]*sUnit[2] - norm[2]*sUnit[1];
   tUnit[1] = norm[2]*sUnit[0] - norm[0]*sUnit[2];
   tUnit[2] = norm[0]*sUnit[1] - norm[1]*sUnit[0];
   Normalize( tUnit );
 #else
   /* Project perpendicular to a coordinate axis -- better numerically */
   sUnit[i] = 0;
   sUnit[(i+1)%3] = S_UNIT_X;
   sUnit[(i+2)%3] = S_UNIT_Y;
 
   tUnit[i] = 0;
   tUnit[(i+1)%3] = (norm[i] > 0) ? -S_UNIT_Y : S_UNIT_Y;
   tUnit[(i+2)%3] = (norm[i] > 0) ? S_UNIT_X : -S_UNIT_X;
 #endif
 
   /* Project the vertices onto the sweep plane */
   for( v = vHead->next; v != vHead; v = v->next ) {
     v->s = Dot( v->coords, sUnit );
     v->t = Dot( v->coords, tUnit );
   }
   if( computedNormal ) {
     static_CheckOrientation( tess );
   }
 }
 
 #endif