src/glu/sgi/libtess/normal.c

   1 /*
   2  * SGI FREE SOFTWARE LICENSE B (Version 2.0, Sept. 18, 2008)
   3  * Copyright (C) 1991-2000 Silicon Graphics, Inc. All Rights Reserved.
   4  *
   5  * Permission is hereby granted, free of charge, to any person obtaining a
   6  * copy of this software and associated documentation files (the "Software"),
   7  * to deal in the Software without restriction, including without limitation
   8  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
   9  * and/or sell copies of the Software, and to permit persons to whom the
  10  * Software is furnished to do so, subject to the following conditions:
  11  *
  12  * The above copyright notice including the dates of first publication and
  13  * either this permission notice or a reference to
  14  * http://oss.sgi.com/projects/FreeB/
  15  * shall be included in all copies or substantial portions of the Software.
  16  *
  17  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  18  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  19  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  20  * SILICON GRAPHICS, INC. BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
  21  * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF
  22  * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  23  * SOFTWARE.
  24  *
  25  * Except as contained in this notice, the name of Silicon Graphics, Inc.
  26  * shall not be used in advertising or otherwise to promote the sale, use or
  27  * other dealings in this Software without prior written authorization from
  28  * Silicon Graphics, Inc.
  29  */
  30 /*
  31 ** Author: Eric Veach, July 1994.
  32 **
  33 */
  34
  35 #include "gluos.h"
  36 #include "mesh.h"
  37 #include "tess.h"
  38 #include "normal.h"
  39 #include <math.h>
  40 #include <assert.h>
  41
  42 #define TRUE 1
  43 #define FALSE 0
  44
  45 #define Dot(u,v)        (u[0]*v[0] + u[1]*v[1] + u[2]*v[2])
  46
  47 #if 0
  48 static void Normalize( GLdouble v[3] )
  49 {
  50   GLdouble len = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
  51
  52   assert( len > 0 );
  53   len = sqrt( len );
  54   v[0] /= len;
  55   v[1] /= len;
  56   v[2] /= len;
  57 }
  58 #endif
  59
  60 #undef  ABS
  61 #define ABS(x)  ((x) < 0 ? -(x) : (x))
  62
  63 static int LongAxis( GLdouble v[3] )
  64 {
  65   int i = 0;
  66
  67   if( ABS(v[1]) > ABS(v[0]) ) { i = 1; }
  68   if( ABS(v[2]) > ABS(v[i]) ) { i = 2; }
  69   return i;
  70 }
  71
  72 static void ComputeNormal( GLUtesselator *tess, GLdouble norm[3] )
  73 {
  74   GLUvertex *v, *v1, *v2;
  75   GLdouble c, tLen2, maxLen2;
  76   GLdouble maxVal[3], minVal[3], d1[3], d2[3], tNorm[3];
  77   GLUvertex *maxVert[3], *minVert[3];
  78   GLUvertex *vHead = &tess->mesh->vHead;
  79   int i;
  80
  81   maxVal[0] = maxVal[1] = maxVal[2] = -2 * GLU_TESS_MAX_COORD;
  82   minVal[0] = minVal[1] = minVal[2] = 2 * GLU_TESS_MAX_COORD;
  83
  84   for( v = vHead->next; v != vHead; v = v->next ) {
  85     for( i = 0; i < 3; ++i ) {
  86       c = v->coords[i];
  87       if( c < minVal[i] ) { minVal[i] = c; minVert[i] = v; }
  88       if( c > maxVal[i] ) { maxVal[i] = c; maxVert[i] = v; }
  89     }
  90   }
  91
  92   /* Find two vertices separated by at least 1/sqrt(3) of the maximum
  93    * distance between any two vertices
  94    */
  95   i = 0;
  96   if( maxVal[1] - minVal[1] > maxVal[0] - minVal[0] ) { i = 1; }
  97   if( maxVal[2] - minVal[2] > maxVal[i] - minVal[i] ) { i = 2; }
  98   if( minVal[i] >= maxVal[i] ) {
  99     /* All vertices are the same -- normal doesn't matter */
 100     norm[0] = 0; norm[1] = 0; norm[2] = 1;
 101     return;
 102   }
 103
 104   /* Look for a third vertex which forms the triangle with maximum area
 105    * (Length of normal == twice the triangle area)
 106    */
 107   maxLen2 = 0;
 108   v1 = minVert[i];
 109   v2 = maxVert[i];
 110   d1[0] = v1->coords[0] - v2->coords[0];
 111   d1[1] = v1->coords[1] - v2->coords[1];
 112   d1[2] = v1->coords[2] - v2->coords[2];
 113   for( v = vHead->next; v != vHead; v = v->next ) {
 114     d2[0] = v->coords[0] - v2->coords[0];
 115     d2[1] = v->coords[1] - v2->coords[1];
 116     d2[2] = v->coords[2] - v2->coords[2];
 117     tNorm[0] = d1[1]*d2[2] - d1[2]*d2[1];
 118     tNorm[1] = d1[2]*d2[0] - d1[0]*d2[2];
 119     tNorm[2] = d1[0]*d2[1] - d1[1]*d2[0];
 120     tLen2 = tNorm[0]*tNorm[0] + tNorm[1]*tNorm[1] + tNorm[2]*tNorm[2];
 121     if( tLen2 > maxLen2 ) {
 122       maxLen2 = tLen2;
 123       norm[0] = tNorm[0];
 124       norm[1] = tNorm[1];
 125       norm[2] = tNorm[2];
 126     }
 127   }
 128
 129   if( maxLen2 <= 0 ) {
 130     /* All points lie on a single line -- any decent normal will do */
 131     norm[0] = norm[1] = norm[2] = 0;
 132     norm[LongAxis(d1)] = 1;
 133   }
 134 }
 135
 136
 137 static void CheckOrientation( GLUtesselator *tess )
 138 {
 139   GLdouble area;
 140   GLUface *f, *fHead = &tess->mesh->fHead;
 141   GLUvertex *v, *vHead = &tess->mesh->vHead;
 142   GLUhalfEdge *e;
 143
 144   /* When we compute the normal automatically, we choose the orientation
 145    * so that the the sum of the signed areas of all contours is non-negative.
 146    */
 147   area = 0;
 148   for( f = fHead->next; f != fHead; f = f->next ) {
 149     e = f->anEdge;
 150     if( e->winding <= 0 ) continue;
 151     do {
 152       area += (e->Org->s - e->Dst->s) * (e->Org->t + e->Dst->t);
 153       e = e->Lnext;
 154     } while( e != f->anEdge );
 155   }
 156   if( area < 0 ) {
 157     /* Reverse the orientation by flipping all the t-coordinates */
 158     for( v = vHead->next; v != vHead; v = v->next ) {
 159       v->t = - v->t;
 160     }
 161     tess->tUnit[0] = - tess->tUnit[0];
 162     tess->tUnit[1] = - tess->tUnit[1];
 163     tess->tUnit[2] = - tess->tUnit[2];
 164   }
 165 }
 166
 167 #ifdef FOR_TRITE_TEST_PROGRAM
 168 #include <stdlib.h>
 169 extern int RandomSweep;
 170 #define S_UNIT_X        (RandomSweep ? (2*drand48()-1) : 1.0)
 171 #define S_UNIT_Y        (RandomSweep ? (2*drand48()-1) : 0.0)
 172 #else
 173 #if defined(SLANTED_SWEEP)
 174 /* The "feature merging" is not intended to be complete.  There are
 175  * special cases where edges are nearly parallel to the sweep line
 176  * which are not implemented.  The algorithm should still behave
 177  * robustly (ie. produce a reasonable tesselation) in the presence
 178  * of such edges, however it may miss features which could have been
 179  * merged.  We could minimize this effect by choosing the sweep line
 180  * direction to be something unusual (ie. not parallel to one of the
 181  * coordinate axes).
 182  */
 183 #define S_UNIT_X        0.50941539564955385     /* Pre-normalized */
 184 #define S_UNIT_Y        0.86052074622010633
 185 #else
 186 #define S_UNIT_X        1.0
 187 #define S_UNIT_Y        0.0
 188 #endif
 189 #endif
 190
 191 /* Determine the polygon normal and project vertices onto the plane
 192  * of the polygon.
 193  */
 194 void __gl_projectPolygon( GLUtesselator *tess )
 195 {
 196   GLUvertex *v, *vHead = &tess->mesh->vHead;
 197   GLdouble norm[3];
 198   GLdouble *sUnit, *tUnit;
 199   int i, computedNormal = FALSE;
 200
 201   norm[0] = tess->normal[0];
 202   norm[1] = tess->normal[1];
 203   norm[2] = tess->normal[2];
 204   if( norm[0] == 0 && norm[1] == 0 && norm[2] == 0 ) {
 205     ComputeNormal( tess, norm );
 206     computedNormal = TRUE;
 207   }
 208   sUnit = tess->sUnit;
 209   tUnit = tess->tUnit;
 210   i = LongAxis( norm );
 211
 212 #if defined(FOR_TRITE_TEST_PROGRAM) || defined(TRUE_PROJECT)
 213   /* Choose the initial sUnit vector to be approximately perpendicular
 214    * to the normal.
 215    */
 216   Normalize( norm );
 217
 218   sUnit[i] = 0;
 219   sUnit[(i+1)%3] = S_UNIT_X;
 220   sUnit[(i+2)%3] = S_UNIT_Y;
 221
 222   /* Now make it exactly perpendicular */
 223   w = Dot( sUnit, norm );
 224   sUnit[0] -= w * norm[0];
 225   sUnit[1] -= w * norm[1];
 226   sUnit[2] -= w * norm[2];
 227   Normalize( sUnit );
 228
 229   /* Choose tUnit so that (sUnit,tUnit,norm) form a right-handed frame */
 230   tUnit[0] = norm[1]*sUnit[2] - norm[2]*sUnit[1];
 231   tUnit[1] = norm[2]*sUnit[0] - norm[0]*sUnit[2];
 232   tUnit[2] = norm[0]*sUnit[1] - norm[1]*sUnit[0];
 233   Normalize( tUnit );
 234 #else
 235   /* Project perpendicular to a coordinate axis -- better numerically */
 236   sUnit[i] = 0;
 237   sUnit[(i+1)%3] = S_UNIT_X;
 238   sUnit[(i+2)%3] = S_UNIT_Y;
 239
 240   tUnit[i] = 0;
 241   tUnit[(i+1)%3] = (norm[i] > 0) ? -S_UNIT_Y : S_UNIT_Y;
 242   tUnit[(i+2)%3] = (norm[i] > 0) ? S_UNIT_X : -S_UNIT_X;
 243 #endif
 244
 245   /* Project the vertices onto the sweep plane */
 246   for( v = vHead->next; v != vHead; v = v->next ) {
 247     v->s = Dot( v->coords, sUnit );
 248     v->t = Dot( v->coords, tUnit );
 249   }
 250   if( computedNormal ) {
 251     CheckOrientation( tess );
 252   }
 253 }