Back to index

radiance  4R0+20100331
rholo2.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: rholo2.c,v 3.28 2004/01/01 11:21:55 schorsch Exp $";
00003 #endif
00004 /*
00005  * Rtrace support routines for holodeck rendering
00006  */
00007 
00008 #include <time.h>
00009 
00010 #include "rholo.h"
00011 #include "paths.h"
00012 #include "random.h"
00013 
00014 
00015 VIEWPOINT     myeye;        /* target view position */
00016 
00017 struct gclim {
00018        HOLO   *hp;                 /* holodeck pointer */
00019        GCOORD gc;                  /* grid cell */
00020        FVECT  egp;                 /* eye grid point */
00021        double erg2;                /* mean square eye grid range */
00022        double gmin[2], gmax[2];    /* grid coordinate limits */
00023 };                          /* a grid coordinate range */
00024 
00025 static void initeyelim(struct gclim       *gcl, HOLO    *hp, GCOORD   *gc);
00026 static void groweyelim(struct gclim *gcl, GCOORD *gc,
00027               double r0, double r1, int tight);
00028 static int clipeyelim(short rrng[2][2], struct gclim    *gcl);
00029 
00030 
00031 static void
00032 initeyelim(          /* initialize grid coordinate limits */
00033        register struct gclim       *gcl,
00034        register HOLO *hp,
00035        GCOORD *gc
00036 )
00037 {
00038        register RREAL       *v;
00039        register int  i;
00040 
00041        if (hp != NULL) {
00042               hdgrid(gcl->egp, gcl->hp = hp, myeye.vpt);
00043               gcl->erg2 = 0;
00044               for (i = 0, v = hp->wg[0]; i < 3; i++, v += 3)
00045                      gcl->erg2 += DOT(v,v);
00046               gcl->erg2 *= (1./3.) * myeye.rng*myeye.rng;
00047        }
00048        if (gc != NULL)
00049               gcl->gc = *gc;
00050        gcl->gmin[0] = gcl->gmin[1] = FHUGE;
00051        gcl->gmax[0] = gcl->gmax[1] = -FHUGE;
00052 }
00053 
00054 
00055 static void
00056 groweyelim(   /* grow grid limits about eye point */
00057        register struct gclim       *gcl,
00058        GCOORD *gc,
00059        double r0,
00060        double r1,
00061        int    tight
00062 )
00063 {
00064        FVECT  gp, ab;
00065        double ab2, od, cfact;
00066        double sqcoef[3], ctcoef[3], licoef[3], cnst;
00067        int    gw, gi[2];
00068        double wallpos, a, b, c, d, e, f;
00069        double root[2], yex;
00070        int    n, i, j, nex;
00071                                           /* point/view cone */
00072        i = gc->w>>1;
00073        gp[i] = gc->w&1 ? gcl->hp->grid[i] : 0;
00074        gp[hdwg0[gc->w]] = gc->i[0] + r0;
00075        gp[hdwg1[gc->w]] = gc->i[1] + r1;
00076        VSUB(ab, gcl->egp, gp);
00077        ab2 = DOT(ab, ab);
00078        gw = gcl->gc.w>>1;
00079        if ((i==gw ? ab[gw]*ab[gw] : ab2)  <= gcl->erg2 + FTINY) {
00080               gcl->gmin[0] = gcl->gmin[1] = -FHUGE;
00081               gcl->gmax[0] = gcl->gmax[1] = FHUGE;
00082               return;                     /* too close (to wall) */
00083        }
00084        ab2 = 1./ab2;                      /* 1/norm2(ab) */
00085        od = DOT(gp, ab);                  /* origin dot direction */
00086        cfact = 1./(1. - ab2*gcl->erg2);   /* tan^2 + 1 of cone angle */
00087        for (i = 0; i < 3; i++) {          /* compute cone equation */
00088               sqcoef[i] = ab[i]*ab[i]*cfact*ab2 - 1.;
00089               ctcoef[i] = 2.*ab[i]*ab[(i+1)%3]*cfact*ab2;
00090               licoef[i] = 2.*(gp[i] - ab[i]*cfact*od*ab2);
00091        }
00092        cnst = cfact*od*od*ab2 - DOT(gp,gp);
00093        /*
00094         * CONE:      sqcoef[0]*x*x + sqcoef[1]*y*y + sqcoef[2]*z*z
00095         *            + ctcoef[0]*x*y + ctcoef[1]*y*z + ctcoef[2]*z*x
00096         *            + licoef[0]*x + licoef[1]*y + licoef[2]*z + cnst == 0
00097         */
00098                             /* equation for conic section in plane */
00099        gi[0] = hdwg0[gcl->gc.w];
00100        gi[1] = hdwg1[gcl->gc.w];
00101        wallpos = gcl->gc.w&1 ? gcl->hp->grid[gw] : 0;
00102        a = sqcoef[gi[0]];                               /* x2 */
00103        b = ctcoef[gi[0]];                               /* xy */
00104        c = sqcoef[gi[1]];                               /* y2 */
00105        d = ctcoef[gw]*wallpos + licoef[gi[0]];                 /* x */
00106        e = ctcoef[gi[1]]*wallpos + licoef[gi[1]];              /* y */
00107        f = wallpos*(wallpos*sqcoef[gw] + licoef[gw]) + cnst;
00108        for (i = 0; i < 2; i++) {
00109               if (i) {             /* swap x and y coefficients */
00110                      register double      t;
00111                      t = a; a = c; c = t;
00112                      t = d; d = e; e = t;
00113               }
00114               nex = 0;             /* check global extrema */
00115               n = quadratic(root, a*(4.*a*c-b*b), 2.*a*(2.*c*d-b*e),
00116                             d*(c*d-b*e) + f*b*b);
00117               while (n-- > 0) {
00118                      if (gc->w>>1 == gi[i] &&
00119                                    (gc->w&1) ^ (root[n] < gp[gc->w>>1])) {
00120                             if (gc->w&1)
00121                                    gcl->gmin[i] = -FHUGE;
00122                             else
00123                                    gcl->gmax[i] = FHUGE;
00124                             nex++;
00125                             continue;            /* hyperbolic */
00126                      }
00127                      if (tight) {
00128                             yex = (-2.*a*root[n] - d)/b;
00129                             if (yex < gcl->gc.i[1-i] ||
00130                                           yex > gcl->gc.i[1-i]+1)
00131                                    continue;     /* outside cell */
00132                      }
00133                      if (root[n] < gcl->gmin[i])
00134                             gcl->gmin[i] = root[n];
00135                      if (root[n] > gcl->gmax[i])
00136                             gcl->gmax[i] = root[n];
00137                      nex++;
00138               }
00139                                    /* check local extrema */
00140               for (j = nex < 2 ? 2 : 0; j--; ) {
00141                      yex = gcl->gc.i[1-i] + j;
00142                      n = quadratic(root, a, b*yex+d, yex*(yex*c+e)+f);
00143                      while (n-- > 0) {
00144                             if (gc->w>>1 == gi[i] &&
00145                                    (gc->w&1) ^ (root[n] < gp[gc->w>>1]))
00146                                    continue;
00147                             if (root[n] < gcl->gmin[i])
00148                                    gcl->gmin[i] = root[n];
00149                             if (root[n] > gcl->gmax[i])
00150                                    gcl->gmax[i] = root[n];
00151                      }
00152               }
00153        }
00154 }
00155 
00156 
00157 static int
00158 clipeyelim(          /* clip eye limits to grid cell */
00159        register short       rrng[2][2],
00160        register struct gclim       *gcl
00161 )
00162 {
00163        int    incell = 1;
00164        register int  i;
00165 
00166        for (i = 0; i < 2; i++) {
00167               if (gcl->gmin[i] < gcl->gc.i[i])
00168                      gcl->gmin[i] = gcl->gc.i[i];
00169               if (gcl->gmax[i] > gcl->gc.i[i]+1)
00170                      gcl->gmax[i] = gcl->gc.i[i]+1;
00171               if (gcl->gmax[i] > gcl->gmin[i]) {
00172                      rrng[i][0] = 256.*(gcl->gmin[i] - gcl->gc.i[i]) +
00173                                    (1.-FTINY);
00174                      rrng[i][1] = 256.*(gcl->gmax[i] - gcl->gc.i[i]) +
00175                                    (1.-FTINY) - rrng[i][0];
00176               } else
00177                      rrng[i][0] = rrng[i][1] = 0;
00178               incell &= rrng[i][1] > 0;
00179        }
00180        return(incell);
00181 }
00182 
00183 
00184 extern void
00185 packrays(            /* pack ray origins and directions */
00186        register float       *rod,
00187        register PACKET      *p
00188 )
00189 {
00190 #if 0
00191        double dist2sum = 0.;
00192        FVECT  vt;
00193 #endif
00194        int    nretries = p->nr + 2;
00195        struct gclim  eyelim;
00196        short  rrng0[2][2], rrng1[2][2];
00197        int    useyelim;
00198        GCOORD gc[2];
00199        FVECT  ro, rd;
00200        double d;
00201        register int  i;
00202 
00203        if (!hdbcoord(gc, hdlist[p->hd], p->bi))
00204               error(CONSISTENCY, "bad beam index in packrays");
00205        if ((useyelim = myeye.rng > FTINY)) {
00206               initeyelim(&eyelim, hdlist[p->hd], gc);
00207               groweyelim(&eyelim, gc+1, 0., 0., 0);
00208               groweyelim(&eyelim, gc+1, 1., 1., 0);
00209               useyelim = clipeyelim(rrng0, &eyelim);
00210 #ifdef DEBUG
00211               if (!useyelim)
00212                      error(WARNING, "no eye overlap in packrays");
00213 #endif
00214        }
00215        for (i = 0; i < p->nr; i++) {
00216        retry:
00217               if (useyelim) {
00218                      initeyelim(&eyelim, NULL, gc+1);
00219                      p->ra[i].r[0][0] = (int)(frandom()*rrng0[0][1])
00220                                           + rrng0[0][0];
00221                      p->ra[i].r[0][1] = (int)(frandom()*rrng0[1][1])
00222                                           + rrng0[1][0];
00223                      groweyelim(&eyelim, gc,
00224                                    (1./256.)*(p->ra[i].r[0][0]+.5),
00225                                    (1./256.)*(p->ra[i].r[0][1]+.5), 1);
00226                      if (!clipeyelim(rrng1, &eyelim)) {
00227                             useyelim = nretries-- > 0;
00228 #ifdef DEBUG
00229                             if (!useyelim)
00230                                    error(WARNING,
00231                                    "exceeded retry limit in packrays");
00232 #endif
00233                             goto retry;
00234                      }
00235                      p->ra[i].r[1][0] = (int)(frandom()*rrng1[0][1])
00236                                           + rrng1[0][0];
00237                      p->ra[i].r[1][1] = (int)(frandom()*rrng1[1][1])
00238                                           + rrng1[1][0];
00239               } else {
00240                      p->ra[i].r[0][0] = frandom() * 256.;
00241                      p->ra[i].r[0][1] = frandom() * 256.;
00242                      p->ra[i].r[1][0] = frandom() * 256.;
00243                      p->ra[i].r[1][1] = frandom() * 256.;
00244               }
00245               d = hdray(ro, rd, hdlist[p->hd], gc, p->ra[i].r);
00246 #if 0
00247               VSUM(vt, ro, rd, d);
00248               dist2sum += dist2line(myeye.vpt, ro, vt);
00249 #endif
00250               if (p->offset != NULL) {
00251                      if (!vdef(OBSTRUCTIONS))
00252                             d *= frandom();             /* random offset */
00253                      VSUM(ro, ro, rd, d);        /* advance ray */
00254                      p->offset[i] = d;
00255               }
00256               VCOPY(rod, ro);
00257               rod += 3;
00258               VCOPY(rod, rd);
00259               rod += 3;
00260        }
00261 #if 0
00262        fprintf(stderr, "%f RMS (%d retries)\t", sqrt(dist2sum/p->nr),
00263                      p->nr + 2 - nretries);
00264 #endif
00265 }
00266 
00267 
00268 extern void
00269 donerays(            /* encode finished ray computations */
00270        register PACKET      *p,
00271        register float       *rvl
00272 )
00273 {
00274        double d;
00275        register int  i;
00276 
00277        for (i = 0; i < p->nr; i++) {
00278               setcolr(p->ra[i].v, rvl[0], rvl[1], rvl[2]);
00279               d = rvl[3];
00280               if (p->offset != NULL)
00281                      d += p->offset[i];
00282               p->ra[i].d = hdcode(hdlist[p->hd], d);
00283               rvl += 4;
00284        }
00285        p->nc += p->nr;
00286 }
00287 
00288 
00289 extern int
00290 done_rtrace(void)                  /* clean up and close rtrace calculation */
00291 {
00292        int    status;
00293                                    /* already closed? */
00294        if (!nprocs)
00295               return(0);
00296                                    /* flush beam queue */
00297        done_packets(flush_queue());
00298                                    /* sync holodeck */
00299        hdsync(NULL, 1);
00300                                    /* close rtrace */
00301        if ((status = end_rtrace()))
00302               error(WARNING, "bad exit status from rtrace");
00303        if (vdef(REPORT)) {         /* report time */
00304               eputs("rtrace process closed\n");
00305               report(0);
00306        }
00307        return(status);                    /* return status */
00308 }
00309 
00310 
00311 extern void
00312 new_rtrace(void)                   /* restart rtrace calculation */
00313 {
00314        char   combuf[128];
00315 
00316        if (nprocs > 0)                    /* already running? */
00317               return;
00318        starttime = time(NULL);            /* reset start time and counts */
00319        npacksdone = nraysdone = 0L;
00320        if (vdef(TIME))                    /* reset end time */
00321               endtime = starttime + vflt(TIME)*3600. + .5;
00322        if (vdef(RIF)) {            /* rerun rad to update octree */
00323               sprintf(combuf, "rad -v 0 -s -w %s", vval(RIF));
00324               if (system(combuf))
00325                      error(WARNING, "error running rad");
00326        }
00327        if (start_rtrace() < 1)            /* start rtrace */
00328               error(WARNING, "cannot restart rtrace");
00329        else if (vdef(REPORT)) {
00330               eputs("rtrace process restarted\n");
00331               report(0);
00332        }
00333 }
00334 
00335 
00336 extern int
00337 getradfile(void)                   /* run rad and get needed variables */
00338 {
00339        static short  mvar[] = {OCTREE,EYESEP,-1};
00340        static char   tf1[] = TEMPLATE;
00341        char   tf2[64];
00342        char   combuf[256];
00343        char   *pippt = NULL;
00344        register int  i;
00345        register char *cp;
00346                                    /* check if rad file specified */
00347        if (!vdef(RIF))
00348               return(0);
00349                                    /* create rad command */
00350        mktemp(tf1);
00351        sprintf(tf2, "%s.rif", tf1);
00352        sprintf(combuf,
00353               "rad -v 0 -s -e -w %s OPTFILE=%s | egrep '^[ \t]*(NOMATCH",
00354                      vval(RIF), tf1);
00355        cp = combuf;
00356        while (*cp){
00357               if (*cp == '|') pippt = cp;
00358               cp++;
00359        }                           /* match unset variables */
00360        for (i = 0; mvar[i] >= 0; i++)
00361               if (!vdef(mvar[i])) {
00362                      *cp++ = '|';
00363                      strcpy(cp, vnam(mvar[i]));
00364                      while (*cp) cp++;
00365                      pippt = NULL;
00366               }
00367        if (pippt != NULL)
00368               strcpy(pippt, "> " NULL_DEVICE);   /* nothing to match */
00369        else
00370               sprintf(cp, ")[ \t]*=' > %s", tf2);
00371 #ifdef DEBUG
00372        wputs(combuf); wputs("\n");
00373 #endif
00374        system(combuf);                           /* ignore exit code */
00375        if (pippt == NULL) {
00376               loadvars(tf2);                     /* load variables */
00377               unlink(tf2);
00378        }
00379        rtargc += wordfile(rtargv+rtargc, tf1);   /* get rtrace options */
00380        unlink(tf1);                /* clean up */
00381        return(1);
00382 }
00383 
00384 
00385 extern void
00386 report(                     /* report progress so far */
00387        time_t t
00388 )
00389 {
00390        static time_t seconds2go = 1000000;
00391 
00392        if (t == 0L)
00393               t = time(NULL);
00394        sprintf(errmsg, "%ld packets (%ld rays) done after %.2f hours\n",
00395                      npacksdone, nraysdone, (t-starttime)/3600.);
00396        eputs(errmsg);
00397        if (seconds2go == 1000000)
00398               seconds2go = vdef(REPORT) ? (long)(vflt(REPORT)*60. + .5) : 0L;
00399        if (seconds2go)
00400               reporttime = t + seconds2go;
00401 }