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

# Content
1 /* Copyright (c) 1998 Silicon Graphics, Inc. */
2
3 #ifndef lint
4 static char SCCSid[] = "$SunId$ SGI";
5 #endif
6
7 /*
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 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
21
22 int
23 beamcmp(b0, b1) /* comparison for compute order */
24 register PACKHEAD *b0, *b1;
25 {
26 return( b1->nr*(b0->nc+1) - b0->nr*(b1->nc+1) );
27 }
28
29
30 int
31 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 register BEAM *b;
44 register HDBEAMI *hb;
45 {
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 for (p->hd = 0; hdlist[p->hd] != hb->h; p->hd++)
62 if (hdlist[p->hd] == NULL)
63 error(CONSISTENCY, "unregistered holodeck in dispbeam");
64 p->bi = hb->b;
65 disp_packet(p); /* display it */
66 }
67
68
69 bundle_set(op, clist, nents) /* bundle set operation */
70 int op;
71 PACKHEAD *clist;
72 int nents;
73 {
74 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 }
108 /* 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 /* complete list operations */
113 switch (op) {
114 case BS_NEW: /* new computation set */
115 listpos = 0; lastin = -1;
116 if (complen) /* free old list */
117 free((char *)complist);
118 complist = NULL;
119 if (!(complen = nents))
120 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 case BS_ADJ: /* adjust set quantities */
128 if (nents <= 0)
129 return;
130 sortcomplist(); /* sort updated list & new entries */
131 qsort((char *)clist, nents, sizeof(PACKHEAD), beamcmp);
132 /* what can't we satisfy? */
133 for (i = nents, csm = clist; i-- && csm->nr > csm->nc; csm++)
134 ;
135 n = csm - clist;
136 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 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 return; /* already done */
159 default:
160 error(CONSISTENCY, "bundle_set called with unknown operation");
161 }
162 if (outdev == NULL) /* nothing to display? */
163 return;
164 /* load and display beams we have */
165 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 }
170 hdloadbeams(hbarr, nents, dispbeam);
171 free((char *)hbarr);
172 return;
173 memerr:
174 error(SYSTEM, "out of memory in bundle_set");
175 }
176
177
178 double
179 beamvolume(hp, bi) /* compute approximate volume of a beam */
180 HOLO *hp;
181 int bi;
182 {
183 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 }
201 /* compute difference vector */
202 VSUM(diffv, cent[1], cent[0], -1.0);
203 for (i = 0; i < 2; i++) { /* compute volume contributions */
204 vol[i] = 0.5*DOT(crossp[i], diffv);
205 if (vol[i] < 0.) vol[i] = -vol[i];
206 }
207 return(vol[0] + vol[1]); /* return total volume */
208 }
209
210
211 init_global() /* initialize global ray computation */
212 {
213 long wtotal = 0;
214 double frac;
215 int i;
216 register int j, k;
217 /* free old list and empty queue */
218 if (complen > 0) {
219 free((char *)complist);
220 done_packets(flush_queue());
221 }
222 /* 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 for (j = 0; hdlist[j] != NULL; j++) {
232 frac = 512. * VLEN(hdlist[j]->wg[0]) *
233 VLEN(hdlist[j]->wg[1]) *
234 VLEN(hdlist[j]->wg[2]);
235 for (i = nbeams(hdlist[j]); i > 0; i--) {
236 complist[k].hd = j;
237 complist[k].bi = i;
238 complist[k].nr = frac*beamvolume(hdlist[j], i) + 0.5;
239 complist[k].nc = bnrays(hdlist[j], i);
240 wtotal += complist[k++].nr;
241 }
242 }
243 /* adjust weights */
244 if (vdef(DISKSPACE))
245 frac = 1024.*1024.*vflt(DISKSPACE) / (wtotal*sizeof(RAYVAL));
246 else
247 frac = 1024.*1024.*16384. / (wtotal*sizeof(RAYVAL));
248 while (k--)
249 complist[k].nr = frac * complist[k].nr;
250 listpos = 0; lastin = -1; /* perform initial sort */
251 sortcomplist();
252 }
253
254
255 mergeclists(cdest, cl1, n1, cl2, n2) /* merge two sorted lists */
256 register PACKHEAD *cdest;
257 register PACKHEAD *cl1, *cl2;
258 int n1, n2;
259 {
260 register int cmp;
261
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 register int i;
282
283 if (complen <= 0) /* check to see if there is even a list */
284 return;
285 if (lastin < 0 || listpos*4 >= complen*3)
286 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 /* drop satisfied requests */
298 for (i = complen; i-- && complist[i].nr <= complist[i].nc; )
299 ;
300 if (i < 0) {
301 free((char *)complist);
302 complist = NULL;
303 complen = 0;
304 } else if (i < complen-1) {
305 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 }
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 * list and start again from the beginning. Since
323 * a merge sort is used, the sorting costs are minimal.
324 */
325 next_packet(p, n) /* prepare packet for computation */
326 register PACKET *p;
327 int n;
328 {
329 register int i;
330
331 if (listpos > lastin) /* time to sort the list */
332 sortcomplist();
333 if (complen <= 0)
334 return(0);
335 p->hd = complist[listpos].hd;
336 p->bi = complist[listpos].bi;
337 p->nc = complist[listpos].nc;
338 p->nr = complist[listpos].nr - p->nc;
339 if (p->nr <= 0)
340 return(0);
341 #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 complist[listpos].nc += p->nr; /* find where this one would go */
348 while (lastin > listpos &&
349 beamcmp(complist+lastin, complist+listpos) > 0)
350 lastin--;
351 listpos++;
352 return(1);
353 }