Add option to perform Euclidean cluster extraction in 2D (#3084) (#3086)

* Add option to perform Euclidean cluster extraction in 2D (#3084)

Changed to return clusters as deque of PointIdLists (rather than
vector). Default behavior will continue to be to cluster in 3D.

* Templatize extractClusters

* Update Cluster filter docs for new is3d option

* Add tests for 2D clustering

* Addressing PR comments

doc/stages/filters.cluster.rst

@@ -42,3 +42,7 @@ tolerance
   Cluster tolerance - maximum Euclidean distance for a point to be added to the
   cluster. [Default: 1.0]
 
+is3d
+  By default, clusters are formed by considering neighbors in a 3D sphere, but
+  if ``is3d`` is set to false, it will instead consider neighbors in a 2D
+  cylinder (XY plane only). [Default: true]

filters/ClusterFilter.cpp

@@ -34,6 +34,8 @@
 
 #include "ClusterFilter.hpp"
 
+#include <pdal/KDIndex.hpp>
+
 #include "private/Segmentation.hpp"
 
 #include <string>
@@ -62,6 +64,7 @@ void ClusterFilter::addArgs(ProgramArgs& args)
     args.add("max_points", "Max points per cluster", m_maxPoints,
         (std::numeric_limits<uint64_t>::max)());
     args.add("tolerance", "Radius", m_tolerance, 1.0);
+    args.add("is3d", "Perform cluster extraction in 3D?", m_is3d, true);
 }
 
 void ClusterFilter::addDimensions(PointLayoutPtr layout)
@@ -71,8 +74,13 @@ void ClusterFilter::addDimensions(PointLayoutPtr layout)
 
 void ClusterFilter::filter(PointView& view)
 {
-    auto clusters = Segmentation::extractClusters(view, m_minPoints,
-        m_maxPoints, m_tolerance);
+    std::deque<PointIdList> clusters;
+    if (m_is3d)
+        clusters = Segmentation::extractClusters<KD3Index>(view, m_minPoints,
+            m_maxPoints, m_tolerance);
+    else
+        clusters = Segmentation::extractClusters<KD2Index>(view, m_minPoints,
+            m_maxPoints, m_tolerance);
 
     uint64_t id = 1;
     for (auto const& c : clusters)

filters/ClusterFilter.hpp

@@ -57,6 +57,7 @@ class PDAL_DLL ClusterFilter : public Filter
     uint64_t m_minPoints;
     uint64_t m_maxPoints;
     double m_tolerance;
+    bool m_is3d;
 
     virtual void addArgs(ProgramArgs& args);
     virtual void addDimensions(PointLayoutPtr layout);

filters/private/Segmentation.cpp

@@ -88,68 +88,6 @@ std::ostream& operator<<(std::ostream& out, const PointClasses& classes)
     return out;
 }
 
-std::vector<PointIdList> extractClusters(PointView& view, uint64_t min_points,
-                                         uint64_t max_points, double tolerance)
-{
-    // Index the incoming PointView for subsequent radius searches.
-    KD3Index kdi(view);
-    kdi.build();
-
-    // Create variables to track PointIds that have already been added to
-    // clusters and to build the list of cluster indices.
-    PointIdList processed(view.size(), 0);
-    std::vector<PointIdList> clusters;
-
-    for (PointId i = 0; i < view.size(); ++i)
-    {
-        // Points can only belong to a single cluster.
-        if (processed[i])
-            continue;
-
-        // Initialize list of indices belonging to current cluster, marking the
-        // seed point as processed.
-        PointIdList seed_queue;
-        size_t sq_idx = 0;
-        seed_queue.push_back(i);
-        processed[i] = 1;
-
-        // Check each point in the cluster for additional neighbors within the
-        // given tolerance, remembering that the list can grow if we add points
-        // to the cluster.
-        while (sq_idx < seed_queue.size())
-        {
-            // Find neighbors of the next cluster point.
-            PointId j = seed_queue[sq_idx];
-            PointIdList ids = kdi.radius(j, tolerance);
-
-            // The case where the only neighbor is the query point.
-            if (ids.size() == 1)
-            {
-                sq_idx++;
-                continue;
-            }
-
-            // Skip neighbors that already belong to a cluster and add the rest
-            // to this cluster.
-            for (auto const& k : ids)
-            {
-                if (processed[k])
-                    continue;
-                seed_queue.push_back(k);
-                processed[k] = 1;
-            }
-
-            sq_idx++;
-        }
-
-        // Keep clusters that are within the min/max number of points.
-        if (seed_queue.size() >= min_points && seed_queue.size() <= max_points)
-            clusters.push_back(seed_queue);
-    }
-
-    return clusters;
-}
-
 void ignoreDimRange(DimRange dr, PointViewPtr input, PointViewPtr keep,
                     PointViewPtr ignore)
 {

filters/private/Segmentation.hpp

@@ -89,12 +89,70 @@ std::ostream& operator<<(std::ostream& out, const PointClasses& classes);
   \param[in] min_points the minimum number of points in a cluster.
   \param[in] max_points the maximum number of points in a cluster.
   \param[in] tolerance the tolerance for adding points to a cluster.
-  \returns a vector of clusters (themselves vectors of PointIds).
+  \returns a deque of clusters (themselves vectors of PointIds).
 */
-PDAL_DLL std::vector<PointIdList> extractClusters(PointView& view,
-                                                  uint64_t min_points,
-                                                  uint64_t max_points,
-                                                  double tolerance);
+template <class KDINDEX>
+PDAL_DLL std::deque<PointIdList> extractClusters(PointView& view, uint64_t min_points,
+                                                 uint64_t max_points, double tolerance)
+{
+    // Index the incoming PointView for subsequent radius searches.
+    KDINDEX kdi(view);
+    kdi.build();
+
+    // Create variables to track PointIds that have already been added to
+    // clusters and to build the list of cluster indices.
+    PointIdList processed(view.size(), 0);
+    std::deque<PointIdList> clusters;
+
+    for (PointId i = 0; i < view.size(); ++i)
+    {
+        // Points can only belong to a single cluster.
+        if (processed[i])
+            continue;
+
+        // Initialize list of indices belonging to current cluster, marking the
+        // seed point as processed.
+        PointIdList seed_queue;
+        size_t sq_idx = 0;
+        seed_queue.push_back(i);
+        processed[i] = 1;
+
+        // Check each point in the cluster for additional neighbors within the
+        // given tolerance, remembering that the list can grow if we add points
+        // to the cluster.
+        while (sq_idx < seed_queue.size())
+        {
+            // Find neighbors of the next cluster point.
+            PointId j = seed_queue[sq_idx];
+            PointIdList ids = kdi.radius(j, tolerance);
+
+            // The case where the only neighbor is the query point.
+            if (ids.size() == 1)
+            {
+                sq_idx++;
+                continue;
+            }
+
+            // Skip neighbors that already belong to a cluster and add the rest
+            // to this cluster.
+            for (auto const& k : ids)
+            {
+                if (processed[k])
+                    continue;
+                seed_queue.push_back(k);
+                processed[k] = 1;
+            }
+
+            sq_idx++;
+        }
+
+        // Keep clusters that are within the min/max number of points.
+        if (seed_queue.size() >= min_points && seed_queue.size() <= max_points)
+            clusters.push_back(seed_queue);
+    }
+
+    return clusters;
+}
 
 PDAL_DLL void ignoreDimRange(DimRange dr, PointViewPtr input, PointViewPtr keep,
                              PointViewPtr ignore);

test/unit/SegmentationTest.cpp

@@ -33,6 +33,7 @@
 ****************************************************************************/
 
 #include <pdal/Dimension.hpp>
+#include <pdal/KDIndex.hpp>
 #include <pdal/pdal_test_main.hpp>
 #include <pdal/PointTable.hpp>
 #include <pdal/PointView.hpp>
@@ -47,7 +48,7 @@ TEST(SegmentationTest, BasicClustering)
 {
     using namespace Segmentation;
 
-    std::vector<PointIdList> clusters;
+    std::deque<PointIdList> clusters;
 
     PointTable table;
     PointLayoutPtr layout(table.layout());
@@ -62,15 +63,15 @@ TEST(SegmentationTest, BasicClustering)
     src->setField(Dimension::Id::X, 0, 0.0);
     src->setField(Dimension::Id::Y, 0, 0.0);
     src->setField(Dimension::Id::Z, 0, 0.0);
-    clusters = extractClusters(*src, 1, 10, 1.0);
+    clusters = extractClusters<KD3Index>(*src, 1, 10, 1.0);
     EXPECT_EQ(1u, clusters.size());
     EXPECT_EQ(1u, clusters[0].size());
 
     // Two separate clusters, both with single point
     src->setField(Dimension::Id::X, 1, 10.0);
     src->setField(Dimension::Id::Y, 1, 10.0);
     src->setField(Dimension::Id::Z, 1, 10.0);
-    clusters = extractClusters(*src, 1, 10, 1.0);
+    clusters = extractClusters<KD3Index>(*src, 1, 10, 1.0);
     EXPECT_EQ(2u, clusters.size());
     EXPECT_EQ(1u, clusters[0].size());
     EXPECT_EQ(1u, clusters[1].size());
@@ -79,18 +80,77 @@ TEST(SegmentationTest, BasicClustering)
     src->setField(Dimension::Id::X, 2, 0.5);
     src->setField(Dimension::Id::Y, 2, 0.5);
     src->setField(Dimension::Id::Z, 2, 0.5);
-    clusters = extractClusters(*src, 1, 10, 1.0);
+    clusters = extractClusters<KD3Index>(*src, 1, 10, 1.0);
     EXPECT_EQ(2u, clusters.size());
     EXPECT_EQ(2u, clusters[0].size());
     EXPECT_EQ(1u, clusters[1].size());
 
     // Reject the cluster with only one point
-    clusters = extractClusters(*src, 2, 10, 1.0);
+    clusters = extractClusters<KD3Index>(*src, 2, 10, 1.0);
     EXPECT_EQ(1u, clusters.size());
     EXPECT_EQ(2u, clusters[0].size());
 
     // Reject the cluster with two points
-    clusters = extractClusters(*src, 1, 1, 1.0);
+    clusters = extractClusters<KD3Index>(*src, 1, 1, 1.0);
+    EXPECT_EQ(1u, clusters.size());
+    EXPECT_EQ(1u, clusters[0].size());
+}
+
+TEST(SegmentationTest, Clustering2D)
+{
+    using namespace Segmentation;
+
+    std::deque<PointIdList> clusters;
+
+    PointTable table;
+    PointLayoutPtr layout(table.layout());
+
+    layout->registerDim(Dimension::Id::X);
+    layout->registerDim(Dimension::Id::Y);
+    layout->registerDim(Dimension::Id::Z);
+
+    PointViewPtr src(new PointView(table));
+
+    // Single point, single cluster
+    src->setField(Dimension::Id::X, 0, 0.0);
+    src->setField(Dimension::Id::Y, 0, 0.0);
+    src->setField(Dimension::Id::Z, 0, 0.0);
+    clusters = extractClusters<KD2Index>(*src, 1, 10, 1.0);
+    EXPECT_EQ(1u, clusters.size());
+    EXPECT_EQ(1u, clusters[0].size());
+
+    // Two separate clusters, both with single point
+    src->setField(Dimension::Id::X, 1, 10.0);
+    src->setField(Dimension::Id::Y, 1, 10.0);
+    src->setField(Dimension::Id::Z, 1, 10.0);
+    clusters = extractClusters<KD2Index>(*src, 1, 10, 1.0);
+    EXPECT_EQ(2u, clusters.size());
+    EXPECT_EQ(1u, clusters[0].size());
+    EXPECT_EQ(1u, clusters[1].size());
+
+    // Still two clusters, one with two points
+    src->setField(Dimension::Id::X, 2, 0.0);
+    src->setField(Dimension::Id::Y, 2, 0.0);
+    src->setField(Dimension::Id::Z, 2, 10.0);
+    clusters = extractClusters<KD2Index>(*src, 1, 10, 1.0);
+    EXPECT_EQ(2u, clusters.size());
+    EXPECT_EQ(2u, clusters[0].size());
+    EXPECT_EQ(1u, clusters[1].size());
+
+    // In 3D, should be three clusters
+    clusters = extractClusters<KD3Index>(*src, 1, 10, 1.0);
+    EXPECT_EQ(3u, clusters.size());
+    EXPECT_EQ(1u, clusters[0].size());
+    EXPECT_EQ(1u, clusters[1].size());
+    EXPECT_EQ(1u, clusters[2].size());
+
+    // Reject the cluster with only one point
+    clusters = extractClusters<KD2Index>(*src, 2, 10, 1.0);
+    EXPECT_EQ(1u, clusters.size());
+    EXPECT_EQ(2u, clusters[0].size());
+
+    // Reject the cluster with two points
+    clusters = extractClusters<KD2Index>(*src, 1, 1, 1.0);
     EXPECT_EQ(1u, clusters.size());
     EXPECT_EQ(1u, clusters[0].size());
 }