src/rt/o_cone.c

#ifndef lint
static const char       RCSid[] = "$Id$";
#endif
/*
 *  o_cone.c - routine to determine ray intersection with cones.
 */

/* ====================================================================
 * The Radiance Software License, Version 1.0
 *
 * Copyright (c) 1990 - 2002 The Regents of the University of California,
 * through Lawrence Berkeley National Laboratory.   All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *         notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *       notice, this list of conditions and the following disclaimer in
 *       the documentation and/or other materials provided with the
 *       distribution.
 *
 * 3. The end-user documentation included with the redistribution,
 *           if any, must include the following acknowledgment:
 *             "This product includes Radiance software
 *                 (http://radsite.lbl.gov/)
 *                 developed by the Lawrence Berkeley National Laboratory
 *               (http://www.lbl.gov/)."
 *       Alternately, this acknowledgment may appear in the software itself,
 *       if and wherever such third-party acknowledgments normally appear.
 *
 * 4. The names "Radiance," "Lawrence Berkeley National Laboratory"
 *       and "The Regents of the University of California" must
 *       not be used to endorse or promote products derived from this
 *       software without prior written permission. For written
 *       permission, please contact [email protected].
 *
 * 5. Products derived from this software may not be called "Radiance",
 *       nor may "Radiance" appear in their name, without prior written
 *       permission of Lawrence Berkeley National Laboratory.
 *
 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED.   IN NO EVENT SHALL Lawrence Berkeley National Laboratory OR
 * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
 * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
 * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 * ====================================================================
 *
 * This software consists of voluntary contributions made by many
 * individuals on behalf of Lawrence Berkeley National Laboratory.   For more
 * information on Lawrence Berkeley National Laboratory, please see
 * <http://www.lbl.gov/>.
 */

#include  "ray.h"

#include  "otypes.h"

#include  "cone.h"


o_cone(o, r)                    /* intersect ray with cone */
OBJREC  *o;
register RAY  *r;
{
        FVECT  rox, rdx;
        double  a, b, c;
        double  root[2];
        int  nroots, rn;
        register CONE  *co;
        register int  i;

                                                /* get cone structure */
        co = getcone(o, 1);

        /*
         *     To intersect a ray with a cone, we transform the
         *  ray into the cone's normalized space.  This greatly
         *  simplifies the computation.
         *     For a cone or cup, normalization results in the
         *  equation:
         *
         *              x*x + y*y - z*z == 0
         *
         *     For a cylinder or tube, the normalized equation is:
         *
         *              x*x + y*y - r*r == 0
         *
         *     A normalized ring obeys the following set of equations:
         *
         *              z == 0                  &&
         *              x*x + y*y >= r0*r0      &&
         *              x*x + y*y <= r1*r1
         */

                                        /* transform ray */
        multp3(rox, r->rorg, co->tm);
        multv3(rdx, r->rdir, co->tm);
                                        /* compute intersection */

        if (o->otype == OBJ_CONE || o->otype == OBJ_CUP) {

                a = rdx[0]*rdx[0] + rdx[1]*rdx[1] - rdx[2]*rdx[2];
                b = 2.0*(rdx[0]*rox[0] + rdx[1]*rox[1] - rdx[2]*rox[2]);
                c = rox[0]*rox[0] + rox[1]*rox[1] - rox[2]*rox[2];

        } else if (o->otype == OBJ_CYLINDER || o->otype == OBJ_TUBE) {

                a = rdx[0]*rdx[0] + rdx[1]*rdx[1];
                b = 2.0*(rdx[0]*rox[0] + rdx[1]*rox[1]);
                c = rox[0]*rox[0] + rox[1]*rox[1] - CO_R0(co)*CO_R0(co);

        } else { /* OBJ_RING */

                if (rdx[2] <= FTINY && rdx[2] >= -FTINY)
                        return(0);                      /* parallel */
                root[0] = -rox[2]/rdx[2];
                if (root[0] <= FTINY || root[0] >= r->rot)
                        return(0);                      /* distance check */
                b = root[0]*rdx[0] + rox[0];
                c = root[0]*rdx[1] + rox[1];
                a = b*b + c*c;
                if (a < CO_R0(co)*CO_R0(co) || a > CO_R1(co)*CO_R1(co))
                        return(0);                      /* outside radii */
                r->ro = o;
                r->rot = root[0];
                for (i = 0; i < 3; i++)
                        r->rop[i] = r->rorg[i] + r->rdir[i]*r->rot;
                VCOPY(r->ron, co->ad);
                r->rod = -rdx[2];
                r->rox = NULL;
                return(1);                              /* good */
        }
                                        /* roots for cone, cup, cyl., tube */
        nroots = quadratic(root, a, b, c);

        for (rn = 0; rn < nroots; rn++) {       /* check real roots */
                if (root[rn] <= FTINY)
                        continue;               /* too small */
                if (root[rn] >= r->rot)
                        break;                  /* too big */
                                                /* check endpoints */
                for (i = 0; i < 3; i++) {
                        rox[i] = r->rorg[i] + root[rn]*r->rdir[i];
                        rdx[i] = rox[i] - CO_P0(co)[i];
                }
                b = DOT(rdx, co->ad); 
                if (b < 0.0)
                        continue;               /* before p0 */
                if (b > co->al)
                        continue;               /* after p1 */
                r->ro = o;
                r->rot = root[rn];
                VCOPY(r->rop, rox);
                                                /* get normal */
                if (o->otype == OBJ_CYLINDER)
                        a = CO_R0(co);
                else if (o->otype == OBJ_TUBE)
                        a = -CO_R0(co);
                else { /* OBJ_CONE || OBJ_CUP */
                        c = CO_R1(co) - CO_R0(co);
                        a = CO_R0(co) + b*c/co->al;
                        if (o->otype == OBJ_CUP) {
                                c = -c;
                                a = -a;
                        }
                }
                for (i = 0; i < 3; i++)
                        r->ron[i] = (rdx[i] - b*co->ad[i])/a;
                if (o->otype == OBJ_CONE || o->otype == OBJ_CUP)
                        for (i = 0; i < 3; i++)
                                r->ron[i] = (co->al*r->ron[i] - c*co->ad[i])
                                                /co->sl;
                r->rod = -DOT(r->rdir, r->ron);
                r->rox = NULL;
                return(1);                      /* good */
        }
        return(0);
}
Revision:	2.2
Committed:	Sat Feb 22 02:07:29 2003 UTC (21 years, 2 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.1:	+57 -5 lines
Log Message:	Changes and check-in for 3.5 release Includes new source files and modifications not recorded for many years See ray/doc/notes/ReleaseNotes for notes between 3.1 and 3.5 release
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	2.2	static const char RCSid[] = "$Id$";
3	greg	1.1	#endif
4			/*
5			* o_cone.c - routine to determine ray intersection with cones.
6	greg	2.2	*/
7
8			/* ====================================================================
9			* The Radiance Software License, Version 1.0
10			*
11			* Copyright (c) 1990 - 2002 The Regents of the University of California,
12			* through Lawrence Berkeley National Laboratory. All rights reserved.
13			*
14			* Redistribution and use in source and binary forms, with or without
15			* modification, are permitted provided that the following conditions
16			* are met:
17			*
18			* 1. Redistributions of source code must retain the above copyright
19			* notice, this list of conditions and the following disclaimer.
20			*
21			* 2. Redistributions in binary form must reproduce the above copyright
22			* notice, this list of conditions and the following disclaimer in
23			* the documentation and/or other materials provided with the
24			* distribution.
25			*
26			* 3. The end-user documentation included with the redistribution,
27			* if any, must include the following acknowledgment:
28			* "This product includes Radiance software
29			* (http://radsite.lbl.gov/)
30			* developed by the Lawrence Berkeley National Laboratory
31			* (http://www.lbl.gov/)."
32			* Alternately, this acknowledgment may appear in the software itself,
33			* if and wherever such third-party acknowledgments normally appear.
34			*
35			* 4. The names "Radiance," "Lawrence Berkeley National Laboratory"
36			* and "The Regents of the University of California" must
37			* not be used to endorse or promote products derived from this
38			* software without prior written permission. For written
39			* permission, please contact [email protected].
40			*
41			* 5. Products derived from this software may not be called "Radiance",
42			* nor may "Radiance" appear in their name, without prior written
43			* permission of Lawrence Berkeley National Laboratory.
44			*
45			* THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
46			* WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
47			* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
48			* DISCLAIMED. IN NO EVENT SHALL Lawrence Berkeley National Laboratory OR
49			* ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
50			* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
51			* LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
52			* USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
53			* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
54			* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
55			* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
56			* SUCH DAMAGE.
57			* ====================================================================
58	greg	1.1	*
59	greg	2.2	* This software consists of voluntary contributions made by many
60			* individuals on behalf of Lawrence Berkeley National Laboratory. For more
61			* information on Lawrence Berkeley National Laboratory, please see
62			* <http://www.lbl.gov/>.
63	greg	1.1	*/
64
65			#include "ray.h"
66
67			#include "otypes.h"
68
69			#include "cone.h"
70
71
72			o_cone(o, r) /* intersect ray with cone */
73			OBJREC *o;
74			register RAY *r;
75			{
76			FVECT rox, rdx;
77			double a, b, c;
78			double root[2];
79			int nroots, rn;
80			register CONE *co;
81			register int i;
82
83			/* get cone structure */
84			co = getcone(o, 1);
85
86			/*
87			* To intersect a ray with a cone, we transform the
88			* ray into the cone's normalized space. This greatly
89			* simplifies the computation.
90			* For a cone or cup, normalization results in the
91			* equation:
92			*
93			* xx + yy - z*z == 0
94			*
95			* For a cylinder or tube, the normalized equation is:
96			*
97			* xx + yy - r*r == 0
98			*
99			* A normalized ring obeys the following set of equations:
100			*
101			* z == 0 &&
102			* xx + yy >= r0*r0 &&
103			* xx + yy <= r1*r1
104			*/
105
106			/* transform ray */
107			multp3(rox, r->rorg, co->tm);
108			multv3(rdx, r->rdir, co->tm);
109			/* compute intersection */
110
111			if (o->otype == OBJ_CONE \|\| o->otype == OBJ_CUP) {
112
113			a = rdx[0]rdx[0] + rdx[1]rdx[1] - rdx[2]*rdx[2];
114			b = 2.0(rdx[0]rox[0] + rdx[1]rox[1] - rdx[2]rox[2]);
115			c = rox[0]rox[0] + rox[1]rox[1] - rox[2]*rox[2];
116
117			} else if (o->otype == OBJ_CYLINDER \|\| o->otype == OBJ_TUBE) {
118
119			a = rdx[0]rdx[0] + rdx[1]rdx[1];
120			b = 2.0(rdx[0]rox[0] + rdx[1]*rox[1]);
121			c = rox[0]rox[0] + rox[1]rox[1] - CO_R0(co)*CO_R0(co);
122
123			} else { /* OBJ_RING */
124
125			if (rdx[2] <= FTINY && rdx[2] >= -FTINY)
126			return(0); /* parallel */
127			root[0] = -rox[2]/rdx[2];
128			if (root[0] <= FTINY \|\| root[0] >= r->rot)
129			return(0); /* distance check */
130			b = root[0]*rdx[0] + rox[0];
131			c = root[0]*rdx[1] + rox[1];
132			a = bb + cc;
133			if (a < CO_R0(co)CO_R0(co) \|\| a > CO_R1(co)CO_R1(co))
134			return(0); /* outside radii */
135			r->ro = o;
136			r->rot = root[0];
137			for (i = 0; i < 3; i++)
138			r->rop[i] = r->rorg[i] + r->rdir[i]*r->rot;
139			VCOPY(r->ron, co->ad);
140			r->rod = -rdx[2];
141	greg	1.3	r->rox = NULL;
142	greg	1.1	return(1); /* good */
143			}
144			/* roots for cone, cup, cyl., tube */
145			nroots = quadratic(root, a, b, c);
146
147			for (rn = 0; rn < nroots; rn++) { /* check real roots */
148			if (root[rn] <= FTINY)
149			continue; /* too small */
150			if (root[rn] >= r->rot)
151			break; /* too big */
152			/* check endpoints */
153			for (i = 0; i < 3; i++) {
154			rox[i] = r->rorg[i] + root[rn]*r->rdir[i];
155			rdx[i] = rox[i] - CO_P0(co)[i];
156			}
157			b = DOT(rdx, co->ad);
158			if (b < 0.0)
159			continue; /* before p0 */
160			if (b > co->al)
161			continue; /* after p1 */
162			r->ro = o;
163			r->rot = root[rn];
164			VCOPY(r->rop, rox);
165			/* get normal */
166			if (o->otype == OBJ_CYLINDER)
167			a = CO_R0(co);
168			else if (o->otype == OBJ_TUBE)
169			a = -CO_R0(co);
170			else { /* OBJ_CONE \|\| OBJ_CUP */
171			c = CO_R1(co) - CO_R0(co);
172			a = CO_R0(co) + b*c/co->al;
173			if (o->otype == OBJ_CUP) {
174			c = -c;
175			a = -a;
176			}
177			}
178			for (i = 0; i < 3; i++)
179			r->ron[i] = (rdx[i] - b*co->ad[i])/a;
180			if (o->otype == OBJ_CONE \|\| o->otype == OBJ_CUP)
181			for (i = 0; i < 3; i++)
182			r->ron[i] = (co->alr->ron[i] - cco->ad[i])
183			/co->sl;
184			r->rod = -DOT(r->rdir, r->ron);
185	greg	1.3	r->rox = NULL;
186	greg	1.1	return(1); /* good */
187			}
188			return(0);
189			}