--- ray/src/hd/rholo3.c 1997/12/01 16:35:02 3.12 +++ ray/src/hd/rholo3.c 1998/01/05 16:46:34 3.18 @@ -1,4 +1,4 @@ -/* Copyright (c) 1997 Silicon Graphics, Inc. */ +/* Copyright (c) 1998 Silicon Graphics, Inc. */ #ifndef lint static char SCCSid[] = "$SunId$ SGI"; @@ -27,13 +27,40 @@ register PACKHEAD *b0, *b1; } +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) { @@ -121,29 +148,17 @@ int nents; default: error(CONSISTENCY, "bundle_set called with unknown operation"); } - if (outdev == NULL) - return; - n = 32*RPACKSIZ; /* allocate packet holder */ - p = (PACKHEAD *)malloc(packsiz(n)); - 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) { - if (b->nrm > n) { - n = b->nrm; - p = (PACKHEAD *)realloc((char *)p, packsiz(n)); - if (p == NULL) - goto memerr; - } - bcopy((char *)hdbray(b), (char *)packra(p), - b->nrm*sizeof(RAYVAL)); - p->hd = clist[i].hd; - p->bi = clist[i].bi; - p->nr = p->nc = b->nrm; - disp_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 */ @@ -158,77 +173,50 @@ memerr: } -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); + /* 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.5*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; - /* free old list */ - if (complen > 0) + int i; + register int j, k; + /* free old list and empty queue */ + if (complen > 0) { free((char *)complist); + done_packets(flush_queue()); + } /* allocate beam list */ complen = 0; for (j = 0; hdlist[j] != NULL; j++) @@ -238,18 +226,18 @@ 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]; + if (frac < 0.) frac = -frac; 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; + complist[k].nc = bnrays(hdlist[j], i); wtotal += complist[k++].nr; } + } /* adjust weights */ if (vdef(DISKSPACE)) frac = 1024.*1024.*vflt(DISKSPACE) / (wtotal*sizeof(RAYVAL)); @@ -262,11 +250,11 @@ init_global() /* initialize global ray computation * mergeclists(cdest, cl1, n1, cl2, n2) /* merge two sorted lists */ -PACKHEAD *cdest; -PACKHEAD *cl1, *cl2; +register PACKHEAD *cdest; +register PACKHEAD *cl1, *cl2; int n1, n2; { - int cmp; + register int cmp; while (n1 | n2) { if (!n1) cmp = 1; @@ -331,8 +319,9 @@ sortcomplist() /* fix our list order */ * 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 */ +next_packet(p, n) /* prepare packet for computation */ register PACKET *p; +int n; { register int i; @@ -346,8 +335,12 @@ register PACKET *p; p->nr = complist[listpos].nr - p->nc; if (p->nr <= 0) return(0); - if (p->nr > RPACKSIZ) - p->nr = RPACKSIZ; +#ifdef DEBUG + if (n < 1 | n > RPACKSIZ) + error(CONSISTENCY, "next_packet called with bad n value"); +#endif + if (p->nr > n) + p->nr = n; complist[listpos].nc += p->nr; /* find where this one would go */ while (lastin > listpos && beamcmp(complist+lastin, complist+listpos) > 0)