ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rholo3.c
Revision: 3.34
Committed: Sat Feb 22 02:07:25 2003 UTC (21 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R5
Changes since 3.33: +9 -12 lines
Log Message:
Changes and check-in for 3.5 release
Includes new source files and modifications not recorded for many years
See ray/doc/notes/ReleaseNotes for notes between 3.1 and 3.5 release

File Contents

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