ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rholo3.c
Revision: 3.18
Committed: Mon Jan 5 16:46:34 1998 UTC (26 years, 9 months ago) by gregl
Content type: text/plain
Branch: MAIN
Changes since 3.17: +6 -3 lines
Log Message:
fixed dire bug in init_global() that caused it to recompute rays

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 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 and empty queue */
216 if (complen > 0) {
217 free((char *)complist);
218 done_packets(flush_queue());
219 }
220 /* allocate beam list */
221 complen = 0;
222 for (j = 0; hdlist[j] != NULL; j++)
223 complen += nbeams(hdlist[j]);
224 complist = (PACKHEAD *)malloc(complen*sizeof(PACKHEAD));
225 if (complist == NULL)
226 error(SYSTEM, "out of memory in init_global");
227 /* compute beam weights */
228 k = 0;
229 for (j = 0; hdlist[j] != NULL; j++) {
230 frac = 512. * hdlist[j]->wg[0] *
231 hdlist[j]->wg[1] * hdlist[j]->wg[2];
232 if (frac < 0.) frac = -frac;
233 for (i = nbeams(hdlist[j]); i > 0; i--) {
234 complist[k].hd = j;
235 complist[k].bi = i;
236 complist[k].nr = frac*beamvolume(hdlist[j], i) + 0.5;
237 complist[k].nc = bnrays(hdlist[j], i);
238 wtotal += complist[k++].nr;
239 }
240 }
241 /* adjust weights */
242 if (vdef(DISKSPACE))
243 frac = 1024.*1024.*vflt(DISKSPACE) / (wtotal*sizeof(RAYVAL));
244 else
245 frac = 1024.*1024.*16384. / (wtotal*sizeof(RAYVAL));
246 while (k--)
247 complist[k].nr = frac * complist[k].nr;
248 listpos = 0; lastin = -1; /* flag initial sort */
249 }
250
251
252 mergeclists(cdest, cl1, n1, cl2, n2) /* merge two sorted lists */
253 register PACKHEAD *cdest;
254 register PACKHEAD *cl1, *cl2;
255 int n1, n2;
256 {
257 register int cmp;
258
259 while (n1 | n2) {
260 if (!n1) cmp = 1;
261 else if (!n2) cmp = -1;
262 else cmp = beamcmp(cl1, cl2);
263 if (cmp > 0) {
264 copystruct(cdest, cl2);
265 cl2++; n2--;
266 } else {
267 copystruct(cdest, cl1);
268 cl1++; n1--;
269 }
270 cdest++;
271 }
272 }
273
274
275 sortcomplist() /* fix our list order */
276 {
277 PACKHEAD *list2;
278 register int i;
279
280 if (complen <= 0) /* check to see if there is even a list */
281 return;
282 if (lastin < 0 || listpos*4 >= complen*3)
283 qsort((char *)complist, complen, sizeof(PACKHEAD), beamcmp);
284 else if (listpos) { /* else sort and merge sublist */
285 list2 = (PACKHEAD *)malloc(listpos*sizeof(PACKHEAD));
286 if (list2 == NULL)
287 error(SYSTEM, "out of memory in sortcomplist");
288 bcopy((char *)complist,(char *)list2,listpos*sizeof(PACKHEAD));
289 qsort((char *)list2, listpos, sizeof(PACKHEAD), beamcmp);
290 mergeclists(complist, list2, listpos,
291 complist+listpos, complen-listpos);
292 free((char *)list2);
293 }
294 /* drop satisfied requests */
295 for (i = complen; i-- && complist[i].nr <= complist[i].nc; )
296 ;
297 if (i < 0) {
298 free((char *)complist);
299 complist = NULL;
300 complen = 0;
301 } else if (i < complen-1) {
302 list2 = (PACKHEAD *)realloc((char *)complist,
303 (i+1)*sizeof(PACKHEAD));
304 if (list2 != NULL) {
305 complist = list2;
306 complen = i+1;
307 }
308 }
309 listpos = 0; lastin = i;
310 }
311
312
313 /*
314 * The following routine works on the assumption that the bundle weights are
315 * more or less evenly distributed, such that computing a packet causes
316 * a given bundle to move way down in the computation order. We keep
317 * track of where the computed bundle with the highest priority would end
318 * up, and if we get further in our compute list than this, we resort the
319 * list and start again from the beginning. Since
320 * a merge sort is used, the sorting costs are minimal.
321 */
322 next_packet(p, n) /* prepare packet for computation */
323 register PACKET *p;
324 int n;
325 {
326 register int i;
327
328 if (listpos > lastin) /* time to sort the list */
329 sortcomplist();
330 if (complen <= 0)
331 return(0);
332 p->hd = complist[listpos].hd;
333 p->bi = complist[listpos].bi;
334 p->nc = complist[listpos].nc;
335 p->nr = complist[listpos].nr - p->nc;
336 if (p->nr <= 0)
337 return(0);
338 #ifdef DEBUG
339 if (n < 1 | n > RPACKSIZ)
340 error(CONSISTENCY, "next_packet called with bad n value");
341 #endif
342 if (p->nr > n)
343 p->nr = n;
344 complist[listpos].nc += p->nr; /* find where this one would go */
345 while (lastin > listpos &&
346 beamcmp(complist+lastin, complist+listpos) > 0)
347 lastin--;
348 listpos++;
349 return(1);
350 }