Add Cloth Simulation Filter (CSF)

CMakeLists.txt

@@ -125,6 +125,13 @@ include(${PDAL_CMAKE_DIR}/arbiter.cmake)
 include(${PDAL_CMAKE_DIR}/nlohmann.cmake)
 include(${PDAL_CMAKE_DIR}/openssl.cmake) # Optional
 
+find_package(OpenMP)
+if (OPENMP_FOUND)
+    set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
+    set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+    set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}")
+endif()
+
 #------------------------------------------------------------------------------
 # generate the pdal_features.hpp header
 #------------------------------------------------------------------------------

doc/references.rst

@@ -65,3 +65,5 @@ Reference
 .. [Weyrich2004] Weyrich, T et al. “Post-Processing of Scanned 3D Surface Data.” Proceedings of Eurographics Symposium on Point-Based Graphics 2004 (2004): 85–94. Print.
 
 .. [Zhang2003] Zhang, Keqi, et al. "A progressive morphological filter for removing nonground measurements from airborne LIDAR data." Geoscience and Remote Sensing, IEEE Transactions on 41.4 (2003): 872-882.
+
+.. [Zhang2016] Zhang, Wuming, et al. "An easy-to-use airborne LiDAR data filtering method based on cloth simulation." Remote Sensing 8.6 (2016): 501.

doc/stages/filters.csf.rst

@@ -0,0 +1,57 @@
+.. _filters.csf:
+
+filters.csf
+===============================================================================
+
+The **Cloth Simulation Filter (CSF)** classifies ground points based on the
+approach outlined in [Zhang2016]_.
+
+.. embed::
+
+Example
+-------
+
+The sample pipeline below uses CSF to segment ground and non-ground returns,
+using default options, and writing only the ground returns to the output file.
+
+.. code-block:: json
+
+  [
+      "input.las",
+      {
+          "type":"filters.csf"
+      },
+      {
+          "type":"filters.range",
+          "limits":"Classification[2:2]"
+      },
+      "output.laz"
+  ]
+
+Options
+-------------------------------------------------------------------------------
+
+resolution
+  Cloth resolution. [Default: **1.0**]
+
+ignore
+  A :ref:`range <ranges>` of values of a dimension to ignore.
+
+returns
+  Return types to include in output.  Valid values are "first", "last",
+  "intermediate" and "only". [Default: **"last, only"**]
+
+threshold
+  Classification threshold. [Default: **0.5**]
+
+smooth
+  Perform slope post-processing? [Default: **true**]
+
+step
+  Time step. [Default: **0.65**]
+
+rigidness
+  Rigidness. [Default: **3**]
+
+iterations
+  Maximum number of iterations. [Default: **500**]

doc/stages/filters.rst

@@ -37,6 +37,7 @@ invalidate an existing KD-tree.
    filters.colorinterp
    filters.colorization
    filters.covariancefeatures
+   filters.csf
    filters.dem
    filters.eigenvalues
    filters.estimaterank
@@ -81,6 +82,9 @@ invalidate an existing KD-tree.
     Filter that calculates local features based on the covariance matrix of a
     point's neighborhood.
 
+:ref:`filters.csf`
+    Label ground/non-ground returns using [Zhang2016]_.
+
 :ref:`filters.eigenvalues`
     Compute pointwise eigenvalues, based on k-nearest neighbors.
 

filters/CSFilter.cpp

@@ -0,0 +1,251 @@
+/******************************************************************************
+ * Copyright (c) 2019, Bradley J Chambers (brad.chambers@gmail.com)
+ *
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following
+ * conditions are met:
+ *
+ *     * Redistributions of source code must retain the above copyright
+ *       notice, this list of conditions and the following disclaimer.
+ *     * Redistributions in binary form must reproduce the above copyright
+ *       notice, this list of conditions and the following disclaimer in
+ *       the documentation and/or other materials provided
+ *       with the distribution.
+ *     * Neither the name of Hobu, Inc. or Flaxen Geo Consulting nor the
+ *       names of its contributors may be used to endorse or promote
+ *       products derived from this software without specific prior
+ *       written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
+ * OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
+ * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+ * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+ * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
+ * OF SUCH DAMAGE.
+ ****************************************************************************/
+
+// PDAL implementation of T. J. Pingel, K. C. Clarke, and W. A. McBride, “An
+// improved simple morphological filter for the terrain classification of
+// airborne LIDAR data,” ISPRS J. Photogramm. Remote Sens., vol. 77, pp. 21–30,
+// 2013.
+
+#include "CSFilter.hpp"
+
+#include <pdal/EigenUtils.hpp>
+#include <pdal/KDIndex.hpp>
+#include <pdal/util/FileUtils.hpp>
+#include <pdal/util/ProgramArgs.hpp>
+
+#include "private/DimRange.hpp"
+#include "private/Segmentation.hpp"
+#include "private/csf/CSF.h"
+
+#include <Eigen/Dense>
+
+#include <algorithm>
+#include <cmath>
+#include <iterator>
+#include <limits>
+#include <numeric>
+#include <string>
+#include <vector>
+
+namespace pdal
+{
+
+using namespace Dimension;
+using namespace Eigen;
+
+static StaticPluginInfo const s_info{
+    "filters.csf", "Cloth Simulation Filter (Zhang et al., 2016)",
+    "http://pdal.io/stages/filters.csf.html"};
+
+CREATE_STATIC_STAGE(CSFilter, s_info)
+
+struct CSArgs
+{
+    bool m_smooth;
+    double m_step;
+    double m_threshold;
+    double m_resolution;
+    int m_rigid;
+    int m_iterations;
+    std::vector<DimRange> m_ignored;
+    StringList m_returns;
+};
+
+CSFilter::CSFilter() : m_args(new CSArgs)
+{
+}
+
+std::string CSFilter::getName() const
+{
+    return s_info.name;
+}
+
+void CSFilter::addArgs(ProgramArgs& args)
+{
+    args.add("smooth", "Slope postprocessing?", m_args->m_smooth, true);
+    args.add("step", "Time step", m_args->m_step, 0.65);
+    args.add("threshold", "Classification threshold", m_args->m_threshold, 0.5);
+    args.add("resolution", "Cloth resolution", m_args->m_resolution, 1.0);
+    args.add("rigidness", "Rigidness", m_args->m_rigid, 3);
+    args.add("iterations", "Max iterations", m_args->m_iterations, 500);
+    args.add("ignore", "Ignore values", m_args->m_ignored);
+    args.add("returns", "Include last returns?", m_args->m_returns,
+             {"last", "only"});
+}
+
+void CSFilter::addDimensions(PointLayoutPtr layout)
+{
+    layout->registerDim(Id::Classification);
+}
+
+void CSFilter::prepared(PointTableRef table)
+{
+    const PointLayoutPtr layout(table.layout());
+
+    for (auto& r : m_args->m_ignored)
+    {
+        r.m_id = layout->findDim(r.m_name);
+        if (r.m_id == Dimension::Id::Unknown)
+            throwError("Invalid dimension name in 'ignored' option: '" +
+                       r.m_name + "'.");
+    }
+    if (m_args->m_returns.size())
+    {
+        for (auto& r : m_args->m_returns)
+        {
+            Utils::trim(r);
+            if ((r != "first") && (r != "intermediate") && (r != "last") &&
+                (r != "only"))
+            {
+                throwError("Unrecognized 'returns' value: '" + r + "'.");
+            }
+        }
+
+        if (!layout->hasDim(Dimension::Id::ReturnNumber) ||
+            !layout->hasDim(Dimension::Id::NumberOfReturns))
+        {
+            log()->get(LogLevel::Warning) << "Could not find ReturnNumber and "
+                                             "NumberOfReturns. Skipping "
+                                             "segmentation of last returns and "
+                                             "proceeding with all returns.\n";
+            m_args->m_returns.clear();
+        }
+    }
+}
+
+PointViewSet CSFilter::run(PointViewPtr view)
+{
+    PointViewSet viewSet;
+    if (!view->size())
+        return viewSet;
+
+    // Segment input view into ignored/kept views.
+    PointViewPtr ignoredView = view->makeNew();
+    PointViewPtr keptView = view->makeNew();
+    if (m_args->m_ignored.empty())
+        keptView->append(*view);
+    else
+        Segmentation::ignoreDimRanges(m_args->m_ignored, view, keptView,
+                                      ignoredView);
+
+    // Check for 0's in ReturnNumber and NumberOfReturns
+    bool nrOneZero(false);
+    bool rnOneZero(false);
+    bool nrAllZero(true);
+    bool rnAllZero(true);
+    for (PointId i = 0; i < keptView->size(); ++i)
+    {
+        uint8_t nr =
+            keptView->getFieldAs<uint8_t>(Dimension::Id::NumberOfReturns, i);
+        uint8_t rn =
+            keptView->getFieldAs<uint8_t>(Dimension::Id::ReturnNumber, i);
+        if ((nr == 0) && !nrOneZero)
+            nrOneZero = true;
+        if ((rn == 0) && !rnOneZero)
+            rnOneZero = true;
+        if (nr != 0)
+            nrAllZero = false;
+        if (rn != 0)
+            rnAllZero = false;
+    }
+
+    if ((nrOneZero || rnOneZero) && !(nrAllZero && rnAllZero))
+        throwError("Some NumberOfReturns or ReternNumber values were 0, but "
+                   "not all. Check that all values in the input file are >= "
+                   "1.");
+
+    // Segment kept view into two views
+    PointViewPtr firstView = keptView->makeNew();
+    PointViewPtr secondView = keptView->makeNew();
+
+    if (nrAllZero && rnAllZero)
+    {
+        log()->get(LogLevel::Warning)
+            << "Both NumberOfReturns and ReturnNumber are filled with 0's. "
+               "Proceeding without any further return filtering.\n";
+        firstView->append(*keptView);
+    }
+    else
+    {
+        Segmentation::segmentReturns(keptView, firstView, secondView,
+                                     m_args->m_returns);
+    }
+
+    if (!firstView->size())
+    {
+        throwError("No returns to process.");
+    }
+
+    for (PointId i = 0; i < secondView->size(); ++i)
+        secondView->setField(Dimension::Id::Classification, i, 1);
+
+    csf::PointCloud csfPC;
+    for (PointId i = 0; i < firstView->size(); ++i)
+    {
+        csf::Point p;
+        p.x = firstView->getFieldAs<double>(Dimension::Id::X, i);
+        p.y = firstView->getFieldAs<double>(Dimension::Id::Y, i);
+        p.z = firstView->getFieldAs<double>(Dimension::Id::Z, i);
+        csfPC.push_back(p);
+    }
+    std::cerr << "CSF: " << firstView->size() << std::endl;
+
+    CSF c(0);
+    c.params.bSloopSmooth = m_args->m_smooth;
+    c.params.time_step = m_args->m_step;
+    c.params.class_threshold = m_args->m_threshold;
+    c.params.cloth_resolution = m_args->m_resolution;
+    c.params.rigidness = m_args->m_rigid;
+    c.params.interations = m_args->m_iterations;
+    std::vector<int> groundIdx, offGroundIdx;
+    c.setPointCloud(csfPC);
+    c.do_filtering(groundIdx, offGroundIdx, true);
+    std::cerr << "CSF: " << groundIdx.size() << ", " << offGroundIdx.size()
+              << std::endl;
+
+    for (auto const& i : groundIdx)
+        firstView->setField(Dimension::Id::Classification, i, 2);
+    for (auto const& i : offGroundIdx)
+        firstView->setField(Dimension::Id::Classification, i, 1);
+
+    PointViewPtr outView = view->makeNew();
+    outView->append(*ignoredView);
+    outView->append(*secondView);
+    outView->append(*firstView);
+    viewSet.insert(outView);
+
+    return viewSet;
+}
+
+} // namespace pdal