ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rholo3.c
Revision: 3.38
Committed: Mon Jun 30 14:59:12 2003 UTC (21 years, 3 months ago) by schorsch
Content type: text/plain
Branch: MAIN
Changes since 3.37: +6 -4 lines
Log Message:
Replaced most outdated BSD function calls with their posix equivalents, and cleaned up a few other platform dependencies.

File Contents

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