/************************************************************************ * * * Copyright (c) 1988, David B. Wecker * * All Rights Reserved * * * * This file is part of DBW_uRAY * * * * DBW_uRAY is distributed in the hope that it will be useful, but * * WITHOUT ANY WARRANTY. No author or distributor accepts * * responsibility to anyone for the consequences of using it or for * * whether it serves any particular purpose or works at all, unless * * he says so in writing. Refer to the DBW_uRAY General Public * * License for full details. * * * * Everyone is granted permission to copy, modify and redistribute * * DBW_uRAY, but only under the conditions described in the * * DBW_uRAY General Public License. A copy of this license is * * supposed to have been given to you along with DBW_uRAY so you * * can know your rights and responsibilities. It should be in a file * * named COPYING. Among other things, the copyright notice and this * * notice must be preserved on all copies. * ************************************************************************ * * * Authors: * * DBW - David B. Wecker * * * * Versions: * * V1.0 881023 DBW - First released version * * V1.1 881110 DBW - Fixed scan coherence code * * V1.2 881125 DBW - Removed ALL scan coherence code (useless) * * added "fat" extent boxes * * * ************************************************************************/ #include "uray.h" /*************************************************************************/ /******** This module computes all ray object intersections **************/ /*************************************************************************/ /* see if we hit a bounding box (n1 = near side, f1 = far side if good) */ /* near point = from + n1 * direction */ /* far point = from + f1 * direction */ NODE *hitextent(s,n,f) EXTENT *s; FLTDBL *n,*f; { /* Use registers to be as efficient as possible */ register FLTDBL n1 = -MYHUGE,f1 = MYHUGE,n2,f2,tval,sval; register int i; /* check each dimension of the bounding box */ for (i=0; i<3; i++) { /* see which side of the box that we're on */ if ((tval=dirinv[i]) < 0.0) { n2 = (s->max[i] - (sval=from[i])) * tval; f2 = (s->min[i] - sval) * tval; } else { n2 = (s->min[i] - (sval=from[i])) * tval; f2 = (s->max[i] - sval) * tval; } /* if new near is worse than the old near (and old far) then bad */ if (n2 > n1 && (n1=n2) > f1) return NULL; /* if new far is worse than old far and (and old near) or if we're too close to 0 then return bad */ if (f2 < f1 && ((f1=f2) < n1 || f1 < TOL)) return NULL; } /* are we worse than the far we ENTERED the routine with? */ if (n1 > *f) return NULL; /* set the return values */ *n = n1; *f = f1; /* return the object hit */ return (NODE *)s; } /* see if we hit the sphere in question */ NODE *hitsphere(s,n1) SPHERE *s; FLTDBL *n1; { register FLTDBL r_r,d_r,radical,t1,t2; VEC r; vcomb(-1.0, from, s->cen, r); r_r = vdot(r,r); d_r = vdot(direction,r); if ((radical = (d_r*d_r) + s->rad - r_r) < 0.0) return NULL; radical = sqrt(radical); if (d_r < radical) { t1 = d_r + radical; t2 = d_r - radical; } else { t1 = d_r - radical; t2 = d_r + radical; } /* since there are 2 possible intersection points, pick the best */ if (t1 <= TOL) t1 = t2; /* are we too close to the intersection point */ if (t1 <= TOL) return NULL; /* are we further than the bounding box intersection */ if (t1 > *n1) return NULL; /* return the intersection point */ *n1 = t1; /* return the object that we hit */ return (NODE *)s; } /* take care of all planar intersections (quad, triangle, ring) */ NODE *hitplane(q,n1) RING *q; FLTDBL *n1; { VEC sfs,h; FLTDBL det,r; /* First get the determinant of the interseection plane */ det = q->v1[0] * ((q->v2[1] * direction[2])-(q->v2[2] * direction[1])); det -= q->v2[0] * ((q->v1[1] * direction[2])-(q->v1[2] * direction[1])); det += direction[0] * ((q->v1[1] * q->v2[2]) -(q->v1[2] * q->v2[1])); /* 0 determinant says that we missed the plane of intersection */ if (det == 0.0) return NULL; vcomb(-1.0,q->p0,from,h); /* now build the exact intersection point */ sfs[2] = h[0] * ((q->v1[1] * q->v2[2]) - (q->v1[2] * q->v2[1])); sfs[2] -= h[1] * ((q->v1[0] * q->v2[2]) - (q->v1[2] * q->v2[0])); sfs[2] += h[2] * ((q->v1[0] * q->v2[1]) - (q->v1[1] * q->v2[0])); sfs[2] /= det; sfs[2] = -sfs[2]; /* Too close to 0? */ if (sfs[2] <= TOL) return NULL; /* Further than bounding box? */ if (sfs[2] > *n1) return NULL; /* Now build comparison points for each object type */ sfs[0] = h[0] * ((q->v2[1] * direction[2]) - (q->v2[2] * direction[1])); sfs[0] -= h[1] * ((q->v2[0] * direction[2]) - (q->v2[2] * direction[0])); sfs[0] += h[2] * ((q->v2[0] * direction[1]) - (q->v2[1] * direction[0])); sfs[0] /= det; sfs[1] = h[0] * ((q->v1[2] * direction[1]) - (q->v1[1] * direction[2])); sfs[1] += h[1] * ((q->v1[0] * direction[2]) - (q->v1[2] * direction[0])); sfs[1] -= h[2] * ((q->v1[0] * direction[1]) - (q->v1[1] * direction[0])); sfs[1] /= det; switch (q->typ) { /* Make sure that we aren't off the triangles surface */ case TYP_T: if (sfs[0] < 0.0 || sfs[1] < 0.0 || sfs[0]+sfs[1] > 1.0) return NULL; break; /* make sure that we're within the bounds of the quad */ case TYP_Q: if (sfs[0] < 0.0 || sfs[1] < 0.0 || sfs[0] > 1.0 || sfs[1] > 1.0) return NULL; break; /* make sure that we're between the ring's radii */ case TYP_R: r = sfs[0] * sfs[0] + sfs[1] * sfs[1]; if (r > q->rad2 || r < q->rad1) return NULL; break; } /* return the intersection point */ *n1 = sfs[2]; return (NODE *)q; } /* recursive part of the intersection() routine */ NODE *intersection2(cur,n1,f1) NODE *cur; FLTDBL *n1,*f1; { register NODE *best = NULL,*hit; register int i; FLTDBL ln1,lf1,ln2,lf2; /* save copies of the near and far intersection points (for extents) */ ln1 = *n1; lf1 = *f1; /* check each object in the hierarchy */ while (cur) { /* re-init the intersection points */ hit = NULL; ln2 = ln1; lf2 = lf1; /* do the specific object */ switch (cur->typ) { /* if we hit an extent, call ourselves recursively */ case TYP_E: if (hitextent(cur,&ln2,&lf2)) { gehits++; hit = intersection2(((EXTENT *)cur)->child,&ln2,&lf2); } else fehits++; break; /* if we hit a sphere, save the intersection info */ case TYP_S: hit = hitsphere(cur,&lf2); ln2 = lf2; break; /* if we hit a planar, save the intersection info */ case TYP_T: case TYP_Q: case TYP_R: hit = hitplane(cur,&lf2); ln2 = lf2; break; } /* if NOT an extent, update the node statistics */ if (cur->typ != TYP_E) { if (hit) gnhits++; else fnhits++; } /* if hit a non extent, then save this is the best so far */ if (hit && hit->typ != TYP_E) { best = hit; ln1 = ln2; lf1 = lf2; } /* go get next object */ cur = cur->next; } /* if we got something... return in */ if (best) { *n1 = ln1; *f1 = lf1; } return best; } /* return closest node (frm + n1 * dir) or NULL */ NODE *intersection(frm, dir, n1) VEC frm, dir; FLTDBL *n1; { register int i; NODE *hit = NULL; FLTDBL f1; /* initialize near and far */ *n1 = -MYHUGE; f1 = MYHUGE; /* set up global vectors FROM, DIRECTION and DIRINV (inverse of dir) */ for (i=0; i<3; i++) { from[i] = frm[i]; direction[i] = dir[i]; if (FABS(direction[i]) < TOL) { if (direction[i] < 0.0) dirinv[i] = -MYHUGE; else dirinv[i] = MYHUGE; } else dirinv[i] = 1.0 / direction[i]; } /* call the recursive part of this routine */ hit = intersection2(nodes,n1,&f1); /* return the result */ return hit; }