src/rt/sphere.c

#ifndef lint
static const char RCSid[] = "$Id: sphere.c,v 2.5 2003/03/12 04:59:05 greg Exp $";
#endif
/*
 *  sphere.c - compute ray intersection with spheres.
 */

#include "copyright.h"

#include  "ray.h"

#include  "otypes.h"


o_sphere(so, r)                 /* compute intersection with sphere */
OBJREC  *so;
register RAY  *r;
{
        double  a, b, c;        /* coefficients for quadratic equation */
        double  root[2];        /* quadratic roots */
        int  nroots;
        double  t;
        register RREAL  *ap;
        register int  i;

        if (so->oargs.nfargs != 4)
                objerror(so, USER, "bad # arguments");
        ap = so->oargs.farg;
        if (ap[3] < -FTINY) {
                objerror(so, WARNING, "negative radius");
                so->otype = so->otype == OBJ_SPHERE ?
                                OBJ_BUBBLE : OBJ_SPHERE;
                ap[3] = -ap[3];
        } else if (ap[3] <= FTINY)
                objerror(so, USER, "zero radius");

        /*
         *      We compute the intersection by substituting into
         *  the surface equation for the sphere.  The resulting
         *  quadratic equation in t is then solved for the
         *  smallest positive root, which is our point of
         *  intersection.
         *      Since the ray is normalized, a should always be
         *  one.  We compute it here to prevent instability in the
         *  intersection calculation.
         */
                                /* compute quadratic coefficients */
        a = b = c = 0.0;
        for (i = 0; i < 3; i++) {
                a += r->rdir[i]*r->rdir[i];
                t = r->rorg[i] - ap[i];
                b += 2.0*r->rdir[i]*t;
                c += t*t;
        }
        c -= ap[3] * ap[3];

        nroots = quadratic(root, a, b, c);      /* solve quadratic */
        
        for (i = 0; i < nroots; i++)            /* get smallest positive */
                if ((t = root[i]) > FTINY)
                        break;
        if (i >= nroots)
                return(0);                      /* no positive root */

        if (t >= r->rot)
                return(0);                      /* other is closer */

        r->ro = so;
        r->rot = t;
                                        /* compute normal */
        a = ap[3];
        if (so->otype == OBJ_BUBBLE)
                a = -a;                 /* reverse */
        for (i = 0; i < 3; i++) {
                r->rop[i] = r->rorg[i] + r->rdir[i]*t;
                r->ron[i] = (r->rop[i] - ap[i]) / a;
        }
        r->rod = -DOT(r->rdir, r->ron);
        r->rox = NULL;
        r->pert[0] = r->pert[1] = r->pert[2] = 0.0;
        r->uv[0] = r->uv[1] = 0.0;

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