Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-09-26 08:01:29

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008 
0009 #pragma once
0010 
0011 #include "Acts/Utilities/BoundingBox.hpp"
0012 
0013 #include <algorithm>
0014 
0015 template <typename entity_t, typename value_t, std::size_t DIM>
0016 Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::AxisAlignedBoundingBox(
0017     const entity_t* entity, const VertexType& vmin, const VertexType& vmax)
0018     : m_entity(entity),
0019       m_vmin(vmin),
0020       m_vmax(vmax),
0021       m_center((vmin + vmax) / 2.),
0022       m_width(vmax - vmin),
0023       m_iwidth(1 / m_width) {}
0024 
0025 template <typename entity_t, typename value_t, std::size_t DIM>
0026 Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::AxisAlignedBoundingBox(
0027     const entity_t* entity, const VertexType& center, const Size& size)
0028     : m_entity(entity),
0029       m_vmin(center - size.get() * 0.5),
0030       m_vmax(center + size.get() * 0.5),
0031       m_center(center),
0032       m_width(size.get()),
0033       m_iwidth(1 / m_width) {}
0034 
0035 template <typename entity_t, typename value_t, std::size_t DIM>
0036 Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::AxisAlignedBoundingBox(
0037     const std::vector<self_t*>& boxes, vertex_array_type envelope)
0038     : m_entity(nullptr) {
0039   assert(boxes.size() > 1);
0040 
0041   for (std::size_t i = 0; i < boxes.size(); i++) {
0042     if (i < boxes.size() - 1) {
0043       // set next on i to i+1
0044       boxes[i]->setSkip(boxes[i + 1]);
0045     } else {
0046       // make sure last is set to nullptr, this marks end
0047       // boxes[i]->m_next = nullptr;
0048       boxes[i]->setSkip(nullptr);
0049     }
0050   }
0051 
0052   m_leftChild = boxes.front();
0053   m_rightChild = boxes.back();
0054   m_skip = nullptr;
0055 
0056   std::tie(m_vmin, m_vmax) = wrap(boxes, envelope);
0057 
0058   m_center = (m_vmin + m_vmax) / 2.;
0059   m_width = m_vmax - m_vmin;
0060   m_iwidth = 1 / m_width;
0061 }
0062 
0063 template <typename entity_t, typename value_t, std::size_t DIM>
0064 std::pair<
0065     typename Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::VertexType,
0066     typename Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::VertexType>
0067 Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::wrap(
0068     const std::vector<const self_t*>& boxes, vertex_array_type envelope) {
0069   assert(boxes.size() > 1);
0070   // figure out extent of boxes
0071   // use array for Eigen coefficient wise min/max
0072   vertex_array_type vmax(
0073       vertex_array_type::Constant(std::numeric_limits<value_type>::lowest()));
0074   vertex_array_type vmin(
0075       vertex_array_type::Constant(std::numeric_limits<value_type>::max()));
0076 
0077   for (std::size_t i = 0; i < boxes.size(); i++) {
0078     vmin = vmin.min(boxes[i]->min().array());
0079     vmax = vmax.max(boxes[i]->max().array());
0080   }
0081 
0082   vmax += envelope;
0083   vmin -= envelope;
0084 
0085   return {vmin, vmax};
0086 }
0087 
0088 template <typename entity_t, typename value_t, std::size_t DIM>
0089 std::pair<
0090     typename Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::VertexType,
0091     typename Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::VertexType>
0092 Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::wrap(
0093     const std::vector<self_t*>& boxes, vertex_array_type envelope) {
0094   assert(boxes.size() > 1);
0095   std::vector<const self_t*> box_ptrs;
0096   box_ptrs.reserve(boxes.size());
0097   std::ranges::transform(boxes, std::back_inserter(box_ptrs),
0098                          [](const auto* box) { return box; });
0099   return wrap(box_ptrs, envelope);
0100 }
0101 
0102 template <typename entity_t, typename value_t, std::size_t DIM>
0103 std::pair<
0104     typename Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::VertexType,
0105     typename Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::VertexType>
0106 Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::wrap(
0107     const std::vector<self_t>& boxes, vertex_array_type envelope) {
0108   assert(boxes.size() > 1);
0109   std::vector<const self_t*> box_ptrs;
0110   box_ptrs.reserve(boxes.size());
0111   std::ranges::transform(boxes, std::back_inserter(box_ptrs),
0112                          [](auto& box) { return &box; });
0113   return wrap(box_ptrs, envelope);
0114 }
0115 
0116 template <typename entity_t, typename value_t, std::size_t DIM>
0117 bool Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::intersect(
0118     const VertexType& point) const {
0119   vertex_array_type t = (point - m_vmin).array() * m_iwidth;
0120   return t.minCoeff() >= 0 && t.maxCoeff() < 1;
0121 }
0122 
0123 template <typename entity_t, typename value_t, std::size_t DIM>
0124 bool Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::intersect(
0125     const Ray<value_type, DIM>& ray) const {
0126   const VertexType& origin = ray.origin();
0127   const vertex_array_type& idir = ray.idir();
0128 
0129   // Calculate the intersect distances with the min and max planes along the ray
0130   // direction, from the ray origin. See Ch VII.5 Fig.1 in [1].
0131   // This is done in all dimensions at the same time:
0132   vertex_array_type t0s = (m_vmin - origin).array() * idir;
0133   vertex_array_type t1s = (m_vmax - origin).array() * idir;
0134 
0135   // Calculate the component wise min/max between the t0s and t1s
0136   // this is non-compliant with IEEE-754-2008, NaN gets propagated through
0137   // https://eigen.tuxfamily.org/bz/show_bug.cgi?id=564
0138   // this means that rays parallel to boundaries might not be considered
0139   // to intersect.
0140   vertex_array_type tsmaller = t0s.min(t1s);
0141   vertex_array_type tbigger = t0s.max(t1s);
0142 
0143   // extract largest and smallest component of the component wise extrema
0144   value_type tmin = tsmaller.maxCoeff();
0145   value_type tmax = tbigger.minCoeff();
0146 
0147   // If tmin is smaller than tmax and tmax is positive, then the box is in
0148   // positive ray direction, and the ray intersects the box.
0149   return tmin < tmax && tmax > 0.0;
0150 }
0151 
0152 template <typename entity_t, typename value_t, std::size_t DIM>
0153 template <std::size_t sides>
0154 bool Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::intersect(
0155     const Frustum<value_type, DIM, sides>& fr) const {
0156   const auto& normals = fr.normals();
0157   // Transform vmin and vmax into the coordinate system, at which the frustum is
0158   // located at the coordinate origin.
0159   const vertex_array_type fr_vmin = m_vmin - fr.origin();
0160   const vertex_array_type fr_vmax = m_vmax - fr.origin();
0161 
0162   // For each plane, find the p-vertex, which is the vertex that is at the
0163   // furthest distance from the plane *along* its normal direction.
0164   // See Fig. 2 in [2].
0165   VertexType p_vtx;
0166   // for loop, we could eliminate this, probably,
0167   // but sides+1 is known at compile time, so the compiler
0168   // will most likely unroll the loop
0169   for (std::size_t i = 0; i < sides + 1; i++) {
0170     const VertexType& normal = normals[i];
0171 
0172     // for AABBs, take the component from the min vertex, if the normal
0173     // component is negative, else take the component from the max vertex.
0174     p_vtx = (normal.array() < 0).template cast<value_type>() * fr_vmin +
0175             (normal.array() >= 0).template cast<value_type>() * fr_vmax;
0176 
0177     // Check if the p-vertex is at positive or negative direction along the
0178     // If the p vertex is along negative normal direction *once*, the box is
0179     // outside the frustum, and we can terminate early.
0180     if (p_vtx.dot(normal) < 0) {
0181       // p vertex is outside on this plane, box must be outside
0182       return false;
0183     }
0184   }
0185 
0186   // If we get here, no p-vertex was outside, so box intersects or is
0187   // contained. We don't care, so report 'intersect'
0188   return true;
0189 }
0190 
0191 template <typename entity_t, typename value_t, std::size_t DIM>
0192 void Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::setSkip(
0193     self_t* skip) {
0194   // set next on this
0195   m_skip = skip;
0196   // find last child and set its skip
0197   if (m_rightChild != nullptr) {
0198     m_rightChild->setSkip(skip);
0199   }
0200 }
0201 
0202 template <typename entity_t, typename value_t, std::size_t DIM>
0203 const Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>*
0204 Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::getLeftChild() const {
0205   return m_leftChild;
0206 }
0207 
0208 template <typename entity_t, typename value_t, std::size_t DIM>
0209 const Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>*
0210 Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::getSkip() const {
0211   return m_skip;
0212 }
0213 
0214 template <typename entity_t, typename value_t, std::size_t DIM>
0215 bool Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::hasEntity() const {
0216   return m_entity != nullptr;
0217 }
0218 
0219 template <typename entity_t, typename value_t, std::size_t DIM>
0220 const entity_t* Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::entity()
0221     const {
0222   return m_entity;
0223 }
0224 
0225 template <typename entity_t, typename value_t, std::size_t DIM>
0226 void Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::setEntity(
0227     const entity_t* entity) {
0228   m_entity = entity;
0229 }
0230 
0231 template <typename entity_t, typename value_t, std::size_t DIM>
0232 const typename Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::VertexType&
0233 Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::center() const {
0234   return m_center;
0235 }
0236 
0237 template <typename entity_t, typename value_t, std::size_t DIM>
0238 const typename Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::VertexType&
0239 Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::min() const {
0240   return m_vmin;
0241 }
0242 
0243 template <typename entity_t, typename value_t, std::size_t DIM>
0244 const typename Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::VertexType&
0245 Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::max() const {
0246   return m_vmax;
0247 }
0248 
0249 template <typename entity_t, typename value_t, std::size_t DIM>
0250 std::ostream& Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::toStream(
0251     std::ostream& os) const {
0252   os << "AABB(ctr=(";
0253 
0254   for (std::size_t i = 0; i < DIM; i++) {
0255     if (i > 0) {
0256       os << ", ";
0257     }
0258     os << m_center[i];
0259   }
0260 
0261   os << ") vmin=(";
0262   for (std::size_t i = 0; i < DIM; i++) {
0263     if (i > 0) {
0264       os << ", ";
0265     }
0266     os << m_vmin[i];
0267   }
0268 
0269   os << ") vmax=(";
0270 
0271   for (std::size_t i = 0; i < DIM; i++) {
0272     if (i > 0) {
0273       os << ", ";
0274     }
0275     os << m_vmax[i];
0276   }
0277 
0278   os << "))";
0279 
0280   return os;
0281 }
0282 
0283 template <typename entity_t, typename value_t, std::size_t DIM>
0284 std::pair<
0285     typename Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::VertexType,
0286     typename Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::VertexType>
0287 Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::transformVertices(
0288     const transform_type& trf) const
0289   requires(DIM == 3)
0290 {
0291   // we need to enumerate all the vertices, transform,
0292   // and then recalculate min and max
0293 
0294   std::array<VertexType, 8> vertices({{
0295       {m_vmin.x(), m_vmin.y(), m_vmin.z()},
0296       {m_vmin.x(), m_vmax.y(), m_vmin.z()},
0297       {m_vmax.x(), m_vmax.y(), m_vmin.z()},
0298       {m_vmax.x(), m_vmin.y(), m_vmin.z()},
0299       {m_vmin.x(), m_vmin.y(), m_vmax.z()},
0300       {m_vmin.x(), m_vmax.y(), m_vmax.z()},
0301       {m_vmax.x(), m_vmax.y(), m_vmax.z()},
0302       {m_vmax.x(), m_vmin.y(), m_vmax.z()},
0303   }});
0304 
0305   VertexType vmin = trf * vertices[0];
0306   VertexType vmax = trf * vertices[0];
0307 
0308   for (std::size_t i = 1; i < 8; i++) {
0309     const VertexType vtx = trf * vertices[i];
0310     vmin = vmin.cwiseMin(vtx);
0311     vmax = vmax.cwiseMax(vtx);
0312   }
0313 
0314   return {vmin, vmax};
0315 }
0316 
0317 template <typename entity_t, typename value_t, std::size_t DIM>
0318 std::pair<
0319     typename Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::VertexType,
0320     typename Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::VertexType>
0321 Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::transformVertices(
0322     const transform_type& trf) const
0323   requires(DIM == 2)
0324 {
0325   // we need to enumerate all the vertices, transform,
0326   // and then recalculate min and max
0327 
0328   std::array<VertexType, 4> vertices({{{m_vmin.x(), m_vmin.y()},
0329                                        {m_vmin.x(), m_vmax.y()},
0330                                        {m_vmax.x(), m_vmax.y()},
0331                                        {m_vmax.x(), m_vmin.y()}}});
0332 
0333   VertexType vmin = trf * vertices[0];
0334   VertexType vmax = trf * vertices[0];
0335 
0336   for (std::size_t i = 1; i < 4; i++) {
0337     const VertexType vtx = trf * vertices[i];
0338     vmin = vmin.cwiseMin(vtx);
0339     vmax = vmax.cwiseMax(vtx);
0340   }
0341 
0342   return {vmin, vmax};
0343 }
0344 
0345 template <typename entity_t, typename value_t, std::size_t DIM>
0346 void Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::transform(
0347     const transform_type& trf) {
0348   std::tie(m_vmin, m_vmax) = transformVertices(trf);
0349 }
0350 
0351 template <typename entity_t, typename value_t, std::size_t DIM>
0352 Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>
0353 Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::transformed(
0354     const transform_type& trf) const {
0355   const auto [vmin, vmax] = transformVertices(trf);
0356   return self_t(m_entity, vmin, vmax);
0357 }
0358 
0359 template <typename entity_t, typename value_t, std::size_t DIM>
0360 void Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::draw(
0361     IVisualization3D& helper, Color color, const transform_type& trf) const
0362   requires(DIM == 3)
0363 {
0364   static_assert(DIM == 3, "PLY output only supported in 3D");
0365 
0366   const VertexType& vmin = m_vmin;
0367   const VertexType& vmax = m_vmax;
0368 
0369   auto write = [&](const VertexType& a, const VertexType& b,
0370                    const VertexType& c, const VertexType& d) {
0371     helper.face(std::vector<VertexType>({trf * a, trf * b, trf * c, trf * d}),
0372                 color);
0373   };
0374 
0375   write({vmin.x(), vmin.y(), vmin.z()}, {vmin.x(), vmax.y(), vmin.z()},
0376         {vmin.x(), vmax.y(), vmax.z()}, {vmin.x(), vmin.y(), vmax.z()});
0377 
0378   write({vmax.x(), vmin.y(), vmin.z()}, {vmax.x(), vmax.y(), vmin.z()},
0379         {vmax.x(), vmax.y(), vmax.z()}, {vmax.x(), vmin.y(), vmax.z()});
0380 
0381   write({vmin.x(), vmin.y(), vmin.z()}, {vmax.x(), vmin.y(), vmin.z()},
0382         {vmax.x(), vmin.y(), vmax.z()}, {vmin.x(), vmin.y(), vmax.z()});
0383 
0384   write({vmin.x(), vmax.y(), vmin.z()}, {vmax.x(), vmax.y(), vmin.z()},
0385         {vmax.x(), vmax.y(), vmax.z()}, {vmin.x(), vmax.y(), vmax.z()});
0386 
0387   write({vmin.x(), vmin.y(), vmin.z()}, {vmax.x(), vmin.y(), vmin.z()},
0388         {vmax.x(), vmax.y(), vmin.z()}, {vmin.x(), vmax.y(), vmin.z()});
0389 
0390   write({vmin.x(), vmin.y(), vmax.z()}, {vmax.x(), vmin.y(), vmax.z()},
0391         {vmax.x(), vmax.y(), vmax.z()}, {vmin.x(), vmax.y(), vmax.z()});
0392 }
0393 
0394 template <typename entity_t, typename value_t, std::size_t DIM>
0395 std::ostream& Acts::AxisAlignedBoundingBox<entity_t, value_t, DIM>::svg(
0396     std::ostream& os, value_type w, value_type h, value_type unit,
0397     const std::string& label, const std::string& fillcolor) const
0398   requires(DIM == 2)
0399 {
0400   static_assert(DIM == 2, "SVG is only supported in 2D");
0401 
0402   VertexType mid(w / 2., h / 2.);
0403 
0404   using transform_t = Eigen::Transform<value_t, DIM, Eigen::Affine>;
0405 
0406   transform_t trf = transform_t::Identity();
0407   trf.translate(mid);
0408   trf = trf * Eigen::Scaling(VertexType(1, -1));
0409   trf.scale(unit);
0410 
0411   auto draw_point = [&](const VertexType& p_, const std::string& color,
0412                         std::size_t r) {
0413     VertexType p = trf * p_;
0414     os << "<circle ";
0415     os << "cx=\"" << p.x() << "\" cy=\"" << p.y() << "\" r=\"" << r << "\"";
0416     os << " fill=\"" << color << "\"";
0417     os << "/>\n";
0418   };
0419 
0420   auto draw_rect = [&](const VertexType& center_, const VertexType& size_,
0421                        const std::string& color) {
0422     VertexType size = size_ * unit;
0423     VertexType center = trf * center_ - size * 0.5;
0424 
0425     os << "<rect ";
0426     os << "x=\"" << center.x() << "\" y=\"" << center.y() << "\" ";
0427     os << "width=\"" << size.x() << "\" height=\"" << size.y() << "\"";
0428     os << " fill=\"" << color << "\"";
0429     os << "/>\n";
0430   };
0431 
0432   auto draw_text = [&](const VertexType& center_, const std::string& text,
0433                        const std::string& color, std::size_t size) {
0434     VertexType center = trf * center_;
0435     os << "<text dominant-baseline=\"middle\" text-anchor=\"middle\" ";
0436     os << "fill=\"" << color << "\" font-size=\"" << size << "\" ";
0437     os << "x=\"" << center.x() << "\" y=\"" << center.y() << "\">";
0438     os << text << "</text>\n";
0439   };
0440 
0441   draw_rect(m_center, m_width, fillcolor);
0442   draw_point(m_vmin, "black", 2);
0443   draw_point(m_vmax, "black", 2);
0444   draw_text(m_center, label, "white", 10);
0445 
0446   return os;
0447 }
0448 
0449 template <typename object_t, typename box_t, typename visitor_t>
0450 void Acts::BoundingBoxHierarchy::visitIntersecting(const object_t& object,
0451                                                    const box_t* root,
0452                                                    visitor_t&& visitor) {
0453   const box_t* node = root;
0454   while (node != nullptr) {
0455     if (!node->intersect(object)) {
0456       node = node->getSkip();
0457       continue;
0458     }
0459 
0460     if (!node->hasEntity()) {
0461       node = node->getLeftChild();
0462       continue;
0463     }
0464 
0465     std::invoke(std::forward<visitor_t>(visitor), *node->entity());
0466     node = node->getSkip();
0467   }
0468 }
0469 
0470 namespace Acts::BoundingBoxHierarchy::detail {
0471 
0472 template <typename box_t>
0473 box_t* octreeInner(std::vector<std::unique_ptr<box_t>>& store,
0474                    std::size_t maxDepth,
0475                    typename box_t::vertex_array_type envelope,
0476                    const std::vector<box_t*>& localPrimitives,
0477                    std::size_t depth) {
0478   using VertexType = typename box_t::VertexType;
0479 
0480   assert(!localPrimitives.empty());
0481   if (localPrimitives.size() == 1) {
0482     // just return
0483     return localPrimitives.front();
0484   }
0485 
0486   if (depth >= maxDepth) {
0487     // just wrap them all up
0488     auto bb = std::make_unique<box_t>(localPrimitives, envelope);
0489     store.push_back(std::move(bb));
0490     return store.back().get();
0491   }
0492 
0493   std::array<std::vector<box_t*>, 8> octants;
0494   // calc center of boxes
0495   const auto [vmin, vmax] = box_t::wrap(localPrimitives);
0496   VertexType globalCenter = (vmin + vmax) / 2.;
0497 
0498   for (auto* box : localPrimitives) {
0499     VertexType center = box->center() - globalCenter;
0500     if (center.x() < 0 && center.y() < 0 && center.z() < 0) {
0501       octants[0].push_back(box);
0502       continue;
0503     }
0504     if (center.x() > 0 && center.y() < 0 && center.z() < 0) {
0505       octants[1].push_back(box);
0506       continue;
0507     }
0508     if (center.x() < 0 && center.y() > 0 && center.z() < 0) {
0509       octants[2].push_back(box);
0510       continue;
0511     }
0512     if (center.x() > 0 && center.y() > 0 && center.z() < 0) {
0513       octants[3].push_back(box);
0514       continue;
0515     }
0516 
0517     if (center.x() < 0 && center.y() < 0 && center.z() > 0) {
0518       octants[4].push_back(box);
0519       continue;
0520     }
0521     if (center.x() > 0 && center.y() < 0 && center.z() > 0) {
0522       octants[5].push_back(box);
0523       continue;
0524     }
0525     if (center.x() < 0 && center.y() > 0 && center.z() > 0) {
0526       octants[6].push_back(box);
0527       continue;
0528     }
0529     if (center.x() > 0 && center.y() > 0 && center.z() > 0) {
0530       octants[7].push_back(box);
0531       continue;
0532     }
0533 
0534     // not in any quadrant (numerics probably)
0535     octants[0].push_back(box);
0536   }
0537 
0538   std::vector<box_t*> subOctants;
0539   for (const auto& subPrimitives : octants) {
0540     if (subPrimitives.size() <= 8) {
0541       if (subPrimitives.empty()) {
0542         // done
0543       } else if (subPrimitives.size() == 1) {
0544         subOctants.push_back(subPrimitives.front());
0545       } else {
0546         store.push_back(std::make_unique<box_t>(subPrimitives, envelope));
0547         subOctants.push_back(store.back().get());
0548       }
0549     } else {
0550       // recurse
0551       subOctants.push_back(
0552           octreeInner(store, maxDepth, envelope, subPrimitives, depth + 1));
0553     }
0554   }
0555 
0556   if (subOctants.size() == 1) {
0557     return subOctants.front();
0558   }
0559 
0560   auto bb = std::make_unique<box_t>(subOctants, envelope);
0561   store.push_back(std::move(bb));
0562   return store.back().get();
0563 }
0564 
0565 }  // namespace Acts::BoundingBoxHierarchy::detail
0566 
0567 template <typename box_t>
0568 box_t* Acts::BoundingBoxHierarchy::makeOctree(
0569     std::vector<std::unique_ptr<box_t>>& store,
0570     const std::vector<box_t*>& prims, std::size_t maxDepth,
0571     typename box_t::value_type envelopeValue) {
0572   static_assert(box_t::dim == 3, "Octree can only be created in 3D");
0573 
0574   using vertex_array_type = typename box_t::vertex_array_type;
0575 
0576   vertex_array_type envelope(vertex_array_type::Constant(envelopeValue));
0577 
0578   box_t* top = detail::octreeInner(store, maxDepth, envelope, prims, 0);
0579   return top;
0580 }
0581 
0582 template <typename box_t>
0583 box_t* Acts::make_octree(std::vector<std::unique_ptr<box_t>>& store,
0584                          const std::vector<box_t*>& prims, std::size_t maxDepth,
0585                          typename box_t::value_type envelope) {
0586   return BoundingBoxHierarchy::makeOctree(store, prims, maxDepth, envelope);
0587 }
0588 
0589 template <typename T, typename U, std::size_t V>
0590 std::ostream& Acts::operator<<(
0591     std::ostream& os, const Acts::AxisAlignedBoundingBox<T, U, V>& box) {
0592   return box.toStream(os);
0593 }