Improve SHAPE_POLY_SET fracture performance

Refactors `SHAPE_POLY_SET::fractureSingle()` to be more efficient, while
not changing the actual algorithm:

* increase cache locality by using contiguous arrays instead of what was
effectively a linked list
* reduce latency and jitter by replacing per-edge allocator calls with
ahead-of-time std::vector reserves
* increase cache efficiency by making the vertex struct smaller
* replace O(n^2) leftmost edge search with O(n log n) std::sort
* sort the polygons instead of the edges
* cut iteration count in half in the remaining O(polygons * edges) part

common/advanced_config.cpp

@@ -109,6 +109,7 @@ static const wxChar OcePluginLinearDeflection[] = wxT( "OcePluginLinearDeflectio
 static const wxChar OcePluginAngularDeflection[] = wxT( "OcePluginAngularDeflection" );
 static const wxChar TriangulateSimplificationLevel[] = wxT( "TriangulateSimplificationLevel" );
 static const wxChar TriangulateMinimumArea[] = wxT( "TriangulateMinimumArea" );
+static const wxChar EnableCacheFriendlyFracture[] = wxT( "EnableCacheFriendlyFracture" );
 } // namespace KEYS
 
 
@@ -260,6 +261,8 @@ ADVANCED_CFG::ADVANCED_CFG()
     m_TriangulateSimplificationLevel = 50;
     m_TriangulateMinimumArea = 1000;
 
+    m_EnableCacheFriendlyFracture = true;
+
     loadFromConfigFile();
 }
 
@@ -475,6 +478,10 @@ void ADVANCED_CFG::loadSettings( wxConfigBase& aCfg )
                                                     &m_TriangulateMinimumArea,
                                                     m_TriangulateMinimumArea, 0, 100000 ) );
 
+    configParams.push_back( new PARAM_CFG_BOOL( true, AC_KEYS::EnableCacheFriendlyFracture,
+                                                &m_EnableCacheFriendlyFracture,
+                                                m_EnableCacheFriendlyFracture ) );
+
     // Special case for trace mask setting...we just grab them and set them immediately
     // Because we even use wxLogTrace inside of advanced config
     wxString traceMasks;

include/advanced_config.h

@@ -566,7 +566,17 @@ public:
      */
     int m_TriangulateMinimumArea;
 
-///@}
+    /**
+     * Enable the use of a cache-friendlier and therefore faster version of the
+     * polygon fracture algorithm.
+     *
+     * Setting name: "EnableCacheFriendlyFracture"
+     * Valid values: 0 or 1
+     * Default value: 1
+     */
+    bool m_EnableCacheFriendlyFracture;
+
+    ///@}
 
 
 private:

libs/kimath/src/geometry/shape_poly_set.cpp

@@ -37,9 +37,9 @@
 #include <map>
 #include <memory>
 #include <set>
-#include <string>                            // for char_traits, operator!=
-#include <type_traits>                       // for swap, move
+#include <string> // for char_traits, operator!=
 #include <unordered_set>
+#include <utility> // for swap, move
 #include <vector>
 
 #include <clipper.hpp>                       // for Clipper, PolyNode, Clipp...
@@ -1443,18 +1443,243 @@ void SHAPE_POLY_SET::importPaths( Clipper2Lib::Paths64&              aPath,
 
 struct FractureEdge
 {
-    FractureEdge( int y = 0 ) :
-        m_connected( false ),
-        m_next( nullptr )
+    using Index = int;
+
+    FractureEdge() = default;
+
+    FractureEdge( const VECTOR2I& p1, const VECTOR2I& p2, Index next ) :
+            m_p1( p1 ), m_p2( p2 ), m_next( next )
+    {
+    }
+
+    bool matches( int y ) const
+    {
+        return ( y >= m_p1.y || y >= m_p2.y ) && ( y <= m_p1.y || y <= m_p2.y );
+    }
+
+    VECTOR2I m_p1;
+    VECTOR2I m_p2;
+    Index    m_next;
+};
+
+
+typedef std::vector<FractureEdge> FractureEdgeSet;
+
+
+static FractureEdge* processHole( FractureEdgeSet& edges, FractureEdge::Index provokingIndex,
+                                  FractureEdge::Index edgeIndex, FractureEdge::Index bridgeIndex )
+{
+    FractureEdge& edge = edges[edgeIndex];
+    int           x = edge.m_p1.x;
+    int           y = edge.m_p1.y;
+    int           min_dist = std::numeric_limits<int>::max();
+    int           x_nearest = 0;
+
+    FractureEdge* e_nearest = nullptr;
+
+    // Since this function is run for all holes left to right, no need to
+    // check for any edge beyond the provoking one because they will always be
+    // further to the right, and unconnected to the outline anyway.
+    for( FractureEdge::Index i = 0; i < provokingIndex; i++ )
+    {
+        FractureEdge& e = edges[i];
+        // Don't consider this edge if it can't be bridged to, or faces left.
+        if( !e.matches( y ) )
+            continue;
+
+        int x_intersect;
+
+        if( e.m_p1.y == e.m_p2.y ) // horizontal edge
+        {
+            x_intersect = std::max( e.m_p1.x, e.m_p2.x );
+        }
+        else
+        {
+            x_intersect =
+                    e.m_p1.x + rescale( e.m_p2.x - e.m_p1.x, y - e.m_p1.y, e.m_p2.y - e.m_p1.y );
+        }
+
+        int dist = ( x - x_intersect );
+
+        if( dist >= 0 && dist < min_dist )
+        {
+            min_dist = dist;
+            x_nearest = x_intersect;
+            e_nearest = &e;
+        }
+    }
+
+    if( e_nearest )
+    {
+        const FractureEdge::Index outline2hole_index = bridgeIndex;
+        const FractureEdge::Index hole2outline_index = bridgeIndex + 1;
+        const FractureEdge::Index split_index = bridgeIndex + 2;
+        // Make an edge between the split outline edge and the hole...
+        edges[outline2hole_index] = FractureEdge( VECTOR2I( x_nearest, y ), edge.m_p1, edgeIndex );
+        // ...between the hole and the edge...
+        edges[hole2outline_index] =
+                FractureEdge( edge.m_p1, VECTOR2I( x_nearest, y ), split_index );
+        // ...and between the split outline edge and the rest.
+        edges[split_index] =
+                FractureEdge( VECTOR2I( x_nearest, y ), e_nearest->m_p2, e_nearest->m_next );
+
+        // Perform the actual outline edge split
+        e_nearest->m_p2 = VECTOR2I( x_nearest, y );
+        e_nearest->m_next = outline2hole_index;
+
+        FractureEdge* last = &edge;
+        for( ; last->m_next != edgeIndex; last = &edges[last->m_next] )
+            ;
+        last->m_next = hole2outline_index;
+    }
+
+    return e_nearest;
+}
+
+
+static void fractureSingleCacheFriendly( SHAPE_POLY_SET::POLYGON& paths )
+{
+    FractureEdgeSet edges;
+    bool            outline = true;
+
+    if( paths.size() == 1 )
+        return;
+
+    size_t total_point_count = 0;
+
+    for( const SHAPE_LINE_CHAIN& path : paths )
+    {
+        total_point_count += path.PointCount();
+    }
+
+    if( total_point_count > std::numeric_limits<FractureEdge::Index>::max() )
+    {
+        wxLogWarning( wxT( "Polygon has more points than int limit" ) );
+        return;
+    }
+
+    // Reserve space in the edge set so pointers don't get invalidated during
+    // the whole fracture process; one for each original edge, plus 3 per
+    // path to join it to the outline.
+    edges.reserve( total_point_count + paths.size() * 3 );
+
+    // Sort the paths by their lowest X bound before processing them.
+    // This ensures the processing order for processEdge() is correct.
+    struct PathInfo
+    {
+        int                 path_or_provoking_index;
+        FractureEdge::Index leftmost;
+        int                 x;
+        int                 y_or_bridge;
+    };
+    std::vector<PathInfo> sorted_paths;
+    const int             paths_count = static_cast<int>( paths.size() );
+    sorted_paths.reserve( paths_count );
+
+    for( int path_index = 0; path_index < paths_count; path_index++ )
+    {
+        const SHAPE_LINE_CHAIN&      path = paths[path_index];
+        const std::vector<VECTOR2I>& points = path.CPoints();
+        const int                    point_count = static_cast<int>( points.size() );
+        int                          x_min = std::numeric_limits<int>::max();
+        int                          y_min = std::numeric_limits<int>::max();
+        int                          leftmost = -1;
+
+        for( int point_index = 0; point_index < point_count; point_index++ )
+        {
+            const VECTOR2I& point = points[point_index];
+            if( point.x < x_min )
+            {
+                x_min = point.x;
+                leftmost = point_index;
+            }
+            if( point.y < y_min )
+                y_min = point.y;
+        }
+
+        sorted_paths.emplace_back( PathInfo{ path_index, leftmost, x_min, y_min } );
+    }
+
+    std::sort( sorted_paths.begin() + 1, sorted_paths.end(),
+               []( const PathInfo& a, const PathInfo& b )
+               {
+                   if( a.x == b.x )
+                       return a.y_or_bridge < b.y_or_bridge;
+                   return a.x < b.x;
+               } );
+
+    FractureEdge::Index edge_index = 0;
+
+    for( PathInfo& path_info : sorted_paths )
+    {
+        const SHAPE_LINE_CHAIN&      path = paths[path_info.path_or_provoking_index];
+        const std::vector<VECTOR2I>& points = path.CPoints();
+        const size_t                 point_count = points.size();
+
+        // Index of the provoking (first) edge for this path
+        const FractureEdge::Index provoking_edge = edge_index;
+
+        for( size_t i = 0; i < point_count - 1; i++ )
+        {
+            edges.emplace_back( points[i], points[i + 1], edge_index + 1 );
+            edge_index++;
+        }
+
+        // Create last edge looping back to the provoking one.
+        edges.emplace_back( points[point_count - 1], points[0], provoking_edge );
+        edge_index++;
+
+        if( !outline )
+        {
+            // Repurpose the path sorting data structure to schedule the leftmost edge
+            // for merging to the outline, which will in turn merge the rest of the path.
+            path_info.path_or_provoking_index = provoking_edge;
+            path_info.y_or_bridge = edge_index;
+
+            // Reserve 3 additional edges to bridge with the outline.
+            edge_index += 3;
+            edges.resize( edge_index );
+        }
+
+        outline = false; // first path is always the outline
+    }
+
+    for( auto it = sorted_paths.begin() + 1; it != sorted_paths.end(); it++ )
     {
-        m_p1.x = m_p2.y = y;
+        auto edge = processHole( edges, it->path_or_provoking_index,
+                                 it->path_or_provoking_index + it->leftmost, it->y_or_bridge );
+
+        // If we can't handle the hole, the zone is broken (maybe)
+        if( !edge )
+        {
+            wxLogWarning( wxT( "Broken polygon, dropping path" ) );
+
+            return;
+        }
     }
 
-    FractureEdge( bool connected, const VECTOR2I& p1, const VECTOR2I& p2 ) :
-        m_connected( connected ),
-        m_p1( p1 ),
-        m_p2( p2 ),
-        m_next( nullptr )
+    paths.resize( 1 );
+    SHAPE_LINE_CHAIN& newPath = paths[0];
+
+    newPath.Clear();
+    newPath.SetClosed( true );
+
+    // Root edge is always at index 0
+    FractureEdge* e = &edges[0];
+
+    for( ; e->m_next != 0; e = &edges[e->m_next] )
+        newPath.Append( e->m_p1 );
+
+    newPath.Append( e->m_p1 );
+}
+
+
+struct FractureEdgeSlow
+{
+    FractureEdgeSlow( int y = 0 ) : m_connected( false ), m_next( nullptr ) { m_p1.x = m_p2.y = y; }
+
+    FractureEdgeSlow( bool connected, const VECTOR2I& p1, const VECTOR2I& p2 ) :
+            m_connected( connected ), m_p1( p1 ), m_p2( p2 ), m_next( nullptr )
     {
     }
 
@@ -1463,26 +1688,26 @@ struct FractureEdge
         return ( y >= m_p1.y || y >= m_p2.y ) && ( y <= m_p1.y || y <= m_p2.y );
     }
 
-    bool          m_connected;
-    VECTOR2I      m_p1;
-    VECTOR2I      m_p2;
-    FractureEdge* m_next;
+    bool              m_connected;
+    VECTOR2I          m_p1;
+    VECTOR2I          m_p2;
+    FractureEdgeSlow* m_next;
 };
 
 
-typedef std::vector<FractureEdge*> FractureEdgeSet;
+typedef std::vector<FractureEdgeSlow*> FractureEdgeSetSlow;
 
 
-static int processEdge( FractureEdgeSet& edges, FractureEdge* edge )
+static int processEdge( FractureEdgeSetSlow& edges, FractureEdgeSlow* edge )
 {
-    int x   = edge->m_p1.x;
-    int y   = edge->m_p1.y;
-    int min_dist    = std::numeric_limits<int>::max();
-    int x_nearest   = 0;
+    int x = edge->m_p1.x;
+    int y = edge->m_p1.y;
+    int min_dist = std::numeric_limits<int>::max();
+    int x_nearest = 0;
 
-    FractureEdge* e_nearest = nullptr;
+    FractureEdgeSlow* e_nearest = nullptr;
 
-    for( FractureEdge* e : edges )
+    for( FractureEdgeSlow* e : edges )
     {
         if( !e->matches( y ) )
             continue;
@@ -1495,17 +1720,17 @@ static int processEdge( FractureEdgeSet& edges, FractureEdge* edge )
         }
         else
         {
-            x_intersect = e->m_p1.x + rescale( e->m_p2.x - e->m_p1.x, y - e->m_p1.y,
-                                               e->m_p2.y - e->m_p1.y );
+            x_intersect = e->m_p1.x
+                          + rescale( e->m_p2.x - e->m_p1.x, y - e->m_p1.y, e->m_p2.y - e->m_p1.y );
         }
 
         int dist = ( x - x_intersect );
 
         if( dist >= 0 && dist < min_dist && e->m_connected )
         {
-            min_dist    = dist;
-            x_nearest   = x_intersect;
-            e_nearest   = e;
+            min_dist = dist;
+            x_nearest = x_intersect;
+            e_nearest = e;
         }
     }
 
@@ -1513,21 +1738,24 @@ static int processEdge( FractureEdgeSet& edges, FractureEdge* edge )
     {
         int count = 0;
 
-        FractureEdge* lead1 = new FractureEdge( true, VECTOR2I( x_nearest, y ), VECTOR2I( x, y ) );
-        FractureEdge* lead2 = new FractureEdge( true, VECTOR2I( x, y ), VECTOR2I( x_nearest, y ) );
-        FractureEdge* split_2 = new FractureEdge( true, VECTOR2I( x_nearest, y ), e_nearest->m_p2 );
+        FractureEdgeSlow* lead1 =
+                new FractureEdgeSlow( true, VECTOR2I( x_nearest, y ), VECTOR2I( x, y ) );
+        FractureEdgeSlow* lead2 =
+                new FractureEdgeSlow( true, VECTOR2I( x, y ), VECTOR2I( x_nearest, y ) );
+        FractureEdgeSlow* split_2 =
+                new FractureEdgeSlow( true, VECTOR2I( x_nearest, y ), e_nearest->m_p2 );
 
         edges.push_back( split_2 );
         edges.push_back( lead1 );
         edges.push_back( lead2 );
 
-        FractureEdge* link = e_nearest->m_next;
+        FractureEdgeSlow* link = e_nearest->m_next;
 
         e_nearest->m_p2 = VECTOR2I( x_nearest, y );
         e_nearest->m_next = lead1;
         lead1->m_next = edge;
 
-        FractureEdge* last;
+        FractureEdgeSlow* last;
 
         for( last = edge; last->m_next != edge; last = last->m_next )
         {
@@ -1536,8 +1764,8 @@ static int processEdge( FractureEdgeSet& edges, FractureEdge* edge )
         }
 
         last->m_connected = true;
-        last->m_next    = lead2;
-        lead2->m_next   = split_2;
+        last->m_next = lead2;
+        lead2->m_next = split_2;
         split_2->m_next = link;
 
         return count + 1;
@@ -1547,11 +1775,11 @@ static int processEdge( FractureEdgeSet& edges, FractureEdge* edge )
 }
 
 
-void SHAPE_POLY_SET::fractureSingle( POLYGON& paths )
+static void fractureSingleSlow( SHAPE_POLY_SET::POLYGON& paths )
 {
-    FractureEdgeSet edges;
-    FractureEdgeSet border_edges;
-    FractureEdge*   root = nullptr;
+    FractureEdgeSetSlow edges;
+    FractureEdgeSetSlow border_edges;
+    FractureEdgeSlow*   root = nullptr;
 
     bool first = true;
 
@@ -1563,9 +1791,9 @@ void SHAPE_POLY_SET::fractureSingle( POLYGON& paths )
     for( const SHAPE_LINE_CHAIN& path : paths )
     {
         const std::vector<VECTOR2I>& points = path.CPoints();
-        int pointCount = points.size();
+        int                          pointCount = points.size();
 
-        FractureEdge* prev = nullptr, * first_edge = nullptr;
+        FractureEdgeSlow *prev = nullptr, *first_edge = nullptr;
 
         int x_min = std::numeric_limits<int>::max();
 
@@ -1576,8 +1804,8 @@ void SHAPE_POLY_SET::fractureSingle( POLYGON& paths )
 
             // Do not use path.CPoint() here; open-coding it using the local variables "points"
             // and "pointCount" gives a non-trivial performance boost to zone fill times.
-            FractureEdge* fe = new FractureEdge( first, points[ i ],
-                                                 points[ i+1 == pointCount ? 0 : i+1 ] );
+            FractureEdgeSlow* fe = new FractureEdgeSlow( first, points[i],
+                                                         points[i + 1 == pointCount ? 0 : i + 1] );
 
             if( !root )
                 root = fe;
@@ -1604,22 +1832,22 @@ void SHAPE_POLY_SET::fractureSingle( POLYGON& paths )
                 num_unconnected++;
         }
 
-        first = false;    // first path is always the outline
+        first = false; // first path is always the outline
     }
 
     // keep connecting holes to the main outline, until there's no holes left...
     while( num_unconnected > 0 )
     {
-        int x_min = std::numeric_limits<int>::max();
+        int  x_min = std::numeric_limits<int>::max();
         auto it = border_edges.begin();
 
-        FractureEdge* smallestX = nullptr;
+        FractureEdgeSlow* smallestX = nullptr;
 
         // find the left-most hole edge and merge with the outline
         for( ; it != border_edges.end(); ++it )
         {
-            FractureEdge* border_edge = *it;
-            int xt = border_edge->m_p1.x;
+            FractureEdgeSlow* border_edge = *it;
+            int               xt = border_edge->m_p1.x;
 
             if( ( xt <= x_min ) && !border_edge->m_connected )
             {
@@ -1635,7 +1863,7 @@ void SHAPE_POLY_SET::fractureSingle( POLYGON& paths )
         {
             wxLogWarning( wxT( "Broken polygon, dropping path" ) );
 
-            for( FractureEdge* edge : edges )
+            for( FractureEdgeSlow* edge : edges )
                 delete edge;
 
             return;
@@ -1649,20 +1877,28 @@ void SHAPE_POLY_SET::fractureSingle( POLYGON& paths )
 
     newPath.SetClosed( true );
 
-    FractureEdge* e;
+    FractureEdgeSlow* e;
 
     for( e = root; e->m_next != root; e = e->m_next )
         newPath.Append( e->m_p1 );
 
     newPath.Append( e->m_p1 );
 
-    for( FractureEdge* edge : edges )
+    for( FractureEdgeSlow* edge : edges )
         delete edge;
 
     paths.push_back( std::move( newPath ) );
 }
 
 
+void SHAPE_POLY_SET::fractureSingle( POLYGON& paths )
+{
+    if( ADVANCED_CFG::GetCfg().m_EnableCacheFriendlyFracture )
+        return fractureSingleCacheFriendly( paths );
+    fractureSingleSlow( paths );
+}
+
+
 void SHAPE_POLY_SET::Fracture( POLYGON_MODE aFastMode )
 {
     Simplify( aFastMode );    // remove overlapping holes/degeneracy