DOLFIN
DOLFIN C++ interface
BoundingBoxTree3D.h
1 // Copyright (C) 2013 Anders Logg
2 //
3 // This file is part of DOLFIN.
4 //
5 // DOLFIN is free software: you can redistribute it and/or modify
6 // it under the terms of the GNU Lesser General Public License as published by
7 // the Free Software Foundation, either version 3 of the License, or
8 // (at your option) any later version.
9 //
10 // DOLFIN is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 // GNU Lesser General Public License for more details.
14 //
15 // You should have received a copy of the GNU Lesser General Public License
16 // along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
17 //
18 // First added: 2013-04-09
19 // Last changed: 2013-11-30
20 
21 #ifndef __BOUNDING_BOX_TREE_3D_H
22 #define __BOUNDING_BOX_TREE_3D_H
23 
24 #include <algorithm>
25 #include <vector>
26 #include <dolfin/common/constants.h>
27 #include "GenericBoundingBoxTree.h"
28 
29 namespace dolfin
30 {
31 
33 
35  {
36  protected:
37 
39 
41 
43  struct less_x_bbox
44  {
46  const std::vector<double>& bboxes;
47 
49  less_x_bbox(const std::vector<double>& bboxes): bboxes(bboxes) {}
50 
52  inline bool operator()(unsigned int i, unsigned int j)
53  {
54  const double* bi = bboxes.data() + 6*i;
55  const double* bj = bboxes.data() + 6*j;
56  return bi[0] + bi[3] < bj[0] + bj[3];
57  }
58  };
59 
61  struct less_y_bbox
62  {
64  const std::vector<double>& bboxes;
65 
67  less_y_bbox(const std::vector<double>& bboxes): bboxes(bboxes) {}
68 
70  inline bool operator()(unsigned int i, unsigned int j)
71  {
72  const double* bi = bboxes.data() + 6*i;
73  const double* bj = bboxes.data() + 6*j;
74  return bi[1] + bi[4] < bj[1] + bj[4];
75  }
76  };
77 
79  struct less_z_bbox
80  {
82  const std::vector<double>& bboxes;
83 
85  less_z_bbox(const std::vector<double>& bboxes): bboxes(bboxes) {}
86 
88  inline bool operator()(unsigned int i, unsigned int j)
89  {
90  const double* bi = bboxes.data() + 6*i;
91  const double* bj = bboxes.data() + 6*j;
92  return bi[2] + bi[5] < bj[2] + bj[5];
93  }
94  };
95 
97  std::size_t gdim() const { return 3; }
98 
100  const double* get_bbox_coordinates(unsigned int node) const
101  {
102  return _bbox_coordinates.data() + 6*node;
103  }
104 
106  bool point_in_bbox(const double* x, const unsigned int node) const
107  {
108  const double* b = _bbox_coordinates.data() + 6*node;
109  const double eps0 = DOLFIN_EPS_LARGE*(b[3] - b[0]);
110  const double eps1 = DOLFIN_EPS_LARGE*(b[4] - b[1]);
111  const double eps2 = DOLFIN_EPS_LARGE*(b[5] - b[2]);
112  return (b[0] - eps0 <= x[0] && x[0] <= b[3] + eps0 &&
113  b[1] - eps1 <= x[1] && x[1] <= b[4] + eps1 &&
114  b[2] - eps2 <= x[2] && x[2] <= b[5] + eps2);
115  }
116 
118  bool bbox_in_bbox(const double* a, unsigned int node) const
119  {
120  const double* b = _bbox_coordinates.data() + 6*node;
121  const double eps0 = DOLFIN_EPS_LARGE*(b[3] - b[0]);
122  const double eps1 = DOLFIN_EPS_LARGE*(b[4] - b[1]);
123  const double eps2 = DOLFIN_EPS_LARGE*(b[5] - b[2]);
124  return (b[0] - eps0 <= a[3] && a[0] <= b[3] + eps0 &&
125  b[1] - eps1 <= a[4] && a[1] <= b[4] + eps1 &&
126  b[2] - eps2 <= a[5] && a[2] <= b[5] + eps2);
127  }
128 
130  double compute_squared_distance_bbox(const double* x,
131  unsigned int node) const
132  {
133  // Note: Some else-if might be in order here but I assume the
134  // compiler can do a better job at optimizing/parallelizing this
135  // version. This is also the way the algorithm is presented in
136  // Ericsson.
137 
138  const double* b = _bbox_coordinates.data() + 6*node;
139  double r2 = 0.0;
140 
141  if (x[0] < b[0]) r2 += (x[0] - b[0])*(x[0] - b[0]);
142  if (x[0] > b[3]) r2 += (x[0] - b[3])*(x[0] - b[3]);
143  if (x[1] < b[1]) r2 += (x[1] - b[1])*(x[1] - b[1]);
144  if (x[1] > b[4]) r2 += (x[1] - b[4])*(x[1] - b[4]);
145  if (x[2] < b[2]) r2 += (x[2] - b[2])*(x[2] - b[2]);
146  if (x[2] > b[5]) r2 += (x[2] - b[5])*(x[2] - b[5]);
147 
148  return r2;
149  }
150 
152  double compute_squared_distance_point(const double* x,
153  unsigned int node) const
154  {
155  const double* p = _bbox_coordinates.data() + 6*node;
156  return ((x[0] - p[0])*(x[0] - p[0]) +
157  (x[1] - p[1])*(x[1] - p[1]) +
158  (x[2] - p[2])*(x[2] - p[2]));
159  }
160 
162  void compute_bbox_of_bboxes(double* bbox,
163  std::size_t& axis,
164  const std::vector<double>& leaf_bboxes,
165  const std::vector<unsigned int>::iterator& begin,
166  const std::vector<unsigned int>::iterator& end)
167  {
168  typedef std::vector<unsigned int>::const_iterator iterator;
169 
170  // Get coordinates for first box
171  iterator it = begin;
172  const double* b = leaf_bboxes.data() + 6*(*it);
173  bbox[0] = b[0];
174  bbox[1] = b[1];
175  bbox[2] = b[2];
176  bbox[3] = b[3];
177  bbox[4] = b[4];
178  bbox[5] = b[5];
179 
180  // Compute min and max over remaining boxes
181  for (; it != end; ++it)
182  {
183  const double* b = leaf_bboxes.data() + 6*(*it);
184  if (b[0] < bbox[0]) bbox[0] = b[0];
185  if (b[1] < bbox[1]) bbox[1] = b[1];
186  if (b[2] < bbox[2]) bbox[2] = b[2];
187  if (b[3] > bbox[3]) bbox[3] = b[3];
188  if (b[4] > bbox[4]) bbox[4] = b[4];
189  if (b[5] > bbox[5]) bbox[5] = b[5];
190  }
191 
192  // Compute longest axis
193  const double x = bbox[3] - bbox[0];
194  const double y = bbox[4] - bbox[1];
195  const double z = bbox[5] - bbox[2];
196 
197  if (x > y && x > z)
198  axis = 0;
199  else if (y > z)
200  axis = 1;
201  else
202  axis = 2;
203  }
204 
206  void compute_bbox_of_points(double* bbox,
207  std::size_t& axis,
208  const std::vector<Point>& points,
209  const std::vector<unsigned int>::iterator& begin,
210  const std::vector<unsigned int>::iterator& end)
211  {
212  typedef std::vector<unsigned int>::const_iterator iterator;
213 
214  // Get coordinates for first point
215  iterator it = begin;
216  const double* p = points[*it].coordinates();
217  bbox[0] = p[0];
218  bbox[1] = p[1];
219  bbox[2] = p[2];
220  bbox[3] = p[0];
221  bbox[4] = p[1];
222  bbox[5] = p[2];
223 
224  // Compute min and max over remaining points
225  for (++it; it != end; ++it)
226  {
227  const double* p = points[*it].coordinates();
228  if (p[0] < bbox[0]) bbox[0] = p[0];
229  if (p[1] < bbox[1]) bbox[1] = p[1];
230  if (p[2] < bbox[2]) bbox[2] = p[2];
231  if (p[0] > bbox[3]) bbox[3] = p[0];
232  if (p[1] > bbox[4]) bbox[4] = p[1];
233  if (p[2] > bbox[5]) bbox[5] = p[2];
234  }
235 
236  // Compute longest axis
237  const double x = bbox[3] - bbox[0];
238  const double y = bbox[4] - bbox[1];
239  const double z = bbox[5] - bbox[2];
240 
241  if (x > y && x > z)
242  axis = 0;
243  else if (y > z)
244  axis = 1;
245  else
246  axis = 2;
247  }
248 
250  void sort_bboxes(std::size_t axis,
251  const std::vector<double>& leaf_bboxes,
252  const std::vector<unsigned int>::iterator& begin,
253  const std::vector<unsigned int>::iterator& middle,
254  const std::vector<unsigned int>::iterator& end)
255  {
256  switch (axis)
257  {
258  case 0:
259  std::nth_element(begin, middle, end, less_x_bbox(leaf_bboxes));
260  break;
261  case 1:
262  std::nth_element(begin, middle, end, less_y_bbox(leaf_bboxes));
263  break;
264  default:
265  std::nth_element(begin, middle, end, less_z_bbox(leaf_bboxes));
266  }
267  }
268 
269  };
270 
271 }
272 
273 #endif
double compute_squared_distance_point(const double *x, unsigned int node) const
Compute squared distance between point and point.
Definition: BoundingBoxTree3D.h:152
void compute_bbox_of_bboxes(double *bbox, std::size_t &axis, const std::vector< double > &leaf_bboxes, const std::vector< unsigned int >::iterator &begin, const std::vector< unsigned int >::iterator &end)
Compute bounding box of bounding boxes.
Definition: BoundingBoxTree3D.h:162
std::size_t gdim() const
Return geometric dimension.
Definition: BoundingBoxTree3D.h:97
Definition: adapt.h:29
Definition: GenericBoundingBoxTree.h:40
void begin(std::string msg,...)
Begin task (increase indentation level)
Definition: log.cpp:153
less_x_bbox(const std::vector< double > &bboxes)
Constructor.
Definition: BoundingBoxTree3D.h:49
bool operator()(unsigned int i, unsigned int j)
Comparison operator.
Definition: BoundingBoxTree3D.h:88
Less than operator in z-direction.
Definition: BoundingBoxTree3D.h:79
less_y_bbox(const std::vector< double > &bboxes)
Constructor.
Definition: BoundingBoxTree3D.h:67
const std::vector< double > & bboxes
Bounding boxes.
Definition: BoundingBoxTree3D.h:64
const double * get_bbox_coordinates(unsigned int node) const
Return bounding box coordinates for node.
Definition: BoundingBoxTree3D.h:100
const std::vector< double > & bboxes
Bounding boxes.
Definition: BoundingBoxTree3D.h:82
void end()
End task (decrease indentation level)
Definition: log.cpp:168
Less than operator in y-direction.
Definition: BoundingBoxTree3D.h:61
bool operator()(unsigned int i, unsigned int j)
Comparison operator.
Definition: BoundingBoxTree3D.h:52
double compute_squared_distance_bbox(const double *x, unsigned int node) const
Compute squared distance between point and bounding box.
Definition: BoundingBoxTree3D.h:130
bool point_in_bbox(const double *x, const unsigned int node) const
Check whether point (x) is in bounding box (node)
Definition: BoundingBoxTree3D.h:106
void sort_bboxes(std::size_t axis, const std::vector< double > &leaf_bboxes, const std::vector< unsigned int >::iterator &begin, const std::vector< unsigned int >::iterator &middle, const std::vector< unsigned int >::iterator &end)
Sort leaf bounding boxes along given axis.
Definition: BoundingBoxTree3D.h:250
less_z_bbox(const std::vector< double > &bboxes)
Constructor.
Definition: BoundingBoxTree3D.h:85
bool operator()(unsigned int i, unsigned int j)
Comparison operator.
Definition: BoundingBoxTree3D.h:70
Comparison operators for sorting of bounding boxes.
Definition: BoundingBoxTree3D.h:43
bool bbox_in_bbox(const double *a, unsigned int node) const
Check whether bounding box (a) collides with bounding box (node)
Definition: BoundingBoxTree3D.h:118
std::vector< double > _bbox_coordinates
List of bounding box coordinates.
Definition: GenericBoundingBoxTree.h:118
void compute_bbox_of_points(double *bbox, std::size_t &axis, const std::vector< Point > &points, const std::vector< unsigned int >::iterator &begin, const std::vector< unsigned int >::iterator &end)
Compute bounding box of points.
Definition: BoundingBoxTree3D.h:206
const std::vector< double > & bboxes
Bounding boxes.
Definition: BoundingBoxTree3D.h:46
Specialization of bounding box implementation to 3D.
Definition: BoundingBoxTree3D.h:34