src/common/face.c

#ifndef lint
static const char RCSid[] = "$Id: face.c,v 2.17 2024/08/18 00:06:10 greg Exp $";
#endif
/*
 *  face.c - routines dealing with polygonal faces.
 */

#include "copyright.h"

#include  "standard.h"

#include  "object.h"

#include  "face.h"

/*
 *      A face is given as a list of 3D vertices.  The normal
 *  direction and therefore the surface orientation is determined
 *  by the ordering of the vertices.  Looking in the direction opposite
 *  the normal (at the front of the face), the vertices will be
 *  listed in counter-clockwise order.
 *      There is no checking done to insure that the edges do not cross
 *  one another.  This was considered too expensive and should be unnecessary.
 *  The last vertex is automatically connected to the first.
 */

#ifdef  SMLFLT
#define  VERTEPS        1e-3            /* allowed vertex error */
#else
#define  VERTEPS        1e-5            /* allowed vertex error */
#endif


FACE *
getface(                                /* get arguments for a face */
        OBJREC  *o
)
{
        double  d1;
        int  smalloff, badvert;
        FVECT  v1, v2, v3;
        FACE  *f;
        int  i;

        if ((f = (FACE *)o->os) != NULL)
                return(f);                      /* already done */

        f = (FACE *)malloc(sizeof(FACE));
        if (f == NULL)
                error(SYSTEM, "out of memory in makeface");

        if (o->oargs.nfargs < 9 || o->oargs.nfargs % 3)
                objerror(o, USER, "bad # arguments");

        o->os = (char *)f;                      /* save face */

        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] = 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];
                fcross(v3, v1, v2);
                f->norm[0] += v3[0];
                f->norm[1] += v3[1];
                f->norm[2] += v3[2];
                VCOPY(v1, v2);
        }
        f->area = normalize(f->norm);
        if (f->area == 0.0) {
                objerror(o, WARNING, "zero area");      /* used to be fatal */
                f->offset = 0.0;
                f->ax = 0;
                return(f);
        }
        f->area *= 0.5;
                                                /* 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));
                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 (f->nv > 3 && badvert)
                objerror(o, WARNING, "non-planar vertex");
                                                /* find axis */
        f->ax = (fabs(f->norm[1]) > fabs(f->norm[0]));
        if (fabs(f->norm[2]) > fabs(f->norm[f->ax]))
                f->ax = 2;

        return(f);
}


void
freeface(                       /* free memory associated with face */
        OBJREC  *o
)
{
        if (o->os == NULL)
                return;
        free(o->os);
        o->os = NULL;
}

int
inface(                         /* determine if point is in face */
        FVECT  p,
        FACE  *f
)
{
        int  ncross, n;
        double  x, y;
        int  tst;
        int  xi, yi;
        RREAL  *p0, *p1;

        if ((xi = f->ax + 1) >= 3) xi -= 3;
        if ((yi = xi + 1) >= 3) yi -= 3;
        x = p[xi];
        y = p[yi];
        n = f->nv;
        p0 = f->va + 3*(n-1);           /* connect last to first */
        p1 = f->va;
        ncross = 0;
                                        /* positive x axis cross test */
        while (n--) {
                if (FABSEQ(p0[yi], y) && FABSEQ(p1[yi], y) &&
                                (p0[xi] > x) ^ (p1[xi] > x))
                        return(1);                      /* edge case #1 */
                if ((p0[yi] > y) ^ (p1[yi] > y)) {
                        tst = (p0[xi] > x) + (p1[xi] > x);
                        if (tst == 2)
                                ncross++;
                        else if (tst) {
                                double  prodA = (p0[yi]-y)*(p1[xi]-x);
                                double  prodB = (p0[xi]-x)*(p1[yi]-y);
                                if (FABSEQ(prodA, prodB))
                                        return(1);      /* edge case #2 */
                                ncross += (p1[yi] > p0[yi]) ^ (prodA > prodB);
                        } else if (FABSEQ(p0[xi], x) && FABSEQ(p1[xi], x))
                                return(1);              /* edge case #3 */
                }
                p0 = p1;
                p1 += 3;
        }
        return(ncross & 01);
}
Revision:	2.18
Committed:	Mon Sep 9 02:48:03 2024 UTC (7 months, 3 weeks ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 2.17:	+4 -5 lines
Log Message:	fix: Thanks to Peter A-B for spotting errors caused by FRELEQ()
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	2.18	static const char RCSid[] = "$Id: face.c,v 2.17 2024/08/18 00:06:10 greg Exp $";
3	greg	1.1	#endif
4			/*
5			* face.c - routines dealing with polygonal faces.
6	greg	2.6	*/
7
8	greg	2.7	#include "copyright.h"
9	greg	1.1
10			#include "standard.h"
11
12			#include "object.h"
13
14			#include "face.h"
15
16			/*
17			* A face is given as a list of 3D vertices. The normal
18			* direction and therefore the surface orientation is determined
19			* by the ordering of the vertices. Looking in the direction opposite
20			* the normal (at the front of the face), the vertices will be
21			* listed in counter-clockwise order.
22			* There is no checking done to insure that the edges do not cross
23			* one another. This was considered too expensive and should be unnecessary.
24			* The last vertex is automatically connected to the first.
25			*/
26
27	greg	2.3	#ifdef SMLFLT
28	greg	2.9	#define VERTEPS 1e-3 /* allowed vertex error */
29	greg	2.3	#else
30	greg	2.9	#define VERTEPS 1e-5 /* allowed vertex error */
31	greg	2.3	#endif
32	greg	1.1
33
34			FACE *
35	greg	2.13	getface( /* get arguments for a face */
36			OBJREC *o
37			)
38	greg	1.1	{
39			double d1;
40	greg	2.10	int smalloff, badvert;
41	greg	1.1	FVECT v1, v2, v3;
42	greg	2.13	FACE *f;
43			int i;
44	greg	1.1
45			if ((f = (FACE *)o->os) != NULL)
46			return(f); /* already done */
47
48			f = (FACE *)malloc(sizeof(FACE));
49			if (f == NULL)
50			error(SYSTEM, "out of memory in makeface");
51
52			if (o->oargs.nfargs < 9 \|\| o->oargs.nfargs % 3)
53			objerror(o, USER, "bad # arguments");
54
55	greg	1.3	o->os = (char )f; / save face */
56
57	greg	1.1	f->va = o->oargs.farg;
58			f->nv = o->oargs.nfargs / 3;
59	greg	2.4	/* check for last==first */
60	greg	2.14	if (f->nv > 3 && dist2(VERTEX(f,0),VERTEX(f,f->nv-1)) <= FTINY*FTINY)
61	greg	2.4	f->nv--;
62	greg	1.1	/* compute area and normal */
63			f->norm[0] = f->norm[1] = f->norm[2] = 0.0;
64	greg	2.5	v1[0] = VERTEX(f,1)[0] - VERTEX(f,0)[0];
65			v1[1] = VERTEX(f,1)[1] - VERTEX(f,0)[1];
66			v1[2] = VERTEX(f,1)[2] - VERTEX(f,0)[2];
67			for (i = 2; i < f->nv; i++) {
68	greg	1.1	v2[0] = VERTEX(f,i)[0] - VERTEX(f,0)[0];
69			v2[1] = VERTEX(f,i)[1] - VERTEX(f,0)[1];
70			v2[2] = VERTEX(f,i)[2] - VERTEX(f,0)[2];
71			fcross(v3, v1, v2);
72			f->norm[0] += v3[0];
73			f->norm[1] += v3[1];
74			f->norm[2] += v3[2];
75			VCOPY(v1, v2);
76			}
77			f->area = normalize(f->norm);
78			if (f->area == 0.0) {
79			objerror(o, WARNING, "zero area"); /* used to be fatal */
80	greg	1.2	f->offset = 0.0;
81	greg	1.1	f->ax = 0;
82			return(f);
83			}
84			f->area *= 0.5;
85	greg	1.2	/* compute offset */
86	greg	1.1	badvert = 0;
87	greg	1.2	f->offset = DOT(f->norm, VERTEX(f,0));
88	greg	2.10	smalloff = fabs(f->offset) <= VERTEPS;
89	greg	1.1	for (i = 1; i < f->nv; i++) {
90			d1 = DOT(f->norm, VERTEX(f,i));
91	greg	2.10	if (smalloff)
92			badvert += fabs(d1 - f->offset/i) > VERTEPS;
93			else
94			badvert += fabs(1.0 - d1*i/f->offset) > VERTEPS;
95	greg	1.2	f->offset += d1;
96	greg	1.1	}
97	greg	1.5	f->offset /= (double)f->nv;
98	greg	2.9	if (f->nv > 3 && badvert)
99	greg	1.1	objerror(o, WARNING, "non-planar vertex");
100			/* find axis */
101	greg	2.16	f->ax = (fabs(f->norm[1]) > fabs(f->norm[0]));
102	greg	1.1	if (fabs(f->norm[2]) > fabs(f->norm[f->ax]))
103			f->ax = 2;
104
105			return(f);
106			}
107
108
109	greg	2.6	void
110	greg	2.13	freeface( /* free memory associated with face */
111			OBJREC *o
112			)
113	greg	1.1	{
114	greg	1.4	if (o->os == NULL)
115			return;
116	greg	1.1	free(o->os);
117			o->os = NULL;
118			}
119
120	greg	2.6	int
121	greg	2.13	inface( /* determine if point is in face */
122			FVECT p,
123			FACE *f
124			)
125	greg	1.1	{
126			int ncross, n;
127			double x, y;
128	greg	2.8	int tst;
129	greg	2.13	int xi, yi;
130			RREAL p0, p1;
131	greg	1.1
132	greg	2.12	if ((xi = f->ax + 1) >= 3) xi -= 3;
133			if ((yi = xi + 1) >= 3) yi -= 3;
134	greg	1.1	x = p[xi];
135			y = p[yi];
136			n = f->nv;
137			p0 = f->va + 3(n-1); / connect last to first */
138			p1 = f->va;
139			ncross = 0;
140			/* positive x axis cross test */
141			while (n--) {
142	greg	2.18	if (FABSEQ(p0[yi], y) && FABSEQ(p1[yi], y) &&
143	greg	2.17	(p0[xi] > x) ^ (p1[xi] > x))
144			return(1); /* edge case #1 */
145	greg	2.8	if ((p0[yi] > y) ^ (p1[yi] > y)) {
146			tst = (p0[xi] > x) + (p1[xi] > x);
147			if (tst == 2)
148	greg	1.1	ncross++;
149	greg	2.15	else if (tst) {
150			double prodA = (p0[yi]-y)*(p1[xi]-x);
151			double prodB = (p0[xi]-x)*(p1[yi]-y);
152	greg	2.18	if (FABSEQ(prodA, prodB))
153	greg	2.17	return(1); /* edge case #2 */
154	greg	2.15	ncross += (p1[yi] > p0[yi]) ^ (prodA > prodB);
155	greg	2.18	} else if (FABSEQ(p0[xi], x) && FABSEQ(p1[xi], x))
156	greg	2.17	return(1); /* edge case #3 */
157			}
158	greg	1.1	p0 = p1;
159			p1 += 3;
160			}
161			return(ncross & 01);
162			}