src/common/tmesh.c

#ifndef lint
static const char RCSid[] = "$Id: tmesh.c,v 2.4 2003/11/14 17:22:06 schorsch Exp $";
#endif
/*
 * Compute and print barycentric coordinates for triangle meshes
 */

#include <stdio.h>

#include "fvect.h"

#include "tmesh.h"

#define ABS(x)          ((x) >= 0 ? (x) : -(x))


int
flat_tri(v1, v2, v3, n1, n2, n3)        /* determine if triangle is flat */
FVECT   v1, v2, v3, n1, n2, n3;
{
        double  d1, d2, d3;
        FVECT   vt1, vt2, vn;
                                        /* compute default normal */
        VSUB(vt1, v2, v1);
        VSUB(vt2, v3, v2);
        VCROSS(vn, vt1, vt2);
        if (normalize(vn) == 0.0)
                return(DEGEN);
                                        /* compare to supplied normals */
        d1 = DOT(vn, n1); d2 = DOT(vn, n2); d3 = DOT(vn, n3);
        if (d1 < 0 && d2 < 0 && d3 < 0) {
                if (d1 > -COSTOL || d2 > -COSTOL || d3 > -COSTOL)
                        return(RVBENT);
                return(RVFLAT);
        }
        if (d1 < COSTOL || d2 < COSTOL || d3 < COSTOL)
                return(ISBENT);
        return(ISFLAT);
}


int
comp_baryc(bcm, v1, v2, v3)             /* compute barycentric vectors */
register BARYCCM        *bcm;
RREAL                   *v1, *v2, *v3;
{
        RREAL   *vt;
        FVECT   va, vab, vcb;
        double  d;
        int     ax0, ax1;
        register int    i;
                                        /* compute major axis */
        VSUB(vab, v1, v2);
        VSUB(vcb, v3, v2);
        VCROSS(va, vab, vcb);
        bcm->ax = ABS(va[0]) > ABS(va[1]) ? 0 : 1;
        bcm->ax = ABS(va[bcm->ax]) > ABS(va[2]) ? bcm->ax : 2;
        if ((ax0 = bcm->ax + 1) >= 3) ax0 -= 3;
        if ((ax1 = ax0 + 1) >= 3) ax1 -= 3;
        for (i = 0; i < 2; i++) {
                vab[0] = v1[ax0] - v2[ax0];
                vcb[0] = v3[ax0] - v2[ax0];
                vab[1] = v1[ax1] - v2[ax1];
                vcb[1] = v3[ax1] - v2[ax1];
                d = vcb[0]*vcb[0] + vcb[1]*vcb[1];
                if (d <= FTINY*FTINY)
                        return(-1);
                d = (vcb[0]*vab[0]+vcb[1]*vab[1])/d;
                va[0] = vab[0] - vcb[0]*d;
                va[1] = vab[1] - vcb[1]*d;
                d = va[0]*va[0] + va[1]*va[1];
                if (d <= FTINY*FTINY)
                        return(-1);
                d = 1.0/d;
                bcm->tm[i][0] = va[0] *= d;
                bcm->tm[i][1] = va[1] *= d;
                bcm->tm[i][2] = -(v2[ax0]*va[0]+v2[ax1]*va[1]);
                                        /* rotate vertices */
                vt = v1;
                v1 = v2;
                v2 = v3;
                v3 = vt;
        }
        return(0);
}


void
eval_baryc(wt, p, bcm)          /* evaluate barycentric weights at p */
RREAL   wt[3];
FVECT   p;
register BARYCCM        *bcm;
{
        double  u, v;
        int     i;
        
        if ((i = bcm->ax + 1) >= 3) i -= 3;
        u = p[i];
        if (++i >= 3) i -= 3;
        v = p[i];
        wt[0] = u*bcm->tm[0][0] + v*bcm->tm[0][1] + bcm->tm[0][2];
        wt[1] = u*bcm->tm[1][0] + v*bcm->tm[1][1] + bcm->tm[1][2];
        wt[2] = 1. - wt[1] - wt[0];
}


int
get_baryc(wt, p, v1, v2, v3)    /* compute barycentric weights at p */
RREAL   wt[3];
FVECT   p;
FVECT   v1, v2, v3;
{
        BARYCCM bcm;
        
        if (comp_baryc(&bcm, v1, v2, v3) < 0)
                return(-1);
        eval_baryc(wt, p, &bcm);
        return(0);
}


#if 0
int
get_baryc(wt, p, v1, v2, v3)    /* compute barycentric weights at p */
RREAL   wt[3];
FVECT   p;
FVECT   v1, v2, v3;
{
        FVECT   ac, bc, pc, cros;
        double  normf;
                                /* area formula w/o 2-D optimization */
        VSUB(ac, v1, v3);
        VSUB(bc, v2, v3);
        VSUB(pc, p, v3);
        VCROSS(cros, ac, bc);
        normf = DOT(cros,cros)
        if (normf <= 0.0)
                return(-1);
        normf = 1./sqrt(normf);
        VCROSS(cros, bc, pc);
        wt[0] = VLEN(cros) * normf;
        VCROSS(cros, ac, pc);
        wt[1] = VLEN(cros) * normf;
        wt[2] = 1. - wt[1] - wt[0];
        return(0);
}
#endif


void
put_baryc(bcm, com, n)          /* put barycentric coord. vectors */
register BARYCCM        *bcm;
register RREAL          com[][3];
int                     n;
{
        double  a, b;
        register int    i;

        printf("%d\t%d\n", 1+3*n, bcm->ax);
        for (i = 0; i < n; i++) {
                a = com[i][0] - com[i][2];
                b = com[i][1] - com[i][2];
                printf("%14.8f %14.8f %14.8f\n",
                        bcm->tm[0][0]*a + bcm->tm[1][0]*b,
                        bcm->tm[0][1]*a + bcm->tm[1][1]*b,
                        bcm->tm[0][2]*a + bcm->tm[1][2]*b + com[i][2]);
        }
}
Revision:	2.5
Committed:	Thu Mar 2 17:16:56 2006 UTC (18 years, 2 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R2, rad4R2P2, rad5R0, rad5R1, rad4R2, rad4R1, rad4R0, rad3R8, rad3R9, rad4R2P1, rad5R3
Changes since 2.4:	+8 -5 lines
Log Message:	Replaced % operations with if statements (hopefully faster)
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.5	static const char RCSid[] = "$Id: tmesh.c,v 2.4 2003/11/14 17:22:06 schorsch Exp $";
3	greg	2.1	#endif
4			/*
5			* Compute and print barycentric coordinates for triangle meshes
6			*/
7
8			#include <stdio.h>
9
10			#include "fvect.h"
11
12			#include "tmesh.h"
13
14			#define ABS(x) ((x) >= 0 ? (x) : -(x))
15
16
17			int
18			flat_tri(v1, v2, v3, n1, n2, n3) /* determine if triangle is flat */
19			FVECT v1, v2, v3, n1, n2, n3;
20			{
21			double d1, d2, d3;
22			FVECT vt1, vt2, vn;
23			/* compute default normal */
24	greg	2.2	VSUB(vt1, v2, v1);
25			VSUB(vt2, v3, v2);
26			VCROSS(vn, vt1, vt2);
27	greg	2.1	if (normalize(vn) == 0.0)
28			return(DEGEN);
29			/* compare to supplied normals */
30			d1 = DOT(vn, n1); d2 = DOT(vn, n2); d3 = DOT(vn, n3);
31			if (d1 < 0 && d2 < 0 && d3 < 0) {
32			if (d1 > -COSTOL \|\| d2 > -COSTOL \|\| d3 > -COSTOL)
33			return(RVBENT);
34			return(RVFLAT);
35			}
36			if (d1 < COSTOL \|\| d2 < COSTOL \|\| d3 < COSTOL)
37			return(ISBENT);
38			return(ISFLAT);
39			}
40
41
42			int
43			comp_baryc(bcm, v1, v2, v3) /* compute barycentric vectors */
44			register BARYCCM *bcm;
45	schorsch	2.3	RREAL v1, v2, *v3;
46	greg	2.1	{
47	schorsch	2.3	RREAL *vt;
48	greg	2.1	FVECT va, vab, vcb;
49			double d;
50			int ax0, ax1;
51			register int i;
52			/* compute major axis */
53	greg	2.2	VSUB(vab, v1, v2);
54			VSUB(vcb, v3, v2);
55			VCROSS(va, vab, vcb);
56	greg	2.1	bcm->ax = ABS(va[0]) > ABS(va[1]) ? 0 : 1;
57			bcm->ax = ABS(va[bcm->ax]) > ABS(va[2]) ? bcm->ax : 2;
58	greg	2.5	if ((ax0 = bcm->ax + 1) >= 3) ax0 -= 3;
59			if ((ax1 = ax0 + 1) >= 3) ax1 -= 3;
60	greg	2.1	for (i = 0; i < 2; i++) {
61			vab[0] = v1[ax0] - v2[ax0];
62			vcb[0] = v3[ax0] - v2[ax0];
63			vab[1] = v1[ax1] - v2[ax1];
64			vcb[1] = v3[ax1] - v2[ax1];
65			d = vcb[0]vcb[0] + vcb[1]vcb[1];
66			if (d <= FTINY*FTINY)
67			return(-1);
68			d = (vcb[0]vab[0]+vcb[1]vab[1])/d;
69			va[0] = vab[0] - vcb[0]*d;
70			va[1] = vab[1] - vcb[1]*d;
71			d = va[0]va[0] + va[1]va[1];
72			if (d <= FTINY*FTINY)
73			return(-1);
74			d = 1.0/d;
75			bcm->tm[i][0] = va[0] *= d;
76			bcm->tm[i][1] = va[1] *= d;
77			bcm->tm[i][2] = -(v2[ax0]va[0]+v2[ax1]va[1]);
78			/* rotate vertices */
79			vt = v1;
80			v1 = v2;
81			v2 = v3;
82			v3 = vt;
83			}
84			return(0);
85			}
86
87
88	greg	2.2	void
89			eval_baryc(wt, p, bcm) /* evaluate barycentric weights at p */
90	schorsch	2.3	RREAL wt[3];
91	greg	2.2	FVECT p;
92			register BARYCCM *bcm;
93			{
94			double u, v;
95	greg	2.5	int i;
96	greg	2.2
97	greg	2.5	if ((i = bcm->ax + 1) >= 3) i -= 3;
98			u = p[i];
99			if (++i >= 3) i -= 3;
100			v = p[i];
101	greg	2.2	wt[0] = ubcm->tm[0][0] + vbcm->tm[0][1] + bcm->tm[0][2];
102			wt[1] = ubcm->tm[1][0] + vbcm->tm[1][1] + bcm->tm[1][2];
103			wt[2] = 1. - wt[1] - wt[0];
104			}
105
106
107			int
108			get_baryc(wt, p, v1, v2, v3) /* compute barycentric weights at p */
109	schorsch	2.3	RREAL wt[3];
110	greg	2.2	FVECT p;
111			FVECT v1, v2, v3;
112			{
113			BARYCCM bcm;
114
115			if (comp_baryc(&bcm, v1, v2, v3) < 0)
116			return(-1);
117			eval_baryc(wt, p, &bcm);
118			return(0);
119			}
120
121
122			#if 0
123			int
124			get_baryc(wt, p, v1, v2, v3) /* compute barycentric weights at p */
125	schorsch	2.3	RREAL wt[3];
126	greg	2.2	FVECT p;
127			FVECT v1, v2, v3;
128			{
129			FVECT ac, bc, pc, cros;
130			double normf;
131			/* area formula w/o 2-D optimization */
132			VSUB(ac, v1, v3);
133			VSUB(bc, v2, v3);
134			VSUB(pc, p, v3);
135			VCROSS(cros, ac, bc);
136			normf = DOT(cros,cros)
137			if (normf <= 0.0)
138			return(-1);
139			normf = 1./sqrt(normf);
140			VCROSS(cros, bc, pc);
141			wt[0] = VLEN(cros) * normf;
142			VCROSS(cros, ac, pc);
143			wt[1] = VLEN(cros) * normf;
144			wt[2] = 1. - wt[1] - wt[0];
145			return(0);
146			}
147			#endif
148
149
150			void
151			put_baryc(bcm, com, n) /* put barycentric coord. vectors */
152	greg	2.1	register BARYCCM *bcm;
153	schorsch	2.3	register RREAL com[][3];
154	greg	2.2	int n;
155	greg	2.1	{
156			double a, b;
157	schorsch	2.4	register int i;
158	greg	2.1
159			printf("%d\t%d\n", 1+3*n, bcm->ax);
160			for (i = 0; i < n; i++) {
161			a = com[i][0] - com[i][2];
162			b = com[i][1] - com[i][2];
163			printf("%14.8f %14.8f %14.8f\n",
164			bcm->tm[0][0]a + bcm->tm[1][0]b,
165			bcm->tm[0][1]a + bcm->tm[1][1]b,
166			bcm->tm[0][2]a + bcm->tm[1][2]b + com[i][2]);
167			}
168			}