ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rholo3.c
Revision: 3.23
Committed: Wed Aug 12 17:56:06 1998 UTC (25 years, 7 months ago) by gwlarson
Content type: text/plain
Branch: MAIN
Changes since 3.22: +68 -47 lines
Log Message:
changed parameters to function passed to hdloadbeams()
improved speed of common member match in bundle_set()

File Contents

# User Rev Content
1 gwlarson 3.23 /* Copyright (c) 1998 Silicon Graphics, Inc. */
2    
3     #ifndef lint
4     static char SCCSid[] = "$SunId$ SGI";
5     #endif
6    
7 gregl 3.1 /*
8     * Routines for tracking beam compuatations
9     */
10    
11     #include "rholo.h"
12    
13     #define abs(x) ((x) > 0 ? (x) : -(x))
14     #define sgn(x) ((x) > 0 ? 1 : (x) < 0 ? -1 : 0)
15    
16 gregl 3.4 static PACKHEAD *complist=NULL; /* list of beams to compute */
17     static int complen=0; /* length of complist */
18     static int listpos=0; /* current list position for next_packet */
19     static int lastin= -1; /* last ordered position in list */
20 gregl 3.1
21    
22     int
23 gwlarson 3.23 beamcmp(b0, b1) /* comparison for compute order */
24 gregl 3.2 register PACKHEAD *b0, *b1;
25     {
26 gregl 3.11 return( b1->nr*(b0->nc+1) - b0->nr*(b1->nc+1) );
27 gregl 3.2 }
28    
29    
30 gregl 3.14 int
31 gwlarson 3.23 beamidcmp(b0, b1) /* comparison for beam searching */
32     register PACKHEAD *b0, *b1;
33     {
34     register int c = b0->hd - b1->hd;
35    
36     if (c) return(c);
37     return(b0->bi - b1->bi);
38     }
39    
40    
41     int
42     dispbeam(b, hb) /* display a holodeck beam */
43 gregl 3.14 register BEAM *b;
44 gwlarson 3.23 register HDBEAMI *hb;
45 gregl 3.14 {
46     static int n = 0;
47     static PACKHEAD *p = NULL;
48    
49     if (b == NULL)
50     return;
51     if (b->nrm > n) { /* (re)allocate packet holder */
52     n = b->nrm;
53     if (p == NULL) p = (PACKHEAD *)malloc(packsiz(n));
54     else p = (PACKHEAD *)realloc((char *)p, packsiz(n));
55     if (p == NULL)
56     error(SYSTEM, "out of memory in dispbeam");
57     }
58     /* assign packet fields */
59     bcopy((char *)hdbray(b), (char *)packra(p), b->nrm*sizeof(RAYVAL));
60     p->nr = p->nc = b->nrm;
61 gwlarson 3.23 for (p->hd = 0; hdlist[p->hd] != hb->h; p->hd++)
62 gregl 3.14 if (hdlist[p->hd] == NULL)
63     error(CONSISTENCY, "unregistered holodeck in dispbeam");
64 gwlarson 3.23 p->bi = hb->b;
65 gregl 3.14 disp_packet(p); /* display it */
66     }
67    
68    
69 gregl 3.2 bundle_set(op, clist, nents) /* bundle set operation */
70     int op;
71 gwlarson 3.23 PACKHEAD *clist;
72 gregl 3.2 int nents;
73     {
74 gwlarson 3.23 int oldnr, n;
75     HDBEAMI *hbarr;
76     register PACKHEAD *csm;
77     register int i;
78     /* search for common members */
79     qsort((char *)clist, nents, sizeof(PACKHEAD), beamidcmp);
80     for (csm = clist+nents; csm-- > clist; )
81     csm->nc = -1;
82     for (i = 0; i < complen; i++) {
83     csm = (PACKHEAD *)bsearch((char *)(complist+i), (char *)clist,
84     nents, sizeof(PACKHEAD), beamidcmp);
85     if (csm == NULL)
86     continue;
87     oldnr = complist[i].nr;
88     csm->nc = complist[i].nc;
89     switch (op) {
90     case BS_ADD: /* add to count */
91     complist[i].nr += csm->nr;
92     csm->nr = 0;
93     break;
94     case BS_ADJ: /* reset count */
95     complist[i].nr = csm->nr;
96     csm->nr = 0;
97     break;
98     case BS_DEL: /* delete count */
99     if (csm->nr == 0 || csm->nr >= complist[i].nr)
100     complist[i].nr = 0;
101     else
102     complist[i].nr -= csm->nr;
103     break;
104     }
105     if (complist[i].nr != oldnr)
106     lastin = -1; /* flag sort */
107 gregl 3.20 }
108 gwlarson 3.23 /* computed rays for each uncommon beams */
109     for (csm = clist+nents; csm-- > clist; )
110     if (csm->nc < 0)
111     csm->nc = bnrays(hdlist[csm->hd], csm->bi);
112 gregl 3.20 /* complete list operations */
113 gregl 3.2 switch (op) {
114     case BS_NEW: /* new computation set */
115 gwlarson 3.21 listpos = 0; lastin = -1;
116 gregl 3.20 if (complen) /* free old list */
117 gregl 3.2 free((char *)complist);
118 gregl 3.20 complist = NULL;
119     if (!(complen = nents))
120 gregl 3.2 return;
121     complist = (PACKHEAD *)malloc(nents*sizeof(PACKHEAD));
122     if (complist == NULL)
123     goto memerr;
124     bcopy((char *)clist, (char *)complist, nents*sizeof(PACKHEAD));
125     break;
126     case BS_ADD: /* add to computation set */
127 gregl 3.11 case BS_ADJ: /* adjust set quantities */
128 gregl 3.2 if (nents <= 0)
129     return;
130 gregl 3.20 sortcomplist(); /* sort updated list & new entries */
131 gregl 3.2 qsort((char *)clist, nents, sizeof(PACKHEAD), beamcmp);
132     /* what can't we satisfy? */
133 gwlarson 3.23 for (i = nents, csm = clist; i-- && csm->nr > csm->nc; csm++)
134 gregl 3.2 ;
135 gwlarson 3.23 n = csm - clist;
136 gwlarson 3.22 if (op == BS_ADJ) { /* don't regenerate adjusted beams */
137     for (i = n; i < nents && clist[i].nr > 0; i++)
138     ;
139     nents = i;
140     }
141 gregl 3.2 if (n) { /* allocate space for merged list */
142     PACKHEAD *newlist;
143     newlist = (PACKHEAD *)malloc(
144     (complen+n)*sizeof(PACKHEAD) );
145     if (newlist == NULL)
146     goto memerr;
147     /* merge lists */
148     mergeclists(newlist, clist, n, complist, complen);
149     if (complen)
150     free((char *)complist);
151     complist = newlist;
152     complen += n;
153     }
154     listpos = 0;
155     lastin = complen-1; /* list is now sorted */
156     break;
157     case BS_DEL: /* delete from computation set */
158 gregl 3.20 return; /* already done */
159 gregl 3.2 default:
160     error(CONSISTENCY, "bundle_set called with unknown operation");
161     }
162 gwlarson 3.22 if (outdev == NULL) /* nothing to display? */
163     return;
164     /* load and display beams we have */
165 gwlarson 3.23 hbarr = (HDBEAMI *)malloc(nents*sizeof(HDBEAMI));
166     for (i = nents; i--; ) {
167     hbarr[i].h = hdlist[clist[i].hd];
168     hbarr[i].b = clist[i].bi;
169 gregl 3.11 }
170 gwlarson 3.23 hdloadbeams(hbarr, nents, dispbeam);
171     free((char *)hbarr);
172 gregl 3.2 return;
173     memerr:
174     error(SYSTEM, "out of memory in bundle_set");
175     }
176    
177    
178 gregl 3.13 double
179     beamvolume(hp, bi) /* compute approximate volume of a beam */
180 gregl 3.1 HOLO *hp;
181 gregl 3.13 int bi;
182 gregl 3.1 {
183 gregl 3.13 GCOORD gc[2];
184     FVECT cp[4], edgeA, edgeB, cent[2];
185     FVECT v, crossp[2], diffv;
186     double vol[2];
187     register int i;
188     /* get grid coordinates */
189     if (!hdbcoord(gc, hp, bi))
190     error(CONSISTENCY, "bad beam index in beamvolume");
191     for (i = 0; i < 2; i++) { /* compute cell area vectors */
192     hdcell(cp, hp, gc+i);
193     VSUM(edgeA, cp[1], cp[0], -1.0);
194     VSUM(edgeB, cp[3], cp[1], -1.0);
195     fcross(crossp[i], edgeA, edgeB);
196     /* compute center */
197     cent[i][0] = 0.5*(cp[0][0] + cp[2][0]);
198     cent[i][1] = 0.5*(cp[0][1] + cp[2][1]);
199     cent[i][2] = 0.5*(cp[0][2] + cp[2][2]);
200 gregl 3.1 }
201 gregl 3.13 /* compute difference vector */
202     VSUM(diffv, cent[1], cent[0], -1.0);
203     for (i = 0; i < 2; i++) { /* compute volume contributions */
204 gregl 3.15 vol[i] = 0.5*DOT(crossp[i], diffv);
205 gregl 3.13 if (vol[i] < 0.) vol[i] = -vol[i];
206     }
207     return(vol[0] + vol[1]); /* return total volume */
208 gregl 3.1 }
209    
210    
211     init_global() /* initialize global ray computation */
212     {
213     long wtotal = 0;
214     double frac;
215 gregl 3.13 int i;
216     register int j, k;
217 gregl 3.18 /* free old list and empty queue */
218     if (complen > 0) {
219 gregl 3.3 free((char *)complist);
220 gregl 3.18 done_packets(flush_queue());
221     }
222 gregl 3.1 /* allocate beam list */
223     complen = 0;
224     for (j = 0; hdlist[j] != NULL; j++)
225     complen += nbeams(hdlist[j]);
226     complist = (PACKHEAD *)malloc(complen*sizeof(PACKHEAD));
227     if (complist == NULL)
228     error(SYSTEM, "out of memory in init_global");
229     /* compute beam weights */
230     k = 0;
231 gregl 3.13 for (j = 0; hdlist[j] != NULL; j++) {
232 gregl 3.19 frac = 512. * VLEN(hdlist[j]->wg[0]) *
233     VLEN(hdlist[j]->wg[1]) *
234     VLEN(hdlist[j]->wg[2]);
235 gregl 3.1 for (i = nbeams(hdlist[j]); i > 0; i--) {
236     complist[k].hd = j;
237     complist[k].bi = i;
238 gregl 3.13 complist[k].nr = frac*beamvolume(hdlist[j], i) + 0.5;
239 gregl 3.18 complist[k].nc = bnrays(hdlist[j], i);
240 gregl 3.1 wtotal += complist[k++].nr;
241     }
242 gregl 3.13 }
243 gregl 3.1 /* adjust weights */
244 gregl 3.12 if (vdef(DISKSPACE))
245 gregl 3.1 frac = 1024.*1024.*vflt(DISKSPACE) / (wtotal*sizeof(RAYVAL));
246 gregl 3.12 else
247     frac = 1024.*1024.*16384. / (wtotal*sizeof(RAYVAL));
248     while (k--)
249     complist[k].nr = frac * complist[k].nr;
250 gwlarson 3.21 listpos = 0; lastin = -1; /* perform initial sort */
251     sortcomplist();
252 gregl 3.1 }
253    
254    
255     mergeclists(cdest, cl1, n1, cl2, n2) /* merge two sorted lists */
256 gregl 3.16 register PACKHEAD *cdest;
257     register PACKHEAD *cl1, *cl2;
258 gregl 3.1 int n1, n2;
259     {
260 gregl 3.16 register int cmp;
261 gregl 3.1
262     while (n1 | n2) {
263     if (!n1) cmp = 1;
264     else if (!n2) cmp = -1;
265     else cmp = beamcmp(cl1, cl2);
266     if (cmp > 0) {
267     copystruct(cdest, cl2);
268     cl2++; n2--;
269     } else {
270     copystruct(cdest, cl1);
271     cl1++; n1--;
272     }
273     cdest++;
274     }
275     }
276    
277    
278     sortcomplist() /* fix our list order */
279     {
280     PACKHEAD *list2;
281 gregl 3.2 register int i;
282    
283 gregl 3.3 if (complen <= 0) /* check to see if there is even a list */
284 gregl 3.2 return;
285 gregl 3.6 if (lastin < 0 || listpos*4 >= complen*3)
286 gregl 3.1 qsort((char *)complist, complen, sizeof(PACKHEAD), beamcmp);
287     else if (listpos) { /* else sort and merge sublist */
288     list2 = (PACKHEAD *)malloc(listpos*sizeof(PACKHEAD));
289     if (list2 == NULL)
290     error(SYSTEM, "out of memory in sortcomplist");
291     bcopy((char *)complist,(char *)list2,listpos*sizeof(PACKHEAD));
292     qsort((char *)list2, listpos, sizeof(PACKHEAD), beamcmp);
293     mergeclists(complist, list2, listpos,
294     complist+listpos, complen-listpos);
295     free((char *)list2);
296     }
297 gregl 3.2 /* drop satisfied requests */
298 gregl 3.11 for (i = complen; i-- && complist[i].nr <= complist[i].nc; )
299 gregl 3.2 ;
300 gregl 3.4 if (i < 0) {
301     free((char *)complist);
302     complist = NULL;
303     complen = 0;
304     } else if (i < complen-1) {
305 gregl 3.2 list2 = (PACKHEAD *)realloc((char *)complist,
306     (i+1)*sizeof(PACKHEAD));
307     if (list2 != NULL) {
308     complist = list2;
309     complen = i+1;
310     }
311     }
312     listpos = 0; lastin = i;
313 gregl 3.1 }
314    
315    
316     /*
317     * The following routine works on the assumption that the bundle weights are
318     * more or less evenly distributed, such that computing a packet causes
319     * a given bundle to move way down in the computation order. We keep
320     * track of where the computed bundle with the highest priority would end
321     * up, and if we get further in our compute list than this, we resort the
322 gregl 3.11 * list and start again from the beginning. Since
323     * a merge sort is used, the sorting costs are minimal.
324 gregl 3.1 */
325 gregl 3.17 next_packet(p, n) /* prepare packet for computation */
326 gregl 3.1 register PACKET *p;
327 gregl 3.17 int n;
328 gregl 3.1 {
329     register int i;
330    
331 gregl 3.10 if (listpos > lastin) /* time to sort the list */
332     sortcomplist();
333 gregl 3.1 if (complen <= 0)
334     return(0);
335     p->hd = complist[listpos].hd;
336     p->bi = complist[listpos].bi;
337 gregl 3.11 p->nc = complist[listpos].nc;
338     p->nr = complist[listpos].nr - p->nc;
339 gregl 3.1 if (p->nr <= 0)
340     return(0);
341 gregl 3.17 #ifdef DEBUG
342     if (n < 1 | n > RPACKSIZ)
343     error(CONSISTENCY, "next_packet called with bad n value");
344     #endif
345     if (p->nr > n)
346     p->nr = n;
347 gregl 3.11 complist[listpos].nc += p->nr; /* find where this one would go */
348     while (lastin > listpos &&
349     beamcmp(complist+lastin, complist+listpos) > 0)
350 gregl 3.1 lastin--;
351     listpos++;
352     return(1);
353     }