]> git.tdb.fi Git - libs/math.git/blob - source/geometry/extrudedshape.h
Fix a case in ExtrudedShape with certain ray parameters
[libs/math.git] / 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<=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>(LinAl::Vector<T, D>(base_bbox.get_minimum_point(), -half_length),
82                 LinAl::Vector<T, D>(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(LinAl::Vector<T, D-1>(point));
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);
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.get_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(LinAl::Vector<T, D-1>(ray_start+ray_direction*offset),
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 array
138                                 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 = LinAl::Vector<T, D>(base_points[i].normal, T(0));
149                                         points[n].distance = x;
150                                 }
151
152                                 ++n;
153                         }
154                 }
155         }
156
157         /* If the ray is not parallel to the base space, it may pass through the
158         caps. */
159         if(n<size && ray_direction[D-1])
160         {
161                 for(int i=-1; (n<size && i<=1); i+=2)
162                 {
163                         T x = (half_length*i-ray_start[D-1])/ray_direction[D-1];
164                         if(!ray.check_limits(x))
165                                 continue;
166
167                         LinAl::Vector<T, D> p = ray_start+ray_direction*x;
168                         if(base->contains(LinAl::Vector<T, D-1>(p)))
169                         {
170                                 if(points)
171                                 {
172                                         points[n].position = p;
173                                         points[n].normal = LinAl::Vector<T, D>();
174                                         points[n].normal[D-1] = i;
175                                         points[n].distance = x;
176                                 }
177
178                                 ++n;
179                         }
180                 }
181
182                 sort_points(points, n);
183         }
184
185         return n;
186 }
187
188 } // namespace Geometry
189 } // namespace Msp
190
191 #endif