ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rholo3.c
Revision: 3.15
Committed: Fri Dec 19 09:55:36 1997 UTC (26 years, 3 months ago) by gregl
Content type: text/plain
Branch: MAIN
Changes since 3.14: +2 -5 lines
Log Message:
simplified beam volume calculation slightly and fixed bug in init_global()

File Contents

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