ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rholo3.c
Revision: 3.40
Committed: Mon Jul 21 22:30:18 2003 UTC (20 years, 8 months ago) by schorsch
Content type: text/plain
Branch: MAIN
Changes since 3.39: +4 -4 lines
Log Message:
Eliminated copystruct() macro, which is unnecessary in ANSI.
Reduced ambiguity warnings for nested if/if/else clauses.

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rholo3.c,v 3.39 2003/07/07 17:21:51 greg Exp $";
3 #endif
4 /*
5 * Routines for tracking beam compuatations
6 */
7
8 #include <string.h>
9
10 #include "rholo.h"
11 #include "view.h"
12
13 #ifndef NFRAG2CHUNK
14 #define NFRAG2CHUNK 4096 /* number of fragments to start chunking */
15 #endif
16
17 #ifndef abs
18 #define abs(x) ((x) > 0 ? (x) : -(x))
19 #endif
20 #ifndef sgn
21 #define sgn(x) ((x) > 0 ? 1 : (x) < 0 ? -1 : 0)
22 #endif
23
24 #define rchunk(n) (((n)+(RPACKSIZ/2))/RPACKSIZ)
25
26 extern time_t time();
27
28 int chunkycmp = 0; /* clump beams together on disk */
29
30 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
35
36 int
37 beamcmp(b0, b1) /* comparison for compute order */
38 register PACKHEAD *b0, *b1;
39 {
40 BEAMI *bip0, *bip1;
41 register long c;
42 /* first check desired quantities */
43 if (chunkycmp)
44 c = rchunk(b1->nr)*(rchunk(b0->nc)+1L) -
45 rchunk(b0->nr)*(rchunk(b1->nc)+1L);
46 else
47 c = b1->nr*(b0->nc+1L) - b0->nr*(b1->nc+1L);
48 if (c > 0) return(1);
49 if (c < 0) return(-1);
50 /* 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 }
67
68
69 int
70 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 void
81 dispbeam(b, hb) /* display a holodeck beam */
82 register BEAM *b;
83 register HDBEAMI *hb;
84 {
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 else p = (PACKHEAD *)realloc((void *)p, packsiz(n));
94 CHECK(p==NULL, SYSTEM, "out of memory in dispbeam");
95 }
96 /* assign packet fields */
97 memcpy((void *)packra(p), (void *)hdbray(b), b->nrm*sizeof(RAYVAL));
98 p->nr = p->nc = b->nrm;
99 for (p->hd = 0; hdlist[p->hd] != hb->h; p->hd++)
100 if (hdlist[p->hd] == NULL)
101 error(CONSISTENCY, "unregistered holodeck in dispbeam");
102 p->bi = hb->b;
103 disp_packet(p); /* display it */
104 if (n >= 1024) { /* free ridiculous packets */
105 free((void *)p);
106 p = NULL; n = 0;
107 }
108 }
109
110
111 bundle_set(op, clist, nents) /* bundle set operation */
112 int op;
113 PACKHEAD *clist;
114 int nents;
115 {
116 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 qsort((void *)clist, nents, sizeof(PACKHEAD), beamidcmp);
124 for (i = 0; i < complen; i++) {
125 csm = (PACKHEAD *)bsearch((void *)(complist+i), (void *)clist,
126 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 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 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 }
155 /* record computed rays for uncommon beams */
156 for (csm = clist+nents; csm-- > clist; )
157 if (csm->nc < 0)
158 csm->nc = bnrays(hdlist[csm->hd], csm->bi);
159 /* complete list operations */
160 switch (op) {
161 case BS_NEW: /* new computation set */
162 listpos = 0; lastin = -1;
163 if (complen) /* free old list */
164 free((void *)complist);
165 complist = NULL;
166 if (!(complen = nents))
167 return;
168 complist = (PACKHEAD *)malloc(nents*sizeof(PACKHEAD));
169 if (complist == NULL)
170 goto memerr;
171 memcpy((void *)complist, (void *)clist, nents*sizeof(PACKHEAD));
172 break;
173 case BS_ADD: /* add to computation set */
174 case BS_MAX: /* maximum of quantities */
175 case BS_ADJ: /* adjust set quantities */
176 if (nents <= 0)
177 return;
178 sortcomplist(); /* sort updated list & new entries */
179 qsort((void *)clist, nents, sizeof(PACKHEAD), beamcmp);
180 /* what can't we satisfy? */
181 for (i = nents, csm = clist; i-- && csm->nr > csm->nc; csm++)
182 ;
183 n = csm - clist;
184 if (op != BS_ADD) { /* don't regenerate adjusted beams */
185 for (++i; i-- && csm->nr > 0; csm++)
186 ;
187 nents = csm - clist;
188 }
189 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 free((void *)complist);
199 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 return; /* already done */
207 default:
208 error(CONSISTENCY, "bundle_set called with unknown operation");
209 }
210 if (outdev == NULL || !nents) /* nothing to display? */
211 return;
212 /* load and display beams we have */
213 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 }
218 hdloadbeams(hbarr, nents, dispbeam);
219 free((void *)hbarr);
220 if (hdfragflags&FF_READ) {
221 listpos = 0;
222 lastin = -1; /* need to re-sort list */
223 }
224 return;
225 memerr:
226 error(SYSTEM, "out of memory in bundle_set");
227 }
228
229
230 double
231 beamvolume(hp, bi) /* compute approximate volume of a beam */
232 HOLO *hp;
233 int bi;
234 {
235 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 }
253 /* compute difference vector */
254 VSUM(diffv, cent[1], cent[0], -1.0);
255 for (i = 0; i < 2; i++) { /* compute volume contributions */
256 vol[i] = 0.5*DOT(crossp[i], diffv);
257 if (vol[i] < 0.) vol[i] = -vol[i];
258 }
259 return(vol[0] + vol[1]); /* return total volume */
260 }
261
262
263 ambient_list() /* compute ambient beam list */
264 {
265 int32 wtotal, minrt;
266 double frac;
267 int i;
268 register int j, k;
269
270 complen = 0;
271 for (j = 0; hdlist[j] != NULL; j++)
272 complen += nbeams(hdlist[j]);
273 complist = (PACKHEAD *)malloc(complen*sizeof(PACKHEAD));
274 CHECK(complist==NULL, SYSTEM, "out of memory in ambient_list");
275 /* compute beam weights */
276 k = 0; wtotal = 0;
277 for (j = 0; hdlist[j] != NULL; j++) {
278 frac = 512. * VLEN(hdlist[j]->wg[0]) *
279 VLEN(hdlist[j]->wg[1]) *
280 VLEN(hdlist[j]->wg[2]);
281 for (i = nbeams(hdlist[j]); i > 0; i--) {
282 complist[k].hd = j;
283 complist[k].bi = i;
284 complist[k].nr = frac*beamvolume(hdlist[j], i) + 0.5;
285 complist[k].nc = bnrays(hdlist[j], i);
286 wtotal += complist[k++].nr;
287 }
288 }
289 /* adjust sample weights */
290 if (vdef(DISKSPACE))
291 frac = 1024.*1024.*vflt(DISKSPACE) / (wtotal*sizeof(RAYVAL));
292 else
293 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 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 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 free((void *)blist.bl);
324 }
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 free((void *)complist);
334 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 /* no view vicinity */
344 myeye.rng = 0;
345 }
346
347
348 mergeclists(cdest, cl1, n1, cl2, n2) /* merge two sorted lists */
349 register PACKHEAD *cdest;
350 register PACKHEAD *cl1, *cl2;
351 int n1, n2;
352 {
353 register int cmp;
354
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 *cdest = *cl2;
361 cl2++; n2--;
362 } else {
363 *cdest = *cl1;
364 cl1++; n1--;
365 }
366 cdest++;
367 }
368 }
369
370
371 sortcomplist() /* fix our list order */
372 {
373 PACKHEAD *list2;
374 int listlen;
375 register int i;
376
377 if (complen <= 0) /* check to see if there is even a list */
378 return;
379 if (!chunkycmp) /* check to see if fragment list is full */
380 if (!hdfragOK(hdlist[0]->fd, &listlen, NULL)
381 #if NFRAG2CHUNK
382 || listlen >= NFRAG2CHUNK
383 #endif
384 ) {
385 chunkycmp++; /* use "chunky" comparison */
386 lastin = -1; /* need to re-sort list */
387 #ifdef DEBUG
388 error(WARNING, "using chunky comparison mode");
389 #endif
390 }
391 if (lastin < 0 || listpos*4 >= complen*3)
392 qsort((void *)complist, complen, sizeof(PACKHEAD), beamcmp);
393 else if (listpos) { /* else sort and merge sublist */
394 list2 = (PACKHEAD *)malloc(listpos*sizeof(PACKHEAD));
395 CHECK(list2==NULL, SYSTEM, "out of memory in sortcomplist");
396 memcpy((void *)list2,(void *)complist,listpos*sizeof(PACKHEAD));
397 qsort((void *)list2, listpos, sizeof(PACKHEAD), beamcmp);
398 mergeclists(complist, list2, listpos,
399 complist+listpos, complen-listpos);
400 free((void *)list2);
401 }
402 /* drop satisfied requests */
403 for (i = complen; i-- && complist[i].nr <= complist[i].nc; )
404 ;
405 if (i < 0) {
406 free((void *)complist);
407 complist = NULL;
408 complen = 0;
409 } else if (i < complen-1) {
410 list2 = (PACKHEAD *)realloc((void *)complist,
411 (i+1)*sizeof(PACKHEAD));
412 if (list2 != NULL)
413 complist = list2;
414 complen = i+1;
415 }
416 listpos = 0; lastin = i;
417 }
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 * up, and if we get further in our compute list than this, we re-sort the
426 * list and start again from the beginning. Since
427 * a merge sort is used, the sorting costs are minimal.
428 */
429 next_packet(p, n) /* prepare packet for computation */
430 register PACKET *p;
431 int n;
432 {
433 register int i;
434
435 if (listpos > lastin) /* time to sort the list */
436 sortcomplist();
437 if (complen <= 0)
438 return(0);
439 p->hd = complist[listpos].hd;
440 p->bi = complist[listpos].bi;
441 p->nc = complist[listpos].nc;
442 p->nr = complist[listpos].nr - p->nc;
443 if (p->nr <= 0)
444 return(0);
445 DCHECK(n < 1 | n > RPACKSIZ,
446 CONSISTENCY, "next_packet called with bad n value");
447 if (p->nr > n)
448 p->nr = n;
449 complist[listpos].nc += p->nr; /* find where this one would go */
450 if (hdgetbeam(hdlist[p->hd], p->bi) != NULL)
451 hdfreefrag(hdlist[p->hd], p->bi);
452 while (lastin > listpos &&
453 beamcmp(complist+lastin, complist+listpos) > 0)
454 lastin--;
455 listpos++;
456 return(1);
457 }