Flow123d  master-f44eb46
bounding_box.hh
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file bounding_box.hh
15  * @brief
16  */
17 
18 #ifndef BOX_ELEMENT_HH_
19 #define BOX_ELEMENT_HH_
20 
21 #include <string.h> // for memcpy
22 #include <algorithm> // for max, min
23 
24 
25 #include <new> // for operator new[]
26 #include <ostream> // for operator<<
27 #include <string> // for basic_string
28 #include <vector> // for vector
29 #include <armadillo>
30 #include "mesh/point.hh" // for Space, Space<>...
31 #include "system/exceptions.hh" // for ExcStream, ope...
32 #include "system/global_defs.h" // for msg, rank, ss
33 
34 using namespace std;
35 
36 /**
37  * @brief Bounding box in 3d ambient space.
38  *
39  * Primary intention is usage in BIHTree and various speedups of non-compatible
40  * intersections.
41  *
42  * Copy constructor and assignment are default provided by compiler.
43  * These can be used to set bounds latter on without particular method
44  * to this end:
45  *
46  * @code
47  * BoundingBox box; // non-initialized box
48  * box=BoundingBox( arma::vec3("0 1 2"), arma::vec3("4 5 6") );
49  * @endcode
50  *
51  * Don;t worry about performance, all is inlined.
52  */
53 class BoundingBox {
54 public:
55  TYPEDEF_ERR_INFO( EI_split_point, double);
56  TYPEDEF_ERR_INFO( EI_interval_left, double);
57  TYPEDEF_ERR_INFO( EI_interval_right, double);
58  DECLARE_EXCEPTION(ExcSplitting, << "Split point " << EI_split_point::val
59  << "out of bounds: <" << EI_interval_left::val
60  << ", " << EI_interval_right::val << ">\n");
61 
62  /// Currently we set dimension to 3.
63  static const unsigned int dimension = 3;
64  /// stabilization parameter
65  static const double epsilon;
66  /// Currently we assume
68 
69 
70  /**
71  * Default constructor.
72  * No initialization of vertices. Be very careful using this.
73  * One necessary usage is vector of BoundigBox.
74  */
76  {
77  // set undefined vertices
78  min_vertex_.fill( std::numeric_limits<double>::signaling_NaN() );
79  max_vertex_.fill( std::numeric_limits<double>::signaling_NaN() );
80  }
81 
82  /**
83  * Constructor for point box.
84  */
85  BoundingBox(const Point &min)
86  : min_vertex_(min), max_vertex_(min)
87  {};
88 
89  /**
90  * Constructor.
91  *
92  * From given minimal and maximal vertex.
93  */
94  BoundingBox(const Point &min, const Point &max)
95  : min_vertex_(min), max_vertex_(max)
96  {
97  ASSERT( arma::min( min <= max ) ).error("Wrong coordinates in constructor.");
98  };
99 
100  /**
101  * Constructor.
102  *
103  * Make bounding box for set of points.
104  */
105  BoundingBox(const vector<Point> &points);
106 
107  /**
108  * Set maximum in given axis.
109  */
110  void set_max(unsigned int axis, double max) {
111  ASSERT_LT( axis , dimension);
112  ASSERT_LE( min(axis) , max);
113  max_vertex_[axis] = max;
114  }
115 
116  /**
117  * Set minimum on given axis.
118  */
119  void set_min(unsigned int axis, double min) {
120  ASSERT_LT(axis, dimension);
121  ASSERT_LE(min , max(axis));
122  min_vertex_[axis] = min;
123  }
124 
125  /**
126  * Return minimal vertex of the bounding box.
127  */
128  const Point &min() const {
129  return min_vertex_;
130  }
131 
132  /**
133  * Return maximal vertex of the bounding box.
134  */
135  const Point &max() const {
136  return max_vertex_;
137  }
138 
139  /**
140  * Return minimal value on given axis.
141  */
142  double min(unsigned int axis) const {
143  return min()[axis];
144  }
145 
146  /**
147  * Return maximal value on given axis.
148  */
149  double max(unsigned int axis) const {
150  return max()[axis];
151  }
152 
153 
154  /**
155  * Return size of the box in given axis.
156  */
157  double size(unsigned int axis) const {
158  return max()[axis] - min()[axis];
159  }
160 
161 
162  /**
163  * Return center of the bounding box.
164  */
165  Point center() const {
166  return (max_vertex_ + min_vertex_) / 2.0;
167  }
168 
169  /**
170  * Return center of projection of the bounding box to given @p axis.
171  * Axis coding is: 0 - axis x, 1 - axis y, 2 - axis z.
172  */
173  double projection_center(unsigned int axis) const {
174  ASSERT_LT(axis, dimension).error("Invalid axis!\n");
175  return (max_vertex_[axis] + min_vertex_[axis])/2;
176  }
177 
178  /**
179  * Returns true is the box element contains @p point
180  *
181  * @param point Testing point
182  * @return True if box element contains point
183  */
184  bool contains_point(const Point &point) const
185  {
186  for (unsigned int i=0; i<dimension; i++) {
187  if ((point(i) + epsilon < min_vertex_(i)) ||
188  (point(i) > epsilon + max_vertex_(i))) return false;
189  }
190  return true;
191  }
192 
193  /**
194  * Returns true if two bounding boxes have intersection.
195  *
196  * This serves as an estimate of intersection of elements.
197  * To make it safe (do not exclude possible intersection) for
198  * 1d and 2d elements aligned with axes, we use some tolerance.
199  * Since this tolerance is fixed, there could be problem with
200  * highly refined meshes (get false positive result).
201  */
202  bool intersect(const BoundingBox &b2) const
203  {
204  for (unsigned int i=0; i<dimension; i++) {
205  if ( (min_vertex_(i) > b2.max_vertex_(i) + epsilon) ||
206  (b2.min_vertex_(i) > max_vertex_(i) + epsilon ) ) return false;
207  }
208  return true;
209  }
210 
211  /**
212  * Returns true if projection of the box to @p axis is an interval
213  * less then (with tolerance) to given @p value.
214  */
215  bool projection_lt(unsigned int axis, double value) const
216  {
217  return max_vertex_(axis) + epsilon < value;
218  }
219 
220  /**
221  * Returns true if projection of the box to @p axis is an interval
222  * greater then (with tolerance) to given @p value.
223  */
224  bool projection_gt(unsigned int axis, double value) const
225  {
226  return min_vertex_(axis) - epsilon > value;
227  }
228 
229  /**
230  * Split box into two boxes along @p axis by the
231  * plane going through @p splitting_point on the axis.
232  */
233  void split(unsigned int axis, double splitting_point,
234  BoundingBox &left, BoundingBox &right ) const
235  {
236  ASSERT_LT(axis , dimension);
237  if (min_vertex_[axis] <= splitting_point && splitting_point <= max_vertex_[axis] ) {
238  left = *this;
239  right = *this;
240  left.max_vertex_[axis] = splitting_point;
241  right.min_vertex_[axis] = splitting_point;
242  } else {
243  THROW( ExcSplitting() << EI_interval_left(min_vertex_[axis])
244  << EI_interval_right(max_vertex_[axis])
245  << EI_split_point(splitting_point) );
246  }
247  }
248 
249  /**
250  * Expand bounding box to contain also given @p point.
251  */
252  void expand(const Point &point) {
253  for(unsigned int j=0; j<dimension; j++) {
254  // parameters of min and max functions must be in correct order, vertices can be set to NaN (see default constructor)
255  min_vertex_(j) = std::min( point[j], min_vertex_[j] );
256  max_vertex_(j) = std::max( point[j], max_vertex_[j] );
257  }
258  }
259 
260  /**
261  * Expand bounding box to contain also given @p box.
262  */
263  void expand(const BoundingBox &box) {
264  for(unsigned int j=0; j<dimension; j++) {
265  // parameters of min and max functions must be in correct order, vertices can be set to NaN (see default constructor)
266  min_vertex_[j] = std::min( box.min_vertex_[j], min_vertex_[j] );
267  max_vertex_[j] = std::max( box.max_vertex_[j], max_vertex_[j] );
268  }
269  }
270 
271  /**
272  * Return index of the axis in which the box has longest projection.
273  */
274  unsigned char longest_axis() const {
275  auto diff=max_vertex_ - min_vertex_;
276  return (diff[1] > diff[0])
277  ? ( diff[2] > diff[1] ? 2 : 1 )
278  : ( diff[2] > diff[0] ? 2 : 0 );
279  }
280 
281  /**
282  * Project point to bounding box.
283  *
284  * If point is in bounding box, returns its.
285  */
286  Point project_point(const Point &point) const {
287  Point projected_point;
288  for (unsigned int i=0; i<dimension; ++i) {
289  if ( projection_gt(i, point[i]) ) projected_point[i] = min_vertex_(i);
290  else if ( projection_lt(i, point[i]) ) projected_point[i] = max_vertex_(i);
291  else projected_point[i] = point[i];
292  }
293 
294  return projected_point;
295  }
296 
297  /**
298  * Return size of the axis in which the box has longest projection.
299  */
300  inline double longest_size() const {
301  return size( longest_axis() );
302  }
303 
304 
305 private:
306  /// minimal coordinates of bounding box
308  /// maximal coordinates of bounding box
310 };
311 
312 /// Overloads output operator for box.
313 inline ostream &operator<<(ostream &stream, const BoundingBox &box) {
314  stream << "Box("
315  << box.min(0) << " "
316  << box.min(1) << " "
317  << box.min(2) << "; "
318  << box.max(0) << " "
319  << box.max(1) << " "
320  << box.max(2) << ")";
321  return stream;
322 }
323 
324 #endif /* BOX_ELEMENT_HH_ */
#define ASSERT(expr)
Definition: asserts.hh:351
#define ASSERT_LT(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
Definition: asserts.hh:301
#define ASSERT_LE(a, b)
Definition of comparative assert macro (Less or Equal) only for debug mode.
Definition: asserts.hh:309
ostream & operator<<(ostream &stream, const BoundingBox &box)
Overloads output operator for box.
Bounding box in 3d ambient space.
Definition: bounding_box.hh:53
const Point & min() const
bool projection_lt(unsigned int axis, double value) const
double longest_size() const
unsigned char longest_axis() const
Point center() const
Point min_vertex_
minimal coordinates of bounding box
TYPEDEF_ERR_INFO(EI_interval_right, double)
bool projection_gt(unsigned int axis, double value) const
void set_min(unsigned int axis, double min)
void split(unsigned int axis, double splitting_point, BoundingBox &left, BoundingBox &right) const
bool contains_point(const Point &point) const
double min(unsigned int axis) const
BoundingBox(const Point &min)
Definition: bounding_box.hh:85
TYPEDEF_ERR_INFO(EI_split_point, double)
double projection_center(unsigned int axis) const
const Point & max() const
double max(unsigned int axis) const
Point max_vertex_
maximal coordinates of bounding box
TYPEDEF_ERR_INFO(EI_interval_left, double)
Point project_point(const Point &point) const
static const double epsilon
stabilization parameter
Definition: bounding_box.hh:65
DECLARE_EXCEPTION(ExcSplitting,<< "Split point "<< EI_split_point::val<< "out of bounds: <"<< EI_interval_left::val<< ", "<< EI_interval_right::val<< ">\n")
double size(unsigned int axis) const
void set_max(unsigned int axis, double max)
void expand(const BoundingBox &box)
Space< dimension >::Point Point
Currently we assume.
Definition: bounding_box.hh:67
void expand(const Point &point)
BoundingBox(const Point &min, const Point &max)
Definition: bounding_box.hh:94
bool intersect(const BoundingBox &b2) const
Definition: point.hh:40
Global macros to enhance readability and debugging, general constants.
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
Definition: exceptions.hh:53
static constexpr bool value
Definition: json.hpp:87