--- ray/src/hd/rholo2.c 1997/12/18 09:33:35 3.10 +++ ray/src/hd/rholo2.c 1998/12/03 15:21:05 3.20 @@ -1,4 +1,4 @@ -/* Copyright (c) 1997 Silicon Graphics, Inc. */ +/* Copyright (c) 1998 Silicon Graphics, Inc. */ #ifndef lint static char SCCSid[] = "$SunId$ SGI"; @@ -13,30 +13,172 @@ static char SCCSid[] = "$SunId$ SGI"; #include "random.h" +VIEWPOINT myeye; /* target view position */ + +struct gclim { + HOLO *hp; /* holodeck pointer */ + GCOORD gc; /* grid cell */ + FVECT egp; /* eye grid point */ + double erg2; /* mean square eye grid range */ + double gmin[2], gmax[2]; /* grid coordinate limits */ +}; /* a grid coordinate range */ + + +static +initeyelim(gcl, hp, gc) /* initialize grid coordinate limits */ +register struct gclim *gcl; +register HOLO *hp; +GCOORD *gc; +{ + register FLOAT *v; + register int i; + + if (hp != NULL) { + hdgrid(gcl->egp, gcl->hp = hp, myeye.vpt); + gcl->erg2 = 0; + for (i = 0, v = hp->wg[0]; i < 3; i++, v += 3) + gcl->erg2 += DOT(v,v); + gcl->erg2 *= (1./3.) * myeye.rng*myeye.rng; + } + if (gc != NULL) + copystruct(&gcl->gc, gc); + gcl->gmin[0] = gcl->gmin[1] = FHUGE; + gcl->gmax[0] = gcl->gmax[1] = -FHUGE; +} + + +static +groweyelim(gcl, gc, r0, r1) /* grow grid limits about eye point */ +register struct gclim *gcl; +GCOORD *gc; +double r0, r1; +{ + FVECT gp, ab; + double vlen, plen, dv0, dv1; + double rd2, dwall, gpos; + int eyeout; + register int i, g0, g1; + + i = gc->w>>1; + if (gc->w&1) + eyeout = (gp[i] = gcl->hp->grid[i]) < gcl->egp[i]; + else + eyeout = (gp[i] = 0) > gcl->egp[i]; + gp[hdwg0[gc->w]] = gc->i[0] + r0; + gp[hdwg1[gc->w]] = gc->i[1] + r1; + VSUB(ab, gcl->egp, gp); + rd2 = DOT(ab,ab); + if (rd2 <= gcl->erg2) { + gcl->gmin[0] = gcl->gmin[1] = -FHUGE; + gcl->gmax[0] = gcl->gmax[1] = FHUGE; + return; + } + rd2 = gcl->erg2 / rd2; + vlen = 1. - rd2; + plen = sqrt(rd2 * vlen); + g0 = gcl->gc.w>>1; + dwall = (gcl->gc.w&1 ? gcl->hp->grid[g0] : 0) - gp[g0]; + for (i = 0; i < 4; i++) { + if (i == 2) + plen = -plen; + g1 = (g0+(i&1)+1)%3; + dv0 = vlen*ab[g0] + plen*ab[g1]; + dv1 = vlen*ab[g1] - plen*ab[g0]; + if ((dv0 < 0 ^ dwall < 0 ^ eyeout) || + (dv0 <= FTINY && dv0 >= -FTINY)) { + if (eyeout) + dv1 = -dv1; + if (dv1 > FTINY) + gcl->gmax[i&1] = FHUGE; + else if (dv1 < -FTINY) + gcl->gmin[i&1] = -FHUGE; + } else { + gpos = gp[g1] + dv1*dwall/dv0; + if (gpos < gcl->gmin[i&1]) + gcl->gmin[i&1] = gpos; + if (gpos > gcl->gmax[i&1]) + gcl->gmax[i&1] = gpos; + } + } +} + + +static int +clipeyelim(rrng, gcl) /* clip eye limits to grid cell */ +register short rrng[2][2]; +register struct gclim *gcl; +{ + int incell = 1; + register int i; + + for (i = 0; i < 2; i++) { + if (gcl->gmin[i] < gcl->gc.i[i]) + gcl->gmin[i] = gcl->gc.i[i]; + if (gcl->gmax[i] > gcl->gc.i[i]+1) + gcl->gmax[i] = gcl->gc.i[i]+1; + if (gcl->gmax[i] > gcl->gmin[i]) { + rrng[i][0] = 256.*(gcl->gmin[i] - gcl->gc.i[i]) + + (1.-FTINY); + rrng[i][1] = 256.*(gcl->gmax[i] - gcl->gc.i[i]) + + (1.-FTINY) - rrng[i][0]; + } else + rrng[i][0] = rrng[i][1] = 0; + incell &= rrng[i][1] > 0; + } + return(incell); +} + + packrays(rod, p) /* pack ray origins and directions */ register float *rod; register PACKET *p; { - static FVECT ro, rd; + int nretries = p->nr + 2; + struct gclim eyelim; + short rrng0[2][2], rrng1[2][2]; + int useyelim; GCOORD gc[2]; - int ila[2], hsh; - double d, sl[4]; + FVECT ro, rd; + double d; register int i; if (!hdbcoord(gc, hdlist[p->hd], p->bi)) error(CONSISTENCY, "bad beam index in packrays"); - ila[0] = p->hd; ila[1] = p->bi; - hsh = ilhash(ila,2) + p->nc; + if ((useyelim = myeye.rng > FTINY)) { + initeyelim(&eyelim, hdlist[p->hd], gc); + groweyelim(&eyelim, gc+1, 0., 0.); + groweyelim(&eyelim, gc+1, 1., 1.); + useyelim &= clipeyelim(rrng0, &eyelim); + } for (i = 0; i < p->nr; i++) { - multisamp(sl, 4, urand(hsh+i)); - p->ra[i].r[0][0] = sl[0] * 256.; - p->ra[i].r[0][1] = sl[1] * 256.; - p->ra[i].r[1][0] = sl[2] * 256.; - p->ra[i].r[1][1] = sl[3] * 256.; + retry: + if (useyelim) { + initeyelim(&eyelim, NULL, gc+1); + p->ra[i].r[0][0] = (int)(frandom()*rrng0[0][1]) + + rrng0[0][0]; + p->ra[i].r[0][1] = (int)(frandom()*rrng0[1][1]) + + rrng0[1][0]; + groweyelim(&eyelim, gc, + (1./256.)*(p->ra[i].r[0][0]+.5), + (1./256.)*(p->ra[i].r[0][1]+.5)); + if (!clipeyelim(rrng1, &eyelim)) { + useyelim &= nretries-- > 0; + goto retry; + } + p->ra[i].r[1][0] = (int)(frandom()*rrng1[0][1]) + + rrng1[0][0]; + p->ra[i].r[1][1] = (int)(frandom()*rrng1[1][1]) + + rrng1[1][0]; + } else { + p->ra[i].r[0][0] = frandom() * 256.; + p->ra[i].r[0][1] = frandom() * 256.; + p->ra[i].r[1][0] = frandom() * 256.; + p->ra[i].r[1][1] = frandom() * 256.; + } d = hdray(ro, rd, hdlist[p->hd], gc, p->ra[i].r); - if (!vdef(OBSTRUCTIONS)) - d *= frandom(); /* random offset */ if (p->offset != NULL) { + if (!vdef(OBSTRUCTIONS)) + d *= frandom(); /* random offset */ VSUM(ro, ro, rd, d); /* advance ray */ p->offset[i] = d; } @@ -115,7 +257,7 @@ new_rtrace() /* restart rtrace calculation */ getradfile() /* run rad and get needed variables */ { - static short mvar[] = {OCTREE,-1}; + static short mvar[] = {OCTREE,EYESEP,-1}; static char tf1[] = TEMPLATE; char tf2[64]; char combuf[256]; @@ -124,7 +266,7 @@ getradfile() /* run rad and get needed variables */ register char *cp; /* check if rad file specified */ if (!vdef(RIF)) - return; + return(0); /* create rad command */ mktemp(tf1); sprintf(tf2, "%s.rif", tf1); @@ -147,17 +289,17 @@ getradfile() /* run rad and get needed variables */ strcpy(pippt, "> /dev/null"); /* nothing to match */ else sprintf(cp, ")[ \t]*=' > %s", tf2); - if (system(combuf)) { - unlink(tf2); /* clean up */ - unlink(tf1); - error(USER, "error executing rad command"); - } +#ifdef DEBUG + wputs(combuf); wputs("\n"); +#endif + system(combuf); /* ignore exit code */ if (pippt == NULL) { loadvars(tf2); /* load variables */ unlink(tf2); } rtargc += wordfile(rtargv+rtargc, tf1); /* get rtrace options */ unlink(tf1); /* clean up */ + return(1); }