orange/detail/BIHTraverser.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2022-2024 UT-Battelle, LLC, and other Celeritas developers.
0003 // See the top-level COPYRIGHT file for details.
0004 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0005 //---------------------------------------------------------------------------//
0006 //! \file orange/detail/BIHTraverser.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/math/Algorithms.hh"
0011
0012 #include "../BoundingBoxUtils.hh"
0013 #include "../OrangeData.hh"
0014
0015 namespace celeritas
0016 {
0017 namespace detail
0018 {
0019 //---------------------------------------------------------------------------//
0020 /*!
0021  * Traverse BIH tree using a depth-first search.
0022  *
0023  * \todo move to top-level orange directory out of detail namespace
0024  */
0025 class BIHTraverser
0026 {
0027   public:
0028     //!@{
0029     //! \name Type aliases
0030     using Storage = NativeCRef<BIHTreeData>;
0031     //!@}
0032
0033     // Construct from vector of bounding boxes and storage for LocalVolumeIds
0034     inline CELER_FUNCTION
0035     BIHTraverser(BIHTree const& tree, Storage const& storage);
0036
0037     // Point-in-volume operation
0038     template<class F>
0039     inline CELER_FUNCTION LocalVolumeId operator()(Real3 const& point,
0040                                                    F&& is_inside) const;
0041
0042   private:
0043     //// DATA ////
0044     BIHTree const& tree_;
0045     Storage const& storage_;
0046     size_type leaf_offset_;
0047
0048     //// HELPER FUNCTIONS ////
0049
0050     // Get the ID of the next node in the traversal sequence
0051     inline CELER_FUNCTION BIHNodeId next_node(BIHNodeId const& current_id,
0052                                               BIHNodeId const& previous_id,
0053                                               Real3 const& point) const;
0054
0055     // Determine if traversal shall proceed down a given edge
0056     inline CELER_FUNCTION bool visit_edge(BIHInnerNode const& node,
0057                                           BIHInnerNode::Edge edge,
0058                                           Real3 const& point) const;
0059
0060     // Determine if a node is inner, i.e., not a leaf
0061     inline CELER_FUNCTION bool is_inner(BIHNodeId id) const;
0062
0063     // Get an inner node for a given BIHNodeId
0064     inline CELER_FUNCTION BIHInnerNode const&
0065     get_inner_node(BIHNodeId id) const;
0066
0067     // Get a leaf node for a given BIHNodeId
0068     inline CELER_FUNCTION BIHLeafNode const& get_leaf_node(BIHNodeId id) const;
0069
0070     // Determine if any leaf node volumes contain the point
0071     template<class F>
0072     inline CELER_FUNCTION LocalVolumeId visit_leaf(BIHLeafNode const& leaf_node,
0073                                                    Real3 const& point,
0074                                                    F&& is_inside) const;
0075
0076     // Determine if any inf_vols contain the point
0077     template<class F>
0078     inline CELER_FUNCTION LocalVolumeId visit_inf_vols(F&& is_inside) const;
0079
0080     // Determine if a single bbox contains the point
0081     inline CELER_FUNCTION bool
0082     visit_bbox(LocalVolumeId const& id, Real3 const& point) const;
0083 };
0084
0085 //---------------------------------------------------------------------------//
0086 // INLINE DEFINITIONS
0087 //---------------------------------------------------------------------------//
0088 /*!
0089  * Construct from vector of bounding boxes and storage.
0090  */
0091 CELER_FUNCTION
0092 BIHTraverser::BIHTraverser(BIHTree const& tree,
0093                            BIHTraverser::Storage const& storage)
0094     : tree_(tree), storage_(storage), leaf_offset_(tree.inner_nodes.size())
0095 {
0096     CELER_EXPECT(tree);
0097 }
0098
0099 //---------------------------------------------------------------------------//
0100 /*!
0101  * Point-in-volume operation.
0102  */
0103 template<class F>
0104 CELER_FUNCTION LocalVolumeId BIHTraverser::operator()(Real3 const& point,
0105                                                       F&& is_inside) const
0106 {
0107     BIHNodeId previous_node;
0108     BIHNodeId current_node{0};
0109     LocalVolumeId id;
0110
0111     do
0112     {
0113         if (!this->is_inner(current_node))
0114         {
0115             id = this->visit_leaf(
0116                 this->get_leaf_node(current_node), point, is_inside);
0117
0118             if (id)
0119             {
0120                 return id;
0121             }
0122         }
0123
0124         previous_node = exchange(
0125             current_node, this->next_node(current_node, previous_node, point));
0126
0127     } while (current_node);
0128
0129     if (!id)
0130     {
0131         id = this->visit_inf_vols(is_inside);
0132     }
0133
0134     return id;
0135 }
0136
0137 //---------------------------------------------------------------------------//
0138 // HELPER FUNCTIONS
0139 //---------------------------------------------------------------------------//
0140 /*!
0141  *  Get the ID of the next node in the traversal sequence.
0142  */
0143 CELER_FUNCTION
0144 BIHNodeId BIHTraverser::next_node(BIHNodeId const& current_id,
0145                                   BIHNodeId const& previous_id,
0146                                   Real3 const& point) const
0147 {
0148     using Edge = BIHInnerNode::Edge;
0149
0150     BIHNodeId next_id;
0151
0152     if (this->is_inner(current_id))
0153     {
0154         auto const& current_node = this->get_inner_node(current_id);
0155         if (previous_id == current_node.parent)
0156         {
0157             // Visiting this inner node for the first time; go down either left
0158             // or right edge
0159             if (this->visit_edge(current_node, Edge::left, point))
0160             {
0161                 next_id = current_node.bounding_planes[Edge::left].child;
0162             }
0163             else
0164             {
0165                 next_id = current_node.bounding_planes[Edge::right].child;
0166             }
0167         }
0168         else if (previous_id == current_node.bounding_planes[Edge::left].child)
0169         {
0170             // Visiting this inner node for the second time; go down right edge
0171             // or return to parent
0172             if (this->visit_edge(current_node, Edge::right, point))
0173             {
0174                 next_id = current_node.bounding_planes[Edge::right].child;
0175             }
0176             else
0177             {
0178                 next_id = current_node.parent;
0179             }
0180         }
0181         else
0182         {
0183             // Visiting this inner node for the third time; return to parent
0184             CELER_EXPECT(previous_id
0185                          == current_node.bounding_planes[Edge::right].child);
0186             next_id = current_node.parent;
0187         }
0188     }
0189     else
0190     {
0191         // Leaf node; return to parent
0192         CELER_EXPECT(previous_id == this->get_leaf_node(current_id).parent);
0193         next_id = previous_id;
0194     }
0195
0196     return next_id;
0197 }
0198
0199 //---------------------------------------------------------------------------//
0200 /*!
0201  * Determine if traversal shall proceed down a given edge.
0202  */
0203 CELER_FUNCTION
0204 bool BIHTraverser::visit_edge(BIHInnerNode const& node,
0205                               BIHInnerNode::Edge edge,
0206                               Real3 const& point) const
0207 {
0208     CELER_EXPECT(edge < BIHInnerNode::Edge::size_);
0209
0210     auto pos = node.bounding_planes[edge].position;
0211     auto point_pos = point[to_int(node.axis)];
0212
0213     return (edge == BIHInnerNode::Edge::left) ? (point_pos < pos)
0214                                               : (pos < point_pos);
0215 }
0216
0217 //---------------------------------------------------------------------------//
0218 /*!
0219  *  Determine if a node is inner, i.e., not a leaf.
0220  */
0221 CELER_FUNCTION
0222 bool BIHTraverser::is_inner(BIHNodeId id) const
0223 {
0224     return id.unchecked_get() < leaf_offset_;
0225 }
0226
0227 //---------------------------------------------------------------------------//
0228 /*!
0229  *  Get an inner node for a given BIHNodeId.
0230  */
0231 CELER_FUNCTION
0232 BIHInnerNode const& BIHTraverser::get_inner_node(BIHNodeId id) const
0233 {
0234     CELER_EXPECT(this->is_inner(id));
0235     return storage_.inner_nodes[tree_.inner_nodes[id.unchecked_get()]];
0236 }
0237
0238 //---------------------------------------------------------------------------//
0239 /*!
0240  *  Get a leaf node for a given BIHNodeId.
0241  */
0242 CELER_FUNCTION
0243 BIHLeafNode const& BIHTraverser::get_leaf_node(BIHNodeId id) const
0244 {
0245     CELER_EXPECT(!this->is_inner(id));
0246     return storage_
0247         .leaf_nodes[tree_.leaf_nodes[id.unchecked_get() - leaf_offset_]];
0248 }
0249
0250 //---------------------------------------------------------------------------//
0251 /*!
0252  * Determine if any leaf node volumes contain the point.
0253  */
0254 template<class F>
0255 CELER_FUNCTION LocalVolumeId BIHTraverser::visit_leaf(
0256     BIHLeafNode const& leaf_node, Real3 const& point, F&& is_inside) const
0257 {
0258     for (auto i : range(leaf_node.vol_ids.size()))
0259     {
0260         auto id = storage_.local_volume_ids[leaf_node.vol_ids[i]];
0261         if (this->visit_bbox(id, point) && is_inside(id))
0262         {
0263             return id;
0264         }
0265     }
0266     return LocalVolumeId{};
0267 }
0268
0269 //---------------------------------------------------------------------------//
0270 /*!
0271  * Determine if any volumes in inf_vols contain the point.
0272  */
0273 template<class F>
0274 CELER_FUNCTION LocalVolumeId BIHTraverser::visit_inf_vols(F&& is_inside) const
0275 {
0276     for (auto i : range(tree_.inf_volids.size()))
0277     {
0278         auto id = storage_.local_volume_ids[tree_.inf_volids[i]];
0279         if (is_inside(id))
0280         {
0281             return id;
0282         }
0283     }
0284     return LocalVolumeId{};
0285 }
0286
0287 //---------------------------------------------------------------------------//
0288 /*!
0289  * Determinate if a single bbox contains the point.
0290  */
0291 CELER_FUNCTION
0292 bool BIHTraverser::visit_bbox(LocalVolumeId const& id, Real3 const& point) const
0293 {
0294     return is_inside(storage_.bboxes[tree_.bboxes[id]], point);
0295 }
0296
0297 //---------------------------------------------------------------------------//
0298 }  // namespace detail
0299 }  // namespace celeritas