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.12 by gwlarson, Thu Jun 10 15:22:22 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;
14 < {
15 < /* Average three triangle vertices to give centroid: return in c */
16 <  c[0] = (v0[0] + v1[0] + v2[0])/3.0;
17 <  c[1] = (v0[1] + v1[1] + v2[1])/3.0;
18 <  c[2] = (v0[2] + v1[2] + v2[2])/3.0;
19 < }
20 <
21 <
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 < vec3_equal(v1,v2)
32 < FVECT v1,v2;
31 > pt_in_cone(p,a,b,c)
32 > double p[3],a[3],b[3],c[3];
33   {
34 <   return(EQUAL(v1[0],v2[0]) && EQUAL(v1[1],v2[1])&& EQUAL(v1[2],v2[2]));
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   }
30 – #if 0
31 – extern FVECT Norm[500];
32 – extern int Ncnt;
55   #endif
56 + #endif
57  
58 < int
59 < convex_angle(v0,v1,v2)
60 < FVECT v0,v1,v2;
38 < {
39 <    double dp,dp1;
40 <    int test,test1;
41 <    FVECT nv0,nv1,nv2;
42 <    FVECT cp01,cp12,cp;
58 >  VSUB(ab,b,a);
59 >  VSUB(ac,c,a);
60 >  VCROSS(r,ab,ac);
61  
62 <    /* test sign of (v0Xv1)X(v1Xv2). v1 */
63 <    /* un-Simplified: */
64 <    VCOPY(nv0,v0);
65 <    normalize(nv0);
66 <    VCOPY(nv1,v1);
67 <    normalize(nv1);
68 <    VCOPY(nv2,v2);
69 <    normalize(nv2);
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 <    VCROSS(cp01,nv0,nv1);
72 <    VCROSS(cp12,nv1,nv2);
73 <    VCROSS(cp,cp01,cp12);
74 <    normalize(cp);
75 <    dp1 = DOT(cp,nv1);
76 <    if(dp1 <= 1e-8 || dp1 >= (1-1e-8)) /* Test if on other side,or colinear*/
77 <      test1 = FALSE;
78 <    else
79 <      test1 = TRUE;
80 <
81 <    dp =   nv0[2]*nv1[0]*nv2[1] -  nv0[2]*nv1[1]*nv2[0] - nv0[0]*nv1[2]*nv2[1]
64 <         + nv0[0]*nv1[1]*nv2[2] +  nv0[1]*nv1[2]*nv2[0] - nv0[1]*nv1[0]*nv2[2];
65 <    
66 <    if(dp <= 1e-8 || dp1 >= (1-1e-8))
67 <      test = FALSE;
68 <    else
69 <      test = TRUE;
70 <
71 <    if(test != test1)
72 <      fprintf(stderr,"test %f simplified  %f\n",dp1,dp);
73 <    return(test1);
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 >  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  
76 – /* 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 89 | 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]);
92 –  
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]);
96 –  
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]);
100 –
113    if(!norm)
114       return(0);
103 –  
115    mag = normalize(n);
105 –
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;
129     FPEQ *peqptr;
130     int norm;
131   {
132      tri_normal(v0,v1,v2,FP_N(*peqptr),norm);
116 –
133      FP_D(*peqptr) = -(DOT(FP_N(*peqptr),v0));
134   }
135  
136 <
137 < /* returns TRUE if ray from origin in direction v intersects plane defined
138 <  * by normal plane_n, and plane_d. If plane is not parallel- returns
139 <  * intersection point if r != NULL. If tptr!= NULL returns value of
140 <  * t, if parallel, returns t=FHUGE
141 <  */
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
127 – intersect_vector_plane(v,peq,tptr,r)
128 –   FVECT v;
129 –   FPEQ peq;
130 –   double *tptr;
131 –   FVECT r;
132 – {
133 –  double t,d;
134 –  int hit;
135 –    /*
136 –      Plane is Ax + By + Cz +D = 0:
137 –      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
138 –    */
139 –
140 –    /* line is  l = p1 + (p2-p1)t, p1=origin */
141 –
142 –    /* Solve for t: */
143 –  d = -(DOT(FP_N(peq),v));
144 –  if(ZERO(d))
145 –  {
146 –      t = FHUGE;
147 –      hit = 0;
148 –  }
149 –  else
150 –  {
151 –      t =  FP_D(peq)/d;
152 –      if(t < 0 )
153 –         hit = 0;
154 –      else
155 –         hit = 1;
156 –      if(r)
157 –         {
158 –             r[0] = v[0]*t;
159 –             r[1] = v[1]*t;
160 –             r[2] = v[2]*t;
161 –         }
162 –  }
163 –    if(tptr)
164 –       *tptr = t;
165 –  return(hit);
166 – }
167 –
168 – int
155   intersect_ray_plane(orig,dir,peq,pd,r)
156     FVECT orig,dir;
157     FPEQ peq;
# Line 174 | Line 160 | intersect_ray_plane(orig,dir,peq,pd,r)
160   {
161    double t,d;
162    int hit;
163 <    /*
178 <      Plane is Ax + By + Cz +D = 0:
179 <      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
180 <    */
181 <     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
182 <         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
183 <       line is  l = p1 + (p2-p1)t
184 <     */
185 <    /* Solve for t: */
163 >
164    d = DOT(FP_N(peq),dir);
165    if(ZERO(d))
166       return(0);
# Line 192 | Line 170 | intersect_ray_plane(orig,dir,peq,pd,r)
170         hit = 0;
171      else
172         hit = 1;
195 –
173    if(r)
174       VSUM(r,orig,dir,t);
175  
# Line 201 | 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;
210 <   FVECT r;
211 < {
212 <  double t,d;
213 <  int hit;
214 <    /*
215 <      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 <    d= DOT(n,dir);
224 <    if(ZERO(d))
225 <       return(0);
226 <    t =  -(DOT(n,orig))/d;
227 <    if(t < 0)
228 <       hit = 0;
229 <    else
230 <       hit = 1;
231 <
232 <  if(r)
233 <     VSUM(r,orig,dir,t);
234 <
235 <    if(pd)
236 <       *pd = t;
237 <  return(hit);
238 < }
239 <
240 < /* Assumption: know crosses plane:dont need to check for 'on' case */
241 < intersect_edge_coord_plane(v0,v1,w,r)
242 < FVECT v0,v1;
243 < int w;
244 < FVECT r;
245 < {
246 <  FVECT dv;
247 <  int wnext;
248 <  double t;
249 <
250 <  VSUB(dv,v1,v0);
251 <  t = -v0[w]/dv[w];
252 <  r[w] = 0.0;
253 <  wnext = (w+1)%3;
254 <  r[wnext] = v0[wnext] + dv[wnext]*t;
255 <  wnext = (w+2)%3;
256 <  r[wnext] = v0[wnext] + dv[wnext]*t;
257 < }
258 <
259 < int
260 < intersect_edge_plane(e0,e1,peq,pd,r)
261 <   FVECT e0,e1;
262 <   FPEQ peq;
263 <   double *pd;
264 <   FVECT r;
265 < {
266 <  double t;
267 <  int hit;
268 <  FVECT d;
269 <  /*
270 <      Plane is Ax + By + Cz +D = 0:
271 <      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
272 <    */
273 <     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
274 <         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
275 <       line is  l = p1 + (p2-p1)t
276 <     */
277 <    /* Solve for t: */
278 <  VSUB(d,e1,e0);
279 <  t =  -(DOT(FP_N(peq),e0) + FP_D(peq))/(DOT(FP_N(peq),d));
280 <    if(t < 0)
281 <       hit = 0;
282 <    else
283 <       hit = 1;
284 <
285 <  VSUM(r,e0,d,t);
286 <
287 <    if(pd)
288 <       *pd = t;
289 <  return(hit);
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],v0,v1);
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],v1,v2);
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],v2,v0);
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 < set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,n)
357 <   FVECT t0,t1,t2,p0,p1,p2;
358 <   int test[3];
359 <   int sides[3][3];
360 <   int nset;
361 <   FVECT n[3];
362 < {
363 <    int t;
364 <    double d;
365 <
366 <    
367 <    /* p=0 */
368 <    test[0] = 0;
369 <    if(sides[0][0] == GT_INVALID)
370 <    {
371 <      if(!NTH_BIT(nset,0))
372 <        VCROSS(n[0],t0,t1);
373 <      /* Test the point for sidedness */
374 <      d  = DOT(n[0],p0);
375 <      if(d >= 0.0)
376 <        SET_NTH_BIT(test[0],0);
377 <    }
378 <    else
379 <      if(sides[0][0] != GT_INTERIOR)
380 <        SET_NTH_BIT(test[0],0);
381 <
382 <    if(sides[0][1] == GT_INVALID)
383 <    {
384 <      if(!NTH_BIT(nset,1))
385 <        VCROSS(n[1],t1,t2);
386 <        /* Test the point for sidedness */
387 <        d  = DOT(n[1],p0);
388 <        if(d >= 0.0)
389 <          SET_NTH_BIT(test[0],1);
390 <    }
391 <    else
392 <      if(sides[0][1] != GT_INTERIOR)
393 <        SET_NTH_BIT(test[0],1);
394 <
395 <    if(sides[0][2] == GT_INVALID)
396 <    {
397 <      if(!NTH_BIT(nset,2))
398 <        VCROSS(n[2],t2,t0);
399 <      /* Test the point for sidedness */
400 <      d  = DOT(n[2],p0);
401 <      if(d >= 0.0)
402 <        SET_NTH_BIT(test[0],2);
403 <    }
404 <    else
405 <      if(sides[0][2] != GT_INTERIOR)
406 <        SET_NTH_BIT(test[0],2);
407 <
408 <    /* p=1 */
409 <    test[1] = 0;
410 <    /* t=0*/
411 <    if(sides[1][0] == GT_INVALID)
412 <    {
413 <      if(!NTH_BIT(nset,0))
414 <        VCROSS(n[0],t0,t1);
415 <      /* Test the point for sidedness */
416 <      d  = DOT(n[0],p1);
417 <      if(d >= 0.0)
418 <        SET_NTH_BIT(test[1],0);
419 <    }
420 <    else
421 <      if(sides[1][0] != GT_INTERIOR)
422 <        SET_NTH_BIT(test[1],0);
423 <    
424 <    /* t=1 */
425 <    if(sides[1][1] == GT_INVALID)
426 <    {
427 <      if(!NTH_BIT(nset,1))
428 <        VCROSS(n[1],t1,t2);
429 <      /* Test the point for sidedness */
430 <      d  = DOT(n[1],p1);
431 <      if(d >= 0.0)
432 <        SET_NTH_BIT(test[1],1);
433 <    }
434 <    else
435 <      if(sides[1][1] != GT_INTERIOR)
436 <        SET_NTH_BIT(test[1],1);
437 <      
438 <    /* t=2 */
439 <    if(sides[1][2] == GT_INVALID)
440 <    {
441 <      if(!NTH_BIT(nset,2))
442 <        VCROSS(n[2],t2,t0);
443 <      /* Test the point for sidedness */
444 <      d  = DOT(n[2],p1);
445 <      if(d >= 0.0)
446 <        SET_NTH_BIT(test[1],2);
447 <    }
448 <    else
449 <      if(sides[1][2] != GT_INTERIOR)
450 <        SET_NTH_BIT(test[1],2);
451 <
452 <    /* p=2 */
453 <    test[2] = 0;
454 <    /* t = 0 */
455 <    if(sides[2][0] == GT_INVALID)
456 <    {
457 <      if(!NTH_BIT(nset,0))
458 <        VCROSS(n[0],t0,t1);
459 <      /* Test the point for sidedness */
460 <      d  = DOT(n[0],p2);
461 <      if(d >= 0.0)
462 <        SET_NTH_BIT(test[2],0);
463 <    }
464 <    else
465 <      if(sides[2][0] != GT_INTERIOR)
466 <        SET_NTH_BIT(test[2],0);
467 <    /* t=1 */
468 <    if(sides[2][1] == GT_INVALID)
469 <    {
470 <      if(!NTH_BIT(nset,1))
471 <        VCROSS(n[1],t1,t2);
472 <      /* Test the point for sidedness */
473 <      d  = DOT(n[1],p2);
474 <      if(d >= 0.0)
475 <        SET_NTH_BIT(test[2],1);
476 <    }
477 <    else
478 <      if(sides[2][1] != GT_INTERIOR)
479 <        SET_NTH_BIT(test[2],1);
480 <    /* t=2 */
481 <    if(sides[2][2] == GT_INVALID)
482 <    {
483 <      if(!NTH_BIT(nset,2))
484 <        VCROSS(n[2],t2,t0);
485 <      /* Test the point for sidedness */
486 <      d  = DOT(n[2],p2);
487 <      if(d >= 0.0)
488 <        SET_NTH_BIT(test[2],2);
489 <    }
490 <    else
491 <      if(sides[2][2] != GT_INTERIOR)
492 <        SET_NTH_BIT(test[2],2);
493 < }
494 <
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 d;
192 <    VSUB(ps,p,c);    
193 <    d= normalize(ps);
194 <    return(d);
192 >  
193 >  VSUB(ps,p,c);    
194 >  d = normalize(ps);
195 >  return(d);
196   }
197  
198 + /*
199 + * int
200 + * point_in_stri(v0,v1,v2,p) : Return TRUE if p is in pyramid defined by
201 + *                            tri v0,v1,v2 and origin
202 + * FVECT v0,v1,v2,p;     :Triangle vertices(v0,v1,v2) and point in question(p)
203 + *
204 + *   Tests orientation of p relative to each edge (v0v1,v1v2,v2v0), if it is
205 + *   inside of all 3 edges, returns TRUE, else FALSE.
206 + */
207   int
208   point_in_stri(v0,v1,v2,p)
209   FVECT v0,v1,v2,p;
# Line 514 | Line 216 | FVECT v0,v1,v2,p;
216      d  = DOT(n,p);
217      if(d > 0.0)
218        return(FALSE);
517 –
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);
524 –
225      /* Test next edge */
226      VCROSS(n,v2,v0);
227      /* Test the point for sidedness */
# Line 529 | 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 <
235 > /*
236 > * int
237 > * ray_intersect_tri(orig,dir,v0,v1,v2,pt)
238 > *    : test if ray orig-dir intersects triangle v0v1v2, result in pt
239 > * FVECT orig,dir;   : Vectors defining ray origin and direction
240 > * FVECT v0,v1,v2;   : Triangle vertices
241 > * FVECT pt;         : Intersection point (if any)
242 > */
243   int
244   ray_intersect_tri(orig,dir,v0,v1,v2,pt)
245   FVECT orig,dir;
# Line 556 | 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 566 | 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 */
569 –    /* hv and vv are the horizontal and vertical vectors in the
570 –       view frame-the magnitude is the dimension of the front frustum
571 –       face at z =1
572 –    */
288      VCOPY(nhv,hv);
289      VCOPY(nvv,vv);
290      w2 = normalize(nhv);
# Line 612 | Line 327 | FVECT fnear[4],ffar[4];
327      ffar[3][2] =  width*nhv[2] - height*nvv[2] + far*ndv[2] + vp[2] ;
328   }
329  
615 – int
616 – max_index(v,r)
617 – FVECT v;
618 – double *r;
619 – {
620 –  double p[3];
621 –  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);
630 < }
631 <
632 < int
633 < closest_point_in_tri(p0,p1,p2,p,p0id,p1id,p2id)
634 < FVECT p0,p1,p2,p;
635 < int p0id,p1id,p2id;
636 < {
637 <    double d,d1;
638 <    int i;
639 <    
640 <    d =  DIST_SQ(p,p0);
641 <    d1 = DIST_SQ(p,p1);
642 <    if(d < d1)
643 <    {
644 <      d1 = DIST_SQ(p,p2);
645 <      i = (d1 < d)?p2id:p0id;
646 <    }
647 <    else
648 <    {
649 <      d = DIST_SQ(p,p2);
650 <      i = (d < d1)? p2id:p1id;
651 <    }
652 <    return(i);
653 < }
654 <
655 < /* Find the normalized barycentric coordinates of p relative to
656 < * 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 672 | Line 353 | double coord[3];
353  
354  
355  
675 – bary_parent(coord,i)
676 – BCOORD coord[3];
677 – int i;
678 – {
679 –  switch(i) {
680 –  case 0:
681 –    /* update bary for child */
682 –    coord[0] = (coord[0] >> 1) + MAXBCOORD4;
683 –    coord[1] >>= 1;
684 –    coord[2] >>= 1;
685 –    break;
686 –  case 1:
687 –    coord[0] >>= 1;
688 –    coord[1]  = (coord[1] >> 1) + MAXBCOORD4;
689 –    coord[2] >>= 1;
690 –    break;
691 –    
692 –  case 2:
693 –    coord[0] >>= 1;
694 –    coord[1] >>= 1;
695 –    coord[2] = (coord[2] >> 1) + MAXBCOORD4;
696 –    break;
697 –    
698 –  case 3:
699 –    coord[0] = MAXBCOORD4 - (coord[0] >> 1);
700 –    coord[1] = MAXBCOORD4 - (coord[1] >> 1);
701 –    coord[2] = MAXBCOORD4 - (coord[2] >> 1);
702 –    break;
703 – #ifdef DEBUG
704 –  default:
705 –    eputs("bary_parent():Invalid child\n");
706 –    break;
707 – #endif
708 –  }
709 – }
356  
711 – bary_from_child(coord,child,next)
712 – BCOORD coord[3];
713 – int child,next;
714 – {
715 – #ifdef DEBUG
716 –  if(child <0 || child > 3)
717 –  {
718 –    eputs("bary_from_child():Invalid child\n");
719 –    return;
720 –  }
721 –  if(next <0 || next > 3)
722 –  {
723 –    eputs("bary_from_child():Invalid next\n");
724 –    return;
725 –  }
726 – #endif
727 –  if(next == child)
728 –    return;
357  
730 –  switch(child){
731 –  case 0:
732 –      coord[0] = 0;
733 –      coord[1] = MAXBCOORD2 - coord[1];
734 –      coord[2] = MAXBCOORD2 - coord[2];
735 –      break;
736 –  case 1:
737 –      coord[0] = MAXBCOORD2 - coord[0];
738 –      coord[1] = 0;
739 –      coord[2] = MAXBCOORD2 - coord[2];
740 –      break;
741 –  case 2:
742 –      coord[0] = MAXBCOORD2 - coord[0];
743 –      coord[1] = MAXBCOORD2 - coord[1];
744 –      coord[2] = 0;
745 –    break;
746 –  case 3:
747 –    switch(next){
748 –    case 0:
749 –      coord[0] = 0;
750 –      coord[1] = MAXBCOORD2 - coord[1];
751 –      coord[2] = MAXBCOORD2 - coord[2];
752 –      break;
753 –    case 1:
754 –      coord[0] = MAXBCOORD2 - coord[0];
755 –      coord[1] = 0;
756 –      coord[2] = MAXBCOORD2 - coord[2];
757 –      break;
758 –    case 2:
759 –      coord[0] = MAXBCOORD2 - coord[0];
760 –      coord[1] = MAXBCOORD2 - coord[1];
761 –      coord[2] = 0;
762 –      break;
763 –    }
764 –    break;
765 –  }
766 – }
358  
768 – int
769 – bary_child(coord)
770 – BCOORD coord[3];
771 – {
772 –  int i;
359  
774 –  if(coord[0] > MAXBCOORD4)
775 –  {
776 –      /* update bary for child */
777 –      coord[0] = (coord[0]<< 1) - MAXBCOORD2;
778 –      coord[1] <<= 1;
779 –      coord[2] <<= 1;
780 –      return(0);
781 –  }
782 –  else
783 –    if(coord[1] > MAXBCOORD4)
784 –    {
785 –      coord[0] <<= 1;
786 –      coord[1] = (coord[1] << 1) - MAXBCOORD2;
787 –      coord[2] <<= 1;
788 –      return(1);
789 –    }
790 –    else
791 –      if(coord[2] > MAXBCOORD4)
792 –      {
793 –        coord[0] <<= 1;
794 –        coord[1] <<= 1;
795 –        coord[2] = (coord[2] << 1) - MAXBCOORD2;
796 –        return(2);
797 –      }
798 –      else
799 –         {
800 –           coord[0] = MAXBCOORD2 - (coord[0] << 1);
801 –           coord[1] = MAXBCOORD2 - (coord[1] << 1);
802 –           coord[2] = MAXBCOORD2 - (coord[2] << 1);
803 –           return(3);
804 –         }
805 – }
360  
807 – int
808 – bary_nth_child(coord,i)
809 – BCOORD coord[3];
810 – int i;
811 – {
361  
362 <  switch(i){
363 <  case 0:
815 <    /* update bary for child */
816 <    coord[0] = (coord[0]<< 1) - MAXBCOORD2;
817 <    coord[1] <<= 1;
818 <    coord[2] <<= 1;
819 <    break;
820 <  case 1:
821 <    coord[0] <<= 1;
822 <    coord[1] = (coord[1] << 1) - MAXBCOORD2;
823 <    coord[2] <<= 1;
824 <    break;
825 <  case 2:
826 <    coord[0] <<= 1;
827 <    coord[1] <<= 1;
828 <    coord[2] = (coord[2] << 1) - MAXBCOORD2;
829 <    break;
830 <  case 3:
831 <    coord[0] = MAXBCOORD2 - (coord[0] << 1);
832 <    coord[1] = MAXBCOORD2 - (coord[1] << 1);
833 <    coord[2] = MAXBCOORD2 - (coord[2] << 1);
834 <    break;
835 <  }
836 < }
362 >
363 >
364  
365  
366  

Diff Legend

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