ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rholo3.c
Revision: 3.33
Committed: Mon Mar 8 17:31:49 1999 UTC (25 years ago) by gwlarson
Content type: text/plain
Branch: MAIN
Changes since 3.32: +1 -2 lines
Log Message:
removed unnecessary inclusion of sys/types.h

File Contents

# Content
1 /* Copyright (c) 1999 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 #include "view.h"
13
14 #ifndef NFRAG2CHUNK
15 #define NFRAG2CHUNK 4096 /* number of fragments to start chunking */
16 #endif
17
18 #ifndef abs
19 #define abs(x) ((x) > 0 ? (x) : -(x))
20 #endif
21 #ifndef sgn
22 #define sgn(x) ((x) > 0 ? 1 : (x) < 0 ? -1 : 0)
23 #endif
24
25 #define rchunk(n) (((n)+(RPACKSIZ/2))/RPACKSIZ)
26
27 extern time_t time();
28
29 int chunkycmp = 0; /* clump beams together on disk */
30
31 static PACKHEAD *complist=NULL; /* list of beams to compute */
32 static int complen=0; /* length of complist */
33 static int listpos=0; /* current list position for next_packet */
34 static int lastin= -1; /* last ordered position in list */
35
36
37 int
38 beamcmp(b0, b1) /* comparison for compute order */
39 register PACKHEAD *b0, *b1;
40 {
41 BEAMI *bip0, *bip1;
42 register long c;
43 /* first check desired quantities */
44 if (chunkycmp)
45 c = rchunk(b1->nr)*(rchunk(b0->nc)+1L) -
46 rchunk(b0->nr)*(rchunk(b1->nc)+1L);
47 else
48 c = b1->nr*(b0->nc+1L) - b0->nr*(b1->nc+1L);
49 if (c > 0) return(1);
50 if (c < 0) return(-1);
51 /* only one file, so skip the following: */
52 #if 0
53 /* next, check file descriptors */
54 c = hdlist[b0->hd]->fd - hdlist[b1->hd]->fd;
55 if (c) return(c);
56 #endif
57 /* finally, check file positions */
58 bip0 = &hdlist[b0->hd]->bi[b0->bi];
59 bip1 = &hdlist[b1->hd]->bi[b1->bi];
60 /* put diskless beams last */
61 if (!bip0->nrd)
62 return(bip1->nrd > 0);
63 if (!bip1->nrd)
64 return(-1);
65 c = bip0->fo - bip1->fo;
66 return(c < 0 ? -1 : c > 0);
67 }
68
69
70 int
71 beamidcmp(b0, b1) /* comparison for beam searching */
72 register PACKHEAD *b0, *b1;
73 {
74 register int c = b0->hd - b1->hd;
75
76 if (c) return(c);
77 return(b0->bi - b1->bi);
78 }
79
80
81 int
82 dispbeam(b, hb) /* display a holodeck beam */
83 register BEAM *b;
84 register HDBEAMI *hb;
85 {
86 static int n = 0;
87 static PACKHEAD *p = NULL;
88
89 if (b == NULL)
90 return;
91 if (b->nrm > n) { /* (re)allocate packet holder */
92 n = b->nrm;
93 if (p == NULL) p = (PACKHEAD *)malloc(packsiz(n));
94 else p = (PACKHEAD *)realloc((char *)p, packsiz(n));
95 CHECK(p==NULL, SYSTEM, "out of memory in dispbeam");
96 }
97 /* assign packet fields */
98 bcopy((char *)hdbray(b), (char *)packra(p), b->nrm*sizeof(RAYVAL));
99 p->nr = p->nc = b->nrm;
100 for (p->hd = 0; hdlist[p->hd] != hb->h; p->hd++)
101 if (hdlist[p->hd] == NULL)
102 error(CONSISTENCY, "unregistered holodeck in dispbeam");
103 p->bi = hb->b;
104 disp_packet(p); /* display it */
105 if (n >= 1024) { /* free ridiculous packets */
106 free((char *)p);
107 p = NULL; n = 0;
108 }
109 }
110
111
112 bundle_set(op, clist, nents) /* bundle set operation */
113 int op;
114 PACKHEAD *clist;
115 int nents;
116 {
117 int oldnr, n;
118 HDBEAMI *hbarr;
119 register PACKHEAD *csm;
120 register int i;
121 /* search for common members */
122 for (csm = clist+nents; csm-- > clist; )
123 csm->nc = -1;
124 qsort((char *)clist, nents, sizeof(PACKHEAD), beamidcmp);
125 for (i = 0; i < complen; i++) {
126 csm = (PACKHEAD *)bsearch((char *)(complist+i), (char *)clist,
127 nents, sizeof(PACKHEAD), beamidcmp);
128 if (csm == NULL)
129 continue;
130 oldnr = complist[i].nr;
131 csm->nc = complist[i].nc;
132 switch (op) {
133 case BS_ADD: /* add to count */
134 complist[i].nr += csm->nr;
135 csm->nr = 0;
136 break;
137 case BS_MAX: /* maximum of counts */
138 if (csm->nr > complist[i].nr)
139 complist[i].nr = csm->nr;
140 csm->nr = 0;
141 break;
142 case BS_ADJ: /* reset count */
143 complist[i].nr = csm->nr;
144 csm->nr = 0;
145 break;
146 case BS_DEL: /* delete count */
147 if (csm->nr == 0 || csm->nr >= complist[i].nr)
148 complist[i].nr = 0;
149 else
150 complist[i].nr -= csm->nr;
151 break;
152 }
153 if (complist[i].nr != oldnr)
154 lastin = -1; /* flag sort */
155 }
156 /* record computed rays for uncommon beams */
157 for (csm = clist+nents; csm-- > clist; )
158 if (csm->nc < 0)
159 csm->nc = bnrays(hdlist[csm->hd], csm->bi);
160 /* complete list operations */
161 switch (op) {
162 case BS_NEW: /* new computation set */
163 listpos = 0; lastin = -1;
164 if (complen) /* free old list */
165 free((char *)complist);
166 complist = NULL;
167 if (!(complen = nents))
168 return;
169 complist = (PACKHEAD *)malloc(nents*sizeof(PACKHEAD));
170 if (complist == NULL)
171 goto memerr;
172 bcopy((char *)clist, (char *)complist, nents*sizeof(PACKHEAD));
173 break;
174 case BS_ADD: /* add to computation set */
175 case BS_MAX: /* maximum of quantities */
176 case BS_ADJ: /* adjust set quantities */
177 if (nents <= 0)
178 return;
179 sortcomplist(); /* sort updated list & new entries */
180 qsort((char *)clist, nents, sizeof(PACKHEAD), beamcmp);
181 /* what can't we satisfy? */
182 for (i = nents, csm = clist; i-- && csm->nr > csm->nc; csm++)
183 ;
184 n = csm - clist;
185 if (op != BS_ADD) { /* don't regenerate adjusted beams */
186 for (++i; i-- && csm->nr > 0; csm++)
187 ;
188 nents = csm - clist;
189 }
190 if (n) { /* allocate space for merged list */
191 PACKHEAD *newlist;
192 newlist = (PACKHEAD *)malloc(
193 (complen+n)*sizeof(PACKHEAD) );
194 if (newlist == NULL)
195 goto memerr;
196 /* merge lists */
197 mergeclists(newlist, clist, n, complist, complen);
198 if (complen)
199 free((char *)complist);
200 complist = newlist;
201 complen += n;
202 }
203 listpos = 0;
204 lastin = complen-1; /* list is now sorted */
205 break;
206 case BS_DEL: /* delete from computation set */
207 return; /* already done */
208 default:
209 error(CONSISTENCY, "bundle_set called with unknown operation");
210 }
211 if (outdev == NULL || !nents) /* nothing to display? */
212 return;
213 /* load and display beams we have */
214 hbarr = (HDBEAMI *)malloc(nents*sizeof(HDBEAMI));
215 for (i = nents; i--; ) {
216 hbarr[i].h = hdlist[clist[i].hd];
217 hbarr[i].b = clist[i].bi;
218 }
219 hdloadbeams(hbarr, nents, dispbeam);
220 free((char *)hbarr);
221 if (hdfragflags&FF_READ) {
222 listpos = 0;
223 lastin = -1; /* need to re-sort list */
224 }
225 return;
226 memerr:
227 error(SYSTEM, "out of memory in bundle_set");
228 }
229
230
231 double
232 beamvolume(hp, bi) /* compute approximate volume of a beam */
233 HOLO *hp;
234 int bi;
235 {
236 GCOORD gc[2];
237 FVECT cp[4], edgeA, edgeB, cent[2];
238 FVECT v, crossp[2], diffv;
239 double vol[2];
240 register int i;
241 /* get grid coordinates */
242 if (!hdbcoord(gc, hp, bi))
243 error(CONSISTENCY, "bad beam index in beamvolume");
244 for (i = 0; i < 2; i++) { /* compute cell area vectors */
245 hdcell(cp, hp, gc+i);
246 VSUM(edgeA, cp[1], cp[0], -1.0);
247 VSUM(edgeB, cp[3], cp[1], -1.0);
248 fcross(crossp[i], edgeA, edgeB);
249 /* compute center */
250 cent[i][0] = 0.5*(cp[0][0] + cp[2][0]);
251 cent[i][1] = 0.5*(cp[0][1] + cp[2][1]);
252 cent[i][2] = 0.5*(cp[0][2] + cp[2][2]);
253 }
254 /* compute difference vector */
255 VSUM(diffv, cent[1], cent[0], -1.0);
256 for (i = 0; i < 2; i++) { /* compute volume contributions */
257 vol[i] = 0.5*DOT(crossp[i], diffv);
258 if (vol[i] < 0.) vol[i] = -vol[i];
259 }
260 return(vol[0] + vol[1]); /* return total volume */
261 }
262
263
264 ambient_list() /* compute ambient beam list */
265 {
266 int4 wtotal, minrt;
267 double frac;
268 int i;
269 register int j, k;
270
271 complen = 0;
272 for (j = 0; hdlist[j] != NULL; j++)
273 complen += nbeams(hdlist[j]);
274 complist = (PACKHEAD *)malloc(complen*sizeof(PACKHEAD));
275 CHECK(complist==NULL, SYSTEM, "out of memory in ambient_list");
276 /* compute beam weights */
277 k = 0; wtotal = 0;
278 for (j = 0; hdlist[j] != NULL; j++) {
279 frac = 512. * VLEN(hdlist[j]->wg[0]) *
280 VLEN(hdlist[j]->wg[1]) *
281 VLEN(hdlist[j]->wg[2]);
282 for (i = nbeams(hdlist[j]); i > 0; i--) {
283 complist[k].hd = j;
284 complist[k].bi = i;
285 complist[k].nr = frac*beamvolume(hdlist[j], i) + 0.5;
286 complist[k].nc = bnrays(hdlist[j], i);
287 wtotal += complist[k++].nr;
288 }
289 }
290 /* adjust sample weights */
291 if (vdef(DISKSPACE))
292 frac = 1024.*1024.*vflt(DISKSPACE) / (wtotal*sizeof(RAYVAL));
293 else
294 frac = 1024.*1024.*2048. / (wtotal*sizeof(RAYVAL));
295 minrt = .02*frac*wtotal/complen + .5; /* heuristic mimimum */
296 if (minrt > RPACKSIZ)
297 minrt = RPACKSIZ;
298 for (k = complen; k--; )
299 if ((complist[k].nr = frac*complist[k].nr + 0.5) < minrt)
300 complist[k].nr = minrt;
301 listpos = 0; lastin = -1; /* flag initial sort */
302 }
303
304
305 view_list(fp) /* assign beam priority from view list */
306 FILE *fp;
307 {
308 double pa = 1.;
309 VIEW curview;
310 int xr, yr;
311 char *err;
312 BEAMLIST blist;
313
314 copystruct(&curview, &stdview);
315 while (nextview(&curview, fp) != EOF) {
316 if ((err = setview(&curview)) != NULL) {
317 error(WARNING, err);
318 continue;
319 }
320 xr = yr = 1024;
321 normaspect(viewaspect(&curview), &pa, &xr, &yr);
322 viewbeams(&curview, xr, yr, &blist);
323 bundle_set(BS_MAX, blist.bl, blist.nb);
324 free((char *)blist.bl);
325 }
326 }
327
328
329 init_global() /* initialize global ray computation */
330 {
331 register int k;
332 /* free old list and empty queue */
333 if (complen > 0) {
334 free((char *)complist);
335 done_packets(flush_queue());
336 }
337 /* reseed random number generator */
338 srandom(time(NULL));
339 /* allocate beam list */
340 if (readinp)
341 view_list(stdin);
342 else
343 ambient_list();
344 /* no view vicinity */
345 myeye.rng = 0;
346 }
347
348
349 mergeclists(cdest, cl1, n1, cl2, n2) /* merge two sorted lists */
350 register PACKHEAD *cdest;
351 register PACKHEAD *cl1, *cl2;
352 int n1, n2;
353 {
354 register int cmp;
355
356 while (n1 | n2) {
357 if (!n1) cmp = 1;
358 else if (!n2) cmp = -1;
359 else cmp = beamcmp(cl1, cl2);
360 if (cmp > 0) {
361 copystruct(cdest, cl2);
362 cl2++; n2--;
363 } else {
364 copystruct(cdest, cl1);
365 cl1++; n1--;
366 }
367 cdest++;
368 }
369 }
370
371
372 sortcomplist() /* fix our list order */
373 {
374 PACKHEAD *list2;
375 int listlen;
376 register int i;
377
378 if (complen <= 0) /* check to see if there is even a list */
379 return;
380 if (!chunkycmp) /* check to see if fragment list is full */
381 if (!hdfragOK(hdlist[0]->fd, &listlen, NULL)
382 #if NFRAG2CHUNK
383 || listlen >= NFRAG2CHUNK
384 #endif
385 ) {
386 chunkycmp++; /* use "chunky" comparison */
387 lastin = -1; /* need to re-sort list */
388 #ifdef DEBUG
389 error(WARNING, "using chunky comparison mode");
390 #endif
391 }
392 if (lastin < 0 || listpos*4 >= complen*3)
393 qsort((char *)complist, complen, sizeof(PACKHEAD), beamcmp);
394 else if (listpos) { /* else sort and merge sublist */
395 list2 = (PACKHEAD *)malloc(listpos*sizeof(PACKHEAD));
396 CHECK(list2==NULL, SYSTEM, "out of memory in sortcomplist");
397 bcopy((char *)complist,(char *)list2,listpos*sizeof(PACKHEAD));
398 qsort((char *)list2, listpos, sizeof(PACKHEAD), beamcmp);
399 mergeclists(complist, list2, listpos,
400 complist+listpos, complen-listpos);
401 free((char *)list2);
402 }
403 /* drop satisfied requests */
404 for (i = complen; i-- && complist[i].nr <= complist[i].nc; )
405 ;
406 if (i < 0) {
407 free((char *)complist);
408 complist = NULL;
409 complen = 0;
410 } else if (i < complen-1) {
411 list2 = (PACKHEAD *)realloc((char *)complist,
412 (i+1)*sizeof(PACKHEAD));
413 if (list2 != NULL)
414 complist = list2;
415 complen = i+1;
416 }
417 listpos = 0; lastin = i;
418 }
419
420
421 /*
422 * The following routine works on the assumption that the bundle weights are
423 * more or less evenly distributed, such that computing a packet causes
424 * a given bundle to move way down in the computation order. We keep
425 * track of where the computed bundle with the highest priority would end
426 * up, and if we get further in our compute list than this, we re-sort the
427 * list and start again from the beginning. Since
428 * a merge sort is used, the sorting costs are minimal.
429 */
430 next_packet(p, n) /* prepare packet for computation */
431 register PACKET *p;
432 int n;
433 {
434 register int i;
435
436 if (listpos > lastin) /* time to sort the list */
437 sortcomplist();
438 if (complen <= 0)
439 return(0);
440 p->hd = complist[listpos].hd;
441 p->bi = complist[listpos].bi;
442 p->nc = complist[listpos].nc;
443 p->nr = complist[listpos].nr - p->nc;
444 if (p->nr <= 0)
445 return(0);
446 DCHECK(n < 1 | n > RPACKSIZ,
447 CONSISTENCY, "next_packet called with bad n value");
448 if (p->nr > n)
449 p->nr = n;
450 complist[listpos].nc += p->nr; /* find where this one would go */
451 if (hdgetbeam(hdlist[p->hd], p->bi) != NULL)
452 hdfreefrag(hdlist[p->hd], p->bi);
453 while (lastin > listpos &&
454 beamcmp(complist+lastin, complist+listpos) > 0)
455 lastin--;
456 listpos++;
457 return(1);
458 }