ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rholo3.c
Revision: 3.32
Committed: Thu Mar 4 15:44:10 1999 UTC (25 years, 1 month ago) by gwlarson
Content type: text/plain
Branch: MAIN
Changes since 3.31: +10 -8 lines
Log Message:
put in mimimum ray count for ambient_list()

File Contents

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