src/rt/o_cone.c

#ifndef lint
static const char RCSid[] = "$Id: o_cone.c,v 2.9 2016/04/21 00:40:35 greg Exp $";
#endif
/*
 *  o_cone.c - routine to determine ray intersection with cones.
 */

#include "copyright.h"

#include  "ray.h"
#include  "otypes.h"
#include  "rtotypes.h"
#include  "cone.h"


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

                                                /* get cone structure */
        co = getcone(o, 1);
        if (co == NULL)
                objerror(o, INTERNAL, "unexpected illegal");

        /*
         *     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 (rayreject(o, r, root[0]))
                        return(0);                      /* have better */
                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];
                VSUM(r->rop, r->rorg, r->rdir, 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 + FTINY)
                        break;                  /* too big */
                                                /* check endpoints */
                VSUM(rox, r->rorg, r->rdir, root[rn]);
                VSUB(rdx, rox, CO_P0(co));
                b = DOT(rdx, co->ad); 
                if (b < 0.0)
                        continue;               /* before p0 */
                if (b > co->al)
                        continue;               /* after p1 */
                if (rayreject(o, r, root[rn]))
                        break;                  /* previous hit better */
                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;
                a = DOT(r->ron, r->ron);
                if (a > 1.+FTINY || a < 1.-FTINY) {
                        c = 1./(.5 + .5*a);     /* avoid numerical error */
                        r->ron[0] *= c; r->ron[1] *= c; r->ron[2] *= c;
                }
                r->rod = -DOT(r->rdir, r->ron);
                r->pert[0] = r->pert[1] = r->pert[2] = 0.0;
                r->uv[0] = r->uv[1] = 0.0;
                r->rox = NULL;
                return(1);                      /* good */
        }
        return(0);
}
Revision:	2.10
Committed:	Sun Jan 31 18:08:04 2021 UTC (3 years, 3 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.9:	+7 -5 lines
Log Message:	feat: added consistent handling of coincident surfaces, thanks to Peter A-B
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	2.10	static const char RCSid[] = "$Id: o_cone.c,v 2.9 2016/04/21 00:40:35 greg Exp $";
3	greg	1.1	#endif
4			/*
5			* o_cone.c - routine to determine ray intersection with cones.
6	greg	2.2	*/
7
8	greg	2.3	#include "copyright.h"
9	greg	1.1
10			#include "ray.h"
11			#include "otypes.h"
12	schorsch	2.5	#include "rtotypes.h"
13	greg	1.1	#include "cone.h"
14
15
16	greg	2.8	int
17	schorsch	2.5	o_cone( /* intersect ray with cone */
18			OBJREC *o,
19	greg	2.8	RAY *r
20	schorsch	2.5	)
21	greg	1.1	{
22			FVECT rox, rdx;
23			double a, b, c;
24			double root[2];
25			int nroots, rn;
26	greg	2.8	CONE *co;
27			int i;
28	greg	1.1
29			/* get cone structure */
30			co = getcone(o, 1);
31	greg	2.9	if (co == NULL)
32			objerror(o, INTERNAL, "unexpected illegal");
33	greg	1.1
34			/*
35			* To intersect a ray with a cone, we transform the
36			* ray into the cone's normalized space. This greatly
37			* simplifies the computation.
38			* For a cone or cup, normalization results in the
39			* equation:
40			*
41			* xx + yy - z*z == 0
42			*
43			* For a cylinder or tube, the normalized equation is:
44			*
45			* xx + yy - r*r == 0
46			*
47			* A normalized ring obeys the following set of equations:
48			*
49			* z == 0 &&
50			* xx + yy >= r0*r0 &&
51			* xx + yy <= r1*r1
52			*/
53
54			/* transform ray */
55			multp3(rox, r->rorg, co->tm);
56			multv3(rdx, r->rdir, co->tm);
57			/* compute intersection */
58
59			if (o->otype == OBJ_CONE \|\| o->otype == OBJ_CUP) {
60
61			a = rdx[0]rdx[0] + rdx[1]rdx[1] - rdx[2]*rdx[2];
62			b = 2.0(rdx[0]rox[0] + rdx[1]rox[1] - rdx[2]rox[2]);
63			c = rox[0]rox[0] + rox[1]rox[1] - rox[2]*rox[2];
64
65			} else if (o->otype == OBJ_CYLINDER \|\| o->otype == OBJ_TUBE) {
66
67			a = rdx[0]rdx[0] + rdx[1]rdx[1];
68			b = 2.0(rdx[0]rox[0] + rdx[1]*rox[1]);
69			c = rox[0]rox[0] + rox[1]rox[1] - CO_R0(co)*CO_R0(co);
70
71			} else { /* OBJ_RING */
72
73	greg	2.10	if ((rdx[2] <= FTINY) & (rdx[2] >= -FTINY))
74	greg	1.1	return(0); /* parallel */
75			root[0] = -rox[2]/rdx[2];
76	greg	2.10	if (rayreject(o, r, root[0]))
77			return(0); /* have better */
78	greg	1.1	b = root[0]*rdx[0] + rox[0];
79			c = root[0]*rdx[1] + rox[1];
80			a = bb + cc;
81			if (a < CO_R0(co)CO_R0(co) \|\| a > CO_R1(co)CO_R1(co))
82			return(0); /* outside radii */
83			r->ro = o;
84			r->rot = root[0];
85	greg	2.7	VSUM(r->rop, r->rorg, r->rdir, r->rot);
86	greg	1.1	VCOPY(r->ron, co->ad);
87			r->rod = -rdx[2];
88	greg	1.3	r->rox = NULL;
89	greg	1.1	return(1); /* good */
90			}
91			/* roots for cone, cup, cyl., tube */
92			nroots = quadratic(root, a, b, c);
93
94			for (rn = 0; rn < nroots; rn++) { /* check real roots */
95			if (root[rn] <= FTINY)
96			continue; /* too small */
97	greg	2.10	if (root[rn] > r->rot + FTINY)
98	greg	1.1	break; /* too big */
99			/* check endpoints */
100	greg	2.7	VSUM(rox, r->rorg, r->rdir, root[rn]);
101			VSUB(rdx, rox, CO_P0(co));
102	greg	1.1	b = DOT(rdx, co->ad);
103			if (b < 0.0)
104			continue; /* before p0 */
105			if (b > co->al)
106			continue; /* after p1 */
107	greg	2.10	if (rayreject(o, r, root[rn]))
108			break; /* previous hit better */
109	greg	1.1	r->ro = o;
110			r->rot = root[rn];
111			VCOPY(r->rop, rox);
112			/* get normal */
113			if (o->otype == OBJ_CYLINDER)
114			a = CO_R0(co);
115			else if (o->otype == OBJ_TUBE)
116			a = -CO_R0(co);
117			else { /* OBJ_CONE \|\| OBJ_CUP */
118			c = CO_R1(co) - CO_R0(co);
119			a = CO_R0(co) + b*c/co->al;
120			if (o->otype == OBJ_CUP) {
121			c = -c;
122			a = -a;
123			}
124			}
125			for (i = 0; i < 3; i++)
126			r->ron[i] = (rdx[i] - b*co->ad[i])/a;
127			if (o->otype == OBJ_CONE \|\| o->otype == OBJ_CUP)
128			for (i = 0; i < 3; i++)
129			r->ron[i] = (co->alr->ron[i] - cco->ad[i])
130	greg	2.7	/ co->sl;
131	greg	2.6	a = DOT(r->ron, r->ron);
132			if (a > 1.+FTINY \|\| a < 1.-FTINY) {
133			c = 1./(.5 + .5a); / avoid numerical error */
134			r->ron[0] = c; r->ron[1] = c; r->ron[2] *= c;
135			}
136	greg	1.1	r->rod = -DOT(r->rdir, r->ron);
137	greg	2.4	r->pert[0] = r->pert[1] = r->pert[2] = 0.0;
138			r->uv[0] = r->uv[1] = 0.0;
139	greg	1.3	r->rox = NULL;
140	greg	1.1	return(1); /* good */
141			}
142			return(0);
143			}