@@ -3660,48 +3660,57 @@ namespace OrthoTree
3660
3660
3661
3661
private: // K Nearest Neighbor helpers
3662
3662
static inline void AddEntityDistance (
3663
- auto const & entities, TVector const & searchPoint, TContainer const & points, std::vector<EntityDistance>& neighborEntities, TGeometry maxDistance) noexcept
3663
+ auto const & entities,
3664
+ TVector const & searchPoint,
3665
+ TContainer const & points,
3666
+ std::vector<EntityDistance>& neighborEntities,
3667
+ std::size_t neighborNo,
3668
+ TGeometry& maxDistanceWithin) noexcept
3664
3669
{
3665
3670
for (auto const entityID : entities)
3666
3671
{
3667
3672
const auto distance = AD::Distance (searchPoint, detail::at (points, entityID));
3668
- if (distance < maxDistance)
3669
- {
3673
+
3674
+ // maxDistanceWithin is implemented for tolerance handling: distance should be less than maxDistanceWithin
3675
+ if (distance >= maxDistanceWithin)
3676
+ continue ;
3677
+
3678
+ if (neighborEntities.size () < neighborNo - 1 )
3670
3679
neighborEntities.push_back ({ { distance }, entityID });
3680
+ else
3681
+ {
3682
+ if (neighborEntities.size () < neighborNo)
3683
+ {
3684
+ std::make_heap (neighborEntities.begin (), neighborEntities.end ());
3685
+ neighborEntities.push_back ({ { distance }, entityID });
3686
+ }
3687
+ else
3688
+ {
3689
+ std::pop_heap (neighborEntities.begin (), neighborEntities.end ());
3690
+ neighborEntities.back () = { { distance }, entityID };
3691
+ }
3692
+ std::push_heap (neighborEntities.begin (), neighborEntities.end ());
3693
+ maxDistanceWithin = neighborEntities.front ().Distance ;
3671
3694
}
3672
3695
}
3673
3696
}
3674
3697
3675
- static inline IGM::Geometry GetFarestDistance (std::vector<EntityDistance>& neighborEntities, std::size_t neighborNo, IGM::Geometry maxDistance) noexcept
3676
- {
3677
- if (neighborEntities.size () < neighborNo)
3678
- {
3679
- return maxDistance;
3680
- }
3681
-
3682
- auto const farestNeighborID = neighborNo - 1 ;
3683
- std::nth_element (neighborEntities.begin (), std::next (neighborEntities.begin (), farestNeighborID), neighborEntities.end ());
3684
- return neighborEntities[farestNeighborID].Distance ;
3685
- }
3686
-
3687
- static std::vector<TEntityID> ConvertEntityDistanceToList (std::vector<EntityDistance>& neighborEntities, std::size_t neighborNo) noexcept
3698
+ static inline constexpr std::vector<TEntityID> ConvertEntityDistanceToList (std::vector<EntityDistance>& neighborEntities, std::size_t neighborNo) noexcept
3688
3699
{
3689
3700
auto entityIDs = std::vector<TEntityID>();
3690
3701
if (neighborEntities.empty ())
3691
- {
3692
3702
return entityIDs;
3693
- }
3694
3703
3695
- if (neighborNo < neighborEntities.size ())
3696
- {
3697
- std::nth_element (neighborEntities. begin (), std::next (neighborEntities. begin (), neighborNo - 1 ), neighborEntities. end ());
3698
- }
3704
+ if (neighborEntities.size () < neighborNo )
3705
+ std::sort (neighborEntities. begin (), neighborEntities. end ());
3706
+ else
3707
+ std::sort_heap (neighborEntities. begin (), neighborEntities. end ());
3699
3708
3700
- auto const entityNo = std::min (neighborNo, neighborEntities.size ());
3701
- auto const lastIt = std::next (neighborEntities.begin (), entityNo);
3709
+ auto const entityNo = neighborEntities.size ();
3710
+ entityIDs.resize (entityNo);
3711
+ for (std::size_t i = 0 ; i < entityNo; ++i)
3712
+ entityIDs[i] = neighborEntities[i].EntityID ;
3702
3713
3703
- entityIDs.reserve (entityNo);
3704
- std::transform (neighborEntities.begin (), lastIt, std::back_inserter (entityIDs), [](auto const & ed) { return ed.EntityID ; });
3705
3714
return entityIDs;
3706
3715
}
3707
3716
@@ -3714,83 +3723,90 @@ namespace OrthoTree
3714
3723
return IGM::GetBoxWallDistanceAD (searchPoint, centerPoint, halfSize, isInsideConsideredAsZero);
3715
3724
}
3716
3725
3717
- void VisitNodesInDFSWithChildrenEdit (MortonNodeIDCR key, auto const & procedure, auto const & childNodeKeyEditor) const noexcept
3726
+ void VisitNodesInDFSWithChildrenEdit (
3727
+ depth_t stackID, std::pair<Node const *, TGeometry> const & nodeWithDistance, auto const & procedure, auto const & childNodeKeyEditor) const noexcept
3718
3728
{
3719
- auto const & node = this ->GetNode (key);
3720
- if (!procedure (node))
3729
+ if (!procedure (nodeWithDistance))
3721
3730
return ;
3722
3731
3723
- for (auto const [childKey, _] : childNodeKeyEditor (node.GetChildren ()))
3724
- this ->VisitNodesInDFSWithChildrenEdit (childKey, procedure, childNodeKeyEditor);
3732
+ auto const childStackID = stackID + 1 ;
3733
+ for (auto const & childNodeWithDistance : childNodeKeyEditor (nodeWithDistance.first ->GetChildren (), stackID))
3734
+ this ->VisitNodesInDFSWithChildrenEdit (childStackID, childNodeWithDistance, procedure, childNodeKeyEditor);
3725
3735
}
3726
3736
3727
3737
public:
3728
- // K Nearest Neighbor
3729
- std::vector<TEntityID> GetNearestNeighbors (TVector const & searchPoint, std::size_t neighborNo, TGeometry maxDistance, TContainer const & points) const noexcept
3738
+ // Get K Nearest Neighbor sorted by distance (point distance should be less than maxDistanceWithin, it is used as a Tolerance check)
3739
+ std::vector<TEntityID> GetNearestNeighbors (
3740
+ TVector const & searchPoint, std::size_t neighborNo, TGeometry maxDistanceWithin, TContainer const & points) const noexcept
3730
3741
{
3731
3742
auto neighborEntities = std::vector<EntityDistance>();
3732
- auto [smallestNodeKey, smallesDepthID] = this ->FindSmallestNodeKeyWithDepth (this ->template GetNodeID <true >(searchPoint));
3743
+ neighborEntities.reserve (neighborNo);
3744
+
3745
+ auto smallestNodeKey = this ->FindSmallestNodeKey (this ->template GetNodeID <true >(searchPoint));
3733
3746
if (!SI::IsValidKey (smallestNodeKey))
3734
- {
3735
3747
smallestNodeKey = SI::GetRootKey ();
3736
- smallesDepthID = 0 ;
3737
- }
3738
3748
3739
- auto const & smallestNode = this ->GetNode (smallestNodeKey);
3740
-
3741
- auto farestEntityDistance = maxDistance;
3749
+ auto farestEntityDistance = maxDistanceWithin;
3742
3750
// Parent checks (in a usual case parents do not have entities)
3743
3751
for (auto nodeKey = smallestNodeKey; SI::IsValidKey (nodeKey); nodeKey = SI::GetParentKey (nodeKey))
3744
- {
3745
- auto const & node = this ->GetNode (nodeKey);
3746
- AddEntityDistance (this ->GetNodeEntities (node), searchPoint, points, neighborEntities, farestEntityDistance);
3747
- farestEntityDistance = GetFarestDistance (neighborEntities, neighborNo, farestEntityDistance);
3748
- }
3752
+ AddEntityDistance (this ->GetNodeEntities (nodeKey), searchPoint, points, neighborEntities, neighborNo, farestEntityDistance);
3749
3753
3750
3754
// Search in itself and the children
3755
+ auto childrenDistanceStack = std::vector<std::vector<std::pair<Node const *, TGeometry>>>(this ->GetMaxDepthID ());
3751
3756
for (auto nodeKey = smallestNodeKey, prevNodeKey = MortonNodeID{}; SI::IsValidKey (nodeKey);
3752
3757
prevNodeKey = nodeKey, nodeKey = SI::GetParentKey (nodeKey))
3753
3758
{
3754
- auto const wallDistance = this ->GetNodeWallDistance (searchPoint, smallestNodeKey, smallestNode, false );
3759
+ auto const node = this ->GetNode (nodeKey);
3760
+ auto const wallDistance = this ->GetNodeWallDistance (searchPoint, nodeKey, node, false );
3755
3761
this ->VisitNodesInDFSWithChildrenEdit (
3756
- nodeKey,
3757
- [&](Node const & node) {
3758
- if (node.GetKey () == nodeKey)
3762
+ 0 ,
3763
+ { &node, wallDistance },
3764
+ [&](std::pair<Node const *, TGeometry> const & nodeDistance) {
3765
+ auto const & [childNode, childNodeDistance] = nodeDistance;
3766
+ auto const childNodeKey = childNode->GetKey ();
3767
+ if (childNodeKey == nodeKey)
3759
3768
return true ; // Self check was already done.
3760
3769
3761
- if (node. GetKey () == prevNodeKey)
3770
+ if (childNodeKey == prevNodeKey)
3762
3771
return false ; // Previous subtree was already checked.
3763
3772
3764
- AddEntityDistance (this ->GetNodeEntities (node), searchPoint, points, neighborEntities, farestEntityDistance);
3765
- farestEntityDistance = GetFarestDistance (neighborEntities, neighborNo, farestEntityDistance);
3773
+ if (childNodeDistance > farestEntityDistance)
3774
+ return false ;
3775
+
3776
+ AddEntityDistance (this ->GetNodeEntities (*childNode), searchPoint, points, neighborEntities, neighborNo, farestEntityDistance);
3766
3777
return true ;
3767
3778
},
3768
- [&](auto const & children) {
3769
- auto childrenDistance = std::vector<std::pair<MortonNodeID, TGeometry>>();
3779
+ [&](auto const & children, depth_t stackID) -> std::span<std::pair<Node const *, TGeometry>> {
3780
+ if (children.empty ())
3781
+ return {};
3782
+
3783
+ auto & childrenDistance = childrenDistanceStack[stackID];
3784
+ childrenDistance.clear ();
3770
3785
for (MortonNodeIDCR childNodeKey : children)
3771
3786
{
3772
3787
auto const & childNode = this ->m_nodes .at (childNodeKey);
3773
3788
auto const wallDistance = this ->GetNodeWallDistance (searchPoint, childNodeKey, childNode, true );
3774
3789
if (wallDistance > farestEntityDistance)
3775
3790
continue ;
3776
3791
3777
- childrenDistance.emplace_back (childNodeKey , wallDistance);
3792
+ childrenDistance.emplace_back (&childNode , wallDistance);
3778
3793
}
3779
3794
3780
3795
std::sort (childrenDistance.begin (), childrenDistance.end (), [&](auto const & leftDistance, auto const & rightDistance) {
3781
3796
return leftDistance.second < rightDistance.second ;
3782
3797
});
3783
3798
3784
- return childrenDistance;
3799
+ return std::span ( childrenDistance) ;
3785
3800
});
3786
3801
3787
- if (farestEntityDistance < wallDistance || maxDistance < wallDistance )
3802
+ if (farestEntityDistance < wallDistance)
3788
3803
break ;
3789
3804
}
3790
3805
3791
3806
return ConvertEntityDistanceToList (neighborEntities, neighborNo);
3792
3807
}
3793
3808
3809
+ // Get K Nearest Neighbor sorted by distance
3794
3810
inline std::vector<TEntityID> GetNearestNeighbors (TVector const & searchPoint, std::size_t neighborNo, TContainer const & points) const noexcept
3795
3811
{
3796
3812
return this ->GetNearestNeighbors (searchPoint, neighborNo, std::numeric_limits<TGeometry>::max (), points);
0 commit comments