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.11 by gwlarson, Sun Jan 10 10:27:47 1999 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   }
30 – #if 0
31 – extern FVECT Norm[500];
32 – extern int Ncnt;
33 – #endif
70  
71 < int
72 < convex_angle(v0,v1,v2)
73 < 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 cp,cp01,cp12,v10,v02;
80 <    double dp;
81 <
42 <    /* test sign of (v0Xv1)X(v1Xv2). v1 */
43 <    VCROSS(cp01,v0,v1);
44 <    VCROSS(cp12,v1,v2);
45 <    VCROSS(cp,cp01,cp12);
46 <        
47 <    dp = DOT(cp,v1);
48 < #if 0
49 <    VCOPY(Norm[Ncnt++],cp01);
50 <    VCOPY(Norm[Ncnt++],cp12);
51 <    VCOPY(Norm[Ncnt++],cp);
52 < #endif
53 <    if(ZERO(dp) || dp < 0.0)
54 <      return(FALSE);
55 <    return(TRUE);
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  
58 – /* 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 71 | 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]);
74 –  
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]);
78 –  
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]);
82 –
113    if(!norm)
114       return(0);
85 –  
115    mag = normalize(n);
87 –
116    return(mag);
117   }
118  
119 <
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;
128 >   FVECT v0,v1,v2;
129     FPEQ *peqptr;
130     int norm;
131   {
132      tri_normal(v0,v1,v2,FP_N(*peqptr),norm);
98 –
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
104 – point_in_circle_thru_origin(p,p0,p1)
105 – FVECT p;
106 – FVECT p0,p1;
107 – {
108 –
109 –    double dp0,dp1;
110 –    double dp,det;
111 –    
112 –    dp0 = DOT_VEC2(p0,p0);
113 –    dp1 = DOT_VEC2(p1,p1);
114 –    dp  = DOT_VEC2(p,p);    
115 –
116 –    det = -dp0*CROSS_VEC2(p1,p) + dp1*CROSS_VEC2(p0,p) - dp*CROSS_VEC2(p0,p1);
117 –    
118 –    return (det > (1e-10));
119 – }
120 –
121 –
122 – double
123 – point_on_sphere(ps,p,c)
124 – FVECT ps,p,c;
125 – {
126 –  double d;
127 –    VSUB(ps,p,c);    
128 –    d= normalize(ps);
129 –    return(d);
130 – }
131 –
132 –
133 – /* returns TRUE if ray from origin in direction v intersects plane defined
134 –  * by normal plane_n, and plane_d. If plane is not parallel- returns
135 –  * intersection point if r != NULL. If tptr!= NULL returns value of
136 –  * t, if parallel, returns t=FHUGE
137 –  */
138 – int
139 – intersect_vector_plane(v,peq,tptr,r)
140 –   FVECT v;
141 –   FPEQ peq;
142 –   double *tptr;
143 –   FVECT r;
144 – {
145 –  double t,d;
146 –  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 –    */
151 –
152 –    /* line is  l = p1 + (p2-p1)t, p1=origin */
153 –
154 –    /* Solve for t: */
155 –  d = -(DOT(FP_N(peq),v));
156 –  if(ZERO(d))
157 –  {
158 –      t = FHUGE;
159 –      hit = 0;
160 –  }
161 –  else
162 –  {
163 –      t =  FP_D(peq)/d;
164 –      if(t < 0 )
165 –         hit = 0;
166 –      else
167 –         hit = 1;
168 –      if(r)
169 –         {
170 –             r[0] = v[0]*t;
171 –             r[1] = v[1]*t;
172 –             r[2] = v[2]*t;
173 –         }
174 –  }
175 –    if(tptr)
176 –       *tptr = t;
177 –  return(hit);
178 – }
179 –
180 – int
155   intersect_ray_plane(orig,dir,peq,pd,r)
156     FVECT orig,dir;
157     FPEQ peq;
# Line 186 | Line 160 | intersect_ray_plane(orig,dir,peq,pd,r)
160   {
161    double t,d;
162    int hit;
163 <    /*
190 <      Plane is Ax + By + Cz +D = 0:
191 <      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
192 <    */
193 <     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
194 <         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
195 <       line is  l = p1 + (p2-p1)t
196 <     */
197 <    /* Solve for t: */
163 >
164    d = DOT(FP_N(peq),dir);
165    if(ZERO(d))
166       return(0);
# Line 204 | Line 170 | intersect_ray_plane(orig,dir,peq,pd,r)
170         hit = 0;
171      else
172         hit = 1;
207 –
173    if(r)
174       VSUM(r,orig,dir,t);
175  
# Line 213 | Line 178 | intersect_ray_plane(orig,dir,peq,pd,r)
178    return(hit);
179   }
180  
181 <
182 < int
183 < intersect_ray_oplane(orig,dir,n,pd,r)
184 <   FVECT orig,dir;
185 <   FVECT n;
186 <   double *pd;
187 <   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,d;
192 <  int hit;
193 <    /*
194 <      Plane is Ax + By + Cz +D = 0:
195 <      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
229 <    */
230 <     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
231 <         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
232 <       line is  l = p1 + (p2-p1)t
233 <     */
234 <    /* Solve for t: */
235 <    d= DOT(n,dir);
236 <    if(ZERO(d))
237 <       return(0);
238 <    t =  -(DOT(n,orig))/d;
239 <    if(t < 0)
240 <       hit = 0;
241 <    else
242 <       hit = 1;
243 <
244 <  if(r)
245 <     VSUM(r,orig,dir,t);
246 <
247 <    if(pd)
248 <       *pd = t;
249 <  return(hit);
191 >  double d;
192 >  
193 >  VSUB(ps,p,c);    
194 >  d = normalize(ps);
195 >  return(d);
196   }
197  
198 < /* Assumption: know crosses plane:dont need to check for 'on' case */
199 < intersect_edge_coord_plane(v0,v1,w,r)
200 < FVECT v0,v1;
201 < int w;
202 < FVECT r;
203 < {
204 <  FVECT dv;
205 <  int wnext;
206 <  double t;
261 <
262 <  VSUB(dv,v1,v0);
263 <  t = -v0[w]/dv[w];
264 <  r[w] = 0.0;
265 <  wnext = (w+1)%3;
266 <  r[wnext] = v0[wnext] + dv[wnext]*t;
267 <  wnext = (w+2)%3;
268 <  r[wnext] = v0[wnext] + dv[wnext]*t;
269 < }
270 <
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
272 – intersect_edge_plane(e0,e1,peq,pd,r)
273 –   FVECT e0,e1;
274 –   FPEQ peq;
275 –   double *pd;
276 –   FVECT r;
277 – {
278 –  double t;
279 –  int hit;
280 –  FVECT d;
281 –  /*
282 –      Plane is Ax + By + Cz +D = 0:
283 –      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
284 –    */
285 –     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
286 –         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
287 –       line is  l = p1 + (p2-p1)t
288 –     */
289 –    /* Solve for t: */
290 –  VSUB(d,e1,e0);
291 –  t =  -(DOT(FP_N(peq),e0) + FP_D(peq))/(DOT(FP_N(peq),d));
292 –    if(t < 0)
293 –       hit = 0;
294 –    else
295 –       hit = 1;
296 –
297 –  VSUM(r,e0,d,t);
298 –
299 –    if(pd)
300 –       *pd = t;
301 –  return(hit);
302 – }
303 –
304 –
305 – int
306 – point_in_cone(p,p0,p1,p2)
307 – FVECT p;
308 – FVECT p0,p1,p2;
309 – {
310 –    FVECT np,x_axis,y_axis;
311 –    double d1,d2;
312 –    FPEQ peq;
313 –    
314 –    /* Find the equation of the circle defined by the intersection
315 –       of the cone with the plane defined by p1,p2,p3- project p into
316 –       that plane and do an in-circle test in the plane
317 –     */
318 –    
319 –    /* find the equation of the plane defined by p0-p2 */
320 –    tri_plane_equation(p0,p1,p2,&peq,FALSE);
321 –
322 –    /* define a coordinate system on the plane: the x axis is in
323 –       the direction of np2-np1, and the y axis is calculated from
324 –       n cross x-axis
325 –     */
326 –    /* Project p onto the plane */
327 –    /* NOTE: check this: does sideness check?*/
328 –    if(!intersect_vector_plane(p,peq,NULL,np))
329 –        return(FALSE);
330 –
331 –    /* create coordinate system on  plane: p1-p0 defines the x_axis*/
332 –    VSUB(x_axis,p1,p0);
333 –    normalize(x_axis);
334 –    /* The y axis is  */
335 –    VCROSS(y_axis,FP_N(peq),x_axis);
336 –    normalize(y_axis);
337 –
338 –    VSUB(p1,p1,p0);
339 –    VSUB(p2,p2,p0);
340 –    VSUB(np,np,p0);
341 –    
342 –    p1[0] = DOT(p1,x_axis);
343 –    p1[1] = 0;
344 –
345 –    d1 = DOT(p2,x_axis);
346 –    d2 = DOT(p2,y_axis);
347 –    p2[0] = d1;
348 –    p2[1] = d2;
349 –    
350 –    d1 = DOT(np,x_axis);
351 –    d2 = DOT(np,y_axis);
352 –    np[0] = d1;
353 –    np[1] = d2;
354 –
355 –    /* perform the in-circle test in the new coordinate system */
356 –    return(point_in_circle_thru_origin(np,p1,p2));
357 – }
358 –
359 – int
360 – point_set_in_stri(v0,v1,v2,p,n,nset,sides)
361 – FVECT v0,v1,v2,p;
362 – FVECT n[3];
363 – int *nset;
364 – int sides[3];
365 –
366 – {
367 –    double d;
368 –    /* Find the normal to the triangle ORIGIN:v0:v1 */
369 –    if(!NTH_BIT(*nset,0))
370 –    {
371 –        VCROSS(n[0],v0,v1);
372 –        SET_NTH_BIT(*nset,0);
373 –    }
374 –    /* Test the point for sidedness */
375 –    d  = DOT(n[0],p);
376 –
377 –    if(d > 0.0)
378 –     {
379 –       sides[0] =  GT_OUT;
380 –       sides[1] = sides[2] = GT_INVALID;
381 –       return(FALSE);
382 –      }
383 –    else
384 –       sides[0] = GT_INTERIOR;
385 –      
386 –    /* Test next edge */
387 –    if(!NTH_BIT(*nset,1))
388 –    {
389 –        VCROSS(n[1],v1,v2);
390 –        SET_NTH_BIT(*nset,1);
391 –    }
392 –    /* Test the point for sidedness */
393 –    d  = DOT(n[1],p);
394 –    if(d > 0.0)
395 –    {
396 –        sides[1] = GT_OUT;
397 –        sides[2] = GT_INVALID;
398 –        return(FALSE);
399 –    }
400 –    else
401 –       sides[1] = GT_INTERIOR;
402 –    /* Test next edge */
403 –    if(!NTH_BIT(*nset,2))
404 –    {
405 –        VCROSS(n[2],v2,v0);
406 –        SET_NTH_BIT(*nset,2);
407 –    }
408 –    /* Test the point for sidedness */
409 –    d  = DOT(n[2],p);
410 –    if(d > 0.0)
411 –    {
412 –      sides[2] = GT_OUT;
413 –      return(FALSE);
414 –    }
415 –    else
416 –       sides[2] = GT_INTERIOR;
417 –    /* Must be interior to the pyramid */
418 –    return(GT_INTERIOR);
419 – }
420 –
421 –
422 –
423 –
424 – int
208   point_in_stri(v0,v1,v2,p)
209   FVECT v0,v1,v2,p;
210   {
# Line 433 | Line 216 | FVECT v0,v1,v2,p;
216      d  = DOT(n,p);
217      if(d > 0.0)
218        return(FALSE);
436 –
219      /* Test next edge */
220      VCROSS(n,v1,v2);
221      /* Test the point for sidedness */
222      d  = DOT(n,p);
223      if(d > 0.0)
224         return(FALSE);
443 –
225      /* Test next edge */
226      VCROSS(n,v2,v0);
227      /* Test the point for sidedness */
# Line 448 | Line 229 | FVECT v0,v1,v2,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
455 – vertices_in_stri(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides)
456 – FVECT t0,t1,t2,p0,p1,p2;
457 – int *nset;
458 – FVECT n[3];
459 – FVECT avg;
460 – int pt_sides[3][3];
461 –
462 – {
463 –    int below_plane[3],test;
464 –
465 –    SUM_3VEC3(avg,t0,t1,t2);
466 –    *nset = 0;
467 –    /* Test vertex v[i] against triangle j*/
468 –    /* Check if v[i] lies below plane defined by avg of 3 vectors
469 –       defining triangle
470 –       */
471 –
472 –    /* test point 0 */
473 –    if(DOT(avg,p0) < 0.0)
474 –      below_plane[0] = 1;
475 –    else
476 –      below_plane[0] = 0;
477 –    /* Test if b[i] lies in or on triangle a */
478 –    test = point_set_in_stri(t0,t1,t2,p0,n,nset,pt_sides[0]);
479 –    /* If pts[i] is interior: done */
480 –    if(!below_plane[0])
481 –      {
482 –        if(test == GT_INTERIOR)
483 –          return(TRUE);
484 –      }
485 –    /* Now test point 1*/
486 –
487 –    if(DOT(avg,p1) < 0.0)
488 –      below_plane[1] = 1;
489 –    else
490 –      below_plane[1]=0;
491 –    /* Test if b[i] lies in or on triangle a */
492 –    test = point_set_in_stri(t0,t1,t2,p1,n,nset,pt_sides[1]);
493 –    /* If pts[i] is interior: done */
494 –    if(!below_plane[1])
495 –    {
496 –      if(test == GT_INTERIOR)
497 –        return(TRUE);
498 –    }
499 –    
500 –    /* Now test point 2 */
501 –    if(DOT(avg,p2) < 0.0)
502 –      below_plane[2] = 1;
503 –    else
504 –      below_plane[2] = 0;
505 –        /* Test if b[i] lies in or on triangle a */
506 –    test = point_set_in_stri(t0,t1,t2,p2,n,nset,pt_sides[2]);
507 –
508 –    /* If pts[i] is interior: done */
509 –    if(!below_plane[2])
510 –      {
511 –        if(test == GT_INTERIOR)
512 –          return(TRUE);
513 –      }
514 –
515 –    /* If all three points below separating plane: trivial reject */
516 –    if(below_plane[0] && below_plane[1] && below_plane[2])
517 –       return(FALSE);
518 –    /* Now check vertices in a against triangle b */
519 –    return(FALSE);
520 – }
521 –
522 –
523 – set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,n)
524 –   FVECT t0,t1,t2,p0,p1,p2;
525 –   int test[3];
526 –   int sides[3][3];
527 –   int nset;
528 –   FVECT n[3];
529 – {
530 –    int t;
531 –    double d;
532 –
533 –    
534 –    /* p=0 */
535 –    test[0] = 0;
536 –    if(sides[0][0] == GT_INVALID)
537 –    {
538 –      if(!NTH_BIT(nset,0))
539 –        VCROSS(n[0],t0,t1);
540 –      /* Test the point for sidedness */
541 –      d  = DOT(n[0],p0);
542 –      if(d >= 0.0)
543 –        SET_NTH_BIT(test[0],0);
544 –    }
545 –    else
546 –      if(sides[0][0] != GT_INTERIOR)
547 –        SET_NTH_BIT(test[0],0);
548 –
549 –    if(sides[0][1] == GT_INVALID)
550 –    {
551 –      if(!NTH_BIT(nset,1))
552 –        VCROSS(n[1],t1,t2);
553 –        /* Test the point for sidedness */
554 –        d  = DOT(n[1],p0);
555 –        if(d >= 0.0)
556 –          SET_NTH_BIT(test[0],1);
557 –    }
558 –    else
559 –      if(sides[0][1] != GT_INTERIOR)
560 –        SET_NTH_BIT(test[0],1);
561 –
562 –    if(sides[0][2] == GT_INVALID)
563 –    {
564 –      if(!NTH_BIT(nset,2))
565 –        VCROSS(n[2],t2,t0);
566 –      /* Test the point for sidedness */
567 –      d  = DOT(n[2],p0);
568 –      if(d >= 0.0)
569 –        SET_NTH_BIT(test[0],2);
570 –    }
571 –    else
572 –      if(sides[0][2] != GT_INTERIOR)
573 –        SET_NTH_BIT(test[0],2);
574 –
575 –    /* p=1 */
576 –    test[1] = 0;
577 –    /* t=0*/
578 –    if(sides[1][0] == GT_INVALID)
579 –    {
580 –      if(!NTH_BIT(nset,0))
581 –        VCROSS(n[0],t0,t1);
582 –      /* Test the point for sidedness */
583 –      d  = DOT(n[0],p1);
584 –      if(d >= 0.0)
585 –        SET_NTH_BIT(test[1],0);
586 –    }
587 –    else
588 –      if(sides[1][0] != GT_INTERIOR)
589 –        SET_NTH_BIT(test[1],0);
590 –    
591 –    /* t=1 */
592 –    if(sides[1][1] == GT_INVALID)
593 –    {
594 –      if(!NTH_BIT(nset,1))
595 –        VCROSS(n[1],t1,t2);
596 –      /* Test the point for sidedness */
597 –      d  = DOT(n[1],p1);
598 –      if(d >= 0.0)
599 –        SET_NTH_BIT(test[1],1);
600 –    }
601 –    else
602 –      if(sides[1][1] != GT_INTERIOR)
603 –        SET_NTH_BIT(test[1],1);
604 –      
605 –    /* t=2 */
606 –    if(sides[1][2] == GT_INVALID)
607 –    {
608 –      if(!NTH_BIT(nset,2))
609 –        VCROSS(n[2],t2,t0);
610 –      /* Test the point for sidedness */
611 –      d  = DOT(n[2],p1);
612 –      if(d >= 0.0)
613 –        SET_NTH_BIT(test[1],2);
614 –    }
615 –    else
616 –      if(sides[1][2] != GT_INTERIOR)
617 –        SET_NTH_BIT(test[1],2);
618 –
619 –    /* p=2 */
620 –    test[2] = 0;
621 –    /* t = 0 */
622 –    if(sides[2][0] == GT_INVALID)
623 –    {
624 –      if(!NTH_BIT(nset,0))
625 –        VCROSS(n[0],t0,t1);
626 –      /* Test the point for sidedness */
627 –      d  = DOT(n[0],p2);
628 –      if(d >= 0.0)
629 –        SET_NTH_BIT(test[2],0);
630 –    }
631 –    else
632 –      if(sides[2][0] != GT_INTERIOR)
633 –        SET_NTH_BIT(test[2],0);
634 –    /* t=1 */
635 –    if(sides[2][1] == GT_INVALID)
636 –    {
637 –      if(!NTH_BIT(nset,1))
638 –        VCROSS(n[1],t1,t2);
639 –      /* Test the point for sidedness */
640 –      d  = DOT(n[1],p2);
641 –      if(d >= 0.0)
642 –        SET_NTH_BIT(test[2],1);
643 –    }
644 –    else
645 –      if(sides[2][1] != GT_INTERIOR)
646 –        SET_NTH_BIT(test[2],1);
647 –    /* t=2 */
648 –    if(sides[2][2] == GT_INVALID)
649 –    {
650 –      if(!NTH_BIT(nset,2))
651 –        VCROSS(n[2],t2,t0);
652 –      /* Test the point for sidedness */
653 –      d  = DOT(n[2],p2);
654 –      if(d >= 0.0)
655 –        SET_NTH_BIT(test[2],2);
656 –    }
657 –    else
658 –      if(sides[2][2] != GT_INTERIOR)
659 –        SET_NTH_BIT(test[2],2);
660 – }
661 –
662 –
663 – int
664 – stri_intersect(a1,a2,a3,b1,b2,b3)
665 – FVECT a1,a2,a3,b1,b2,b3;
666 – {
667 –  int which,test,n_set[2];
668 –  int sides[2][3][3],i,j,inext,jnext;
669 –  int tests[2][3];
670 –  FVECT n[2][3],p,avg[2],t1,t2,t3;
671 –
672 – #ifdef DEBUG
673 –    tri_normal(b1,b2,b3,p,FALSE);
674 –    if(DOT(p,b1) > 0)
675 –      {
676 –       VCOPY(t3,b1);
677 –       VCOPY(t2,b2);
678 –       VCOPY(t1,b3);
679 –      }
680 –    else
681 – #endif
682 –      {
683 –       VCOPY(t1,b1);
684 –       VCOPY(t2,b2);
685 –       VCOPY(t3,b3);
686 –      }
687 –
688 –  /* Test the vertices of triangle a against the pyramid formed by triangle
689 –     b and the origin. If any vertex of a is interior to triangle b, or
690 –     if all 3 vertices of a are ON the edges of b,return TRUE. Remember
691 –     the results of the edge normal and sidedness tests for later.
692 –   */
693 – if(vertices_in_stri(a1,a2,a3,t1,t2,t3,&(n_set[0]),n[0],avg[0],sides[1]))
694 –     return(TRUE);
695 –  
696 – if(vertices_in_stri(t1,t2,t3,a1,a2,a3,&(n_set[1]),n[1],avg[1],sides[0]))
697 –     return(TRUE);
698 –
699 –
700 –  set_sidedness_tests(t1,t2,t3,a1,a2,a3,tests[0],sides[0],n_set[1],n[1]);
701 –  if(tests[0][0]&tests[0][1]&tests[0][2])
702 –    return(FALSE);
703 –
704 –  set_sidedness_tests(a1,a2,a3,t1,t2,t3,tests[1],sides[1],n_set[0],n[0]);
705 –  if(tests[1][0]&tests[1][1]&tests[1][2])
706 –    return(FALSE);
707 –
708 –  for(j=0; j < 3;j++)
709 –  {
710 –    jnext = (j+1)%3;
711 –    /* IF edge b doesnt cross any great circles of a, punt */
712 –    if(tests[1][j] & tests[1][jnext])
713 –      continue;
714 –    for(i=0;i<3;i++)
715 –    {
716 –      inext = (i+1)%3;
717 –      /* IF edge a doesnt cross any great circles of b, punt */
718 –      if(tests[0][i] & tests[0][inext])
719 –        continue;
720 –      /* Now find the great circles that cross and test */
721 –      if((NTH_BIT(tests[0][i],j)^(NTH_BIT(tests[0][inext],j)))
722 –          && (NTH_BIT(tests[1][j],i)^NTH_BIT(tests[1][jnext],i)))
723 –      {
724 –        VCROSS(p,n[0][i],n[1][j]);
725 –                    
726 –        /* If zero cp= done */
727 –        if(ZERO_VEC3(p))
728 –          continue;
729 –        /* check above both planes */
730 –        if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0)
731 –          {
732 –            NEGATE_VEC3(p);
733 –            if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0)
734 –              continue;
735 –          }
736 –        return(TRUE);
737 –      }
738 –    }
739 –  }
740 –  return(FALSE);
741 – }
742 –
743 – int
244   ray_intersect_tri(orig,dir,v0,v1,v2,pt)
245   FVECT orig,dir;
246   FVECT v0,v1,v2;
# Line 763 | Line 263 | FVECT pt;
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 773 | 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 */
776 –    /* hv and vv are the horizontal and vertical vectors in the
777 –       view frame-the magnitude is the dimension of the front frustum
778 –       face at z =1
779 –    */
288      VCOPY(nhv,hv);
289      VCOPY(nvv,vv);
290      w2 = normalize(nhv);
# Line 819 | Line 327 | FVECT fnear[4],ffar[4];
327      ffar[3][2] =  width*nhv[2] - height*nvv[2] + far*ndv[2] + vp[2] ;
328   }
329  
822 – int
823 – max_index(v,r)
824 – FVECT v;
825 – double *r;
826 – {
827 –  double p[3];
828 –  int i;
330  
331 <  p[0] = fabs(v[0]);
332 <  p[1] = fabs(v[1]);
333 <  p[2] = fabs(v[2]);
334 <  i = (p[0]>=p[1])?((p[0]>=p[2])?0:2):((p[1]>=p[2])?1:2);  
335 <  if(r)
336 <    *r = p[i];
337 <  return(i);
837 < }
838 <
839 < int
840 < closest_point_in_tri(p0,p1,p2,p,p0id,p1id,p2id)
841 < FVECT p0,p1,p2,p;
842 < int p0id,p1id,p2id;
843 < {
844 <    double d,d1;
845 <    int i;
846 <    
847 <    d =  DIST_SQ(p,p0);
848 <    d1 = DIST_SQ(p,p1);
849 <    if(d < d1)
850 <    {
851 <      d1 = DIST_SQ(p,p2);
852 <      i = (d1 < d)?p2id:p0id;
853 <    }
854 <    else
855 <    {
856 <      d = DIST_SQ(p,p2);
857 <      i = (d < d1)? p2id:p1id;
858 <    }
859 <    return(i);
860 < }
861 <
862 <
863 < int
864 < sedge_intersect(a0,a1,b0,b1)
865 < FVECT a0,a1,b0,b1;
866 < {
867 <    FVECT na,nb,avga,avgb,p;
868 <    double d;
869 <    int sb0,sb1,sa0,sa1;
870 <
871 <    /* First test if edge b straddles great circle of a */
872 <    VCROSS(na,a0,a1);
873 <    d = DOT(na,b0);
874 <    sb0 = ZERO(d)?0:(d<0)? -1: 1;
875 <    d = DOT(na,b1);
876 <    sb1 = ZERO(d)?0:(d<0)? -1: 1;
877 <    /* edge b entirely on one side of great circle a: edges cannot intersect*/
878 <    if(sb0*sb1 > 0)
879 <       return(FALSE);
880 <    /* test if edge a straddles great circle of b */
881 <    VCROSS(nb,b0,b1);
882 <    d = DOT(nb,a0);
883 <    sa0 = ZERO(d)?0:(d<0)? -1: 1;
884 <    d = DOT(nb,a1);
885 <    sa1 = ZERO(d)?0:(d<0)? -1: 1;
886 <    /* edge a entirely on one side of great circle b: edges cannot intersect*/
887 <    if(sa0*sa1 > 0)
888 <       return(FALSE);
889 <
890 <    /* Find one of intersection points of the great circles */
891 <    VCROSS(p,na,nb);
892 <    /* If they lie on same great circle: call an intersection */
893 <    if(ZERO_VEC3(p))
894 <       return(TRUE);
895 <
896 <    VADD(avga,a0,a1);
897 <    VADD(avgb,b0,b1);
898 <    if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
899 <    {
900 <      NEGATE_VEC3(p);
901 <      if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
902 <        return(FALSE);
903 <    }
904 <    if((!sb0 || !sb1) && (!sa0 || !sa1))
905 <      return(FALSE);
906 <    return(TRUE);
907 < }
908 <
909 <
910 < /* Find the normalized barycentric coordinates of p relative to
911 < * triangle v0,v1,v2. Return result in coord
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 927 | Line 353 | double coord[3];
353  
354  
355  
930 – bary_parent(coord,i)
931 – BCOORD coord[3];
932 – int i;
933 – {
934 –  switch(i) {
935 –  case 0:
936 –    /* update bary for child */
937 –    coord[0] = (coord[0] >> 1) + MAXBCOORD4;
938 –    coord[1] >>= 1;
939 –    coord[2] >>= 1;
940 –    break;
941 –  case 1:
942 –    coord[0] >>= 1;
943 –    coord[1]  = (coord[1] >> 1) + MAXBCOORD4;
944 –    coord[2] >>= 1;
945 –    break;
946 –    
947 –  case 2:
948 –    coord[0] >>= 1;
949 –    coord[1] >>= 1;
950 –    coord[2] = (coord[2] >> 1) + MAXBCOORD4;
951 –    break;
952 –    
953 –  case 3:
954 –    coord[0] = MAXBCOORD4 - (coord[0] >> 1);
955 –    coord[1] = MAXBCOORD4 - (coord[1] >> 1);
956 –    coord[2] = MAXBCOORD4 - (coord[2] >> 1);
957 –    break;
958 – #ifdef DEBUG
959 –  default:
960 –    eputs("bary_parent():Invalid child\n");
961 –    break;
962 – #endif
963 –  }
964 – }
356  
966 – bary_from_child(coord,child,next)
967 – BCOORD coord[3];
968 – int child,next;
969 – {
970 – #ifdef DEBUG
971 –  if(child <0 || child > 3)
972 –  {
973 –    eputs("bary_from_child():Invalid child\n");
974 –    return;
975 –  }
976 –  if(next <0 || next > 3)
977 –  {
978 –    eputs("bary_from_child():Invalid next\n");
979 –    return;
980 –  }
981 – #endif
982 –  if(next == child)
983 –    return;
357  
985 –  switch(child){
986 –  case 0:
987 –      coord[0] = 0;
988 –      coord[1] = MAXBCOORD2 - coord[1];
989 –      coord[2] = MAXBCOORD2 - coord[2];
990 –      break;
991 –  case 1:
992 –      coord[0] = MAXBCOORD2 - coord[0];
993 –      coord[1] = 0;
994 –      coord[2] = MAXBCOORD2 - coord[2];
995 –      break;
996 –  case 2:
997 –      coord[0] = MAXBCOORD2 - coord[0];
998 –      coord[1] = MAXBCOORD2 - coord[1];
999 –      coord[2] = 0;
1000 –    break;
1001 –  case 3:
1002 –    switch(next){
1003 –    case 0:
1004 –      coord[0] = 0;
1005 –      coord[1] = MAXBCOORD2 - coord[1];
1006 –      coord[2] = MAXBCOORD2 - coord[2];
1007 –      break;
1008 –    case 1:
1009 –      coord[0] = MAXBCOORD2 - coord[0];
1010 –      coord[1] = 0;
1011 –      coord[2] = MAXBCOORD2 - coord[2];
1012 –      break;
1013 –    case 2:
1014 –      coord[0] = MAXBCOORD2 - coord[0];
1015 –      coord[1] = MAXBCOORD2 - coord[1];
1016 –      coord[2] = 0;
1017 –      break;
1018 –    }
1019 –    break;
1020 –  }
1021 – }
358  
1023 – int
1024 – bary_child(coord)
1025 – BCOORD coord[3];
1026 – {
1027 –  int i;
359  
1029 –  if(coord[0] > MAXBCOORD4)
1030 –  {
1031 –      /* update bary for child */
1032 –      coord[0] = (coord[0]<< 1) - MAXBCOORD2;
1033 –      coord[1] <<= 1;
1034 –      coord[2] <<= 1;
1035 –      return(0);
1036 –  }
1037 –  else
1038 –    if(coord[1] > MAXBCOORD4)
1039 –    {
1040 –      coord[0] <<= 1;
1041 –      coord[1] = (coord[1] << 1) - MAXBCOORD2;
1042 –      coord[2] <<= 1;
1043 –      return(1);
1044 –    }
1045 –    else
1046 –      if(coord[2] > MAXBCOORD4)
1047 –      {
1048 –        coord[0] <<= 1;
1049 –        coord[1] <<= 1;
1050 –        coord[2] = (coord[2] << 1) - MAXBCOORD2;
1051 –        return(2);
1052 –      }
1053 –      else
1054 –         {
1055 –           coord[0] = MAXBCOORD2 - (coord[0] << 1);
1056 –           coord[1] = MAXBCOORD2 - (coord[1] << 1);
1057 –           coord[2] = MAXBCOORD2 - (coord[2] << 1);
1058 –           return(3);
1059 –         }
1060 – }
360  
1062 – int
1063 – bary_nth_child(coord,i)
1064 – BCOORD coord[3];
1065 – int i;
1066 – {
361  
362 <  switch(i){
363 <  case 0:
1070 <    /* update bary for child */
1071 <    coord[0] = (coord[0]<< 1) - MAXBCOORD2;
1072 <    coord[1] <<= 1;
1073 <    coord[2] <<= 1;
1074 <    break;
1075 <  case 1:
1076 <    coord[0] <<= 1;
1077 <    coord[1] = (coord[1] << 1) - MAXBCOORD2;
1078 <    coord[2] <<= 1;
1079 <    break;
1080 <  case 2:
1081 <    coord[0] <<= 1;
1082 <    coord[1] <<= 1;
1083 <    coord[2] = (coord[2] << 1) - MAXBCOORD2;
1084 <    break;
1085 <  case 3:
1086 <    coord[0] = MAXBCOORD2 - (coord[0] << 1);
1087 <    coord[1] = MAXBCOORD2 - (coord[1] << 1);
1088 <    coord[2] = MAXBCOORD2 - (coord[2] << 1);
1089 <    break;
1090 <  }
1091 < }
362 >
363 >
364  
365  
366  

Diff Legend

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