ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rholo3.c
Revision: 3.42
Committed: Tue Jun 8 19:48:30 2004 UTC (20 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R7P2, rad3R7P1, rad4R0, rad3R6, rad3R6P1, rad3R8, rad3R9
Changes since 3.41: +1 -6 lines
Log Message:
Removed redundant #include's and fixed ordering on some headers

File Contents

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