ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BoundingBox.ipp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file BoundingBox.ipp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2018-2019 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 template <typename entity_t, typename value_t, size_t DIM>
11  const entity_t* entity, const VertexType& vmin, const VertexType& vmax)
12  : m_entity(entity),
13  m_vmin(vmin),
14  m_vmax(vmax),
15  m_center((vmin + vmax) / 2.),
16  m_width(vmax - vmin),
17  m_iwidth(1 / m_width) {}
18 
19 template <typename entity_t, typename value_t, size_t DIM>
21  const entity_t* entity, const VertexType& center, const Size& size)
22  : m_entity(entity),
23  m_vmin(center - size.get() * 0.5),
24  m_vmax(center + size.get() * 0.5),
25  m_center(center),
26  m_width(size.get()),
27  m_iwidth(1 / m_width) {}
28 
29 template <typename entity_t, typename value_t, size_t DIM>
31  const std::vector<self_t*>& boxes, vertex_array_type envelope)
32  : m_entity(nullptr) {
33  assert(boxes.size() > 1);
34 
35  for (size_t i = 0; i < boxes.size(); i++) {
36  if (i < boxes.size() - 1) {
37  // set next on i to i+1
38  boxes[i]->setSkip(boxes[i + 1]);
39  } else {
40  // make sure last is set to nullptr, this marks end
41  // boxes[i]->m_next = nullptr;
42  boxes[i]->setSkip(nullptr);
43  }
44  }
45 
46  m_left_child = boxes.front();
47  m_right_child = boxes.back();
48  m_skip = nullptr;
49 
50  std::tie(m_vmin, m_vmax) = wrap(boxes, envelope);
51 
52  m_center = (m_vmin + m_vmax) / 2.;
53  m_width = m_vmax - m_vmin;
54  m_iwidth = 1 / m_width;
55 }
56 
57 template <typename entity_t, typename value_t, size_t DIM>
58 std::pair<
62  const std::vector<const self_t*>& boxes, vertex_array_type envelope) {
63  assert(boxes.size() > 1);
64  // figure out extent of boxes
65  // use array for Eigen coefficient wise min/max
66  vertex_array_type vmax(
67  vertex_array_type::Constant(std::numeric_limits<value_type>::lowest()));
68  vertex_array_type vmin(
69  vertex_array_type::Constant(std::numeric_limits<value_type>::max()));
70 
71  for (size_t i = 0; i < boxes.size(); i++) {
72  vmin = vmin.min(boxes[i]->min().array());
73  vmax = vmax.max(boxes[i]->max().array());
74  }
75 
76  vmax += envelope;
77  vmin -= envelope;
78 
79  return {vmin, vmax};
80 }
81 
82 template <typename entity_t, typename value_t, size_t DIM>
83 std::pair<
87  const std::vector<self_t*>& boxes, vertex_array_type envelope) {
88  assert(boxes.size() > 1);
89  std::vector<const self_t*> box_ptrs;
90  box_ptrs.reserve(boxes.size());
91  std::transform(boxes.begin(), boxes.end(), std::back_inserter(box_ptrs),
92  [](const auto* box) { return box; });
93  return wrap(box_ptrs, envelope);
94 }
95 
96 template <typename entity_t, typename value_t, size_t DIM>
97 std::pair<
101  const std::vector<self_t>& boxes, vertex_array_type envelope) {
102  assert(boxes.size() > 1);
103  std::vector<const self_t*> box_ptrs;
104  box_ptrs.reserve(boxes.size());
105  std::transform(boxes.begin(), boxes.end(), std::back_inserter(box_ptrs),
106  [](auto& box) { return &box; });
107  return wrap(box_ptrs, envelope);
108 }
109 
110 template <typename entity_t, typename value_t, size_t DIM>
112  const VertexType& point) const {
113  vertex_array_type t = (point - m_vmin).array() * m_iwidth;
114  return t.minCoeff() >= 0 && t.maxCoeff() < 1;
115 }
116 
117 template <typename entity_t, typename value_t, size_t DIM>
119  const Ray<value_type, DIM>& ray) const {
120  const VertexType& origin = ray.origin();
121  const vertex_array_type& idir = ray.idir();
122 
123  // Calculate the intersect distances with the min and max planes along the ray
124  // direction, from the ray origin. See Ch VII.5 Fig.1 in [1].
125  // This is done in all dimensions at the same time:
126  vertex_array_type t0s = (m_vmin - origin).array() * idir;
127  vertex_array_type t1s = (m_vmax - origin).array() * idir;
128 
129  // Calculate the component wise min/max between the t0s and t1s
130  // this is non-compliant with IEEE-754-2008, NaN gets propagated through
131  // http://eigen.tuxfamily.org/bz/show_bug.cgi?id=564
132  // this means that rays parallel to boundaries might not be considered
133  // to intersect.
134  vertex_array_type tsmaller = t0s.min(t1s);
135  vertex_array_type tbigger = t0s.max(t1s);
136 
137  // extract largest and smallest component of the component wise extrema
138  value_type tmin = tsmaller.maxCoeff();
139  value_type tmax = tbigger.minCoeff();
140 
141  // If tmin is smaller than tmax and tmax is positive, then the box is in
142  // positive ray direction, and the ray intersects the box.
143  return tmin < tmax && tmax > 0.0;
144 }
145 
146 template <typename entity_t, typename value_t, size_t DIM>
147 template <size_t sides>
149  const Frustum<value_type, DIM, sides>& fr) const {
150  const auto& normals = fr.normals();
151  // Transform vmin and vmax into the coordinate system, at which the frustum is
152  // located at the coordinate origin.
153  const vertex_array_type fr_vmin = m_vmin - fr.origin();
154  const vertex_array_type fr_vmax = m_vmax - fr.origin();
155 
156  // For each plane, find the p-vertex, which is the vertex that is at the
157  // furthest distance from the plane *along* it's normal direction.
158  // See Fig. 2 in [2].
159  VertexType p_vtx;
160  // for loop, we could eliminate this, probably,
161  // but sides+1 is known at compile time, so the compiler
162  // will most likely unroll the loop
163  for (size_t i = 0; i < sides + 1; i++) {
164  const VertexType& normal = normals[i];
165 
166  // for AABBs, take the component from the min vertex, if the normal
167  // component is negative, else take the component from the max vertex.
168  p_vtx = (normal.array() < 0).template cast<value_type>() * fr_vmin +
169  (normal.array() >= 0).template cast<value_type>() * fr_vmax;
170 
171  // Check if the p-vertex is at positive or negative direction along the
172  // If the p vertex is along negative normal direction *once*, the box is
173  // outside the frustum, and we can terminate early.
174  if (p_vtx.dot(normal) < 0) {
175  // p vertex is outside on this plane, box must be outside
176  return false;
177  }
178  }
179 
180  // If we get here, no p-vertex was outside, so box intersects or is
181  // contained. We don't care, so report 'intersect'
182  return true;
183 }
184 
185 template <typename entity_t, typename value_t, size_t DIM>
187  self_t* skip) {
188  // set next on this
189  m_skip = skip;
190  // find last child and set its skip
191  if (m_right_child != nullptr) {
192  m_right_child->setSkip(skip);
193  }
194 }
195 
196 template <typename entity_t, typename value_t, size_t DIM>
199  return m_left_child;
200 }
201 
202 template <typename entity_t, typename value_t, size_t DIM>
205  return m_skip;
206 }
207 
208 template <typename entity_t, typename value_t, size_t DIM>
210  return m_entity != nullptr;
211 }
212 
213 template <typename entity_t, typename value_t, size_t DIM>
215  const {
216  return m_entity;
217 }
218 
219 template <typename entity_t, typename value_t, size_t DIM>
221  const entity_t* entity) {
222  m_entity = entity;
223 }
224 
225 template <typename entity_t, typename value_t, size_t DIM>
228  return m_center;
229 }
230 
231 template <typename entity_t, typename value_t, size_t DIM>
234  return m_vmin;
235 }
236 
237 template <typename entity_t, typename value_t, size_t DIM>
240  return m_vmax;
241 }
242 
243 template <typename entity_t, typename value_t, size_t DIM>
245  std::ostream& os) const {
246  os << "AABB(ctr=(";
247 
248  for (size_t i = 0; i < DIM; i++) {
249  if (i > 0) {
250  os << ", ";
251  }
252  os << m_center[i];
253  }
254 
255  os << ") vmin=(";
256  for (size_t i = 0; i < DIM; i++) {
257  if (i > 0) {
258  os << ", ";
259  }
260  os << m_vmin[i];
261  }
262 
263  os << ") vmax=(";
264 
265  for (size_t i = 0; i < DIM; i++) {
266  if (i > 0) {
267  os << ", ";
268  }
269  os << m_vmax[i];
270  }
271 
272  os << "))";
273 
274  return os;
275 }
276 
277 template <typename entity_t, typename value_t, size_t DIM>
278 template <size_t D, std::enable_if_t<D == 3, int>>
279 std::pair<
283  const transform_type& trf) const {
284  // we need to enumerate all the vertices, transform,
285  // and then recalculate min and max
286 
287  std::array<VertexType, 8> vertices({{
288  {m_vmin.x(), m_vmin.y(), m_vmin.z()},
289  {m_vmin.x(), m_vmax.y(), m_vmin.z()},
290  {m_vmax.x(), m_vmax.y(), m_vmin.z()},
291  {m_vmax.x(), m_vmin.y(), m_vmin.z()},
292  {m_vmin.x(), m_vmin.y(), m_vmax.z()},
293  {m_vmin.x(), m_vmax.y(), m_vmax.z()},
294  {m_vmax.x(), m_vmax.y(), m_vmax.z()},
295  {m_vmax.x(), m_vmin.y(), m_vmax.z()},
296  }});
297 
298  VertexType vmin = trf * vertices[0];
299  VertexType vmax = trf * vertices[0];
300 
301  for (size_t i = 1; i < 8; i++) {
302  const VertexType vtx = trf * vertices[i];
303  vmin = vmin.cwiseMin(vtx);
304  vmax = vmax.cwiseMax(vtx);
305  }
306 
307  return {vmin, vmax};
308 }
309 
310 template <typename entity_t, typename value_t, size_t DIM>
311 template <size_t D, std::enable_if_t<D == 2, int>>
312 std::pair<
316  const transform_type& trf) const {
317  // we need to enumerate all the vertices, transform,
318  // and then recalculate min and max
319 
320  std::array<VertexType, 4> vertices({{{m_vmin.x(), m_vmin.y()},
321  {m_vmin.x(), m_vmax.y()},
322  {m_vmax.x(), m_vmax.y()},
323  {m_vmax.x(), m_vmin.y()}}});
324 
325  VertexType vmin = trf * vertices[0];
326  VertexType vmax = trf * vertices[0];
327 
328  for (size_t i = 1; i < 4; i++) {
329  const VertexType vtx = trf * vertices[i];
330  vmin = vmin.cwiseMin(vtx);
331  vmax = vmax.cwiseMax(vtx);
332  }
333 
334  return {vmin, vmax};
335 }
336 
337 template <typename entity_t, typename value_t, size_t DIM>
339  const transform_type& trf) {
340  std::tie(m_vmin, m_vmax) = transformVertices(trf);
341 }
342 
343 template <typename entity_t, typename value_t, size_t DIM>
346  const transform_type& trf) const {
347  VertexType vmin, vmax;
348  std::tie(vmin, vmax) = transformVertices(trf);
349  return self_t(m_entity, vmin, vmax);
350 }
351 
352 template <typename entity_t, typename value_t, size_t DIM>
353 template <size_t D, std::enable_if_t<D == 3, int>>
355  IVisualization& helper, std::array<int, 3> color,
356  const transform_type& trf) const {
357  static_assert(DIM == 3, "PLY output only supported in 3D");
358 
359  const VertexType& vmin = m_vmin;
360  const VertexType& vmax = m_vmax;
361 
362  auto write = [&](const VertexType& a, const VertexType& b,
363  const VertexType& c, const VertexType& d) {
364  helper.face(std::vector<VertexType>({trf * a, trf * b, trf * c, trf * d}),
365  color);
366  };
367 
368  write({vmin.x(), vmin.y(), vmin.z()}, {vmin.x(), vmax.y(), vmin.z()},
369  {vmin.x(), vmax.y(), vmax.z()}, {vmin.x(), vmin.y(), vmax.z()});
370 
371  write({vmax.x(), vmin.y(), vmin.z()}, {vmax.x(), vmax.y(), vmin.z()},
372  {vmax.x(), vmax.y(), vmax.z()}, {vmax.x(), vmin.y(), vmax.z()});
373 
374  write({vmin.x(), vmin.y(), vmin.z()}, {vmax.x(), vmin.y(), vmin.z()},
375  {vmax.x(), vmin.y(), vmax.z()}, {vmin.x(), vmin.y(), vmax.z()});
376 
377  write({vmin.x(), vmax.y(), vmin.z()}, {vmax.x(), vmax.y(), vmin.z()},
378  {vmax.x(), vmax.y(), vmax.z()}, {vmin.x(), vmax.y(), vmax.z()});
379 
380  write({vmin.x(), vmin.y(), vmin.z()}, {vmax.x(), vmin.y(), vmin.z()},
381  {vmax.x(), vmax.y(), vmin.z()}, {vmin.x(), vmax.y(), vmin.z()});
382 
383  write({vmin.x(), vmin.y(), vmax.z()}, {vmax.x(), vmin.y(), vmax.z()},
384  {vmax.x(), vmax.y(), vmax.z()}, {vmin.x(), vmax.y(), vmax.z()});
385 }
386 
387 template <typename entity_t, typename value_t, size_t DIM>
388 template <size_t D, std::enable_if_t<D == 2, int>>
390  std::ostream& os, value_type w, value_type h, value_type unit,
391  std::string label, std::string fillcolor) const {
392  static_assert(DIM == 2, "SVG is only supported in 2D");
393 
394  VertexType mid(w / 2., h / 2.);
395 
396  using transform_t = Eigen::Transform<value_t, DIM, Eigen::Affine>;
397 
398  transform_t trf = transform_t::Identity();
399  trf.translate(mid);
400  trf = trf * Eigen::Scaling(VertexType(1, -1));
401  trf.scale(unit);
402 
403  auto draw_point = [&](const VertexType& p_, std::string color, size_t r) {
404  VertexType p = trf * p_;
405  os << "<circle ";
406  os << "cx=\"" << p.x() << "\" cy=\"" << p.y() << "\" r=\"" << r << "\"";
407  os << " fill=\"" << color << "\"";
408  os << "/>\n";
409  };
410 
411  auto draw_rect = [&](const VertexType& center_, const VertexType& size_,
412  std::string color) {
413  VertexType size = size_ * unit;
414  VertexType center = trf * center_ - size * 0.5;
415 
416  os << "<rect ";
417  os << "x=\"" << center.x() << "\" y=\"" << center.y() << "\" ";
418  os << "width=\"" << size.x() << "\" height=\"" << size.y() << "\"";
419  os << " fill=\"" << color << "\"";
420  os << "/>\n";
421  };
422 
423  auto draw_text = [&](const VertexType& center_, std::string text,
424  std::string color, size_t size) {
425  VertexType center = trf * center_;
426  os << "<text dominant-baseline=\"middle\" text-anchor=\"middle\" ";
427  os << "fill=\"" << color << "\" font-size=\"" << size << "\" ";
428  os << "x=\"" << center.x() << "\" y=\"" << center.y() << "\">";
429  os << text << "</text>\n";
430  };
431 
432  draw_rect(m_center, m_width, fillcolor);
433  draw_point(m_vmin, "black", 2);
434  draw_point(m_vmax, "black", 2);
435  draw_text(m_center, label, "white", 10);
436 
437  return os;
438 }
439 
440 template <typename box_t>
441 box_t* octree_inner(std::vector<std::unique_ptr<box_t>>& store,
442  size_t max_depth,
443  typename box_t::vertex_array_type envelope,
444  const std::vector<box_t*>& lprims, size_t depth) {
445  using VertexType = typename box_t::VertexType;
446 
447  assert(lprims.size() > 0);
448  if (lprims.size() == 1) {
449  // just return
450  return lprims.front();
451  }
452 
453  if (depth >= max_depth) {
454  // just wrap them all up
455  auto bb = std::make_unique<box_t>(lprims, envelope);
456  store.push_back(std::move(bb));
457  return store.back().get();
458  }
459 
460  std::array<std::vector<box_t*>, 8> octants;
461  // calc center of boxes
462  VertexType vmin, vmax;
463  std::tie(vmin, vmax) = box_t::wrap(lprims);
464  VertexType glob_ctr = (vmin + vmax) / 2.;
465 
466  for (auto* box : lprims) {
467  VertexType ctr = box->center() - glob_ctr;
468  if (ctr.x() < 0 && ctr.y() < 0 && ctr.z() < 0) {
469  octants[0].push_back(box);
470  continue;
471  }
472  if (ctr.x() > 0 && ctr.y() < 0 && ctr.z() < 0) {
473  octants[1].push_back(box);
474  continue;
475  }
476  if (ctr.x() < 0 && ctr.y() > 0 && ctr.z() < 0) {
477  octants[2].push_back(box);
478  continue;
479  }
480  if (ctr.x() > 0 && ctr.y() > 0 && ctr.z() < 0) {
481  octants[3].push_back(box);
482  continue;
483  }
484 
485  if (ctr.x() < 0 && ctr.y() < 0 && ctr.z() > 0) {
486  octants[4].push_back(box);
487  continue;
488  }
489  if (ctr.x() > 0 && ctr.y() < 0 && ctr.z() > 0) {
490  octants[5].push_back(box);
491  continue;
492  }
493  if (ctr.x() < 0 && ctr.y() > 0 && ctr.z() > 0) {
494  octants[6].push_back(box);
495  continue;
496  }
497  if (ctr.x() > 0 && ctr.y() > 0 && ctr.z() > 0) {
498  octants[7].push_back(box);
499  continue;
500  }
501 
502  // not in any quadrant (numerics probably)
503  octants[0].push_back(box);
504  }
505 
506  std::vector<box_t*> sub_octs;
507  for (const auto& sub_prims : octants) {
508  if (sub_prims.size() <= 8) {
509  if (sub_prims.empty()) {
510  // done
511  } else if (sub_prims.size() == 1) {
512  sub_octs.push_back(sub_prims.front());
513  } else {
514  store.push_back(std::make_unique<box_t>(sub_prims, envelope));
515  sub_octs.push_back(store.back().get());
516  }
517  } else {
518  // recurse
519  sub_octs.push_back(
520  octree_inner(store, max_depth, envelope, sub_prims, depth + 1));
521  }
522  }
523 
524  if (sub_octs.size() == 1) {
525  return sub_octs.front();
526  }
527 
528  auto bb = std::make_unique<box_t>(sub_octs, envelope);
529  store.push_back(std::move(bb));
530  return store.back().get();
531 }
532 
533 template <typename box_t>
534 box_t* Acts::make_octree(std::vector<std::unique_ptr<box_t>>& store,
535  const std::vector<box_t*>& prims, size_t max_depth,
536  typename box_t::value_type envelope1) {
537  static_assert(box_t::dim == 3, "Octree can only be created in 3D");
538 
539  using vertex_array_type = typename box_t::vertex_array_type;
540 
541  vertex_array_type envelope(vertex_array_type::Constant(envelope1));
542 
543  box_t* top = octree_inner(store, max_depth, envelope, prims, 0);
544  return top;
545 }
546 
547 template <typename T, typename U, size_t V>
548 std::ostream& operator<<(std::ostream& os,
550  box.dump(os);
551  return os;
552 }