ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/o_cone.c
Revision: 2.6
Committed: Mon Jun 28 10:07:17 2004 UTC (19 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R7P2, rad3R7P1, rad4R0, rad3R6, rad3R6P1, rad3R8, rad3R9
Changes since 2.5: +6 -1 lines
Log Message:
Added numerical correction for needle-like cones

File Contents

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