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.1 by gwlarson, Wed Aug 19 17:45:24 1998 UTC vs.
Revision 3.12 by gwlarson, Thu Jun 10 15:22:22 1999 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 +    double dp,dp1;
40 +    int test,test1;
41 +    FVECT nv0,nv1,nv2;
42      FVECT cp01,cp12,cp;
43 <    
43 >
44      /* test sign of (v0Xv1)X(v1Xv2). v1 */
45 <    VCROSS(cp01,v0,v1);
46 <    VCROSS(cp12,v1,v2);
45 >    /* un-Simplified: */
46 >    VCOPY(nv0,v0);
47 >    normalize(nv0);
48 >    VCOPY(nv1,v1);
49 >    normalize(nv1);
50 >    VCOPY(nv2,v2);
51 >    normalize(nv2);
52 >
53 >    VCROSS(cp01,nv0,nv1);
54 >    VCROSS(cp12,nv1,nv2);
55      VCROSS(cp,cp01,cp12);
56 <    if(DOT(cp,v1) < 0)
57 <       return(FALSE);
58 <    return(TRUE);
56 >    normalize(cp);
57 >    dp1 = DOT(cp,nv1);
58 >    if(dp1 <= 1e-8 || dp1 >= (1-1e-8)) /* Test if on other side,or colinear*/
59 >      test1 = FALSE;
60 >    else
61 >      test1 = TRUE;
62 >
63 >    dp =   nv0[2]*nv1[0]*nv2[1] -  nv0[2]*nv1[1]*nv2[0] - nv0[0]*nv1[2]*nv2[1]
64 >         + nv0[0]*nv1[1]*nv2[2] +  nv0[1]*nv1[2]*nv2[0] - nv0[1]*nv1[0]*nv2[2];
65 >    
66 >    if(dp <= 1e-8 || dp1 >= (1-1e-8))
67 >      test = FALSE;
68 >    else
69 >      test = TRUE;
70 >
71 >    if(test != test1)
72 >      fprintf(stderr,"test %f simplified  %f\n",dp1,dp);
73 >    return(test1);
74   }
75  
76   /* calculates the normal of a face contour using Newell's formula. e
77  
78 <               a =  SUMi (yi - yi+1)(zi + zi+1)
78 >               a =  SUMi (yi - yi+1)(zi + zi+1);
79                 b =  SUMi (zi - zi+1)(xi + xi+1)
80                 c =  SUMi (xi - xi+1)(yi + yi+1)
81   */
82   double
83   tri_normal(v0,v1,v2,n,norm)
84   FVECT v0,v1,v2,n;
85 < char norm;
85 > int norm;
86   {
87    double mag;
88  
# Line 63 | Line 92 | char norm;
92    
93    n[1] = (v0[2] - v1[2]) * (v0[0] + v1[0]) +
94             (v1[2] - v2[2]) * (v1[0] + v2[0]) +
95 <           (v2[2] - v0[2]) * (v2[0] + v0[0]);
67 <
95 >            (v2[2] - v0[2]) * (v2[0] + v0[0]);
96    
97    n[2] = (v0[1] + v1[1]) * (v0[0] - v1[0]) +
98           (v1[1] + v2[1]) * (v1[0] - v2[0]) +
# Line 72 | Line 100 | char norm;
100  
101    if(!norm)
102       return(0);
75 –
103    
104    mag = normalize(n);
105  
# Line 80 | Line 107 | char norm;
107   }
108  
109  
110 < tri_plane_equation(v0,v1,v2,n,nd,norm)
111 <   FVECT v0,v1,v2,n;
112 <   double *nd;
113 <   char norm;
110 > tri_plane_equation(v0,v1,v2,peqptr,norm)
111 >   FVECT v0,v1,v2;
112 >   FPEQ *peqptr;
113 >   int norm;
114   {
115 <    tri_normal(v0,v1,v2,n,norm);
115 >    tri_normal(v0,v1,v2,FP_N(*peqptr),norm);
116  
117 <    *nd = -(DOT(n,v0));
117 >    FP_D(*peqptr) = -(DOT(FP_N(*peqptr),v0));
118   }
119  
93 – int
94 – point_relative_to_plane(p,n,nd)
95 –   FVECT p,n;
96 –   double nd;
97 – {
98 –    double d;
99 –    
100 –    d = p[0]*n[0] + p[1]*n[1] + p[2]*n[2] + nd;
101 –    if(d < 0)
102 –       return(-1);
103 –    if(ZERO(d))
104 –       return(0);
105 –    else
106 –       return(1);
107 – }
120  
121 < /* From quad_edge-code */
121 > /* returns TRUE if ray from origin in direction v intersects plane defined
122 >  * by normal plane_n, and plane_d. If plane is not parallel- returns
123 >  * intersection point if r != NULL. If tptr!= NULL returns value of
124 >  * t, if parallel, returns t=FHUGE
125 >  */
126   int
127 < point_in_circle_thru_origin(p,p0,p1)
128 < FVECT p;
129 < FVECT p0,p1;
127 > intersect_vector_plane(v,peq,tptr,r)
128 >   FVECT v;
129 >   FPEQ peq;
130 >   double *tptr;
131 >   FVECT r;
132   {
133 +  double t,d;
134 +  int hit;
135 +    /*
136 +      Plane is Ax + By + Cz +D = 0:
137 +      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
138 +    */
139  
140 <    double dp0,dp1;
117 <    double dp,det;
118 <    
119 <    dp0 = DOT_VEC2(p0,p0);
120 <    dp1 = DOT_VEC2(p1,p1);
121 <    dp  = DOT_VEC2(p,p);    
140 >    /* line is  l = p1 + (p2-p1)t, p1=origin */
141  
142 <    det = -dp0*CROSS_VEC2(p1,p) + dp1*CROSS_VEC2(p0,p) - dp*CROSS_VEC2(p0,p1);
143 <    
144 <    return (det > 0);
142 >    /* Solve for t: */
143 >  d = -(DOT(FP_N(peq),v));
144 >  if(ZERO(d))
145 >  {
146 >      t = FHUGE;
147 >      hit = 0;
148 >  }
149 >  else
150 >  {
151 >      t =  FP_D(peq)/d;
152 >      if(t < 0 )
153 >         hit = 0;
154 >      else
155 >         hit = 1;
156 >      if(r)
157 >         {
158 >             r[0] = v[0]*t;
159 >             r[1] = v[1]*t;
160 >             r[2] = v[2]*t;
161 >         }
162 >  }
163 >    if(tptr)
164 >       *tptr = t;
165 >  return(hit);
166   }
167  
128 –
129 –
130 – point_on_sphere(ps,p,c)
131 – FVECT ps,p,c;
132 – {
133 –    VSUB(ps,p,c);    
134 –    normalize(ps);
135 – }
136 –
137 –
168   int
169 < intersect_vector_plane(v,plane_n,plane_d,pd,r)
170 <   FVECT v,plane_n;
171 <   double plane_d;
169 > intersect_ray_plane(orig,dir,peq,pd,r)
170 >   FVECT orig,dir;
171 >   FPEQ peq;
172     double *pd;
173     FVECT r;
174   {
175 <  double t;
175 >  double t,d;
176    int hit;
177      /*
178        Plane is Ax + By + Cz +D = 0:
179        plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
180      */
181 <
182 <    /* line is  l = p1 + (p2-p1)t, p1=origin */
183 <
181 >     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
182 >         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
183 >       line is  l = p1 + (p2-p1)t
184 >     */
185      /* Solve for t: */
186 <    t =  plane_d/-(DOT(plane_n,v));
187 <    if(t >0 || ZERO(t))
188 <       hit = 1;
189 <    else
186 >  d = DOT(FP_N(peq),dir);
187 >  if(ZERO(d))
188 >     return(0);
189 >  t =  -(DOT(FP_N(peq),orig) + FP_D(peq))/d;
190 >
191 >  if(t < 0)
192         hit = 0;
193 <    r[0] = v[0]*t;
194 <    r[1] = v[1]*t;
195 <    r[2] = v[2]*t;
193 >    else
194 >       hit = 1;
195 >
196 >  if(r)
197 >     VSUM(r,orig,dir,t);
198 >
199      if(pd)
200         *pd = t;
201    return(hit);
202   }
203  
204 +
205   int
206 < intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
206 > intersect_ray_oplane(orig,dir,n,pd,r)
207     FVECT orig,dir;
208 <   FVECT plane_n;
172 <   double plane_d;
208 >   FVECT n;
209     double *pd;
210     FVECT r;
211   {
212 <  double t;
212 >  double t,d;
213    int hit;
214      /*
215        Plane is Ax + By + Cz +D = 0:
216        plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
217      */
218 <     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0 */
219 <    /* t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
220 <    /* line is  l = p1 + (p2-p1)t */
218 >     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
219 >         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
220 >       line is  l = p1 + (p2-p1)t
221 >     */
222      /* Solve for t: */
223 <    t =  -(DOT(plane_n,orig) + plane_d)/(DOT(plane_n,dir));
224 <    if(ZERO(t) || t >0)
225 <       hit = 1;
226 <    else
223 >    d= DOT(n,dir);
224 >    if(ZERO(d))
225 >       return(0);
226 >    t =  -(DOT(n,orig))/d;
227 >    if(t < 0)
228         hit = 0;
229 +    else
230 +       hit = 1;
231  
232 <  VSUM(r,orig,dir,t);
232 >  if(r)
233 >     VSUM(r,orig,dir,t);
234  
235      if(pd)
236         *pd = t;
237    return(hit);
238   }
239  
240 + /* Assumption: know crosses plane:dont need to check for 'on' case */
241 + intersect_edge_coord_plane(v0,v1,w,r)
242 + FVECT v0,v1;
243 + int w;
244 + FVECT r;
245 + {
246 +  FVECT dv;
247 +  int wnext;
248 +  double t;
249  
250 +  VSUB(dv,v1,v0);
251 +  t = -v0[w]/dv[w];
252 +  r[w] = 0.0;
253 +  wnext = (w+1)%3;
254 +  r[wnext] = v0[wnext] + dv[wnext]*t;
255 +  wnext = (w+2)%3;
256 +  r[wnext] = v0[wnext] + dv[wnext]*t;
257 + }
258 +
259   int
260 < point_in_cone(p,p0,p1,p2)
261 < FVECT p;
262 < FVECT p0,p1,p2;
260 > intersect_edge_plane(e0,e1,peq,pd,r)
261 >   FVECT e0,e1;
262 >   FPEQ peq;
263 >   double *pd;
264 >   FVECT r;
265   {
266 <    FVECT n;
267 <    FVECT np,x_axis,y_axis;
268 <    double d1,d2,d;
269 <    
270 <    /* Find the equation of the circle defined by the intersection
271 <       of the cone with the plane defined by p1,p2,p3- project p into
272 <       that plane and do an in-circle test in the plane
266 >  double t;
267 >  int hit;
268 >  FVECT d;
269 >  /*
270 >      Plane is Ax + By + Cz +D = 0:
271 >      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
272 >    */
273 >     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
274 >         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
275 >       line is  l = p1 + (p2-p1)t
276       */
277 <    
278 <    /* find the equation of the plane defined by p1-p3 */
279 <    tri_plane_equation(p0,p1,p2,n,&d,FALSE);
277 >    /* Solve for t: */
278 >  VSUB(d,e1,e0);
279 >  t =  -(DOT(FP_N(peq),e0) + FP_D(peq))/(DOT(FP_N(peq),d));
280 >    if(t < 0)
281 >       hit = 0;
282 >    else
283 >       hit = 1;
284  
285 <    /* define a coordinate system on the plane: the x axis is in
218 <       the direction of np2-np1, and the y axis is calculated from
219 <       n cross x-axis
220 <     */
221 <    /* Project p onto the plane */
222 <    if(!intersect_vector_plane(p,n,d,NULL,np))
223 <        return(FALSE);
285 >  VSUM(r,e0,d,t);
286  
287 <    /* create coordinate system on  plane: p2-p1 defines the x_axis*/
288 <    VSUB(x_axis,p1,p0);
289 <    normalize(x_axis);
228 <    /* The y axis is  */
229 <    VCROSS(y_axis,n,x_axis);
230 <    normalize(y_axis);
231 <
232 <    VSUB(p1,p1,p0);
233 <    VSUB(p2,p2,p0);
234 <    VSUB(np,np,p0);
235 <    
236 <    p1[0] = VLEN(p1);
237 <    p1[1] = 0;
238 <
239 <    d1 = DOT(p2,x_axis);
240 <    d2 = DOT(p2,y_axis);
241 <    p2[0] = d1;
242 <    p2[1] = d2;
243 <    
244 <    d1 = DOT(np,x_axis);
245 <    d2 = DOT(np,y_axis);
246 <    np[0] = d1;
247 <    np[1] = d2;
248 <
249 <    /* perform the in-circle test in the new coordinate system */
250 <    return(point_in_circle_thru_origin(np,p1,p2));
287 >    if(pd)
288 >       *pd = t;
289 >  return(hit);
290   }
291  
292   int
293 < test_point_against_spherical_tri(v0,v1,v2,p,n,nset,which,sides)
293 > point_set_in_stri(v0,v1,v2,p,n,nset,sides)
294   FVECT v0,v1,v2,p;
295   FVECT n[3];
296 < char *nset;
297 < char *which;
259 < char sides[3];
296 > int *nset;
297 > int sides[3];
298  
299   {
300 <    float d;
263 <
300 >    double d;
301      /* Find the normal to the triangle ORIGIN:v0:v1 */
302      if(!NTH_BIT(*nset,0))
303      {
304 <        VCROSS(n[0],v1,v0);
304 >        VCROSS(n[0],v0,v1);
305          SET_NTH_BIT(*nset,0);
306      }
307      /* Test the point for sidedness */
308      d  = DOT(n[0],p);
309  
310 <    if(ZERO(d))
311 <       sides[0] = GT_EDGE;
312 <    else
313 <       if(d > 0)
314 <      {
278 <          sides[0] =  GT_OUT;
279 <          sides[1] = sides[2] = GT_INVALID;
280 <          return(FALSE);
310 >    if(d > 0.0)
311 >     {
312 >       sides[0] =  GT_OUT;
313 >       sides[1] = sides[2] = GT_INVALID;
314 >       return(FALSE);
315        }
316      else
317         sides[0] = GT_INTERIOR;
# Line 285 | Line 319 | char sides[3];
319      /* Test next edge */
320      if(!NTH_BIT(*nset,1))
321      {
322 <        VCROSS(n[1],v2,v1);
322 >        VCROSS(n[1],v1,v2);
323          SET_NTH_BIT(*nset,1);
324      }
325      /* Test the point for sidedness */
326      d  = DOT(n[1],p);
327 <    if(ZERO(d))
327 >    if(d > 0.0)
328      {
295 –        sides[1] = GT_EDGE;
296 –        /* If on plane 0-and on plane 1: lies on edge */
297 –        if(sides[0] == GT_EDGE)
298 –        {
299 –            *which = 1;
300 –            sides[2] = GT_INVALID;
301 –            return(GT_EDGE);
302 –        }
303 –    }
304 –    else if(d > 0)
305 –    {
329          sides[1] = GT_OUT;
330          sides[2] = GT_INVALID;
331          return(FALSE);
# Line 312 | Line 335 | char sides[3];
335      /* Test next edge */
336      if(!NTH_BIT(*nset,2))
337      {
338 <
316 <        VCROSS(n[2],v0,v2);
338 >        VCROSS(n[2],v2,v0);
339          SET_NTH_BIT(*nset,2);
340      }
341      /* Test the point for sidedness */
342      d  = DOT(n[2],p);
343 <    if(ZERO(d))
343 >    if(d > 0.0)
344      {
345 <        sides[2] = GT_EDGE;
346 <
325 <        /* If on plane 0 and 2: lies on edge 0*/
326 <        if(sides[0] == GT_EDGE)
327 <           {
328 <               *which = 0;
329 <               return(GT_EDGE);
330 <           }
331 <        /* If on plane 1 and 2: lies on edge  2*/
332 <        if(sides[1] == GT_EDGE)
333 <           {
334 <               *which = 2;
335 <               return(GT_EDGE);
336 <           }
337 <        /* otherwise: on face 2 */
338 <        else
339 <           {
340 <               *which = 2;
341 <               return(GT_FACE);
342 <           }
345 >      sides[2] = GT_OUT;
346 >      return(FALSE);
347      }
344 –    else if(d > 0)
345 –      {
346 –        sides[2] = GT_OUT;
347 –        return(FALSE);
348 –      }
349 –    /* If on edge */
348      else
349         sides[2] = GT_INTERIOR;
352 –    
353 –    /* If on plane 0 only: on face 0 */
354 –    if(sides[0] == GT_EDGE)
355 –    {
356 –        *which = 0;
357 –        return(GT_FACE);
358 –    }
359 –    /* If on plane 1 only: on face 1 */
360 –    if(sides[1] == GT_EDGE)
361 –    {
362 –        *which = 1;
363 –        return(GT_FACE);
364 –    }
350      /* Must be interior to the pyramid */
351      return(GT_INTERIOR);
352   }
353  
354  
355  
371 –
372 – int
373 – test_single_point_against_spherical_tri(v0,v1,v2,p,which)
374 – FVECT v0,v1,v2,p;
375 – char *which;
376 – {
377 –    float d;
378 –    FVECT n;  
379 –    char sides[3];
380 –
381 –    /* First test if point coincides with any of the vertices */
382 –    if(EQUAL_VEC3(p,v0))
383 –    {
384 –        *which = 0;
385 –        return(GT_VERTEX);
386 –    }
387 –    if(EQUAL_VEC3(p,v1))
388 –    {
389 –        *which = 1;
390 –        return(GT_VERTEX);
391 –    }
392 –    if(EQUAL_VEC3(p,v2))
393 –    {
394 –        *which = 2;
395 –        return(GT_VERTEX);
396 –    }
397 –    VCROSS(n,v1,v0);
398 –    /* Test the point for sidedness */
399 –    d  = DOT(n,p);
400 –    if(ZERO(d))
401 –       sides[0] = GT_EDGE;
402 –    else
403 –       if(d > 0)
404 –          return(FALSE);
405 –       else
406 –          sides[0] = GT_INTERIOR;
407 –    /* Test next edge */
408 –    VCROSS(n,v2,v1);
409 –    /* Test the point for sidedness */
410 –    d  = DOT(n,p);
411 –    if(ZERO(d))
412 –    {
413 –        sides[1] = GT_EDGE;
414 –        /* If on plane 0-and on plane 1: lies on edge */
415 –        if(sides[0] == GT_EDGE)
416 –        {
417 –            *which = 1;
418 –            return(GT_VERTEX);
419 –        }
420 –    }
421 –    else if(d > 0)
422 –       return(FALSE);
423 –    else
424 –       sides[1] = GT_INTERIOR;
425 –
426 –    /* Test next edge */
427 –    VCROSS(n,v0,v2);
428 –    /* Test the point for sidedness */
429 –    d  = DOT(n,p);
430 –    if(ZERO(d))
431 –    {
432 –        sides[2] = GT_EDGE;
433 –        
434 –        /* If on plane 0 and 2: lies on edge 0*/
435 –        if(sides[0] == GT_EDGE)
436 –        {
437 –            *which = 0;
438 –            return(GT_VERTEX);
439 –        }
440 –        /* If on plane 1 and 2: lies on edge  2*/
441 –        if(sides[1] == GT_EDGE)
442 –        {
443 –            *which = 2;
444 –            return(GT_VERTEX);
445 –        }
446 –        /* otherwise: on face 2 */
447 –        else
448 –       {
449 –           return(GT_FACE);
450 –       }
451 –    }
452 –    else if(d > 0)
453 –       return(FALSE);
454 –    /* Must be interior to the pyramid */
455 –    return(GT_FACE);
456 – }
457 –
458 – int
459 – test_vertices_for_tri_inclusion(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides)
460 – FVECT t0,t1,t2,p0,p1,p2;
461 – char *nset;
462 – FVECT n[3];
463 – FVECT avg;
464 – char pt_sides[3][3];
465 –
466 – {
467 –    char below_plane[3],on_edge,test;
468 –    char which;
469 –
470 –    SUM_3VEC3(avg,t0,t1,t2);
471 –    on_edge = 0;
472 –    *nset = 0;
473 –    /* Test vertex v[i] against triangle j*/
474 –    /* Check if v[i] lies below plane defined by avg of 3 vectors
475 –       defining triangle
476 –       */
477 –
478 –    /* test point 0 */
479 –    if(DOT(avg,p0) < 0)
480 –      below_plane[0] = 1;
481 –    else
482 –      below_plane[0]=0;
483 –    /* Test if b[i] lies in or on triangle a */
484 –    test = test_point_against_spherical_tri(t0,t1,t2,p0,
485 –                                                 n,nset,&which,pt_sides[0]);
486 –    /* If pts[i] is interior: done */
487 –    if(!below_plane[0])
488 –      {
489 –        if(test == GT_INTERIOR)
490 –          return(TRUE);
491 –        /* Remember if b[i] fell on one of the 3 defining planes */
492 –        if(test)
493 –          on_edge++;
494 –      }
495 –    /* Now test point 1*/
496 –
497 –    if(DOT(avg,p1) < 0)
498 –      below_plane[1] = 1;
499 –    else
500 –      below_plane[1]=0;
501 –    /* Test if b[i] lies in or on triangle a */
502 –    test = test_point_against_spherical_tri(t0,t1,t2,p1,
503 –                                                 n,nset,&which,pt_sides[1]);
504 –    /* If pts[i] is interior: done */
505 –    if(!below_plane[1])
506 –    {
507 –      if(test == GT_INTERIOR)
508 –        return(TRUE);
509 –      /* Remember if b[i] fell on one of the 3 defining planes */
510 –      if(test)
511 –        on_edge++;
512 –    }
513 –    
514 –    /* Now test point 2 */
515 –    if(DOT(avg,p2) < 0)
516 –      below_plane[2] = 1;
517 –    else
518 –      below_plane[2]=0;
519 –        /* Test if b[i] lies in or on triangle a */
520 –    test = test_point_against_spherical_tri(t0,t1,t2,p2,
521 –                                                 n,nset,&which,pt_sides[2]);
522 –
523 –    /* If pts[i] is interior: done */
524 –    if(!below_plane[2])
525 –      {
526 –        if(test == GT_INTERIOR)
527 –          return(TRUE);
528 –        /* Remember if b[i] fell on one of the 3 defining planes */
529 –        if(test)
530 –          on_edge++;
531 –      }
532 –
533 –    /* If all three points below separating plane: trivial reject */
534 –    if(below_plane[0] && below_plane[1] && below_plane[2])
535 –       return(FALSE);
536 –    /* Accept if all points lie on a triangle vertex/edge edge- accept*/
537 –    if(on_edge == 3)
538 –       return(TRUE);
539 –    /* Now check vertices in a against triangle b */
540 –    return(FALSE);
541 – }
542 –
543 –
356   set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,n)
357     FVECT t0,t1,t2,p0,p1,p2;
358 <   char test[3];
359 <   char sides[3][3];
360 <   char nset;
358 >   int test[3];
359 >   int sides[3][3];
360 >   int nset;
361     FVECT n[3];
362   {
363 <    char t;
363 >    int t;
364      double d;
365  
366      
# Line 557 | Line 369 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
369      if(sides[0][0] == GT_INVALID)
370      {
371        if(!NTH_BIT(nset,0))
372 <        VCROSS(n[0],t1,t0);
372 >        VCROSS(n[0],t0,t1);
373        /* Test the point for sidedness */
374        d  = DOT(n[0],p0);
375 <      if(d >= 0)
375 >      if(d >= 0.0)
376          SET_NTH_BIT(test[0],0);
377      }
378      else
# Line 570 | Line 382 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
382      if(sides[0][1] == GT_INVALID)
383      {
384        if(!NTH_BIT(nset,1))
385 <        VCROSS(n[1],t2,t1);
385 >        VCROSS(n[1],t1,t2);
386          /* Test the point for sidedness */
387          d  = DOT(n[1],p0);
388 <        if(d >= 0)
388 >        if(d >= 0.0)
389            SET_NTH_BIT(test[0],1);
390      }
391      else
# Line 583 | Line 395 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
395      if(sides[0][2] == GT_INVALID)
396      {
397        if(!NTH_BIT(nset,2))
398 <        VCROSS(n[2],t0,t2);
398 >        VCROSS(n[2],t2,t0);
399        /* Test the point for sidedness */
400        d  = DOT(n[2],p0);
401 <      if(d >= 0)
401 >      if(d >= 0.0)
402          SET_NTH_BIT(test[0],2);
403      }
404      else
# Line 599 | Line 411 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
411      if(sides[1][0] == GT_INVALID)
412      {
413        if(!NTH_BIT(nset,0))
414 <        VCROSS(n[0],t1,t0);
414 >        VCROSS(n[0],t0,t1);
415        /* Test the point for sidedness */
416        d  = DOT(n[0],p1);
417 <      if(d >= 0)
417 >      if(d >= 0.0)
418          SET_NTH_BIT(test[1],0);
419      }
420      else
# Line 613 | Line 425 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
425      if(sides[1][1] == GT_INVALID)
426      {
427        if(!NTH_BIT(nset,1))
428 <        VCROSS(n[1],t2,t1);
428 >        VCROSS(n[1],t1,t2);
429        /* Test the point for sidedness */
430        d  = DOT(n[1],p1);
431 <      if(d >= 0)
431 >      if(d >= 0.0)
432          SET_NTH_BIT(test[1],1);
433      }
434      else
# Line 627 | Line 439 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
439      if(sides[1][2] == GT_INVALID)
440      {
441        if(!NTH_BIT(nset,2))
442 <        VCROSS(n[2],t0,t2);
442 >        VCROSS(n[2],t2,t0);
443        /* Test the point for sidedness */
444        d  = DOT(n[2],p1);
445 <      if(d >= 0)
445 >      if(d >= 0.0)
446          SET_NTH_BIT(test[1],2);
447      }
448      else
# Line 643 | Line 455 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
455      if(sides[2][0] == GT_INVALID)
456      {
457        if(!NTH_BIT(nset,0))
458 <        VCROSS(n[0],t1,t0);
458 >        VCROSS(n[0],t0,t1);
459        /* Test the point for sidedness */
460        d  = DOT(n[0],p2);
461 <      if(d >= 0)
461 >      if(d >= 0.0)
462          SET_NTH_BIT(test[2],0);
463      }
464      else
# Line 656 | Line 468 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
468      if(sides[2][1] == GT_INVALID)
469      {
470        if(!NTH_BIT(nset,1))
471 <        VCROSS(n[1],t2,t1);
471 >        VCROSS(n[1],t1,t2);
472        /* Test the point for sidedness */
473        d  = DOT(n[1],p2);
474 <      if(d >= 0)
474 >      if(d >= 0.0)
475          SET_NTH_BIT(test[2],1);
476      }
477      else
# Line 669 | Line 481 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
481      if(sides[2][2] == GT_INVALID)
482      {
483        if(!NTH_BIT(nset,2))
484 <        VCROSS(n[2],t0,t2);
484 >        VCROSS(n[2],t2,t0);
485        /* Test the point for sidedness */
486        d  = DOT(n[2],p2);
487 <      if(d >= 0)
487 >      if(d >= 0.0)
488          SET_NTH_BIT(test[2],2);
489      }
490      else
# Line 680 | Line 492 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
492          SET_NTH_BIT(test[2],2);
493   }
494  
495 + double
496 + point_on_sphere(ps,p,c)
497 + FVECT ps,p,c;
498 + {
499 +  double d;
500 +    VSUB(ps,p,c);    
501 +    d= normalize(ps);
502 +    return(d);
503 + }
504  
505   int
506 < spherical_tri_intersect(a1,a2,a3,b1,b2,b3)
507 < FVECT a1,a2,a3,b1,b2,b3;
506 > point_in_stri(v0,v1,v2,p)
507 > FVECT v0,v1,v2,p;
508   {
509 <  char which,test,n_set[2];
510 <  char sides[2][3][3],i,j,inext,jnext;
690 <  char tests[2][3];
691 <  FVECT n[2][3],p,avg[2];
692 <
693 <  /* Test the vertices of triangle a against the pyramid formed by triangle
694 <     b and the origin. If any vertex of a is interior to triangle b, or
695 <     if all 3 vertices of a are ON the edges of b,return TRUE. Remember
696 <     the results of the edge normal and sidedness tests for later.
697 <   */
698 < if(test_vertices_for_tri_inclusion(a1,a2,a3,b1,b2,b3,
699 <                                    &(n_set[0]),n[0],avg[0],sides[1]))
700 <     return(TRUE);
701 <  
702 < if(test_vertices_for_tri_inclusion(b1,b2,b3,a1,a2,a3,
703 <                                    &(n_set[1]),n[1],avg[1],sides[0]))
704 <     return(TRUE);
509 >    double d;
510 >    FVECT n;  
511  
512 +    VCROSS(n,v0,v1);
513 +    /* Test the point for sidedness */
514 +    d  = DOT(n,p);
515 +    if(d > 0.0)
516 +      return(FALSE);
517  
518 <  set_sidedness_tests(b1,b2,b3,a1,a2,a3,tests[0],sides[0],n_set[1],n[1]);
519 <  if(tests[0][0]&tests[0][1]&tests[0][2])
520 <    return(FALSE);
518 >    /* Test next edge */
519 >    VCROSS(n,v1,v2);
520 >    /* Test the point for sidedness */
521 >    d  = DOT(n,p);
522 >    if(d > 0.0)
523 >       return(FALSE);
524  
525 <  set_sidedness_tests(a1,a2,a3,b1,b2,b3,tests[1],sides[1],n_set[0],n[0]);
526 <  if(tests[1][0]&tests[1][1]&tests[1][2])
527 <    return(FALSE);
528 <
529 <  for(j=0; j < 3;j++)
530 <  {
531 <    jnext = (j+1)%3;
532 <    /* IF edge b doesnt cross any great circles of a, punt */
719 <    if(tests[1][j] & tests[1][jnext])
720 <      continue;
721 <    for(i=0;i<3;i++)
722 <    {
723 <      inext = (i+1)%3;
724 <      /* IF edge a doesnt cross any great circles of b, punt */
725 <      if(tests[0][i] & tests[0][inext])
726 <        continue;
727 <      /* Now find the great circles that cross and test */
728 <      if((NTH_BIT(tests[0][i],j)^(NTH_BIT(tests[0][inext],j)))
729 <          && (NTH_BIT(tests[1][j],i)^NTH_BIT(tests[1][jnext],i)))
730 <      {
731 <        VCROSS(p,n[0][i],n[1][j]);
732 <                    
733 <        /* If zero cp= done */
734 <        if(ZERO_VEC3(p))
735 <          continue;
736 <        /* check above both planes */
737 <        if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0)
738 <          {
739 <            NEGATE_VEC3(p);
740 <            if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0)
741 <              continue;
742 <          }
743 <        return(TRUE);
744 <      }
745 <    }
746 <  }
747 <  return(FALSE);
525 >    /* Test next edge */
526 >    VCROSS(n,v2,v0);
527 >    /* Test the point for sidedness */
528 >    d  = DOT(n,p);
529 >    if(d > 0.0)
530 >       return(FALSE);
531 >    /* Must be interior to the pyramid */
532 >    return(GT_INTERIOR);
533   }
534  
535 +
536   int
537 < ray_intersect_tri(orig,dir,v0,v1,v2,pt,wptr)
537 > ray_intersect_tri(orig,dir,v0,v1,v2,pt)
538   FVECT orig,dir;
539   FVECT v0,v1,v2;
540   FVECT pt;
755 – char *wptr;
541   {
542 <  FVECT p0,p1,p2,p,n;
543 <  char type,which;
544 <  double pd;
760 <  
761 <  point_on_sphere(p0,v0,orig);
762 <  point_on_sphere(p1,v1,orig);
763 <  point_on_sphere(p2,v2,orig);
764 <  type = test_single_point_against_spherical_tri(p0,p1,p2,dir,&which);
542 >  FVECT p0,p1,p2,p;
543 >  FPEQ peq;
544 >  int type;
545  
546 <  if(type)
546 >  VSUB(p0,v0,orig);
547 >  VSUB(p1,v1,orig);
548 >  VSUB(p2,v2,orig);
549 >
550 >  if(point_in_stri(p0,p1,p2,dir))
551    {
552        /* Intersect the ray with the triangle plane */
553 <      tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
554 <      intersect_ray_plane(orig,dir,n,pd,NULL,pt);        
553 >      tri_plane_equation(v0,v1,v2,&peq,FALSE);
554 >      return(intersect_ray_plane(orig,dir,peq,NULL,pt));
555    }
556 <  if(wptr)
773 <    *wptr = which;
774 <
775 <  return(type);
556 >  return(FALSE);
557   }
558  
559  
# Line 831 | Line 612 | FVECT fnear[4],ffar[4];
612      ffar[3][2] =  width*nhv[2] - height*nvv[2] + far*ndv[2] + vp[2] ;
613   }
614  
615 + int
616 + max_index(v,r)
617 + FVECT v;
618 + double *r;
619 + {
620 +  double p[3];
621 +  int i;
622  
623 +  p[0] = fabs(v[0]);
624 +  p[1] = fabs(v[1]);
625 +  p[2] = fabs(v[2]);
626 +  i = (p[0]>=p[1])?((p[0]>=p[2])?0:2):((p[1]>=p[2])?1:2);  
627 +  if(r)
628 +    *r = p[i];
629 +  return(i);
630 + }
631  
836 –
632   int
633 < spherical_polygon_edge_intersect(a0,a1,b0,b1)
634 < FVECT a0,a1,b0,b1;
633 > closest_point_in_tri(p0,p1,p2,p,p0id,p1id,p2id)
634 > FVECT p0,p1,p2,p;
635 > int p0id,p1id,p2id;
636   {
637 <    FVECT na,nb,avga,avgb,p;
638 <    double d;
639 <    int sb0,sb1,sa0,sa1;
637 >    double d,d1;
638 >    int i;
639 >    
640 >    d =  DIST_SQ(p,p0);
641 >    d1 = DIST_SQ(p,p1);
642 >    if(d < d1)
643 >    {
644 >      d1 = DIST_SQ(p,p2);
645 >      i = (d1 < d)?p2id:p0id;
646 >    }
647 >    else
648 >    {
649 >      d = DIST_SQ(p,p2);
650 >      i = (d < d1)? p2id:p1id;
651 >    }
652 >    return(i);
653 > }
654  
655 <    /* First test if edge b straddles great circle of a */
656 <    VCROSS(na,a0,a1);
657 <    d = DOT(na,b0);
658 <    sb0 = ZERO(d)?0:(d<0)? -1: 1;
659 <    d = DOT(na,b1);
660 <    sb1 = ZERO(d)?0:(d<0)? -1: 1;
661 <    /* edge b entirely on one side of great circle a: edges cannot intersect*/
662 <    if(sb0*sb1 > 0)
663 <       return(FALSE);
854 <    /* test if edge a straddles great circle of b */
855 <    VCROSS(nb,b0,b1);
856 <    d = DOT(nb,a0);
857 <    sa0 = ZERO(d)?0:(d<0)? -1: 1;
858 <    d = DOT(nb,a1);
859 <    sa1 = ZERO(d)?0:(d<0)? -1: 1;
860 <    /* edge a entirely on one side of great circle b: edges cannot intersect*/
861 <    if(sa0*sa1 > 0)
862 <       return(FALSE);
655 > /* Find the normalized barycentric coordinates of p relative to
656 > * triangle v0,v1,v2. Return result in coord
657 > */
658 > bary2d(x1,y1,x2,y2,x3,y3,px,py,coord)
659 > double x1,y1,x2,y2,x3,y3;
660 > double px,py;
661 > double coord[3];
662 > {
663 >  double a;
664  
665 <    /* Find one of intersection points of the great circles */
666 <    VCROSS(p,na,nb);
667 <    /* If they lie on same great circle: call an intersection */
668 <    if(ZERO_VEC3(p))
669 <       return(TRUE);
665 >  a =  (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1);
666 >  coord[0] = ((x2 - px) * (y3 - py) - (x3 - px) * (y2 - py)) / a;
667 >  coord[1] = ((x3 - px) * (y1 - py) - (x1 - px) * (y3 - py)) / a;
668 >  coord[2] = ((x1 - px) * (y2 - py) - (x2 - px) * (y1 - py)) / a;
669 >
670 > }
671  
672 <    VADD(avga,a0,a1);
673 <    VADD(avgb,b0,b1);
674 <    if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
672 >
673 >
674 >
675 > bary_parent(coord,i)
676 > BCOORD coord[3];
677 > int i;
678 > {
679 >  switch(i) {
680 >  case 0:
681 >    /* update bary for child */
682 >    coord[0] = (coord[0] >> 1) + MAXBCOORD4;
683 >    coord[1] >>= 1;
684 >    coord[2] >>= 1;
685 >    break;
686 >  case 1:
687 >    coord[0] >>= 1;
688 >    coord[1]  = (coord[1] >> 1) + MAXBCOORD4;
689 >    coord[2] >>= 1;
690 >    break;
691 >    
692 >  case 2:
693 >    coord[0] >>= 1;
694 >    coord[1] >>= 1;
695 >    coord[2] = (coord[2] >> 1) + MAXBCOORD4;
696 >    break;
697 >    
698 >  case 3:
699 >    coord[0] = MAXBCOORD4 - (coord[0] >> 1);
700 >    coord[1] = MAXBCOORD4 - (coord[1] >> 1);
701 >    coord[2] = MAXBCOORD4 - (coord[2] >> 1);
702 >    break;
703 > #ifdef DEBUG
704 >  default:
705 >    eputs("bary_parent():Invalid child\n");
706 >    break;
707 > #endif
708 >  }
709 > }
710 >
711 > bary_from_child(coord,child,next)
712 > BCOORD coord[3];
713 > int child,next;
714 > {
715 > #ifdef DEBUG
716 >  if(child <0 || child > 3)
717 >  {
718 >    eputs("bary_from_child():Invalid child\n");
719 >    return;
720 >  }
721 >  if(next <0 || next > 3)
722 >  {
723 >    eputs("bary_from_child():Invalid next\n");
724 >    return;
725 >  }
726 > #endif
727 >  if(next == child)
728 >    return;
729 >
730 >  switch(child){
731 >  case 0:
732 >      coord[0] = 0;
733 >      coord[1] = MAXBCOORD2 - coord[1];
734 >      coord[2] = MAXBCOORD2 - coord[2];
735 >      break;
736 >  case 1:
737 >      coord[0] = MAXBCOORD2 - coord[0];
738 >      coord[1] = 0;
739 >      coord[2] = MAXBCOORD2 - coord[2];
740 >      break;
741 >  case 2:
742 >      coord[0] = MAXBCOORD2 - coord[0];
743 >      coord[1] = MAXBCOORD2 - coord[1];
744 >      coord[2] = 0;
745 >    break;
746 >  case 3:
747 >    switch(next){
748 >    case 0:
749 >      coord[0] = 0;
750 >      coord[1] = MAXBCOORD2 - coord[1];
751 >      coord[2] = MAXBCOORD2 - coord[2];
752 >      break;
753 >    case 1:
754 >      coord[0] = MAXBCOORD2 - coord[0];
755 >      coord[1] = 0;
756 >      coord[2] = MAXBCOORD2 - coord[2];
757 >      break;
758 >    case 2:
759 >      coord[0] = MAXBCOORD2 - coord[0];
760 >      coord[1] = MAXBCOORD2 - coord[1];
761 >      coord[2] = 0;
762 >      break;
763 >    }
764 >    break;
765 >  }
766 > }
767 >
768 > int
769 > bary_child(coord)
770 > BCOORD coord[3];
771 > {
772 >  int i;
773 >
774 >  if(coord[0] > MAXBCOORD4)
775 >  {
776 >      /* update bary for child */
777 >      coord[0] = (coord[0]<< 1) - MAXBCOORD2;
778 >      coord[1] <<= 1;
779 >      coord[2] <<= 1;
780 >      return(0);
781 >  }
782 >  else
783 >    if(coord[1] > MAXBCOORD4)
784      {
785 <      NEGATE_VEC3(p);
786 <      if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
787 <        return(FALSE);
785 >      coord[0] <<= 1;
786 >      coord[1] = (coord[1] << 1) - MAXBCOORD2;
787 >      coord[2] <<= 1;
788 >      return(1);
789      }
790 <    if((!sb0 || !sb1) && (!sa0 || !sa1))
791 <      return(FALSE);
792 <    return(TRUE);
790 >    else
791 >      if(coord[2] > MAXBCOORD4)
792 >      {
793 >        coord[0] <<= 1;
794 >        coord[1] <<= 1;
795 >        coord[2] = (coord[2] << 1) - MAXBCOORD2;
796 >        return(2);
797 >      }
798 >      else
799 >         {
800 >           coord[0] = MAXBCOORD2 - (coord[0] << 1);
801 >           coord[1] = MAXBCOORD2 - (coord[1] << 1);
802 >           coord[2] = MAXBCOORD2 - (coord[2] << 1);
803 >           return(3);
804 >         }
805   }
806 +
807 + int
808 + bary_nth_child(coord,i)
809 + BCOORD coord[3];
810 + int i;
811 + {
812 +
813 +  switch(i){
814 +  case 0:
815 +    /* update bary for child */
816 +    coord[0] = (coord[0]<< 1) - MAXBCOORD2;
817 +    coord[1] <<= 1;
818 +    coord[2] <<= 1;
819 +    break;
820 +  case 1:
821 +    coord[0] <<= 1;
822 +    coord[1] = (coord[1] << 1) - MAXBCOORD2;
823 +    coord[2] <<= 1;
824 +    break;
825 +  case 2:
826 +    coord[0] <<= 1;
827 +    coord[1] <<= 1;
828 +    coord[2] = (coord[2] << 1) - MAXBCOORD2;
829 +    break;
830 +  case 3:
831 +    coord[0] = MAXBCOORD2 - (coord[0] << 1);
832 +    coord[1] = MAXBCOORD2 - (coord[1] << 1);
833 +    coord[2] = MAXBCOORD2 - (coord[2] << 1);
834 +    break;
835 +  }
836 + }
837 +
838 +
839  

Diff Legend

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