--- ray/src/hd/rholo3.c 1997/11/04 11:52:21 3.4 +++ ray/src/hd/rholo3.c 1997/12/18 14:46:00 3.14 @@ -10,11 +10,9 @@ static char SCCSid[] = "$SunId$ SGI"; #include "rholo.h" - #define abs(x) ((x) > 0 ? (x) : -(x)) #define sgn(x) ((x) > 0 ? 1 : (x) < 0 ? -1 : 0) - static PACKHEAD *complist=NULL; /* list of beams to compute */ static int complen=0; /* length of complist */ static int listpos=0; /* current list position for next_packet */ @@ -25,18 +23,44 @@ int beamcmp(b0, b1) /* comparison for descending compute order */ register PACKHEAD *b0, *b1; { - return( b1->nr*(bnrays(hdlist[b0->hd],b0->bi)+1) - - b0->nr*(bnrays(hdlist[b1->hd],b1->bi)+1) ); + return( b1->nr*(b0->nc+1) - b0->nr*(b1->nc+1) ); } +int +dispbeam(b, hp, bi) /* display a holodeck beam */ +register BEAM *b; +HOLO *hp; +int bi; +{ + static int n = 0; + static PACKHEAD *p = NULL; + + if (b == NULL) + return; + if (b->nrm > n) { /* (re)allocate packet holder */ + n = b->nrm; + if (p == NULL) p = (PACKHEAD *)malloc(packsiz(n)); + else p = (PACKHEAD *)realloc((char *)p, packsiz(n)); + if (p == NULL) + error(SYSTEM, "out of memory in dispbeam"); + } + /* assign packet fields */ + bcopy((char *)hdbray(b), (char *)packra(p), b->nrm*sizeof(RAYVAL)); + p->nr = p->nc = b->nrm; + for (p->hd = 0; hdlist[p->hd] != hp; p->hd++) + if (hdlist[p->hd] == NULL) + error(CONSISTENCY, "unregistered holodeck in dispbeam"); + p->bi = bi; + disp_packet(p); /* display it */ +} + + bundle_set(op, clist, nents) /* bundle set operation */ int op; -PACKHEAD *clist; +register PACKHEAD *clist; int nents; { - BEAM *b; - PACKHEAD *p; register int i, n; switch (op) { @@ -53,31 +77,41 @@ int nents; if (complist == NULL) goto memerr; bcopy((char *)clist, (char *)complist, nents*sizeof(PACKHEAD)); - complen = nents; - listpos = 0; - lastin = -1; /* flag for initial sort */ + complen = nents; /* finish initialization below */ break; case BS_ADD: /* add to computation set */ + case BS_ADJ: /* adjust set quantities */ if (nents <= 0) return; /* merge any common members */ - for (i = 0; i < complen; i++) - for (n = 0; n < nents; n++) + for (n = 0; n < nents; n++) { + for (i = 0; i < complen; i++) if (clist[n].bi == complist[i].bi && clist[n].hd == complist[i].hd) { - complist[i].nr += clist[n].nr; + int oldnr = complist[i].nr; + if (op == BS_ADD) + complist[i].nr += clist[n].nr; + else /* op == BS_ADJ */ + complist[i].nr = clist[n].nr; clist[n].nr = 0; - lastin = -1; /* flag full sort */ + clist[n].nc = complist[i].nc; + if (complist[i].nr != oldnr) + lastin = -1; /* flag sort */ break; } + if (i >= complen) + clist[n].nc = bnrays(hdlist[clist[n].hd], + clist[n].bi); + } /* sort updated list */ sortcomplist(); /* sort new entries */ qsort((char *)clist, nents, sizeof(PACKHEAD), beamcmp); /* what can't we satisfy? */ - for (n = 0; n < nents && clist[n].nr > - bnrays(hdlist[clist[n].hd],clist[n].bi); n++) + for (n = 0; n < nents && clist[n].nr > clist[n].nc; n++) ; + if (op == BS_ADJ) + nents = n; if (n) { /* allocate space for merged list */ PACKHEAD *newlist; newlist = (PACKHEAD *)malloc( @@ -110,101 +144,78 @@ int nents; lastin = -1; /* flag full sort */ break; } - if (lastin < 0) /* sort updated list */ - sortcomplist(); return; /* no display */ default: error(CONSISTENCY, "bundle_set called with unknown operation"); } - n = 0; /* allocate packet holder */ - for (i = 0; i < nents; i++) - if (clist[i].nr > n) - n = clist[i].nr; - p = (PACKHEAD *)malloc(sizeof(PACKHEAD) + n*sizeof(RAYVAL)); - if (p == NULL) - goto memerr; - /* display what we have */ - for (i = 0; i < nents; i++) - if ((b = hdgetbeam(hdlist[clist[i].hd], clist[i].bi)) != NULL) { - bcopy((char *)hdbray(b), (char *)(p+1), - (p->nr=b->nrm)*sizeof(RAYVAL)); - disp_packet((PACKET *)p); + if (outdev != NULL) { /* load and display beams we have */ + register HDBEAMI *hb; + + hb = (HDBEAMI *)malloc(nents*sizeof(HDBEAMI)); + for (i = 0; i < nents; i++) { + hb[i].h = hdlist[clist[i].hd]; + hb[i].b = clist[i].bi; } - free((char *)p); /* clean up */ + hdloadbeams(hb, nents, dispbeam); + free((char *)hb); + } + if (op == BS_NEW) { + done_packets(flush_queue()); /* empty queue, so we can... */ + for (i = 0; i < complen; i++) /* ...get number computed */ + complist[i].nc = bnrays(hdlist[complist[i].hd], + complist[i].bi); + listpos = 0; + lastin = -1; /* flag for initial sort */ + } return; memerr: error(SYSTEM, "out of memory in bundle_set"); } -int -weightf(hp, x0, x1, x2) /* voxel weighting function */ -register HOLO *hp; -register int x0, x1, x2; -{ - switch (vlet(OCCUPANCY)) { - case 'U': /* uniform weighting */ - return(1); - case 'C': /* center weighting (crude) */ - x0 += x0 - hp->grid[0] + 1; - x0 = abs(x0)*hp->grid[1]*hp->grid[2]; - x1 += x1 - hp->grid[1] + 1; - x1 = abs(x1)*hp->grid[0]*hp->grid[2]; - x2 += x2 - hp->grid[2] + 1; - x2 = abs(x2)*hp->grid[0]*hp->grid[1]; - return(hp->grid[0]*hp->grid[1]*hp->grid[2] - - (x0+x1+x2)/3); - default: - badvalue(OCCUPANCY); - } -} - - -/* The following is by Daniel Cohen, taken from Graphics Gems IV, p. 368 */ -long -lineweight(hp, x, y, z, dx, dy, dz) /* compute weights along a line */ +double +beamvolume(hp, bi) /* compute approximate volume of a beam */ HOLO *hp; -int x, y, z, dx, dy, dz; +int bi; { - long wres = 0; - int n, sx, sy, sz, exy, exz, ezy, ax, ay, az, bx, by, bz; - - sx = sgn(dx); sy = sgn(dy); sz = sgn(dz); - ax = abs(dx); ay = abs(dy); az = abs(dz); - bx = 2*ax; by = 2*ay; bz = 2*az; - exy = ay-ax; exz = az-ax; ezy = ay-az; - n = ax+ay+az + 1; /* added increment to visit last */ - while (n--) { - wres += weightf(hp, x, y, z); - if (exy < 0) { - if (exz < 0) { - x += sx; - exy += by; exz += bz; - } else { - z += sz; - exz -= bx; ezy += by; - } - } else { - if (ezy < 0) { - z += sz; - exz -= bx; ezy += by; - } else { - y += sy; - exy -= bx; ezy -= bz; - } - } + GCOORD gc[2]; + FVECT cp[4], edgeA, edgeB, cent[2]; + FVECT v, crossp[2], diffv; + double vol[2]; + register int i; + /* get grid coordinates */ + if (!hdbcoord(gc, hp, bi)) + error(CONSISTENCY, "bad beam index in beamvolume"); + for (i = 0; i < 2; i++) { /* compute cell area vectors */ + hdcell(cp, hp, gc+i); + VSUM(edgeA, cp[1], cp[0], -1.0); + VSUM(edgeB, cp[3], cp[1], -1.0); + fcross(crossp[i], edgeA, edgeB); + VSUM(edgeA, cp[2], cp[3], -1.0); + VSUM(edgeB, cp[0], cp[2], -1.0); + fcross(v, edgeA, edgeB); + VSUM(crossp[i], crossp[i], v, 1.0); + /* compute center */ + cent[i][0] = 0.5*(cp[0][0] + cp[2][0]); + cent[i][1] = 0.5*(cp[0][1] + cp[2][1]); + cent[i][2] = 0.5*(cp[0][2] + cp[2][2]); } - return(wres); + /* compute difference vector */ + VSUM(diffv, cent[1], cent[0], -1.0); + for (i = 0; i < 2; i++) { /* compute volume contributions */ + vol[i] = 0.25*DOT(crossp[i], diffv); + if (vol[i] < 0.) vol[i] = -vol[i]; + } + return(vol[0] + vol[1]); /* return total volume */ } init_global() /* initialize global ray computation */ { long wtotal = 0; - int i, j; - int lseg[2][3]; double frac; - register int k; + int i; + register int j, k; /* free old list */ if (complen > 0) free((char *)complist); @@ -217,25 +228,23 @@ init_global() /* initialize global ray computation * error(SYSTEM, "out of memory in init_global"); /* compute beam weights */ k = 0; - for (j = 0; hdlist[j] != NULL; j++) + for (j = 0; hdlist[j] != NULL; j++) { + frac = 512. * hdlist[j]->wg[0] * + hdlist[j]->wg[1] * hdlist[j]->wg[2]; for (i = nbeams(hdlist[j]); i > 0; i--) { - hdlseg(lseg, hdlist[j], i); complist[k].hd = j; complist[k].bi = i; - complist[k].nr = lineweight( hdlist[j], - lseg[0][0], lseg[0][1], lseg[0][2], - lseg[1][0] - lseg[0][0], - lseg[1][1] - lseg[0][1], - lseg[1][2] - lseg[0][2] ); + complist[k].nr = frac*beamvolume(hdlist[j], i) + 0.5; wtotal += complist[k++].nr; } + } /* adjust weights */ - if (vdef(DISKSPACE)) { + if (vdef(DISKSPACE)) frac = 1024.*1024.*vflt(DISKSPACE) / (wtotal*sizeof(RAYVAL)); - if (frac < 0.95 | frac > 1.05) - while (k--) - complist[k].nr = frac * complist[k].nr; - } + else + frac = 1024.*1024.*16384. / (wtotal*sizeof(RAYVAL)); + while (k--) + complist[k].nr = frac * complist[k].nr; listpos = 0; lastin = -1; /* flag initial sort */ } @@ -268,11 +277,9 @@ sortcomplist() /* fix our list order */ PACKHEAD *list2; register int i; - /* empty queue */ - done_packets(flush_queue()); if (complen <= 0) /* check to see if there is even a list */ return; - if (lastin < 0) /* flag to sort entire list */ + if (lastin < 0 || listpos*4 >= complen*3) qsort((char *)complist, complen, sizeof(PACKHEAD), beamcmp); else if (listpos) { /* else sort and merge sublist */ list2 = (PACKHEAD *)malloc(listpos*sizeof(PACKHEAD)); @@ -284,15 +291,8 @@ sortcomplist() /* fix our list order */ complist+listpos, complen-listpos); free((char *)list2); } - /* check for all finished */ - if (complist[0].nr <= bnrays(hdlist[complist[0].hd],complist[0].bi)) { - free((char *)complist); - complist = NULL; - complen = 0; - } /* drop satisfied requests */ - for (i = complen; i-- && complist[i].nr <= - bnrays(hdlist[complist[i].hd],complist[i].bi); ) + for (i = complen; i-- && complist[i].nr <= complist[i].nc; ) ; if (i < 0) { free((char *)complist); @@ -316,35 +316,29 @@ sortcomplist() /* fix our list order */ * a given bundle to move way down in the computation order. We keep * track of where the computed bundle with the highest priority would end * up, and if we get further in our compute list than this, we resort the - * list and start again from the beginning. We have to flush the queue - * each time we sort, to ensure that we are not disturbing the order. - * If our major assumption is violated, and we have a very steep - * descent in our weights, then we will end up resorting much more often - * than necessary, resulting in frequent flushing of the queue. Since - * a merge sort is used, the sorting costs will be minimal. + * list and start again from the beginning. Since + * a merge sort is used, the sorting costs are minimal. */ next_packet(p) /* prepare packet for computation */ register PACKET *p; { - int ncomp; register int i; - if (complen <= 0) - return(0); if (listpos > lastin) /* time to sort the list */ sortcomplist(); + if (complen <= 0) + return(0); p->hd = complist[listpos].hd; p->bi = complist[listpos].bi; - ncomp = bnrays(hdlist[p->hd],p->bi); - p->nr = complist[listpos].nr - ncomp; + p->nc = complist[listpos].nc; + p->nr = complist[listpos].nr - p->nc; if (p->nr <= 0) return(0); if (p->nr > RPACKSIZ) p->nr = RPACKSIZ; - ncomp += p->nr; /* find where this one would go */ - while (lastin > listpos && complist[listpos].nr * - (bnrays(hdlist[complist[lastin].hd],complist[lastin].bi)+1) - > complist[lastin].nr * (ncomp+1)) + complist[listpos].nc += p->nr; /* find where this one would go */ + while (lastin > listpos && + beamcmp(complist+lastin, complist+listpos) > 0) lastin--; listpos++; return(1);