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

# User Rev Content
1 gwlarson 3.23 /* Copyright (c) 1998 Silicon Graphics, Inc. */
2    
3     #ifndef lint
4     static char SCCSid[] = "$SunId$ SGI";
5     #endif
6    
7 gregl 3.1 /*
8     * Routines for tracking beam compuatations
9     */
10    
11     #include "rholo.h"
12 gwlarson 3.31 #include "view.h"
13 gwlarson 3.28 #include <sys/types.h>
14 gregl 3.1
15 gwlarson 3.27 #ifndef NFRAG2CHUNK
16     #define NFRAG2CHUNK 4096 /* number of fragments to start chunking */
17     #endif
18    
19     #ifndef abs
20 gregl 3.1 #define abs(x) ((x) > 0 ? (x) : -(x))
21 gwlarson 3.27 #endif
22     #ifndef sgn
23 gregl 3.1 #define sgn(x) ((x) > 0 ? 1 : (x) < 0 ? -1 : 0)
24 gwlarson 3.27 #endif
25 gregl 3.1
26 gwlarson 3.27 #define rchunk(n) (((n)+(RPACKSIZ/2))/RPACKSIZ)
27    
28 gwlarson 3.28 extern time_t time();
29    
30 gwlarson 3.30 int chunkycmp = 0; /* clump beams together on disk */
31    
32 gregl 3.4 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 gregl 3.1
37    
38     int
39 gwlarson 3.23 beamcmp(b0, b1) /* comparison for compute order */
40 gregl 3.2 register PACKHEAD *b0, *b1;
41     {
42 gwlarson 3.27 BEAMI *bip0, *bip1;
43     register long c;
44     /* first check desired quantities */
45 gwlarson 3.30 if (chunkycmp)
46 gwlarson 3.28 c = rchunk(b1->nr)*(rchunk(b0->nc)+1L) -
47     rchunk(b0->nr)*(rchunk(b1->nc)+1L);
48 gwlarson 3.27 else
49 gwlarson 3.28 c = b1->nr*(b0->nc+1L) - b0->nr*(b1->nc+1L);
50     if (c > 0) return(1);
51     if (c < 0) return(-1);
52 gwlarson 3.27 /* 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 gregl 3.2 }
69    
70    
71 gregl 3.14 int
72 gwlarson 3.23 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 gregl 3.14 register BEAM *b;
85 gwlarson 3.23 register HDBEAMI *hb;
86 gregl 3.14 {
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 gwlarson 3.31 CHECK(p==NULL, SYSTEM, "out of memory in dispbeam");
97 gregl 3.14 }
98     /* assign packet fields */
99     bcopy((char *)hdbray(b), (char *)packra(p), b->nrm*sizeof(RAYVAL));
100     p->nr = p->nc = b->nrm;
101 gwlarson 3.23 for (p->hd = 0; hdlist[p->hd] != hb->h; p->hd++)
102 gregl 3.14 if (hdlist[p->hd] == NULL)
103     error(CONSISTENCY, "unregistered holodeck in dispbeam");
104 gwlarson 3.23 p->bi = hb->b;
105 gregl 3.14 disp_packet(p); /* display it */
106 gwlarson 3.29 if (n >= 1024) { /* free ridiculous packets */
107     free((char *)p);
108     p = NULL; n = 0;
109     }
110 gregl 3.14 }
111    
112    
113 gregl 3.2 bundle_set(op, clist, nents) /* bundle set operation */
114     int op;
115 gwlarson 3.23 PACKHEAD *clist;
116 gregl 3.2 int nents;
117     {
118 gwlarson 3.23 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 gwlarson 3.25 qsort((char *)clist, nents, sizeof(PACKHEAD), beamidcmp);
126 gwlarson 3.23 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 gwlarson 3.31 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 gwlarson 3.23 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 gregl 3.20 }
157 gwlarson 3.24 /* record computed rays for uncommon beams */
158 gwlarson 3.23 for (csm = clist+nents; csm-- > clist; )
159     if (csm->nc < 0)
160     csm->nc = bnrays(hdlist[csm->hd], csm->bi);
161 gregl 3.20 /* complete list operations */
162 gregl 3.2 switch (op) {
163     case BS_NEW: /* new computation set */
164 gwlarson 3.21 listpos = 0; lastin = -1;
165 gregl 3.20 if (complen) /* free old list */
166 gregl 3.2 free((char *)complist);
167 gregl 3.20 complist = NULL;
168     if (!(complen = nents))
169 gregl 3.2 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 gwlarson 3.31 case BS_MAX: /* maximum of quantities */
177 gregl 3.11 case BS_ADJ: /* adjust set quantities */
178 gregl 3.2 if (nents <= 0)
179     return;
180 gregl 3.20 sortcomplist(); /* sort updated list & new entries */
181 gregl 3.2 qsort((char *)clist, nents, sizeof(PACKHEAD), beamcmp);
182     /* what can't we satisfy? */
183 gwlarson 3.23 for (i = nents, csm = clist; i-- && csm->nr > csm->nc; csm++)
184 gregl 3.2 ;
185 gwlarson 3.23 n = csm - clist;
186 gwlarson 3.31 if (op != BS_ADD) { /* don't regenerate adjusted beams */
187 gwlarson 3.24 for (++i; i-- && csm->nr > 0; csm++)
188 gwlarson 3.22 ;
189 gwlarson 3.24 nents = csm - clist;
190 gwlarson 3.22 }
191 gregl 3.2 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 gregl 3.20 return; /* already done */
209 gregl 3.2 default:
210     error(CONSISTENCY, "bundle_set called with unknown operation");
211     }
212 gwlarson 3.27 if (outdev == NULL || !nents) /* nothing to display? */
213 gwlarson 3.22 return;
214     /* load and display beams we have */
215 gwlarson 3.23 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 gregl 3.11 }
220 gwlarson 3.23 hdloadbeams(hbarr, nents, dispbeam);
221     free((char *)hbarr);
222 gwlarson 3.27 if (hdfragflags&FF_READ) {
223     listpos = 0;
224     lastin = -1; /* need to re-sort list */
225     }
226 gregl 3.2 return;
227     memerr:
228     error(SYSTEM, "out of memory in bundle_set");
229     }
230    
231    
232 gregl 3.13 double
233     beamvolume(hp, bi) /* compute approximate volume of a beam */
234 gregl 3.1 HOLO *hp;
235 gregl 3.13 int bi;
236 gregl 3.1 {
237 gregl 3.13 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 gregl 3.1 }
255 gregl 3.13 /* compute difference vector */
256     VSUM(diffv, cent[1], cent[0], -1.0);
257     for (i = 0; i < 2; i++) { /* compute volume contributions */
258 gregl 3.15 vol[i] = 0.5*DOT(crossp[i], diffv);
259 gregl 3.13 if (vol[i] < 0.) vol[i] = -vol[i];
260     }
261     return(vol[0] + vol[1]); /* return total volume */
262 gregl 3.1 }
263    
264    
265 gwlarson 3.31 ambient_list() /* compute ambient beam list */
266 gregl 3.1 {
267 gwlarson 3.32 int4 wtotal, minrt;
268 gregl 3.1 double frac;
269 gregl 3.13 int i;
270     register int j, k;
271 gwlarson 3.31
272 gregl 3.1 complen = 0;
273     for (j = 0; hdlist[j] != NULL; j++)
274     complen += nbeams(hdlist[j]);
275     complist = (PACKHEAD *)malloc(complen*sizeof(PACKHEAD));
276 gwlarson 3.31 CHECK(complist==NULL, SYSTEM, "out of memory in ambient_list");
277 gregl 3.1 /* compute beam weights */
278 gwlarson 3.32 k = 0; wtotal = 0;
279 gregl 3.13 for (j = 0; hdlist[j] != NULL; j++) {
280 gregl 3.19 frac = 512. * VLEN(hdlist[j]->wg[0]) *
281     VLEN(hdlist[j]->wg[1]) *
282     VLEN(hdlist[j]->wg[2]);
283 gregl 3.1 for (i = nbeams(hdlist[j]); i > 0; i--) {
284     complist[k].hd = j;
285     complist[k].bi = i;
286 gregl 3.13 complist[k].nr = frac*beamvolume(hdlist[j], i) + 0.5;
287 gregl 3.18 complist[k].nc = bnrays(hdlist[j], i);
288 gregl 3.1 wtotal += complist[k++].nr;
289     }
290 gregl 3.13 }
291 gwlarson 3.31 /* adjust sample weights */
292 gregl 3.12 if (vdef(DISKSPACE))
293 gwlarson 3.32 frac = 1024.*1024.*vflt(DISKSPACE) / (wtotal*sizeof(RAYVAL));
294 gregl 3.12 else
295 gwlarson 3.32 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 gwlarson 3.31 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 gwlarson 3.25 /* no view vicinity */
346     myeye.rng = 0;
347 gregl 3.1 }
348    
349    
350     mergeclists(cdest, cl1, n1, cl2, n2) /* merge two sorted lists */
351 gregl 3.16 register PACKHEAD *cdest;
352     register PACKHEAD *cl1, *cl2;
353 gregl 3.1 int n1, n2;
354     {
355 gregl 3.16 register int cmp;
356 gregl 3.1
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 gwlarson 3.27 int listlen;
377 gregl 3.2 register int i;
378    
379 gregl 3.3 if (complen <= 0) /* check to see if there is even a list */
380 gregl 3.2 return;
381 gwlarson 3.30 if (!chunkycmp) /* check to see if fragment list is full */
382 gwlarson 3.27 if (!hdfragOK(hdlist[0]->fd, &listlen, NULL)
383     #if NFRAG2CHUNK
384     || listlen >= NFRAG2CHUNK
385     #endif
386     ) {
387 gwlarson 3.30 chunkycmp++; /* use "chunky" comparison */
388     lastin = -1; /* need to re-sort list */
389 gwlarson 3.27 #ifdef DEBUG
390     error(WARNING, "using chunky comparison mode");
391     #endif
392     }
393 gregl 3.6 if (lastin < 0 || listpos*4 >= complen*3)
394 gregl 3.1 qsort((char *)complist, complen, sizeof(PACKHEAD), beamcmp);
395     else if (listpos) { /* else sort and merge sublist */
396     list2 = (PACKHEAD *)malloc(listpos*sizeof(PACKHEAD));
397 gwlarson 3.31 CHECK(list2==NULL, SYSTEM, "out of memory in sortcomplist");
398 gregl 3.1 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 gregl 3.2 /* drop satisfied requests */
405 gregl 3.11 for (i = complen; i-- && complist[i].nr <= complist[i].nc; )
406 gregl 3.2 ;
407 gregl 3.4 if (i < 0) {
408     free((char *)complist);
409     complist = NULL;
410     complen = 0;
411     } else if (i < complen-1) {
412 gregl 3.2 list2 = (PACKHEAD *)realloc((char *)complist,
413     (i+1)*sizeof(PACKHEAD));
414 gwlarson 3.26 if (list2 != NULL)
415 gregl 3.2 complist = list2;
416 gwlarson 3.26 complen = i+1;
417 gregl 3.2 }
418     listpos = 0; lastin = i;
419 gregl 3.1 }
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 gwlarson 3.27 * up, and if we get further in our compute list than this, we re-sort the
428 gregl 3.11 * list and start again from the beginning. Since
429     * a merge sort is used, the sorting costs are minimal.
430 gregl 3.1 */
431 gregl 3.17 next_packet(p, n) /* prepare packet for computation */
432 gregl 3.1 register PACKET *p;
433 gregl 3.17 int n;
434 gregl 3.1 {
435     register int i;
436    
437 gregl 3.10 if (listpos > lastin) /* time to sort the list */
438     sortcomplist();
439 gregl 3.1 if (complen <= 0)
440     return(0);
441     p->hd = complist[listpos].hd;
442     p->bi = complist[listpos].bi;
443 gregl 3.11 p->nc = complist[listpos].nc;
444     p->nr = complist[listpos].nr - p->nc;
445 gregl 3.1 if (p->nr <= 0)
446     return(0);
447 gwlarson 3.28 DCHECK(n < 1 | n > RPACKSIZ,
448     CONSISTENCY, "next_packet called with bad n value");
449 gregl 3.17 if (p->nr > n)
450     p->nr = n;
451 gregl 3.11 complist[listpos].nc += p->nr; /* find where this one would go */
452 gwlarson 3.27 if (hdgetbeam(hdlist[p->hd], p->bi) != NULL)
453     hdfreefrag(hdlist[p->hd], p->bi);
454 gregl 3.11 while (lastin > listpos &&
455     beamcmp(complist+lastin, complist+listpos) > 0)
456 gregl 3.1 lastin--;
457     listpos++;
458     return(1);
459     }