ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rholo3.c
Revision: 3.43
Committed: Thu Sep 30 14:10:19 2010 UTC (13 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2P2, rad5R0, rad5R1, rad4R2, rad4R1, rad4R2P1
Changes since 3.42: +7 -2 lines
Log Message:
Increased maximum "ambient" holodeck size from 2 GBytes to 10 GBytes

File Contents

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