| 1 | gregl | 3.1 | /* Copyright (c) 1997 Silicon Graphics, Inc. */ | 
| 2 |  |  |  | 
| 3 |  |  | #ifndef lint | 
| 4 |  |  | static char SCCSid[] = "$SunId$ SGI"; | 
| 5 |  |  | #endif | 
| 6 |  |  |  | 
| 7 |  |  | /* | 
| 8 |  |  | * Routines for tracking beam compuatations | 
| 9 |  |  | */ | 
| 10 |  |  |  | 
| 11 |  |  | #include "rholo.h" | 
| 12 |  |  |  | 
| 13 |  |  |  | 
| 14 |  |  | #define abs(x)          ((x) > 0 ? (x) : -(x)) | 
| 15 |  |  | #define sgn(x)          ((x) > 0 ? 1 : (x) < 0 ? -1 : 0) | 
| 16 |  |  |  | 
| 17 |  |  |  | 
| 18 | gregl | 3.4 | static PACKHEAD *complist=NULL; /* list of beams to compute */ | 
| 19 |  |  | static int      complen=0;      /* length of complist */ | 
| 20 |  |  | static int      listpos=0;      /* current list position for next_packet */ | 
| 21 |  |  | static int      lastin= -1;     /* last ordered position in list */ | 
| 22 | gregl | 3.1 |  | 
| 23 |  |  |  | 
| 24 |  |  | int | 
| 25 | gregl | 3.2 | beamcmp(b0, b1)                 /* comparison for descending compute order */ | 
| 26 |  |  | register PACKHEAD       *b0, *b1; | 
| 27 |  |  | { | 
| 28 |  |  | return( b1->nr*(bnrays(hdlist[b0->hd],b0->bi)+1) - | 
| 29 |  |  | b0->nr*(bnrays(hdlist[b1->hd],b1->bi)+1) ); | 
| 30 |  |  | } | 
| 31 |  |  |  | 
| 32 |  |  |  | 
| 33 |  |  | bundle_set(op, clist, nents)    /* bundle set operation */ | 
| 34 |  |  | int     op; | 
| 35 |  |  | PACKHEAD        *clist; | 
| 36 |  |  | int     nents; | 
| 37 |  |  | { | 
| 38 |  |  | BEAM    *b; | 
| 39 |  |  | PACKHEAD        *p; | 
| 40 |  |  | register int    i, n; | 
| 41 |  |  |  | 
| 42 |  |  | switch (op) { | 
| 43 |  |  | case BS_NEW:                    /* new computation set */ | 
| 44 |  |  | if (complen) | 
| 45 |  |  | free((char *)complist); | 
| 46 |  |  | if (nents <= 0) { | 
| 47 |  |  | complist = NULL; | 
| 48 |  |  | listpos = complen = 0; | 
| 49 |  |  | lastin = -1; | 
| 50 |  |  | return; | 
| 51 |  |  | } | 
| 52 |  |  | complist = (PACKHEAD *)malloc(nents*sizeof(PACKHEAD)); | 
| 53 |  |  | if (complist == NULL) | 
| 54 |  |  | goto memerr; | 
| 55 |  |  | bcopy((char *)clist, (char *)complist, nents*sizeof(PACKHEAD)); | 
| 56 |  |  | complen = nents; | 
| 57 |  |  | listpos = 0; | 
| 58 |  |  | lastin = -1;            /* flag for initial sort */ | 
| 59 |  |  | break; | 
| 60 |  |  | case BS_ADD:                    /* add to computation set */ | 
| 61 |  |  | if (nents <= 0) | 
| 62 |  |  | return; | 
| 63 |  |  | /* merge any common members */ | 
| 64 |  |  | for (i = 0; i < complen; i++) | 
| 65 |  |  | for (n = 0; n < nents; n++) | 
| 66 |  |  | if (clist[n].bi == complist[i].bi && | 
| 67 |  |  | clist[n].hd == complist[i].hd) { | 
| 68 |  |  | complist[i].nr += clist[n].nr; | 
| 69 |  |  | clist[n].nr = 0; | 
| 70 |  |  | lastin = -1;    /* flag full sort */ | 
| 71 |  |  | break; | 
| 72 |  |  | } | 
| 73 |  |  | /* sort updated list */ | 
| 74 |  |  | sortcomplist(); | 
| 75 |  |  | /* sort new entries */ | 
| 76 |  |  | qsort((char *)clist, nents, sizeof(PACKHEAD), beamcmp); | 
| 77 |  |  | /* what can't we satisfy? */ | 
| 78 |  |  | for (n = 0; n < nents && clist[n].nr > | 
| 79 |  |  | bnrays(hdlist[clist[n].hd],clist[n].bi); n++) | 
| 80 |  |  | ; | 
| 81 |  |  | if (n) {                /* allocate space for merged list */ | 
| 82 |  |  | PACKHEAD        *newlist; | 
| 83 |  |  | newlist = (PACKHEAD *)malloc( | 
| 84 |  |  | (complen+n)*sizeof(PACKHEAD) ); | 
| 85 |  |  | if (newlist == NULL) | 
| 86 |  |  | goto memerr; | 
| 87 |  |  | /* merge lists */ | 
| 88 |  |  | mergeclists(newlist, clist, n, complist, complen); | 
| 89 |  |  | if (complen) | 
| 90 |  |  | free((char *)complist); | 
| 91 |  |  | complist = newlist; | 
| 92 |  |  | complen += n; | 
| 93 |  |  | } | 
| 94 |  |  | listpos = 0; | 
| 95 |  |  | lastin = complen-1;     /* list is now sorted */ | 
| 96 |  |  | break; | 
| 97 |  |  | case BS_DEL:                    /* delete from computation set */ | 
| 98 |  |  | if (nents <= 0) | 
| 99 |  |  | return; | 
| 100 |  |  | /* find each member */ | 
| 101 |  |  | for (i = 0; i < complen; i++) | 
| 102 |  |  | for (n = 0; n < nents; n++) | 
| 103 |  |  | if (clist[n].bi == complist[i].bi && | 
| 104 |  |  | clist[n].hd == complist[i].hd) { | 
| 105 |  |  | if (clist[n].nr == 0 || | 
| 106 |  |  | clist[n].nr >= complist[i].nr) | 
| 107 |  |  | complist[i].nr = 0; | 
| 108 |  |  | else | 
| 109 |  |  | complist[i].nr -= clist[n].nr; | 
| 110 |  |  | lastin = -1;    /* flag full sort */ | 
| 111 |  |  | break; | 
| 112 |  |  | } | 
| 113 | gregl | 3.3 | if (lastin < 0)         /* sort updated list */ | 
| 114 |  |  | sortcomplist(); | 
| 115 | gregl | 3.2 | return;                 /* no display */ | 
| 116 |  |  | default: | 
| 117 |  |  | error(CONSISTENCY, "bundle_set called with unknown operation"); | 
| 118 |  |  | } | 
| 119 | gregl | 3.7 | n = RPACKSIZ;                           /* allocate packet holder */ | 
| 120 | gregl | 3.2 | for (i = 0; i < nents; i++) | 
| 121 |  |  | if (clist[i].nr > n) | 
| 122 |  |  | n = clist[i].nr; | 
| 123 | gregl | 3.5 | p = (PACKHEAD *)malloc(packsiz(n)); | 
| 124 | gregl | 3.3 | if (p == NULL) | 
| 125 |  |  | goto memerr; | 
| 126 |  |  | /* display what we have */ | 
| 127 | gregl | 3.2 | for (i = 0; i < nents; i++) | 
| 128 | gregl | 3.3 | if ((b = hdgetbeam(hdlist[clist[i].hd], clist[i].bi)) != NULL) { | 
| 129 | gregl | 3.7 | if (b->nrm > n) { | 
| 130 |  |  | n = b->nrm; | 
| 131 |  |  | p = (PACKHEAD *)realloc((char *)p, packsiz(n)); | 
| 132 |  |  | if (p == NULL) | 
| 133 |  |  | goto memerr; | 
| 134 |  |  | } | 
| 135 | gregl | 3.6 | bcopy((char *)hdbray(b), (char *)packra(p), | 
| 136 | gregl | 3.2 | (p->nr=b->nrm)*sizeof(RAYVAL)); | 
| 137 | gregl | 3.7 | p->hd = clist[i].hd; | 
| 138 |  |  | p->bi = clist[i].bi; | 
| 139 | gregl | 3.5 | disp_packet(p); | 
| 140 | gregl | 3.2 | } | 
| 141 |  |  | free((char *)p);                /* clean up */ | 
| 142 |  |  | return; | 
| 143 |  |  | memerr: | 
| 144 |  |  | error(SYSTEM, "out of memory in bundle_set"); | 
| 145 |  |  | } | 
| 146 |  |  |  | 
| 147 |  |  |  | 
| 148 |  |  | int | 
| 149 | gregl | 3.1 | weightf(hp, x0, x1, x2)         /* voxel weighting function */ | 
| 150 |  |  | register HOLO   *hp; | 
| 151 |  |  | register int    x0, x1, x2; | 
| 152 |  |  | { | 
| 153 |  |  | switch (vlet(OCCUPANCY)) { | 
| 154 |  |  | case 'U':               /* uniform weighting */ | 
| 155 |  |  | return(1); | 
| 156 |  |  | case 'C':               /* center weighting (crude) */ | 
| 157 |  |  | x0 += x0 - hp->grid[0] + 1; | 
| 158 |  |  | x0 = abs(x0)*hp->grid[1]*hp->grid[2]; | 
| 159 |  |  | x1 += x1 - hp->grid[1] + 1; | 
| 160 |  |  | x1 = abs(x1)*hp->grid[0]*hp->grid[2]; | 
| 161 |  |  | x2 += x2 - hp->grid[2] + 1; | 
| 162 |  |  | x2 = abs(x2)*hp->grid[0]*hp->grid[1]; | 
| 163 |  |  | return(hp->grid[0]*hp->grid[1]*hp->grid[2] - | 
| 164 |  |  | (x0+x1+x2)/3); | 
| 165 |  |  | default: | 
| 166 |  |  | badvalue(OCCUPANCY); | 
| 167 |  |  | } | 
| 168 |  |  | } | 
| 169 |  |  |  | 
| 170 |  |  |  | 
| 171 |  |  | /* The following is by Daniel Cohen, taken from Graphics Gems IV, p. 368 */ | 
| 172 |  |  | long | 
| 173 |  |  | lineweight(hp, x, y, z, dx, dy, dz)     /* compute weights along a line */ | 
| 174 |  |  | HOLO    *hp; | 
| 175 |  |  | int     x, y, z, dx, dy, dz; | 
| 176 |  |  | { | 
| 177 |  |  | long    wres = 0; | 
| 178 |  |  | int     n, sx, sy, sz, exy, exz, ezy, ax, ay, az, bx, by, bz; | 
| 179 |  |  |  | 
| 180 |  |  | sx = sgn(dx);   sy = sgn(dy);   sz = sgn(dz); | 
| 181 |  |  | ax = abs(dx);   ay = abs(dy);   az = abs(dz); | 
| 182 |  |  | bx = 2*ax;      by = 2*ay;      bz = 2*az; | 
| 183 |  |  | exy = ay-ax;    exz = az-ax;    ezy = ay-az; | 
| 184 |  |  | n = ax+ay+az + 1;               /* added increment to visit last */ | 
| 185 |  |  | while (n--) { | 
| 186 |  |  | wres += weightf(hp, x, y, z); | 
| 187 |  |  | if (exy < 0) { | 
| 188 |  |  | if (exz < 0) { | 
| 189 |  |  | x += sx; | 
| 190 |  |  | exy += by; exz += bz; | 
| 191 |  |  | } else { | 
| 192 |  |  | z += sz; | 
| 193 |  |  | exz -= bx; ezy += by; | 
| 194 |  |  | } | 
| 195 |  |  | } else { | 
| 196 |  |  | if (ezy < 0) { | 
| 197 |  |  | z += sz; | 
| 198 |  |  | exz -= bx; ezy += by; | 
| 199 |  |  | } else { | 
| 200 |  |  | y += sy; | 
| 201 |  |  | exy -= bx; ezy -= bz; | 
| 202 |  |  | } | 
| 203 |  |  | } | 
| 204 |  |  | } | 
| 205 |  |  | return(wres); | 
| 206 |  |  | } | 
| 207 |  |  |  | 
| 208 |  |  |  | 
| 209 |  |  | init_global()                   /* initialize global ray computation */ | 
| 210 |  |  | { | 
| 211 |  |  | long    wtotal = 0; | 
| 212 |  |  | int     i, j; | 
| 213 |  |  | int     lseg[2][3]; | 
| 214 |  |  | double  frac; | 
| 215 |  |  | register int    k; | 
| 216 | gregl | 3.3 | /* free old list */ | 
| 217 |  |  | if (complen > 0) | 
| 218 |  |  | free((char *)complist); | 
| 219 | gregl | 3.1 | /* allocate beam list */ | 
| 220 |  |  | complen = 0; | 
| 221 |  |  | for (j = 0; hdlist[j] != NULL; j++) | 
| 222 |  |  | complen += nbeams(hdlist[j]); | 
| 223 |  |  | complist = (PACKHEAD *)malloc(complen*sizeof(PACKHEAD)); | 
| 224 |  |  | if (complist == NULL) | 
| 225 |  |  | error(SYSTEM, "out of memory in init_global"); | 
| 226 |  |  | /* compute beam weights */ | 
| 227 |  |  | k = 0; | 
| 228 |  |  | for (j = 0; hdlist[j] != NULL; j++) | 
| 229 |  |  | for (i = nbeams(hdlist[j]); i > 0; i--) { | 
| 230 |  |  | hdlseg(lseg, hdlist[j], i); | 
| 231 |  |  | complist[k].hd = j; | 
| 232 |  |  | complist[k].bi = i; | 
| 233 |  |  | complist[k].nr = lineweight( hdlist[j], | 
| 234 |  |  | lseg[0][0], lseg[0][1], lseg[0][2], | 
| 235 |  |  | lseg[1][0] - lseg[0][0], | 
| 236 |  |  | lseg[1][1] - lseg[0][1], | 
| 237 |  |  | lseg[1][2] - lseg[0][2] ); | 
| 238 |  |  | wtotal += complist[k++].nr; | 
| 239 |  |  | } | 
| 240 |  |  | /* adjust weights */ | 
| 241 |  |  | if (vdef(DISKSPACE)) { | 
| 242 |  |  | frac = 1024.*1024.*vflt(DISKSPACE) / (wtotal*sizeof(RAYVAL)); | 
| 243 |  |  | if (frac < 0.95 | frac > 1.05) | 
| 244 |  |  | while (k--) | 
| 245 |  |  | complist[k].nr = frac * complist[k].nr; | 
| 246 |  |  | } | 
| 247 | gregl | 3.4 | listpos = 0; lastin = -1;       /* flag initial sort */ | 
| 248 | gregl | 3.1 | } | 
| 249 |  |  |  | 
| 250 |  |  |  | 
| 251 |  |  | mergeclists(cdest, cl1, n1, cl2, n2)    /* merge two sorted lists */ | 
| 252 |  |  | PACKHEAD        *cdest; | 
| 253 |  |  | PACKHEAD        *cl1, *cl2; | 
| 254 |  |  | int     n1, n2; | 
| 255 |  |  | { | 
| 256 |  |  | int     cmp; | 
| 257 |  |  |  | 
| 258 |  |  | while (n1 | n2) { | 
| 259 |  |  | if (!n1) cmp = 1; | 
| 260 |  |  | else if (!n2) cmp = -1; | 
| 261 |  |  | else cmp = beamcmp(cl1, cl2); | 
| 262 |  |  | if (cmp > 0) { | 
| 263 |  |  | copystruct(cdest, cl2); | 
| 264 |  |  | cl2++; n2--; | 
| 265 |  |  | } else { | 
| 266 |  |  | copystruct(cdest, cl1); | 
| 267 |  |  | cl1++; n1--; | 
| 268 |  |  | } | 
| 269 |  |  | cdest++; | 
| 270 |  |  | } | 
| 271 |  |  | } | 
| 272 |  |  |  | 
| 273 |  |  |  | 
| 274 |  |  | sortcomplist()                  /* fix our list order */ | 
| 275 |  |  | { | 
| 276 |  |  | PACKHEAD        *list2; | 
| 277 | gregl | 3.2 | register int    i; | 
| 278 |  |  |  | 
| 279 | gregl | 3.1 | /* empty queue */ | 
| 280 |  |  | done_packets(flush_queue()); | 
| 281 | gregl | 3.3 | if (complen <= 0)       /* check to see if there is even a list */ | 
| 282 | gregl | 3.2 | return; | 
| 283 | gregl | 3.6 | if (lastin < 0 || listpos*4 >= complen*3) | 
| 284 | gregl | 3.1 | qsort((char *)complist, complen, sizeof(PACKHEAD), beamcmp); | 
| 285 |  |  | else if (listpos) {     /* else sort and merge sublist */ | 
| 286 |  |  | list2 = (PACKHEAD *)malloc(listpos*sizeof(PACKHEAD)); | 
| 287 |  |  | if (list2 == NULL) | 
| 288 |  |  | error(SYSTEM, "out of memory in sortcomplist"); | 
| 289 |  |  | bcopy((char *)complist,(char *)list2,listpos*sizeof(PACKHEAD)); | 
| 290 |  |  | qsort((char *)list2, listpos, sizeof(PACKHEAD), beamcmp); | 
| 291 |  |  | mergeclists(complist, list2, listpos, | 
| 292 |  |  | complist+listpos, complen-listpos); | 
| 293 |  |  | free((char *)list2); | 
| 294 |  |  | } | 
| 295 | gregl | 3.2 | /* check for all finished */ | 
| 296 |  |  | if (complist[0].nr <= bnrays(hdlist[complist[0].hd],complist[0].bi)) { | 
| 297 |  |  | free((char *)complist); | 
| 298 |  |  | complist = NULL; | 
| 299 |  |  | complen = 0; | 
| 300 |  |  | } | 
| 301 |  |  | /* drop satisfied requests */ | 
| 302 |  |  | for (i = complen; i-- && complist[i].nr <= | 
| 303 |  |  | bnrays(hdlist[complist[i].hd],complist[i].bi); ) | 
| 304 |  |  | ; | 
| 305 | gregl | 3.4 | if (i < 0) { | 
| 306 |  |  | free((char *)complist); | 
| 307 |  |  | complist = NULL; | 
| 308 |  |  | complen = 0; | 
| 309 |  |  | } else if (i < complen-1) { | 
| 310 | gregl | 3.2 | list2 = (PACKHEAD *)realloc((char *)complist, | 
| 311 |  |  | (i+1)*sizeof(PACKHEAD)); | 
| 312 |  |  | if (list2 != NULL) { | 
| 313 |  |  | complist = list2; | 
| 314 |  |  | complen = i+1; | 
| 315 |  |  | } | 
| 316 |  |  | } | 
| 317 |  |  | listpos = 0; lastin = i; | 
| 318 | gregl | 3.1 | } | 
| 319 |  |  |  | 
| 320 |  |  |  | 
| 321 |  |  | /* | 
| 322 |  |  | * The following routine works on the assumption that the bundle weights are | 
| 323 |  |  | * more or less evenly distributed, such that computing a packet causes | 
| 324 |  |  | * a given bundle to move way down in the computation order.  We keep | 
| 325 |  |  | * track of where the computed bundle with the highest priority would end | 
| 326 |  |  | * up, and if we get further in our compute list than this, we resort the | 
| 327 |  |  | * list and start again from the beginning.  We have to flush the queue | 
| 328 |  |  | * each time we sort, to ensure that we are not disturbing the order. | 
| 329 |  |  | *      If our major assumption is violated, and we have a very steep | 
| 330 |  |  | * descent in our weights, then we will end up resorting much more often | 
| 331 |  |  | * than necessary, resulting in frequent flushing of the queue.  Since | 
| 332 |  |  | * a merge sort is used, the sorting costs will be minimal. | 
| 333 |  |  | */ | 
| 334 |  |  | next_packet(p)                  /* prepare packet for computation */ | 
| 335 |  |  | register PACKET *p; | 
| 336 |  |  | { | 
| 337 |  |  | int     ncomp; | 
| 338 |  |  | register int    i; | 
| 339 |  |  |  | 
| 340 |  |  | if (complen <= 0) | 
| 341 |  |  | return(0); | 
| 342 |  |  | if (listpos > lastin)           /* time to sort the list */ | 
| 343 |  |  | sortcomplist(); | 
| 344 |  |  | p->hd = complist[listpos].hd; | 
| 345 |  |  | p->bi = complist[listpos].bi; | 
| 346 |  |  | ncomp = bnrays(hdlist[p->hd],p->bi); | 
| 347 |  |  | p->nr = complist[listpos].nr - ncomp; | 
| 348 |  |  | if (p->nr <= 0) | 
| 349 |  |  | return(0); | 
| 350 |  |  | if (p->nr > RPACKSIZ) | 
| 351 |  |  | p->nr = RPACKSIZ; | 
| 352 |  |  | ncomp += p->nr;                 /* find where this one would go */ | 
| 353 |  |  | while (lastin > listpos && complist[listpos].nr * | 
| 354 |  |  | (bnrays(hdlist[complist[lastin].hd],complist[lastin].bi)+1) | 
| 355 |  |  | > complist[lastin].nr * (ncomp+1)) | 
| 356 |  |  | lastin--; | 
| 357 |  |  | listpos++; | 
| 358 |  |  | return(1); | 
| 359 |  |  | } |