ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rholo3.c
Revision: 3.12
Committed: Mon Dec 1 16:35:02 1997 UTC (26 years, 10 months ago) by gregl
Content type: text/plain
Branch: MAIN
Changes since 3.11: +17 -15 lines
Log Message:
minor bug fixes and enhancements

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