ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rhd_qtree.c
(Generate patch)

Comparing ray/src/hd/rhd_qtree.c (file contents):
Revision 3.10 by gregl, Fri Dec 5 09:40:05 1997 UTC vs.
Revision 3.12 by gregl, Fri Dec 5 16:22:49 1997 UTC

# Line 16 | Line 16 | static char SCCSid[] = "$SunId$ SGI";
16   #endif
17                                  /* maximum allowed angle difference (deg.) */
18   #ifndef MAXANG
19 < #define MAXANG          20.
19 > #define MAXANG          20
20   #endif
21 + #if MAXANG>0
22 + #define MAXDIFF2        ( MAXANG*MAXANG * (PI*PI/180./180.))
23 + #endif
24  
22 #define MAXDIFF2        (long)( MAXANG*MAXANG /90./90.*(1L<<15)*(1L<<15))
23
25   #define abs(i)          ((i) < 0 ? -(i) : (i))
26  
27   RTREE   qtrunk;                 /* our quadtree trunk */
# Line 102 | Line 103 | newleaf()                      /* allocate a leaf from our pile */
103   }
104  
105  
106 < #define LEAFSIZ         (3*sizeof(float)+2*sizeof(short)+\
106 > #define LEAFSIZ         (3*sizeof(float)+sizeof(int4)+\
107                          sizeof(TMbright)+6*sizeof(BYTE))
108  
109   int
# Line 129 | Line 130 | register int   n;
130                  return(0);
131                                  /* assign larger alignment types earlier */
132          qtL.wp = (float (*)[3])qtL.base;
133 <        qtL.wd = (short (*)[2])(qtL.wp + n);
133 >        qtL.wd = (int4 *)(qtL.wp + n);
134          qtL.brt = (TMbright *)(qtL.wd + n);
135          qtL.chr = (BYTE (*)[3])(qtL.brt + n);
136          qtL.rgb = (BYTE (*)[3])(qtL.chr + n);
# Line 202 | Line 203 | int    pct;
203   }
204  
205  
206 < static
207 < encodedir(pa, dv)               /* encode a normalized direction vector */
208 < short   pa[2];
206 > #define DCSCALE         11585.2         /* (1<<13)*sqrt(2) */
207 > #define FXNEG           01
208 > #define FYNEG           02
209 > #define FZNEG           04
210 > #define F1X             010
211 > #define F2Z             020
212 > #define F1SFT           5
213 > #define F2SFT           18
214 > #define FMASK           0x1fff
215 >
216 > static int4
217 > encodedir(dv)           /* encode a normalized direction vector */
218   FVECT   dv;
219   {
220 <        pa[1] = 0;
221 <        if (dv[2] >= 1.)
222 <                pa[0] = (1L<<15)-1;
223 <        else if (dv[2] <= -1.)
224 <                pa[0] = -((1L<<15)-1);
225 <        else {
226 <                pa[0] = ((1L<<15)-1)/(PI/2.) * asin(dv[2]);
227 <                pa[1] = ((1L<<15)-1)/PI * atan2(dv[1], dv[0]);
220 >        register int4   dc = 0;
221 >        int     cd[3], cm;
222 >        register int    i;
223 >
224 >        for (i = 0; i < 3; i++)
225 >                if (dv[i] < 0.) {
226 >                        cd[i] = dv[i] * -DCSCALE;
227 >                        dc |= 1<<i;
228 >                } else
229 >                        cd[i] = dv[i] * DCSCALE;
230 >        if (cd[0] <= cd[1]) {
231 >                dc |= F1X | cd[0] << F1SFT;
232 >                cm = cd[1];
233 >        } else {
234 >                dc |= cd[1] << F1SFT;
235 >                cm = cd[0];
236          }
237 +        if (cd[2] <= cm)
238 +                dc |= F2Z | cd[2] << F2SFT;
239 +        else
240 +                dc |= cm << F2SFT;
241 +        return(dc);
242   }
243  
244  
245 < #define ALTSHFT         5
246 < #define NALT            (1<<ALTSHFT)
247 < #define azisft(alt)     azisftab[abs(alt)>>(15-ALTSHFT)]
245 > static
246 > decodedir(dv, dc)       /* decode a normalized direction vector */
247 > register FVECT  dv;     /* returned */
248 > register int4   dc;
249 > {
250 >        double  d1, d2, der;
251  
252 < static unsigned short   azisftab[NALT];
252 >        d1 = ((dc>>F1SFT & FMASK)+.5)/DCSCALE;
253 >        d2 = ((dc>>F2SFT & FMASK)+.5)/DCSCALE;
254 >        der = sqrt(1. - d1*d1 - d2*d2);
255 >        if (dc & F1X) {
256 >                dv[0] = d1;
257 >                if (dc & F2Z) { dv[1] = der; dv[2] = d2; }
258 >                else { dv[1] = d2; dv[2] = der; }
259 >        } else {
260 >                dv[1] = d1;
261 >                if (dc & F2Z) { dv[0] = der; dv[2] = d2; }
262 >                else { dv[0] = d2; dv[2] = der; }
263 >        }
264 >        if (dc & FXNEG) dv[0] = -dv[0];
265 >        if (dc & FYNEG) dv[1] = -dv[1];
266 >        if (dc & FZNEG) dv[2] = -dv[2];
267 > }
268  
269 < static
270 < azisftinit(alt)         /* initialize azimuth scale factor table */
271 < int     alt;
269 >
270 > static double
271 > dir2diff(dc1, dc2)              /* relative radians^2 between directions */
272 > int4    dc1, dc2;
273   {
274 <        register int    i;
274 >        FVECT   v1, v2;
275  
276 <        for (i = NALT; i--; )
277 <                azisftab[i] = 2.*(1L<<15) * cos(PI/2.*(i+.5)/NALT);
278 <        return(azisft(alt));
276 >        decodedir(v1, dc1);
277 >        decodedir(v2, dc2);
278 >
279 >        return(2. - 2.*DOT(v1,v2));
280   }
281  
239 #define azisf(alt)      (azisftab[0] ? azisft(alt) : azisftinit(alt)) >> 15
282  
283 < static long
284 < dir2diff(pa1, pa2)              /* relative distance^2 between directions */
285 < short   pa1[2], pa2[2];
283 > static double
284 > fdir2diff(dc1, v2)              /* relative radians^2 between directions */
285 > int4    dc1;
286 > register FVECT  v2;
287   {
288 <        long    altd2, azid2;
246 <        int     alt;
288 >        FVECT   v1;
289  
290 <        altd2 = pa1[0] - pa2[0];        /* get altitude difference^2 */
291 <        altd2 *= altd2;
292 <        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);
290 >        decodedir(v1, dc1);
291 >
292 >        return(2. - 2.*DOT(v1,v2));
293   }
294  
295  
# Line 306 | Line 337 | int    li;
337          register RTREE  *tp = &qtrunk;
338          int     x0=0, y0=0, x1=odev.hres, y1=odev.vres;
339          int     lo = -1;
340 <        long    d2;
310 <        short   dc[2];
340 >        double  d2;
341          int     x, y, mx, my;
342          double  z;
343 <        FVECT   ip, wp;
343 >        FVECT   ip, wp, vd;
344          register int    q;
345                                          /* compute leaf location in view */
346          VCOPY(wp, qtL.wp[li]);
347          viewloc(ip, &odev.v, wp);
348          if (ip[2] <= 0. || ip[0] < 0. || ip[0] >= 1.
349                          || ip[1] < 0. || ip[1] >= 1.)
350 <                return(0);                      /* behind or outside view */
350 >                return(-1);                     /* behind or outside view */
351   #ifdef DEBUG
352          if (odev.v.type == VT_PAR | odev.v.vfore > FTINY)
353                  error(INTERNAL, "bad view assumption in addleaf");
354   #endif
355          for (q = 0; q < 3; q++)
356 <                wp[q] = (wp[q] - odev.v.vp[q])/ip[2];
357 <        encodedir(dc, wp);              /* compute pixel direction */
358 <        d2 = dir2diff(dc, qtL.wd[li]);
356 >                vd[q] = (wp[q] - odev.v.vp[q])/ip[2];
357 >        d2 = fdir2diff(qtL.wd[li], vd);
358 > #ifdef MAXDIFF2
359          if (d2 > MAXDIFF2)
360                  return(0);                      /* leaf dir. too far off */
361 + #endif
362          x = ip[0] * odev.hres;
363          y = ip[1] * odev.vres;
364          z = ip[2];
# Line 360 | Line 391 | int    li;
391                          if (z > (1.+qtDepthEps)*ip[2])
392                                  return(0);              /* old one closer */
393                          if (z >= (1.-qtDepthEps)*ip[2] &&
394 <                                        dir2diff(dc, qtL.wd[lo]) < d2)
394 >                                        fdir2diff(qtL.wd[lo], vd) < d2)
395                                  return(0);              /* old one better */
396                          tp->k[q].li = li;               /* else new one is */
397                          tp->flgs |= CHF(q);
# Line 374 | Line 405 | int    li;
405                  my = ip[1] * odev.vres;
406                  if (mx >= (x0 + x1) >> 1) q |= 01;
407                  if (my >= (y0 + y1) >> 1) q |= 02;
408 +                tp->flgs = CH_ANY|LFF(q);       /* all new */
409                  tp->k[q].li = lo;
378                tp->flgs |= LFF(q)|CH_ANY;      /* all new */
410          }
411          return(1);              /* done */
412   }
# Line 389 | Line 420 | FVECT  p, v;
420  
421          li = newleaf();
422          VCOPY(qtL.wp[li], p);
423 <        encodedir(qtL.wd[li], v);
423 >        qtL.wd[li] = encodedir(v);
424          tmCvColrs(&qtL.brt[li], qtL.chr[li], c, 1);
425          if (!addleaf(li))
426                  ungetleaf(li);

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines