Analysis Software
Documentation for sPHENIX simulation 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* its 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  IVisualization3D& 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  const std::string& label, const 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_, const std::string& color,
404  size_t r) {
405  VertexType p = trf * p_;
406  os << "<circle ";
407  os << "cx=\"" << p.x() << "\" cy=\"" << p.y() << "\" r=\"" << r << "\"";
408  os << " fill=\"" << color << "\"";
409  os << "/>\n";
410  };
411 
412  auto draw_rect = [&](const VertexType& center_, const VertexType& size_,
413  const std::string& color) {
414  VertexType size = size_ * unit;
415  VertexType center = trf * center_ - size * 0.5;
416 
417  os << "<rect ";
418  os << "x=\"" << center.x() << "\" y=\"" << center.y() << "\" ";
419  os << "width=\"" << size.x() << "\" height=\"" << size.y() << "\"";
420  os << " fill=\"" << color << "\"";
421  os << "/>\n";
422  };
423 
424  auto draw_text = [&](const VertexType& center_, const std::string& text,
425  const std::string& color, size_t size) {
426  VertexType center = trf * center_;
427  os << "<text dominant-baseline=\"middle\" text-anchor=\"middle\" ";
428  os << "fill=\"" << color << "\" font-size=\"" << size << "\" ";
429  os << "x=\"" << center.x() << "\" y=\"" << center.y() << "\">";
430  os << text << "</text>\n";
431  };
432 
433  draw_rect(m_center, m_width, fillcolor);
434  draw_point(m_vmin, "black", 2);
435  draw_point(m_vmax, "black", 2);
436  draw_text(m_center, label, "white", 10);
437 
438  return os;
439 }
440 
441 template <typename box_t>
442 box_t* octree_inner(std::vector<std::unique_ptr<box_t>>& store,
443  size_t max_depth,
444  typename box_t::vertex_array_type envelope,
445  const std::vector<box_t*>& lprims, size_t depth) {
446  using VertexType = typename box_t::VertexType;
447 
448  assert(lprims.size() > 0);
449  if (lprims.size() == 1) {
450  // just return
451  return lprims.front();
452  }
453 
454  if (depth >= max_depth) {
455  // just wrap them all up
456  auto bb = std::make_unique<box_t>(lprims, envelope);
457  store.push_back(std::move(bb));
458  return store.back().get();
459  }
460 
461  std::array<std::vector<box_t*>, 8> octants;
462  // calc center of boxes
463  VertexType vmin, vmax;
464  std::tie(vmin, vmax) = box_t::wrap(lprims);
465  VertexType glob_ctr = (vmin + vmax) / 2.;
466 
467  for (auto* box : lprims) {
468  VertexType ctr = box->center() - glob_ctr;
469  if (ctr.x() < 0 && ctr.y() < 0 && ctr.z() < 0) {
470  octants[0].push_back(box);
471  continue;
472  }
473  if (ctr.x() > 0 && ctr.y() < 0 && ctr.z() < 0) {
474  octants[1].push_back(box);
475  continue;
476  }
477  if (ctr.x() < 0 && ctr.y() > 0 && ctr.z() < 0) {
478  octants[2].push_back(box);
479  continue;
480  }
481  if (ctr.x() > 0 && ctr.y() > 0 && ctr.z() < 0) {
482  octants[3].push_back(box);
483  continue;
484  }
485 
486  if (ctr.x() < 0 && ctr.y() < 0 && ctr.z() > 0) {
487  octants[4].push_back(box);
488  continue;
489  }
490  if (ctr.x() > 0 && ctr.y() < 0 && ctr.z() > 0) {
491  octants[5].push_back(box);
492  continue;
493  }
494  if (ctr.x() < 0 && ctr.y() > 0 && ctr.z() > 0) {
495  octants[6].push_back(box);
496  continue;
497  }
498  if (ctr.x() > 0 && ctr.y() > 0 && ctr.z() > 0) {
499  octants[7].push_back(box);
500  continue;
501  }
502 
503  // not in any quadrant (numerics probably)
504  octants[0].push_back(box);
505  }
506 
507  std::vector<box_t*> sub_octs;
508  for (const auto& sub_prims : octants) {
509  if (sub_prims.size() <= 8) {
510  if (sub_prims.empty()) {
511  // done
512  } else if (sub_prims.size() == 1) {
513  sub_octs.push_back(sub_prims.front());
514  } else {
515  store.push_back(std::make_unique<box_t>(sub_prims, envelope));
516  sub_octs.push_back(store.back().get());
517  }
518  } else {
519  // recurse
520  sub_octs.push_back(
521  octree_inner(store, max_depth, envelope, sub_prims, depth + 1));
522  }
523  }
524 
525  if (sub_octs.size() == 1) {
526  return sub_octs.front();
527  }
528 
529  auto bb = std::make_unique<box_t>(sub_octs, envelope);
530  store.push_back(std::move(bb));
531  return store.back().get();
532 }
533 
534 template <typename box_t>
535 box_t* Acts::make_octree(std::vector<std::unique_ptr<box_t>>& store,
536  const std::vector<box_t*>& prims, size_t max_depth,
537  typename box_t::value_type envelope1) {
538  static_assert(box_t::dim == 3, "Octree can only be created in 3D");
539 
541 
542  vertex_array_type envelope(vertex_array_type::Constant(envelope1));
543 
544  box_t* top = octree_inner(store, max_depth, envelope, prims, 0);
545  return top;
546 }
547 
548 template <typename T, typename U, size_t V>
549 std::ostream& Acts::operator<<(
550  std::ostream& os, const Acts::AxisAlignedBoundingBox<T, U, V>& box) {
551  return box.toStream(os);
552 }