ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rhd_qtree.c
Revision: 3.11
Committed: Fri Dec 5 15:40:54 1997 UTC (26 years, 3 months ago) by gregl
Content type: text/plain
Branch: MAIN
Changes since 3.10: +82 -54 lines
Log Message:
change to better, faster angle-difference computation

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 * Quadtree driver support routines.
9 */
10
11 #include "standard.h"
12 #include "rhd_qtree.h"
13 /* quantity of leaves to free at a time */
14 #ifndef LFREEPCT
15 #define LFREEPCT 25
16 #endif
17 /* maximum allowed angle difference (deg.) */
18 #ifndef MAXANG
19 #define MAXANG 20.
20 #endif
21
22 #define MAXDIFF2 (PI*PI/180./180.* MAXANG*MAXANG )
23
24 #define abs(i) ((i) < 0 ? -(i) : (i))
25
26 RTREE qtrunk; /* our quadtree trunk */
27 double qtDepthEps = .05; /* epsilon to compare depths (z fraction) */
28 int qtMinNodesiz = 2; /* minimum node dimension (pixels) */
29 struct rleaves qtL; /* our pile of leaves */
30
31 #define TBUNDLESIZ 409 /* number of twigs in a bundle */
32
33 static RTREE **twigbundle; /* free twig blocks (NULL term.) */
34 static int nexttwig; /* next free twig */
35
36 #define ungetleaf(li) (qtL.tl=(li)) /* dangerous if used improperly */
37
38
39 static RTREE *
40 newtwig() /* allocate a twig */
41 {
42 register int bi;
43
44 if (twigbundle == NULL) { /* initialize */
45 twigbundle = (RTREE **)malloc(sizeof(RTREE *));
46 if (twigbundle == NULL)
47 goto memerr;
48 twigbundle[0] = NULL;
49 }
50 bi = nexttwig / TBUNDLESIZ;
51 if (twigbundle[bi] == NULL) { /* new block */
52 twigbundle = (RTREE **)realloc((char *)twigbundle,
53 (bi+2)*sizeof(RTREE *));
54 if (twigbundle == NULL)
55 goto memerr;
56 twigbundle[bi] = (RTREE *)calloc(TBUNDLESIZ, sizeof(RTREE));
57 if (twigbundle[bi] == NULL)
58 goto memerr;
59 twigbundle[bi+1] = NULL;
60 }
61 /* nexttwig++ % TBUNDLESIZ */
62 return(twigbundle[bi] + (nexttwig++ - bi*TBUNDLESIZ));
63 memerr:
64 error(SYSTEM, "out of memory in newtwig");
65 }
66
67
68 qtFreeTree(really) /* free allocated twigs */
69 int really;
70 {
71 register int i;
72
73 qtrunk.flgs = CH_ANY; /* chop down tree */
74 if (twigbundle == NULL)
75 return;
76 i = (TBUNDLESIZ-1+nexttwig)/TBUNDLESIZ;
77 nexttwig = 0;
78 if (!really) { /* just clear allocated blocks */
79 while (i--)
80 bzero((char *)twigbundle[i], TBUNDLESIZ*sizeof(RTREE));
81 return;
82 }
83 /* else "really" means free up memory */
84 for (i = 0; twigbundle[i] != NULL; i++)
85 free((char *)twigbundle[i]);
86 free((char *)twigbundle);
87 twigbundle = NULL;
88 }
89
90
91 static int
92 newleaf() /* allocate a leaf from our pile */
93 {
94 int li;
95
96 li = qtL.tl++;
97 if (qtL.tl >= qtL.nl) /* get next leaf in ring */
98 qtL.tl = 0;
99 if (qtL.tl == qtL.bl) /* need to shake some free */
100 qtCompost(LFREEPCT);
101 return(li);
102 }
103
104
105 #define LEAFSIZ (3*sizeof(float)+sizeof(int4)+\
106 sizeof(TMbright)+6*sizeof(BYTE))
107
108 int
109 qtAllocLeaves(n) /* allocate space for n leaves */
110 register int n;
111 {
112 unsigned nbytes;
113 register unsigned i;
114
115 qtFreeTree(0); /* make sure tree is empty */
116 if (n <= 0)
117 return(0);
118 if (qtL.nl >= n)
119 return(qtL.nl);
120 else if (qtL.nl > 0)
121 free(qtL.base);
122 /* round space up to nearest power of 2 */
123 nbytes = n*LEAFSIZ + 8;
124 for (i = 1024; nbytes > i; i <<= 1)
125 ;
126 n = (i - 8) / LEAFSIZ; /* should we make sure n is even? */
127 qtL.base = (char *)malloc(n*LEAFSIZ);
128 if (qtL.base == NULL)
129 return(0);
130 /* assign larger alignment types earlier */
131 qtL.wp = (float (*)[3])qtL.base;
132 qtL.wd = (int4 *)(qtL.wp + n);
133 qtL.brt = (TMbright *)(qtL.wd + n);
134 qtL.chr = (BYTE (*)[3])(qtL.brt + n);
135 qtL.rgb = (BYTE (*)[3])(qtL.chr + n);
136 qtL.nl = n;
137 qtL.tml = qtL.bl = qtL.tl = 0;
138 return(n);
139 }
140
141 #undef LEAFSIZ
142
143
144 qtFreeLeaves() /* free our allocated leaves and twigs */
145 {
146 qtFreeTree(1); /* free tree also */
147 if (qtL.nl <= 0)
148 return;
149 free(qtL.base);
150 qtL.base = NULL;
151 qtL.nl = 0;
152 }
153
154
155 static
156 shaketree(tp) /* shake dead leaves from tree */
157 register RTREE *tp;
158 {
159 register int i, li;
160
161 for (i = 0; i < 4; i++)
162 if (tp->flgs & BRF(i)) {
163 shaketree(tp->k[i].b);
164 if (is_stump(tp->k[i].b))
165 tp->flgs &= ~BRF(i);
166 } else if (tp->flgs & LFF(i)) {
167 li = tp->k[i].li;
168 if (qtL.bl < qtL.tl ?
169 (li < qtL.bl || li >= qtL.tl) :
170 (li < qtL.bl && li >= qtL.tl))
171 tp->flgs &= ~LFF(i);
172 }
173 }
174
175
176 int
177 qtCompost(pct) /* free up some leaves */
178 int pct;
179 {
180 int nused, nclear, nmapped;
181
182 /* figure out how many leaves to clear */
183 nclear = qtL.nl * pct / 100;
184 nused = qtL.tl - qtL.bl;
185 if (nused <= 0) nused += qtL.nl;
186 nclear -= qtL.nl - nused;
187 if (nclear <= 0)
188 return(0);
189 if (nclear >= nused) { /* clear them all */
190 qtFreeTree(0);
191 qtL.tml = qtL.bl = qtL.tl = 0;
192 return(nused);
193 }
194 /* else clear leaves from bottom */
195 nmapped = qtL.tml - qtL.bl;
196 if (nmapped < 0) nmapped += qtL.nl;
197 qtL.bl += nclear;
198 if (qtL.bl >= qtL.nl) qtL.bl -= qtL.nl;
199 if (nmapped <= nclear) qtL.tml = qtL.bl;
200 shaketree(&qtrunk);
201 return(nclear);
202 }
203
204
205 #define DCSCALE 11585.2 /* (1<<13)*sqrt(2) */
206 #define FXNEG 01
207 #define FYNEG 02
208 #define FZNEG 04
209 #define FXACT 010
210 #define FZACT 020
211 #define F1SFT 5
212 #define F2SFT 18
213 #define FMASK 0x1fff
214
215 static int4
216 encodedir(dv) /* encode a normalized direction vector */
217 FVECT dv;
218 {
219 register int4 dc = 0;
220 int cd[3], cm;
221 register int i;
222
223 for (i = 0; i < 3; i++)
224 if (dv[i] < 0.) {
225 cd[i] = dv[i] * -DCSCALE;
226 dc |= 1<<i;
227 } else
228 cd[i] = dv[i] * DCSCALE;
229 if (cd[0] <= cd[1]) {
230 dc |= FXACT | cd[0] << F1SFT;
231 cm = cd[1];
232 } else {
233 dc |= cd[1] << F1SFT;
234 cm = cd[0];
235 }
236 if (cd[2] <= cm)
237 dc |= FZACT | cd[2] << F2SFT;
238 else
239 dc |= cm << F2SFT;
240 return(dc);
241 }
242
243
244 static
245 decodedir(dv, dc) /* decode a normalized direction vector */
246 FVECT dv; /* returned */
247 register int4 dc;
248 {
249 double d1, d2, der;
250
251 d1 = ((dc>>F1SFT & FMASK)+.5)/DCSCALE;
252 d2 = ((dc>>F2SFT & FMASK)+.5)/DCSCALE;
253 der = sqrt(1. - d1*d1 - d2*d2);
254 if (dc & FXACT) {
255 dv[0] = d1;
256 if (dc & FZACT) { dv[1] = der; dv[2] = d2; }
257 else { dv[1] = d2; dv[2] = der; }
258 } else {
259 dv[1] = d1;
260 if (dc & FZACT) { dv[0] = der; dv[2] = d2; }
261 else { dv[0] = d2; dv[2] = der; }
262 }
263 if (dc & FXNEG) dv[0] = -dv[0];
264 if (dc & FYNEG) dv[1] = -dv[1];
265 if (dc & FZNEG) dv[2] = -dv[2];
266 }
267
268
269 static double
270 dir2diff(dc1, dc2) /* relative radians^2 between directions */
271 int4 dc1, dc2;
272 {
273 FVECT v1, v2;
274
275 decodedir(v1, dc1);
276 decodedir(v2, dc2);
277
278 return(2. - 2.*DOT(v1,v2));
279 }
280
281
282 static double
283 fdir2diff(dc1, v2) /* relative radians^2 between directions */
284 int4 dc1;
285 register FVECT v2;
286 {
287 FVECT v1;
288
289 decodedir(v1, dc1);
290
291 return(2. - 2.*DOT(v1,v2));
292 }
293
294
295 int
296 qtFindLeaf(x, y) /* find closest leaf to (x,y) */
297 int x, y;
298 {
299 register RTREE *tp = &qtrunk;
300 int li = -1;
301 int x0=0, y0=0, x1=odev.hres, y1=odev.vres;
302 int mx, my;
303 register int q;
304 /* check limits */
305 if (x < 0 || x >= odev.hres || y < 0 || y >= odev.vres)
306 return(-1);
307 /* find nearby leaf in our tree */
308 for ( ; ; ) {
309 for (q = 0; q < 4; q++) /* find any leaf this level */
310 if (tp->flgs & LFF(q)) {
311 li = tp->k[q].li;
312 break;
313 }
314 q = 0; /* which quadrant are we? */
315 mx = (x0 + x1) >> 1;
316 my = (y0 + y1) >> 1;
317 if (x < mx) x1 = mx;
318 else {x0 = mx; q |= 01;}
319 if (y < my) y1 = my;
320 else {y0 = my; q |= 02;}
321 if (tp->flgs & BRF(q)) { /* branch down if not a leaf */
322 tp = tp->k[q].b;
323 continue;
324 }
325 if (tp->flgs & LFF(q)) /* good shot! */
326 return(tp->k[q].li);
327 return(li); /* else return what we have */
328 }
329 }
330
331
332 static
333 addleaf(li) /* add a leaf to our tree */
334 int li;
335 {
336 register RTREE *tp = &qtrunk;
337 int x0=0, y0=0, x1=odev.hres, y1=odev.vres;
338 int lo = -1;
339 double d2;
340 int x, y, mx, my;
341 double z;
342 FVECT ip, wp, vd;
343 register int q;
344 /* compute leaf location in view */
345 VCOPY(wp, qtL.wp[li]);
346 viewloc(ip, &odev.v, wp);
347 if (ip[2] <= 0. || ip[0] < 0. || ip[0] >= 1.
348 || ip[1] < 0. || ip[1] >= 1.)
349 return(0); /* behind or outside view */
350 #ifdef DEBUG
351 if (odev.v.type == VT_PAR | odev.v.vfore > FTINY)
352 error(INTERNAL, "bad view assumption in addleaf");
353 #endif
354 for (q = 0; q < 3; q++)
355 vd[q] = (wp[q] - odev.v.vp[q])/ip[2];
356 d2 = fdir2diff(qtL.wd[li], vd);
357 if (d2 > MAXDIFF2)
358 return(0); /* leaf dir. too far off */
359 x = ip[0] * odev.hres;
360 y = ip[1] * odev.vres;
361 z = ip[2];
362 /* find the place for it */
363 for ( ; ; ) {
364 q = 0; /* which quadrant? */
365 mx = (x0 + x1) >> 1;
366 my = (y0 + y1) >> 1;
367 if (x < mx) x1 = mx;
368 else {x0 = mx; q |= 01;}
369 if (y < my) y1 = my;
370 else {y0 = my; q |= 02;}
371 if (tp->flgs & BRF(q)) { /* move to next branch */
372 tp->flgs |= CHF(q); /* not sure; guess */
373 tp = tp->k[q].b;
374 continue;
375 }
376 if (!(tp->flgs & LFF(q))) { /* found stem for leaf */
377 tp->k[q].li = li;
378 tp->flgs |= CHLFF(q);
379 break;
380 }
381 if (lo != tp->k[q].li) { /* check old leaf */
382 lo = tp->k[q].li;
383 VCOPY(wp, qtL.wp[lo]);
384 viewloc(ip, &odev.v, wp);
385 }
386 /* is node minimum size? */
387 if (y1-y0 <= qtMinNodesiz || x1-x0 <= qtMinNodesiz) {
388 if (z > (1.+qtDepthEps)*ip[2])
389 return(0); /* old one closer */
390 if (z >= (1.-qtDepthEps)*ip[2] &&
391 fdir2diff(qtL.wd[lo], vd) < d2)
392 return(0); /* old one better */
393 tp->k[q].li = li; /* else new one is */
394 tp->flgs |= CHF(q);
395 break;
396 }
397 tp->flgs &= ~LFF(q); /* else grow tree */
398 tp->flgs |= CHBRF(q);
399 tp = tp->k[q].b = newtwig();
400 q = 0; /* old leaf -> new branch */
401 mx = ip[0] * odev.hres;
402 my = ip[1] * odev.vres;
403 if (mx >= (x0 + x1) >> 1) q |= 01;
404 if (my >= (y0 + y1) >> 1) q |= 02;
405 tp->flgs = CH_ANY|LFF(q); /* all new */
406 tp->k[q].li = lo;
407 }
408 return(1); /* done */
409 }
410
411
412 dev_value(c, p, v) /* add a pixel value to our quadtree */
413 COLR c;
414 FVECT p, v;
415 {
416 register int li;
417
418 li = newleaf();
419 VCOPY(qtL.wp[li], p);
420 qtL.wd[li] = encodedir(v);
421 tmCvColrs(&qtL.brt[li], qtL.chr[li], c, 1);
422 if (!addleaf(li))
423 ungetleaf(li);
424 }
425
426
427 qtReplant() /* replant our tree using new view */
428 {
429 register int i;
430 /* anything to replant? */
431 if (qtL.bl == qtL.tl)
432 return;
433 qtFreeTree(0); /* blow the old tree away */
434 /* regrow it in new place */
435 for (i = qtL.bl; i != qtL.tl; ) {
436 addleaf(i);
437 if (++i >= qtL.nl) i = 0;
438 }
439 }
440
441
442 qtMapLeaves(redo) /* map our leaves to RGB */
443 int redo;
444 {
445 int aorg, alen, borg, blen;
446 /* recompute mapping? */
447 if (redo)
448 qtL.tml = qtL.bl;
449 /* already done? */
450 if (qtL.tml == qtL.tl)
451 return(1);
452 /* compute segments */
453 aorg = qtL.tml;
454 if (qtL.tl >= aorg) {
455 alen = qtL.tl - aorg;
456 blen = 0;
457 } else {
458 alen = qtL.nl - aorg;
459 borg = 0;
460 blen = qtL.tl;
461 }
462 /* (re)compute tone mapping? */
463 if (qtL.tml == qtL.bl) {
464 tmClearHisto();
465 tmAddHisto(qtL.brt+aorg, alen, 1);
466 if (blen > 0)
467 tmAddHisto(qtL.brt+borg, blen, 1);
468 if (tmComputeMapping(0., 0., 0.) != TM_E_OK)
469 return(0);
470 }
471 if (tmMapPixels(qtL.rgb+aorg, qtL.brt+aorg,
472 qtL.chr+aorg, alen) != TM_E_OK)
473 return(0);
474 if (blen > 0)
475 tmMapPixels(qtL.rgb+borg, qtL.brt+borg,
476 qtL.chr+borg, blen);
477 qtL.tml = qtL.tl;
478 return(1);
479 }