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.2 by gwlarson, Thu Aug 20 16:47:21 1998 UTC vs.
Revision 3.13 by greg, Sat Feb 22 02:07:25 2003 UTC

# Line 1 | Line 1
1 – /* Copyright (c) 1998 Silicon Graphics, Inc. */
2 –
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ SGI";
2 > static const char       RCSid[] = "$Id$";
3   #endif
6 –
4   /*
5   *  sm_geom.c
6 + *    some geometric utility routines
7   */
8  
9   #include "standard.h"
10   #include "sm_geom.h"
11  
12 < tri_centroid(v0,v1,v2,c)
13 < FVECT v0,v1,v2,c;
12 > /*
13 > * int
14 > * pt_in_cone(p,a,b,c)
15 > *             : test if point p lies in cone defined by a,b,c and origin
16 > * double p[3];              : point to test
17 > * double a[3],b[3],c[3];    : points forming triangle
18 > *
19 > *  Assumes apex at origin, a,b,c are unit vectors defining the
20 > *  triangle which the cone circumscribes. Assumes p is also normalized
21 > *  Test is implemented as:
22 > *             r =  (b-a)X(c-a)
23 > *             in = (p.r) > (a.r)
24 > *  The center of the cone is r, and corresponds to the triangle normal.
25 > *  p.r is the proportional to the cosine of the angle between p and the
26 > *  the cone center ray, and a.r to the radius of the cone. If the cosine
27 > *  of the angle for p is greater than that for a, the angle between p
28 > *  and r is smaller, and p lies in the cone.
29 > */
30 > int
31 > pt_in_cone(p,a,b,c)
32 > double p[3],a[3],b[3],c[3];
33   {
34 < /* Average three triangle vertices to give centroid: return in c */
35 <  c[0] = (v0[0] + v1[0] + v2[0])/3.0;
36 <  c[1] = (v0[1] + v1[1] + v2[1])/3.0;
37 <  c[2] = (v0[2] + v1[2] + v2[2])/3.0;
34 >  double r[3];
35 >  double pr,ar;
36 >  double ab[3],ac[3];
37 >
38 > #ifdef DEBUG
39 > #if DEBUG > 1
40 > {
41 >  double l;
42 >  VSUB(ab,b,a);
43 >  normalize(ab);
44 >  VSUB(ac,c,a);
45 >  normalize(ac);
46 >  VCROSS(r,ab,ac);
47 >  l = normalize(r);
48 >  /* l = sin@ between ab,ac - if 0 vectors are colinear */
49 >  if( l <= COLINEAR_EPS)
50 >  {
51 >    eputs("pt in cone: null triangle:returning FALSE\n");
52 >    return(FALSE);
53 >  }
54   }
55 + #endif
56 + #endif
57  
58 +  VSUB(ab,b,a);
59 +  VSUB(ac,c,a);
60 +  VCROSS(r,ab,ac);
61  
62 < int
63 < vec3_equal(v1,v2)
64 < FVECT v1,v2;
65 < {
66 <   return(EQUAL(v1[0],v2[0]) && EQUAL(v1[1],v2[1])&& EQUAL(v1[2],v2[2]));
62 >  pr = DOT(p,r);        
63 >  ar = DOT(a,r);
64 >  /* Need to check for equality for degeneracy of 4 points on circle */
65 >  if( pr > ar *( 1.0 + EQUALITY_EPS))
66 >    return(TRUE);
67 >  else
68 >    return(FALSE);
69   }
70  
71 <
72 < int
73 < convex_angle(v0,v1,v2)
74 < FVECT v0,v1,v2;
71 > /*
72 > * tri_centroid(v0,v1,v2,c)
73 > *             : Average triangle vertices to give centroid: return in c
74 > *FVECT v0,v1,v2,c; : triangle vertices(v0,v1,v2) and vector to hold result(c)
75 > */                  
76 > tri_centroid(v0,v1,v2,c)
77 > FVECT v0,v1,v2,c;
78   {
79 <    FVECT cp01,cp12,cp;
80 <    
81 <    /* test sign of (v0Xv1)X(v1Xv2). v1 */
39 <    VCROSS(cp01,v0,v1);
40 <    VCROSS(cp12,v1,v2);
41 <    VCROSS(cp,cp01,cp12);
42 <    if(DOT(cp,v1) < 0)
43 <       return(FALSE);
44 <    return(TRUE);
79 >  c[0] = (v0[0] + v1[0] + v2[0])/3.0;
80 >  c[1] = (v0[1] + v1[1] + v2[1])/3.0;
81 >  c[2] = (v0[2] + v1[2] + v2[2])/3.0;
82   }
83  
47 – /* calculates the normal of a face contour using Newell's formula. e
84  
85 <               a =  SUMi (yi - yi+1)(zi + zi+1)
86 <               b =  SUMi (zi - zi+1)(xi + xi+1)
87 <               c =  SUMi (xi - xi+1)(yi + yi+1)
85 > /*
86 > * double  
87 > * tri_normal(v0,v1,v2,n,norm) : Calculates the normal of a face contour using
88 > *                            Newell's formula.
89 > * FVECT v0,v1,v2,n;  : Triangle vertices(v0,v1,v2) and vector for result(n)
90 > * int norm;          : If true result is normalized
91 > *
92 > * Triangle normal is calculated using the following:
93 > *    A =  SUMi (yi - yi+1)(zi + zi+1);
94 > *    B =  SUMi (zi - zi+1)(xi + xi+1)
95 > *    C =  SUMi (xi - xi+1)(yi + yi+1)
96   */
97   double
98   tri_normal(v0,v1,v2,n,norm)
99   FVECT v0,v1,v2,n;
100 < char norm;
100 > int norm;
101   {
102    double mag;
103  
104    n[0] = (v0[2] + v1[2]) * (v0[1] - v1[1]) +
105           (v1[2] + v2[2]) * (v1[1] - v2[1])  +
106           (v2[2] + v0[2]) * (v2[1] - v0[1]);
63 –  
107    n[1] = (v0[2] - v1[2]) * (v0[0] + v1[0]) +
108             (v1[2] - v2[2]) * (v1[0] + v2[0]) +
109 <           (v2[2] - v0[2]) * (v2[0] + v0[0]);
67 <
68 <  
109 >            (v2[2] - v0[2]) * (v2[0] + v0[0]);
110    n[2] = (v0[1] + v1[1]) * (v0[0] - v1[0]) +
111           (v1[1] + v2[1]) * (v1[0] - v2[0]) +
112           (v2[1] + v0[1]) * (v2[0] - v0[0]);
72 –
113    if(!norm)
114       return(0);
75 –
76 –  
115    mag = normalize(n);
78 –
116    return(mag);
117   }
118  
119 <
120 < tri_plane_equation(v0,v1,v2,n,nd,norm)
121 <   FVECT v0,v1,v2,n;
122 <   double *nd;
123 <   char norm;
119 > /*
120 > * tri_plane_equation(v0,v1,v2,peqptr,norm)
121 > *             : Calculates the plane equation (A,B,C,D) for triangle
122 > *               v0,v1,v2 ( Ax + By + Cz = D)
123 > * FVECT v0,v1,v2;   : Triangle vertices
124 > * FPEQ *peqptr;     : ptr to structure to hold result
125 > *  int norm;        : if TRUE, return unit normal
126 > */
127 > tri_plane_equation(v0,v1,v2,peqptr,norm)
128 >   FVECT v0,v1,v2;
129 >   FPEQ *peqptr;
130 >   int norm;
131   {
132 <    tri_normal(v0,v1,v2,n,norm);
133 <
90 <    *nd = -(DOT(n,v0));
132 >    tri_normal(v0,v1,v2,FP_N(*peqptr),norm);
133 >    FP_D(*peqptr) = -(DOT(FP_N(*peqptr),v0));
134   }
135  
136 + /*
137 + * int
138 + * intersect_ray_plane(orig,dir,peq,pd,r)
139 + *             : Intersects ray (orig,dir) with plane (peq). Returns TRUE
140 + *             if intersection occurs. If r!=NULL, sets with resulting i
141 + *             intersection point, and pd is set with parametric value of the
142 + *             intersection.
143 + * FVECT orig,dir;      : vectors defining the ray
144 + * FPEQ peq;            : plane equation
145 + * double *pd;          : holds resulting parametric intersection point
146 + * FVECT r;             : holds resulting intersection point
147 + *
148 + *    Plane is Ax + By + Cz +D = 0:
149 + *    A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
150 + *     t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
151 + *      line is  l = p1 + (p2-p1)t
152 + *    Solve for t
153 + */
154   int
155 < point_relative_to_plane(p,n,nd)
156 <   FVECT p,n;
157 <   double nd;
158 < {
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 <
109 < /* From quad_edge-code */
110 < int
111 < point_in_circle_thru_origin(p,p0,p1)
112 < FVECT p;
113 < FVECT p0,p1;
114 < {
115 <
116 <    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);    
122 <
123 <    det = -dp0*CROSS_VEC2(p1,p) + dp1*CROSS_VEC2(p0,p) - dp*CROSS_VEC2(p0,p1);
124 <    
125 <    return (det > 0);
126 < }
127 <
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 <
138 < int
139 < intersect_vector_plane(v,plane_n,plane_d,tptr,r)
140 <   FVECT v,plane_n;
141 <   double plane_d;
142 <   double *tptr;
155 > intersect_ray_plane(orig,dir,peq,pd,r)
156 >   FVECT orig,dir;
157 >   FPEQ peq;
158 >   double *pd;
159     FVECT r;
160   {
161 <  double t;
161 >  double t,d;
162    int hit;
147 –    /*
148 –      Plane is Ax + By + Cz +D = 0:
149 –      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
150 –    */
163  
164 <    /* line is  l = p1 + (p2-p1)t, p1=origin */
164 >  d = DOT(FP_N(peq),dir);
165 >  if(ZERO(d))
166 >     return(0);
167 >  t =  -(DOT(FP_N(peq),orig) + FP_D(peq))/d;
168  
169 <    /* Solve for t: */
155 <    t =  plane_d/-(DOT(plane_n,v));
156 <    if(t >0 || ZERO(t))
157 <       hit = 1;
158 <    else
169 >  if(t < 0)
170         hit = 0;
160 –    r[0] = v[0]*t;
161 –    r[1] = v[1]*t;
162 –    r[2] = v[2]*t;
163 –    if(tptr)
164 –       *tptr = t;
165 –  return(hit);
166 – }
167 –
168 – int
169 – intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
170 –   FVECT orig,dir;
171 –   FVECT plane_n;
172 –   double plane_d;
173 –   double *pd;
174 –   FVECT r;
175 – {
176 –  double t;
177 –  int hit;
178 –    /*
179 –      Plane is Ax + By + Cz +D = 0:
180 –      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
181 –    */
182 –     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
183 –         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
184 –       line is  l = p1 + (p2-p1)t
185 –     */
186 –    /* Solve for t: */
187 –    t =  -(DOT(plane_n,orig) + plane_d)/(DOT(plane_n,dir));
188 –    if(ZERO(t) || t >0)
189 –       hit = 1;
171      else
172 <       hit = 0;
172 >       hit = 1;
173 >  if(r)
174 >     VSUM(r,orig,dir,t);
175  
193 –  VSUM(r,orig,dir,t);
194 –
176      if(pd)
177         *pd = t;
178    return(hit);
179   }
180  
181 <
182 < int
183 < point_in_cone(p,p0,p1,p2)
184 < FVECT p;
185 < FVECT p0,p1,p2;
181 > /*
182 > * double
183 > * point_on_sphere(ps,p,c) : normalize p relative to sphere with center c
184 > * FVECT ps,p,c;       : ps Holds result vector,p is the original point,
185 > *                       and c is the sphere center
186 > */
187 > double
188 > point_on_sphere(ps,p,c)
189 > FVECT ps,p,c;
190   {
191 <    FVECT n;
192 <    FVECT np,x_axis,y_axis;
193 <    double d1,d2,d;
194 <    
195 <    /* Find the equation of the circle defined by the intersection
211 <       of the cone with the plane defined by p1,p2,p3- project p into
212 <       that plane and do an in-circle test in the plane
213 <     */
214 <    
215 <    /* find the equation of the plane defined by p1-p3 */
216 <    tri_plane_equation(p0,p1,p2,n,&d,FALSE);
217 <
218 <    /* define a coordinate system on the plane: the x axis is in
219 <       the direction of np2-np1, and the y axis is calculated from
220 <       n cross x-axis
221 <     */
222 <    /* Project p onto the plane */
223 <    if(!intersect_vector_plane(p,n,d,NULL,np))
224 <        return(FALSE);
225 <
226 <    /* create coordinate system on  plane: p2-p1 defines the x_axis*/
227 <    VSUB(x_axis,p1,p0);
228 <    normalize(x_axis);
229 <    /* The y axis is  */
230 <    VCROSS(y_axis,n,x_axis);
231 <    normalize(y_axis);
232 <
233 <    VSUB(p1,p1,p0);
234 <    VSUB(p2,p2,p0);
235 <    VSUB(np,np,p0);
236 <    
237 <    p1[0] = VLEN(p1);
238 <    p1[1] = 0;
239 <
240 <    d1 = DOT(p2,x_axis);
241 <    d2 = DOT(p2,y_axis);
242 <    p2[0] = d1;
243 <    p2[1] = d2;
244 <    
245 <    d1 = DOT(np,x_axis);
246 <    d2 = DOT(np,y_axis);
247 <    np[0] = d1;
248 <    np[1] = d2;
249 <
250 <    /* perform the in-circle test in the new coordinate system */
251 <    return(point_in_circle_thru_origin(np,p1,p2));
191 >  double d;
192 >  
193 >  VSUB(ps,p,c);    
194 >  d = normalize(ps);
195 >  return(d);
196   }
197  
198 + /*
199 + * int
200 + * point_in_stri(v0,v1,v2,p) : Return TRUE if p is in pyramid defined by
201 + *                            tri v0,v1,v2 and origin
202 + * FVECT v0,v1,v2,p;     :Triangle vertices(v0,v1,v2) and point in question(p)
203 + *
204 + *   Tests orientation of p relative to each edge (v0v1,v1v2,v2v0), if it is
205 + *   inside of all 3 edges, returns TRUE, else FALSE.
206 + */
207   int
208 < test_point_against_spherical_tri(v0,v1,v2,p,n,nset,which,sides)
208 > point_in_stri(v0,v1,v2,p)
209   FVECT v0,v1,v2,p;
257 – FVECT n[3];
258 – char *nset;
259 – char *which;
260 – char sides[3];
261 –
210   {
211 <    float d;
264 <
265 <    /* Find the normal to the triangle ORIGIN:v0:v1 */
266 <    if(!NTH_BIT(*nset,0))
267 <    {
268 <        VCROSS(n[0],v1,v0);
269 <        SET_NTH_BIT(*nset,0);
270 <    }
271 <    /* Test the point for sidedness */
272 <    d  = DOT(n[0],p);
273 <
274 <    if(ZERO(d))
275 <       sides[0] = GT_EDGE;
276 <    else
277 <       if(d > 0)
278 <      {
279 <          sides[0] =  GT_OUT;
280 <          sides[1] = sides[2] = GT_INVALID;
281 <          return(FALSE);
282 <      }
283 <    else
284 <       sides[0] = GT_INTERIOR;
285 <      
286 <    /* Test next edge */
287 <    if(!NTH_BIT(*nset,1))
288 <    {
289 <        VCROSS(n[1],v2,v1);
290 <        SET_NTH_BIT(*nset,1);
291 <    }
292 <    /* Test the point for sidedness */
293 <    d  = DOT(n[1],p);
294 <    if(ZERO(d))
295 <    {
296 <        sides[1] = GT_EDGE;
297 <        /* If on plane 0-and on plane 1: lies on edge */
298 <        if(sides[0] == GT_EDGE)
299 <        {
300 <            *which = 1;
301 <            sides[2] = GT_INVALID;
302 <            return(GT_EDGE);
303 <        }
304 <    }
305 <    else if(d > 0)
306 <    {
307 <        sides[1] = GT_OUT;
308 <        sides[2] = GT_INVALID;
309 <        return(FALSE);
310 <    }
311 <    else
312 <       sides[1] = GT_INTERIOR;
313 <    /* Test next edge */
314 <    if(!NTH_BIT(*nset,2))
315 <    {
316 <
317 <        VCROSS(n[2],v0,v2);
318 <        SET_NTH_BIT(*nset,2);
319 <    }
320 <    /* Test the point for sidedness */
321 <    d  = DOT(n[2],p);
322 <    if(ZERO(d))
323 <    {
324 <        sides[2] = GT_EDGE;
325 <
326 <        /* If on plane 0 and 2: lies on edge 0*/
327 <        if(sides[0] == GT_EDGE)
328 <           {
329 <               *which = 0;
330 <               return(GT_EDGE);
331 <           }
332 <        /* If on plane 1 and 2: lies on edge  2*/
333 <        if(sides[1] == GT_EDGE)
334 <           {
335 <               *which = 2;
336 <               return(GT_EDGE);
337 <           }
338 <        /* otherwise: on face 2 */
339 <        else
340 <           {
341 <               *which = 2;
342 <               return(GT_FACE);
343 <           }
344 <    }
345 <    else if(d > 0)
346 <      {
347 <        sides[2] = GT_OUT;
348 <        return(FALSE);
349 <      }
350 <    /* If on edge */
351 <    else
352 <       sides[2] = GT_INTERIOR;
353 <    
354 <    /* If on plane 0 only: on face 0 */
355 <    if(sides[0] == GT_EDGE)
356 <    {
357 <        *which = 0;
358 <        return(GT_FACE);
359 <    }
360 <    /* If on plane 1 only: on face 1 */
361 <    if(sides[1] == GT_EDGE)
362 <    {
363 <        *which = 1;
364 <        return(GT_FACE);
365 <    }
366 <    /* Must be interior to the pyramid */
367 <    return(GT_INTERIOR);
368 < }
369 <
370 <
371 <
372 <
373 < int
374 < test_single_point_against_spherical_tri(v0,v1,v2,p,which)
375 < FVECT v0,v1,v2,p;
376 < char *which;
377 < {
378 <    float d;
211 >    double d;
212      FVECT n;  
380 –    char sides[3];
213  
214 <    /* First test if point coincides with any of the vertices */
383 <    if(EQUAL_VEC3(p,v0))
384 <    {
385 <        *which = 0;
386 <        return(GT_VERTEX);
387 <    }
388 <    if(EQUAL_VEC3(p,v1))
389 <    {
390 <        *which = 1;
391 <        return(GT_VERTEX);
392 <    }
393 <    if(EQUAL_VEC3(p,v2))
394 <    {
395 <        *which = 2;
396 <        return(GT_VERTEX);
397 <    }
398 <    VCROSS(n,v1,v0);
214 >    VCROSS(n,v0,v1);
215      /* Test the point for sidedness */
216      d  = DOT(n,p);
217 <    if(ZERO(d))
218 <       sides[0] = GT_EDGE;
403 <    else
404 <       if(d > 0)
405 <          return(FALSE);
406 <       else
407 <          sides[0] = GT_INTERIOR;
217 >    if(d > 0.0)
218 >      return(FALSE);
219      /* Test next edge */
220 <    VCROSS(n,v2,v1);
220 >    VCROSS(n,v1,v2);
221      /* Test the point for sidedness */
222      d  = DOT(n,p);
223 <    if(ZERO(d))
413 <    {
414 <        sides[1] = GT_EDGE;
415 <        /* If on plane 0-and on plane 1: lies on edge */
416 <        if(sides[0] == GT_EDGE)
417 <        {
418 <            *which = 1;
419 <            return(GT_VERTEX);
420 <        }
421 <    }
422 <    else if(d > 0)
223 >    if(d > 0.0)
224         return(FALSE);
424 –    else
425 –       sides[1] = GT_INTERIOR;
426 –
225      /* Test next edge */
226 <    VCROSS(n,v0,v2);
226 >    VCROSS(n,v2,v0);
227      /* Test the point for sidedness */
228      d  = DOT(n,p);
229 <    if(ZERO(d))
432 <    {
433 <        sides[2] = GT_EDGE;
434 <        
435 <        /* If on plane 0 and 2: lies on edge 0*/
436 <        if(sides[0] == GT_EDGE)
437 <        {
438 <            *which = 0;
439 <            return(GT_VERTEX);
440 <        }
441 <        /* If on plane 1 and 2: lies on edge  2*/
442 <        if(sides[1] == GT_EDGE)
443 <        {
444 <            *which = 2;
445 <            return(GT_VERTEX);
446 <        }
447 <        /* otherwise: on face 2 */
448 <        else
449 <       {
450 <           return(GT_FACE);
451 <       }
452 <    }
453 <    else if(d > 0)
229 >    if(d > 0.0)
230         return(FALSE);
231      /* Must be interior to the pyramid */
232 <    return(GT_FACE);
232 >    return(TRUE);
233   }
234  
235 + /*
236 + * int
237 + * ray_intersect_tri(orig,dir,v0,v1,v2,pt)
238 + *    : test if ray orig-dir intersects triangle v0v1v2, result in pt
239 + * FVECT orig,dir;   : Vectors defining ray origin and direction
240 + * FVECT v0,v1,v2;   : Triangle vertices
241 + * FVECT pt;         : Intersection point (if any)
242 + */
243   int
244 < test_vertices_for_tri_inclusion(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides)
461 < FVECT t0,t1,t2,p0,p1,p2;
462 < char *nset;
463 < FVECT n[3];
464 < FVECT avg;
465 < char pt_sides[3][3];
466 <
467 < {
468 <    char below_plane[3],on_edge,test;
469 <    char which;
470 <
471 <    SUM_3VEC3(avg,t0,t1,t2);
472 <    on_edge = 0;
473 <    *nset = 0;
474 <    /* Test vertex v[i] against triangle j*/
475 <    /* Check if v[i] lies below plane defined by avg of 3 vectors
476 <       defining triangle
477 <       */
478 <
479 <    /* test point 0 */
480 <    if(DOT(avg,p0) < 0)
481 <      below_plane[0] = 1;
482 <    else
483 <      below_plane[0]=0;
484 <    /* Test if b[i] lies in or on triangle a */
485 <    test = test_point_against_spherical_tri(t0,t1,t2,p0,
486 <                                                 n,nset,&which,pt_sides[0]);
487 <    /* If pts[i] is interior: done */
488 <    if(!below_plane[0])
489 <      {
490 <        if(test == GT_INTERIOR)
491 <          return(TRUE);
492 <        /* Remember if b[i] fell on one of the 3 defining planes */
493 <        if(test)
494 <          on_edge++;
495 <      }
496 <    /* Now test point 1*/
497 <
498 <    if(DOT(avg,p1) < 0)
499 <      below_plane[1] = 1;
500 <    else
501 <      below_plane[1]=0;
502 <    /* Test if b[i] lies in or on triangle a */
503 <    test = test_point_against_spherical_tri(t0,t1,t2,p1,
504 <                                                 n,nset,&which,pt_sides[1]);
505 <    /* If pts[i] is interior: done */
506 <    if(!below_plane[1])
507 <    {
508 <      if(test == GT_INTERIOR)
509 <        return(TRUE);
510 <      /* Remember if b[i] fell on one of the 3 defining planes */
511 <      if(test)
512 <        on_edge++;
513 <    }
514 <    
515 <    /* Now test point 2 */
516 <    if(DOT(avg,p2) < 0)
517 <      below_plane[2] = 1;
518 <    else
519 <      below_plane[2]=0;
520 <        /* Test if b[i] lies in or on triangle a */
521 <    test = test_point_against_spherical_tri(t0,t1,t2,p2,
522 <                                                 n,nset,&which,pt_sides[2]);
523 <
524 <    /* If pts[i] is interior: done */
525 <    if(!below_plane[2])
526 <      {
527 <        if(test == GT_INTERIOR)
528 <          return(TRUE);
529 <        /* Remember if b[i] fell on one of the 3 defining planes */
530 <        if(test)
531 <          on_edge++;
532 <      }
533 <
534 <    /* If all three points below separating plane: trivial reject */
535 <    if(below_plane[0] && below_plane[1] && below_plane[2])
536 <       return(FALSE);
537 <    /* Accept if all points lie on a triangle vertex/edge edge- accept*/
538 <    if(on_edge == 3)
539 <       return(TRUE);
540 <    /* Now check vertices in a against triangle b */
541 <    return(FALSE);
542 < }
543 <
544 <
545 < set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,n)
546 <   FVECT t0,t1,t2,p0,p1,p2;
547 <   char test[3];
548 <   char sides[3][3];
549 <   char nset;
550 <   FVECT n[3];
551 < {
552 <    char t;
553 <    double d;
554 <
555 <    
556 <    /* p=0 */
557 <    test[0] = 0;
558 <    if(sides[0][0] == GT_INVALID)
559 <    {
560 <      if(!NTH_BIT(nset,0))
561 <        VCROSS(n[0],t1,t0);
562 <      /* Test the point for sidedness */
563 <      d  = DOT(n[0],p0);
564 <      if(d >= 0)
565 <        SET_NTH_BIT(test[0],0);
566 <    }
567 <    else
568 <      if(sides[0][0] != GT_INTERIOR)
569 <        SET_NTH_BIT(test[0],0);
570 <
571 <    if(sides[0][1] == GT_INVALID)
572 <    {
573 <      if(!NTH_BIT(nset,1))
574 <        VCROSS(n[1],t2,t1);
575 <        /* Test the point for sidedness */
576 <        d  = DOT(n[1],p0);
577 <        if(d >= 0)
578 <          SET_NTH_BIT(test[0],1);
579 <    }
580 <    else
581 <      if(sides[0][1] != GT_INTERIOR)
582 <        SET_NTH_BIT(test[0],1);
583 <
584 <    if(sides[0][2] == GT_INVALID)
585 <    {
586 <      if(!NTH_BIT(nset,2))
587 <        VCROSS(n[2],t0,t2);
588 <      /* Test the point for sidedness */
589 <      d  = DOT(n[2],p0);
590 <      if(d >= 0)
591 <        SET_NTH_BIT(test[0],2);
592 <    }
593 <    else
594 <      if(sides[0][2] != GT_INTERIOR)
595 <        SET_NTH_BIT(test[0],2);
596 <
597 <    /* p=1 */
598 <    test[1] = 0;
599 <    /* t=0*/
600 <    if(sides[1][0] == GT_INVALID)
601 <    {
602 <      if(!NTH_BIT(nset,0))
603 <        VCROSS(n[0],t1,t0);
604 <      /* Test the point for sidedness */
605 <      d  = DOT(n[0],p1);
606 <      if(d >= 0)
607 <        SET_NTH_BIT(test[1],0);
608 <    }
609 <    else
610 <      if(sides[1][0] != GT_INTERIOR)
611 <        SET_NTH_BIT(test[1],0);
612 <    
613 <    /* t=1 */
614 <    if(sides[1][1] == GT_INVALID)
615 <    {
616 <      if(!NTH_BIT(nset,1))
617 <        VCROSS(n[1],t2,t1);
618 <      /* Test the point for sidedness */
619 <      d  = DOT(n[1],p1);
620 <      if(d >= 0)
621 <        SET_NTH_BIT(test[1],1);
622 <    }
623 <    else
624 <      if(sides[1][1] != GT_INTERIOR)
625 <        SET_NTH_BIT(test[1],1);
626 <      
627 <    /* t=2 */
628 <    if(sides[1][2] == GT_INVALID)
629 <    {
630 <      if(!NTH_BIT(nset,2))
631 <        VCROSS(n[2],t0,t2);
632 <      /* Test the point for sidedness */
633 <      d  = DOT(n[2],p1);
634 <      if(d >= 0)
635 <        SET_NTH_BIT(test[1],2);
636 <    }
637 <    else
638 <      if(sides[1][2] != GT_INTERIOR)
639 <        SET_NTH_BIT(test[1],2);
640 <
641 <    /* p=2 */
642 <    test[2] = 0;
643 <    /* t = 0 */
644 <    if(sides[2][0] == GT_INVALID)
645 <    {
646 <      if(!NTH_BIT(nset,0))
647 <        VCROSS(n[0],t1,t0);
648 <      /* Test the point for sidedness */
649 <      d  = DOT(n[0],p2);
650 <      if(d >= 0)
651 <        SET_NTH_BIT(test[2],0);
652 <    }
653 <    else
654 <      if(sides[2][0] != GT_INTERIOR)
655 <        SET_NTH_BIT(test[2],0);
656 <    /* t=1 */
657 <    if(sides[2][1] == GT_INVALID)
658 <    {
659 <      if(!NTH_BIT(nset,1))
660 <        VCROSS(n[1],t2,t1);
661 <      /* Test the point for sidedness */
662 <      d  = DOT(n[1],p2);
663 <      if(d >= 0)
664 <        SET_NTH_BIT(test[2],1);
665 <    }
666 <    else
667 <      if(sides[2][1] != GT_INTERIOR)
668 <        SET_NTH_BIT(test[2],1);
669 <    /* t=2 */
670 <    if(sides[2][2] == GT_INVALID)
671 <    {
672 <      if(!NTH_BIT(nset,2))
673 <        VCROSS(n[2],t0,t2);
674 <      /* Test the point for sidedness */
675 <      d  = DOT(n[2],p2);
676 <      if(d >= 0)
677 <        SET_NTH_BIT(test[2],2);
678 <    }
679 <    else
680 <      if(sides[2][2] != GT_INTERIOR)
681 <        SET_NTH_BIT(test[2],2);
682 < }
683 <
684 <
685 < int
686 < spherical_tri_intersect(a1,a2,a3,b1,b2,b3)
687 < FVECT a1,a2,a3,b1,b2,b3;
688 < {
689 <  char which,test,n_set[2];
690 <  char sides[2][3][3],i,j,inext,jnext;
691 <  char tests[2][3];
692 <  FVECT n[2][3],p,avg[2];
693 <
694 <  /* Test the vertices of triangle a against the pyramid formed by triangle
695 <     b and the origin. If any vertex of a is interior to triangle b, or
696 <     if all 3 vertices of a are ON the edges of b,return TRUE. Remember
697 <     the results of the edge normal and sidedness tests for later.
698 <   */
699 < if(test_vertices_for_tri_inclusion(a1,a2,a3,b1,b2,b3,
700 <                                    &(n_set[0]),n[0],avg[0],sides[1]))
701 <     return(TRUE);
702 <  
703 < if(test_vertices_for_tri_inclusion(b1,b2,b3,a1,a2,a3,
704 <                                    &(n_set[1]),n[1],avg[1],sides[0]))
705 <     return(TRUE);
706 <
707 <
708 <  set_sidedness_tests(b1,b2,b3,a1,a2,a3,tests[0],sides[0],n_set[1],n[1]);
709 <  if(tests[0][0]&tests[0][1]&tests[0][2])
710 <    return(FALSE);
711 <
712 <  set_sidedness_tests(a1,a2,a3,b1,b2,b3,tests[1],sides[1],n_set[0],n[0]);
713 <  if(tests[1][0]&tests[1][1]&tests[1][2])
714 <    return(FALSE);
715 <
716 <  for(j=0; j < 3;j++)
717 <  {
718 <    jnext = (j+1)%3;
719 <    /* IF edge b doesnt cross any great circles of a, punt */
720 <    if(tests[1][j] & tests[1][jnext])
721 <      continue;
722 <    for(i=0;i<3;i++)
723 <    {
724 <      inext = (i+1)%3;
725 <      /* IF edge a doesnt cross any great circles of b, punt */
726 <      if(tests[0][i] & tests[0][inext])
727 <        continue;
728 <      /* Now find the great circles that cross and test */
729 <      if((NTH_BIT(tests[0][i],j)^(NTH_BIT(tests[0][inext],j)))
730 <          && (NTH_BIT(tests[1][j],i)^NTH_BIT(tests[1][jnext],i)))
731 <      {
732 <        VCROSS(p,n[0][i],n[1][j]);
733 <                    
734 <        /* If zero cp= done */
735 <        if(ZERO_VEC3(p))
736 <          continue;
737 <        /* check above both planes */
738 <        if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0)
739 <          {
740 <            NEGATE_VEC3(p);
741 <            if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0)
742 <              continue;
743 <          }
744 <        return(TRUE);
745 <      }
746 <    }
747 <  }
748 <  return(FALSE);
749 < }
750 <
751 < int
752 < ray_intersect_tri(orig,dir,v0,v1,v2,pt,wptr)
244 > ray_intersect_tri(orig,dir,v0,v1,v2,pt)
245   FVECT orig,dir;
246   FVECT v0,v1,v2;
247   FVECT pt;
756 – char *wptr;
248   {
249 <  FVECT p0,p1,p2,p,n;
250 <  char type,which;
251 <  double pd;
761 <  
762 <  point_on_sphere(p0,v0,orig);
763 <  point_on_sphere(p1,v1,orig);
764 <  point_on_sphere(p2,v2,orig);
765 <  type = test_single_point_against_spherical_tri(p0,p1,p2,dir,&which);
249 >  FVECT p0,p1,p2,p;
250 >  FPEQ peq;
251 >  int type;
252  
253 <  if(type)
253 >  VSUB(p0,v0,orig);
254 >  VSUB(p1,v1,orig);
255 >  VSUB(p2,v2,orig);
256 >
257 >  if(point_in_stri(p0,p1,p2,dir))
258    {
259        /* Intersect the ray with the triangle plane */
260 <      tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
261 <      intersect_ray_plane(orig,dir,n,pd,NULL,pt);        
260 >      tri_plane_equation(v0,v1,v2,&peq,FALSE);
261 >      return(intersect_ray_plane(orig,dir,peq,NULL,pt));
262    }
263 <  if(wptr)
774 <    *wptr = which;
775 <
776 <  return(type);
263 >  return(FALSE);
264   }
265  
266 <
266 > /*
267 > * calculate_view_frustum(vp,hv,vv,horiz,vert,near,far,fnear,ffar)
268 > *    : Calculate vertices defining front and rear clip rectangles of
269 > *      view frustum defined by vp,hv,vv,horiz,vert,near, and far and
270 > *      return in fnear and ffar.
271 > * FVECT vp,hv,vv;        : Viewpoint(vp),hv and vv are the horizontal and
272 > *                          vertical vectors in the view frame-magnitude is
273 > *                          the dimension of the front frustum face at z =1
274 > * double horiz,vert,near,far; : View angle horizontal and vertical(horiz,vert)
275 > *                              and distance to the near,far clipping planes
276 > * FVECT fnear[4],ffar[4];     : holds results
277 > *
278 > */
279   calculate_view_frustum(vp,hv,vv,horiz,vert,near,far,fnear,ffar)
280   FVECT vp,hv,vv;
281   double horiz,vert,near,far;
# Line 786 | Line 285 | FVECT fnear[4],ffar[4];
285      FVECT t,nhv,nvv,ndv;
286      double w2,h2;
287      /* Calculate the x and y dimensions of the near face */
789 –    /* hv and vv are the horizontal and vertical vectors in the
790 –       view frame-the magnitude is the dimension of the front frustum
791 –       face at z =1
792 –    */
288      VCOPY(nhv,hv);
289      VCOPY(nvv,vv);
290      w2 = normalize(nhv);
# Line 833 | Line 328 | FVECT fnear[4],ffar[4];
328   }
329  
330  
331 <
332 <
333 < int
334 < spherical_polygon_edge_intersect(a0,a1,b0,b1)
335 < FVECT a0,a1,b0,b1;
336 < {
337 <    FVECT na,nb,avga,avgb,p;
843 <    double d;
844 <    int sb0,sb1,sa0,sa1;
845 <
846 <    /* First test if edge b straddles great circle of a */
847 <    VCROSS(na,a0,a1);
848 <    d = DOT(na,b0);
849 <    sb0 = ZERO(d)?0:(d<0)? -1: 1;
850 <    d = DOT(na,b1);
851 <    sb1 = ZERO(d)?0:(d<0)? -1: 1;
852 <    /* edge b entirely on one side of great circle a: edges cannot intersect*/
853 <    if(sb0*sb1 > 0)
854 <       return(FALSE);
855 <    /* test if edge a straddles great circle of b */
856 <    VCROSS(nb,b0,b1);
857 <    d = DOT(nb,a0);
858 <    sa0 = ZERO(d)?0:(d<0)? -1: 1;
859 <    d = DOT(nb,a1);
860 <    sa1 = ZERO(d)?0:(d<0)? -1: 1;
861 <    /* edge a entirely on one side of great circle b: edges cannot intersect*/
862 <    if(sa0*sa1 > 0)
863 <       return(FALSE);
864 <
865 <    /* Find one of intersection points of the great circles */
866 <    VCROSS(p,na,nb);
867 <    /* If they lie on same great circle: call an intersection */
868 <    if(ZERO_VEC3(p))
869 <       return(TRUE);
870 <
871 <    VADD(avga,a0,a1);
872 <    VADD(avgb,b0,b1);
873 <    if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
874 <    {
875 <      NEGATE_VEC3(p);
876 <      if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
877 <        return(FALSE);
878 <    }
879 <    if((!sb0 || !sb1) && (!sa0 || !sa1))
880 <      return(FALSE);
881 <    return(TRUE);
882 < }
883 <
884 <
885 <
886 < /* Find the normalized barycentric coordinates of p relative to
887 < * triangle v0,v1,v2. Return result in coord
331 > /*
332 > * bary2d(x1,y1,x2,y2,x3,y3,px,py,coord)
333 > *  : Find the normalized barycentric coordinates of p relative to
334 > *    triangle v0,v1,v2. Return result in coord
335 > * double x1,y1,x2,y2,x3,y3; : defines triangle vertices 1,2,3
336 > * double px,py;             : coordinates of pt
337 > * double coord[3];          : result
338   */
339   bary2d(x1,y1,x2,y2,x3,y3,px,py,coord)
340   double x1,y1,x2,y2,x3,y3;
# Line 896 | Line 346 | double coord[3];
346    a =  (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1);
347    coord[0] = ((x2 - px) * (y3 - py) - (x3 - px) * (y2 - py)) / a;
348    coord[1] = ((x3 - px) * (y1 - py) - (x1 - px) * (y3 - py)) / a;
349 <  coord[2]  = 1.0 - coord[0] - coord[1];
349 >  coord[2] = ((x1 - px) * (y2 - py) - (x2 - px) * (y1 - py)) / a;
350  
351   }
352  
903 – int
904 – bary2d_child(coord)
905 – double coord[3];
906 – {
907 –  int i;
353  
909 –  /* First check if one of the original vertices */
910 –  for(i=0;i<3;i++)
911 –    if(EQUAL(coord[i],1.0))
912 –      return(i);
354  
914 –  /* Check if one of the new vertices: for all return center child */
915 –  if(ZERO(coord[0]) && EQUAL(coord[1],0.5))
916 –  {
917 –    coord[0] = 1.0f;
918 –    coord[1] = 0.0f;
919 –    coord[2] = 0.0f;
920 –    return(3);
921 –  }
922 –  if(ZERO(coord[1]) && EQUAL(coord[0],0.5))
923 –  {
924 –    coord[0] = 0.0f;
925 –    coord[1] = 1.0f;
926 –    coord[2] = 0.0f;
927 –    return(3);
928 –  }
929 –  if(ZERO(coord[2]) && EQUAL(coord[0],0.5))
930 –    {
931 –      coord[0] = 0.0f;
932 –      coord[1] = 0.0f;
933 –      coord[2] = 1.0f;
934 –      return(3);
935 –    }
355  
937 –  /* Otherwise return child */
938 –  if(coord[0] > 0.5)
939 –  {
940 –      /* update bary for child */
941 –      coord[0] = 2.0*coord[0]- 1.0;
942 –      coord[1] *= 2.0;
943 –      coord[2] *= 2.0;
944 –      return(0);
945 –  }
946 –  else
947 –    if(coord[1] > 0.5)
948 –    {
949 –      coord[0] *= 2.0;
950 –      coord[1] = 2.0*coord[1]- 1.0;
951 –      coord[2] *= 2.0;
952 –      return(1);
953 –    }
954 –    else
955 –      if(coord[2] > 0.5)
956 –      {
957 –        coord[0] *= 2.0;
958 –        coord[1] *= 2.0;
959 –        coord[2] = 2.0*coord[2]- 1.0;
960 –        return(2);
961 –      }
962 –      else
963 –         {
964 –           coord[0] = 1.0 - 2.0*coord[0];
965 –           coord[1] = 1.0 - 2.0*coord[1];
966 –           coord[2] = 1.0 - 2.0*coord[2];
967 –           return(3);
968 –         }
969 – }
970 –
971 – int
972 – max_index(v)
973 – FVECT v;
974 – {
975 –  double a,b,c;
976 –  int i;
977 –
978 –  a = fabs(v[0]);
979 –  b = fabs(v[1]);
980 –  c = fabs(v[2]);
981 –  i = (a>=b)?((a>=c)?0:2):((b>=c)?1:2);  
982 –  return(i);
983 – }
356  
357  
358  

Diff Legend

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