perf(AABB): parallel and deterministic self intersections

.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,

include/geode/geometry/aabb.hpp

@@ -29,6 +29,9 @@
 
 #pragma once
 
+#include <utility>
+#include <vector>
+
 #include <absl/types/span.h>
 
 #include <geode/basic/pimpl.hpp>
@@ -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

include/geode/geometry/detail/aabb_impl.hpp

@@ -31,8 +31,13 @@
 
 #include <algorithm>
 #include <cmath>
+#include <utility>
+#include <vector>
 
 #include <absl/algorithm/container.h>
+#include <absl/container/fixed_array.h>
+
+#include <async++.h>
 
 #include <geode/basic/pimpl_impl.hpp>
 
@@ -71,6 +76,7 @@
     {
     public:
         static constexpr index_t ROOT_INDEX{ 0 };
+        static constexpr index_t SELF_INTERSECTION_CHUNK_SIZE{ 128 };
 
         struct Iterator
         {
@@ -91,13 +97,13 @@
     public:
         Impl() = default;
 
         Impl( absl::Span< const BoundingBox< dimension > > bboxes )
             : element_order_( bboxes.size() )
         {
             if( !bboxes.empty() )
             {
                 absl::c_iota( element_order_, index_t{ 0 } );
                 tree_.resize( 2 * bboxes.size() - 1 );
                 index_t next_free_index{ 0 };
                 initialize_tree_recursive(
                     bboxes, next_free_index, 0, bboxes.size() );
@@ -115,13 +121,13 @@
             return element_begin + 1 == element_end;
         }
 
         [[nodiscard]] Iterator get_recursive_iterators( index_t node_index,
             index_t element_begin,
             index_t element_end ) const
         {
             Iterator it;
             it.element_middle =
                 element_begin + ( element_end - element_begin ) / 2;
             it.child_left = node_index + 1;
             it.child_right = tree_[node_index].right_child;
             return it;
@@ -164,7 +170,7 @@
             }
             const auto axis = std::get< 0 >( range_box.largest_length() );
             const auto element_middle =
                 element_begin + ( element_end - element_begin ) / 2;
             std::nth_element( element_order_.begin() + element_begin,
                 element_order_.begin() + element_middle,
                 element_order_.begin() + element_end,
@@ -184,7 +190,7 @@
         }
 
         template < typename ACTION >
         void closest_element_box_recursive( const Point< dimension >& query,
             index_t& nearest_box,
             double& squared_distance,
             index_t node_index,
@@ -212,7 +218,7 @@
                 }
                 return;
             }
             const auto it = get_recursive_iterators(
                 node_index, element_begin, element_end );
             const auto squared_distance_left =
                 node( it.child_left ).squared_signed_distance( query );
@@ -254,77 +260,100 @@
         }
 
         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 @@
 
     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 >

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 >