--- ray/src/hd/rholo2.c 1998/12/01 15:47:05 3.19 +++ ray/src/hd/rholo2.c 1998/12/03 15:21:05 3.20 @@ -19,67 +19,86 @@ struct gclim { HOLO *hp; /* holodeck pointer */ GCOORD gc; /* grid cell */ FVECT egp; /* eye grid point */ - double erg[2]; /* eye range in wall grid coords */ + double erg2; /* mean square eye grid range */ double gmin[2], gmax[2]; /* grid coordinate limits */ }; /* a grid coordinate range */ static -initeyelim(gcl, hd, gc) /* initialize grid coordinate limits */ +initeyelim(gcl, hp, gc) /* initialize grid coordinate limits */ register struct gclim *gcl; -int hd; +register HOLO *hp; GCOORD *gc; { register FLOAT *v; register int i; - gcl->hp = hdlist[hd]; - copystruct(&gcl->gc, gc); - hdgrid(gcl->egp, gcl->hp, myeye.vpt); - for (i = 0; i < 2; i++) { - v = gcl->hp->wg[((gcl->gc.w>>1)+i+1)%3]; - gcl->erg[i] = myeye.rng * VLEN(v); - gcl->gmin[i] = FHUGE; gcl->gmax[i] = -FHUGE; + 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, gp) /* grow grid limits about eye point */ +groweyelim(gcl, gc, r0, r1) /* grow grid limits about eye point */ register struct gclim *gcl; -FVECT gp; +GCOORD *gc; +double r0, r1; { - FVECT ab; - double l2, d, mult, wg; - register int i, g; + 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); - l2 = DOT(ab,ab); - if (l2 <= gcl->erg[0]*gcl->erg[1]) { + rd2 = DOT(ab,ab); + if (rd2 <= gcl->erg2) { gcl->gmin[0] = gcl->gmin[1] = -FHUGE; gcl->gmax[0] = gcl->gmax[1] = FHUGE; return; } - mult = gp[g = gcl->gc.w>>1]; - if (gcl->gc.w&1) - mult -= gcl->hp->grid[g]; - if (ab[g]*ab[g] > gcl->erg[0]*gcl->erg[1]) - mult /= -ab[g]; - else if (fabs(ab[hdwg0[gcl->gc.w]]) > fabs(ab[hdwg1[gcl->gc.w]])) - mult = (gcl->gc.i[0] + .5 - gp[hdwg0[gcl->gc.w]]) / - ab[hdwg0[gcl->gc.w]]; - else - mult = (gcl->gc.i[1] + .5 - gp[hdwg1[gcl->gc.w]]) / - ab[hdwg1[gcl->gc.w]]; - for (i = 0; i < 2; i++) { - g = ((gcl->gc.w>>1)+i+1)%3; - wg = gp[g] + mult*ab[g]; - d = mult*gcl->erg[i]; - if (d < 0.) d = -d; - if (wg - d < gcl->gmin[i]) - gcl->gmin[i] = wg - d; - if (wg + d > gcl->gmax[i]) - gcl->gmax[i] = wg + d; + 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; + } } } @@ -97,13 +116,14 @@ register struct gclim *gcl; 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 ((incell &= gcl->gmax[i] > gcl->gmin[i])) { + 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]; - incell &= rrng[i][1] > 0; - } + } else + rrng[i][0] = rrng[i][1] = 0; + incell &= rrng[i][1] > 0; } return(incell); } @@ -113,11 +133,6 @@ packrays(rod, p) /* pack ray origins and directions * register float *rod; register PACKET *p; { -#define gp ro -#ifdef DEBUG - double dist2sum = 0.; - FVECT vt; -#endif int nretries = p->nr + 2; struct gclim eyelim; short rrng0[2][2], rrng1[2][2]; @@ -130,38 +145,24 @@ register PACKET *p; if (!hdbcoord(gc, hdlist[p->hd], p->bi)) error(CONSISTENCY, "bad beam index in packrays"); if ((useyelim = myeye.rng > FTINY)) { - initeyelim(&eyelim, p->hd, gc); - gp[gc[1].w>>1] = gc[1].w&1 ? - hdlist[p->hd]->grid[gc[1].w>>1] : 0; - gp[hdwg0[gc[1].w]] = gc[1].i[0]; - gp[hdwg1[gc[1].w]] = gc[1].i[1]; - groweyelim(&eyelim, gp); - gp[hdwg0[gc[1].w]]++; - gp[hdwg1[gc[1].w]]++; - groweyelim(&eyelim, gp); + 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++) { 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]; - initeyelim(&eyelim, p->hd, gc+1); - gp[gc[0].w>>1] = gc[0].w&1 ? - hdlist[p->hd]->grid[gc[0].w>>1] : 0; - gp[hdwg0[gc[0].w]] = gc[0].i[0] + - (1./256.)*(p->ra[i].r[0][0]+.5); - gp[hdwg1[gc[0].w]] = gc[0].i[1] + - (1./256.)*(p->ra[i].r[0][1]+.5); - groweyelim(&eyelim, gp); + 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; -#ifdef DEBUG - if (!useyelim) - error(WARNING, "exceeded retry limit in packrays"); -#endif goto retry; } p->ra[i].r[1][0] = (int)(frandom()*rrng1[0][1]) @@ -175,10 +176,6 @@ register PACKET *p; p->ra[i].r[1][1] = frandom() * 256.; } d = hdray(ro, rd, hdlist[p->hd], gc, p->ra[i].r); -#ifdef DEBUG - VSUM(vt, ro, rd, d); - dist2sum += dist2line(myeye.vpt, ro, vt); -#endif if (p->offset != NULL) { if (!vdef(OBSTRUCTIONS)) d *= frandom(); /* random offset */ @@ -190,10 +187,6 @@ register PACKET *p; VCOPY(rod, rd); rod += 3; } -#ifdef DEBUG - fprintf(stderr, "RMS distance = %f\n", sqrt(dist2sum/p->nr)); -#endif -#undef gp }