File indexing completed on 2024-03-24 15:17:33

0001 /*
0002 ** Author: Eric Veach, July 1994.
0003 **
0004 */
0005 
0006 #include "gluos.h"
0007 #include "mesh.h"
0008 #include "tess.h"
0009 #include "normal.h"
0010 #include <math.h>
0011 #include <assert.h>
0012 
0013 #ifndef TRUE
0014 #define TRUE 1
0015 #endif
0016 #ifndef FALSE
0017 #define FALSE 0
0018 #endif
0019 
0020 #define Dot(u, v) (u[0] * v[0] + u[1] * v[1] + u[2] * v[2])
0021 
0022 #if 0
0023 static void Normalize( GLdouble v[3] )
0024 {
0025   GLdouble len = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
0026 
0027   assert( len > 0 );
0028   len = sqrt( len );
0029   v[0] /= len;
0030   v[1] /= len;
0031   v[2] /= len;
0032 }
0033 #endif
0034 
0035 #undef ABS
0036 #define ABS(x) ((x) < 0 ? -(x) : (x))
0037 
0038 static int LongAxis(GLdouble v[3])
0039 {
0040     int i = 0;
0041 
0042     if (ABS(v[1]) > ABS(v[0]))
0043     {
0044         i = 1;
0045     }
0046     if (ABS(v[2]) > ABS(v[i]))
0047     {
0048         i = 2;
0049     }
0050     return i;
0051 }
0052 
0053 static void ComputeNormal(GLUtesselator *tess, GLdouble norm[3])
0054 {
0055     GLUvertex *v, *v1, *v2;
0056     GLdouble c, tLen2, maxLen2;
0057     GLdouble maxVal[3], minVal[3], d1[3], d2[3], tNorm[3];
0058     GLUvertex *maxVert[3], *minVert[3];
0059     GLUvertex *vHead = &tess->mesh->vHead;
0060     int i;
0061 
0062     maxVal[0] = maxVal[1] = maxVal[2] = -2 * GLU_TESS_MAX_COORD;
0063     minVal[0] = minVal[1] = minVal[2] = 2 * GLU_TESS_MAX_COORD;
0064 
0065     for (v = vHead->next; v != vHead; v = v->next)
0066     {
0067         for (i = 0; i < 3; ++i)
0068         {
0069             c = v->coords[i];
0070             if (c < minVal[i])
0071             {
0072                 minVal[i]  = c;
0073                 minVert[i] = v;
0074             }
0075             if (c > maxVal[i])
0076             {
0077                 maxVal[i]  = c;
0078                 maxVert[i] = v;
0079             }
0080         }
0081     }
0082 
0083     /* Find two vertices separated by at least 1/sqrt(3) of the maximum
0084    * distance between any two vertices
0085    */
0086     i = 0;
0087     if (maxVal[1] - minVal[1] > maxVal[0] - minVal[0])
0088     {
0089         i = 1;
0090     }
0091     if (maxVal[2] - minVal[2] > maxVal[i] - minVal[i])
0092     {
0093         i = 2;
0094     }
0095     if (minVal[i] >= maxVal[i])
0096     {
0097         /* All vertices are the same -- normal doesn't matter */
0098         norm[0] = 0;
0099         norm[1] = 0;
0100         norm[2] = 1;
0101         return;
0102     }
0103 
0104     /* Look for a third vertex which forms the triangle with maximum area
0105    * (Length of normal == twice the triangle area)
0106    */
0107     maxLen2 = 0;
0108     v1      = minVert[i];
0109     v2      = maxVert[i];
0110     d1[0]   = v1->coords[0] - v2->coords[0];
0111     d1[1]   = v1->coords[1] - v2->coords[1];
0112     d1[2]   = v1->coords[2] - v2->coords[2];
0113     for (v = vHead->next; v != vHead; v = v->next)
0114     {
0115         d2[0]    = v->coords[0] - v2->coords[0];
0116         d2[1]    = v->coords[1] - v2->coords[1];
0117         d2[2]    = v->coords[2] - v2->coords[2];
0118         tNorm[0] = d1[1] * d2[2] - d1[2] * d2[1];
0119         tNorm[1] = d1[2] * d2[0] - d1[0] * d2[2];
0120         tNorm[2] = d1[0] * d2[1] - d1[1] * d2[0];
0121         tLen2    = tNorm[0] * tNorm[0] + tNorm[1] * tNorm[1] + tNorm[2] * tNorm[2];
0122         if (tLen2 > maxLen2)
0123         {
0124             maxLen2 = tLen2;
0125             norm[0] = tNorm[0];
0126             norm[1] = tNorm[1];
0127             norm[2] = tNorm[2];
0128         }
0129     }
0130 
0131     if (maxLen2 <= 0)
0132     {
0133         /* All points lie on a single line -- any decent normal will do */
0134         norm[0] = norm[1] = norm[2] = 0;
0135         norm[LongAxis(d1)]          = 1;
0136     }
0137 }
0138 
0139 static void CheckOrientation(GLUtesselator *tess)
0140 {
0141     GLdouble area;
0142     GLUface *f, *fHead   = &tess->mesh->fHead;
0143     GLUvertex *v, *vHead = &tess->mesh->vHead;
0144     GLUhalfEdge *e;
0145 
0146     /* When we compute the normal automatically, we choose the orientation
0147    * so that the sum of the signed areas of all contours is non-negative.
0148    */
0149     area = 0;
0150     for (f = fHead->next; f != fHead; f = f->next)
0151     {
0152         e = f->anEdge;
0153         if (e->winding <= 0)
0154             continue;
0155         do
0156         {
0157             area += (e->Org->s - e->Dst->s) * (e->Org->t + e->Dst->t);
0158             e = e->Lnext;
0159         } while (e != f->anEdge);
0160     }
0161     if (area < 0)
0162     {
0163         /* Reverse the orientation by flipping all the t-coordinates */
0164         for (v = vHead->next; v != vHead; v = v->next)
0165         {
0166             v->t = -v->t;
0167         }
0168         tess->tUnit[0] = -tess->tUnit[0];
0169         tess->tUnit[1] = -tess->tUnit[1];
0170         tess->tUnit[2] = -tess->tUnit[2];
0171     }
0172 }
0173 
0174 #ifdef FOR_TRITE_TEST_PROGRAM
0175 #include <stdlib.h>
0176 extern int RandomSweep;
0177 #define S_UNIT_X (RandomSweep ? (2 * drand48() - 1) : 1.0)
0178 #define S_UNIT_Y (RandomSweep ? (2 * drand48() - 1) : 0.0)
0179 #else
0180 #if defined(SLANTED_SWEEP)
0181 /* The "feature merging" is not intended to be complete.  There are
0182  * special cases where edges are nearly parallel to the sweep line
0183  * which are not implemented.  The algorithm should still behave
0184  * robustly (ie. produce a reasonable tesselation) in the presence
0185  * of such edges, however it may miss features which could have been
0186  * merged.  We could minimize this effect by choosing the sweep line
0187  * direction to be something unusual (ie. not parallel to one of the
0188  * coordinate axes).
0189  */
0190 #define S_UNIT_X 0.50941539564955385 /* Pre-normalized */
0191 #define S_UNIT_Y 0.86052074622010633
0192 #else
0193 #define S_UNIT_X 1.0
0194 #define S_UNIT_Y 0.0
0195 #endif
0196 #endif
0197 
0198 /* Determine the polygon normal and project vertices onto the plane
0199  * of the polygon.
0200  */
0201 void __gl_projectPolygon(GLUtesselator *tess)
0202 {
0203     GLUvertex *v, *vHead = &tess->mesh->vHead;
0204     GLdouble norm[3];
0205     GLdouble *sUnit, *tUnit;
0206     int i, computedNormal = FALSE;
0207 
0208     norm[0] = tess->normal[0];
0209     norm[1] = tess->normal[1];
0210     norm[2] = tess->normal[2];
0211     if (norm[0] == 0 && norm[1] == 0 && norm[2] == 0)
0212     {
0213         ComputeNormal(tess, norm);
0214         computedNormal = TRUE;
0215     }
0216     sUnit = tess->sUnit;
0217     tUnit = tess->tUnit;
0218     i     = LongAxis(norm);
0219 
0220 #if defined(FOR_TRITE_TEST_PROGRAM) || defined(TRUE_PROJECT)
0221     /* Choose the initial sUnit vector to be approximately perpendicular
0222    * to the normal.
0223    */
0224     Normalize(norm);
0225 
0226     sUnit[i]           = 0;
0227     sUnit[(i + 1) % 3] = S_UNIT_X;
0228     sUnit[(i + 2) % 3] = S_UNIT_Y;
0229 
0230     /* Now make it exactly perpendicular */
0231     w = Dot(sUnit, norm);
0232     sUnit[0] -= w * norm[0];
0233     sUnit[1] -= w * norm[1];
0234     sUnit[2] -= w * norm[2];
0235     Normalize(sUnit);
0236 
0237     /* Choose tUnit so that (sUnit,tUnit,norm) form a right-handed frame */
0238     tUnit[0] = norm[1] * sUnit[2] - norm[2] * sUnit[1];
0239     tUnit[1] = norm[2] * sUnit[0] - norm[0] * sUnit[2];
0240     tUnit[2] = norm[0] * sUnit[1] - norm[1] * sUnit[0];
0241     Normalize(tUnit);
0242 #else
0243     /* Project perpendicular to a coordinate axis -- better numerically */
0244     sUnit[i]           = 0;
0245     sUnit[(i + 1) % 3] = S_UNIT_X;
0246     sUnit[(i + 2) % 3] = S_UNIT_Y;
0247 
0248     tUnit[i]           = 0;
0249     tUnit[(i + 1) % 3] = (norm[i] > 0) ? -S_UNIT_Y : S_UNIT_Y;
0250     tUnit[(i + 2) % 3] = (norm[i] > 0) ? S_UNIT_X : -S_UNIT_X;
0251 #endif
0252 
0253     /* Project the vertices onto the sweep plane */
0254     for (v = vHead->next; v != vHead; v = v->next)
0255     {
0256         v->s = Dot(v->coords, sUnit);
0257         v->t = Dot(v->coords, tUnit);
0258     }
0259     if (computedNormal)
0260     {
0261         CheckOrientation(tess);
0262     }
0263 }