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
29namespace dolfin
30{
31
33
35 {
36 protected:
37
39
41
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
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
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
Specialization of bounding box implementation to 3D.
Definition: BoundingBoxTree3D.h:35
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
double compute_squared_distance_point(const double *x, unsigned int node) const
Compute squared distance between point and point.
Definition: BoundingBoxTree3D.h:152
const double * get_bbox_coordinates(unsigned int node) const
Return bounding box coordinates for node.
Definition: BoundingBoxTree3D.h:100
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
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
double compute_squared_distance_bbox(const double *x, unsigned int node) const
Compute squared distance between point and bounding box.
Definition: BoundingBoxTree3D.h:130
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
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::size_t gdim() const
Return geometric dimension.
Definition: BoundingBoxTree3D.h:97
Definition: GenericBoundingBoxTree.h:41
std::vector< double > _bbox_coordinates
List of bounding box coordinates.
Definition: GenericBoundingBoxTree.h:118
Definition: adapt.h:30
void begin(std::string msg,...)
Begin task (increase indentation level)
Definition: log.cpp:153
void end()
End task (decrease indentation level)
Definition: log.cpp:168
Comparison operators for sorting of bounding boxes.
Definition: BoundingBoxTree3D.h:44
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:52
const std::vector< double > & bboxes
Bounding boxes.
Definition: BoundingBoxTree3D.h:46
Less than operator in y-direction.
Definition: BoundingBoxTree3D.h:62
const std::vector< double > & bboxes
Bounding boxes.
Definition: BoundingBoxTree3D.h:64
bool operator()(unsigned int i, unsigned int j)
Comparison operator.
Definition: BoundingBoxTree3D.h:70
less_y_bbox(const std::vector< double > &bboxes)
Constructor.
Definition: BoundingBoxTree3D.h:67
Less than operator in z-direction.
Definition: BoundingBoxTree3D.h:80
bool operator()(unsigned int i, unsigned int j)
Comparison operator.
Definition: BoundingBoxTree3D.h:88
less_z_bbox(const std::vector< double > &bboxes)
Constructor.
Definition: BoundingBoxTree3D.h:85
const std::vector< double > & bboxes
Bounding boxes.
Definition: BoundingBoxTree3D.h:82