Back to index

radiance  4R0+20100331
rholo3.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: rholo3.c,v 3.42 2004/06/08 19:48:30 greg Exp $";
00003 #endif
00004 /*
00005  * Routines for tracking beam compuatations
00006  */
00007 
00008 #include "rholo.h"
00009 
00010 #ifndef NFRAG2CHUNK
00011 #define NFRAG2CHUNK  4096   /* number of fragments to start chunking */
00012 #endif
00013 
00014 #ifndef abs
00015 #define abs(x)              ((x) > 0 ? (x) : -(x))
00016 #endif
00017 #ifndef sgn
00018 #define sgn(x)              ((x) > 0 ? 1 : (x) < 0 ? -1 : 0)
00019 #endif
00020 
00021 #define rchunk(n)    (((n)+(RPACKSIZ/2))/RPACKSIZ)
00022 
00023 int    chunkycmp = 0;              /* clump beams together on disk */
00024 
00025 static PACKHEAD      *complist=NULL;      /* list of beams to compute */
00026 static int    complen=0;    /* length of complist */
00027 static int    listpos=0;    /* current list position for next_packet */
00028 static int    lastin= -1;   /* last ordered position in list */
00029 
00030 static void sortcomplist(void);
00031 static void mergeclists(PACKHEAD *cdest, PACKHEAD *cl1, int n1, PACKHEAD *cl2, int n2);
00032 static void view_list(FILE  *fp);
00033 static void ambient_list(void);
00034 static double beamvolume(HOLO      *hp, int      bi);
00035 static void dispbeam(BEAM   *b, HDBEAMI   *hb);
00036 
00037 
00038 
00039 static int
00040 beamcmp(b0, b1)                           /* comparison for compute order */
00041 register PACKHEAD    *b0, *b1;
00042 {
00043        BEAMI  *bip0, *bip1;
00044        register long c;
00045                                    /* first check desired quantities */
00046        if (chunkycmp)
00047               c = rchunk(b1->nr)*(rchunk(b0->nc)+1L) -
00048                             rchunk(b0->nr)*(rchunk(b1->nc)+1L);
00049        else
00050               c = b1->nr*(b0->nc+1L) - b0->nr*(b1->nc+1L);
00051        if (c > 0) return(1);
00052        if (c < 0) return(-1);
00053                             /* only one file, so skip the following: */
00054 #if 0
00055                                    /* next, check file descriptors */
00056        c = hdlist[b0->hd]->fd - hdlist[b1->hd]->fd;
00057        if (c) return(c);
00058 #endif
00059                                    /* finally, check file positions */
00060        bip0 = &hdlist[b0->hd]->bi[b0->bi];
00061        bip1 = &hdlist[b1->hd]->bi[b1->bi];
00062                                    /* put diskless beams last */
00063        if (!bip0->nrd)
00064               return(bip1->nrd > 0);
00065        if (!bip1->nrd)
00066               return(-1);
00067        c = bip0->fo - bip1->fo;
00068        return(c < 0 ? -1 : c > 0);
00069 }
00070 
00071 
00072 int
00073 beamidcmp(b0, b1)                  /* comparison for beam searching */
00074 register PACKHEAD    *b0, *b1;
00075 {
00076        register int  c = b0->hd - b1->hd;
00077 
00078        if (c) return(c);
00079        return(b0->bi - b1->bi);
00080 }
00081 
00082 
00083 static void
00084 dispbeam(                          /* display a holodeck beam */
00085        register BEAM *b,
00086        register HDBEAMI     *hb
00087 )
00088 {
00089        static int    n = 0;
00090        static PACKHEAD      *p = NULL;
00091 
00092        if (b == NULL)
00093               return;
00094        if (b->nrm > n) {           /* (re)allocate packet holder */
00095               n = b->nrm;
00096               if (p == NULL) p = (PACKHEAD *)malloc(packsiz(n));
00097               else p = (PACKHEAD *)realloc((void *)p, packsiz(n));
00098               CHECK(p==NULL, SYSTEM, "out of memory in dispbeam");
00099        }
00100                                    /* assign packet fields */
00101        memcpy((void *)packra(p), (void *)hdbray(b), b->nrm*sizeof(RAYVAL));
00102        p->nr = p->nc = b->nrm;
00103        for (p->hd = 0; hdlist[p->hd] != hb->h; p->hd++)
00104               if (hdlist[p->hd] == NULL)
00105                      error(CONSISTENCY, "unregistered holodeck in dispbeam");
00106        p->bi = hb->b;
00107        disp_packet(p);                    /* display it */
00108        if (n >= 1024) {            /* free ridiculous packets */
00109               free((void *)p);
00110               p = NULL; n = 0;
00111        }
00112 }
00113 
00114 
00115 extern void
00116 bundle_set(   /* bundle set operation */
00117        int    op,
00118        PACKHEAD      *clist,
00119        int    nents
00120 )
00121 {
00122        int    oldnr, n;
00123        HDBEAMI       *hbarr;
00124        register PACKHEAD    *csm;
00125        register int  i;
00126                                    /* search for common members */
00127        for (csm = clist+nents; csm-- > clist; )
00128               csm->nc = -1;
00129        qsort((void *)clist, nents, sizeof(PACKHEAD), beamidcmp);
00130        for (i = 0; i < complen; i++) {
00131               csm = (PACKHEAD *)bsearch((void *)(complist+i), (void *)clist,
00132                             nents, sizeof(PACKHEAD), beamidcmp);
00133               if (csm == NULL)
00134                      continue;
00135               oldnr = complist[i].nr;
00136               csm->nc = complist[i].nc;
00137               switch (op) {
00138               case BS_ADD:         /* add to count */
00139                      complist[i].nr += csm->nr;
00140                      csm->nr = 0;
00141                      break;
00142               case BS_MAX:         /* maximum of counts */
00143                      if (csm->nr > complist[i].nr)
00144                             complist[i].nr = csm->nr;
00145                      csm->nr = 0;
00146                      break;
00147               case BS_ADJ:         /* reset count */
00148                      complist[i].nr = csm->nr;
00149                      csm->nr = 0;
00150                      break;
00151               case BS_DEL:         /* delete count */
00152                      if (csm->nr == 0 || csm->nr >= complist[i].nr)
00153                             complist[i].nr = 0;
00154                      else
00155                             complist[i].nr -= csm->nr;
00156                      break;
00157               }
00158               if (complist[i].nr != oldnr)
00159                      lastin = -1;  /* flag sort */
00160        }
00161                             /* record computed rays for uncommon beams */
00162        for (csm = clist+nents; csm-- > clist; )
00163               if (csm->nc < 0)
00164                      csm->nc = bnrays(hdlist[csm->hd], csm->bi);
00165                             /* complete list operations */
00166        switch (op) {
00167        case BS_NEW:                /* new computation set */
00168               listpos = 0; lastin = -1;
00169               if (complen)         /* free old list */
00170                      free((void *)complist);
00171               complist = NULL;
00172               if (!(complen = nents))
00173                      return;
00174               complist = (PACKHEAD *)malloc(nents*sizeof(PACKHEAD));
00175               if (complist == NULL)
00176                      goto memerr;
00177               memcpy((void *)complist, (void *)clist, nents*sizeof(PACKHEAD));
00178               break;
00179        case BS_ADD:                /* add to computation set */
00180        case BS_MAX:                /* maximum of quantities */
00181        case BS_ADJ:                /* adjust set quantities */
00182               if (nents <= 0)
00183                      return;
00184               sortcomplist();             /* sort updated list & new entries */
00185               qsort((void *)clist, nents, sizeof(PACKHEAD), beamcmp);
00186                                    /* what can't we satisfy? */
00187               for (i = nents, csm = clist; i-- && csm->nr > csm->nc; csm++)
00188                      ;
00189               n = csm - clist;
00190               if (op != BS_ADD) {  /* don't regenerate adjusted beams */
00191                      for (++i; i-- && csm->nr > 0; csm++)
00192                             ;
00193                      nents = csm - clist;
00194               }
00195               if (n) {             /* allocate space for merged list */
00196                      PACKHEAD      *newlist;
00197                      newlist = (PACKHEAD *)malloc(
00198                                    (complen+n)*sizeof(PACKHEAD) );
00199                      if (newlist == NULL)
00200                             goto memerr;
00201                                           /* merge lists */
00202                      mergeclists(newlist, clist, n, complist, complen);
00203                      if (complen)
00204                             free((void *)complist);
00205                      complist = newlist;
00206                      complen += n;
00207               }
00208               listpos = 0;
00209               lastin = complen-1;  /* list is now sorted */
00210               break;
00211        case BS_DEL:                /* delete from computation set */
00212               return;                     /* already done */
00213        default:
00214               error(CONSISTENCY, "bundle_set called with unknown operation");
00215        }
00216        if (outdev == NULL || !nents)      /* nothing to display? */
00217               return;
00218                                    /* load and display beams we have */
00219        hbarr = (HDBEAMI *)malloc(nents*sizeof(HDBEAMI));
00220        for (i = nents; i--; ) {
00221               hbarr[i].h = hdlist[clist[i].hd];
00222               hbarr[i].b = clist[i].bi;
00223        }
00224        hdloadbeams(hbarr, nents, dispbeam);
00225        free((void *)hbarr);
00226        if (hdfragflags&FF_READ) {
00227               listpos = 0;
00228               lastin = -1;         /* need to re-sort list */
00229        }
00230        return;
00231 memerr:
00232        error(SYSTEM, "out of memory in bundle_set");
00233 }
00234 
00235 
00236 static double
00237 beamvolume(   /* compute approximate volume of a beam */
00238        HOLO   *hp,
00239        int    bi
00240 )
00241 {
00242        GCOORD gc[2];
00243        FVECT  cp[4], edgeA, edgeB, cent[2];
00244        FVECT  crossp[2], diffv;
00245        double vol[2];
00246        register int  i;
00247                                    /* get grid coordinates */
00248        if (!hdbcoord(gc, hp, bi))
00249               error(CONSISTENCY, "bad beam index in beamvolume");
00250        for (i = 0; i < 2; i++) {   /* compute cell area vectors */
00251               hdcell(cp, hp, gc+i);
00252               VSUM(edgeA, cp[1], cp[0], -1.0);
00253               VSUM(edgeB, cp[3], cp[1], -1.0);
00254               fcross(crossp[i], edgeA, edgeB);
00255                                    /* compute center */
00256               cent[i][0] = 0.5*(cp[0][0] + cp[2][0]);
00257               cent[i][1] = 0.5*(cp[0][1] + cp[2][1]);
00258               cent[i][2] = 0.5*(cp[0][2] + cp[2][2]);
00259        }
00260                                    /* compute difference vector */
00261        VSUM(diffv, cent[1], cent[0], -1.0);
00262        for (i = 0; i < 2; i++) {   /* compute volume contributions */
00263               vol[i] = 0.5*DOT(crossp[i], diffv);
00264               if (vol[i] < 0.) vol[i] = -vol[i];
00265        }
00266        return(vol[0] + vol[1]);    /* return total volume */
00267 }
00268 
00269 
00270 static void
00271 ambient_list(void)                 /* compute ambient beam list */
00272 {
00273        int32  wtotal, minrt;
00274        double frac;
00275        int    i;
00276        register int  j, k;
00277 
00278        complen = 0;
00279        for (j = 0; hdlist[j] != NULL; j++)
00280               complen += nbeams(hdlist[j]);
00281        complist = (PACKHEAD *)malloc(complen*sizeof(PACKHEAD));
00282        CHECK(complist==NULL, SYSTEM, "out of memory in ambient_list");
00283                                    /* compute beam weights */
00284        k = 0; wtotal = 0;
00285        for (j = 0; hdlist[j] != NULL; j++) {
00286               frac = 512. * VLEN(hdlist[j]->wg[0]) *
00287                             VLEN(hdlist[j]->wg[1]) *
00288                             VLEN(hdlist[j]->wg[2]);
00289               for (i = nbeams(hdlist[j]); i > 0; i--) {
00290                      complist[k].hd = j;
00291                      complist[k].bi = i;
00292                      complist[k].nr = frac*beamvolume(hdlist[j], i) + 0.5;
00293                      complist[k].nc = bnrays(hdlist[j], i);
00294                      wtotal += complist[k++].nr;
00295               }
00296        }
00297                                    /* adjust sample weights */
00298        if (vdef(DISKSPACE))
00299               frac = 1024.*1024.*vflt(DISKSPACE) / (wtotal*sizeof(RAYVAL));
00300        else
00301               frac = 1024.*1024.*2048. / (wtotal*sizeof(RAYVAL));
00302        minrt = .02*frac*wtotal/complen + .5;     /* heuristic mimimum */
00303        if (minrt > RPACKSIZ)
00304               minrt = RPACKSIZ;
00305        for (k = complen; k--; )
00306               if ((complist[k].nr = frac*complist[k].nr + 0.5) < minrt)
00307                      complist[k].nr = minrt;
00308        listpos = 0; lastin = -1;   /* flag initial sort */
00309 }
00310 
00311 
00312 static void
00313 view_list(                  /* assign beam priority from view list */
00314        FILE   *fp
00315 )
00316 {
00317        double pa = 1.;
00318        VIEW   curview;
00319        int    xr, yr;
00320        char   *err;
00321        BEAMLIST      blist;
00322 
00323        curview = stdview;
00324        while (nextview(&curview, fp) != EOF) {
00325               if ((err = setview(&curview)) != NULL) {
00326                      error(WARNING, err);
00327                      continue;
00328               }
00329               xr = yr = 1024;
00330               normaspect(viewaspect(&curview), &pa, &xr, &yr);
00331               viewbeams(&curview, xr, yr, &blist);
00332               bundle_set(BS_MAX, blist.bl, blist.nb);
00333               free((void *)blist.bl);
00334        }
00335 }
00336 
00337 
00338 extern void
00339 init_global(void)                  /* initialize global ray computation */
00340 {
00341                                    /* free old list and empty queue */
00342        if (complen > 0) {
00343               free((void *)complist);
00344               done_packets(flush_queue());
00345        }
00346                                    /* reseed random number generator */
00347        srandom(time(NULL));
00348                                    /* allocate beam list */
00349        if (readinp)
00350               view_list(stdin);
00351        else
00352               ambient_list();
00353                                    /* no view vicinity */
00354        myeye.rng = 0;
00355 }
00356 
00357 
00358 static void
00359 mergeclists(  /* merge two sorted lists */
00360        register PACKHEAD    *cdest,
00361        register PACKHEAD    *cl1,
00362        int    n1,
00363        register PACKHEAD    *cl2,
00364        int    n2
00365 )
00366 {
00367        register int  cmp;
00368 
00369        while (n1 | n2) {
00370               if (!n1) cmp = 1;
00371               else if (!n2) cmp = -1;
00372               else cmp = beamcmp(cl1, cl2);
00373               if (cmp > 0) {
00374                      *cdest = *cl2;
00375                      cl2++; n2--;
00376               } else {
00377                      *cdest = *cl1;
00378                      cl1++; n1--;
00379               }
00380               cdest++;
00381        }
00382 }
00383 
00384 
00385 static void
00386 sortcomplist(void)                 /* fix our list order */
00387 {
00388        PACKHEAD      *list2;
00389        int    listlen;
00390        register int  i;
00391 
00392        if (complen <= 0)    /* check to see if there is even a list */
00393               return;
00394        if (!chunkycmp)             /* check to see if fragment list is full */
00395               if (!hdfragOK(hdlist[0]->fd, &listlen, NULL)
00396 #if NFRAG2CHUNK
00397                             || listlen >= NFRAG2CHUNK
00398 #endif
00399                             ) {
00400                      chunkycmp++;  /* use "chunky" comparison */
00401                      lastin = -1;  /* need to re-sort list */
00402 #ifdef DEBUG
00403                      error(WARNING, "using chunky comparison mode");
00404 #endif
00405               }
00406        if (lastin < 0 || listpos*4 >= complen*3)
00407               qsort((void *)complist, complen, sizeof(PACKHEAD), beamcmp);
00408        else if (listpos) {  /* else sort and merge sublist */
00409               list2 = (PACKHEAD *)malloc(listpos*sizeof(PACKHEAD));
00410               CHECK(list2==NULL, SYSTEM, "out of memory in sortcomplist");
00411               memcpy((void *)list2,(void *)complist,listpos*sizeof(PACKHEAD));
00412               qsort((void *)list2, listpos, sizeof(PACKHEAD), beamcmp);
00413               mergeclists(complist, list2, listpos,
00414                             complist+listpos, complen-listpos);
00415               free((void *)list2);
00416        }
00417                                    /* drop satisfied requests */
00418        for (i = complen; i-- && complist[i].nr <= complist[i].nc; )
00419               ;
00420        if (i < 0) {
00421               free((void *)complist);
00422               complist = NULL;
00423               complen = 0;
00424        } else if (i < complen-1) {
00425               list2 = (PACKHEAD *)realloc((void *)complist,
00426                             (i+1)*sizeof(PACKHEAD));
00427               if (list2 != NULL)
00428                      complist = list2;
00429               complen = i+1;
00430        }
00431        listpos = 0; lastin = i;
00432 }
00433 
00434 
00435 /*
00436  * The following routine works on the assumption that the bundle weights are
00437  * more or less evenly distributed, such that computing a packet causes
00438  * a given bundle to move way down in the computation order.  We keep
00439  * track of where the computed bundle with the highest priority would end
00440  * up, and if we get further in our compute list than this, we re-sort the
00441  * list and start again from the beginning.  Since
00442  * a merge sort is used, the sorting costs are minimal.
00443  */
00444 extern int
00445 next_packet(         /* prepare packet for computation */
00446        register PACKET      *p,
00447        int    n
00448 )
00449 {
00450        if (listpos > lastin)              /* time to sort the list */
00451               sortcomplist();
00452        if (complen <= 0)
00453               return(0);
00454        p->hd = complist[listpos].hd;
00455        p->bi = complist[listpos].bi;
00456        p->nc = complist[listpos].nc;
00457        p->nr = complist[listpos].nr - p->nc;
00458        if (p->nr <= 0)
00459               return(0);
00460        DCHECK(n < 1 | n > RPACKSIZ,
00461                      CONSISTENCY, "next_packet called with bad n value");
00462        if (p->nr > n)
00463               p->nr = n;
00464        complist[listpos].nc += p->nr;     /* find where this one would go */
00465        if (hdgetbeam(hdlist[p->hd], p->bi) != NULL)
00466               hdfreefrag(hdlist[p->hd], p->bi);
00467        while (lastin > listpos && 
00468                      beamcmp(complist+lastin, complist+listpos) > 0)
00469               lastin--;
00470        listpos++;
00471        return(1);
00472 }