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

# User Rev Content
1 gregl 3.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 gregl 3.6 /* quantity of leaves to free at a time */
14     #ifndef LFREEPCT
15     #define LFREEPCT 25
16     #endif
17 gregl 3.10 /* maximum allowed angle difference (deg.) */
18     #ifndef MAXANG
19     #define MAXANG 20.
20     #endif
21 gregl 3.1
22 gregl 3.11 #define MAXDIFF2 (PI*PI/180./180.* MAXANG*MAXANG )
23 gregl 3.10
24     #define abs(i) ((i) < 0 ? -(i) : (i))
25    
26 gregl 3.2 RTREE qtrunk; /* our quadtree trunk */
27 gregl 3.10 double qtDepthEps = .05; /* epsilon to compare depths (z fraction) */
28 gregl 3.2 int qtMinNodesiz = 2; /* minimum node dimension (pixels) */
29 gregl 3.5 struct rleaves qtL; /* our pile of leaves */
30 gregl 3.2
31 gregl 3.1 #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 gregl 3.10 #define ungetleaf(li) (qtL.tl=(li)) /* dangerous if used improperly */
37 gregl 3.1
38 gregl 3.10
39 gregl 3.1 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 gregl 3.3 qtFreeTree(really) /* free allocated twigs */
69 gregl 3.1 int really;
70     {
71     register int i;
72    
73 gregl 3.7 qtrunk.flgs = CH_ANY; /* chop down tree */
74 gregl 3.1 if (twigbundle == NULL)
75     return;
76 gregl 3.7 i = (TBUNDLESIZ-1+nexttwig)/TBUNDLESIZ;
77     nexttwig = 0;
78 gregl 3.1 if (!really) { /* just clear allocated blocks */
79 gregl 3.7 while (i--)
80 gregl 3.1 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 gregl 3.5 static int
92 gregl 3.1 newleaf() /* allocate a leaf from our pile */
93     {
94 gregl 3.5 int li;
95 gregl 3.4
96 gregl 3.5 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 gregl 3.1 qtCompost(LFREEPCT);
101 gregl 3.5 return(li);
102 gregl 3.1 }
103    
104    
105 gregl 3.11 #define LEAFSIZ (3*sizeof(float)+sizeof(int4)+\
106 gregl 3.10 sizeof(TMbright)+6*sizeof(BYTE))
107 gregl 3.5
108 gregl 3.1 int
109     qtAllocLeaves(n) /* allocate space for n leaves */
110 gregl 3.5 register int n;
111 gregl 3.1 {
112     unsigned nbytes;
113     register unsigned i;
114    
115 gregl 3.3 qtFreeTree(0); /* make sure tree is empty */
116 gregl 3.1 if (n <= 0)
117     return(0);
118 gregl 3.5 if (qtL.nl >= n)
119     return(qtL.nl);
120     else if (qtL.nl > 0)
121     free(qtL.base);
122 gregl 3.1 /* round space up to nearest power of 2 */
123 gregl 3.5 nbytes = n*LEAFSIZ + 8;
124 gregl 3.1 for (i = 1024; nbytes > i; i <<= 1)
125     ;
126 gregl 3.5 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 gregl 3.11 qtL.wd = (int4 *)(qtL.wp + n);
133 gregl 3.10 qtL.brt = (TMbright *)(qtL.wd + n);
134 gregl 3.5 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 gregl 3.1 }
140    
141 gregl 3.5 #undef LEAFSIZ
142 gregl 3.1
143 gregl 3.5
144 gregl 3.1 qtFreeLeaves() /* free our allocated leaves and twigs */
145     {
146 gregl 3.3 qtFreeTree(1); /* free tree also */
147 gregl 3.5 if (qtL.nl <= 0)
148 gregl 3.1 return;
149 gregl 3.5 free(qtL.base);
150     qtL.base = NULL;
151     qtL.nl = 0;
152 gregl 3.1 }
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 gregl 3.5 if (tp->flgs & BRF(i)) {
163 gregl 3.2 shaketree(tp->k[i].b);
164 gregl 3.5 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 gregl 3.1 }
173     }
174    
175    
176     int
177     qtCompost(pct) /* free up some leaves */
178     int pct;
179     {
180 gregl 3.5 int nused, nclear, nmapped;
181 gregl 3.4
182 gregl 3.1 /* figure out how many leaves to clear */
183 gregl 3.5 nclear = qtL.nl * pct / 100;
184     nused = qtL.tl - qtL.bl;
185     if (nused <= 0) nused += qtL.nl;
186     nclear -= qtL.nl - nused;
187 gregl 3.1 if (nclear <= 0)
188     return(0);
189     if (nclear >= nused) { /* clear them all */
190 gregl 3.3 qtFreeTree(0);
191 gregl 3.5 qtL.tml = qtL.bl = qtL.tl = 0;
192 gregl 3.1 return(nused);
193     }
194     /* else clear leaves from bottom */
195 gregl 3.5 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 gregl 3.1 shaketree(&qtrunk);
201     return(nclear);
202     }
203    
204    
205 gregl 3.11 #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 gregl 3.10 FVECT dv;
218     {
219 gregl 3.11 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 gregl 3.10 }
236 gregl 3.11 if (cd[2] <= cm)
237     dc |= FZACT | cd[2] << F2SFT;
238     else
239     dc |= cm << F2SFT;
240     return(dc);
241 gregl 3.10 }
242    
243    
244 gregl 3.11 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 gregl 3.10
251 gregl 3.11 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 gregl 3.10
268 gregl 3.11
269     static double
270     dir2diff(dc1, dc2) /* relative radians^2 between directions */
271     int4 dc1, dc2;
272 gregl 3.10 {
273 gregl 3.11 FVECT v1, v2;
274 gregl 3.10
275 gregl 3.11 decodedir(v1, dc1);
276     decodedir(v2, dc2);
277    
278     return(2. - 2.*DOT(v1,v2));
279 gregl 3.10 }
280    
281    
282 gregl 3.11 static double
283     fdir2diff(dc1, v2) /* relative radians^2 between directions */
284     int4 dc1;
285     register FVECT v2;
286 gregl 3.10 {
287 gregl 3.11 FVECT v1;
288 gregl 3.10
289 gregl 3.11 decodedir(v1, dc1);
290    
291     return(2. - 2.*DOT(v1,v2));
292 gregl 3.10 }
293    
294    
295 gregl 3.5 int
296 gregl 3.3 qtFindLeaf(x, y) /* find closest leaf to (x,y) */
297     int x, y;
298     {
299     register RTREE *tp = &qtrunk;
300 gregl 3.5 int li = -1;
301 gregl 3.3 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 gregl 3.5 return(-1);
307 gregl 3.3 /* find nearby leaf in our tree */
308     for ( ; ; ) {
309     for (q = 0; q < 4; q++) /* find any leaf this level */
310 gregl 3.5 if (tp->flgs & LFF(q)) {
311     li = tp->k[q].li;
312 gregl 3.3 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 gregl 3.5 if (tp->flgs & LFF(q)) /* good shot! */
326     return(tp->k[q].li);
327     return(li); /* else return what we have */
328 gregl 3.3 }
329     }
330    
331    
332 gregl 3.1 static
333 gregl 3.5 addleaf(li) /* add a leaf to our tree */
334     int li;
335 gregl 3.1 {
336     register RTREE *tp = &qtrunk;
337     int x0=0, y0=0, x1=odev.hres, y1=odev.vres;
338 gregl 3.5 int lo = -1;
339 gregl 3.11 double d2;
340 gregl 3.1 int x, y, mx, my;
341     double z;
342 gregl 3.11 FVECT ip, wp, vd;
343 gregl 3.1 register int q;
344 gregl 3.10 /* compute leaf location in view */
345 gregl 3.5 VCOPY(wp, qtL.wp[li]);
346 gregl 3.1 viewloc(ip, &odev.v, wp);
347     if (ip[2] <= 0. || ip[0] < 0. || ip[0] >= 1.
348     || ip[1] < 0. || ip[1] >= 1.)
349 gregl 3.10 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 gregl 3.11 vd[q] = (wp[q] - odev.v.vp[q])/ip[2];
356     d2 = fdir2diff(qtL.wd[li], vd);
357 gregl 3.10 if (d2 > MAXDIFF2)
358     return(0); /* leaf dir. too far off */
359 gregl 3.1 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 gregl 3.5 if (!(tp->flgs & LFF(q))) { /* found stem for leaf */
377     tp->k[q].li = li;
378     tp->flgs |= CHLFF(q);
379 gregl 3.1 break;
380     }
381 gregl 3.10 if (lo != tp->k[q].li) { /* check old leaf */
382 gregl 3.5 lo = tp->k[q].li;
383     VCOPY(wp, qtL.wp[lo]);
384 gregl 3.1 viewloc(ip, &odev.v, wp);
385     }
386     /* is node minimum size? */
387 gregl 3.10 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 gregl 3.11 fdir2diff(qtL.wd[lo], vd) < d2)
392 gregl 3.10 return(0); /* old one better */
393     tp->k[q].li = li; /* else new one is */
394 gregl 3.1 tp->flgs |= CHF(q);
395     break;
396     }
397 gregl 3.5 tp->flgs &= ~LFF(q); /* else grow tree */
398     tp->flgs |= CHBRF(q);
399 gregl 3.1 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 gregl 3.11 tp->flgs = CH_ANY|LFF(q); /* all new */
406 gregl 3.5 tp->k[q].li = lo;
407 gregl 3.1 }
408 gregl 3.10 return(1); /* done */
409 gregl 3.1 }
410    
411    
412 gregl 3.10 dev_value(c, p, v) /* add a pixel value to our quadtree */
413 gregl 3.1 COLR c;
414 gregl 3.10 FVECT p, v;
415 gregl 3.1 {
416 gregl 3.5 register int li;
417 gregl 3.1
418 gregl 3.5 li = newleaf();
419     VCOPY(qtL.wp[li], p);
420 gregl 3.11 qtL.wd[li] = encodedir(v);
421 gregl 3.5 tmCvColrs(&qtL.brt[li], qtL.chr[li], c, 1);
422 gregl 3.10 if (!addleaf(li))
423     ungetleaf(li);
424 gregl 3.1 }
425    
426    
427     qtReplant() /* replant our tree using new view */
428     {
429     register int i;
430 gregl 3.5 /* anything to replant? */
431     if (qtL.bl == qtL.tl)
432 gregl 3.1 return;
433 gregl 3.5 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 gregl 3.1 }
439     }
440    
441    
442 gregl 3.5 qtMapLeaves(redo) /* map our leaves to RGB */
443     int redo;
444     {
445     int aorg, alen, borg, blen;
446 gregl 3.6 /* recompute mapping? */
447     if (redo)
448     qtL.tml = qtL.bl;
449 gregl 3.5 /* 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 gregl 3.1 }