1 #ifndef CYLINDER_HXX
  2 #define CYLINDER_HXX
  3 
  4 #include "Primitive.hxx"
  5 #include "InfinitePlane.hxx"
  6 
  7 /**
  8  * axis aligned cylinder
  9  */ 
 10 typedef enum align
 11 {
 12     xAxis,
 13     yAxis,
 14     zAxis
 15 };
 16 
 17 class Cylinder : public Primitive
 18 {
 19     Vec3f _down;
 20     Vec3f _up;
 21     Vec3f _center;
 22 
 23     enum  align _dim;
 24     float radius;
 25     float height;
 26 
 27     /**
 28      * compute the box of a cylinder
 29      * therefore take the center of the sphere upwards and do the same as for the sphere
 30      * and then continue for the center of the sphere downwards
 31      */
 32     Box CalcBounds()
 33     {       
 34         Box BoundingBox;
 35 
 36         Vec3f circle = Vec3f( radius );
 37         circle [ _dim ] = 0;
 38                 
 39         BoundingBox.Extend( _down + circle);
 40         BoundingBox.Extend( _down - circle);
 41 
 42         BoundingBox.Extend( _up + circle );
 43         BoundingBox.Extend( _up - circle );
 44 
 45         return BoundingBox;
 46     }
 47 
 48 public:
 49 
 50     Cylinder(Vec3f point, enum align _dim, float radius, float height)
 51     : _dim(_dim), radius(radius), height(height)
 52     {
 53         // unit vector in _dim direction
 54         Vec3f add = Vec3f(0);
 55         add[_dim] = 1;
 56 
 57         // compute center of the cylinder
 58         _center = point;
 59         // center of the sphere downwards
 60         _down = _center - height * add;
 61         // center of the sphere upwards
 62         _up   = _center + height * add;
 63     };
 64     
 65     virtual ~Cylinder()
 66     {}
 67     
 68     virtual bool Intersect(Ray &ray)
 69     {
 70         /** 
 71          * first: infinite long cylinder
 72          * computation as for sphere, just one dimension missing
 73          */
 74                  
 75         /** 
 76          * x^2 + y^2 = radius^2 
 77          */
 78         Vec3f diff = ray.org - _down;
 79         Vec3f dir  = ray.dir;
 80 
 81         diff[_dim] = 0.0;
 82         dir [_dim] = 0.0;
 83         
 84         float a = Dot(dir,dir);
 85         float b = 2 * Dot(dir,diff);
 86         float c = Dot(diff,diff) - radius * radius;
 87         
 88         float inRoot = b*b - 4*a*c;
 89         if (inRoot < 0) 
 90             return false;
 91         float root = sqrt(inRoot);
 92                 
 93         float dist = (-b - root)/(2*a);
 94         if (dist > ray.t)
 95             return false;
 96         
 97         if (dist < Epsilon) 
 98         {
 99             dist = (-b + root)/(2*a);
100             if (dist < Epsilon || dist > ray.t)
101                 return false;
102         }   
103 
104         /**
105          * side
106          * cut off, so that a 2*height cylinder stays
107          * which is open to both sides
108          */
109         Vec3f hit = ray.org + dist * ray.dir;
110         float t1 = (hit - _center)[_dim];
111                 
112         if ( - height < t1 && t1 < height ) 
113         {
114             ray.t = dist;
115             ray.hit = this;
116 
117             return true;
118         }
119 
120         /** 
121          * compute top and bottom surfaces => close the cylinder
122          * check if ray hits planes located at _down and _up within the radius of the cylinder
123          */
124         Ray ray1 = ray;
125         Ray ray2 = ray;
126                                 
127         Vec3f normal = Vec3f(0);
128         normal[_dim] = 1.0;
129                                 
130         Primitive* plane1 = new InfinitePlane( _down, normal );
131         bool hit_plane1 = plane1->Intersect(ray1);
132         delete plane1;
133         
134         Primitive* plane2 = new InfinitePlane( _up,   normal );
135         bool hit_plane2 = plane2->Intersect(ray2);
136         delete plane2;
137         
138         // bottom surface
139         if (hit_plane1)
140         {
141             Vec3f hit1 = (ray1.org + ray1.t * ray1.dir) - _down;
142             hit1[_dim] = 0.0;
143         
144             if ( Length(hit1) < radius * radius )
145             {                                       
146                 ray.t = ray1.t;
147                 ray.hit = this;
148 
149                 return true;
150             }
151         }
152         
153         // top surface
154         if (hit_plane2)
155         {
156             Vec3f hit2 = (ray2.org + ray2.t * ray2.dir) - _up;
157             hit2[_dim] = 0.0;
158         
159             if ( Length(hit2) < radius * radius )
160             {                                       
161                 ray.t = ray2.t;
162                 ray.hit = this;
163 
164                 return true;
165             }
166         }
167                                 
168         return false;       
169     };
170     
171     virtual Vec3f GetNormal(Ray &ray)
172     {
173         Vec3f hit = (ray.org + ray.t * ray.dir) - _center;
174         Vec3f normal = Vec3f(0);
175         
176         // normal vector of bottom surface
177         if ( fabs(hit[_dim] - height) < Epsilon)
178         {
179             normal[_dim] = +1.0;
180         }
181         // normal vector of top surface
182         else if ( fabs(hit[_dim] + height) < Epsilon)                   
183         {
184             normal[_dim] = -1.0;
185         }
186         // normal vector of side
187         else 
188         {
189             normal[_dim] = hit[_dim];
190             normal = hit - normal;
191             Normalize(normal);                  
192         }
193         
194         return normal;
195     };
196 
197 };
198 
199 #endif


syntax highlighted by Code2HTML, v. 0.9.1