src/rt/o_mesh.c

#ifndef lint
static const char RCSid[] = "$Id: o_mesh.c,v 2.9 2003/07/21 22:30:19 schorsch Exp $";
#endif
/*
 *  Routines for computing ray intersections with meshes.
 *
 *  Intersection with a triangle mesh is based on Segura and Feito's
 *  WSCG 2001 paper, "Algorithms to Test Ray-Triangle Intersection,
 *  Comparative Study."  This method avoids additional storage
 *  requirements, floating divides, and allows some savings by
 *  caching ray-edge comparisons that are otherwise repeated locally
 *  in typical mesh geometries.  (This is our own optimization.)
 *
 *  The code herein is quite similar to that in o_instance.c, the
 *  chief differences being the custom triangle intersection routines
 *  and the fact that an "OBJECT" in the mesh octree is not an index
 *  into the Radiance OBJREC list, but a mesh triangle index.  We still
 *  utilize the standard octree traversal code by setting the hitf
 *  function pointer in the RAY struct to our custom mesh_hit() call.
 */

#include  "copyright.h"

#include <string.h>

#include  "ray.h"
#include  "mesh.h"
#include  "tmesh.h"
#include  "rtotypes.h"


#define  EDGE_CACHE_SIZ         251     /* length of mesh edge cache */

#define  curmi                  (edge_cache.mi)
#define  curmsh                 (curmi->msh)


/* Cache of signed volumes for this ray and this mesh */
struct EdgeCache {
        OBJREC          *o;     /* mesh object */
        MESHINST        *mi;    /* current mesh instance */
        struct EdgeSide {
                int32   v1i, v2i;       /* vertex indices (lowest first) */
                short   signum;         /* signed volume */
        }               cache[EDGE_CACHE_SIZ];
}       edge_cache;


static void
prep_edge_cache(o)              /* get instance and clear edge cache */
OBJREC  *o;
{
                                        /* get mesh instance */
        edge_cache.mi = getmeshinst(edge_cache.o = o, IO_ALL);
                                        /* clear edge cache */
        memset((void *)edge_cache.cache, '\0', sizeof(edge_cache.cache));
}


static int
signed_volume(r, v1, v2)        /* get signed volume for ray and edge */
register RAY    *r;
int32           v1, v2;
{
        int                             reversed = 0;
        register struct EdgeSide        *ecp;
        
        if (v1 > v2) {
                int32   t = v2; v2 = v1; v1 = t;
                reversed = 1;
        }
        ecp = &edge_cache.cache[((v2<<11 ^ v1) & 0x7fffffff) % EDGE_CACHE_SIZ];
        if (ecp->v1i != v1 || ecp->v2i != v2) {
                MESHVERT        tv1, tv2;       /* compute signed volume */
                FVECT           v2d;
                double          vol;
                if (!getmeshvert(&tv1, edge_cache.mi->msh, v1, MT_V) ||
                            !getmeshvert(&tv2, edge_cache.mi->msh, v2, MT_V))
                        objerror(edge_cache.o, INTERNAL,
                                "missing mesh vertex in signed_volume");
                VSUB(v2d, tv2.v, r->rorg);
                vol = (tv1.v[0] - r->rorg[0]) *
                                (v2d[1]*r->rdir[2] - v2d[2]*r->rdir[1]);
                vol += (tv1.v[1] - r->rorg[1]) *
                                (v2d[2]*r->rdir[0] - v2d[0]*r->rdir[2]);
                vol += (tv1.v[2] - r->rorg[2]) *
                                (v2d[0]*r->rdir[1] - v2d[1]*r->rdir[0]);
                if (vol > .0)
                        ecp->signum = 1;
                else if (vol < .0)
                        ecp->signum = -1;
                else
                        ecp->signum = 0;
                ecp->v1i = v1;
                ecp->v2i = v2;
        }
        return(reversed ? -ecp->signum : ecp->signum);
}


static void
mesh_hit(oset, r)               /* intersect ray with mesh triangle(s) */
OBJECT  *oset;
RAY     *r;
{
        int32           tvi[3];
        int             sv1, sv2, sv3;
        MESHVERT        tv[3];
        OBJECT          tmod;
        FVECT           va, vb, nrm;
        double          d;
        int             i;
                                        /* check each triangle */
        for (i = oset[0]; i > 0; i--) {
                if (!getmeshtrivid(tvi, &tmod, curmsh, oset[i]))
                        objerror(edge_cache.o, INTERNAL,
                                "missing triangle vertices in mesh_hit");
                sv1 = signed_volume(r, tvi[0], tvi[1]);
                sv2 = signed_volume(r, tvi[1], tvi[2]);
                sv3 = signed_volume(r, tvi[2], tvi[0]);
                if (sv1 != sv2 || sv2 != sv3)   /* compare volume signs */
                        if (sv1 && sv2 && sv3)
                                continue;
                                                /* compute intersection */
                getmeshvert(&tv[0], curmsh, tvi[0], MT_V);
                getmeshvert(&tv[1], curmsh, tvi[1], MT_V);
                getmeshvert(&tv[2], curmsh, tvi[2], MT_V);
                VSUB(va, tv[0].v, tv[2].v);
                VSUB(vb, tv[1].v, tv[0].v);
                VCROSS(nrm, va, vb);
                d = DOT(r->rdir, nrm);
                if (d == 0.0)
                        continue;               /* ray is tangent */
                VSUB(va, tv[0].v, r->rorg);
                d = DOT(va, nrm) / d;
                if (d <= FTINY || d >= r->rot)
                        continue;               /* not good enough */
                r->robj = oset[i];              /* else record hit */
                r->ro = edge_cache.o;
                r->rot = d;
                VSUM(r->rop, r->rorg, r->rdir, d);
                VCOPY(r->ron, nrm);
                /* normalize(r->ron) called & r->rod set in o_mesh() */
        }
}


extern int
o_mesh(                 /* compute ray intersection with a mesh */
        OBJREC          *o,
        register RAY    *r
)
{
        RAY             rcont;
        int             flags;
        MESHVERT        tv[3];
        OBJECT          tmod;
        RREAL           wt[3];
        int             i;
                                        /* get the mesh instance */
        prep_edge_cache(o);
                                        /* copy and transform ray */
        rcont = *r;
        multp3(rcont.rorg, r->rorg, curmi->x.b.xfm);
        multv3(rcont.rdir, r->rdir, curmi->x.b.xfm);
        for (i = 0; i < 3; i++)
                rcont.rdir[i] /= curmi->x.b.sca;
        rcont.rmax *= curmi->x.b.sca;
                                        /* clear and trace ray */
        rayclear(&rcont);
        rcont.hitf = mesh_hit;
        if (!localhit(&rcont, &curmi->msh->mcube))
                return(0);                      /* missed */
        if (rcont.rot * curmi->x.f.sca >= r->rot)
                return(0);                      /* not close enough */
                                        /* transform ray back */
        r->rot = rcont.rot * curmi->x.f.sca;
        multp3(r->rop, rcont.rop, curmi->x.f.xfm);
        multv3(r->ron, rcont.ron, curmi->x.f.xfm);
        normalize(r->ron);
        r->rod = -DOT(r->rdir, r->ron);
                                        /* get triangle */
        flags = getmeshtri(tv, &tmod, curmsh, rcont.robj, MT_ALL);
        if (!(flags & MT_V))
                objerror(o, INTERNAL, "missing mesh vertices in o_mesh");
        r->robj = objndx(o);            /* set object and material */
        if (o->omod == OVOID && tmod != OVOID) {
                r->ro = getmeshpseudo(curmsh, tmod);
                r->rox = &curmi->x;
        } else
                r->ro = o;
                                        /* compute barycentric weights */
        if (flags & (MT_N|MT_UV))
                if (get_baryc(wt, rcont.rop, tv[0].v, tv[1].v, tv[2].v) < 0) {
                        objerror(o, WARNING, "bad triangle in o_mesh");
                        flags &= ~(MT_N|MT_UV);
                }
        if (flags & MT_N) {             /* interpolate normal */
                for (i = 0; i < 3; i++)
                        rcont.pert[i] = wt[0]*tv[0].n[i] +
                                        wt[1]*tv[1].n[i] +
                                        wt[2]*tv[2].n[i];
                multv3(r->pert, rcont.pert, curmi->x.f.xfm);
                if (normalize(r->pert) != 0.0)
                        for (i = 0; i < 3; i++)
                                r->pert[i] -= r->ron[i];
        } else
                r->pert[0] = r->pert[1] = r->pert[2] = .0;

        if (flags & MT_UV)              /* interpolate uv coordinates */
                for (i = 0; i < 2; i++)
                        r->uv[i] = wt[0]*tv[0].uv[i] +
                                        wt[1]*tv[1].uv[i] +
                                        wt[2]*tv[2].uv[i];
        else
                r->uv[0] = r->uv[1] = .0;

                                        /* return hit */
        return(1);
}
Revision:	2.10
Committed:	Tue Mar 30 16:13:01 2004 UTC (20 years, 1 month ago) by schorsch
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad3R7P2, rad3R7P1, rad3R6, rad3R6P1
Changes since 2.9:	+7 -5 lines
Log Message:	Continued ANSIfication. There are only bits and pieces left now.
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	schorsch	2.10	static const char RCSid[] = "$Id: o_mesh.c,v 2.9 2003/07/21 22:30:19 schorsch Exp $";
3	greg	2.1	#endif
4			/*
5			* Routines for computing ray intersections with meshes.
6			*
7			* Intersection with a triangle mesh is based on Segura and Feito's
8			* WSCG 2001 paper, "Algorithms to Test Ray-Triangle Intersection,
9			* Comparative Study." This method avoids additional storage
10			* requirements, floating divides, and allows some savings by
11			* caching ray-edge comparisons that are otherwise repeated locally
12			* in typical mesh geometries. (This is our own optimization.)
13			*
14			* The code herein is quite similar to that in o_instance.c, the
15			* chief differences being the custom triangle intersection routines
16			* and the fact that an "OBJECT" in the mesh octree is not an index
17			* into the Radiance OBJREC list, but a mesh triangle index. We still
18			* utilize the standard octree traversal code by setting the hitf
19			* function pointer in the RAY struct to our custom mesh_hit() call.
20			*/
21
22			#include "copyright.h"
23
24	schorsch	2.8	#include <string.h>
25
26	greg	2.1	#include "ray.h"
27			#include "mesh.h"
28			#include "tmesh.h"
29	schorsch	2.10	#include "rtotypes.h"
30	greg	2.1
31
32	greg	2.2	#define EDGE_CACHE_SIZ 251 /* length of mesh edge cache */
33	greg	2.1
34			#define curmi (edge_cache.mi)
35			#define curmsh (curmi->msh)
36
37
38			/* Cache of signed volumes for this ray and this mesh */
39			struct EdgeCache {
40			OBJREC o; / mesh object */
41			MESHINST mi; / current mesh instance */
42			struct EdgeSide {
43	greg	2.6	int32 v1i, v2i; /* vertex indices (lowest first) */
44	greg	2.1	short signum; /* signed volume */
45			} cache[EDGE_CACHE_SIZ];
46			} edge_cache;
47
48
49			static void
50			prep_edge_cache(o) /* get instance and clear edge cache */
51			OBJREC *o;
52			{
53			/* get mesh instance */
54			edge_cache.mi = getmeshinst(edge_cache.o = o, IO_ALL);
55			/* clear edge cache */
56	schorsch	2.8	memset((void *)edge_cache.cache, '\0', sizeof(edge_cache.cache));
57	greg	2.1	}
58
59
60			static int
61			signed_volume(r, v1, v2) /* get signed volume for ray and edge */
62			register RAY *r;
63	greg	2.6	int32 v1, v2;
64	greg	2.1	{
65			int reversed = 0;
66			register struct EdgeSide *ecp;
67
68			if (v1 > v2) {
69	greg	2.6	int32 t = v2; v2 = v1; v1 = t;
70	greg	2.1	reversed = 1;
71			}
72	greg	2.2	ecp = &edge_cache.cache[((v2<<11 ^ v1) & 0x7fffffff) % EDGE_CACHE_SIZ];
73	greg	2.1	if (ecp->v1i != v1 \|\| ecp->v2i != v2) {
74			MESHVERT tv1, tv2; /* compute signed volume */
75	greg	2.4	FVECT v2d;
76	greg	2.1	double vol;
77			if (!getmeshvert(&tv1, edge_cache.mi->msh, v1, MT_V) \|\|
78	greg	2.2	!getmeshvert(&tv2, edge_cache.mi->msh, v2, MT_V))
79	greg	2.1	objerror(edge_cache.o, INTERNAL,
80	greg	2.2	"missing mesh vertex in signed_volume");
81	greg	2.4	VSUB(v2d, tv2.v, r->rorg);
82	greg	2.1	vol = (tv1.v[0] - r->rorg[0]) *
83	greg	2.4	(v2d[1]r->rdir[2] - v2d[2]r->rdir[1]);
84	greg	2.1	vol += (tv1.v[1] - r->rorg[1]) *
85	greg	2.4	(v2d[2]r->rdir[0] - v2d[0]r->rdir[2]);
86	greg	2.1	vol += (tv1.v[2] - r->rorg[2]) *
87	greg	2.4	(v2d[0]r->rdir[1] - v2d[1]r->rdir[0]);
88	greg	2.1	if (vol > .0)
89			ecp->signum = 1;
90			else if (vol < .0)
91			ecp->signum = -1;
92			else
93			ecp->signum = 0;
94			ecp->v1i = v1;
95			ecp->v2i = v2;
96			}
97			return(reversed ? -ecp->signum : ecp->signum);
98			}
99
100
101			static void
102			mesh_hit(oset, r) /* intersect ray with mesh triangle(s) */
103			OBJECT *oset;
104			RAY *r;
105			{
106	greg	2.6	int32 tvi[3];
107	greg	2.1	int sv1, sv2, sv3;
108			MESHVERT tv[3];
109	greg	2.5	OBJECT tmod;
110	greg	2.1	FVECT va, vb, nrm;
111			double d;
112			int i;
113			/* check each triangle */
114			for (i = oset[0]; i > 0; i--) {
115	greg	2.5	if (!getmeshtrivid(tvi, &tmod, curmsh, oset[i]))
116	greg	2.1	objerror(edge_cache.o, INTERNAL,
117	greg	2.2	"missing triangle vertices in mesh_hit");
118	greg	2.1	sv1 = signed_volume(r, tvi[0], tvi[1]);
119			sv2 = signed_volume(r, tvi[1], tvi[2]);
120			sv3 = signed_volume(r, tvi[2], tvi[0]);
121			if (sv1 != sv2 \|\| sv2 != sv3) /* compare volume signs */
122			if (sv1 && sv2 && sv3)
123			continue;
124			/* compute intersection */
125			getmeshvert(&tv[0], curmsh, tvi[0], MT_V);
126			getmeshvert(&tv[1], curmsh, tvi[1], MT_V);
127			getmeshvert(&tv[2], curmsh, tvi[2], MT_V);
128			VSUB(va, tv[0].v, tv[2].v);
129			VSUB(vb, tv[1].v, tv[0].v);
130			VCROSS(nrm, va, vb);
131			d = DOT(r->rdir, nrm);
132			if (d == 0.0)
133			continue; /* ray is tangent */
134			VSUB(va, tv[0].v, r->rorg);
135			d = DOT(va, nrm) / d;
136			if (d <= FTINY \|\| d >= r->rot)
137			continue; /* not good enough */
138			r->robj = oset[i]; /* else record hit */
139			r->ro = edge_cache.o;
140			r->rot = d;
141			VSUM(r->rop, r->rorg, r->rdir, d);
142			VCOPY(r->ron, nrm);
143	greg	2.2	/* normalize(r->ron) called & r->rod set in o_mesh() */
144	greg	2.1	}
145			}
146
147
148	schorsch	2.10	extern int
149			o_mesh( /* compute ray intersection with a mesh */
150			OBJREC *o,
151			register RAY *r
152			)
153	greg	2.1	{
154			RAY rcont;
155			int flags;
156			MESHVERT tv[3];
157	greg	2.5	OBJECT tmod;
158	schorsch	2.7	RREAL wt[3];
159	greg	2.1	int i;
160			/* get the mesh instance */
161			prep_edge_cache(o);
162			/* copy and transform ray */
163	schorsch	2.9	rcont = *r;
164	greg	2.1	multp3(rcont.rorg, r->rorg, curmi->x.b.xfm);
165			multv3(rcont.rdir, r->rdir, curmi->x.b.xfm);
166			for (i = 0; i < 3; i++)
167			rcont.rdir[i] /= curmi->x.b.sca;
168			rcont.rmax *= curmi->x.b.sca;
169			/* clear and trace ray */
170			rayclear(&rcont);
171			rcont.hitf = mesh_hit;
172			if (!localhit(&rcont, &curmi->msh->mcube))
173			return(0); /* missed */
174			if (rcont.rot * curmi->x.f.sca >= r->rot)
175			return(0); /* not close enough */
176			/* transform ray back */
177			r->rot = rcont.rot * curmi->x.f.sca;
178			multp3(r->rop, rcont.rop, curmi->x.f.xfm);
179			multv3(r->ron, rcont.ron, curmi->x.f.xfm);
180			normalize(r->ron);
181			r->rod = -DOT(r->rdir, r->ron);
182	greg	2.5	/* get triangle */
183			flags = getmeshtri(tv, &tmod, curmsh, rcont.robj, MT_ALL);
184	greg	2.1	if (!(flags & MT_V))
185			objerror(o, INTERNAL, "missing mesh vertices in o_mesh");
186	greg	2.5	r->robj = objndx(o); /* set object and material */
187			if (o->omod == OVOID && tmod != OVOID) {
188			r->ro = getmeshpseudo(curmsh, tmod);
189			r->rox = &curmi->x;
190			} else
191			r->ro = o;
192			/* compute barycentric weights */
193	greg	2.1	if (flags & (MT_N\|MT_UV))
194			if (get_baryc(wt, rcont.rop, tv[0].v, tv[1].v, tv[2].v) < 0) {
195			objerror(o, WARNING, "bad triangle in o_mesh");
196			flags &= ~(MT_N\|MT_UV);
197			}
198			if (flags & MT_N) { /* interpolate normal */
199			for (i = 0; i < 3; i++)
200			rcont.pert[i] = wt[0]*tv[0].n[i] +
201			wt[1]*tv[1].n[i] +
202			wt[2]*tv[2].n[i];
203			multv3(r->pert, rcont.pert, curmi->x.f.xfm);
204			if (normalize(r->pert) != 0.0)
205			for (i = 0; i < 3; i++)
206			r->pert[i] -= r->ron[i];
207	greg	2.3	} else
208			r->pert[0] = r->pert[1] = r->pert[2] = .0;
209
210	greg	2.1	if (flags & MT_UV) /* interpolate uv coordinates */
211			for (i = 0; i < 2; i++)
212			r->uv[i] = wt[0]*tv[0].uv[i] +
213			wt[1]*tv[1].uv[i] +
214			wt[2]*tv[2].uv[i];
215	greg	2.3	else
216			r->uv[0] = r->uv[1] = .0;
217	greg	2.1
218			/* return hit */
219			return(1);
220			}