First pass at DRC RTree functionality

This implements a copper-layer RTree with functions for iterating over
the elements in a copper layer and providing Nearest Neighbor returns
for BOARD_CONNECTED_ITEMS

libs/kimath/include/geometry/shape.h

@@ -124,7 +124,7 @@ public:
      */
     virtual bool Collide( const VECTOR2I& aP, int aClearance = 0, int* aActual = nullptr ) const
     {
-        return Collide( SEG( aP, aP ), aClearance );
+        return Collide( SEG( aP, aP ), aClearance, aActual );
     }
 
     /**

libs/kimath/src/geometry/shape_compound.cpp

@@ -115,11 +115,21 @@ bool SHAPE_COMPOUND::IsSolid() const
 
 bool SHAPE_COMPOUND::Collide( const SEG& aSeg, int aClearance, int* aActual ) const
 {
+    int dist = std::numeric_limits<int>::max();
+
     for( auto& item : m_shapes )
     {
         if( item->Collide( aSeg, aClearance, aActual ) )
-            return true;
+        {
+            if( !aActual || *aActual == 0 )
+                return true;
+
+            dist = std::min( dist, *aActual );
+        }
     }
 
-    return false;
+    if( aActual )
+        *aActual = dist;
+
+    return dist != std::numeric_limits<int>::max();
 }

qa/drc_proto/CMakeLists.txt

@@ -66,6 +66,7 @@ include_directories(
     ${CMAKE_SOURCE_DIR}/pcbnew/dialogs
     ${CMAKE_SOURCE_DIR}/polygon
     ${CMAKE_SOURCE_DIR}/common/geometry
+    ${CMAKE_SOURCE_DIR}/libs/kimath/include/math
     ${CMAKE_SOURCE_DIR}/qa/common
     ${CMAKE_SOURCE_DIR}/qa
     ${CMAKE_SOURCE_DIR}/qa/qa_utils

qa/drc_proto/drc_rtree.h

@@ -0,0 +1,308 @@
+/*
+ * This program source code file is part of KiCad, a free EDA CAD application.
+ *
+ * Copyright (C) 2019 KiCad Developers, see AUTHORS.txt for contributors.
+ * Copyright (C) 2020 CERN
+ *
+ * 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 DRC_RTREE_H_
+#define DRC_RTREE_H_
+
+#include <eda_rect.h>
+#include <board_connected_item.h>
+#include <set>
+#include <vector>
+
+#include <geometry/rtree.h>
+#include <vector2d.h>
+
+/**
+ * DRC_RTREE -
+ * Implements an R-tree for fast spatial and layer indexing of connectable items.
+ * Non-owning.
+ */
+class DRC_RTREE
+{
+private:
+    using drc_rtree = RTree<BOARD_CONNECTED_ITEM*, int, 2, double>;
+
+public:
+    DRC_RTREE()
+    {
+        for( int layer : LSET::AllCuMask().Seq() )
+            m_tree[layer] = new drc_rtree();
+
+        m_count = 0;
+    }
+
+    ~DRC_RTREE()
+    {
+        for( auto tree : m_tree )
+            delete tree;
+    }
+
+    /**
+     * Function Insert()
+     * Inserts an item into the tree. Item's bounding box is taken via its GetBoundingBox() method.
+     */
+    void insert( BOARD_CONNECTED_ITEM* aItem )
+    {
+        const EDA_RECT& bbox    = aItem->GetBoundingBox();
+        const int       mmin[2] = { bbox.GetX(), bbox.GetY() };
+        const int       mmax[2] = { bbox.GetRight(), bbox.GetBottom() };
+
+        for( int layer : aItem->GetLayerSet().Seq() )
+            m_tree[layer]->Insert( mmin, mmax, aItem );
+
+        m_count++;
+    }
+
+    /**
+     * Function Remove()
+     * Removes an item from the tree. Removal is done by comparing pointers, attempting
+     * to remove a copy of the item will fail.
+     */
+    bool remove( BOARD_CONNECTED_ITEM* aItem )
+    {
+        // First, attempt to remove the item using its given BBox
+        const EDA_RECT& bbox    = aItem->GetBoundingBox();
+        const int       mmin[2] = { bbox.GetX(), bbox.GetY() };
+        const int       mmax[2] = { bbox.GetRight(), bbox.GetBottom() };
+        bool            removed = false;
+
+        for( auto layer : aItem->GetLayerSet().Seq() )
+        {
+            // If we are not successful ( true == not found ), then we expand
+            // the search to the full tree
+            if( m_tree[int( layer )]->Remove( mmin, mmax, aItem ) )
+            {
+                // N.B. We must search the whole tree for the pointer to remove
+                // because the item may have been moved before we have the chance to
+                // delete it from the tree
+                const int mmin2[2] = { INT_MIN, INT_MIN };
+                const int mmax2[2] = { INT_MAX, INT_MAX };
+
+                if( m_tree[int( layer )]->Remove( mmin2, mmax2, aItem ) )
+                    continue;
+            }
+
+            removed = true;
+        }
+
+        m_count -= int( removed );
+
+        return removed;
+    }
+
+    /**
+     * Function RemoveAll()
+     * Removes all items from the RTree
+     */
+    void clear()
+    {
+        for( auto tree : m_tree )
+            tree->RemoveAll();
+
+        m_count = 0;
+    }
+
+    /**
+     * Determine if a given item exists in the tree.  Note that this does not search the full tree
+     * so if the item has been moved, this will return false when it should be true.
+     *
+     * @param aItem Item that may potentially exist in the tree
+     * @param aRobust If true, search the whole tree, not just the bounding box
+     * @return true if the item definitely exists, false if it does not exist within bbox
+     */
+    bool contains( BOARD_CONNECTED_ITEM* aItem, bool aRobust = false )
+    {
+        const EDA_RECT& bbox    = aItem->GetBoundingBox();
+        const int       mmin[2] = { bbox.GetX(), bbox.GetY() };
+        const int       mmax[2] = { bbox.GetRight(), bbox.GetBottom() };
+        bool            found   = false;
+
+        auto search = [&found, &aItem]( const BOARD_CONNECTED_ITEM* aSearchItem ) {
+            if( aSearchItem == aItem )
+            {
+                found = true;
+                return false;
+            }
+
+            return true;
+        };
+
+        for( int layer : aItem->GetLayerSet().Seq() )
+        {
+            m_tree[layer]->Search( mmin, mmax, search );
+
+            if( found )
+                break;
+        }
+
+        if( !found && aRobust )
+        {
+            for( int layer : LSET::AllCuMask().Seq() )
+            {
+                // N.B. We must search the whole tree for the pointer to remove
+                // because the item may have been moved.  We do not expand the item
+                // layer search as this should not change.
+
+                const int mmin2[2] = { INT_MIN, INT_MIN };
+                const int mmax2[2] = { INT_MAX, INT_MAX };
+
+                m_tree[layer]->Search( mmin2, mmax2, search );
+
+                if( found )
+                    break;
+            }
+        }
+
+        return found;
+    }
+
+    std::vector<std::pair<int, BOARD_CONNECTED_ITEM*>> GetNearest( const wxPoint &aPoint,
+                                                                   PCB_LAYER_ID aLayer,
+                                                                   int aLimit )
+    {
+
+        const int point[2] = { aPoint.x, aPoint.y };
+        auto result = m_tree[int( aLayer )]->NearestNeighbors( point,
+                [aLimit]( std::size_t a_count, int a_maxDist ) -> bool
+                {
+                    return a_count >= aLimit;
+                },
+                []( BOARD_CONNECTED_ITEM* aElement) -> bool
+                {
+                    // Don't remove any elements from the list
+                    return false;
+                },
+                [aLayer]( const int* a_point, BOARD_CONNECTED_ITEM* a_data ) -> int
+                {
+                    switch( a_data->Type() )
+                    {
+                    case PCB_TRACE_T:
+                    {
+                        TRACK* track = static_cast<TRACK*>( a_data );
+                        SEG seg( track->GetStart(), track->GetEnd() );
+                        return seg.Distance( VECTOR2I( a_point[0], a_point[1] ) ) -
+                                           ( track->GetWidth() + 1 ) / 2;
+                    }
+                    case PCB_VIA_T:
+                    {
+                        VIA* via = static_cast<VIA*>( a_data );
+                        return ( VECTOR2I( via->GetPosition() ) -
+                                 VECTOR2I( a_point[0], a_point[1] ) ).EuclideanNorm() -
+                               ( via->GetWidth() + 1 ) / 2;
+
+                    }
+                    default:
+                    {
+                        VECTOR2I point( a_point[0], a_point[1] );
+                        int      dist = 0;
+                        auto     shape = a_data->GetEffectiveShape( aLayer );
+
+                        // Here we use a hack to get the distance by colliding with a large area
+                        // However, we can't use just MAX_INT because we will overflow the collision calculations
+                        shape->Collide( point, std::numeric_limits<int>::max() / 2, &dist);
+                        return dist;
+                    }
+                    }
+                    return 0;
+                });
+
+        return result;
+    }
+
+    /**
+     * Returns the number of items in the tree
+     * @return number of elements in the tree;
+     */
+    size_t size()
+    {
+        return m_count;
+    }
+
+    bool empty()
+    {
+        return m_count == 0;
+    }
+
+    using iterator = typename drc_rtree::Iterator;
+
+    /**
+     * The DRC_LAYER struct provides a layer-specific auto-range iterator to the RTree.  Using
+     * this struct, one can write lines like:
+     *
+     * for( auto item : rtree.OnLayer( In1_Cu ) )
+     *
+     * and iterate over only the RTree items that are on In1
+     */
+    struct DRC_LAYER
+    {
+        DRC_LAYER( drc_rtree* aTree ) : layer_tree( aTree )
+        {
+            m_rect = { { INT_MIN, INT_MIN }, { INT_MAX, INT_MAX } };
+        };
+
+        DRC_LAYER( drc_rtree* aTree, const EDA_RECT aRect ) : layer_tree( aTree )
+        {
+            m_rect = { { aRect.GetX(), aRect.GetY() },
+                       { aRect.GetRight(), aRect.GetBottom() } };
+        };
+
+        drc_rtree::Rect m_rect;
+        drc_rtree*      layer_tree;
+
+        iterator begin()
+        {
+            return layer_tree->begin( m_rect );
+        }
+
+        iterator end()
+        {
+            return layer_tree->end( m_rect );
+        }
+    };
+
+    DRC_LAYER OnLayer( PCB_LAYER_ID aLayer )
+    {
+        return DRC_LAYER( m_tree[int( aLayer )] );
+    }
+
+    DRC_LAYER Overlapping( PCB_LAYER_ID aLayer, const wxPoint& aPoint, int aAccuracy = 0 )
+    {
+        EDA_RECT rect( aPoint, wxSize( 0, 0 ) );
+        rect.Inflate( aAccuracy );
+        return DRC_LAYER( m_tree[int( aLayer )], rect );
+    }
+
+    DRC_LAYER Overlapping( PCB_LAYER_ID aLayer, const EDA_RECT& aRect )
+    {
+        return DRC_LAYER( m_tree[int( aLayer )], aRect );
+    }
+
+
+private:
+    drc_rtree*  m_tree[MAX_CU_LAYERS];
+    size_t      m_count;
+};
+
+
+#endif /* DRC_RTREE_H_ */

qa/drc_proto/drc_test_provider_copper_clearance.cpp

@@ -36,6 +36,7 @@
 
 #include <drc_proto/drc_engine.h>
 #include <drc_proto/drc_item.h>
+#include <drc_proto/drc_rtree.h>
 #include <drc_proto/drc_rule.h>
 #include <drc_proto/drc_test_provider_clearance_base.h>
 

thirdparty/rtree/geometry/rtree.h

@@ -15,10 +15,32 @@
 //    * 2004 Templated C++ port by Greg Douglas
 //    * 2013 CERN (www.cern.ch)
 //    * 2020 KiCad Developers - Add std::iterator support for searching
+//    * 2020 KiCad Developers - Add container nearest neighbor based on Hjaltason & Samet
 //
-//LICENSE:
+
-//
+/*
-//    Entirely free for all uses. Enjoy!
+ * This program source code file is part of KiCad, a free EDA CAD application.
+ *
+ * Copyright (C) 2020 KiCad Developers, see AUTHORS.txt for contributors.
+ * Copyright (C) 2013 CERN
+ *
+ * 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 RTREE_H
 #define RTREE_H
@@ -35,6 +57,8 @@
 #include <array>
 #include <functional>
 #include <iterator>
+#include <queue>
+#include <vector>
 
 #ifdef DEBUG
 #define ASSERT assert    // RTree uses ASSERT( condition )
@@ -202,23 +226,19 @@ public:
     /// Save tree contents to stream
     bool    Save( RTFileStream& a_stream );
 
-    /// Find the nearest neighbor of a specific point.
+    /**
-    /// It uses the MINDIST method, simplifying the one from "R-Trees: Theory and Applications" by Yannis Manolopoulos et al.
+     * Gets an ordered vector of the nearest data elements to a specified point
-    /// The bounding rectangle is used to calculate the distance to the DATATYPE.
+     * @param aPoint coordinate to measure against
-    /// \param a_point point to start the search
+     * @param aTerminate Callback routine to check when we have gathered sufficient elements
-    /// \return Returns the DATATYPE located closest to a_point, 0 if the tree is empty.
+     * @param aFilter Callback routine to remove specific elements from the query results
-    DATATYPE NearestNeighbor( const ELEMTYPE a_point[NUMDIMS] );
+     * @param aSquaredDist Callback routine to measure the distance from the point to the data element
-
+     * @return vector of matching elements and their distance to the point
-    /// Find the nearest neighbor of a specific point.
+     */
-    /// It uses the MINDIST method, simplifying the one from "R-Trees: Theory and Applications" by Yannis Manolopoulos et al.
+    std::vector<std::pair<ELEMTYPE, DATATYPE>> NearestNeighbors(
-    /// It receives a callback function to calculate the distance to a DATATYPE object, instead of using the bounding rectangle.
+            const ELEMTYPE aPoint[NUMDIMS],
-    /// \param a_point point to start the search
+            std::function<bool( const std::size_t aNumResults, const ELEMTYPE aMinDist )> aTerminate,
-    /// \param a_squareDistanceCallback function that performs the square distance calculation for the selected DATATYPE.
+            std::function<bool( const DATATYPE aElement )> aFilter,
-    /// \param a_squareDistance Pointer in which the square distance to the nearest neighbour will be returned.
+            std::function<ELEMTYPE( const ELEMTYPE a_point[NUMDIMS], const DATATYPE a_data )> aSquaredDist );
-    /// \return Returns the DATATYPE located closest to a_point, 0 if the tree is empty.
-    DATATYPE NearestNeighbor( const ELEMTYPE a_point[NUMDIMS],
-                              ELEMTYPE a_squareDistanceCallback( const ELEMTYPE a_point[NUMDIMS], DATATYPE a_data ),
-                              ELEMTYPE* a_squareDistance );
 
 public:
     /// Iterator is not remove safe.
@@ -495,6 +515,12 @@ protected:
         Branch m_branch;
         ELEMTYPE minDist;
         bool isLeaf;
+
+        inline bool operator<(const NNNode &other) const
+        {
+            /// This is reversed on purpose to use std::priority_queue
+            return other.minDist < minDist;
+        }
     };
 
     Node*           AllocNode();
@@ -531,8 +557,7 @@ protected:
     void            FreeListNode( ListNode* a_listNode );
     static bool     Overlap( Rect* a_rectA, Rect* a_rectB );
     void            ReInsert( Node* a_node, ListNode** a_listNode );
-    ELEMTYPE        MinDist( const ELEMTYPE a_point[NUMDIMS], Rect* a_rect );
+    ELEMTYPE        MinDist( const ELEMTYPE a_point[NUMDIMS], const Rect& a_rect );
-    void            InsertNNListSorted( std::vector<NNNode*>* nodeList, NNNode* newNode );
 
     bool Search( Node * a_node, Rect * a_rect, int& a_foundCount,
                  std::function<bool (const DATATYPE&)> a_callback ) const;
@@ -815,58 +840,67 @@ int RTREE_QUAL::Search( const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMD
 
 
 RTREE_TEMPLATE
-DATATYPE RTREE_QUAL::NearestNeighbor( const ELEMTYPE a_point[NUMDIMS] )
+std::vector<std::pair<ELEMTYPE, DATATYPE>> RTREE_QUAL::NearestNeighbors(
+        const ELEMTYPE a_point[NUMDIMS],
+        std::function<bool( const std::size_t aNumResults, const ELEMTYPE aMinDist )> aTerminate,
+        std::function<bool( const DATATYPE aElement )> aFilter,
+        std::function<ELEMTYPE( const ELEMTYPE a_point[NUMDIMS], const DATATYPE a_data )> aSquaredDist )
 {
-    return this->NearestNeighbor( a_point, 0, 0 );
+    std::vector<std::pair<ELEMTYPE, DATATYPE>> result;
-}
+    std::priority_queue<NNNode> search_q;
 
-
+    for( int i = 0; i < m_root->m_count; ++i )
-RTREE_TEMPLATE
-DATATYPE RTREE_QUAL::NearestNeighbor( const ELEMTYPE a_point[NUMDIMS],
-                                      ELEMTYPE a_squareDistanceCallback( const ELEMTYPE a_point[NUMDIMS], DATATYPE a_data ),
-                                      ELEMTYPE* a_squareDistance )
-{
-    std::vector<NNNode*> nodeList;
-    Node* node = m_root;
-    NNNode* closestNode = 0;
-    while( !closestNode || !closestNode->isLeaf )
     {
-        //check every node on this level
+        if( m_root->IsLeaf() )
-        for( int index = 0; index < node->m_count; ++index )
         {
-            NNNode* newNode = new NNNode;
+            search_q.push( NNNode{ m_root->m_branch[i],
-            newNode->isLeaf = node->IsLeaf();
+                               aSquaredDist( a_point, m_root->m_branch[i].m_data ),
-            newNode->m_branch = node->m_branch[index];
+                               m_root->IsLeaf() });
-            if( newNode->isLeaf && a_squareDistanceCallback )
-                newNode->minDist = a_squareDistanceCallback( a_point, newNode->m_branch.m_data );
-            else
-                newNode->minDist = this->MinDist( a_point, &(node->m_branch[index].m_rect) );
-
-            //TODO: a custom list could be more efficient than a vector
-            this->InsertNNListSorted( &nodeList, newNode );
         }
-        if( nodeList.size() == 0 )
+        else
         {
-            return 0;
+            search_q.push( NNNode{ m_root->m_branch[i],
+                               MinDist(a_point, m_root->m_branch[i].m_rect),
+                               m_root->IsLeaf() });
         }
-        closestNode = nodeList.back();
-        node = closestNode->m_branch.m_child;
-        nodeList.pop_back();
-        free(closestNode);
     }
 
-    // free memory used for remaining NNNodes in nodeList
+    while( !search_q.empty() )
-    for( auto node_it : nodeList )
     {
-        NNNode* nnode = node_it;
+        const NNNode curNode = search_q.top();
-        free(nnode);
+
+        if( aTerminate( result.size(), curNode.minDist ) )
+            break;
+
+        search_q.pop();
+
+        if( curNode.isLeaf )
+        {
+            if( aFilter( curNode.m_branch.m_data ) )
+                result.emplace_back( curNode.minDist, curNode.m_branch.m_data );
+        }
+        else
+        {
+            Node* node = curNode.m_branch.m_child;
+
+            for( int i = 0; i < node->m_count; ++i )
+            {
+                NNNode newNode;
+                newNode.isLeaf = node->IsLeaf();
+                newNode.m_branch = node->m_branch[i];
+                if( newNode.isLeaf )
+                    newNode.minDist = aSquaredDist( a_point, newNode.m_branch.m_data );
+                else
+                    newNode.minDist = this->MinDist( a_point, node->m_branch[i].m_rect );
+
+                search_q.push( newNode );
+            }
+        }
     }
 
-    *a_squareDistance = closestNode->minDist;
+    return result;
-    return closestNode->m_branch.m_data;
 }
 
-
 RTREE_TEMPLATE
 int RTREE_QUAL::Count()
 {
@@ -1895,7 +1929,7 @@ void RTREE_QUAL::ReInsert( Node* a_node, ListNode** a_listNode )
 }
 
 
-// Search in an index tree or subtree for all data retangles that overlap the argument rectangle.
+// Search in an index tree or subtree for all data rectangles that overlap the argument rectangle.
 RTREE_TEMPLATE
 bool RTREE_QUAL::Search( Node* a_node, Rect* a_rect, int& a_foundCount,
         std::function<bool (const DATATYPE&)> a_callback ) const
@@ -1939,44 +1973,34 @@ bool RTREE_QUAL::Search( Node* a_node, Rect* a_rect, int& a_foundCount,
 
 
 //calculate the minimum distance between a point and a rectangle as defined by Manolopoulos et al.
-//it uses the square distance to avoid the use of ELEMTYPEREAL values, which are slower.
+// returns Euclidean norm to ensure value fits in ELEMTYPE
 RTREE_TEMPLATE
-ELEMTYPE RTREE_QUAL::MinDist( const ELEMTYPE a_point[NUMDIMS], Rect* a_rect )
+ELEMTYPE RTREE_QUAL::MinDist( const ELEMTYPE a_point[NUMDIMS], const Rect& a_rect )
 {
-    ELEMTYPE *q, *s, *t;
+    const ELEMTYPE *q, *s, *t;
-    q = (ELEMTYPE*) a_point;
+    q = a_point;
-    s = a_rect->m_min;
+    s = a_rect.m_min;
-    t = a_rect->m_max;
+    t = a_rect.m_max;
-    int minDist = 0;
+    double minDist = 0.0;
+
     for( int index = 0; index < NUMDIMS; index++ )
     {
         int r = q[index];
+
         if( q[index] < s[index] )
         {
             r = s[index];
         }
-        else if( q[index] >t[index] )
+        else if( q[index] > t[index] )
         {
             r = t[index];
         }
-        int addend = q[index] - r;
+
+        double addend = q[index] - r;
         minDist += addend * addend;
     }
-    return minDist;
-}
 
-
+    return std::lround( std::sqrt( minDist ) );
-//insert a NNNode in a list sorted by its minDist (desc.)
-RTREE_TEMPLATE
-void RTREE_QUAL::InsertNNListSorted( std::vector<NNNode*>* nodeList, NNNode* newNode )
-{
-    auto iter = nodeList->begin();
-    while( iter != nodeList->end() && (*iter)->minDist > newNode->minDist )
-    {
-        ++iter;
-    }
-
-    nodeList->insert(iter, newNode);
 }