diff --git a/.clang-tidy b/.clang-tidy index 725ab31e1..5dacfaeed 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -14,6 +14,7 @@ Checks: > -misc-no-recursion, -misc-include-cleaner, -misc-const-correctness, + -modernize-use-scoped-lock, -modernize-use-trailing-return-type, -portability-avoid-pragma-once, -readability-use-anyofallof, diff --git a/include/geode/geometry/aabb.hpp b/include/geode/geometry/aabb.hpp index ff87029f0..567d40cd1 100644 --- a/include/geode/geometry/aabb.hpp +++ b/include/geode/geometry/aabb.hpp @@ -29,6 +29,9 @@ #pragma once +#include +#include + #include #include @@ -125,24 +128,26 @@ namespace geode EvalIntersection& action ) const; /*! - * @brief Computes the self intersections of the element boxes. + * @brief Computes the self intersections of the element boxes, in + * parallel. * @param[in] action The functor to run when two boxes intersect * @tparam EvalIntersection this functor should have an operator() * defined like this: * bool operator()( index_t cur_element_box1, index_t cur_element_box2 - * ); + * ) const; * @note cur_element_box1 and cur_element_box2 are the element box - * indices that intersect. - * @note the operator defines what to do when two boxes of the - * tree ( \p cur_element_box1 and \p cur_element_box2 ) intersect each - * other (for example: test real intersection between each element in - * boxes and store the result.) - * @note The returned boolean indicates if the search should stop or - * continue. Return true to stop the search, false to continue. + * indices that intersect. Each pair is evaluated once. + * @note the operator is called concurrently from several threads, it + * must be thread-safe. It tests the two elements (for example: real + * intersection between the elements in the boxes) and returns true if + * the pair should be kept. + * @return The kept pairs. Their order is deterministic: it does not + * depend on the number of threads. */ template < class EvalIntersection > - void compute_self_element_bbox_intersections( - EvalIntersection& action ) const; + [[nodiscard]] std::vector< std::pair< index_t, index_t > > + compute_self_element_bbox_intersections( + const EvalIntersection& action ) const; /*! * @brief Computes all the intersections of the element boxes between diff --git a/include/geode/geometry/detail/aabb_impl.hpp b/include/geode/geometry/detail/aabb_impl.hpp index 413ea861a..15e16c1b6 100644 --- a/include/geode/geometry/detail/aabb_impl.hpp +++ b/include/geode/geometry/detail/aabb_impl.hpp @@ -31,8 +31,13 @@ #include #include +#include +#include #include +#include + +#include #include @@ -71,6 +76,7 @@ namespace geode { public: static constexpr index_t ROOT_INDEX{ 0 }; + static constexpr index_t SELF_INTERSECTION_CHUNK_SIZE{ 128 }; struct Iterator { @@ -254,77 +260,100 @@ namespace geode } template < typename ACTION > - bool self_intersect_recursive( index_t node_index1, - index_t element_begin1, - index_t element_end1, - index_t node_index2, - index_t element_begin2, - index_t element_end2, - ACTION& action ) const + [[nodiscard]] std::vector< std::pair< index_t, index_t > > + self_intersect( const ACTION& action ) const { - OpenGeodeGeometryException::check_assertion( - element_end1 != element_begin1, - "No iteration allowed start == end" ); - OpenGeodeGeometryException::check_assertion( - element_end2 != element_begin2, - "No iteration allowed start == end" ); - - // Since we are intersecting the AABBTree with *itself*, - // we can prune half of the cases by skipping the test - // whenever node2's polygon index interval is greated than - // node1's polygon index interval. - if( element_end2 <= element_begin1 ) + const auto nb_chunks = + ( nb_bboxes() + SELF_INTERSECTION_CHUNK_SIZE - 1 ) + / SELF_INTERSECTION_CHUNK_SIZE; + absl::FixedArray< std::vector< std::pair< index_t, index_t > > > + chunk_pairs( nb_chunks ); + async::parallel_for( async::irange( index_t{ 0 }, nb_chunks ), + [this, &action, &chunk_pairs]( index_t chunk ) { + const auto chunk_begin = + chunk * SELF_INTERSECTION_CHUNK_SIZE; + const auto chunk_end = + std::min( chunk_begin + SELF_INTERSECTION_CHUNK_SIZE, + nb_bboxes() ); + for( const auto position : Range{ chunk_begin, chunk_end } ) + { + leaf_self_intersect_recursive( leaf_node( position ), + position, ROOT_INDEX, 0, nb_bboxes(), action, + chunk_pairs[chunk] ); + } + } ); + std::size_t nb_pairs{ 0 }; + for( const auto& pairs : chunk_pairs ) { - return false; + nb_pairs += pairs.size(); } - - // The acceleration is here: - if( !node( node_index1 ).epsilon_intersects( node( node_index2 ) ) ) + std::vector< std::pair< index_t, index_t > > result; + result.reserve( nb_pairs ); + for( const auto& pairs : chunk_pairs ) { - return false; + result.insert( result.end(), pairs.begin(), pairs.end() ); } + return result; + } - // Simple case: leaf - leaf intersection. - if( is_leaf( element_begin1, element_end1 ) - && is_leaf( element_begin2, element_end2 ) ) + [[nodiscard]] index_t leaf_node( index_t position ) const + { + index_t element_begin{ 0 }; + index_t element_end{ nb_bboxes() }; + index_t node_index{ ROOT_INDEX }; + while( !is_leaf( element_begin, element_end ) ) { - if( node_index1 == node_index2 ) + const auto it = get_recursive_iterators( + node_index, element_begin, element_end ); + if( position < it.element_middle ) { - return false; + element_end = it.element_middle; + node_index = it.child_left; + } + else + { + element_begin = it.element_middle; + node_index = it.child_right; } - return action( element_order( element_begin1 ), - element_order( element_begin2 ) ); } + return node_index; + } - // If node2 has more polygons than node1, then - // intersect node2's two children with node1 - // else - // intersect node1's two children with node2 - if( element_end2 - element_begin2 > element_end1 - element_begin1 ) + template < typename ACTION > + void leaf_self_intersect_recursive( index_t leaf_node_index, + index_t leaf_position, + index_t node_index, + index_t element_begin, + index_t element_end, + const ACTION& action, + std::vector< std::pair< index_t, index_t > >& pairs ) const + { + if( element_end <= leaf_position + 1 ) { - const auto it = get_recursive_iterators( - node_index2, element_begin2, element_end2 ); - if( self_intersect_recursive( node_index1, element_begin1, - element_end1, it.child_left, element_begin2, - it.element_middle, action ) ) + return; + } + if( !node( leaf_node_index ) + .epsilon_intersects( node( node_index ) ) ) + { + return; + } + if( is_leaf( element_begin, element_end ) ) + { + const auto element = element_order( leaf_position ); + const auto other_element = element_order( element_begin ); + if( action( element, other_element ) ) { - return true; + pairs.emplace_back( element, other_element ); } - return self_intersect_recursive( node_index1, element_begin1, - element_end1, it.child_right, it.element_middle, - element_end2, action ); + return; } const auto it = get_recursive_iterators( - node_index1, element_begin1, element_end1 ); - if( self_intersect_recursive( it.child_left, element_begin1, - it.element_middle, node_index2, element_begin2, - element_end2, action ) ) - { - return true; - } - return self_intersect_recursive( it.child_right, it.element_middle, - element_end1, node_index2, element_begin2, element_end2, - action ); + node_index, element_begin, element_end ); + leaf_self_intersect_recursive( leaf_node_index, leaf_position, + it.child_left, element_begin, it.element_middle, action, + pairs ); + leaf_self_intersect_recursive( leaf_node_index, leaf_position, + it.child_right, it.element_middle, element_end, action, pairs ); } template < typename ACTION > @@ -532,15 +561,15 @@ namespace geode template < index_t dimension > template < class EvalIntersection > - void AABBTree< dimension >::compute_self_element_bbox_intersections( - EvalIntersection& action ) const + std::vector< std::pair< index_t, index_t > > + AABBTree< dimension >::compute_self_element_bbox_intersections( + const EvalIntersection& action ) const { if( nb_bboxes() == 0 ) { - return; + return {}; } - impl_->self_intersect_recursive( Impl::ROOT_INDEX, 0, nb_bboxes(), - Impl::ROOT_INDEX, 0, nb_bboxes(), action ); + return impl_->self_intersect( action ); } template < index_t dimension > diff --git a/tests/geometry/test-aabb.cpp b/tests/geometry/test-aabb.cpp index e7bbb8fcd..c7a6de68c 100644 --- a/tests/geometry/test-aabb.cpp +++ b/tests/geometry/test-aabb.cpp @@ -181,32 +181,37 @@ class BoxAABBIntersection return false; } - // test box strict inclusion - bool box_contains_box( geode::index_t box1, geode::index_t box2 ) +public: + std::mutex mutex_; + absl::flat_hash_set< geode::index_t > box_intersections_; + +private: + absl::Span< const geode::BoundingBox< dimension > > bounding_boxes_; +}; + +template < geode::index_t dimension > +class BoxAABBInclusion +{ +public: + explicit BoxAABBInclusion( + absl::Span< const geode::BoundingBox< dimension > > bounding_boxes ) + : bounding_boxes_( bounding_boxes ) + { + } + + [[nodiscard]] bool box_contains_box( + geode::index_t box1, geode::index_t box2 ) const { return bounding_boxes_[box1].contains( bounding_boxes_[box2].min() ) && bounding_boxes_[box1].contains( bounding_boxes_[box2].max() ); } - bool operator()( geode::index_t box1, geode::index_t box2 ) + + [[nodiscard]] bool operator()( + geode::index_t box1, geode::index_t box2 ) const { - if( box_contains_box( box1, box2 ) ) - { - std::lock_guard< std::mutex > lock( mutex_ ); - included_box_.emplace_back( box1, box2 ); - } - else if( box_contains_box( box2, box1 ) ) - { - std::lock_guard< std::mutex > lock( mutex_ ); - included_box_.emplace_back( box2, box1 ); - } - return false; + return box_contains_box( box1, box2 ) || box_contains_box( box2, box1 ); } -public: - std::mutex mutex_; - absl::flat_hash_set< geode::index_t > box_intersections_; - std::vector< std::pair< geode::index_t, geode::index_t > > included_box_; - private: absl::Span< const geode::BoundingBox< dimension > > bounding_boxes_; }; @@ -414,22 +419,28 @@ void test_self_intersections() geode::AABBTree< dimension > aabb{ box_vector }; - BoxAABBIntersection< dimension > eval_intersection{ box_vector }; - // investigate box inclusions - eval_intersection.included_box_.clear(); - aabb.compute_self_element_bbox_intersections( eval_intersection ); + const BoxAABBInclusion< dimension > box_inclusion{ box_vector }; + const auto included_boxes = + aabb.compute_self_element_bbox_intersections( box_inclusion ); geode::OpenGeodeGeometryException::test( - eval_intersection.included_box_.size() == nb_boxes * nb_boxes, + included_boxes.size() == nb_boxes * nb_boxes, "Box self intersection - Every box should have one box " "inside" ); - for( const auto& result : eval_intersection.included_box_ ) + for( const auto& [box1, box2] : included_boxes ) { + const auto container = + box_inclusion.box_contains_box( box1, box2 ) ? box1 : box2; + const auto contained = container == box1 ? box2 : box1; geode::OpenGeodeGeometryException::test( - result.first == result.second - ( nb_boxes * nb_boxes ), + container == contained - ( nb_boxes * nb_boxes ), "Box self intersection - Wrong box inclusion result" ); } + geode::OpenGeodeGeometryException::test( + aabb.compute_self_element_bbox_intersections( box_inclusion ) + == included_boxes, + "Box self intersection - Result order should be deterministic" ); } template < geode::index_t dimension >