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.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;
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 <
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,pd,r)
140 <   FVECT v,plane_n;
141 <   double plane_d;
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(pd)
164 –       *pd = 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 –    /* Solve for t: */
186 –    t =  -(DOT(plane_n,orig) + plane_d)/(DOT(plane_n,dir));
187 –    if(ZERO(t) || t >0)
188 –       hit = 1;
171      else
172 <       hit = 0;
172 >       hit = 1;
173 >  if(r)
174 >     VSUM(r,orig,dir,t);
175  
192 –  VSUM(r,orig,dir,t);
193 –
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
210 <       of the cone with the plane defined by p1,p2,p3- project p into
211 <       that plane and do an in-circle test in the plane
212 <     */
213 <    
214 <    /* find the equation of the plane defined by p1-p3 */
215 <    tri_plane_equation(p0,p1,p2,n,&d,FALSE);
216 <
217 <    /* 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);
224 <
225 <    /* create coordinate system on  plane: p2-p1 defines the x_axis*/
226 <    VSUB(x_axis,p1,p0);
227 <    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));
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;
256 – FVECT n[3];
257 – char *nset;
258 – char *which;
259 – char sides[3];
260 –
210   {
211 <    float d;
263 <
264 <    /* Find the normal to the triangle ORIGIN:v0:v1 */
265 <    if(!NTH_BIT(*nset,0))
266 <    {
267 <        VCROSS(n[0],v1,v0);
268 <        SET_NTH_BIT(*nset,0);
269 <    }
270 <    /* Test the point for sidedness */
271 <    d  = DOT(n[0],p);
272 <
273 <    if(ZERO(d))
274 <       sides[0] = GT_EDGE;
275 <    else
276 <       if(d > 0)
277 <      {
278 <          sides[0] =  GT_OUT;
279 <          sides[1] = sides[2] = GT_INVALID;
280 <          return(FALSE);
281 <      }
282 <    else
283 <       sides[0] = GT_INTERIOR;
284 <      
285 <    /* Test next edge */
286 <    if(!NTH_BIT(*nset,1))
287 <    {
288 <        VCROSS(n[1],v2,v1);
289 <        SET_NTH_BIT(*nset,1);
290 <    }
291 <    /* Test the point for sidedness */
292 <    d  = DOT(n[1],p);
293 <    if(ZERO(d))
294 <    {
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 <    {
306 <        sides[1] = GT_OUT;
307 <        sides[2] = GT_INVALID;
308 <        return(FALSE);
309 <    }
310 <    else
311 <       sides[1] = GT_INTERIOR;
312 <    /* Test next edge */
313 <    if(!NTH_BIT(*nset,2))
314 <    {
315 <
316 <        VCROSS(n[2],v0,v2);
317 <        SET_NTH_BIT(*nset,2);
318 <    }
319 <    /* Test the point for sidedness */
320 <    d  = DOT(n[2],p);
321 <    if(ZERO(d))
322 <    {
323 <        sides[2] = GT_EDGE;
324 <
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 <           }
343 <    }
344 <    else if(d > 0)
345 <      {
346 <        sides[2] = GT_OUT;
347 <        return(FALSE);
348 <      }
349 <    /* If on edge */
350 <    else
351 <       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 <    }
365 <    /* Must be interior to the pyramid */
366 <    return(GT_INTERIOR);
367 < }
368 <
369 <
370 <
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;
211 >    double d;
212      FVECT n;  
379 –    char sides[3];
213  
214 <    /* 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);
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;
402 <    else
403 <       if(d > 0)
404 <          return(FALSE);
405 <       else
406 <          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))
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)
223 >    if(d > 0.0)
224         return(FALSE);
423 –    else
424 –       sides[1] = GT_INTERIOR;
425 –
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))
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)
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)
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 <
544 < set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,n)
545 <   FVECT t0,t1,t2,p0,p1,p2;
546 <   char test[3];
547 <   char sides[3][3];
548 <   char nset;
549 <   FVECT n[3];
550 < {
551 <    char t;
552 <    double d;
553 <
554 <    
555 <    /* p=0 */
556 <    test[0] = 0;
557 <    if(sides[0][0] == GT_INVALID)
558 <    {
559 <      if(!NTH_BIT(nset,0))
560 <        VCROSS(n[0],t1,t0);
561 <      /* Test the point for sidedness */
562 <      d  = DOT(n[0],p0);
563 <      if(d >= 0)
564 <        SET_NTH_BIT(test[0],0);
565 <    }
566 <    else
567 <      if(sides[0][0] != GT_INTERIOR)
568 <        SET_NTH_BIT(test[0],0);
569 <
570 <    if(sides[0][1] == GT_INVALID)
571 <    {
572 <      if(!NTH_BIT(nset,1))
573 <        VCROSS(n[1],t2,t1);
574 <        /* Test the point for sidedness */
575 <        d  = DOT(n[1],p0);
576 <        if(d >= 0)
577 <          SET_NTH_BIT(test[0],1);
578 <    }
579 <    else
580 <      if(sides[0][1] != GT_INTERIOR)
581 <        SET_NTH_BIT(test[0],1);
582 <
583 <    if(sides[0][2] == GT_INVALID)
584 <    {
585 <      if(!NTH_BIT(nset,2))
586 <        VCROSS(n[2],t0,t2);
587 <      /* Test the point for sidedness */
588 <      d  = DOT(n[2],p0);
589 <      if(d >= 0)
590 <        SET_NTH_BIT(test[0],2);
591 <    }
592 <    else
593 <      if(sides[0][2] != GT_INTERIOR)
594 <        SET_NTH_BIT(test[0],2);
595 <
596 <    /* p=1 */
597 <    test[1] = 0;
598 <    /* t=0*/
599 <    if(sides[1][0] == GT_INVALID)
600 <    {
601 <      if(!NTH_BIT(nset,0))
602 <        VCROSS(n[0],t1,t0);
603 <      /* Test the point for sidedness */
604 <      d  = DOT(n[0],p1);
605 <      if(d >= 0)
606 <        SET_NTH_BIT(test[1],0);
607 <    }
608 <    else
609 <      if(sides[1][0] != GT_INTERIOR)
610 <        SET_NTH_BIT(test[1],0);
611 <    
612 <    /* t=1 */
613 <    if(sides[1][1] == GT_INVALID)
614 <    {
615 <      if(!NTH_BIT(nset,1))
616 <        VCROSS(n[1],t2,t1);
617 <      /* Test the point for sidedness */
618 <      d  = DOT(n[1],p1);
619 <      if(d >= 0)
620 <        SET_NTH_BIT(test[1],1);
621 <    }
622 <    else
623 <      if(sides[1][1] != GT_INTERIOR)
624 <        SET_NTH_BIT(test[1],1);
625 <      
626 <    /* t=2 */
627 <    if(sides[1][2] == GT_INVALID)
628 <    {
629 <      if(!NTH_BIT(nset,2))
630 <        VCROSS(n[2],t0,t2);
631 <      /* Test the point for sidedness */
632 <      d  = DOT(n[2],p1);
633 <      if(d >= 0)
634 <        SET_NTH_BIT(test[1],2);
635 <    }
636 <    else
637 <      if(sides[1][2] != GT_INTERIOR)
638 <        SET_NTH_BIT(test[1],2);
639 <
640 <    /* p=2 */
641 <    test[2] = 0;
642 <    /* t = 0 */
643 <    if(sides[2][0] == GT_INVALID)
644 <    {
645 <      if(!NTH_BIT(nset,0))
646 <        VCROSS(n[0],t1,t0);
647 <      /* Test the point for sidedness */
648 <      d  = DOT(n[0],p2);
649 <      if(d >= 0)
650 <        SET_NTH_BIT(test[2],0);
651 <    }
652 <    else
653 <      if(sides[2][0] != GT_INTERIOR)
654 <        SET_NTH_BIT(test[2],0);
655 <    /* t=1 */
656 <    if(sides[2][1] == GT_INVALID)
657 <    {
658 <      if(!NTH_BIT(nset,1))
659 <        VCROSS(n[1],t2,t1);
660 <      /* Test the point for sidedness */
661 <      d  = DOT(n[1],p2);
662 <      if(d >= 0)
663 <        SET_NTH_BIT(test[2],1);
664 <    }
665 <    else
666 <      if(sides[2][1] != GT_INTERIOR)
667 <        SET_NTH_BIT(test[2],1);
668 <    /* t=2 */
669 <    if(sides[2][2] == GT_INVALID)
670 <    {
671 <      if(!NTH_BIT(nset,2))
672 <        VCROSS(n[2],t0,t2);
673 <      /* Test the point for sidedness */
674 <      d  = DOT(n[2],p2);
675 <      if(d >= 0)
676 <        SET_NTH_BIT(test[2],2);
677 <    }
678 <    else
679 <      if(sides[2][2] != GT_INTERIOR)
680 <        SET_NTH_BIT(test[2],2);
681 < }
682 <
683 <
684 < int
685 < spherical_tri_intersect(a1,a2,a3,b1,b2,b3)
686 < FVECT a1,a2,a3,b1,b2,b3;
687 < {
688 <  char which,test,n_set[2];
689 <  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);
705 <
706 <
707 <  set_sidedness_tests(b1,b2,b3,a1,a2,a3,tests[0],sides[0],n_set[1],n[1]);
708 <  if(tests[0][0]&tests[0][1]&tests[0][2])
709 <    return(FALSE);
710 <
711 <  set_sidedness_tests(a1,a2,a3,b1,b2,b3,tests[1],sides[1],n_set[0],n[0]);
712 <  if(tests[1][0]&tests[1][1]&tests[1][2])
713 <    return(FALSE);
714 <
715 <  for(j=0; j < 3;j++)
716 <  {
717 <    jnext = (j+1)%3;
718 <    /* 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);
748 < }
749 <
750 < int
751 < 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;
755 – char *wptr;
248   {
249 <  FVECT p0,p1,p2,p,n;
250 <  char type,which;
251 <  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);
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)
773 <    *wptr = which;
774 <
775 <  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 785 | 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 */
788 –    /* hv and vv are the horizontal and vertical vectors in the
789 –       view frame-the magnitude is the dimension of the front frustum
790 –       face at z =1
791 –    */
288      VCOPY(nhv,hv);
289      VCOPY(nvv,vv);
290      w2 = normalize(nhv);
# Line 832 | Line 328 | FVECT fnear[4],ffar[4];
328   }
329  
330  
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;
341 + double px,py;
342 + double coord[3];
343 + {
344 +  double a;
345  
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] = ((x1 - px) * (y2 - py) - (x2 - px) * (y1 - py)) / a;
350 +
351 + }
352  
837 – int
838 – spherical_polygon_edge_intersect(a0,a1,b0,b1)
839 – FVECT a0,a1,b0,b1;
840 – {
841 –    FVECT na,nb,avga,avgb,p;
842 –    double d;
843 –    int sb0,sb1,sa0,sa1;
353  
845 –    /* First test if edge b straddles great circle of a */
846 –    VCROSS(na,a0,a1);
847 –    d = DOT(na,b0);
848 –    sb0 = ZERO(d)?0:(d<0)? -1: 1;
849 –    d = DOT(na,b1);
850 –    sb1 = ZERO(d)?0:(d<0)? -1: 1;
851 –    /* edge b entirely on one side of great circle a: edges cannot intersect*/
852 –    if(sb0*sb1 > 0)
853 –       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);
354  
864 –    /* Find one of intersection points of the great circles */
865 –    VCROSS(p,na,nb);
866 –    /* If they lie on same great circle: call an intersection */
867 –    if(ZERO_VEC3(p))
868 –       return(TRUE);
355  
356 <    VADD(avga,a0,a1);
357 <    VADD(avgb,b0,b1);
358 <    if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
359 <    {
360 <      NEGATE_VEC3(p);
361 <      if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
362 <        return(FALSE);
363 <    }
364 <    if((!sb0 || !sb1) && (!sa0 || !sa1))
365 <      return(FALSE);
880 <    return(TRUE);
881 < }
356 >
357 >
358 >
359 >
360 >
361 >
362 >
363 >
364 >
365 >
366  

Diff Legend

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