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.11 by gwlarson, Sun Jan 10 10:27:47 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 <    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)
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   */
64   double
65   tri_normal(v0,v1,v2,n,norm)
66   FVECT v0,v1,v2,n;
67 < char norm;
67 > int norm;
68   {
69    double mag;
70  
# Line 63 | Line 74 | char 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 | char norm;
82  
83    if(!norm)
84       return(0);
75 –
85    
86    mag = normalize(n);
87  
# Line 80 | Line 89 | char norm;
89   }
90  
91  
92 < tri_plane_equation(v0,v1,v2,n,nd,norm)
93 <   FVECT v0,v1,v2,n;
94 <   double *nd;
95 <   char norm;
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  
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 – }
108 –
102   /* From quad_edge-code */
103   int
104   point_in_circle_thru_origin(p,p0,p1)
# Line 122 | Line 115 | FVECT p0,p1;
115  
116      det = -dp0*CROSS_VEC2(p1,p) + dp1*CROSS_VEC2(p0,p) - dp*CROSS_VEC2(p0,p1);
117      
118 <    return (det > 0);
118 >    return (det > (1e-10));
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  
133 + /* returns TRUE if ray from origin in direction v intersects plane defined
134 +  * by normal plane_n, and plane_d. If plane is not parallel- returns
135 +  * intersection point if r != NULL. If tptr!= NULL returns value of
136 +  * t, if parallel, returns t=FHUGE
137 +  */
138   int
139 < intersect_vector_plane(v,plane_n,plane_d,pd,r)
140 <   FVECT v,plane_n;
141 <   double plane_d;
142 <   double *pd;
139 > intersect_vector_plane(v,peq,tptr,r)
140 >   FVECT v;
141 >   FPEQ peq;
142 >   double *tptr;
143     FVECT r;
144   {
145 <  double t;
145 >  double t,d;
146    int hit;
147      /*
148        Plane is Ax + By + Cz +D = 0:
# Line 152 | Line 152 | intersect_vector_plane(v,plane_n,plane_d,pd,r)
152      /* line is  l = p1 + (p2-p1)t, p1=origin */
153  
154      /* Solve for t: */
155 <    t =  plane_d/-(DOT(plane_n,v));
156 <    if(t >0 || ZERO(t))
157 <       hit = 1;
158 <    else
155 >  d = -(DOT(FP_N(peq),v));
156 >  if(ZERO(d))
157 >  {
158 >      t = FHUGE;
159 >      hit = 0;
160 >  }
161 >  else
162 >  {
163 >      t =  FP_D(peq)/d;
164 >      if(t < 0 )
165 >         hit = 0;
166 >      else
167 >         hit = 1;
168 >      if(r)
169 >         {
170 >             r[0] = v[0]*t;
171 >             r[1] = v[1]*t;
172 >             r[2] = v[2]*t;
173 >         }
174 >  }
175 >    if(tptr)
176 >       *tptr = t;
177 >  return(hit);
178 > }
179 >
180 > int
181 > intersect_ray_plane(orig,dir,peq,pd,r)
182 >   FVECT orig,dir;
183 >   FPEQ peq;
184 >   double *pd;
185 >   FVECT r;
186 > {
187 >  double t,d;
188 >  int hit;
189 >    /*
190 >      Plane is Ax + By + Cz +D = 0:
191 >      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
192 >    */
193 >     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
194 >         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
195 >       line is  l = p1 + (p2-p1)t
196 >     */
197 >    /* Solve for t: */
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 <    r[0] = v[0]*t;
206 <    r[1] = v[1]*t;
207 <    r[2] = v[2]*t;
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_plane(orig,dir,plane_n,plane_d,pd,r)
218 > intersect_ray_oplane(orig,dir,n,pd,r)
219     FVECT orig,dir;
220 <   FVECT plane_n;
172 <   double plane_d;
220 >   FVECT n;
221     double *pd;
222     FVECT r;
223   {
224 <  double t;
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 */
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 <    t =  -(DOT(plane_n,orig) + plane_d)/(DOT(plane_n,dir));
236 <    if(ZERO(t) || t >0)
237 <       hit = 1;
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
242 +       hit = 1;
243 +
244 +  if(r)
245 +     VSUM(r,orig,dir,t);
246 +
247 +    if(pd)
248 +       *pd = t;
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,peq,pd,r)
273 +   FVECT e0,e1;
274 +   FPEQ peq;
275 +   double *pd;
276 +   FVECT r;
277 + {
278 +  double t;
279 +  int hit;
280 +  FVECT d;
281 +  /*
282 +      Plane is Ax + By + Cz +D = 0:
283 +      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
284 +    */
285 +     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
286 +         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
287 +       line is  l = p1 + (p2-p1)t
288 +     */
289 +    /* Solve for t: */
290 +  VSUB(d,e1,e0);
291 +  t =  -(DOT(FP_N(peq),e0) + FP_D(peq))/(DOT(FP_N(peq),d));
292 +    if(t < 0)
293         hit = 0;
294 +    else
295 +       hit = 1;
296  
297 <  VSUM(r,orig,dir,t);
297 >  VSUM(r,e0,d,t);
298  
299      if(pd)
300         *pd = t;
# Line 202 | Line 307 | point_in_cone(p,p0,p1,p2)
307   FVECT p;
308   FVECT p0,p1,p2;
309   {
205 –    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
324         n cross x-axis
325       */
326      /* Project p onto the plane */
327 <    if(!intersect_vector_plane(p,n,d,NULL,np))
327 >    /* NOTE: check this: does sideness check?*/
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 251 | Line 357 | FVECT p0,p1,p2;
357   }
358  
359   int
360 < test_point_against_spherical_tri(v0,v1,v2,p,n,nset,which,sides)
360 > point_set_in_stri(v0,v1,v2,p,n,nset,sides)
361   FVECT v0,v1,v2,p;
362   FVECT n[3];
363 < char *nset;
364 < char *which;
259 < char sides[3];
363 > int *nset;
364 > int sides[3];
365  
366   {
367 <    float d;
263 <
367 >    double d;
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 */
375      d  = DOT(n[0],p);
376  
377 <    if(ZERO(d))
378 <       sides[0] = GT_EDGE;
379 <    else
380 <       if(d > 0)
381 <      {
278 <          sides[0] =  GT_OUT;
279 <          sides[1] = sides[2] = GT_INVALID;
280 <          return(FALSE);
377 >    if(d > 0.0)
378 >     {
379 >       sides[0] =  GT_OUT;
380 >       sides[1] = sides[2] = GT_INVALID;
381 >       return(FALSE);
382        }
383      else
384         sides[0] = GT_INTERIOR;
# Line 285 | Line 386 | char 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 */
393      d  = DOT(n[1],p);
394 <    if(ZERO(d))
394 >    if(d > 0.0)
395      {
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 –    {
396          sides[1] = GT_OUT;
397          sides[2] = GT_INVALID;
398          return(FALSE);
# Line 312 | Line 402 | char sides[3];
402      /* Test next edge */
403      if(!NTH_BIT(*nset,2))
404      {
405 <
316 <        VCROSS(n[2],v0,v2);
405 >        VCROSS(n[2],v2,v0);
406          SET_NTH_BIT(*nset,2);
407      }
408      /* Test the point for sidedness */
409      d  = DOT(n[2],p);
410 <    if(ZERO(d))
410 >    if(d > 0.0)
411      {
412 <        sides[2] = GT_EDGE;
413 <
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 <           }
412 >      sides[2] = GT_OUT;
413 >      return(FALSE);
414      }
344 –    else if(d > 0)
345 –      {
346 –        sides[2] = GT_OUT;
347 –        return(FALSE);
348 –      }
349 –    /* If on edge */
415      else
416         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 –    }
417      /* Must be interior to the pyramid */
418      return(GT_INTERIOR);
419   }
420  
421  
422  
423 <
423 >
424   int
425 < test_single_point_against_spherical_tri(v0,v1,v2,p,which)
425 > point_in_stri(v0,v1,v2,p)
426   FVECT v0,v1,v2,p;
375 – char *which;
427   {
428 <    float d;
428 >    double d;
429      FVECT n;  
379 –    char sides[3];
430  
431 <    /* 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);
431 >    VCROSS(n,v0,v1);
432      /* Test the point for sidedness */
433      d  = DOT(n,p);
434 <    if(ZERO(d))
435 <       sides[0] = GT_EDGE;
436 <    else
403 <       if(d > 0)
404 <          return(FALSE);
405 <       else
406 <          sides[0] = GT_INTERIOR;
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(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)
441 >    if(d > 0.0)
442         return(FALSE);
423 –    else
424 –       sides[1] = GT_INTERIOR;
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(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)
448 >    if(d > 0.0)
449         return(FALSE);
450      /* Must be interior to the pyramid */
451 <    return(GT_FACE);
451 >    return(GT_INTERIOR);
452   }
453  
454   int
455 < test_vertices_for_tri_inclusion(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides)
455 > vertices_in_stri(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides)
456   FVECT t0,t1,t2,p0,p1,p2;
457 < char *nset;
457 > int *nset;
458   FVECT n[3];
459   FVECT avg;
460 < char pt_sides[3][3];
460 > int pt_sides[3][3];
461  
462   {
463 <    char below_plane[3],on_edge,test;
468 <    char which;
463 >    int below_plane[3],test;
464  
465      SUM_3VEC3(avg,t0,t1,t2);
471 –    on_edge = 0;
466      *nset = 0;
467      /* Test vertex v[i] against triangle j*/
468      /* Check if v[i] lies below plane defined by avg of 3 vectors
# Line 476 | Line 470 | char pt_sides[3][3];
470         */
471  
472      /* test point 0 */
473 <    if(DOT(avg,p0) < 0)
473 >    if(DOT(avg,p0) < 0.0)
474        below_plane[0] = 1;
475      else
476 <      below_plane[0]=0;
476 >      below_plane[0] = 0;
477      /* Test if b[i] lies in or on triangle a */
478 <    test = test_point_against_spherical_tri(t0,t1,t2,p0,
485 <                                                 n,nset,&which,pt_sides[0]);
478 >    test = point_set_in_stri(t0,t1,t2,p0,n,nset,pt_sides[0]);
479      /* If pts[i] is interior: done */
480      if(!below_plane[0])
481        {
482          if(test == GT_INTERIOR)
483            return(TRUE);
491 –        /* Remember if b[i] fell on one of the 3 defining planes */
492 –        if(test)
493 –          on_edge++;
484        }
485      /* Now test point 1*/
486  
487 <    if(DOT(avg,p1) < 0)
487 >    if(DOT(avg,p1) < 0.0)
488        below_plane[1] = 1;
489      else
490        below_plane[1]=0;
491      /* Test if b[i] lies in or on triangle a */
492 <    test = test_point_against_spherical_tri(t0,t1,t2,p1,
503 <                                                 n,nset,&which,pt_sides[1]);
492 >    test = point_set_in_stri(t0,t1,t2,p1,n,nset,pt_sides[1]);
493      /* If pts[i] is interior: done */
494      if(!below_plane[1])
495      {
496        if(test == GT_INTERIOR)
497          return(TRUE);
509 –      /* Remember if b[i] fell on one of the 3 defining planes */
510 –      if(test)
511 –        on_edge++;
498      }
499      
500      /* Now test point 2 */
501 <    if(DOT(avg,p2) < 0)
501 >    if(DOT(avg,p2) < 0.0)
502        below_plane[2] = 1;
503      else
504 <      below_plane[2]=0;
504 >      below_plane[2] = 0;
505          /* Test if b[i] lies in or on triangle a */
506 <    test = test_point_against_spherical_tri(t0,t1,t2,p2,
521 <                                                 n,nset,&which,pt_sides[2]);
506 >    test = point_set_in_stri(t0,t1,t2,p2,n,nset,pt_sides[2]);
507  
508      /* If pts[i] is interior: done */
509      if(!below_plane[2])
510        {
511          if(test == GT_INTERIOR)
512            return(TRUE);
528 –        /* Remember if b[i] fell on one of the 3 defining planes */
529 –        if(test)
530 –          on_edge++;
513        }
514  
515      /* If all three points below separating plane: trivial reject */
516      if(below_plane[0] && below_plane[1] && below_plane[2])
517         return(FALSE);
536 –    /* Accept if all points lie on a triangle vertex/edge edge- accept*/
537 –    if(on_edge == 3)
538 –       return(TRUE);
518      /* Now check vertices in a against triangle b */
519      return(FALSE);
520   }
# Line 543 | Line 522 | char pt_sides[3][3];
522  
523   set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,n)
524     FVECT t0,t1,t2,p0,p1,p2;
525 <   char test[3];
526 <   char sides[3][3];
527 <   char nset;
525 >   int test[3];
526 >   int sides[3][3];
527 >   int nset;
528     FVECT n[3];
529   {
530 <    char t;
530 >    int t;
531      double d;
532  
533      
# Line 557 | 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)
542 >      if(d >= 0.0)
543          SET_NTH_BIT(test[0],0);
544      }
545      else
# Line 570 | 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)
555 >        if(d >= 0.0)
556            SET_NTH_BIT(test[0],1);
557      }
558      else
# Line 583 | 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)
568 >      if(d >= 0.0)
569          SET_NTH_BIT(test[0],2);
570      }
571      else
# Line 599 | 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)
584 >      if(d >= 0.0)
585          SET_NTH_BIT(test[1],0);
586      }
587      else
# Line 613 | 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)
598 >      if(d >= 0.0)
599          SET_NTH_BIT(test[1],1);
600      }
601      else
# Line 627 | 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)
612 >      if(d >= 0.0)
613          SET_NTH_BIT(test[1],2);
614      }
615      else
# Line 643 | 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)
628 >      if(d >= 0.0)
629          SET_NTH_BIT(test[2],0);
630      }
631      else
# Line 656 | 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)
641 >      if(d >= 0.0)
642          SET_NTH_BIT(test[2],1);
643      }
644      else
# Line 669 | 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)
654 >      if(d >= 0.0)
655          SET_NTH_BIT(test[2],2);
656      }
657      else
# Line 682 | Line 661 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
661  
662  
663   int
664 < spherical_tri_intersect(a1,a2,a3,b1,b2,b3)
664 > stri_intersect(a1,a2,a3,b1,b2,b3)
665   FVECT a1,a2,a3,b1,b2,b3;
666   {
667 <  char which,test,n_set[2];
668 <  char sides[2][3][3],i,j,inext,jnext;
669 <  char tests[2][3];
670 <  FVECT n[2][3],p,avg[2];
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],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(test_vertices_for_tri_inclusion(a1,a2,a3,b1,b2,b3,
699 <                                    &(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(test_vertices_for_tri_inclusion(b1,b2,b3,a1,a2,a3,
703 <                                    &(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 748 | Line 741 | FVECT a1,a2,a3,b1,b2,b3;
741   }
742  
743   int
744 < ray_intersect_tri(orig,dir,v0,v1,v2,pt,wptr)
744 > ray_intersect_tri(orig,dir,v0,v1,v2,pt)
745   FVECT orig,dir;
746   FVECT v0,v1,v2;
747   FVECT pt;
755 – char *wptr;
748   {
749 <  FVECT p0,p1,p2,p,n;
750 <  char type,which;
751 <  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);
749 >  FVECT p0,p1,p2,p;
750 >  FPEQ peq;
751 >  int type;
752  
753 <  if(type)
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 <      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 <  if(wptr)
773 <    *wptr = which;
774 <
775 <  return(type);
763 >  return(FALSE);
764   }
765  
766  
# Line 831 | 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 < spherical_polygon_edge_intersect(a0,a1,b0,b1)
864 > sedge_intersect(a0,a1,b0,b1)
865   FVECT a0,a1,b0,b1;
866   {
867      FVECT na,nb,avga,avgb,p;
# Line 879 | Line 905 | FVECT a0,a1,b0,b1;
905        return(FALSE);
906      return(TRUE);
907   }
908 +
909 +
910 + /* Find the normalized barycentric coordinates of p relative to
911 + * triangle v0,v1,v2. Return result in coord
912 + */
913 + bary2d(x1,y1,x2,y2,x3,y3,px,py,coord)
914 + double x1,y1,x2,y2,x3,y3;
915 + double px,py;
916 + double coord[3];
917 + {
918 +  double a;
919 +
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] = ((x1 - px) * (y2 - py) - (x2 - px) * (y1 - py)) / a;
924 +
925 + }
926 +
927 +
928 +
929 +
930 + bary_parent(coord,i)
931 + BCOORD coord[3];
932 + int i;
933 + {
934 +  switch(i) {
935 +  case 0:
936 +    /* update bary for child */
937 +    coord[0] = (coord[0] >> 1) + MAXBCOORD4;
938 +    coord[1] >>= 1;
939 +    coord[2] >>= 1;
940 +    break;
941 +  case 1:
942 +    coord[0] >>= 1;
943 +    coord[1]  = (coord[1] >> 1) + MAXBCOORD4;
944 +    coord[2] >>= 1;
945 +    break;
946 +    
947 +  case 2:
948 +    coord[0] >>= 1;
949 +    coord[1] >>= 1;
950 +    coord[2] = (coord[2] >> 1) + MAXBCOORD4;
951 +    break;
952 +    
953 +  case 3:
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:
960 +    eputs("bary_parent():Invalid child\n");
961 +    break;
962 + #endif
963 +  }
964 + }
965 +
966 + bary_from_child(coord,child,next)
967 + BCOORD coord[3];
968 + int child,next;
969 + {
970 + #ifdef DEBUG
971 +  if(child <0 || child > 3)
972 +  {
973 +    eputs("bary_from_child():Invalid child\n");
974 +    return;
975 +  }
976 +  if(next <0 || next > 3)
977 +  {
978 +    eputs("bary_from_child():Invalid next\n");
979 +    return;
980 +  }
981 + #endif
982 +  if(next == child)
983 +    return;
984 +
985 +  switch(child){
986 +  case 0:
987 +      coord[0] = 0;
988 +      coord[1] = MAXBCOORD2 - coord[1];
989 +      coord[2] = MAXBCOORD2 - coord[2];
990 +      break;
991 +  case 1:
992 +      coord[0] = MAXBCOORD2 - coord[0];
993 +      coord[1] = 0;
994 +      coord[2] = MAXBCOORD2 - coord[2];
995 +      break;
996 +  case 2:
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] = 0;
1005 +      coord[1] = MAXBCOORD2 - coord[1];
1006 +      coord[2] = MAXBCOORD2 - coord[2];
1007 +      break;
1008 +    case 1:
1009 +      coord[0] = MAXBCOORD2 - coord[0];
1010 +      coord[1] = 0;
1011 +      coord[2] = MAXBCOORD2 - coord[2];
1012 +      break;
1013 +    case 2:
1014 +      coord[0] = MAXBCOORD2 - coord[0];
1015 +      coord[1] = MAXBCOORD2 - coord[1];
1016 +      coord[2] = 0;
1017 +      break;
1018 +    }
1019 +    break;
1020 +  }
1021 + }
1022 +
1023 + int
1024 + bary_child(coord)
1025 + BCOORD coord[3];
1026 + {
1027 +  int i;
1028 +
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 +    if(coord[1] > MAXBCOORD4)
1039 +    {
1040 +      coord[0] <<= 1;
1041 +      coord[1] = (coord[1] << 1) - MAXBCOORD2;
1042 +      coord[2] <<= 1;
1043 +      return(1);
1044 +    }
1045 +    else
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 +  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