1 /*
2  * cylinder.c - This file contains the functions for dealing with cylinders.
3  *
4  *  $Id: cylinder.c,v 1.26 2007/02/04 05:38:57 johns Exp $
5  */
6 
7 #include "machine.h"
8 #include "types.h"
9 #include "macros.h"
10 #include "vector.h"
11 #include "intersect.h"
12 #include "util.h"
13 
14 #define CYLINDER_PRIVATE
15 #include "cylinder.h"
16 
17 static object_methods cylinder_methods = {
18   (void (*)(const void *, void *))(cylinder_intersect),
19   (void (*)(const void *, const void *, const void *, void *))(cylinder_normal),
20   cylinder_bbox,
21   free
22 };
23 
24 static object_methods fcylinder_methods = {
25   (void (*)(const void *, void *))(fcylinder_intersect),
26   (void (*)(const void *, const void *, const void *, void *))(cylinder_normal),
27   fcylinder_bbox,
28   free
29 };
30 
31 
newcylinder(void * tex,vector ctr,vector axis,flt rad)32 object * newcylinder(void * tex, vector ctr, vector axis, flt rad) {
33   cylinder * c;
34 
35   c=(cylinder *) malloc(sizeof(cylinder));
36   memset(c, 0, sizeof(cylinder));
37   c->methods = &cylinder_methods;
38 
39   c->tex=(texture *) tex;
40   c->ctr=ctr;
41   c->axis=axis;
42   c->rad=rad;
43   return (object *) c;
44 }
45 
cylinder_bbox(void * obj,vector * min,vector * max)46 static int cylinder_bbox(void * obj, vector * min, vector * max) {
47   return 0; /* infinite / unbounded object */
48 }
49 
cylinder_intersect(const cylinder * cyl,ray * ry)50 static void cylinder_intersect(const cylinder * cyl, ray * ry) {
51   vector rc, n, D, O;
52   flt t, s, tin, tout, ln, d;
53 
54   rc.x = ry->o.x - cyl->ctr.x;
55   rc.y = ry->o.y - cyl->ctr.y;
56   rc.z = ry->o.z - cyl->ctr.z;
57 
58   VCross(&ry->d, &cyl->axis, &n);
59 
60   ln=sqrt(n.x*n.x + n.y*n.y + n.z*n.z);    /* finish length calculation */
61 
62   if (ln == 0.0) {  /* ray is parallel to the cylinder.. */
63     VDOT(d, rc, cyl->axis);
64     D.x = rc.x - d * cyl->axis.x;
65     D.y = rc.y - d * cyl->axis.y;
66     D.z = rc.z - d * cyl->axis.z;
67     VDOT(d, D, D);
68     d = sqrt(d);
69     tin = -FHUGE;
70     tout = FHUGE;
71     /* if (d <= cyl->rad) then ray is inside cylinder.. else outside */
72   }
73 
74   n.x /= ln;
75   n.y /= ln;
76   n.z /= ln;
77 
78   VDOT(d, rc, n);
79   d = fabs(d);
80 
81   if (d <= cyl->rad) {  /* ray intersects cylinder.. */
82     VCross(&rc, &cyl->axis, &O);
83     VDOT(t, O, n);
84     t = - t / ln;
85     VCross(&n, &cyl->axis, &O);
86 
87     ln = sqrt(O.x*O.x + O.y*O.y + O.z*O.z);
88     O.x /= ln;
89     O.y /= ln;
90     O.z /= ln;
91 
92     VDOT(s, ry->d, O);
93     s = fabs(sqrt(cyl->rad*cyl->rad - d*d) / s);
94     tin = t - s;
95     ry->add_intersection(tin, (object *) cyl, ry);
96     tout = t + s;
97     ry->add_intersection(tout, (object *) cyl, ry);
98   }
99 }
100 
cylinder_normal(const cylinder * cyl,const vector * pnt,const ray * incident,vector * N)101 static void cylinder_normal(const cylinder * cyl, const vector * pnt, const ray * incident, vector * N) {
102   vector a, b;
103   flt t, invlen, invlen2;
104 
105   a.x = pnt->x - cyl->ctr.x;
106   a.y = pnt->y - cyl->ctr.y;
107   a.z = pnt->z - cyl->ctr.z;
108 
109   b=cyl->axis;
110 
111   invlen = 1.0 / sqrt(b.x*b.x + b.y*b.y + b.z*b.z);
112   b.x *= invlen;
113   b.y *= invlen;
114   b.z *= invlen;
115 
116   VDOT(t, a, b);
117 
118   N->x = pnt->x - (b.x * t + cyl->ctr.x);
119   N->y = pnt->y - (b.y * t + cyl->ctr.y);
120   N->z = pnt->z - (b.z * t + cyl->ctr.z);
121 
122   invlen2 = 1.0 / sqrt(N->x*N->x + N->y*N->y + N->z*N->z);
123   N->x *= invlen2;
124   N->y *= invlen2;
125   N->z *= invlen2;
126 
127   /* Flip surface normal to point toward the viewer if necessary */
128   if (VDot(N, &(incident->d)) > 0.0)  {
129     N->x=-N->x;
130     N->y=-N->y;
131     N->z=-N->z;
132   }
133 }
134 
newfcylinder(void * tex,vector ctr,vector axis,flt rad)135 object * newfcylinder(void * tex, vector ctr, vector axis, flt rad) {
136   cylinder * c;
137 
138   c=(cylinder *) malloc(sizeof(cylinder));
139   memset(c, 0, sizeof(cylinder));
140   c->methods = &fcylinder_methods;
141 
142   c->tex=(texture *) tex;
143   c->ctr=ctr;
144   c->axis=axis;
145   c->rad=rad;
146 
147   return (object *) c;
148 }
149 
fcylinder_bbox(void * obj,vector * min,vector * max)150 static int fcylinder_bbox(void * obj, vector * min, vector * max) {
151   cylinder * c = (cylinder *) obj;
152   vector mintmp, maxtmp;
153 
154   mintmp.x = c->ctr.x;
155   mintmp.y = c->ctr.y;
156   mintmp.z = c->ctr.z;
157   maxtmp.x = c->ctr.x + c->axis.x;
158   maxtmp.y = c->ctr.y + c->axis.y;
159   maxtmp.z = c->ctr.z + c->axis.z;
160 
161   min->x = MYMIN(mintmp.x, maxtmp.x);
162   min->y = MYMIN(mintmp.y, maxtmp.y);
163   min->z = MYMIN(mintmp.z, maxtmp.z);
164   min->x -= c->rad;
165   min->y -= c->rad;
166   min->z -= c->rad;
167 
168   max->x = MYMAX(mintmp.x, maxtmp.x);
169   max->y = MYMAX(mintmp.y, maxtmp.y);
170   max->z = MYMAX(mintmp.z, maxtmp.z);
171   max->x += c->rad;
172   max->y += c->rad;
173   max->z += c->rad;
174 
175   return 1;
176 }
177 
178 
fcylinder_intersect(const cylinder * cyl,ray * ry)179 static void fcylinder_intersect(const cylinder * cyl, ray * ry) {
180   vector rc, n, O, hit, tmp2, ctmp4;
181   flt t, s, tin, tout, ln, d, tmp, tmp3;
182 
183   rc.x = ry->o.x - cyl->ctr.x;
184   rc.y = ry->o.y - cyl->ctr.y;
185   rc.z = ry->o.z - cyl->ctr.z;
186 
187   VCross(&ry->d, &cyl->axis, &n);
188 
189   ln=sqrt(n.x*n.x + n.y*n.y + n.z*n.z);    /* finish length calculation */
190 
191   if (ln == 0.0) {  /* ray is parallel to the cylinder.. */
192     return;       /* in this case, we want to miss or go through the "hole" */
193   }
194 
195   n.x /= ln;
196   n.y /= ln;
197   n.z /= ln;
198 
199   VDOT(d, rc, n);
200   d = fabs(d);
201 
202   if (d <= cyl->rad) {  /* ray intersects cylinder.. */
203     VCross(&rc, &cyl->axis, &O);
204     VDOT(t, O, n);
205     t = - t / ln;
206     VCross(&n, &cyl->axis, &O);
207 
208     ln = sqrt(O.x*O.x + O.y*O.y + O.z*O.z);
209     O.x /= ln;
210     O.y /= ln;
211     O.z /= ln;
212 
213     VDOT(s, ry->d, O);
214     s = fabs(sqrt(cyl->rad*cyl->rad - d*d) / s);
215     tin = t - s;
216 
217     RAYPNT(hit, (*ry), tin);
218 
219     ctmp4=cyl->axis;
220     VNorm(&ctmp4);
221 
222     tmp2.x = hit.x - cyl->ctr.x;
223     tmp2.y = hit.y - cyl->ctr.y;
224     tmp2.z = hit.z - cyl->ctr.z;
225 
226     VDOT(tmp,  tmp2, ctmp4);
227     VDOT(tmp3, cyl->axis, cyl->axis);
228 
229     if ((tmp > 0.0) && (tmp < sqrt(tmp3)))
230       ry->add_intersection(tin, (object *) cyl, ry);
231     tout = t + s;
232 
233     RAYPNT(hit, (*ry), tout);
234 
235     tmp2.x = hit.x - cyl->ctr.x;
236     tmp2.y = hit.y - cyl->ctr.y;
237     tmp2.z = hit.z - cyl->ctr.z;
238 
239     VDOT(tmp,  tmp2, ctmp4);
240     VDOT(tmp3, cyl->axis, cyl->axis);
241 
242     if ((tmp > 0.0) && (tmp < sqrt(tmp3)))
243       ry->add_intersection(tout, (object *) cyl, ry);
244   }
245 }
246 
247