--- ray/src/common/face.c	1992/10/02 16:02:20	2.2
+++ ray/src/common/face.c	2021/02/12 00:47:33	2.15
@@ -1,15 +1,12 @@
-/* Copyright (c) 1991 Regents of the University of California */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
+static const char RCSid[] = "$Id: face.c,v 2.15 2021/02/12 00:47:33 greg Exp $";
 #endif
-
 /*
  *  face.c - routines dealing with polygonal faces.
- *
- *     8/30/85
  */
 
+#include "copyright.h"
+
 #include  "standard.h"
 
 #include  "object.h"
@@ -27,18 +24,23 @@ static char SCCSid[] = "$SunId$ LBL";
  *  The last vertex is automatically connected to the first.
  */
 
-#define  VERTEPS	1e-4		/* allowed vertex error */
+#ifdef  SMLFLT
+#define  VERTEPS	1e-3		/* allowed vertex error */
+#else
+#define  VERTEPS	1e-5		/* allowed vertex error */
+#endif
 
 
 FACE *
-getface(o)			/* get arguments for a face */
-OBJREC  *o;
+getface(				/* get arguments for a face */
+	OBJREC  *o
+)
 {
 	double  d1;
-	int  badvert;
+	int  smalloff, badvert;
 	FVECT  v1, v2, v3;
-	register FACE  *f;
-	register int  i;
+	FACE  *f;
+	int  i;
 
 	if ((f = (FACE *)o->os) != NULL)
 		return(f);			/* already done */
@@ -54,10 +56,15 @@ OBJREC  *o;
 
 	f->va = o->oargs.farg;
 	f->nv = o->oargs.nfargs / 3;
+						/* check for last==first */
+	if (f->nv > 3 && dist2(VERTEX(f,0),VERTEX(f,f->nv-1)) <= FTINY*FTINY)
+		f->nv--;
 						/* compute area and normal */
 	f->norm[0] = f->norm[1] = f->norm[2] = 0.0;
-	v1[0] = v1[1] = v1[2] = 0.0;
-	for (i = 1; i < f->nv; i++) {
+	v1[0] = VERTEX(f,1)[0] - VERTEX(f,0)[0];
+	v1[1] = VERTEX(f,1)[1] - VERTEX(f,0)[1];
+	v1[2] = VERTEX(f,1)[2] - VERTEX(f,0)[2];
+	for (i = 2; i < f->nv; i++) {
 		v2[0] = VERTEX(f,i)[0] - VERTEX(f,0)[0];
 		v2[1] = VERTEX(f,i)[1] - VERTEX(f,0)[1];
 		v2[2] = VERTEX(f,i)[2] - VERTEX(f,0)[2];
@@ -78,16 +85,20 @@ OBJREC  *o;
 						/* compute offset */
 	badvert = 0;
 	f->offset = DOT(f->norm, VERTEX(f,0));
+	smalloff = fabs(f->offset) <= VERTEPS;
 	for (i = 1; i < f->nv; i++) {
 		d1 = DOT(f->norm, VERTEX(f,i));
-		badvert += fabs(d1 - f->offset/i) > VERTEPS;
+		if (smalloff)
+			badvert += fabs(d1 - f->offset/i) > VERTEPS;
+		else
+			badvert += fabs(1.0 - d1*i/f->offset) > VERTEPS;
 		f->offset += d1;
 	}
 	f->offset /= (double)f->nv;
-	if (badvert)
+	if (f->nv > 3 && badvert)
 		objerror(o, WARNING, "non-planar vertex");
 						/* find axis */
-	f->ax = fabs(f->norm[0]) > fabs(f->norm[1]) ? 0 : 1;
+	f->ax = (fabs(f->norm[0]) > fabs(f->norm[1]));
 	if (fabs(f->norm[2]) > fabs(f->norm[f->ax]))
 		f->ax = 2;
 
@@ -95,8 +106,10 @@ OBJREC  *o;
 }
 
 
-freeface(o)			/* free memory associated with face */
-OBJREC  *o;
+void
+freeface(			/* free memory associated with face */
+	OBJREC  *o
+)
 {
 	if (o->os == NULL)
 		return;
@@ -105,17 +118,20 @@ OBJREC  *o;
 }
 
 
-inface(p, f)			/* determine if point is in face */
-FVECT  p;
-FACE  *f;
+int
+inface(				/* determine if point is in face */
+	FVECT  p,
+	FACE  *f
+)
 {
 	int  ncross, n;
 	double  x, y;
-	register int  xi, yi;
-	register FLOAT  *p0, *p1;
+	int  tst;
+	int  xi, yi;
+	RREAL  *p0, *p1;
 
-	xi = (f->ax+1)%3;
-	yi = (f->ax+2)%3;
+	if ((xi = f->ax + 1) >= 3) xi -= 3;
+	if ((yi = xi + 1) >= 3) yi -= 3;
 	x = p[xi];
 	y = p[yi];
 	n = f->nv;
@@ -124,13 +140,21 @@ FACE  *f;
 	ncross = 0;
 					/* positive x axis cross test */
 	while (n--) {
-		if ((p0[yi] > y) ^ (p1[yi] > y))
-			if (p0[xi] > x && p1[xi] > x)
+		if ((p0[yi] > y) ^ (p1[yi] > y)) {
+			tst = (p0[xi] > x) + (p1[xi] > x);
+			if (tst == 2)
 				ncross++;
-			else if (p0[xi] > x || p1[xi] > x)
-				ncross += (p1[yi] > p0[yi]) ^
-						((p0[yi]-y)*(p1[xi]-x) >
-						(p0[xi]-x)*(p1[yi]-y));
+			else if (tst) {
+				double	prodA = (p0[yi]-y)*(p1[xi]-x);
+				double	prodB = (p0[xi]-x)*(p1[yi]-y);
+				if (FRELEQ(prodA, prodB))
+					return(1);	/* edge case #1 */
+				ncross += (p1[yi] > p0[yi]) ^ (prodA > prodB);
+			} else if (FRELEQ(p0[xi], x) && FRELEQ(p1[xi], x))
+				return(1);		/* edge case #2 */
+		} else if (FRELEQ(p0[yi], y) && FRELEQ(p1[yi], y) &&
+				(p0[xi] > x) ^ (p1[xi] > x))
+			return(1);			/* edge case #3 */
 		p0 = p1;
 		p1 += 3;
 	}