Improve zone fill parallelism and indexing

Remove the per-zone mutex that serialized all layer fills of the
same zone. fillSingleZone() is read-only on zone state, so
concurrent fills of different layers are safe. Make m_needRefill
atomic and narrow the CacheTriangulation lock to the map lookup.

Also add an R-tree index to reduce queries over large zones.

Fixes https://gitlab.com/kicad/code/kicad/-/issues/23450

libs/kimath/include/geometry/poly_containment_index.h

@@ -0,0 +1,187 @@
+/*
+ * This program source code file is part of KiCad, a free EDA CAD application.
+ *
+ * Copyright The KiCad Developers, see AUTHORS.txt for contributors.
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 3
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, you may find one here:
+ * http://www.gnu.org/licenses/old-licenses/gpl-3.0.html
+ * or you may search the http://www.gnu.org website for the version 3 license,
+ * or you may write to the Free Software Foundation, Inc.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
+ */
+
+#ifndef POLY_CONTAINMENT_INDEX_H
+#define POLY_CONTAINMENT_INDEX_H
+
+#include <climits>
+#include <cstdint>
+#include <vector>
+
+#include <geometry/rtree.h>
+#include <geometry/seg.h>
+#include <geometry/shape_poly_set.h>
+#include <math/util.h>
+#include <math/vector2d.h>
+
+/**
+ * Spatial index for efficient point-in-polygon containment testing.
+ *
+ * Standard SHAPE_LINE_CHAIN::PointInside() is O(V) per query, ray-casting against every edge.
+ * For large polygons with many containment queries (e.g. testing thousands of via positions
+ * against zone fills with tens of thousands of vertices), this becomes a bottleneck.
+ *
+ * This class builds an R-tree of polygon edges so containment queries become O(log V + K)
+ * where K is the number of edges the horizontal ray actually crosses. The ray-crossing
+ * algorithm matches SHAPE_LINE_CHAIN_BASE::PointInside() exactly, including the accuracy
+ * semantics where aAccuracy > 1 falls back to edge-distance testing.
+ */
+class POLY_CONTAINMENT_INDEX
+{
+public:
+    POLY_CONTAINMENT_INDEX() = default;
+
+    /**
+     * Build the spatial index from a SHAPE_POLY_SET's outlines.
+     *
+     * Indexes every edge of every outline in the polygon set. Must be called before Contains().
+     * Only indexes outlines, not holes (zone fills are fractured and have no holes).
+     */
+    void Build( const SHAPE_POLY_SET& aPolySet )
+    {
+        m_outlineCount = aPolySet.OutlineCount();
+
+        for( int outlineIdx = 0; outlineIdx < m_outlineCount; outlineIdx++ )
+        {
+            const SHAPE_LINE_CHAIN& outline = aPolySet.COutline( outlineIdx );
+            int                     ptCount = outline.PointCount();
+
+            if( ptCount < 3 )
+                continue;
+
+            for( int j = 0; j < ptCount; j++ )
+            {
+                const VECTOR2I& p1 = outline.CPoint( j );
+                const VECTOR2I& p2 = outline.CPoint( ( j + 1 ) % ptCount );
+
+                intptr_t idx = static_cast<intptr_t>( m_segments.size() );
+                m_segments.push_back( { p1, p2, outlineIdx } );
+
+                int min[2] = { std::min( p1.x, p2.x ), std::min( p1.y, p2.y ) };
+                int max[2] = { std::max( p1.x, p2.x ), std::max( p1.y, p2.y ) };
+                m_tree.Insert( min, max, idx );
+            }
+        }
+    }
+
+    /**
+     * Test whether a point is inside the indexed polygon set.
+     *
+     * Uses the same ray-crossing algorithm as SHAPE_LINE_CHAIN_BASE::PointInside(). When
+     * aAccuracy > 1, also checks if the point is within aAccuracy distance of any edge
+     * (matching the PointOnEdge fallback behavior).
+     *
+     * @param aPt       The point to test.
+     * @param aAccuracy Distance threshold for edge-proximity fallback. Values <= 1 skip
+     *                  the edge test for performance.
+     * @return true if the point is inside any outline or (when aAccuracy > 1) within
+     *         aAccuracy distance of any edge.
+     */
+    bool Contains( const VECTOR2I& aPt, int aAccuracy = 0 ) const
+    {
+        if( m_segments.empty() )
+            return false;
+
+        // Most polygon sets have very few outlines, so use a stack buffer to avoid
+        // per-query heap allocation.
+        int  crossingsStack[8] = {};
+        int* crossings = crossingsStack;
+        std::vector<int> crossingsHeap;
+
+        if( m_outlineCount > 8 )
+        {
+            crossingsHeap.resize( m_outlineCount, 0 );
+            crossings = crossingsHeap.data();
+        }
+
+        // Only segments whose X extent reaches past aPt.x can produce a rightward ray crossing.
+        int searchMin[2] = { aPt.x, aPt.y };
+        int searchMax[2] = { INT_MAX, aPt.y };
+
+        m_tree.Search( searchMin, searchMax,
+            [&]( intptr_t idx ) -> bool
+            {
+                const EDGE& seg = m_segments[idx];
+                const VECTOR2I& p1 = seg.p1;
+                const VECTOR2I& p2 = seg.p2;
+
+                if( ( p1.y >= aPt.y ) == ( p2.y >= aPt.y ) )
+                    return true;
+
+                const VECTOR2I diff = p2 - p1;
+                const int d = rescale( diff.x, ( aPt.y - p1.y ), diff.y );
+
+                if( aPt.x - p1.x < d )
+                    crossings[seg.outlineIdx]++;
+
+                return true;
+            } );
+
+        for( int i = 0; i < m_outlineCount; i++ )
+        {
+            if( crossings[i] & 1 )
+                return true;
+        }
+
+        if( aAccuracy > 1 )
+        {
+            int edgeMin[2] = { aPt.x - aAccuracy, aPt.y - aAccuracy };
+            int edgeMax[2] = { aPt.x + aAccuracy, aPt.y + aAccuracy };
+            SEG::ecoord accuracySq = SEG::Square( aAccuracy );
+            bool onEdge = false;
+
+            m_tree.Search( edgeMin, edgeMax,
+                [&]( intptr_t idx ) -> bool
+                {
+                    const EDGE& seg = m_segments[idx];
+                    SEG s( seg.p1, seg.p2 );
+
+                    if( s.SquaredDistance( aPt ) <= accuracySq )
+                    {
+                        onEdge = true;
+                        return false;
+                    }
+
+                    return true;
+                } );
+
+            return onEdge;
+        }
+
+        return false;
+    }
+
+private:
+    struct EDGE
+    {
+        VECTOR2I p1;
+        VECTOR2I p2;
+        int      outlineIdx;
+    };
+
+    std::vector<EDGE>                m_segments;
+    RTree<intptr_t, int, 2, double>  m_tree;
+    int                              m_outlineCount = 0;
+};
+
+#endif // POLY_CONTAINMENT_INDEX_H

pcbnew/zone.cpp

@@ -169,7 +169,7 @@ void ZONE::InitDataFromSrcInCopyCtor( const ZONE& aZone, PCB_LAYER_ID aLayer )
     m_minIslandArea           = aZone.m_minIslandArea;
 
     m_isFilled                = aZone.m_isFilled;
-    m_needRefill              = aZone.m_needRefill;
+    m_needRefill              = aZone.m_needRefill.load();
     m_teardropType            = aZone.m_teardropType;
 
     m_thermalReliefGap        = aZone.m_thermalReliefGap;
@@ -1352,10 +1352,10 @@ void ZONE::swapData( BOARD_ITEM* aImage )
 
 void ZONE::CacheTriangulation( PCB_LAYER_ID aLayer )
 {
-    std::lock_guard<std::mutex> lock( m_filledPolysListMutex );
-
     if( aLayer == UNDEFINED_LAYER )
     {
+        std::lock_guard<std::mutex> lock( m_filledPolysListMutex );
+
         for( auto& [ layer, poly ] : m_FilledPolysList )
             poly->CacheTriangulation();
 
@@ -1363,8 +1363,22 @@ void ZONE::CacheTriangulation( PCB_LAYER_ID aLayer )
     }
     else
     {
-        if( m_FilledPolysList.count( aLayer ) )
-            m_FilledPolysList[ aLayer ]->CacheTriangulation();
+        // Grab a shared_ptr copy under the lock, then triangulate outside it.
+        // Each layer's SHAPE_POLY_SET is independent, so concurrent triangulation
+        // of different layers is safe once we have the shared_ptr.
+        std::shared_ptr<SHAPE_POLY_SET> poly;
+
+        {
+            std::lock_guard<std::mutex> lock( m_filledPolysListMutex );
+
+            auto it = m_FilledPolysList.find( aLayer );
+
+            if( it != m_FilledPolysList.end() )
+                poly = it->second;
+        }
+
+        if( poly )
+            poly->CacheTriangulation();
     }
 }
 

pcbnew/zone.h

@@ -26,6 +26,7 @@
 #define ZONE_H
 
 
+#include <atomic>
 #include <mutex>
 #include <vector>
 #include <map>
@@ -277,11 +278,6 @@ public:
         return m_outlinearea;
     }
 
-    std::mutex& GetLock()
-    {
-        return m_lock;
-    }
-
     int GetFillFlag( PCB_LAYER_ID aLayer )
     {
         std::lock_guard<std::mutex> lock( m_fillFlagsMutex );
@@ -893,7 +889,7 @@ protected:
      * m_needRefill = false does not imply filled areas are up to date, just
      * the zone was refilled after edition, and does not need refilling
      */
-    bool             m_needRefill;
+    std::atomic<bool> m_needRefill;
 
     int              m_thermalReliefGap;        // Width of the gap in thermal reliefs.
     int              m_thermalReliefSpokeWidth; // Width of the copper bridge in thermal reliefs.
@@ -940,8 +936,6 @@ protected:
     double                    m_area;              // The filled zone area
     double                    m_outlinearea;       // The outline zone area
 
-    /// Lock used for multi-threaded filling on multi-layer zones
-    std::mutex                m_lock;
 };
 
 

pcbnew/zone_filler.cpp

@@ -59,6 +59,7 @@ static const wxChar traceZoneFiller[] = wxT( "KICAD_ZONE_FILLER" );
 #include <geometry/convex_hull.h>
 #include <geometry/geometry_utils.h>
 #include <geometry/vertex_set.h>
+#include <geometry/poly_containment_index.h>
 #include <kidialog.h>
 #include <thread_pool.h>
 #include <math/util.h>      // for KiROUND
@@ -680,11 +681,6 @@ bool ZONE_FILLER::Fill( const std::vector<ZONE*>& aZones, bool aCheck, wxWindow*
 
                 // Now we're ready to fill.
                 {
-                    std::unique_lock<std::mutex> zoneLock( zone->GetLock(), std::try_to_lock );
-
-                    if( !zoneLock.owns_lock() )
-                        return 0;
-
                     SHAPE_POLY_SET fillPolys;
 
                     if( !fillSingleZone( zone, layer, fillPolys ) )
@@ -708,15 +704,8 @@ bool ZONE_FILLER::Fill( const std::vector<ZONE*>& aZones, bool aCheck, wxWindow*
                 PCB_LAYER_ID layer = aFillItem.second;
                 ZONE*        zone = aFillItem.first;
 
-                {
-                    std::unique_lock<std::mutex> zoneLock( zone->GetLock(), std::try_to_lock );
-
-                    if( !zoneLock.owns_lock() )
-                        return 0;
-
-                    zone->CacheTriangulation( layer );
-                    zone->SetFillFlag( layer, true );
-                }
+                zone->CacheTriangulation( layer );
+                zone->SetFillFlag( layer, true );
 
                 return 1;
             };
@@ -783,7 +772,6 @@ bool ZONE_FILLER::Fill( const std::vector<ZONE*>& aZones, bool aCheck, wxWindow*
 
         std::this_thread::sleep_for( std::chrono::milliseconds( 100 ) );
 
-
         if( m_progressReporter )
         {
             m_progressReporter->KeepRefreshing();
@@ -1282,6 +1270,70 @@ bool ZONE_FILLER::Fill( const std::vector<ZONE*>& aZones, bool aCheck, wxWindow*
     // The first pass (before filling) marks vias as ZLO_FORCE_FLASHED if they're within the
     // zone outline. However, if the fill doesn't actually reach the via (due to obstacles like
     // tracks), we should not flash the via. See https://gitlab.com/kicad/code/kicad/-/issues/22010
+    //
+    // Build a spatial index per filled zone-layer for O(log V) containment queries instead of
+    // O(V) ray-casting. This is critical for boards with large zone fills (many vertices) and
+    // many vias/pads.
+    struct INDEXED_ZONE
+    {
+        BOX2I                                       bbox;
+        std::unique_ptr<POLY_CONTAINMENT_INDEX>     index;
+    };
+
+    struct NET_LAYER_HASH
+    {
+        size_t operator()( const std::pair<int, PCB_LAYER_ID>& k ) const
+        {
+            return std::hash<int>()( k.first ) ^ ( std::hash<int>()( k.second ) << 16 );
+        }
+    };
+
+    std::unordered_map<std::pair<int, PCB_LAYER_ID>, std::vector<INDEXED_ZONE>, NET_LAYER_HASH>
+            filledZonesByNetLayer;
+
+    for( ZONE* zone : m_board->Zones() )
+    {
+        if( zone->GetIsRuleArea() )
+            continue;
+
+        for( PCB_LAYER_ID layer : zone->GetLayerSet() )
+        {
+            if( !zone->HasFilledPolysForLayer( layer ) )
+                continue;
+
+            const std::shared_ptr<SHAPE_POLY_SET>& fill = zone->GetFilledPolysList( layer );
+
+            if( fill->IsEmpty() )
+                continue;
+
+            INDEXED_ZONE iz;
+            iz.bbox = zone->GetBoundingBox();
+            iz.index = std::make_unique<POLY_CONTAINMENT_INDEX>();
+            iz.index->Build( *fill );
+            filledZonesByNetLayer[{ zone->GetNetCode(), layer }].push_back( std::move( iz ) );
+        }
+    }
+
+    auto zoneReachesPoint =
+            [&]( int aNetcode, PCB_LAYER_ID aLayer, const VECTOR2I& aCenter, int aRadius ) -> bool
+            {
+                auto it = filledZonesByNetLayer.find( { aNetcode, aLayer } );
+
+                if( it == filledZonesByNetLayer.end() )
+                    return false;
+
+                for( const INDEXED_ZONE& iz : it->second )
+                {
+                    if( !iz.bbox.GetInflated( aRadius ).Contains( aCenter ) )
+                        continue;
+
+                    if( iz.index->Contains( aCenter, aRadius ) )
+                        return true;
+                }
+
+                return false;
+            };
+
     for( PCB_TRACK* track : m_board->Tracks() )
     {
         if( track->Type() != PCB_VIA_T )
@@ -1298,42 +1350,11 @@ bool ZONE_FILLER::Fill( const std::vector<ZONE*>& aZones, bool aCheck, wxWindow*
             if( via->GetZoneLayerOverride( layer ) != ZLO_FORCE_FLASHED )
                 continue;
 
-            bool zoneReachesVia = false;
-
-            for( ZONE* zone : m_board->Zones() )
-            {
-                if( zone->GetIsRuleArea() )
-                    continue;
-
-                if( zone->GetNetCode() != netcode )
-                    continue;
-
-                if( !zone->IsOnLayer( layer ) )
-                    continue;
-
-                if( !zone->HasFilledPolysForLayer( layer ) )
-                    continue;
-
-                const std::shared_ptr<SHAPE_POLY_SET>& fill = zone->GetFilledPolysList( layer );
-
-                if( fill->IsEmpty() )
-                    continue;
-
-                // Check if the filled zone reaches the via hole. Use holeRadius as reach distance
-                // to match the pre-fill check logic.
-                if( fill->Contains( center, -1, holeRadius ) )
-                {
-                    zoneReachesVia = true;
-                    break;
-                }
-            }
-
-            if( !zoneReachesVia )
+            if( !zoneReachesPoint( netcode, layer, center, holeRadius ) )
                 via->SetZoneLayerOverride( layer, ZLO_FORCE_NO_ZONE_CONNECTION );
         }
     }
 
-    // Same logic for pads
     for( FOOTPRINT* footprint : m_board->Footprints() )
     {
         for( PAD* pad : footprint->Pads() )
@@ -1342,8 +1363,6 @@ bool ZONE_FILLER::Fill( const std::vector<ZONE*>& aZones, bool aCheck, wxWindow*
             int      netcode = pad->GetNetCode();
             LSET     layers = pad->GetLayerSet() & boardCuMask;
 
-            // For TH pads, use the hole radius as tolerance since the filled zone creates a
-            // thermal relief around the pad hole, similar to vias.
             int holeRadius = 0;
 
             if( pad->HasHole() )
@@ -1354,35 +1373,7 @@ bool ZONE_FILLER::Fill( const std::vector<ZONE*>& aZones, bool aCheck, wxWindow*
                 if( pad->GetZoneLayerOverride( layer ) != ZLO_FORCE_FLASHED )
                     continue;
 
-                bool zoneReachesPad = false;
-
-                for( ZONE* zone : m_board->Zones() )
-                {
-                    if( zone->GetIsRuleArea() )
-                        continue;
-
-                    if( zone->GetNetCode() != netcode )
-                        continue;
-
-                    if( !zone->IsOnLayer( layer ) )
-                        continue;
-
-                    if( !zone->HasFilledPolysForLayer( layer ) )
-                        continue;
-
-                    const std::shared_ptr<SHAPE_POLY_SET>& fill = zone->GetFilledPolysList( layer );
-
-                    if( fill->IsEmpty() )
-                        continue;
-
-                    if( fill->Contains( center, -1, holeRadius ) )
-                    {
-                        zoneReachesPad = true;
-                        break;
-                    }
-                }
-
-                if( !zoneReachesPad )
+                if( !zoneReachesPoint( netcode, layer, center, holeRadius ) )
                     pad->SetZoneLayerOverride( layer, ZLO_FORCE_NO_ZONE_CONNECTION );
             }
         }