ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/sphere.c
Revision: 2.7
Committed: Tue Mar 30 16:13:01 2004 UTC (20 years ago) by schorsch
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R7P2, rad3R7P1, rad4R1, rad4R0, rad3R6, rad3R6P1, rad3R8, rad3R9
Changes since 2.6: +7 -5 lines
Log Message:
Continued ANSIfication. There are only bits and pieces left now.

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: sphere.c,v 2.6 2003/06/26 00:58:10 schorsch Exp $";
3 #endif
4 /*
5 * sphere.c - compute ray intersection with spheres.
6 */
7
8 #include "copyright.h"
9
10 #include "ray.h"
11 #include "otypes.h"
12 #include "rtotypes.h"
13
14
15 extern int
16 o_sphere( /* compute intersection with sphere */
17 OBJREC *so,
18 register RAY *r
19 )
20 {
21 double a, b, c; /* coefficients for quadratic equation */
22 double root[2]; /* quadratic roots */
23 int nroots;
24 double t;
25 register RREAL *ap;
26 register int i;
27
28 if (so->oargs.nfargs != 4)
29 objerror(so, USER, "bad # arguments");
30 ap = so->oargs.farg;
31 if (ap[3] < -FTINY) {
32 objerror(so, WARNING, "negative radius");
33 so->otype = so->otype == OBJ_SPHERE ?
34 OBJ_BUBBLE : OBJ_SPHERE;
35 ap[3] = -ap[3];
36 } else if (ap[3] <= FTINY)
37 objerror(so, USER, "zero radius");
38
39 /*
40 * We compute the intersection by substituting into
41 * the surface equation for the sphere. The resulting
42 * quadratic equation in t is then solved for the
43 * smallest positive root, which is our point of
44 * intersection.
45 * Since the ray is normalized, a should always be
46 * one. We compute it here to prevent instability in the
47 * intersection calculation.
48 */
49 /* compute quadratic coefficients */
50 a = b = c = 0.0;
51 for (i = 0; i < 3; i++) {
52 a += r->rdir[i]*r->rdir[i];
53 t = r->rorg[i] - ap[i];
54 b += 2.0*r->rdir[i]*t;
55 c += t*t;
56 }
57 c -= ap[3] * ap[3];
58
59 nroots = quadratic(root, a, b, c); /* solve quadratic */
60
61 for (i = 0; i < nroots; i++) /* get smallest positive */
62 if ((t = root[i]) > FTINY)
63 break;
64 if (i >= nroots)
65 return(0); /* no positive root */
66
67 if (t >= r->rot)
68 return(0); /* other is closer */
69
70 r->ro = so;
71 r->rot = t;
72 /* compute normal */
73 a = ap[3];
74 if (so->otype == OBJ_BUBBLE)
75 a = -a; /* reverse */
76 for (i = 0; i < 3; i++) {
77 r->rop[i] = r->rorg[i] + r->rdir[i]*t;
78 r->ron[i] = (r->rop[i] - ap[i]) / a;
79 }
80 r->rod = -DOT(r->rdir, r->ron);
81 r->rox = NULL;
82 r->pert[0] = r->pert[1] = r->pert[2] = 0.0;
83 r->uv[0] = r->uv[1] = 0.0;
84
85 return(1); /* hit */
86 }