DOLFIN
DOLFIN C++ interface
BoundingBoxTree2D.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 // Modified by August Johansson 2016
19 //
20 // First added: 2013-05-02
21 // Last changed: 2016-11-15
22 
23 #ifndef __BOUNDING_BOX_TREE_2D_H
24 #define __BOUNDING_BOX_TREE_2D_H
25 
26 #include <algorithm>
27 #include <vector>
28 #include <dolfin/common/constants.h>
29 #include "GenericBoundingBoxTree.h"
30 
31 namespace dolfin
32 {
33 
35 
37  {
38  protected:
39 
42  struct less_x
43  {
45  const std::vector<double>& bboxes;
46 
48  less_x(const std::vector<double>& bboxes): bboxes(bboxes) {}
49 
51  inline bool operator()(unsigned int i, unsigned int j)
52  {
53  const double* bi = bboxes.data() + 4*i;
54  const double* bj = bboxes.data() + 4*j;
55  return bi[0] + bi[2] < bj[0] + bj[2];
56  }
57  };
58 
61  struct less_y
62  {
64  const std::vector<double>& bboxes;
65 
67  less_y(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() + 4*i;
73  const double* bj = bboxes.data() + 4*j;
74  return bi[1] + bi[3] < bj[1] + bj[3];
75  }
76  };
77 
79  std::size_t gdim() const { return 2; }
80 
82  const double* get_bbox_coordinates(unsigned int node) const
83  {
84  return _bbox_coordinates.data() + 4*node;
85  }
86 
88  bool point_in_bbox(const double* x, unsigned int node) const
89  {
90  const double* b = _bbox_coordinates.data() + 4*node;
91  const double eps0 = DOLFIN_EPS_LARGE*(b[2] - b[0]);
92  const double eps1 = DOLFIN_EPS_LARGE*(b[3] - b[1]);
93  return (b[0] - eps0 <= x[0] && x[0] <= b[2] + eps0 &&
94  b[1] - eps1 <= x[1] && x[1] <= b[3] + eps1);
95  }
96 
98  bool bbox_in_bbox(const double* a, unsigned int node) const
99  {
100  const double* b = _bbox_coordinates.data() + 4*node;
101  const double eps0 = DOLFIN_EPS_LARGE*(b[2] - b[0]);
102  const double eps1 = DOLFIN_EPS_LARGE*(b[3] - b[1]);
103  return (b[0] - eps0 <= a[2] && a[0] <= b[2] + eps0 &&
104  b[1] - eps1 <= a[3] && a[1] <= b[3] + eps1);
105  }
106 
108  double compute_squared_distance_bbox(const double* x,
109  unsigned int node) const
110  {
111  // Note: Some else-if might be in order here but I assume the
112  // compiler can do a better job at optimizing/parallelizing this
113  // version. This is also the way the algorithm is presented in
114  // Ericsson.
115 
116  const double* b = _bbox_coordinates.data() + 4*node;
117  double r2 = 0.0;
118 
119  if (x[0] < b[0]) r2 += (x[0] - b[0])*(x[0] - b[0]);
120  if (x[0] > b[2]) r2 += (x[0] - b[2])*(x[0] - b[2]);
121  if (x[1] < b[1]) r2 += (x[1] - b[1])*(x[1] - b[1]);
122  if (x[1] > b[3]) r2 += (x[1] - b[3])*(x[1] - b[3]);
123 
124  return r2;
125  }
126 
128  double compute_squared_distance_point(const double* x,
129  unsigned int node) const
130  {
131  const double* p = _bbox_coordinates.data() + 4*node;
132  return (x[0] - p[0])*(x[0] - p[0]) + (x[1] - p[1])*(x[1] - p[1]);
133  }
134 
136  void compute_bbox_of_bboxes(double* bbox,
137  std::size_t& axis,
138  const std::vector<double>& leaf_bboxes,
139  const std::vector<unsigned int>::iterator& begin,
140  const std::vector<unsigned int>::iterator& end)
141  {
142  typedef std::vector<unsigned int>::const_iterator iterator;
143 
144  // Get coordinates for first box
145  iterator it = begin;
146  const double* b = leaf_bboxes.data() + 4*(*it);
147  bbox[0] = b[0];
148  bbox[1] = b[1];
149  bbox[2] = b[2];
150  bbox[3] = b[3];
151 
153  for (++it; it != end; ++it)
154  {
155  const double* b = leaf_bboxes.data() + 4*(*it);
156  if (b[0] < bbox[0]) bbox[0] = b[0];
157  if (b[1] < bbox[1]) bbox[1] = b[1];
158  if (b[2] > bbox[2]) bbox[2] = b[2];
159  if (b[3] > bbox[3]) bbox[3] = b[3];
160  }
161 
162  // Compute longest axis
163  const double x = bbox[2] - bbox[0];
164  const double y = bbox[3] - bbox[1];
165 
166  if (x > y)
167  axis = 0;
168  else
169  axis = 1;
170  }
171 
173  void compute_bbox_of_points(double* bbox,
174  std::size_t& axis,
175  const std::vector<Point>& points,
176  const std::vector<unsigned int>::iterator& begin,
177  const std::vector<unsigned int>::iterator& end)
178  {
179  typedef std::vector<unsigned int>::const_iterator iterator;
180 
181  // Get coordinates for first point
182  iterator it = begin;
183  const double* p = points[*it].coordinates();
184  bbox[0] = p[0];
185  bbox[1] = p[1];
186  bbox[2] = p[0];
187  bbox[3] = p[1];
188 
189  // Compute min and max over remaining points
190  for (; it != end; ++it)
191  {
192  const double* p = points[*it].coordinates();
193  if (p[0] < bbox[0]) bbox[0] = p[0];
194  if (p[1] < bbox[1]) bbox[1] = p[1];
195  if (p[0] > bbox[2]) bbox[2] = p[0];
196  if (p[1] > bbox[3]) bbox[3] = p[1];
197  }
198 
199  // Compute longest axis
200  const double x = bbox[2] - bbox[0];
201  const double y = bbox[3] - bbox[1];
202 
203  if (x > y)
204  axis = 0;
205  else
206  axis = 1;
207  }
208 
210  void sort_bboxes(std::size_t axis,
211  const std::vector<double>& leaf_bboxes,
212  const std::vector<unsigned int>::iterator& begin,
213  const std::vector<unsigned int>::iterator& middle,
214  const std::vector<unsigned int>::iterator& end)
215  {
216  if (axis == 0)
217  std::nth_element(begin, middle, end, less_x(leaf_bboxes));
218  else
219  std::nth_element(begin, middle, end, less_y(leaf_bboxes));
220  }
221 
222  };
223 
224 }
225 
226 #endif
bool operator()(unsigned int i, unsigned int j)
Comparison operator.
Definition: BoundingBoxTree2D.h:51
bool operator()(unsigned int i, unsigned int j)
Comparison operator.
Definition: BoundingBoxTree2D.h:70
Definition: adapt.h:29
Definition: BoundingBoxTree2D.h:42
Definition: GenericBoundingBoxTree.h:40
std::size_t gdim() const
Return geometric dimension.
Definition: BoundingBoxTree2D.h:79
void begin(std::string msg,...)
Begin task (increase indentation level)
Definition: log.cpp:153
bool point_in_bbox(const double *x, unsigned int node) const
Check whether point (x) is in bounding box (node)
Definition: BoundingBoxTree2D.h:88
const double * get_bbox_coordinates(unsigned int node) const
Return bounding box coordinates for node.
Definition: BoundingBoxTree2D.h:82
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: BoundingBoxTree2D.h:210
bool bbox_in_bbox(const double *a, unsigned int node) const
Check whether bounding box (a) collides with bounding box (node)
Definition: BoundingBoxTree2D.h:98
void end()
End task (decrease indentation level)
Definition: log.cpp:168
double compute_squared_distance_point(const double *x, unsigned int node) const
Compute squared distance between point and point.
Definition: BoundingBoxTree2D.h:128
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: BoundingBoxTree2D.h:173
Specialization of bounding box implementation to 2D.
Definition: BoundingBoxTree2D.h:36
double compute_squared_distance_bbox(const double *x, unsigned int node) const
Compute squared distance between point and bounding box.
Definition: BoundingBoxTree2D.h:108
const std::vector< double > & bboxes
Bounding boxes.
Definition: BoundingBoxTree2D.h:45
less_x(const std::vector< double > &bboxes)
Constructor.
Definition: BoundingBoxTree2D.h:48
less_y(const std::vector< double > &bboxes)
Constructor.
Definition: BoundingBoxTree2D.h:67
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: BoundingBoxTree2D.h:136
std::vector< double > _bbox_coordinates
List of bounding box coordinates.
Definition: GenericBoundingBoxTree.h:118
const std::vector< double > & bboxes
Bounding boxes.
Definition: BoundingBoxTree2D.h:64
Definition: BoundingBoxTree2D.h:61