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

Comparing ray/src/hd/sm_geom.c (file contents):
Revision 3.4 by gwlarson, Fri Sep 11 11:52:25 1998 UTC vs.
Revision 3.9 by gwlarson, Mon Dec 28 18:07:35 1998 UTC

# Line 27 | Line 27 | FVECT v1,v2;
27   {
28     return(EQUAL(v1[0],v2[0]) && EQUAL(v1[1],v2[1])&& EQUAL(v1[2],v2[2]));
29   }
30 + #if 0
31 + extern FVECT Norm[500];
32 + extern int Ncnt;
33 + #endif
34  
31 –
35   int
36   convex_angle(v0,v1,v2)
37   FVECT v0,v1,v2;
38   {
39 <    FVECT cp01,cp12,cp;
40 <    
39 >    FVECT cp,cp01,cp12,v10,v02;
40 >    double dp;
41 >
42      /* test sign of (v0Xv1)X(v1Xv2). v1 */
43      VCROSS(cp01,v0,v1);
44      VCROSS(cp12,v1,v2);
45      VCROSS(cp,cp01,cp12);
46 <    if(DOT(cp,v1) < 0.0)
47 <       return(FALSE);
46 >        
47 >    dp = DOT(cp,v1);
48 > #if 0
49 >    VCOPY(Norm[Ncnt++],cp01);
50 >    VCOPY(Norm[Ncnt++],cp12);
51 >    VCOPY(Norm[Ncnt++],cp);
52 > #endif
53 >    if(ZERO(dp) || dp < 0.0)
54 >      return(FALSE);
55      return(TRUE);
56   }
57  
58   /* calculates the normal of a face contour using Newell's formula. e
59  
60 <               a =  SUMi (yi - yi+1)(zi + zi+1)
60 >               a =  SUMi (yi - yi+1)(zi + zi+1);
61                 b =  SUMi (zi - zi+1)(xi + xi+1)
62                 c =  SUMi (xi - xi+1)(yi + yi+1)
63   */
# Line 63 | Line 74 | int norm;
74    
75    n[1] = (v0[2] - v1[2]) * (v0[0] + v1[0]) +
76             (v1[2] - v2[2]) * (v1[0] + v2[0]) +
77 <           (v2[2] - v0[2]) * (v2[0] + v0[0]);
67 <
77 >            (v2[2] - v0[2]) * (v2[0] + v0[0]);
78    
79    n[2] = (v0[1] + v1[1]) * (v0[0] - v1[0]) +
80           (v1[1] + v2[1]) * (v1[0] - v2[0]) +
# Line 72 | Line 82 | int norm;
82  
83    if(!norm)
84       return(0);
75 –
85    
86    mag = normalize(n);
87  
# Line 80 | Line 89 | int norm;
89   }
90  
91  
92 < tri_plane_equation(v0,v1,v2,n,nd,norm)
93 <   FVECT v0,v1,v2,n;
94 <   double *nd;
92 > tri_plane_equation(v0,v1,v2,peqptr,norm)
93 >   FVECT v0,v1,v2;
94 >   FPEQ *peqptr;
95     int norm;
96   {
97 <    tri_normal(v0,v1,v2,n,norm);
97 >    tri_normal(v0,v1,v2,FP_N(*peqptr),norm);
98  
99 <    *nd = -(DOT(n,v0));
99 >    FP_D(*peqptr) = -(DOT(FP_N(*peqptr),v0));
100   }
101  
102   /* From quad_edge-code */
# Line 110 | Line 119 | FVECT p0,p1;
119   }
120  
121  
122 <
122 > double
123   point_on_sphere(ps,p,c)
124   FVECT ps,p,c;
125   {
126 +  double d;
127      VSUB(ps,p,c);    
128 <    normalize(ps);
128 >    d= normalize(ps);
129 >    return(d);
130   }
131  
132  
# Line 125 | Line 136 | FVECT ps,p,c;
136    * t, if parallel, returns t=FHUGE
137    */
138   int
139 < intersect_vector_plane(v,plane_n,plane_d,tptr,r)
140 <   FVECT v,plane_n;
141 <   double plane_d;
139 > intersect_vector_plane(v,peq,tptr,r)
140 >   FVECT v;
141 >   FPEQ peq;
142     double *tptr;
143     FVECT r;
144   {
# Line 141 | Line 152 | intersect_vector_plane(v,plane_n,plane_d,tptr,r)
152      /* line is  l = p1 + (p2-p1)t, p1=origin */
153  
154      /* Solve for t: */
155 <  d = -(DOT(plane_n,v));
155 >  d = -(DOT(FP_N(peq),v));
156    if(ZERO(d))
157    {
158        t = FHUGE;
# Line 149 | Line 160 | intersect_vector_plane(v,plane_n,plane_d,tptr,r)
160    }
161    else
162    {
163 <      t =  plane_d/d;
163 >      t =  FP_D(peq)/d;
164        if(t < 0 )
165           hit = 0;
166        else
# Line 167 | Line 178 | intersect_vector_plane(v,plane_n,plane_d,tptr,r)
178   }
179  
180   int
181 < intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
181 > intersect_ray_plane(orig,dir,peq,pd,r)
182     FVECT orig,dir;
183 <   FVECT plane_n;
173 <   double plane_d;
183 >   FPEQ peq;
184     double *pd;
185     FVECT r;
186   {
187 <  double t;
187 >  double t,d;
188    int hit;
189      /*
190        Plane is Ax + By + Cz +D = 0:
# Line 185 | Line 195 | intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
195         line is  l = p1 + (p2-p1)t
196       */
197      /* Solve for t: */
198 <    t =  -(DOT(plane_n,orig) + plane_d)/(DOT(plane_n,dir));
198 >  d = DOT(FP_N(peq),dir);
199 >  if(ZERO(d))
200 >     return(0);
201 >  t =  -(DOT(FP_N(peq),orig) + FP_D(peq))/d;
202 >
203 >  if(t < 0)
204 >       hit = 0;
205 >    else
206 >       hit = 1;
207 >
208 >  if(r)
209 >     VSUM(r,orig,dir,t);
210 >
211 >    if(pd)
212 >       *pd = t;
213 >  return(hit);
214 > }
215 >
216 >
217 > int
218 > intersect_ray_oplane(orig,dir,n,pd,r)
219 >   FVECT orig,dir;
220 >   FVECT n;
221 >   double *pd;
222 >   FVECT r;
223 > {
224 >  double t,d;
225 >  int hit;
226 >    /*
227 >      Plane is Ax + By + Cz +D = 0:
228 >      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
229 >    */
230 >     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
231 >         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
232 >       line is  l = p1 + (p2-p1)t
233 >     */
234 >    /* Solve for t: */
235 >    d= DOT(n,dir);
236 >    if(ZERO(d))
237 >       return(0);
238 >    t =  -(DOT(n,orig))/d;
239      if(t < 0)
240         hit = 0;
241      else
# Line 199 | Line 249 | intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
249    return(hit);
250   }
251  
252 + /* Assumption: know crosses plane:dont need to check for 'on' case */
253 + intersect_edge_coord_plane(v0,v1,w,r)
254 + FVECT v0,v1;
255 + int w;
256 + FVECT r;
257 + {
258 +  FVECT dv;
259 +  int wnext;
260 +  double t;
261  
262 +  VSUB(dv,v1,v0);
263 +  t = -v0[w]/dv[w];
264 +  r[w] = 0.0;
265 +  wnext = (w+1)%3;
266 +  r[wnext] = v0[wnext] + dv[wnext]*t;
267 +  wnext = (w+2)%3;
268 +  r[wnext] = v0[wnext] + dv[wnext]*t;
269 + }
270 +
271   int
272 < intersect_edge_plane(e0,e1,plane_n,plane_d,pd,r)
272 > intersect_edge_plane(e0,e1,peq,pd,r)
273     FVECT e0,e1;
274 <   FVECT plane_n;
207 <   double plane_d;
274 >   FPEQ peq;
275     double *pd;
276     FVECT r;
277   {
# Line 221 | Line 288 | intersect_edge_plane(e0,e1,plane_n,plane_d,pd,r)
288       */
289      /* Solve for t: */
290    VSUB(d,e1,e0);
291 <  t =  -(DOT(plane_n,e0) + plane_d)/(DOT(plane_n,d));
291 >  t =  -(DOT(FP_N(peq),e0) + FP_D(peq))/(DOT(FP_N(peq),d));
292      if(t < 0)
293         hit = 0;
294      else
# Line 240 | Line 307 | point_in_cone(p,p0,p1,p2)
307   FVECT p;
308   FVECT p0,p1,p2;
309   {
243 –    FVECT n;
310      FVECT np,x_axis,y_axis;
311 <    double d1,d2,d;
311 >    double d1,d2;
312 >    FPEQ peq;
313      
314      /* Find the equation of the circle defined by the intersection
315         of the cone with the plane defined by p1,p2,p3- project p into
316         that plane and do an in-circle test in the plane
317       */
318      
319 <    /* find the equation of the plane defined by p1-p3 */
320 <    tri_plane_equation(p0,p1,p2,n,&d,FALSE);
319 >    /* find the equation of the plane defined by p0-p2 */
320 >    tri_plane_equation(p0,p1,p2,&peq,FALSE);
321  
322      /* define a coordinate system on the plane: the x axis is in
323         the direction of np2-np1, and the y axis is calculated from
# Line 258 | Line 325 | FVECT p0,p1,p2;
325       */
326      /* Project p onto the plane */
327      /* NOTE: check this: does sideness check?*/
328 <    if(!intersect_vector_plane(p,n,d,NULL,np))
328 >    if(!intersect_vector_plane(p,peq,NULL,np))
329          return(FALSE);
330  
331 <    /* create coordinate system on  plane: p2-p1 defines the x_axis*/
331 >    /* create coordinate system on  plane: p1-p0 defines the x_axis*/
332      VSUB(x_axis,p1,p0);
333      normalize(x_axis);
334      /* The y axis is  */
335 <    VCROSS(y_axis,n,x_axis);
335 >    VCROSS(y_axis,FP_N(peq),x_axis);
336      normalize(y_axis);
337  
338      VSUB(p1,p1,p0);
339      VSUB(p2,p2,p0);
340      VSUB(np,np,p0);
341      
342 <    p1[0] = VLEN(p1);
342 >    p1[0] = DOT(p1,x_axis);
343      p1[1] = 0;
344  
345      d1 = DOT(p2,x_axis);
# Line 301 | Line 368 | int sides[3];
368      /* Find the normal to the triangle ORIGIN:v0:v1 */
369      if(!NTH_BIT(*nset,0))
370      {
371 <        VCROSS(n[0],v1,v0);
371 >        VCROSS(n[0],v0,v1);
372          SET_NTH_BIT(*nset,0);
373      }
374      /* Test the point for sidedness */
# Line 319 | Line 386 | int sides[3];
386      /* Test next edge */
387      if(!NTH_BIT(*nset,1))
388      {
389 <        VCROSS(n[1],v2,v1);
389 >        VCROSS(n[1],v1,v2);
390          SET_NTH_BIT(*nset,1);
391      }
392      /* Test the point for sidedness */
# Line 335 | Line 402 | int sides[3];
402      /* Test next edge */
403      if(!NTH_BIT(*nset,2))
404      {
405 <        VCROSS(n[2],v0,v2);
405 >        VCROSS(n[2],v2,v0);
406          SET_NTH_BIT(*nset,2);
407      }
408      /* Test the point for sidedness */
# Line 353 | Line 420 | int sides[3];
420  
421  
422  
423 <
423 >
424   int
425   point_in_stri(v0,v1,v2,p)
426   FVECT v0,v1,v2,p;
# Line 361 | Line 428 | FVECT v0,v1,v2,p;
428      double d;
429      FVECT n;  
430  
431 <    VCROSS(n,v1,v0);
431 >    VCROSS(n,v0,v1);
432      /* Test the point for sidedness */
433      d  = DOT(n,p);
434      if(d > 0.0)
435        return(FALSE);
436  
437      /* Test next edge */
438 <    VCROSS(n,v2,v1);
438 >    VCROSS(n,v1,v2);
439      /* Test the point for sidedness */
440      d  = DOT(n,p);
441      if(d > 0.0)
442         return(FALSE);
443  
444      /* Test next edge */
445 <    VCROSS(n,v0,v2);
445 >    VCROSS(n,v2,v0);
446      /* Test the point for sidedness */
447      d  = DOT(n,p);
448      if(d > 0.0)
# Line 469 | Line 536 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
536      if(sides[0][0] == GT_INVALID)
537      {
538        if(!NTH_BIT(nset,0))
539 <        VCROSS(n[0],t1,t0);
539 >        VCROSS(n[0],t0,t1);
540        /* Test the point for sidedness */
541        d  = DOT(n[0],p0);
542        if(d >= 0.0)
# Line 482 | Line 549 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
549      if(sides[0][1] == GT_INVALID)
550      {
551        if(!NTH_BIT(nset,1))
552 <        VCROSS(n[1],t2,t1);
552 >        VCROSS(n[1],t1,t2);
553          /* Test the point for sidedness */
554          d  = DOT(n[1],p0);
555          if(d >= 0.0)
# Line 495 | Line 562 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
562      if(sides[0][2] == GT_INVALID)
563      {
564        if(!NTH_BIT(nset,2))
565 <        VCROSS(n[2],t0,t2);
565 >        VCROSS(n[2],t2,t0);
566        /* Test the point for sidedness */
567        d  = DOT(n[2],p0);
568        if(d >= 0.0)
# Line 511 | Line 578 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
578      if(sides[1][0] == GT_INVALID)
579      {
580        if(!NTH_BIT(nset,0))
581 <        VCROSS(n[0],t1,t0);
581 >        VCROSS(n[0],t0,t1);
582        /* Test the point for sidedness */
583        d  = DOT(n[0],p1);
584        if(d >= 0.0)
# Line 525 | Line 592 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
592      if(sides[1][1] == GT_INVALID)
593      {
594        if(!NTH_BIT(nset,1))
595 <        VCROSS(n[1],t2,t1);
595 >        VCROSS(n[1],t1,t2);
596        /* Test the point for sidedness */
597        d  = DOT(n[1],p1);
598        if(d >= 0.0)
# Line 539 | Line 606 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
606      if(sides[1][2] == GT_INVALID)
607      {
608        if(!NTH_BIT(nset,2))
609 <        VCROSS(n[2],t0,t2);
609 >        VCROSS(n[2],t2,t0);
610        /* Test the point for sidedness */
611        d  = DOT(n[2],p1);
612        if(d >= 0.0)
# Line 555 | Line 622 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
622      if(sides[2][0] == GT_INVALID)
623      {
624        if(!NTH_BIT(nset,0))
625 <        VCROSS(n[0],t1,t0);
625 >        VCROSS(n[0],t0,t1);
626        /* Test the point for sidedness */
627        d  = DOT(n[0],p2);
628        if(d >= 0.0)
# Line 568 | Line 635 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
635      if(sides[2][1] == GT_INVALID)
636      {
637        if(!NTH_BIT(nset,1))
638 <        VCROSS(n[1],t2,t1);
638 >        VCROSS(n[1],t1,t2);
639        /* Test the point for sidedness */
640        d  = DOT(n[1],p2);
641        if(d >= 0.0)
# Line 581 | Line 648 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
648      if(sides[2][2] == GT_INVALID)
649      {
650        if(!NTH_BIT(nset,2))
651 <        VCROSS(n[2],t0,t2);
651 >        VCROSS(n[2],t2,t0);
652        /* Test the point for sidedness */
653        d  = DOT(n[2],p2);
654        if(d >= 0.0)
# Line 600 | Line 667 | FVECT a1,a2,a3,b1,b2,b3;
667    int which,test,n_set[2];
668    int sides[2][3][3],i,j,inext,jnext;
669    int tests[2][3];
670 <  FVECT n[2][3],p,avg[2];
670 >  FVECT n[2][3],p,avg[2],t1,t2,t3;
671  
672 + #ifdef DEBUG
673 +    tri_normal(b1,b2,b3,p,FALSE);
674 +    if(DOT(p,b1) > 0)
675 +      {
676 +       VCOPY(t3,b1);
677 +       VCOPY(t2,b2);
678 +       VCOPY(t1,b3);
679 +      }
680 +    else
681 + #endif
682 +      {
683 +       VCOPY(t1,b1);
684 +       VCOPY(t2,b2);
685 +       VCOPY(t3,b3);
686 +      }
687 +
688    /* Test the vertices of triangle a against the pyramid formed by triangle
689       b and the origin. If any vertex of a is interior to triangle b, or
690       if all 3 vertices of a are ON the edges of b,return TRUE. Remember
691       the results of the edge normal and sidedness tests for later.
692     */
693 < if(vertices_in_stri(a1,a2,a3,b1,b2,b3,&(n_set[0]),n[0],avg[0],sides[1]))
693 > if(vertices_in_stri(a1,a2,a3,t1,t2,t3,&(n_set[0]),n[0],avg[0],sides[1]))
694       return(TRUE);
695    
696 < if(vertices_in_stri(b1,b2,b3,a1,a2,a3,&(n_set[1]),n[1],avg[1],sides[0]))
696 > if(vertices_in_stri(t1,t2,t3,a1,a2,a3,&(n_set[1]),n[1],avg[1],sides[0]))
697       return(TRUE);
698  
699  
700 <  set_sidedness_tests(b1,b2,b3,a1,a2,a3,tests[0],sides[0],n_set[1],n[1]);
700 >  set_sidedness_tests(t1,t2,t3,a1,a2,a3,tests[0],sides[0],n_set[1],n[1]);
701    if(tests[0][0]&tests[0][1]&tests[0][2])
702      return(FALSE);
703  
704 <  set_sidedness_tests(a1,a2,a3,b1,b2,b3,tests[1],sides[1],n_set[0],n[0]);
704 >  set_sidedness_tests(a1,a2,a3,t1,t2,t3,tests[1],sides[1],n_set[0],n[0]);
705    if(tests[1][0]&tests[1][1]&tests[1][2])
706      return(FALSE);
707  
# Line 663 | Line 746 | FVECT orig,dir;
746   FVECT v0,v1,v2;
747   FVECT pt;
748   {
749 <  FVECT p0,p1,p2,p,n;
750 <  double pd;
749 >  FVECT p0,p1,p2,p;
750 >  FPEQ peq;
751    int type;
752  
753 <  point_on_sphere(p0,v0,orig);
754 <  point_on_sphere(p1,v1,orig);
755 <  point_on_sphere(p2,v2,orig);
756 <  
753 >  VSUB(p0,v0,orig);
754 >  VSUB(p1,v1,orig);
755 >  VSUB(p2,v2,orig);
756 >
757    if(point_in_stri(p0,p1,p2,dir))
758    {
759        /* Intersect the ray with the triangle plane */
760 <      tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
761 <      return(intersect_ray_plane(orig,dir,n,pd,NULL,pt));
760 >      tri_plane_equation(v0,v1,v2,&peq,FALSE);
761 >      return(intersect_ray_plane(orig,dir,peq,NULL,pt));
762    }
763    return(FALSE);
764   }
# Line 736 | Line 819 | FVECT fnear[4],ffar[4];
819      ffar[3][2] =  width*nhv[2] - height*nvv[2] + far*ndv[2] + vp[2] ;
820   }
821  
822 + int
823 + max_index(v,r)
824 + FVECT v;
825 + double *r;
826 + {
827 +  double p[3];
828 +  int i;
829  
830 +  p[0] = fabs(v[0]);
831 +  p[1] = fabs(v[1]);
832 +  p[2] = fabs(v[2]);
833 +  i = (p[0]>=p[1])?((p[0]>=p[2])?0:2):((p[1]>=p[2])?1:2);  
834 +  if(r)
835 +    *r = p[i];
836 +  return(i);
837 + }
838  
839 + int
840 + closest_point_in_tri(p0,p1,p2,p,p0id,p1id,p2id)
841 + FVECT p0,p1,p2,p;
842 + int p0id,p1id,p2id;
843 + {
844 +    double d,d1;
845 +    int i;
846 +    
847 +    d =  DIST_SQ(p,p0);
848 +    d1 = DIST_SQ(p,p1);
849 +    if(d < d1)
850 +    {
851 +      d1 = DIST_SQ(p,p2);
852 +      i = (d1 < d)?p2id:p0id;
853 +    }
854 +    else
855 +    {
856 +      d = DIST_SQ(p,p2);
857 +      i = (d < d1)? p2id:p1id;
858 +    }
859 +    return(i);
860 + }
861  
862 +
863   int
864   sedge_intersect(a0,a1,b0,b1)
865   FVECT a0,a1,b0,b1;
# Line 786 | Line 907 | FVECT a0,a1,b0,b1;
907   }
908  
909  
789 –
910   /* Find the normalized barycentric coordinates of p relative to
911   * triangle v0,v1,v2. Return result in coord
912   */
# Line 800 | Line 920 | double coord[3];
920    a =  (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1);
921    coord[0] = ((x2 - px) * (y3 - py) - (x3 - px) * (y2 - py)) / a;
922    coord[1] = ((x3 - px) * (y1 - py) - (x1 - px) * (y3 - py)) / a;
923 <  coord[2]  = 1.0 - coord[0] - coord[1];
923 >  coord[2] = ((x1 - px) * (y2 - py) - (x2 - px) * (y1 - py)) / a;
924  
925   }
926  
807 – bary_ith_child(coord,i)
808 – double coord[3];
809 – int i;
810 – {
811 –
812 –  switch(i){
813 –  case 0:
814 –      /* update bary for child */
815 –      coord[0] = 2.0*coord[0]- 1.0;
816 –      coord[1] *= 2.0;
817 –      coord[2] *= 2.0;
818 –      break;
819 –  case 1:
820 –    coord[0] *= 2.0;
821 –    coord[1] = 2.0*coord[1]- 1.0;
822 –    coord[2] *= 2.0;
823 –    break;
824 –  case 2:
825 –    coord[0] *= 2.0;
826 –    coord[1] *= 2.0;
827 –    coord[2] = 2.0*coord[2]- 1.0;
828 –    break;
829 –  case 3:
830 –    coord[0] = 1.0 - 2.0*coord[0];
831 –    coord[1] = 1.0 - 2.0*coord[1];
832 –    coord[2] = 1.0 - 2.0*coord[2];
833 –    break;
834 – #ifdef DEBUG
835 –  default:
836 –    eputs("bary_ith_child():Invalid child\n");
837 –    break;
838 – #endif
839 –  }
840 – }
927  
928  
843 – int
844 – bary_child(coord)
845 – double coord[3];
846 – {
847 –  int i;
929  
849 –  if(coord[0] > 0.5)
850 –  {
851 –      /* update bary for child */
852 –      coord[0] = 2.0*coord[0]- 1.0;
853 –      coord[1] *= 2.0;
854 –      coord[2] *= 2.0;
855 –      return(0);
856 –  }
857 –  else
858 –    if(coord[1] > 0.5)
859 –    {
860 –      coord[0] *= 2.0;
861 –      coord[1] = 2.0*coord[1]- 1.0;
862 –      coord[2] *= 2.0;
863 –      return(1);
864 –    }
865 –    else
866 –      if(coord[2] > 0.5)
867 –      {
868 –        coord[0] *= 2.0;
869 –        coord[1] *= 2.0;
870 –        coord[2] = 2.0*coord[2]- 1.0;
871 –        return(2);
872 –      }
873 –      else
874 –         {
875 –           coord[0] = 1.0 - 2.0*coord[0];
876 –           coord[1] = 1.0 - 2.0*coord[1];
877 –           coord[2] = 1.0 - 2.0*coord[2];
878 –           return(3);
879 –         }
880 – }
881 –
882 – /* Coord was the ith child of the parent: set the coordinate
883 –   relative to the parent
884 – */
930   bary_parent(coord,i)
931 < double coord[3];
931 > BCOORD coord[3];
932   int i;
933   {
889 –
934    switch(i) {
935    case 0:
936      /* update bary for child */
937 <    coord[0] = coord[0]*0.5 + 0.5;
938 <    coord[1] *= 0.5;
939 <    coord[2] *= 0.5;
937 >    coord[0] = (coord[0] >> 1) + MAXBCOORD4;
938 >    coord[1] >>= 1;
939 >    coord[2] >>= 1;
940      break;
941    case 1:
942 <    coord[0] *= 0.5;
943 <    coord[1]  = coord[1]*0.5 + 0.5;
944 <    coord[2] *= 0.5;
942 >    coord[0] >>= 1;
943 >    coord[1]  = (coord[1] >> 1) + MAXBCOORD4;
944 >    coord[2] >>= 1;
945      break;
946      
947    case 2:
948 <    coord[0] *= 0.5;
949 <    coord[1] *= 0.5;
950 <    coord[2] = coord[2]*0.5 + 0.5;
948 >    coord[0] >>= 1;
949 >    coord[1] >>= 1;
950 >    coord[2] = (coord[2] >> 1) + MAXBCOORD4;
951      break;
952      
953    case 3:
954 <    coord[0] = 0.5 - 0.5*coord[0];
955 <    coord[1] = 0.5 - 0.5*coord[1];
956 <    coord[2] = 0.5 - 0.5*coord[2];
954 >    coord[0] = MAXBCOORD4 - (coord[0] >> 1);
955 >    coord[1] = MAXBCOORD4 - (coord[1] >> 1);
956 >    coord[2] = MAXBCOORD4 - (coord[2] >> 1);
957      break;
958   #ifdef DEBUG
959    default:
# Line 920 | Line 964 | int i;
964   }
965  
966   bary_from_child(coord,child,next)
967 < double coord[3];
967 > BCOORD coord[3];
968   int child,next;
969   {
970   #ifdef DEBUG
# Line 940 | Line 984 | int child,next;
984  
985    switch(child){
986    case 0:
987 <    switch(next){
988 <    case 1:
989 <      coord[0] += 1.0;
946 <      coord[1] -= 1.0;
987 >      coord[0] = 0;
988 >      coord[1] = MAXBCOORD2 - coord[1];
989 >      coord[2] = MAXBCOORD2 - coord[2];
990        break;
948 –    case 2:
949 –      coord[0] += 1.0;
950 –      coord[2] -= 1.0;
951 –      break;
952 –    case 3:
953 –      coord[0] *= -1.0;
954 –      coord[1] = 1 - coord[1];
955 –      coord[2] = 1 - coord[2];
956 –      break;
957 –
958 –    }
959 –    break;
991    case 1:
992 <    switch(next){
993 <    case 0:
994 <      coord[0] -= 1.0;
964 <      coord[1] += 1.0;
992 >      coord[0] = MAXBCOORD2 - coord[0];
993 >      coord[1] = 0;
994 >      coord[2] = MAXBCOORD2 - coord[2];
995        break;
966 –    case 2:
967 –      coord[1] += 1.0;
968 –      coord[2] -= 1.0;
969 –      break;
970 –    case 3:
971 –      coord[0] = 1 - coord[0];
972 –      coord[1] *= -1.0;
973 –      coord[2] = 1 - coord[2];
974 –      break;
975 –    }
976 –    break;
996    case 2:
997 <    switch(next){
998 <    case 0:
999 <      coord[0] -= 1.0;
981 <      coord[2] += 1.0;
982 <      break;
983 <    case 1:
984 <      coord[1] -= 1.0;
985 <      coord[2] += 1.0;
986 <      break;
987 <    case 3:
988 <      coord[0] = 1 - coord[0];
989 <      coord[1] = 1 - coord[1];
990 <      coord[2] *= -1.0;
991 <      break;
992 <    }
997 >      coord[0] = MAXBCOORD2 - coord[0];
998 >      coord[1] = MAXBCOORD2 - coord[1];
999 >      coord[2] = 0;
1000      break;
1001    case 3:
1002      switch(next){
1003      case 0:
1004 <      coord[0] *= -1.0;
1005 <      coord[1] = 1 - coord[1];
1006 <      coord[2] = 1 - coord[2];
1004 >      coord[0] = 0;
1005 >      coord[1] = MAXBCOORD2 - coord[1];
1006 >      coord[2] = MAXBCOORD2 - coord[2];
1007        break;
1008      case 1:
1009 <      coord[0] = 1 - coord[0];
1010 <      coord[1] *= -1.0;
1011 <      coord[2] = 1 - coord[2];
1009 >      coord[0] = MAXBCOORD2 - coord[0];
1010 >      coord[1] = 0;
1011 >      coord[2] = MAXBCOORD2 - coord[2];
1012        break;
1013      case 2:
1014 <      coord[0] = 1 - coord[0];
1015 <      coord[1] = 1 - coord[1];
1016 <      coord[2] *= -1.0;
1014 >      coord[0] = MAXBCOORD2 - coord[0];
1015 >      coord[1] = MAXBCOORD2 - coord[1];
1016 >      coord[2] = 0;
1017        break;
1018      }
1019      break;
# Line 1014 | Line 1021 | int child,next;
1021   }
1022  
1023   int
1024 < max_index(v)
1025 < FVECT v;
1024 > bary_child(coord)
1025 > BCOORD coord[3];
1026   {
1020 –  double a,b,c;
1027    int i;
1028  
1029 <  a = fabs(v[0]);
1030 <  b = fabs(v[1]);
1031 <  c = fabs(v[2]);
1032 <  i = (a>=b)?((a>=c)?0:2):((b>=c)?1:2);  
1033 <  return(i);
1034 < }
1035 <
1030 <
1031 <
1032 < /*
1033 < * int
1034 < * traceRay(FVECT orig, FVECT dir,FVECT v0,FVECT v1,FVECT v2,FVECT r)
1035 < *
1036 < *   Intersect the ray with triangle v0v1v2, return intersection point in r
1037 < *
1038 < *    Assumes orig,v0,v1,v2 are in spherical coordinates, and orig is
1039 < *    unit
1040 < */
1041 < int
1042 < traceRay(orig,dir,v0,v1,v2,r)
1043 <  FVECT orig,dir;
1044 <  FVECT v0,v1,v2;
1045 <  FVECT r;
1046 < {
1047 <  FVECT n,p[3],d;
1048 <  double pt[3],r_eps;
1049 <  int i;
1050 <  int which;
1051 <
1052 <  /* Find the plane equation for the triangle defined by the edge v0v1 and
1053 <   the view center*/
1054 <  VCROSS(n,v0,v1);
1055 <  /* Intersect the ray with this plane */
1056 <  i = intersect_ray_plane(orig,dir,n,0.0,&(pt[0]),p[0]);
1057 <  if(i)
1058 <    which = 0;
1059 <  else
1060 <    which = -1;
1061 <
1062 <  VCROSS(n,v1,v2);
1063 <  i = intersect_ray_plane(orig,dir,n,0.0,&(pt[1]),p[1]);
1064 <  if(i)
1065 <    if((which==-1) || pt[1] < pt[0])
1066 <      which = 1;
1067 <
1068 <  VCROSS(n,v2,v0);
1069 <  i = intersect_ray_plane(orig,dir,n,0.0,&(pt[2]),p[2]);
1070 <  if(i)
1071 <    if((which==-1) || pt[2] < pt[which])
1072 <      which = 2;
1073 <
1074 <  if(which != -1)
1075 <  {
1076 <      /* Return point that is K*eps along projection of the ray on
1077 <         the sphere to push intersection point p[which] into next cell
1078 <      */
1079 <      normalize(p[which]);
1080 <      /* Calculate the ray perpendicular to the intersection point: approx
1081 <       the direction moved along the sphere a distance of K*epsilon*/
1082 <      r_eps = -DOT(p[which],dir);
1083 <      VSUM(n,dir,p[which],r_eps);
1084 <     /* Calculate the length along ray p[which]-dir needed to move to
1085 <         cause a move along the sphere of k*epsilon
1086 <       */
1087 <       r_eps = DOT(n,dir);
1088 <      VSUM(r,p[which],dir,(20*FTINY)/r_eps);
1089 <      normalize(r);
1090 <      return(TRUE);
1029 >  if(coord[0] > MAXBCOORD4)
1030 >  {
1031 >      /* update bary for child */
1032 >      coord[0] = (coord[0]<< 1) - MAXBCOORD2;
1033 >      coord[1] <<= 1;
1034 >      coord[2] <<= 1;
1035 >      return(0);
1036    }
1037    else
1038 <  {
1094 <      /* Unable to find intersection: move along ray and try again */
1095 <      r_eps = -DOT(orig,dir);
1096 <      VSUM(n,dir,orig,r_eps);
1097 <      r_eps = DOT(n,dir);
1098 <      VSUM(r,orig,dir,(20*FTINY)/r_eps);
1099 <      normalize(r);
1100 < #ifdef DEBUG
1101 <      eputs("traceRay:Ray does not intersect triangle\n");
1102 < #endif
1103 <      return(FALSE);
1104 <  }
1105 < }
1106 <
1107 <
1108 < int
1109 < closest_point_in_tri(p0,p1,p2,p,p0id,p1id,p2id)
1110 < FVECT p0,p1,p2,p;
1111 < int p0id,p1id,p2id;
1112 < {
1113 <    double d,d1;
1114 <    int i;
1115 <    
1116 <    d =  DIST_SQ(p,p0);
1117 <    d1 = DIST_SQ(p,p1);
1118 <    if(d < d1)
1038 >    if(coord[1] > MAXBCOORD4)
1039      {
1040 <      d1 = DIST_SQ(p,p2);
1041 <      i = (d1 < d)?p2id:p0id;
1040 >      coord[0] <<= 1;
1041 >      coord[1] = (coord[1] << 1) - MAXBCOORD2;
1042 >      coord[2] <<= 1;
1043 >      return(1);
1044      }
1045      else
1046 <    {
1047 <      d = DIST_SQ(p,p2);
1048 <      i = (d < d1)? p2id:p1id;
1049 <    }
1050 <    return(i);
1046 >      if(coord[2] > MAXBCOORD4)
1047 >      {
1048 >        coord[0] <<= 1;
1049 >        coord[1] <<= 1;
1050 >        coord[2] = (coord[2] << 1) - MAXBCOORD2;
1051 >        return(2);
1052 >      }
1053 >      else
1054 >         {
1055 >           coord[0] = MAXBCOORD2 - (coord[0] << 1);
1056 >           coord[1] = MAXBCOORD2 - (coord[1] << 1);
1057 >           coord[2] = MAXBCOORD2 - (coord[2] << 1);
1058 >           return(3);
1059 >         }
1060   }
1061  
1062 + int
1063 + bary_nth_child(coord,i)
1064 + BCOORD coord[3];
1065 + int i;
1066 + {
1067  
1068 <
1069 <
1068 >  switch(i){
1069 >  case 0:
1070 >    /* update bary for child */
1071 >    coord[0] = (coord[0]<< 1) - MAXBCOORD2;
1072 >    coord[1] <<= 1;
1073 >    coord[2] <<= 1;
1074 >    break;
1075 >  case 1:
1076 >    coord[0] <<= 1;
1077 >    coord[1] = (coord[1] << 1) - MAXBCOORD2;
1078 >    coord[2] <<= 1;
1079 >    break;
1080 >  case 2:
1081 >    coord[0] <<= 1;
1082 >    coord[1] <<= 1;
1083 >    coord[2] = (coord[2] << 1) - MAXBCOORD2;
1084 >    break;
1085 >  case 3:
1086 >    coord[0] = MAXBCOORD2 - (coord[0] << 1);
1087 >    coord[1] = MAXBCOORD2 - (coord[1] << 1);
1088 >    coord[2] = MAXBCOORD2 - (coord[2] << 1);
1089 >    break;
1090 >  }
1091 > }
1092  
1093  
1094  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines