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.5 by gwlarson, Mon Sep 14 10:33:46 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.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)
# Line 60 | Line 104 | int norm;
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;
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 < /* From quad_edge-code */
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_in_circle_thru_origin(p,p0,p1)
156 < FVECT p;
157 < FVECT p0,p1;
158 < {
99 <
100 <    double dp0,dp1;
101 <    double dp,det;
102 <    
103 <    dp0 = DOT_VEC2(p0,p0);
104 <    dp1 = DOT_VEC2(p1,p1);
105 <    dp  = DOT_VEC2(p,p);    
106 <
107 <    det = -dp0*CROSS_VEC2(p1,p) + dp1*CROSS_VEC2(p0,p) - dp*CROSS_VEC2(p0,p1);
108 <    
109 <    return (det > 0);
110 < }
111 <
112 <
113 <
114 < point_on_sphere(ps,p,c)
115 < FVECT ps,p,c;
116 < {
117 <    VSUB(ps,p,c);    
118 <    normalize(ps);
119 < }
120 <
121 <
122 < /* returns TRUE if ray from origin in direction v intersects plane defined
123 <  * by normal plane_n, and plane_d. If plane is not parallel- returns
124 <  * intersection point if r != NULL. If tptr!= NULL returns value of
125 <  * t, if parallel, returns t=FHUGE
126 <  */
127 < int
128 < intersect_vector_plane(v,plane_n,plane_d,tptr,r)
129 <   FVECT v,plane_n;
130 <   double plane_d;
131 <   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,d;
162    int hit;
136 –    /*
137 –      Plane is Ax + By + Cz +D = 0:
138 –      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
139 –    */
163  
164 <    /* line is  l = p1 + (p2-p1)t, p1=origin */
142 <
143 <    /* Solve for t: */
144 <  d = -(DOT(plane_n,v));
164 >  d = DOT(FP_N(peq),dir);
165    if(ZERO(d))
166 <  {
167 <      t = FHUGE;
148 <      hit = 0;
149 <  }
150 <  else
151 <  {
152 <      t =  plane_d/d;
153 <      if(t < 0 )
154 <         hit = 0;
155 <      else
156 <         hit = 1;
157 <      if(r)
158 <         {
159 <             r[0] = v[0]*t;
160 <             r[1] = v[1]*t;
161 <             r[2] = v[2]*t;
162 <         }
163 <  }
164 <    if(tptr)
165 <       *tptr = t;
166 <  return(hit);
167 < }
166 >     return(0);
167 >  t =  -(DOT(FP_N(peq),orig) + FP_D(peq))/d;
168  
169 < int
170 < intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
171 <   FVECT orig,dir;
172 <   FVECT plane_n;
173 <   double plane_d;
174 <   double *pd;
175 <   FVECT r;
176 < {
177 <  double t;
178 <  int hit;
179 <    /*
180 <      Plane is Ax + By + Cz +D = 0:
181 <      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
182 <    */
183 <     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
184 <         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
185 <       line is  l = p1 + (p2-p1)t
186 <     */
187 <    /* Solve for t: */
188 <    t =  -(DOT(plane_n,orig) + plane_d)/(DOT(plane_n,dir));
189 <    if(t < 0)
169 >  if(t < 0)
170         hit = 0;
171      else
172         hit = 1;
193 –
173    if(r)
174       VSUM(r,orig,dir,t);
175  
# Line 199 | Line 178 | intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
178    return(hit);
179   }
180  
181 <
182 < int
183 < intersect_edge_plane(e0,e1,plane_n,plane_d,pd,r)
184 <   FVECT e0,e1;
185 <   FVECT plane_n;
186 <   double plane_d;
187 <   double *pd;
188 <   FVECT r;
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 <  double t;
192 <  int hit;
193 <  FVECT d;
194 <  /*
195 <      Plane is Ax + By + Cz +D = 0:
216 <      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
217 <    */
218 <     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
219 <         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
220 <       line is  l = p1 + (p2-p1)t
221 <     */
222 <    /* Solve for t: */
223 <  VSUB(d,e1,e0);
224 <  t =  -(DOT(plane_n,e0) + plane_d)/(DOT(plane_n,d));
225 <    if(t < 0)
226 <       hit = 0;
227 <    else
228 <       hit = 1;
229 <
230 <  VSUM(r,e0,d,t);
231 <
232 <    if(pd)
233 <       *pd = t;
234 <  return(hit);
191 >  double d;
192 >  
193 >  VSUB(ps,p,c);    
194 >  d = normalize(ps);
195 >  return(d);
196   }
197  
198 <
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
239 – point_in_cone(p,p0,p1,p2)
240 – FVECT p;
241 – FVECT p0,p1,p2;
242 – {
243 –    FVECT n;
244 –    FVECT np,x_axis,y_axis;
245 –    double d1,d2,d;
246 –    
247 –    /* Find the equation of the circle defined by the intersection
248 –       of the cone with the plane defined by p1,p2,p3- project p into
249 –       that plane and do an in-circle test in the plane
250 –     */
251 –    
252 –    /* find the equation of the plane defined by p1-p3 */
253 –    tri_plane_equation(p0,p1,p2,n,&d,FALSE);
254 –
255 –    /* define a coordinate system on the plane: the x axis is in
256 –       the direction of np2-np1, and the y axis is calculated from
257 –       n cross x-axis
258 –     */
259 –    /* Project p onto the plane */
260 –    /* NOTE: check this: does sideness check?*/
261 –    if(!intersect_vector_plane(p,n,d,NULL,np))
262 –        return(FALSE);
263 –
264 –    /* create coordinate system on  plane: p2-p1 defines the x_axis*/
265 –    VSUB(x_axis,p1,p0);
266 –    normalize(x_axis);
267 –    /* The y axis is  */
268 –    VCROSS(y_axis,n,x_axis);
269 –    normalize(y_axis);
270 –
271 –    VSUB(p1,p1,p0);
272 –    VSUB(p2,p2,p0);
273 –    VSUB(np,np,p0);
274 –    
275 –    p1[0] = VLEN(p1);
276 –    p1[1] = 0;
277 –
278 –    d1 = DOT(p2,x_axis);
279 –    d2 = DOT(p2,y_axis);
280 –    p2[0] = d1;
281 –    p2[1] = d2;
282 –    
283 –    d1 = DOT(np,x_axis);
284 –    d2 = DOT(np,y_axis);
285 –    np[0] = d1;
286 –    np[1] = d2;
287 –
288 –    /* perform the in-circle test in the new coordinate system */
289 –    return(point_in_circle_thru_origin(np,p1,p2));
290 – }
291 –
292 – int
293 – point_set_in_stri(v0,v1,v2,p,n,nset,sides)
294 – FVECT v0,v1,v2,p;
295 – FVECT n[3];
296 – int *nset;
297 – int sides[3];
298 –
299 – {
300 –    double d;
301 –    /* Find the normal to the triangle ORIGIN:v0:v1 */
302 –    if(!NTH_BIT(*nset,0))
303 –    {
304 –        VCROSS(n[0],v1,v0);
305 –        SET_NTH_BIT(*nset,0);
306 –    }
307 –    /* Test the point for sidedness */
308 –    d  = DOT(n[0],p);
309 –
310 –    if(d > 0.0)
311 –     {
312 –       sides[0] =  GT_OUT;
313 –       sides[1] = sides[2] = GT_INVALID;
314 –       return(FALSE);
315 –      }
316 –    else
317 –       sides[0] = GT_INTERIOR;
318 –      
319 –    /* Test next edge */
320 –    if(!NTH_BIT(*nset,1))
321 –    {
322 –        VCROSS(n[1],v2,v1);
323 –        SET_NTH_BIT(*nset,1);
324 –    }
325 –    /* Test the point for sidedness */
326 –    d  = DOT(n[1],p);
327 –    if(d > 0.0)
328 –    {
329 –        sides[1] = GT_OUT;
330 –        sides[2] = GT_INVALID;
331 –        return(FALSE);
332 –    }
333 –    else
334 –       sides[1] = GT_INTERIOR;
335 –    /* Test next edge */
336 –    if(!NTH_BIT(*nset,2))
337 –    {
338 –        VCROSS(n[2],v0,v2);
339 –        SET_NTH_BIT(*nset,2);
340 –    }
341 –    /* Test the point for sidedness */
342 –    d  = DOT(n[2],p);
343 –    if(d > 0.0)
344 –    {
345 –      sides[2] = GT_OUT;
346 –      return(FALSE);
347 –    }
348 –    else
349 –       sides[2] = GT_INTERIOR;
350 –    /* Must be interior to the pyramid */
351 –    return(GT_INTERIOR);
352 – }
353 –
354 –
355 –
356 –
357 – int
208   point_in_stri(v0,v1,v2,p)
209   FVECT v0,v1,v2,p;
210   {
211      double d;
212      FVECT n;  
213  
214 <    VCROSS(n,v1,v0);
214 >    VCROSS(n,v0,v1);
215      /* Test the point for sidedness */
216      d  = DOT(n,p);
217      if(d > 0.0)
218        return(FALSE);
369 –
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(d > 0.0)
224         return(FALSE);
376 –
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(d > 0.0)
230         return(FALSE);
231      /* Must be interior to the pyramid */
232 <    return(GT_INTERIOR);
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
388 – vertices_in_stri(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides)
389 – FVECT t0,t1,t2,p0,p1,p2;
390 – int *nset;
391 – FVECT n[3];
392 – FVECT avg;
393 – int pt_sides[3][3];
394 –
395 – {
396 –    int below_plane[3],test;
397 –
398 –    SUM_3VEC3(avg,t0,t1,t2);
399 –    *nset = 0;
400 –    /* Test vertex v[i] against triangle j*/
401 –    /* Check if v[i] lies below plane defined by avg of 3 vectors
402 –       defining triangle
403 –       */
404 –
405 –    /* test point 0 */
406 –    if(DOT(avg,p0) < 0.0)
407 –      below_plane[0] = 1;
408 –    else
409 –      below_plane[0] = 0;
410 –    /* Test if b[i] lies in or on triangle a */
411 –    test = point_set_in_stri(t0,t1,t2,p0,n,nset,pt_sides[0]);
412 –    /* If pts[i] is interior: done */
413 –    if(!below_plane[0])
414 –      {
415 –        if(test == GT_INTERIOR)
416 –          return(TRUE);
417 –      }
418 –    /* Now test point 1*/
419 –
420 –    if(DOT(avg,p1) < 0.0)
421 –      below_plane[1] = 1;
422 –    else
423 –      below_plane[1]=0;
424 –    /* Test if b[i] lies in or on triangle a */
425 –    test = point_set_in_stri(t0,t1,t2,p1,n,nset,pt_sides[1]);
426 –    /* If pts[i] is interior: done */
427 –    if(!below_plane[1])
428 –    {
429 –      if(test == GT_INTERIOR)
430 –        return(TRUE);
431 –    }
432 –    
433 –    /* Now test point 2 */
434 –    if(DOT(avg,p2) < 0.0)
435 –      below_plane[2] = 1;
436 –    else
437 –      below_plane[2] = 0;
438 –        /* Test if b[i] lies in or on triangle a */
439 –    test = point_set_in_stri(t0,t1,t2,p2,n,nset,pt_sides[2]);
440 –
441 –    /* If pts[i] is interior: done */
442 –    if(!below_plane[2])
443 –      {
444 –        if(test == GT_INTERIOR)
445 –          return(TRUE);
446 –      }
447 –
448 –    /* If all three points below separating plane: trivial reject */
449 –    if(below_plane[0] && below_plane[1] && below_plane[2])
450 –       return(FALSE);
451 –    /* Now check vertices in a against triangle b */
452 –    return(FALSE);
453 – }
454 –
455 –
456 – set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,n)
457 –   FVECT t0,t1,t2,p0,p1,p2;
458 –   int test[3];
459 –   int sides[3][3];
460 –   int nset;
461 –   FVECT n[3];
462 – {
463 –    int t;
464 –    double d;
465 –
466 –    
467 –    /* p=0 */
468 –    test[0] = 0;
469 –    if(sides[0][0] == GT_INVALID)
470 –    {
471 –      if(!NTH_BIT(nset,0))
472 –        VCROSS(n[0],t1,t0);
473 –      /* Test the point for sidedness */
474 –      d  = DOT(n[0],p0);
475 –      if(d >= 0.0)
476 –        SET_NTH_BIT(test[0],0);
477 –    }
478 –    else
479 –      if(sides[0][0] != GT_INTERIOR)
480 –        SET_NTH_BIT(test[0],0);
481 –
482 –    if(sides[0][1] == GT_INVALID)
483 –    {
484 –      if(!NTH_BIT(nset,1))
485 –        VCROSS(n[1],t2,t1);
486 –        /* Test the point for sidedness */
487 –        d  = DOT(n[1],p0);
488 –        if(d >= 0.0)
489 –          SET_NTH_BIT(test[0],1);
490 –    }
491 –    else
492 –      if(sides[0][1] != GT_INTERIOR)
493 –        SET_NTH_BIT(test[0],1);
494 –
495 –    if(sides[0][2] == GT_INVALID)
496 –    {
497 –      if(!NTH_BIT(nset,2))
498 –        VCROSS(n[2],t0,t2);
499 –      /* Test the point for sidedness */
500 –      d  = DOT(n[2],p0);
501 –      if(d >= 0.0)
502 –        SET_NTH_BIT(test[0],2);
503 –    }
504 –    else
505 –      if(sides[0][2] != GT_INTERIOR)
506 –        SET_NTH_BIT(test[0],2);
507 –
508 –    /* p=1 */
509 –    test[1] = 0;
510 –    /* t=0*/
511 –    if(sides[1][0] == GT_INVALID)
512 –    {
513 –      if(!NTH_BIT(nset,0))
514 –        VCROSS(n[0],t1,t0);
515 –      /* Test the point for sidedness */
516 –      d  = DOT(n[0],p1);
517 –      if(d >= 0.0)
518 –        SET_NTH_BIT(test[1],0);
519 –    }
520 –    else
521 –      if(sides[1][0] != GT_INTERIOR)
522 –        SET_NTH_BIT(test[1],0);
523 –    
524 –    /* t=1 */
525 –    if(sides[1][1] == GT_INVALID)
526 –    {
527 –      if(!NTH_BIT(nset,1))
528 –        VCROSS(n[1],t2,t1);
529 –      /* Test the point for sidedness */
530 –      d  = DOT(n[1],p1);
531 –      if(d >= 0.0)
532 –        SET_NTH_BIT(test[1],1);
533 –    }
534 –    else
535 –      if(sides[1][1] != GT_INTERIOR)
536 –        SET_NTH_BIT(test[1],1);
537 –      
538 –    /* t=2 */
539 –    if(sides[1][2] == GT_INVALID)
540 –    {
541 –      if(!NTH_BIT(nset,2))
542 –        VCROSS(n[2],t0,t2);
543 –      /* Test the point for sidedness */
544 –      d  = DOT(n[2],p1);
545 –      if(d >= 0.0)
546 –        SET_NTH_BIT(test[1],2);
547 –    }
548 –    else
549 –      if(sides[1][2] != GT_INTERIOR)
550 –        SET_NTH_BIT(test[1],2);
551 –
552 –    /* p=2 */
553 –    test[2] = 0;
554 –    /* t = 0 */
555 –    if(sides[2][0] == GT_INVALID)
556 –    {
557 –      if(!NTH_BIT(nset,0))
558 –        VCROSS(n[0],t1,t0);
559 –      /* Test the point for sidedness */
560 –      d  = DOT(n[0],p2);
561 –      if(d >= 0.0)
562 –        SET_NTH_BIT(test[2],0);
563 –    }
564 –    else
565 –      if(sides[2][0] != GT_INTERIOR)
566 –        SET_NTH_BIT(test[2],0);
567 –    /* t=1 */
568 –    if(sides[2][1] == GT_INVALID)
569 –    {
570 –      if(!NTH_BIT(nset,1))
571 –        VCROSS(n[1],t2,t1);
572 –      /* Test the point for sidedness */
573 –      d  = DOT(n[1],p2);
574 –      if(d >= 0.0)
575 –        SET_NTH_BIT(test[2],1);
576 –    }
577 –    else
578 –      if(sides[2][1] != GT_INTERIOR)
579 –        SET_NTH_BIT(test[2],1);
580 –    /* t=2 */
581 –    if(sides[2][2] == GT_INVALID)
582 –    {
583 –      if(!NTH_BIT(nset,2))
584 –        VCROSS(n[2],t0,t2);
585 –      /* Test the point for sidedness */
586 –      d  = DOT(n[2],p2);
587 –      if(d >= 0.0)
588 –        SET_NTH_BIT(test[2],2);
589 –    }
590 –    else
591 –      if(sides[2][2] != GT_INTERIOR)
592 –        SET_NTH_BIT(test[2],2);
593 – }
594 –
595 –
596 – int
597 – stri_intersect(a1,a2,a3,b1,b2,b3)
598 – FVECT a1,a2,a3,b1,b2,b3;
599 – {
600 –  int which,test,n_set[2];
601 –  int sides[2][3][3],i,j,inext,jnext;
602 –  int tests[2][3];
603 –  FVECT n[2][3],p,avg[2];
604 –
605 –  /* Test the vertices of triangle a against the pyramid formed by triangle
606 –     b and the origin. If any vertex of a is interior to triangle b, or
607 –     if all 3 vertices of a are ON the edges of b,return TRUE. Remember
608 –     the results of the edge normal and sidedness tests for later.
609 –   */
610 – if(vertices_in_stri(a1,a2,a3,b1,b2,b3,&(n_set[0]),n[0],avg[0],sides[1]))
611 –     return(TRUE);
612 –  
613 – if(vertices_in_stri(b1,b2,b3,a1,a2,a3,&(n_set[1]),n[1],avg[1],sides[0]))
614 –     return(TRUE);
615 –
616 –
617 –  set_sidedness_tests(b1,b2,b3,a1,a2,a3,tests[0],sides[0],n_set[1],n[1]);
618 –  if(tests[0][0]&tests[0][1]&tests[0][2])
619 –    return(FALSE);
620 –
621 –  set_sidedness_tests(a1,a2,a3,b1,b2,b3,tests[1],sides[1],n_set[0],n[0]);
622 –  if(tests[1][0]&tests[1][1]&tests[1][2])
623 –    return(FALSE);
624 –
625 –  for(j=0; j < 3;j++)
626 –  {
627 –    jnext = (j+1)%3;
628 –    /* IF edge b doesnt cross any great circles of a, punt */
629 –    if(tests[1][j] & tests[1][jnext])
630 –      continue;
631 –    for(i=0;i<3;i++)
632 –    {
633 –      inext = (i+1)%3;
634 –      /* IF edge a doesnt cross any great circles of b, punt */
635 –      if(tests[0][i] & tests[0][inext])
636 –        continue;
637 –      /* Now find the great circles that cross and test */
638 –      if((NTH_BIT(tests[0][i],j)^(NTH_BIT(tests[0][inext],j)))
639 –          && (NTH_BIT(tests[1][j],i)^NTH_BIT(tests[1][jnext],i)))
640 –      {
641 –        VCROSS(p,n[0][i],n[1][j]);
642 –                    
643 –        /* If zero cp= done */
644 –        if(ZERO_VEC3(p))
645 –          continue;
646 –        /* check above both planes */
647 –        if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0)
648 –          {
649 –            NEGATE_VEC3(p);
650 –            if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0)
651 –              continue;
652 –          }
653 –        return(TRUE);
654 –      }
655 –    }
656 –  }
657 –  return(FALSE);
658 – }
659 –
660 – int
244   ray_intersect_tri(orig,dir,v0,v1,v2,pt)
245   FVECT orig,dir;
246   FVECT v0,v1,v2;
247   FVECT pt;
248   {
249 <  FVECT p0,p1,p2,p,n;
250 <  double pd;
249 >  FVECT p0,p1,p2,p;
250 >  FPEQ peq;
251    int type;
252  
253    VSUB(p0,v0,orig);
# Line 674 | Line 257 | FVECT pt;
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 <      return(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    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 690 | 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 */
693 –    /* hv and vv are the horizontal and vertical vectors in the
694 –       view frame-the magnitude is the dimension of the front frustum
695 –       face at z =1
696 –    */
288      VCOPY(nhv,hv);
289      VCOPY(nvv,vv);
290      w2 = normalize(nhv);
# Line 737 | Line 328 | FVECT fnear[4],ffar[4];
328   }
329  
330  
331 <
332 <
333 < int
334 < sedge_intersect(a0,a1,b0,b1)
335 < FVECT a0,a1,b0,b1;
336 < {
337 <    FVECT na,nb,avga,avgb,p;
747 <    double d;
748 <    int sb0,sb1,sa0,sa1;
749 <
750 <    /* First test if edge b straddles great circle of a */
751 <    VCROSS(na,a0,a1);
752 <    d = DOT(na,b0);
753 <    sb0 = ZERO(d)?0:(d<0)? -1: 1;
754 <    d = DOT(na,b1);
755 <    sb1 = ZERO(d)?0:(d<0)? -1: 1;
756 <    /* edge b entirely on one side of great circle a: edges cannot intersect*/
757 <    if(sb0*sb1 > 0)
758 <       return(FALSE);
759 <    /* test if edge a straddles great circle of b */
760 <    VCROSS(nb,b0,b1);
761 <    d = DOT(nb,a0);
762 <    sa0 = ZERO(d)?0:(d<0)? -1: 1;
763 <    d = DOT(nb,a1);
764 <    sa1 = ZERO(d)?0:(d<0)? -1: 1;
765 <    /* edge a entirely on one side of great circle b: edges cannot intersect*/
766 <    if(sa0*sa1 > 0)
767 <       return(FALSE);
768 <
769 <    /* Find one of intersection points of the great circles */
770 <    VCROSS(p,na,nb);
771 <    /* If they lie on same great circle: call an intersection */
772 <    if(ZERO_VEC3(p))
773 <       return(TRUE);
774 <
775 <    VADD(avga,a0,a1);
776 <    VADD(avgb,b0,b1);
777 <    if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
778 <    {
779 <      NEGATE_VEC3(p);
780 <      if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
781 <        return(FALSE);
782 <    }
783 <    if((!sb0 || !sb1) && (!sa0 || !sa1))
784 <      return(FALSE);
785 <    return(TRUE);
786 < }
787 <
788 <
789 <
790 < /* Find the normalized barycentric coordinates of p relative to
791 < * 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 800 | 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  
807 – bary_ith_child(coord,i)
808 – double coord[3];
809 – int i;
810 – {
811 –
812 –  switch(i){
813 –  case 0:
814 –      /* update bary for child */
815 –      coord[0] = 2.0*coord[0]- 1.0;
816 –      coord[1] *= 2.0;
817 –      coord[2] *= 2.0;
818 –      break;
819 –  case 1:
820 –    coord[0] *= 2.0;
821 –    coord[1] = 2.0*coord[1]- 1.0;
822 –    coord[2] *= 2.0;
823 –    break;
824 –  case 2:
825 –    coord[0] *= 2.0;
826 –    coord[1] *= 2.0;
827 –    coord[2] = 2.0*coord[2]- 1.0;
828 –    break;
829 –  case 3:
830 –    coord[0] = 1.0 - 2.0*coord[0];
831 –    coord[1] = 1.0 - 2.0*coord[1];
832 –    coord[2] = 1.0 - 2.0*coord[2];
833 –    break;
834 – #ifdef DEBUG
835 –  default:
836 –    eputs("bary_ith_child():Invalid child\n");
837 –    break;
838 – #endif
839 –  }
840 – }
353  
354  
843 – int
844 – bary_child(coord)
845 – double coord[3];
846 – {
847 –  int i;
355  
849 –  if(coord[0] > 0.5)
850 –  {
851 –      /* update bary for child */
852 –      coord[0] = 2.0*coord[0]- 1.0;
853 –      coord[1] *= 2.0;
854 –      coord[2] *= 2.0;
855 –      return(0);
856 –  }
857 –  else
858 –    if(coord[1] > 0.5)
859 –    {
860 –      coord[0] *= 2.0;
861 –      coord[1] = 2.0*coord[1]- 1.0;
862 –      coord[2] *= 2.0;
863 –      return(1);
864 –    }
865 –    else
866 –      if(coord[2] > 0.5)
867 –      {
868 –        coord[0] *= 2.0;
869 –        coord[1] *= 2.0;
870 –        coord[2] = 2.0*coord[2]- 1.0;
871 –        return(2);
872 –      }
873 –      else
874 –         {
875 –           coord[0] = 1.0 - 2.0*coord[0];
876 –           coord[1] = 1.0 - 2.0*coord[1];
877 –           coord[2] = 1.0 - 2.0*coord[2];
878 –           return(3);
879 –         }
880 – }
356  
882 – /* Coord was the ith child of the parent: set the coordinate
883 –   relative to the parent
884 – */
885 – bary_parent(coord,i)
886 – double coord[3];
887 – int i;
888 – {
357  
890 –  switch(i) {
891 –  case 0:
892 –    /* update bary for child */
893 –    coord[0] = coord[0]*0.5 + 0.5;
894 –    coord[1] *= 0.5;
895 –    coord[2] *= 0.5;
896 –    break;
897 –  case 1:
898 –    coord[0] *= 0.5;
899 –    coord[1]  = coord[1]*0.5 + 0.5;
900 –    coord[2] *= 0.5;
901 –    break;
902 –    
903 –  case 2:
904 –    coord[0] *= 0.5;
905 –    coord[1] *= 0.5;
906 –    coord[2] = coord[2]*0.5 + 0.5;
907 –    break;
908 –    
909 –  case 3:
910 –    coord[0] = 0.5 - 0.5*coord[0];
911 –    coord[1] = 0.5 - 0.5*coord[1];
912 –    coord[2] = 0.5 - 0.5*coord[2];
913 –    break;
914 – #ifdef DEBUG
915 –  default:
916 –    eputs("bary_parent():Invalid child\n");
917 –    break;
918 – #endif
919 –  }
920 – }
358  
922 – bary_from_child(coord,child,next)
923 – double coord[3];
924 – int child,next;
925 – {
926 – #ifdef DEBUG
927 –  if(child <0 || child > 3)
928 –  {
929 –    eputs("bary_from_child():Invalid child\n");
930 –    return;
931 –  }
932 –  if(next <0 || next > 3)
933 –  {
934 –    eputs("bary_from_child():Invalid next\n");
935 –    return;
936 –  }
937 – #endif
938 –  if(next == child)
939 –    return;
359  
941 –  switch(child){
942 –  case 0:
943 –    switch(next){
944 –    case 1:
945 –      coord[0] += 1.0;
946 –      coord[1] -= 1.0;
947 –      break;
948 –    case 2:
949 –      coord[0] += 1.0;
950 –      coord[2] -= 1.0;
951 –      break;
952 –    case 3:
953 –      coord[0] *= -1.0;
954 –      coord[1] = 1 - coord[1];
955 –      coord[2] = 1 - coord[2];
956 –      break;
957 –
958 –    }
959 –    break;
960 –  case 1:
961 –    switch(next){
962 –    case 0:
963 –      coord[0] -= 1.0;
964 –      coord[1] += 1.0;
965 –      break;
966 –    case 2:
967 –      coord[1] += 1.0;
968 –      coord[2] -= 1.0;
969 –      break;
970 –    case 3:
971 –      coord[0] = 1 - coord[0];
972 –      coord[1] *= -1.0;
973 –      coord[2] = 1 - coord[2];
974 –      break;
975 –    }
976 –    break;
977 –  case 2:
978 –    switch(next){
979 –    case 0:
980 –      coord[0] -= 1.0;
981 –      coord[2] += 1.0;
982 –      break;
983 –    case 1:
984 –      coord[1] -= 1.0;
985 –      coord[2] += 1.0;
986 –      break;
987 –    case 3:
988 –      coord[0] = 1 - coord[0];
989 –      coord[1] = 1 - coord[1];
990 –      coord[2] *= -1.0;
991 –      break;
992 –    }
993 –    break;
994 –  case 3:
995 –    switch(next){
996 –    case 0:
997 –      coord[0] *= -1.0;
998 –      coord[1] = 1 - coord[1];
999 –      coord[2] = 1 - coord[2];
1000 –      break;
1001 –    case 1:
1002 –      coord[0] = 1 - coord[0];
1003 –      coord[1] *= -1.0;
1004 –      coord[2] = 1 - coord[2];
1005 –      break;
1006 –    case 2:
1007 –      coord[0] = 1 - coord[0];
1008 –      coord[1] = 1 - coord[1];
1009 –      coord[2] *= -1.0;
1010 –      break;
1011 –    }
1012 –    break;
1013 –  }
1014 – }
1015 –
1016 – int
1017 – max_index(v)
1018 – FVECT v;
1019 – {
1020 –  double a,b,c;
1021 –  int i;
1022 –
1023 –  a = fabs(v[0]);
1024 –  b = fabs(v[1]);
1025 –  c = fabs(v[2]);
1026 –  i = (a>=b)?((a>=c)?0:2):((b>=c)?1:2);  
1027 –  return(i);
1028 – }
1029 –
1030 – int
1031 – closest_point_in_tri(p0,p1,p2,p,p0id,p1id,p2id)
1032 – FVECT p0,p1,p2,p;
1033 – int p0id,p1id,p2id;
1034 – {
1035 –    double d,d1;
1036 –    int i;
1037 –    
1038 –    d =  DIST_SQ(p,p0);
1039 –    d1 = DIST_SQ(p,p1);
1040 –    if(d < d1)
1041 –    {
1042 –      d1 = DIST_SQ(p,p2);
1043 –      i = (d1 < d)?p2id:p0id;
1044 –    }
1045 –    else
1046 –    {
1047 –      d = DIST_SQ(p,p2);
1048 –      i = (d < d1)? p2id:p1id;
1049 –    }
1050 –    return(i);
1051 – }
360  
361  
362  

Diff Legend

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