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.8 by gwlarson, Wed Nov 11 12:05:38 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   }
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)smMesh->samples->max_samp+4);
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 <  
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]);
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 > 0);
119 – }
120 –
121 –
122 –
123 – point_on_sphere(ps,p,c)
124 – FVECT ps,p,c;
125 – {
126 –    VSUB(ps,p,c);    
127 –    normalize(ps);
128 – }
129 –
130 –
131 – /* returns TRUE if ray from origin in direction v intersects plane defined
132 –  * by normal plane_n, and plane_d. If plane is not parallel- returns
133 –  * intersection point if r != NULL. If tptr!= NULL returns value of
134 –  * t, if parallel, returns t=FHUGE
135 –  */
136 – int
137 – intersect_vector_plane(v,peq,tptr,r)
138 –   FVECT v;
139 –   FPEQ peq;
140 –   double *tptr;
141 –   FVECT r;
142 – {
143 –  double t,d;
144 –  int hit;
145 –    /*
146 –      Plane is Ax + By + Cz +D = 0:
147 –      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
148 –    */
149 –
150 –    /* line is  l = p1 + (p2-p1)t, p1=origin */
151 –
152 –    /* Solve for t: */
153 –  d = -(DOT(FP_N(peq),v));
154 –  if(ZERO(d))
155 –  {
156 –      t = FHUGE;
157 –      hit = 0;
158 –  }
159 –  else
160 –  {
161 –      t =  FP_D(peq)/d;
162 –      if(t < 0 )
163 –         hit = 0;
164 –      else
165 –         hit = 1;
166 –      if(r)
167 –         {
168 –             r[0] = v[0]*t;
169 –             r[1] = v[1]*t;
170 –             r[2] = v[2]*t;
171 –         }
172 –  }
173 –    if(tptr)
174 –       *tptr = t;
175 –  return(hit);
176 – }
177 –
178 – int
155   intersect_ray_plane(orig,dir,peq,pd,r)
156     FVECT orig,dir;
157     FPEQ peq;
# Line 184 | Line 160 | intersect_ray_plane(orig,dir,peq,pd,r)
160   {
161    double t,d;
162    int hit;
163 <    /*
188 <      Plane is Ax + By + Cz +D = 0:
189 <      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
190 <    */
191 <     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
192 <         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
193 <       line is  l = p1 + (p2-p1)t
194 <     */
195 <    /* Solve for t: */
163 >
164    d = DOT(FP_N(peq),dir);
165    if(ZERO(d))
166       return(0);
# Line 202 | Line 170 | intersect_ray_plane(orig,dir,peq,pd,r)
170         hit = 0;
171      else
172         hit = 1;
205 –
173    if(r)
174       VSUM(r,orig,dir,t);
175  
# Line 211 | 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
227 <    */
228 <     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
229 <         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
230 <       line is  l = p1 + (p2-p1)t
231 <     */
232 <    /* Solve for t: */
233 <    d= DOT(n,dir);
234 <    if(ZERO(d))
235 <       return(0);
236 <    t =  -(DOT(n,orig))/d;
237 <    if(t < 0)
238 <       hit = 0;
239 <    else
240 <       hit = 1;
241 <
242 <  if(r)
243 <     VSUM(r,orig,dir,t);
244 <
245 <    if(pd)
246 <       *pd = t;
247 <  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
252 – intersect_edge_plane(e0,e1,peq,pd,r)
253 –   FVECT e0,e1;
254 –   FPEQ peq;
255 –   double *pd;
256 –   FVECT r;
257 – {
258 –  double t;
259 –  int hit;
260 –  FVECT d;
261 –  /*
262 –      Plane is Ax + By + Cz +D = 0:
263 –      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
264 –    */
265 –     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
266 –         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
267 –       line is  l = p1 + (p2-p1)t
268 –     */
269 –    /* Solve for t: */
270 –  VSUB(d,e1,e0);
271 –  t =  -(DOT(FP_N(peq),e0) + FP_D(peq))/(DOT(FP_N(peq),d));
272 –    if(t < 0)
273 –       hit = 0;
274 –    else
275 –       hit = 1;
276 –
277 –  VSUM(r,e0,d,t);
278 –
279 –    if(pd)
280 –       *pd = t;
281 –  return(hit);
282 – }
283 –
284 –
285 – int
286 – point_in_cone(p,p0,p1,p2)
287 – FVECT p;
288 – FVECT p0,p1,p2;
289 – {
290 –    FVECT np,x_axis,y_axis;
291 –    double d1,d2;
292 –    FPEQ peq;
293 –    
294 –    /* Find the equation of the circle defined by the intersection
295 –       of the cone with the plane defined by p1,p2,p3- project p into
296 –       that plane and do an in-circle test in the plane
297 –     */
298 –    
299 –    /* find the equation of the plane defined by p0-p2 */
300 –    tri_plane_equation(p0,p1,p2,&peq,FALSE);
301 –
302 –    /* define a coordinate system on the plane: the x axis is in
303 –       the direction of np2-np1, and the y axis is calculated from
304 –       n cross x-axis
305 –     */
306 –    /* Project p onto the plane */
307 –    /* NOTE: check this: does sideness check?*/
308 –    if(!intersect_vector_plane(p,peq,NULL,np))
309 –        return(FALSE);
310 –
311 –    /* create coordinate system on  plane: p1-p0 defines the x_axis*/
312 –    VSUB(x_axis,p1,p0);
313 –    normalize(x_axis);
314 –    /* The y axis is  */
315 –    VCROSS(y_axis,FP_N(peq),x_axis);
316 –    normalize(y_axis);
317 –
318 –    VSUB(p1,p1,p0);
319 –    VSUB(p2,p2,p0);
320 –    VSUB(np,np,p0);
321 –    
322 –    p1[0] = VLEN(p1);
323 –    p1[1] = 0;
324 –
325 –    d1 = DOT(p2,x_axis);
326 –    d2 = DOT(p2,y_axis);
327 –    p2[0] = d1;
328 –    p2[1] = d2;
329 –    
330 –    d1 = DOT(np,x_axis);
331 –    d2 = DOT(np,y_axis);
332 –    np[0] = d1;
333 –    np[1] = d2;
334 –
335 –    /* perform the in-circle test in the new coordinate system */
336 –    return(point_in_circle_thru_origin(np,p1,p2));
337 – }
338 –
339 – int
340 – point_set_in_stri(v0,v1,v2,p,n,nset,sides)
341 – FVECT v0,v1,v2,p;
342 – FVECT n[3];
343 – int *nset;
344 – int sides[3];
345 –
346 – {
347 –    double d;
348 –    /* Find the normal to the triangle ORIGIN:v0:v1 */
349 –    if(!NTH_BIT(*nset,0))
350 –    {
351 –        VCROSS(n[0],v0,v1);
352 –        SET_NTH_BIT(*nset,0);
353 –    }
354 –    /* Test the point for sidedness */
355 –    d  = DOT(n[0],p);
356 –
357 –    if(d > 0.0)
358 –     {
359 –       sides[0] =  GT_OUT;
360 –       sides[1] = sides[2] = GT_INVALID;
361 –       return(FALSE);
362 –      }
363 –    else
364 –       sides[0] = GT_INTERIOR;
365 –      
366 –    /* Test next edge */
367 –    if(!NTH_BIT(*nset,1))
368 –    {
369 –        VCROSS(n[1],v1,v2);
370 –        SET_NTH_BIT(*nset,1);
371 –    }
372 –    /* Test the point for sidedness */
373 –    d  = DOT(n[1],p);
374 –    if(d > 0.0)
375 –    {
376 –        sides[1] = GT_OUT;
377 –        sides[2] = GT_INVALID;
378 –        return(FALSE);
379 –    }
380 –    else
381 –       sides[1] = GT_INTERIOR;
382 –    /* Test next edge */
383 –    if(!NTH_BIT(*nset,2))
384 –    {
385 –        VCROSS(n[2],v2,v0);
386 –        SET_NTH_BIT(*nset,2);
387 –    }
388 –    /* Test the point for sidedness */
389 –    d  = DOT(n[2],p);
390 –    if(d > 0.0)
391 –    {
392 –      sides[2] = GT_OUT;
393 –      return(FALSE);
394 –    }
395 –    else
396 –       sides[2] = GT_INTERIOR;
397 –    /* Must be interior to the pyramid */
398 –    return(GT_INTERIOR);
399 – }
400 –
401 –
402 –
403 –
404 – int
208   point_in_stri(v0,v1,v2,p)
209   FVECT v0,v1,v2,p;
210   {
# Line 413 | Line 216 | FVECT v0,v1,v2,p;
216      d  = DOT(n,p);
217      if(d > 0.0)
218        return(FALSE);
416 –
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);
423 –
225      /* Test next edge */
226      VCROSS(n,v2,v0);
227      /* Test the point for sidedness */
# Line 428 | 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
435 – vertices_in_stri(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides)
436 – FVECT t0,t1,t2,p0,p1,p2;
437 – int *nset;
438 – FVECT n[3];
439 – FVECT avg;
440 – int pt_sides[3][3];
441 –
442 – {
443 –    int below_plane[3],test;
444 –
445 –    SUM_3VEC3(avg,t0,t1,t2);
446 –    *nset = 0;
447 –    /* Test vertex v[i] against triangle j*/
448 –    /* Check if v[i] lies below plane defined by avg of 3 vectors
449 –       defining triangle
450 –       */
451 –
452 –    /* test point 0 */
453 –    if(DOT(avg,p0) < 0.0)
454 –      below_plane[0] = 1;
455 –    else
456 –      below_plane[0] = 0;
457 –    /* Test if b[i] lies in or on triangle a */
458 –    test = point_set_in_stri(t0,t1,t2,p0,n,nset,pt_sides[0]);
459 –    /* If pts[i] is interior: done */
460 –    if(!below_plane[0])
461 –      {
462 –        if(test == GT_INTERIOR)
463 –          return(TRUE);
464 –      }
465 –    /* Now test point 1*/
466 –
467 –    if(DOT(avg,p1) < 0.0)
468 –      below_plane[1] = 1;
469 –    else
470 –      below_plane[1]=0;
471 –    /* Test if b[i] lies in or on triangle a */
472 –    test = point_set_in_stri(t0,t1,t2,p1,n,nset,pt_sides[1]);
473 –    /* If pts[i] is interior: done */
474 –    if(!below_plane[1])
475 –    {
476 –      if(test == GT_INTERIOR)
477 –        return(TRUE);
478 –    }
479 –    
480 –    /* Now test point 2 */
481 –    if(DOT(avg,p2) < 0.0)
482 –      below_plane[2] = 1;
483 –    else
484 –      below_plane[2] = 0;
485 –        /* Test if b[i] lies in or on triangle a */
486 –    test = point_set_in_stri(t0,t1,t2,p2,n,nset,pt_sides[2]);
487 –
488 –    /* If pts[i] is interior: done */
489 –    if(!below_plane[2])
490 –      {
491 –        if(test == GT_INTERIOR)
492 –          return(TRUE);
493 –      }
494 –
495 –    /* If all three points below separating plane: trivial reject */
496 –    if(below_plane[0] && below_plane[1] && below_plane[2])
497 –       return(FALSE);
498 –    /* Now check vertices in a against triangle b */
499 –    return(FALSE);
500 – }
501 –
502 –
503 – set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,n)
504 –   FVECT t0,t1,t2,p0,p1,p2;
505 –   int test[3];
506 –   int sides[3][3];
507 –   int nset;
508 –   FVECT n[3];
509 – {
510 –    int t;
511 –    double d;
512 –
513 –    
514 –    /* p=0 */
515 –    test[0] = 0;
516 –    if(sides[0][0] == GT_INVALID)
517 –    {
518 –      if(!NTH_BIT(nset,0))
519 –        VCROSS(n[0],t0,t1);
520 –      /* Test the point for sidedness */
521 –      d  = DOT(n[0],p0);
522 –      if(d >= 0.0)
523 –        SET_NTH_BIT(test[0],0);
524 –    }
525 –    else
526 –      if(sides[0][0] != GT_INTERIOR)
527 –        SET_NTH_BIT(test[0],0);
528 –
529 –    if(sides[0][1] == GT_INVALID)
530 –    {
531 –      if(!NTH_BIT(nset,1))
532 –        VCROSS(n[1],t1,t2);
533 –        /* Test the point for sidedness */
534 –        d  = DOT(n[1],p0);
535 –        if(d >= 0.0)
536 –          SET_NTH_BIT(test[0],1);
537 –    }
538 –    else
539 –      if(sides[0][1] != GT_INTERIOR)
540 –        SET_NTH_BIT(test[0],1);
541 –
542 –    if(sides[0][2] == GT_INVALID)
543 –    {
544 –      if(!NTH_BIT(nset,2))
545 –        VCROSS(n[2],t2,t0);
546 –      /* Test the point for sidedness */
547 –      d  = DOT(n[2],p0);
548 –      if(d >= 0.0)
549 –        SET_NTH_BIT(test[0],2);
550 –    }
551 –    else
552 –      if(sides[0][2] != GT_INTERIOR)
553 –        SET_NTH_BIT(test[0],2);
554 –
555 –    /* p=1 */
556 –    test[1] = 0;
557 –    /* t=0*/
558 –    if(sides[1][0] == GT_INVALID)
559 –    {
560 –      if(!NTH_BIT(nset,0))
561 –        VCROSS(n[0],t0,t1);
562 –      /* Test the point for sidedness */
563 –      d  = DOT(n[0],p1);
564 –      if(d >= 0.0)
565 –        SET_NTH_BIT(test[1],0);
566 –    }
567 –    else
568 –      if(sides[1][0] != GT_INTERIOR)
569 –        SET_NTH_BIT(test[1],0);
570 –    
571 –    /* t=1 */
572 –    if(sides[1][1] == GT_INVALID)
573 –    {
574 –      if(!NTH_BIT(nset,1))
575 –        VCROSS(n[1],t1,t2);
576 –      /* Test the point for sidedness */
577 –      d  = DOT(n[1],p1);
578 –      if(d >= 0.0)
579 –        SET_NTH_BIT(test[1],1);
580 –    }
581 –    else
582 –      if(sides[1][1] != GT_INTERIOR)
583 –        SET_NTH_BIT(test[1],1);
584 –      
585 –    /* t=2 */
586 –    if(sides[1][2] == GT_INVALID)
587 –    {
588 –      if(!NTH_BIT(nset,2))
589 –        VCROSS(n[2],t2,t0);
590 –      /* Test the point for sidedness */
591 –      d  = DOT(n[2],p1);
592 –      if(d >= 0.0)
593 –        SET_NTH_BIT(test[1],2);
594 –    }
595 –    else
596 –      if(sides[1][2] != GT_INTERIOR)
597 –        SET_NTH_BIT(test[1],2);
598 –
599 –    /* p=2 */
600 –    test[2] = 0;
601 –    /* t = 0 */
602 –    if(sides[2][0] == GT_INVALID)
603 –    {
604 –      if(!NTH_BIT(nset,0))
605 –        VCROSS(n[0],t0,t1);
606 –      /* Test the point for sidedness */
607 –      d  = DOT(n[0],p2);
608 –      if(d >= 0.0)
609 –        SET_NTH_BIT(test[2],0);
610 –    }
611 –    else
612 –      if(sides[2][0] != GT_INTERIOR)
613 –        SET_NTH_BIT(test[2],0);
614 –    /* t=1 */
615 –    if(sides[2][1] == GT_INVALID)
616 –    {
617 –      if(!NTH_BIT(nset,1))
618 –        VCROSS(n[1],t1,t2);
619 –      /* Test the point for sidedness */
620 –      d  = DOT(n[1],p2);
621 –      if(d >= 0.0)
622 –        SET_NTH_BIT(test[2],1);
623 –    }
624 –    else
625 –      if(sides[2][1] != GT_INTERIOR)
626 –        SET_NTH_BIT(test[2],1);
627 –    /* t=2 */
628 –    if(sides[2][2] == GT_INVALID)
629 –    {
630 –      if(!NTH_BIT(nset,2))
631 –        VCROSS(n[2],t2,t0);
632 –      /* Test the point for sidedness */
633 –      d  = DOT(n[2],p2);
634 –      if(d >= 0.0)
635 –        SET_NTH_BIT(test[2],2);
636 –    }
637 –    else
638 –      if(sides[2][2] != GT_INTERIOR)
639 –        SET_NTH_BIT(test[2],2);
640 – }
641 –
642 –
643 – int
644 – stri_intersect(a1,a2,a3,b1,b2,b3)
645 – FVECT a1,a2,a3,b1,b2,b3;
646 – {
647 –  int which,test,n_set[2];
648 –  int sides[2][3][3],i,j,inext,jnext;
649 –  int tests[2][3];
650 –  FVECT n[2][3],p,avg[2];
651 –
652 –  /* Test the vertices of triangle a against the pyramid formed by triangle
653 –     b and the origin. If any vertex of a is interior to triangle b, or
654 –     if all 3 vertices of a are ON the edges of b,return TRUE. Remember
655 –     the results of the edge normal and sidedness tests for later.
656 –   */
657 – if(vertices_in_stri(a1,a2,a3,b1,b2,b3,&(n_set[0]),n[0],avg[0],sides[1]))
658 –     return(TRUE);
659 –  
660 – if(vertices_in_stri(b1,b2,b3,a1,a2,a3,&(n_set[1]),n[1],avg[1],sides[0]))
661 –     return(TRUE);
662 –
663 –
664 –  set_sidedness_tests(b1,b2,b3,a1,a2,a3,tests[0],sides[0],n_set[1],n[1]);
665 –  if(tests[0][0]&tests[0][1]&tests[0][2])
666 –    return(FALSE);
667 –
668 –  set_sidedness_tests(a1,a2,a3,b1,b2,b3,tests[1],sides[1],n_set[0],n[0]);
669 –  if(tests[1][0]&tests[1][1]&tests[1][2])
670 –    return(FALSE);
671 –
672 –  for(j=0; j < 3;j++)
673 –  {
674 –    jnext = (j+1)%3;
675 –    /* IF edge b doesnt cross any great circles of a, punt */
676 –    if(tests[1][j] & tests[1][jnext])
677 –      continue;
678 –    for(i=0;i<3;i++)
679 –    {
680 –      inext = (i+1)%3;
681 –      /* IF edge a doesnt cross any great circles of b, punt */
682 –      if(tests[0][i] & tests[0][inext])
683 –        continue;
684 –      /* Now find the great circles that cross and test */
685 –      if((NTH_BIT(tests[0][i],j)^(NTH_BIT(tests[0][inext],j)))
686 –          && (NTH_BIT(tests[1][j],i)^NTH_BIT(tests[1][jnext],i)))
687 –      {
688 –        VCROSS(p,n[0][i],n[1][j]);
689 –                    
690 –        /* If zero cp= done */
691 –        if(ZERO_VEC3(p))
692 –          continue;
693 –        /* check above both planes */
694 –        if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0)
695 –          {
696 –            NEGATE_VEC3(p);
697 –            if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0)
698 –              continue;
699 –          }
700 –        return(TRUE);
701 –      }
702 –    }
703 –  }
704 –  return(FALSE);
705 – }
706 –
707 – int
244   ray_intersect_tri(orig,dir,v0,v1,v2,pt)
245   FVECT orig,dir;
246   FVECT v0,v1,v2;
# Line 727 | 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 737 | 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 */
740 –    /* hv and vv are the horizontal and vertical vectors in the
741 –       view frame-the magnitude is the dimension of the front frustum
742 –       face at z =1
743 –    */
288      VCOPY(nhv,hv);
289      VCOPY(nvv,vv);
290      w2 = normalize(nhv);
# Line 783 | Line 327 | FVECT fnear[4],ffar[4];
327      ffar[3][2] =  width*nhv[2] - height*nvv[2] + far*ndv[2] + vp[2] ;
328   }
329  
786 – int
787 – max_index(v,r)
788 – FVECT v;
789 – double *r;
790 – {
791 –  double p[3];
792 –  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);
801 < }
802 <
803 < int
804 < closest_point_in_tri(p0,p1,p2,p,p0id,p1id,p2id)
805 < FVECT p0,p1,p2,p;
806 < int p0id,p1id,p2id;
807 < {
808 <    double d,d1;
809 <    int i;
810 <    
811 <    d =  DIST_SQ(p,p0);
812 <    d1 = DIST_SQ(p,p1);
813 <    if(d < d1)
814 <    {
815 <      d1 = DIST_SQ(p,p2);
816 <      i = (d1 < d)?p2id:p0id;
817 <    }
818 <    else
819 <    {
820 <      d = DIST_SQ(p,p2);
821 <      i = (d < d1)? p2id:p1id;
822 <    }
823 <    return(i);
824 < }
825 <
826 <
827 < int
828 < sedge_intersect(a0,a1,b0,b1)
829 < FVECT a0,a1,b0,b1;
830 < {
831 <    FVECT na,nb,avga,avgb,p;
832 <    double d;
833 <    int sb0,sb1,sa0,sa1;
834 <
835 <    /* First test if edge b straddles great circle of a */
836 <    VCROSS(na,a0,a1);
837 <    d = DOT(na,b0);
838 <    sb0 = ZERO(d)?0:(d<0)? -1: 1;
839 <    d = DOT(na,b1);
840 <    sb1 = ZERO(d)?0:(d<0)? -1: 1;
841 <    /* edge b entirely on one side of great circle a: edges cannot intersect*/
842 <    if(sb0*sb1 > 0)
843 <       return(FALSE);
844 <    /* test if edge a straddles great circle of b */
845 <    VCROSS(nb,b0,b1);
846 <    d = DOT(nb,a0);
847 <    sa0 = ZERO(d)?0:(d<0)? -1: 1;
848 <    d = DOT(nb,a1);
849 <    sa1 = ZERO(d)?0:(d<0)? -1: 1;
850 <    /* edge a entirely on one side of great circle b: edges cannot intersect*/
851 <    if(sa0*sa1 > 0)
852 <       return(FALSE);
853 <
854 <    /* Find one of intersection points of the great circles */
855 <    VCROSS(p,na,nb);
856 <    /* If they lie on same great circle: call an intersection */
857 <    if(ZERO_VEC3(p))
858 <       return(TRUE);
859 <
860 <    VADD(avga,a0,a1);
861 <    VADD(avgb,b0,b1);
862 <    if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
863 <    {
864 <      NEGATE_VEC3(p);
865 <      if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
866 <        return(FALSE);
867 <    }
868 <    if((!sb0 || !sb1) && (!sa0 || !sa1))
869 <      return(FALSE);
870 <    return(TRUE);
871 < }
872 <
873 <
874 <
875 < /* Find the normalized barycentric coordinates of p relative to
876 < * 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 889 | Line 350 | double coord[3];
350  
351   }
352  
892 – bary_ith_child(coord,i)
893 – double coord[3];
894 – int i;
895 – {
896 –
897 –  switch(i){
898 –  case 0:
899 –      /* update bary for child */
900 –      coord[0] = 2.0*coord[0]- 1.0;
901 –      coord[1] *= 2.0;
902 –      coord[2] *= 2.0;
903 –      break;
904 –  case 1:
905 –    coord[0] *= 2.0;
906 –    coord[1] = 2.0*coord[1]- 1.0;
907 –    coord[2] *= 2.0;
908 –    break;
909 –  case 2:
910 –    coord[0] *= 2.0;
911 –    coord[1] *= 2.0;
912 –    coord[2] = 2.0*coord[2]- 1.0;
913 –    break;
914 –  case 3:
915 –    coord[0] = 1.0 - 2.0*coord[0];
916 –    coord[1] = 1.0 - 2.0*coord[1];
917 –    coord[2] = 1.0 - 2.0*coord[2];
918 –    break;
919 – #ifdef DEBUG
920 –  default:
921 –    eputs("bary_ith_child():Invalid child\n");
922 –    break;
923 – #endif
924 –  }
925 – }
353  
354  
928 – int
929 – bary_child(coord)
930 – double coord[3];
931 – {
932 –  int i;
355  
934 –  if(coord[0] > 0.5)
935 –  {
936 –      /* update bary for child */
937 –      coord[0] = 2.0*coord[0]- 1.0;
938 –      coord[1] *= 2.0;
939 –      coord[2] *= 2.0;
940 –      return(0);
941 –  }
942 –  else
943 –    if(coord[1] > 0.5)
944 –    {
945 –      coord[0] *= 2.0;
946 –      coord[1] = 2.0*coord[1]- 1.0;
947 –      coord[2] *= 2.0;
948 –      return(1);
949 –    }
950 –    else
951 –      if(coord[2] > 0.5)
952 –      {
953 –        coord[0] *= 2.0;
954 –        coord[1] *= 2.0;
955 –        coord[2] = 2.0*coord[2]- 1.0;
956 –        return(2);
957 –      }
958 –      else
959 –         {
960 –           coord[0] = 1.0 - 2.0*coord[0];
961 –           coord[1] = 1.0 - 2.0*coord[1];
962 –           coord[2] = 1.0 - 2.0*coord[2];
963 –           return(3);
964 –         }
965 – }
966 –
967 – /* Coord was the ith child of the parent: set the coordinate
968 –   relative to the parent
969 – */
970 – bary_parent(coord,i)
971 – double coord[3];
972 – int i;
973 – {
974 –
975 –  switch(i) {
976 –  case 0:
977 –    /* update bary for child */
978 –    coord[0] = coord[0]*0.5 + 0.5;
979 –    coord[1] *= 0.5;
980 –    coord[2] *= 0.5;
981 –    break;
982 –  case 1:
983 –    coord[0] *= 0.5;
984 –    coord[1]  = coord[1]*0.5 + 0.5;
985 –    coord[2] *= 0.5;
986 –    break;
987 –    
988 –  case 2:
989 –    coord[0] *= 0.5;
990 –    coord[1] *= 0.5;
991 –    coord[2] = coord[2]*0.5 + 0.5;
992 –    break;
993 –    
994 –  case 3:
995 –    coord[0] = 0.5 - 0.5*coord[0];
996 –    coord[1] = 0.5 - 0.5*coord[1];
997 –    coord[2] = 0.5 - 0.5*coord[2];
998 –    break;
999 – #ifdef DEBUG
1000 –  default:
1001 –    eputs("bary_parent():Invalid child\n");
1002 –    break;
1003 – #endif
1004 –  }
1005 – }
1006 –
1007 – bary_from_child(coord,child,next)
1008 – double coord[3];
1009 – int child,next;
1010 – {
1011 – #ifdef DEBUG
1012 –  if(child <0 || child > 3)
1013 –  {
1014 –    eputs("bary_from_child():Invalid child\n");
1015 –    return;
1016 –  }
1017 –  if(next <0 || next > 3)
1018 –  {
1019 –    eputs("bary_from_child():Invalid next\n");
1020 –    return;
1021 –  }
1022 – #endif
1023 –  if(next == child)
1024 –    return;
1025 –
1026 –  switch(child){
1027 –  case 0:
1028 –    switch(next){
1029 –    case 1:
1030 –      coord[0] += 1.0;
1031 –      coord[1] -= 1.0;
1032 –      break;
1033 –    case 2:
1034 –      coord[0] += 1.0;
1035 –      coord[2] -= 1.0;
1036 –      break;
1037 –    case 3:
1038 –      coord[0] *= -1.0;
1039 –      coord[1] = 1 - coord[1];
1040 –      coord[2] = 1 - coord[2];
1041 –      break;
1042 –
1043 –    }
1044 –    break;
1045 –  case 1:
1046 –    switch(next){
1047 –    case 0:
1048 –      coord[0] -= 1.0;
1049 –      coord[1] += 1.0;
1050 –      break;
1051 –    case 2:
1052 –      coord[1] += 1.0;
1053 –      coord[2] -= 1.0;
1054 –      break;
1055 –    case 3:
1056 –      coord[0] = 1 - coord[0];
1057 –      coord[1] *= -1.0;
1058 –      coord[2] = 1 - coord[2];
1059 –      break;
1060 –    }
1061 –    break;
1062 –  case 2:
1063 –    switch(next){
1064 –    case 0:
1065 –      coord[0] -= 1.0;
1066 –      coord[2] += 1.0;
1067 –      break;
1068 –    case 1:
1069 –      coord[1] -= 1.0;
1070 –      coord[2] += 1.0;
1071 –      break;
1072 –    case 3:
1073 –      coord[0] = 1 - coord[0];
1074 –      coord[1] = 1 - coord[1];
1075 –      coord[2] *= -1.0;
1076 –      break;
1077 –    }
1078 –    break;
1079 –  case 3:
1080 –    switch(next){
1081 –    case 0:
1082 –      coord[0] *= -1.0;
1083 –      coord[1] = 1 - coord[1];
1084 –      coord[2] = 1 - coord[2];
1085 –      break;
1086 –    case 1:
1087 –      coord[0] = 1 - coord[0];
1088 –      coord[1] *= -1.0;
1089 –      coord[2] = 1 - coord[2];
1090 –      break;
1091 –    case 2:
1092 –      coord[0] = 1 - coord[0];
1093 –      coord[1] = 1 - coord[1];
1094 –      coord[2] *= -1.0;
1095 –      break;
1096 –    }
1097 –    break;
1098 –  }
1099 – }
1100 –
1101 –
1102 – baryi_parent(coord,i)
1103 – BCOORD coord[3];
1104 – int i;
1105 – {
1106 –
1107 –  switch(i) {
1108 –  case 0:
1109 –    /* update bary for child */
1110 –    coord[0] = (coord[0] >> 1) + MAXBCOORD2;
1111 –    coord[1] >>= 1;
1112 –    coord[2] >>= 1;
1113 –    break;
1114 –  case 1:
1115 –    coord[0] >>= 1;
1116 –    coord[1]  = (coord[1] >> 1) + MAXBCOORD2;
1117 –    coord[2] >>= 1;
1118 –    break;
1119 –    
1120 –  case 2:
1121 –    coord[0] >>= 1;
1122 –    coord[1] >>= 1;
1123 –    coord[2] = (coord[2] >> 1) + MAXBCOORD2;
1124 –    break;
1125 –    
1126 –  case 3:
1127 –    coord[0] = MAXBCOORD2 - (coord[0] >> 1);
1128 –    coord[1] = MAXBCOORD2 - (coord[1] >> 1);
1129 –    coord[2] = MAXBCOORD2 - (coord[2] >> 1);
1130 –    break;
1131 – #ifdef DEBUG
1132 –  default:
1133 –    eputs("baryi_parent():Invalid child\n");
1134 –    break;
1135 – #endif
1136 –  }
1137 – }
1138 –
1139 – baryi_from_child(coord,child,next)
1140 – BCOORD coord[3];
1141 – int child,next;
1142 – {
1143 – #ifdef DEBUG
1144 –  if(child <0 || child > 3)
1145 –  {
1146 –    eputs("baryi_from_child():Invalid child\n");
1147 –    return;
1148 –  }
1149 –  if(next <0 || next > 3)
1150 –  {
1151 –    eputs("baryi_from_child():Invalid next\n");
1152 –    return;
1153 –  }
1154 – #endif
1155 –  if(next == child)
1156 –    return;
1157 –
1158 –  switch(child){
1159 –  case 0:
1160 –      coord[0] = 0;
1161 –      coord[1] = MAXBCOORD - coord[1];
1162 –      coord[2] = MAXBCOORD - coord[2];
1163 –      break;
1164 –  case 1:
1165 –      coord[0] = MAXBCOORD - coord[0];
1166 –      coord[1] = 0;
1167 –      coord[2] = MAXBCOORD - coord[2];
1168 –      break;
1169 –  case 2:
1170 –      coord[0] = MAXBCOORD - coord[0];
1171 –      coord[1] = MAXBCOORD - coord[1];
1172 –      coord[2] = 0;
1173 –    break;
1174 –  case 3:
1175 –    switch(next){
1176 –    case 0:
1177 –      coord[0] = 0;
1178 –      coord[1] = MAXBCOORD - coord[1];
1179 –      coord[2] = MAXBCOORD - coord[2];
1180 –      break;
1181 –    case 1:
1182 –      coord[0] = MAXBCOORD - coord[0];
1183 –      coord[1] = 0;
1184 –      coord[2] = MAXBCOORD - coord[2];
1185 –      break;
1186 –    case 2:
1187 –      coord[0] = MAXBCOORD - coord[0];
1188 –      coord[1] = MAXBCOORD - coord[1];
1189 –      coord[2] = 0;
1190 –      break;
1191 –    }
1192 –    break;
1193 –  }
1194 – }
1195 –
1196 – int
1197 – baryi_child(coord)
1198 – BCOORD coord[3];
1199 – {
1200 –  int i;
1201 –
1202 –  if(coord[0] > MAXBCOORD2)
1203 –  {
1204 –      /* update bary for child */
1205 –      coord[0] = (coord[0]<< 1) - MAXBCOORD;
1206 –      coord[1] <<= 1;
1207 –      coord[2] <<= 1;
1208 –      return(0);
1209 –  }
1210 –  else
1211 –    if(coord[1] > MAXBCOORD2)
1212 –    {
1213 –      coord[0] <<= 1;
1214 –      coord[1] = (coord[1] << 1) - MAXBCOORD;
1215 –      coord[2] <<= 1;
1216 –      return(1);
1217 –    }
1218 –    else
1219 –      if(coord[2] > MAXBCOORD2)
1220 –      {
1221 –        coord[0] <<= 1;
1222 –        coord[1] <<= 1;
1223 –        coord[2] = (coord[2] << 1) - MAXBCOORD;
1224 –        return(2);
1225 –      }
1226 –      else
1227 –         {
1228 –           coord[0] = MAXBCOORD - (coord[0] << 1);
1229 –           coord[1] = MAXBCOORD - (coord[1] << 1);
1230 –           coord[2] = MAXBCOORD - (coord[2] << 1);
1231 –           return(3);
1232 –         }
1233 – }
1234 –
1235 – int
1236 – baryi_nth_child(coord,i)
1237 – BCOORD coord[3];
1238 – int i;
1239 – {
1240 –
1241 –  switch(i){
1242 –  case 0:
1243 –    /* update bary for child */
1244 –    coord[0] = (coord[0]<< 1) - MAXBCOORD;
1245 –    coord[1] <<= 1;
1246 –    coord[2] <<= 1;
1247 –    break;
1248 –  case 1:
1249 –    coord[0] <<= 1;
1250 –    coord[1] = (coord[1] << 1) - MAXBCOORD;
1251 –    coord[2] <<= 1;
1252 –    break;
1253 –  case 2:
1254 –    coord[0] <<= 1;
1255 –    coord[1] <<= 1;
1256 –    coord[2] = (coord[2] << 1) - MAXBCOORD;
1257 –    break;
1258 –  case 3:
1259 –    coord[0] = MAXBCOORD - (coord[0] << 1);
1260 –    coord[1] = MAXBCOORD - (coord[1] << 1);
1261 –    coord[2] = MAXBCOORD - (coord[2] << 1);
1262 –    break;
1263 –  }
1264 – }
1265 –
1266 –
1267 – baryi_children(coord,i,in_tri,rcoord,rin_tri)
1268 – BCOORD coord[3][3];
1269 – int i;
1270 – int in_tri[3];
1271 – BCOORD rcoord[3][3];
1272 – int rin_tri[3];
1273 – {
1274 –  int j;
1275 –
1276 –  for(j=0; j< 3; j++)
1277 –  {
1278 –    if(!in_tri[j])
1279 –    {
1280 –      rin_tri[j]=0;
1281 –      continue;
1282 –    }
1283 –    
1284 –    if(i != 3)
1285 –    {
1286 –      if(coord[j][i] < MAXBCOORD2)
1287 –        {
1288 –          rin_tri[j] = 0;
1289 –          continue;
1290 –        }
1291 –    }
1292 –    else
1293 –      if( !(coord[j][0] <= MAXBCOORD2 && coord[j][1] <= MAXBCOORD2 &&
1294 –            coord[j][2] <= MAXBCOORD2))
1295 –        {
1296 –          rin_tri[j] = 0;
1297 –          continue;
1298 –        }
1299 –      
1300 –    rin_tri[j]=1;
1301 –    VCOPY(rcoord[j],coord[j]);
1302 –    baryi_nth_child(rcoord[j],i);
1303 –  }
1304 –
1305 – }
1306 –
1307 – convert_dtol(b,bi)
1308 – double b[3];
1309 – BCOORD bi[3];
1310 – {
1311 –  int i;
1312 –
1313 –  for(i=0; i < 2;i++)
1314 –  {
1315 –
1316 –    if(b[i] <= 0.0)
1317 –    {
1318 –      bi[i]= 0;
1319 –    }
1320 –    else
1321 –      if(b[i] >= 1.0)
1322 –      {
1323 –        bi[i]= MAXBCOORD;
1324 –      }
1325 –      else
1326 –        bi[i] = (BCOORD)(b[i]*MAXBCOORD);
1327 –  }
1328 –  bi[2] = bi[0] +  bi[1];
1329 –  if(bi[2] > MAXBCOORD)
1330 –  {
1331 –      bi[2] = 0;
1332 –      bi[1] = MAXBCOORD - bi[0];
1333 –  }
1334 –  else
1335 –    bi[2] = MAXBCOORD - bi[2];
1336 –
1337 – }
1338 –
1339 – /* convert barycentric coordinate b in [-eps,1+eps] to [0,MAXLONG],
1340 –   dir unbounded to [-MAXLONG,MAXLONG]
1341 – */
1342 – bary_dtol(b,db,bi,dbi,t,w)
1343 – double b[3],db[3][3];
1344 – BCOORD bi[3];
1345 – BDIR dbi[3][3];
1346 – TINT t[3];
1347 – int w;
1348 – {
1349 –  int i,id,j,k;
1350 –  double d;
1351 –
1352 –  convert_dtol(b,bi);
1353 –
1354 –  for(j=w; j< 3; j++)
1355 –  {
1356 –    if(t[j] == HUGET)
1357 –    {
1358 –      max_index(db[j],&d);
1359 –      for(i=0; i< 3; i++)
1360 –        dbi[j][i] = (BDIR)(db[j][i]/d*MAXBDIR);
1361 –      break;
1362 –    }
1363 –    else
1364 –    {
1365 –      for(i=0; i< 3; i++)
1366 –          dbi[j][i] = (BDIR)(db[j][i]*MAXBDIR);
1367 –    }
1368 –  }
1369 – }
1370 –
1371 –
1372 – /* convert barycentric coordinate b in [-eps,1+eps] to [0,MAXLONG],
1373 –   dir unbounded to [-MAXLONG,MAXLONG]
1374 – */
1375 – bary_dtol_new(b,db,bi,boi,dbi,t)
1376 – double b[4][3],db[3][3];
1377 – BCOORD bi[3],boi[3][3];
1378 – BDIR dbi[3][3];
1379 – int t[3];
1380 – {
1381 –  int i,id,j,k;
1382 –  double d;
1383 –
1384 –  convert_dtol(b[3],bi);
1385 –
1386 –  for(j=0; j<3;j++)
1387 –  {
1388 –    if(t[j] != 1)
1389 –      continue;
1390 –    convert_dtol(b[j],boi[j]);
1391 –  }
1392 –  for(j=0; j< 3; j++)
1393 –  {
1394 –    k = (j+1)%3;
1395 –    if(t[k]==0)
1396 –      continue;
1397 –    if(t[k] == -1)
1398 –      {
1399 –        max_index(db[j],&d);
1400 –        for(i=0; i< 3; i++)
1401 –          dbi[j][i] = (BDIR)(db[j][i]/d*MAXBDIR);
1402 –        t[k] = 0;
1403 –      }
1404 –    else
1405 –      if(t[j] != 1)
1406 –        for(i=0; i< 3; i++)
1407 –          dbi[j][i] = (BDIR)(db[j][i]*MAXBDIR);
1408 –    else
1409 –      for(i=0; i< 3; i++)
1410 –        dbi[j][i] = boi[k][i] - boi[j][i];
1411 –    
1412 –  }
1413 – }
1414 –
1415 –
1416 – bary_dtolb(b,bi,in_tri)
1417 – double b[3][3];
1418 – BCOORD bi[3][3];
1419 – int in_tri[3];
1420 – {
1421 –  int i,j;
1422 –
1423 –  for(j=0; j<3;j++)
1424 –  {
1425 –    if(!in_tri[j])
1426 –      continue;
1427 –    for(i=0; i < 2;i++)
1428 –    {
1429 –    if(b[j][i] <= 0.0)
1430 –    {
1431 –      bi[j][i]= 0;
1432 –    }
1433 –    else
1434 –      if(b[j][i] >= 1.0)
1435 –      {
1436 –        bi[j][i]= MAXBCOORD;
1437 –      }
1438 –      else
1439 –        bi[j][i] = (BCOORD)(b[j][i]*MAXBCOORD);
1440 –    }
1441 –    bi[j][2] = MAXBCOORD - bi[j][0] - bi[j][1];
1442 –    if(bi[j][2] < 0)
1443 –      {
1444 –        bi[j][2] = 0;
1445 –        bi[j][1] = MAXBCOORD - bi[j][0];
1446 –      }
1447 –  }
1448 – }
356  
357  
358  

Diff Legend

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