src/rt/o_mesh.c

#ifndef lint
static const char RCSid[] = "$Id: o_mesh.c,v 2.10 2004/03/30 16:13:01 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
volume_sign(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 volume_sign");
                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]);
                                                /* don't generate 0 */
                ecp->signum = vol > .0 ? 1 : -1;
                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 = volume_sign(r, tvi[0], tvi[1]);
                sv2 = volume_sign(r, tvi[1], tvi[2]);
                sv3 = volume_sign(r, tvi[2], tvi[0]);
                                                /* compare volume signs */
                if ((sv1 != sv2) | (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.11
Committed:	Sat Feb 25 19:49:16 2006 UTC (18 years, 2 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad4R0, rad3R8, rad3R9
Changes since 2.10:	+13 -17 lines
Log Message:	Minor optimizations in ray-triangle intersection test
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.11	static const char RCSid[] = "$Id: o_mesh.c,v 2.10 2004/03/30 16:13:01 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	greg	2.11	volume_sign(r, v1, v2) /* get signed volume for ray and edge */
62	greg	2.1	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.11	if ((ecp->v1i != v1) \| (ecp->v2i != v2)) {
74	greg	2.1	MESHVERT tv1, tv2; /* compute signed volume */
75	greg	2.4	FVECT v2d;
76	greg	2.1	double vol;
77	greg	2.11	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.11	"missing mesh vertex in volume_sign");
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.11	/* don't generate 0 */
89			ecp->signum = vol > .0 ? 1 : -1;
90	greg	2.1	ecp->v1i = v1;
91			ecp->v2i = v2;
92			}
93			return(reversed ? -ecp->signum : ecp->signum);
94			}
95
96
97			static void
98			mesh_hit(oset, r) /* intersect ray with mesh triangle(s) */
99			OBJECT *oset;
100			RAY *r;
101			{
102	greg	2.6	int32 tvi[3];
103	greg	2.1	int sv1, sv2, sv3;
104			MESHVERT tv[3];
105	greg	2.5	OBJECT tmod;
106	greg	2.1	FVECT va, vb, nrm;
107			double d;
108			int i;
109			/* check each triangle */
110			for (i = oset[0]; i > 0; i--) {
111	greg	2.5	if (!getmeshtrivid(tvi, &tmod, curmsh, oset[i]))
112	greg	2.1	objerror(edge_cache.o, INTERNAL,
113	greg	2.2	"missing triangle vertices in mesh_hit");
114	greg	2.11	sv1 = volume_sign(r, tvi[0], tvi[1]);
115			sv2 = volume_sign(r, tvi[1], tvi[2]);
116			sv3 = volume_sign(r, tvi[2], tvi[0]);
117			/* compare volume signs */
118			if ((sv1 != sv2) \| (sv2 != sv3))
119			continue;
120	greg	2.1	/* compute intersection */
121			getmeshvert(&tv[0], curmsh, tvi[0], MT_V);
122			getmeshvert(&tv[1], curmsh, tvi[1], MT_V);
123			getmeshvert(&tv[2], curmsh, tvi[2], MT_V);
124			VSUB(va, tv[0].v, tv[2].v);
125			VSUB(vb, tv[1].v, tv[0].v);
126			VCROSS(nrm, va, vb);
127			d = DOT(r->rdir, nrm);
128			if (d == 0.0)
129			continue; /* ray is tangent */
130			VSUB(va, tv[0].v, r->rorg);
131			d = DOT(va, nrm) / d;
132			if (d <= FTINY \|\| d >= r->rot)
133			continue; /* not good enough */
134			r->robj = oset[i]; /* else record hit */
135			r->ro = edge_cache.o;
136			r->rot = d;
137			VSUM(r->rop, r->rorg, r->rdir, d);
138			VCOPY(r->ron, nrm);
139	greg	2.2	/* normalize(r->ron) called & r->rod set in o_mesh() */
140	greg	2.1	}
141			}
142
143
144	schorsch	2.10	extern int
145			o_mesh( /* compute ray intersection with a mesh */
146			OBJREC *o,
147			register RAY *r
148			)
149	greg	2.1	{
150			RAY rcont;
151			int flags;
152			MESHVERT tv[3];
153	greg	2.5	OBJECT tmod;
154	schorsch	2.7	RREAL wt[3];
155	greg	2.1	int i;
156			/* get the mesh instance */
157			prep_edge_cache(o);
158			/* copy and transform ray */
159	schorsch	2.9	rcont = *r;
160	greg	2.1	multp3(rcont.rorg, r->rorg, curmi->x.b.xfm);
161			multv3(rcont.rdir, r->rdir, curmi->x.b.xfm);
162			for (i = 0; i < 3; i++)
163			rcont.rdir[i] /= curmi->x.b.sca;
164			rcont.rmax *= curmi->x.b.sca;
165			/* clear and trace ray */
166			rayclear(&rcont);
167			rcont.hitf = mesh_hit;
168			if (!localhit(&rcont, &curmi->msh->mcube))
169			return(0); /* missed */
170			if (rcont.rot * curmi->x.f.sca >= r->rot)
171			return(0); /* not close enough */
172			/* transform ray back */
173			r->rot = rcont.rot * curmi->x.f.sca;
174			multp3(r->rop, rcont.rop, curmi->x.f.xfm);
175			multv3(r->ron, rcont.ron, curmi->x.f.xfm);
176			normalize(r->ron);
177			r->rod = -DOT(r->rdir, r->ron);
178	greg	2.5	/* get triangle */
179			flags = getmeshtri(tv, &tmod, curmsh, rcont.robj, MT_ALL);
180	greg	2.1	if (!(flags & MT_V))
181			objerror(o, INTERNAL, "missing mesh vertices in o_mesh");
182	greg	2.5	r->robj = objndx(o); /* set object and material */
183			if (o->omod == OVOID && tmod != OVOID) {
184			r->ro = getmeshpseudo(curmsh, tmod);
185			r->rox = &curmi->x;
186			} else
187			r->ro = o;
188			/* compute barycentric weights */
189	greg	2.1	if (flags & (MT_N\|MT_UV))
190			if (get_baryc(wt, rcont.rop, tv[0].v, tv[1].v, tv[2].v) < 0) {
191			objerror(o, WARNING, "bad triangle in o_mesh");
192			flags &= ~(MT_N\|MT_UV);
193			}
194			if (flags & MT_N) { /* interpolate normal */
195			for (i = 0; i < 3; i++)
196			rcont.pert[i] = wt[0]*tv[0].n[i] +
197			wt[1]*tv[1].n[i] +
198			wt[2]*tv[2].n[i];
199			multv3(r->pert, rcont.pert, curmi->x.f.xfm);
200			if (normalize(r->pert) != 0.0)
201			for (i = 0; i < 3; i++)
202			r->pert[i] -= r->ron[i];
203	greg	2.3	} else
204			r->pert[0] = r->pert[1] = r->pert[2] = .0;
205
206	greg	2.1	if (flags & MT_UV) /* interpolate uv coordinates */
207			for (i = 0; i < 2; i++)
208			r->uv[i] = wt[0]*tv[0].uv[i] +
209			wt[1]*tv[1].uv[i] +
210			wt[2]*tv[2].uv[i];
211	greg	2.3	else
212			r->uv[0] = r->uv[1] = .0;
213	greg	2.1
214			/* return hit */
215			return(1);
216			}