ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rhd_qtree.c
Revision: 3.10
Committed: Fri Dec 5 09:40:05 1997 UTC (26 years, 10 months ago) by gregl
Content type: text/plain
Branch: MAIN
Changes since 3.9: +103 -14 lines
Log Message:
added code for determining closest ray direction in addleaf()

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.10 #define MAXDIFF2 (long)( MAXANG*MAXANG /90./90.*(1L<<15)*(1L<<15))
23    
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.10 #define LEAFSIZ (3*sizeof(float)+2*sizeof(short)+\
106     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.10 qtL.wd = (short (*)[2])(qtL.wp + n);
133     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.10 static
206     encodedir(pa, dv) /* encode a normalized direction vector */
207     short pa[2];
208     FVECT dv;
209     {
210     pa[1] = 0;
211     if (dv[2] >= 1.)
212     pa[0] = (1L<<15)-1;
213     else if (dv[2] <= -1.)
214     pa[0] = -((1L<<15)-1);
215     else {
216     pa[0] = ((1L<<15)-1)/(PI/2.) * asin(dv[2]);
217     pa[1] = ((1L<<15)-1)/PI * atan2(dv[1], dv[0]);
218     }
219     }
220    
221    
222     #define ALTSHFT 5
223     #define NALT (1<<ALTSHFT)
224     #define azisft(alt) azisftab[abs(alt)>>(15-ALTSHFT)]
225    
226     static unsigned short azisftab[NALT];
227    
228     static
229     azisftinit(alt) /* initialize azimuth scale factor table */
230     int alt;
231     {
232     register int i;
233    
234     for (i = NALT; i--; )
235     azisftab[i] = 2.*(1L<<15) * cos(PI/2.*(i+.5)/NALT);
236     return(azisft(alt));
237     }
238    
239     #define azisf(alt) (azisftab[0] ? azisft(alt) : azisftinit(alt)) >> 15
240    
241     static long
242     dir2diff(pa1, pa2) /* relative distance^2 between directions */
243     short pa1[2], pa2[2];
244     {
245     long altd2, azid2;
246     int alt;
247    
248     altd2 = pa1[0] - pa2[0]; /* get altitude difference^2 */
249     altd2 *= altd2;
250     if (altd2 > MAXDIFF2)
251     return(altd2); /* too large already */
252     azid2 = pa1[1] - pa2[1]; /* get adjusted azimuth difference^2 */
253     if (azid2 < 0) azid2 = -azid2;
254     if (azid2 >= 1L<<15) { /* wrap sphere */
255     azid2 -= 1L<<16;
256     if (azid2 < 0) azid2 = -azid2;
257     }
258     alt = (pa1[0]+pa2[0])/2;
259     azid2 = azid2*azisf(alt); /* evaluation order is important */
260     azid2 *= azid2;
261     return(altd2 + azid2);
262     }
263    
264    
265 gregl 3.5 int
266 gregl 3.3 qtFindLeaf(x, y) /* find closest leaf to (x,y) */
267     int x, y;
268     {
269     register RTREE *tp = &qtrunk;
270 gregl 3.5 int li = -1;
271 gregl 3.3 int x0=0, y0=0, x1=odev.hres, y1=odev.vres;
272     int mx, my;
273     register int q;
274     /* check limits */
275     if (x < 0 || x >= odev.hres || y < 0 || y >= odev.vres)
276 gregl 3.5 return(-1);
277 gregl 3.3 /* find nearby leaf in our tree */
278     for ( ; ; ) {
279     for (q = 0; q < 4; q++) /* find any leaf this level */
280 gregl 3.5 if (tp->flgs & LFF(q)) {
281     li = tp->k[q].li;
282 gregl 3.3 break;
283     }
284     q = 0; /* which quadrant are we? */
285     mx = (x0 + x1) >> 1;
286     my = (y0 + y1) >> 1;
287     if (x < mx) x1 = mx;
288     else {x0 = mx; q |= 01;}
289     if (y < my) y1 = my;
290     else {y0 = my; q |= 02;}
291     if (tp->flgs & BRF(q)) { /* branch down if not a leaf */
292     tp = tp->k[q].b;
293     continue;
294     }
295 gregl 3.5 if (tp->flgs & LFF(q)) /* good shot! */
296     return(tp->k[q].li);
297     return(li); /* else return what we have */
298 gregl 3.3 }
299     }
300    
301    
302 gregl 3.1 static
303 gregl 3.5 addleaf(li) /* add a leaf to our tree */
304     int li;
305 gregl 3.1 {
306     register RTREE *tp = &qtrunk;
307     int x0=0, y0=0, x1=odev.hres, y1=odev.vres;
308 gregl 3.5 int lo = -1;
309 gregl 3.10 long d2;
310     short dc[2];
311 gregl 3.1 int x, y, mx, my;
312     double z;
313     FVECT ip, wp;
314     register int q;
315 gregl 3.10 /* compute leaf location in view */
316 gregl 3.5 VCOPY(wp, qtL.wp[li]);
317 gregl 3.1 viewloc(ip, &odev.v, wp);
318     if (ip[2] <= 0. || ip[0] < 0. || ip[0] >= 1.
319     || ip[1] < 0. || ip[1] >= 1.)
320 gregl 3.10 return(0); /* behind or outside view */
321     #ifdef DEBUG
322     if (odev.v.type == VT_PAR | odev.v.vfore > FTINY)
323     error(INTERNAL, "bad view assumption in addleaf");
324     #endif
325     for (q = 0; q < 3; q++)
326     wp[q] = (wp[q] - odev.v.vp[q])/ip[2];
327     encodedir(dc, wp); /* compute pixel direction */
328     d2 = dir2diff(dc, qtL.wd[li]);
329     if (d2 > MAXDIFF2)
330     return(0); /* leaf dir. too far off */
331 gregl 3.1 x = ip[0] * odev.hres;
332     y = ip[1] * odev.vres;
333     z = ip[2];
334     /* find the place for it */
335     for ( ; ; ) {
336     q = 0; /* which quadrant? */
337     mx = (x0 + x1) >> 1;
338     my = (y0 + y1) >> 1;
339     if (x < mx) x1 = mx;
340     else {x0 = mx; q |= 01;}
341     if (y < my) y1 = my;
342     else {y0 = my; q |= 02;}
343     if (tp->flgs & BRF(q)) { /* move to next branch */
344     tp->flgs |= CHF(q); /* not sure; guess */
345     tp = tp->k[q].b;
346     continue;
347     }
348 gregl 3.5 if (!(tp->flgs & LFF(q))) { /* found stem for leaf */
349     tp->k[q].li = li;
350     tp->flgs |= CHLFF(q);
351 gregl 3.1 break;
352     }
353 gregl 3.10 if (lo != tp->k[q].li) { /* check old leaf */
354 gregl 3.5 lo = tp->k[q].li;
355     VCOPY(wp, qtL.wp[lo]);
356 gregl 3.1 viewloc(ip, &odev.v, wp);
357     }
358     /* is node minimum size? */
359 gregl 3.10 if (y1-y0 <= qtMinNodesiz || x1-x0 <= qtMinNodesiz) {
360     if (z > (1.+qtDepthEps)*ip[2])
361     return(0); /* old one closer */
362     if (z >= (1.-qtDepthEps)*ip[2] &&
363     dir2diff(dc, qtL.wd[lo]) < d2)
364     return(0); /* old one better */
365     tp->k[q].li = li; /* else new one is */
366 gregl 3.1 tp->flgs |= CHF(q);
367     break;
368     }
369 gregl 3.5 tp->flgs &= ~LFF(q); /* else grow tree */
370     tp->flgs |= CHBRF(q);
371 gregl 3.1 tp = tp->k[q].b = newtwig();
372     q = 0; /* old leaf -> new branch */
373     mx = ip[0] * odev.hres;
374     my = ip[1] * odev.vres;
375     if (mx >= (x0 + x1) >> 1) q |= 01;
376     if (my >= (y0 + y1) >> 1) q |= 02;
377 gregl 3.5 tp->k[q].li = lo;
378     tp->flgs |= LFF(q)|CH_ANY; /* all new */
379 gregl 3.1 }
380 gregl 3.10 return(1); /* done */
381 gregl 3.1 }
382    
383    
384 gregl 3.10 dev_value(c, p, v) /* add a pixel value to our quadtree */
385 gregl 3.1 COLR c;
386 gregl 3.10 FVECT p, v;
387 gregl 3.1 {
388 gregl 3.5 register int li;
389 gregl 3.1
390 gregl 3.5 li = newleaf();
391     VCOPY(qtL.wp[li], p);
392 gregl 3.10 encodedir(qtL.wd[li], v);
393 gregl 3.5 tmCvColrs(&qtL.brt[li], qtL.chr[li], c, 1);
394 gregl 3.10 if (!addleaf(li))
395     ungetleaf(li);
396 gregl 3.1 }
397    
398    
399     qtReplant() /* replant our tree using new view */
400     {
401     register int i;
402 gregl 3.5 /* anything to replant? */
403     if (qtL.bl == qtL.tl)
404 gregl 3.1 return;
405 gregl 3.5 qtFreeTree(0); /* blow the old tree away */
406     /* regrow it in new place */
407     for (i = qtL.bl; i != qtL.tl; ) {
408     addleaf(i);
409     if (++i >= qtL.nl) i = 0;
410 gregl 3.1 }
411     }
412    
413    
414 gregl 3.5 qtMapLeaves(redo) /* map our leaves to RGB */
415     int redo;
416     {
417     int aorg, alen, borg, blen;
418 gregl 3.6 /* recompute mapping? */
419     if (redo)
420     qtL.tml = qtL.bl;
421 gregl 3.5 /* already done? */
422     if (qtL.tml == qtL.tl)
423     return(1);
424     /* compute segments */
425     aorg = qtL.tml;
426     if (qtL.tl >= aorg) {
427     alen = qtL.tl - aorg;
428     blen = 0;
429     } else {
430     alen = qtL.nl - aorg;
431     borg = 0;
432     blen = qtL.tl;
433     }
434     /* (re)compute tone mapping? */
435     if (qtL.tml == qtL.bl) {
436     tmClearHisto();
437     tmAddHisto(qtL.brt+aorg, alen, 1);
438     if (blen > 0)
439     tmAddHisto(qtL.brt+borg, blen, 1);
440     if (tmComputeMapping(0., 0., 0.) != TM_E_OK)
441     return(0);
442     }
443     if (tmMapPixels(qtL.rgb+aorg, qtL.brt+aorg,
444     qtL.chr+aorg, alen) != TM_E_OK)
445     return(0);
446     if (blen > 0)
447     tmMapPixels(qtL.rgb+borg, qtL.brt+borg,
448     qtL.chr+borg, blen);
449     qtL.tml = qtL.tl;
450     return(1);
451 gregl 3.1 }