1 #ifndef CONE_HXX
  2 #define CONE_HXX
  3 
  4 #include "Primitive.hxx"
  5 
  6 class Cone : public Primitive
  7 {
  8     Vec3f center;
  9     float radius;
 10     float height;
 11     int dim;
 12 
 13     /** 
 14      * bounding box of a cone
 15      */
 16     Box CalcBounds()
 17     {
 18         Box bb;
 19         // downwards
 20         Vec3f sphere = Vec3f(radius);
 21         sphere[dim] = 0;
 22         bb.Extend( center + sphere );
 23         bb.Extend( center - sphere );
 24         // upwards
 25         Vec3f normal = Vec3f(0);
 26         normal[dim] = height;
 27         bb.Extend ( center + normal + sphere );
 28         bb.Extend ( center + normal - sphere );
 29         
 30         return bb;
 31     }
 32 
 33 public:
 34     /**
 35      * axis aligned cone
 36      * the alignment is not given as parameter, thus: only cone in y-direction!
 37      * since I have a translation class which supports rotation
 38      */
 39     Cone(Vec3f center, float radius, float height)
 40     : center(center), radius(radius), height(height)
 41     {
 42         dim = 1;
 43     };
 44     
 45     virtual ~Cone()
 46     {}
 47     
 48     virtual bool Intersect(Ray &ray)
 49     {   
 50         /**
 51          * formula for the cone:
 52          * (x^2 + y^2) = (h/r)^2 * (z - h)^2
 53          * with ray = org + t * dir, we obtain a quadratic equation
 54          */
 55     
 56         // side
 57         Vec3f diff = ray.org - center;
 58         Vec3f dir  = ray.dir;
 59                 
 60         float ratio = radius / height;
 61         
 62         Vec3f help = Vec3f(1.0);
 63         help[dim] = -powf(ratio,2.0);
 64         
 65         float a = Dot(dir, dir * help );
 66         
 67         float b = Dot( diff, dir * help ) + height * powf(ratio,2.0) * dir[dim];
 68         b *= 2.0;       
 69         
 70         float c = Dot( diff, diff * help );
 71         c += 2.0 * height * powf(ratio,2.0) * diff[dim];
 72         c -= powf(ratio,2.0) * powf(height,2.0);
 73         
 74         /**
 75          * use 'abc'-formula for finding root t_1,2 = (-b +/- sqrt(b^2-4ac))/(2a)
 76          */
 77         float inRoot = b*b - 4*a*c;
 78         if (inRoot < 0) 
 79             return false;
 80         float root = sqrt(inRoot);
 81         
 82         float dist = (-b - root)/(2*a);
 83         if (dist > ray.t)
 84             return false;
 85         
 86         if (dist < Epsilon) 
 87         {
 88             dist = (-b + root)/(2*a);
 89             if (dist < Epsilon || dist > ray.t)
 90                 return false;
 91         }
 92 
 93         /** 
 94          * compute where the hitpoint is lying
 95          * side of cone is hit between 0 and height in dim direction
 96          */
 97         Vec3f hitpoint = ray.org + dist * ray.dir;
 98         float t1 = (hitpoint - center)[dim];
 99         
100         if ( 0 < t1 && t1 < height )
101         {
102             ray.t = dist;
103             ray.hit = this;
104             return true;
105         }
106 
107         /**
108          * and the base is hit at y = 0
109          * thus: compute the plane going through 'center' with dim.th
110          * unit vector as normal vector, to close the cone
111          */
112         Vec3f normal = Vec3f(0);
113         normal[dim] = 1.0;
114 
115         Ray clone = ray;
116         
117         Primitive* plane = new InfinitePlane( center, normal );
118             bool hit = plane->Intersect(clone);
119         delete plane;
120         
121         if ( hit )
122         {
123             /**
124              * is the distance of the hitpoint on the plane to the center smaller than radius^2 ?
125              * then part of the base found
126              * otherwise the ray neither hits the cone's side nor its base 
127              */             
128             Vec3f hitpoint = clone.org + clone.t * clone.dir;
129 
130             if ( Length(hitpoint - center) < powf(radius,2) )
131             {       
132                 ray.t = clone.t;
133                 ray.hit = this;
134 
135                 return true;
136             }
137             else 
138                 return false;
139         }
140         
141         return false;
142     };
143     
144     /**
145      * GetNormal: computes the normal vector of a cone
146      * Therefore determine whether the ray hits the object at the base
147      * (then the normal vector is only one of the unit vector, since we have axis aligned cones)
148      * or whether the ray hits the cone at the side (then the position of the hitpoint determines the normal vector)
149      */
150     virtual Vec3f GetNormal(Ray &ray)
151     {
152         Vec3f hit = (ray.org + ray.t * ray.dir) - center;
153         Vec3f normal = Vec3f(0);
154         
155         if ( fabs(hit[dim]) < Epsilon )
156         {
157             normal[dim] = -1.0;
158         }       
159         else
160         {
161             normal[dim] = hit[dim];
162             normal = hit - normal;
163             Normalize(normal);
164         }
165         
166         return normal;
167     };
168 };
169 
170 #endif


syntax highlighted by Code2HTML, v. 0.9.1