Back to index

radiance  4R0+20100331
holo.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: holo.c,v 3.21 2004/01/01 11:21:55 schorsch Exp $";
00003 #endif
00004 /*
00005  * Routines for converting holodeck coordinates, etc.
00006  *
00007  *     10/22/97      GWLarson
00008  */
00009 
00010 #include "holo.h"
00011 
00012 float  hd_depthmap[DCINF-DCLIN];
00013 
00014 int    hdwg0[6] = {1,1,2,2,0,0};
00015 int    hdwg1[6] = {2,2,0,0,1,1};
00016 
00017 static double logstep;
00018 
00019 
00020 extern void
00021 hdcompgrid(                 /* compute derived grid vector and index */
00022        register HOLO *hp
00023 )
00024 {
00025        double d;
00026        register int  i, j;
00027                             /* initialize depth map */
00028        if (hd_depthmap[0] < 1.) {
00029               d = 1. + .5/DCLIN;
00030               for (i = 0; i < DCINF-DCLIN; i++) {
00031                      hd_depthmap[i] = d;
00032                      d *= 1. + 1./DCLIN;
00033               }
00034               logstep = log(1. + 1./DCLIN);
00035        }
00036                             /* compute grid coordinate vectors */
00037        for (i = 0; i < 3; i++) {
00038               fcross(hp->wg[i], hp->xv[(i+1)%3], hp->xv[(i+2)%3]);
00039               d = DOT(hp->wg[i],hp->xv[i]);
00040               if ((d <= FTINY) & (d >= -FTINY))
00041                      error(USER, "degenerate holodeck section");
00042               d = hp->grid[i] / d;
00043               hp->wg[i][0] *= d; hp->wg[i][1] *= d; hp->wg[i][2] *= d;
00044        }
00045                             /* compute linear depth range */
00046        hp->tlin = VLEN(hp->xv[0]) + VLEN(hp->xv[1]) + VLEN(hp->xv[2]);
00047                             /* compute wall super-indices from grid */
00048        hp->wi[0] = 1;              /**** index values begin at 1 ****/
00049        for (i = 1; i < 6; i++) {
00050               hp->wi[i] = 0;
00051               for (j = i; j < 6; j++)
00052                      hp->wi[i] += hp->grid[hdwg0[j]] * hp->grid[hdwg1[j]];
00053               hp->wi[i] *= hp->grid[hdwg0[i-1]] * hp->grid[hdwg1[i-1]];
00054               hp->wi[i] += hp->wi[i-1];
00055        }
00056 }
00057 
00058 
00059 extern int
00060 hdbcoord(            /* compute beam coordinates from index */
00061        GCOORD gc[2],        /* returned */
00062        register HOLO *hp,
00063        register int  i
00064 )
00065 {
00066        register int  j, n;
00067        int    n2, reverse;
00068        GCOORD g2[2];
00069                                    /* check range */
00070        if ((i < 1) | (i > nbeams(hp)))
00071               return(0);
00072        if ( (reverse = i >= hp->wi[5]) )
00073               i -= hp->wi[5] - 1;
00074        for (j = 0; j < 5; j++)            /* find w0 */
00075               if (hp->wi[j+1] > i)
00076                      break;
00077        i -= hp->wi[gc[0].w=j];
00078                                    /* find w1 */
00079        n2 = hp->grid[hdwg0[j]] * hp->grid[hdwg1[j]];
00080        while (++j < 5)      {
00081               n = n2 * hp->grid[hdwg0[j]] * hp->grid[hdwg1[j]];
00082               if (n > i)
00083                      break;
00084               i -= n;
00085        }
00086        gc[1].w = j;
00087                                    /* find position on w0 */
00088        n2 = hp->grid[hdwg0[j]] * hp->grid[hdwg1[j]];
00089        n = i / n2;
00090        gc[0].i[1] = n / hp->grid[hdwg0[gc[0].w]];
00091        gc[0].i[0] = n - gc[0].i[1]*hp->grid[hdwg0[gc[0].w]];
00092        i -= n*n2;
00093                                    /* find position on w1 */
00094        gc[1].i[1] = i / hp->grid[hdwg0[gc[1].w]];
00095        gc[1].i[0] = i - gc[1].i[1]*hp->grid[hdwg0[gc[1].w]];
00096        if (reverse) {
00097               *g2 = *(gc+1);
00098               *(gc+1) = *gc;
00099               *gc = *g2;
00100        }
00101        return(1);                  /* we're done */
00102 }
00103 
00104 
00105 extern int
00106 hdbindex(            /* compute index from beam coordinates */
00107        register HOLO *hp,
00108        register GCOORD      gc[2]
00109 )
00110 {
00111        GCOORD g2[2];
00112        int    reverse;
00113        register int  i, j;
00114                                    /* check ordering and limits */
00115        if ( (reverse = gc[0].w > gc[1].w) ) {
00116               *g2 = *(gc+1);
00117               *(g2+1) = *gc;
00118               gc = g2;
00119        } else if (gc[0].w == gc[1].w)
00120               return(0);
00121        if ((gc[0].w < 0) | (gc[1].w > 5))
00122               return(0);
00123        i = 0;                      /* compute index */
00124        for (j = gc[0].w+1; j < gc[1].w; j++)
00125               i += hp->grid[hdwg0[j]] * hp->grid[hdwg1[j]];
00126        i *= hp->grid[hdwg0[gc[0].w]] * hp->grid[hdwg1[gc[0].w]];
00127        i += hp->wi[gc[0].w];
00128        i += (hp->grid[hdwg0[gc[0].w]]*gc[0].i[1] + gc[0].i[0]) *
00129                      hp->grid[hdwg0[gc[1].w]] * hp->grid[hdwg1[gc[1].w]] ;
00130        i += hp->grid[hdwg0[gc[1].w]]*gc[1].i[1] + gc[1].i[0];
00131        if (reverse)
00132               i += hp->wi[5] - 1;
00133        return(i);
00134 }
00135 
00136 
00137 extern void
00138 hdcell(              /* compute cell coordinates */
00139        register FVECT       cp[4], /* returned (may be passed as FVECT cp[2][2]) */
00140        register HOLO *hp,
00141        register GCOORD      *gc
00142 )
00143 {
00144        register RREAL       *v;
00145        double d;
00146                                    /* compute common component */
00147        VCOPY(cp[0], hp->orig);
00148        if (gc->w & 1) {
00149               v = hp->xv[gc->w>>1];
00150               cp[0][0] += v[0]; cp[0][1] += v[1]; cp[0][2] += v[2];
00151        }
00152        v = hp->xv[hdwg0[gc->w]];
00153        d = (double)gc->i[0] / hp->grid[hdwg0[gc->w]];
00154        VSUM(cp[0], cp[0], v, d);
00155        v = hp->xv[hdwg1[gc->w]];
00156        d = (double)gc->i[1] / hp->grid[hdwg1[gc->w]];
00157        VSUM(cp[0], cp[0], v, d);
00158                                    /* compute x1 sums */
00159        v = hp->xv[hdwg0[gc->w]];
00160        d = 1.0 / hp->grid[hdwg0[gc->w]];
00161        VSUM(cp[1], cp[0], v, d);
00162        VSUM(cp[3], cp[0], v, d);
00163                                    /* compute y1 sums */
00164        v = hp->xv[hdwg1[gc->w]];
00165        d = 1.0 / hp->grid[hdwg1[gc->w]];
00166        VSUM(cp[2], cp[0], v, d);
00167        VSUM(cp[3], cp[3], v, d);
00168 }
00169 
00170 
00171 extern int
00172 hdlseg(                     /* compute line segment for beam */
00173        register int  lseg[2][3],
00174        register HOLO *hp,
00175        GCOORD gc[2]
00176 )
00177 {
00178        register int  k;
00179 
00180        for (k = 0; k < 2; k++) {          /* compute end points */
00181               lseg[k][gc[k].w>>1] = gc[k].w&1    ? hp->grid[gc[k].w>>1]-1 : 0 ;
00182               lseg[k][hdwg0[gc[k].w]] = gc[k].i[0];
00183               lseg[k][hdwg1[gc[k].w]] = gc[k].i[1];
00184        }
00185        return(1);
00186 }
00187 
00188 
00189 extern unsigned int
00190 hdcode(                     /* compute depth code for d */
00191        HOLO   *hp,
00192        double d
00193 )
00194 {
00195        double tl = hp->tlin;
00196        register long c;
00197 
00198        if (d <= 0.)
00199               return(0);
00200        if (d >= .99*FHUGE)
00201               return(DCINF);
00202        if (d < tl)
00203               return((unsigned)(d*DCLIN/tl));
00204        c = (long)(log(d/tl)/logstep) + DCLIN;
00205        return(c > DCINF ? (unsigned)DCINF : (unsigned)c);
00206 }
00207 
00208 
00209 extern void
00210 hdgrid(              /* compute grid coordinates */
00211        FVECT  gp,           /* returned */
00212        register HOLO *hp,
00213        FVECT  wp
00214 )
00215 {
00216        FVECT  vt;
00217 
00218        VSUB(vt, wp, hp->orig);
00219        gp[0] = DOT(vt, hp->wg[0]);
00220        gp[1] = DOT(vt, hp->wg[1]);
00221        gp[2] = DOT(vt, hp->wg[2]);
00222 }
00223 
00224 
00225 extern void
00226 hdworld(             /* compute world coordinates */
00227        register FVECT       wp,
00228        register HOLO *hp,
00229        FVECT  gp
00230 )
00231 {
00232        register double      d;
00233 
00234        d = gp[0]/hp->grid[0];
00235        VSUM(wp, hp->orig, hp->xv[0], d);
00236 
00237        d = gp[1]/hp->grid[1];
00238        VSUM(wp, wp, hp->xv[1], d);
00239 
00240        d = gp[2]/hp->grid[2];
00241        VSUM(wp, wp, hp->xv[2], d);
00242 }
00243 
00244 
00245 extern double
00246 hdray( /* compute ray within a beam */
00247        FVECT  ro,
00248        FVECT  rd,           /* returned */
00249        HOLO   *hp,
00250        GCOORD gc[2],
00251        BYTE   r[2][2]
00252 )
00253 {
00254        FVECT  cp[4], p[2];
00255        register int  i, j;
00256        double d0, d1;
00257                                    /* compute entry and exit points */
00258        for (i = 0; i < 2; i++) {
00259               hdcell(cp, hp, gc+i);
00260               d0 = (1./256.)*(r[i][0]+.5);
00261               d1 = (1./256.)*(r[i][1]+.5);
00262               for (j = 0; j < 3; j++)
00263                      p[i][j] = (1.-d0-d1)*cp[0][j] +
00264                                    d0*cp[1][j] + d1*cp[2][j];
00265        }
00266        VCOPY(ro, p[0]);            /* assign ray origin and direction */
00267        VSUB(rd, p[1], p[0]);
00268        return(normalize(rd));             /* return maximum inside distance */
00269 }
00270 
00271 
00272 extern double
00273 hdinter(      /* compute ray intersection with section */
00274        register GCOORD      gc[2], /* returned */
00275        BYTE   r[2][2],      /* returned (optional) */
00276        double *ed,          /* returned (optional) */
00277        register HOLO *hp,
00278        FVECT  ro,
00279        FVECT  rd            /* normalization of rd affects distances */
00280 )
00281 {
00282        FVECT  p[2], vt;
00283        double d, t0, t1, d0, d1;
00284        register RREAL       *v;
00285        register int  i;
00286                                    /* first, intersect walls */
00287        gc[0].w = gc[1].w = -1;
00288        t0 = -FHUGE; t1 = FHUGE;
00289        VSUB(vt, ro, hp->orig);
00290        for (i = 0; i < 3; i++) {          /* for each wall pair */
00291               d = -DOT(rd, hp->wg[i]);    /* plane distance */
00292               if (d <= FTINY && d >= -FTINY)     /* check for parallel */
00293                      continue;
00294               d1 = DOT(vt, hp->wg[i]);    /* ray distances */
00295               d0 = d1 / d;
00296               d1 = (d1 - hp->grid[i]) / d;
00297               if (d < 0) {                /* check against best */
00298                      if (d0 > t0) {
00299                             t0 = d0;
00300                             gc[0].w = i<<1;
00301                      }
00302                      if (d1 < t1) {
00303                             t1 = d1;
00304                             gc[1].w = i<<1 | 1;
00305                      }
00306               } else {
00307                      if (d1 > t0) {
00308                             t0 = d1;
00309                             gc[0].w = i<<1 | 1;
00310                      }
00311                      if (d0 < t1) {
00312                             t1 = d0;
00313                             gc[1].w = i<<1;
00314                      }
00315               }
00316        }
00317        if ((gc[0].w < 0) | (gc[1].w < 0))        /* paranoid check */
00318               return(FHUGE);
00319                                           /* compute intersections */
00320        VSUM(p[0], ro, rd, t0);
00321        VSUM(p[1], ro, rd, t1);
00322                                    /* now, compute grid coordinates */
00323        for (i = 0; i < 2; i++) {
00324               VSUB(vt, p[i], hp->orig);
00325               v = hp->wg[hdwg0[gc[i].w]];
00326               d = DOT(vt, v);
00327               if (d < 0 || d >= hp->grid[hdwg0[gc[i].w]])
00328                      return(FHUGE);              /* outside wall */
00329               gc[i].i[0] = d;
00330               if (r != NULL)
00331                      r[i][0] = 256. * (d - gc[i].i[0]);
00332               v = hp->wg[hdwg1[gc[i].w]];
00333               d = DOT(vt, v);
00334               if (d < 0 || d >= hp->grid[hdwg1[gc[i].w]])
00335                      return(FHUGE);              /* outside wall */
00336               gc[i].i[1] = d;
00337               if (r != NULL)
00338                      r[i][1] = 256. * (d - gc[i].i[1]);
00339        }
00340        if (ed != NULL)                    /* assign distance to exit point */
00341               *ed = t1;
00342        return(t0);                 /* return distance to entry point */
00343 }