source/geometry/extrudedshape.h

   1 #ifndef MSP_GEOMETRY_EXTRUDEDSHAPE_H_
   2 #define MSP_GEOMETRY_EXTRUDEDSHAPE_H_
   3
   4 #include <algorithm>
   5 #include <cmath>
   6 #include "shape.h"
   7
   8 namespace Msp {
   9 namespace Geometry {
  10
  11 /**
  12 A shape embedded in space of dimension higher by one and extruded towards the
  13 highest dimension.  As an example, extruding a circle creates a cylinder.  The
  14 base shape's orientation is not changed.
  15 */
  16 template<typename T, unsigned D>
  17 class ExtrudedShape: public Shape<T, D>
  18 {
  19 private:
  20         Shape<T, D-1> *base;
  21         T length;
  22
  23 public:
  24         ExtrudedShape(const Shape<T, D-1> &, T);
  25         ExtrudedShape(const ExtrudedShape &);
  26         ExtrudedShape &operator=(const ExtrudedShape &);
  27         virtual ~ExtrudedShape();
  28
  29         virtual ExtrudedShape *clone() const;
  30
  31         const Shape<T, D-1> &get_base() const { return *base; }
  32         T get_length() const { return length; }
  33
  34         virtual BoundingBox<T, D> get_axis_aligned_bounding_box() const;
  35         virtual bool contains(const LinAl::Vector<T, D> &) const;
  36         virtual unsigned get_max_ray_intersections() const;
  37         virtual unsigned get_intersections(const Ray<T, D> &, SurfacePoint<T, D> *, unsigned) const;
  38 };
  39
  40 template<typename T, unsigned D>
  41 inline ExtrudedShape<T, D>::ExtrudedShape(const Shape<T, D-1> &b, T l):
  42         length(l)
  43 {
  44         if(l<=T(0))
  45                 throw std::invalid_argument("ExtrudedShape::ExtrudedShape");
  46
  47         base = b.clone();
  48 }
  49
  50 template<typename T, unsigned D>
  51 inline ExtrudedShape<T, D>::ExtrudedShape(const ExtrudedShape<T, D> &other):
  52         base(other.base.clone()),
  53         length(other.length)
  54 { }
  55
  56 template<typename T, unsigned D>
  57 inline ExtrudedShape<T, D> &ExtrudedShape<T, D>::operator=(const ExtrudedShape<T, D> &other)
  58 {
  59         delete base;
  60         base = other.base.clone();
  61         length = other.length;
  62 }
  63
  64 template<typename T, unsigned D>
  65 inline ExtrudedShape<T, D>::~ExtrudedShape()
  66 {
  67         delete base;
  68 }
  69
  70 template<typename T, unsigned D>
  71 inline ExtrudedShape<T, D> *ExtrudedShape<T, D>::clone() const
  72 {
  73         return new ExtrudedShape<T, D>(*base, length);
  74 }
  75
  76 template<typename T, unsigned D>
  77 inline BoundingBox<T, D> ExtrudedShape<T, D>::get_axis_aligned_bounding_box() const
  78 {
  79         BoundingBox<T, D-1> base_bbox = base->get_axis_aligned_bounding_box();
  80         T half_length = length/T(2);
  81         return BoundingBox<T, D>(compose(base_bbox.get_minimum_point(), -half_length),
  82                 compose(base_bbox.get_maximum_point(), half_length));
  83 }
  84
  85 template<typename T, unsigned D>
  86 inline bool ExtrudedShape<T, D>::contains(const LinAl::Vector<T, D> &point) const
  87 {
  88         using std::abs;
  89
  90         if(abs(point[D-1])>length/T(2))
  91                 return false;
  92
  93         return base->contains(point.template slice<D-1>(0));
  94 }
  95
  96 template<typename T, unsigned D>
  97 inline unsigned ExtrudedShape<T, D>::get_max_ray_intersections() const
  98 {
  99         return std::max(base->get_max_ray_intersections(), 2U);
 100 }
 101
 102 template<typename T, unsigned D>
 103 inline unsigned ExtrudedShape<T, D>::get_intersections(const Ray<T, D> &ray, SurfacePoint<T, D> *points, unsigned size) const
 104 {
 105         using std::swap;
 106
 107         unsigned n = 0;
 108         T half_length = length/T(2);
 109         const LinAl::Vector<T, D> &ray_start = ray.get_start();
 110         const LinAl::Vector<T, D> &ray_direction = ray.get_direction();
 111         LinAl::Vector<T, D-1> base_dir = ray_direction.template slice<D-1>(0);
 112
 113         /* If the ray does not degenerate to a point in the base space, it could
 114         intersect the base shape. */
 115         if(inner_product(base_dir, base_dir)!=T(0))
 116         {
 117                 T offset = T();
 118                 T limit = T();
 119                 if(ray_direction[D-1]!=T(0))
 120                 {
 121                         offset = (half_length-ray_start[D-1])/ray_direction[D-1];
 122                         limit = (-half_length-ray_start[D-1])/ray_direction[D-1];
 123                         if(offset>limit)
 124                                 swap(offset, limit);
 125                         if(offset<T(0))
 126                                 offset = T(0);
 127                 }
 128
 129                 if(limit>=offset)
 130                 {
 131                         T distortion = base_dir.norm();
 132                         Ray<T, D-1> base_ray((ray_start+ray_direction*offset).template slice<D-1>(0),
 133                                 base_dir, (limit-offset)*distortion);
 134
 135                         SurfacePoint<T, D-1> *base_points = 0;
 136                         if(points)
 137                                 /* Shamelessly reuse the provided storage.  Align to the end of the
 138                                 array so processing can start from the first (nearest) point. */
 139                                 base_points = reinterpret_cast<SurfacePoint<T, D-1> *>(points+size)-size;
 140
 141                         unsigned count = base->get_intersections(base_ray, base_points, size);
 142                         for(unsigned i=0; (n<size && i<count); ++i)
 143                         {
 144                                 if(points)
 145                                 {
 146                                         T x = offset+base_points[i].distance/distortion;
 147                                         points[n].position = ray_start+ray_direction*x;
 148                                         points[n].normal = compose(base_points[i].normal, T(0));
 149                                         points[n].distance = x;
 150                                         points[n].entry = base_points[i].entry;
 151                                 }
 152
 153                                 ++n;
 154                         }
 155                 }
 156         }
 157
 158         /* If the ray is not parallel to the base space, it may pass through the
 159         caps. */
 160         if(n<size && ray_direction[D-1]!=T(0))
 161         {
 162                 for(int i=-1; (n<size && i<=1); i+=2)
 163                 {
 164                         T x = (half_length*i-ray_start[D-1])/ray_direction[D-1];
 165                         if(!ray.check_limits(x))
 166                                 continue;
 167
 168                         LinAl::Vector<T, D> p = ray_start+ray_direction*x;
 169                         if(base->contains(p.template slice<D-1>(0)))
 170                         {
 171                                 if(points)
 172                                 {
 173                                         points[n].position = p;
 174                                         points[n].normal = LinAl::Vector<T, D>();
 175                                         points[n].normal[D-1] = i;
 176                                         points[n].distance = x;
 177                                         points[n].entry = (T(i)*ray_direction[D-1]<T(0));
 178                                 }
 179
 180                                 ++n;
 181                         }
 182                 }
 183
 184                 sort_points(points, n);
 185         }
 186
 187         return n;
 188 }
 189
 190 } // namespace Geometry
 191 } // namespace Msp
 192
 193 #endif