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
 #------------------------------------------------------------------------------

filters/CSFilter.cpp

@@ -0,0 +1,256 @@
+/******************************************************************************
+ * 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
+{
+    double m_cell;
+    double m_slope;
+    double m_window;
+    double m_scalar;
+    double m_threshold;
+    double m_cut;
+    std::string m_dir;
+    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("cell", "Cell size?", m_args->m_cell, 1.0);
+    args.add("slope", "Percent slope?", m_args->m_slope, 0.15);
+    args.add("window", "Max window size?", m_args->m_window, 18.0);
+    args.add("scalar", "Elevation scalar?", m_args->m_scalar, 1.25);
+    args.add("threshold", "Elevation threshold?", m_args->m_threshold, 0.5);
+    args.add("cut", "Cut net size?", m_args->m_cut, 0.0);
+    args.add("dir", "Optional output directory for debugging", m_args->m_dir);
+    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();
+        }
+    }
+}
+
+void CSFilter::ready(PointTableRef table)
+{
+    if (m_args->m_dir.empty())
+        return;
+
+    if (!FileUtils::directoryExists(m_args->m_dir))
+        throwError("Output directory '" + m_args->m_dir + "' does not exist");
+}
+
+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);
+    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

filters/CSFilter.hpp

@@ -0,0 +1,66 @@
+/******************************************************************************
+ * 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.
+ ****************************************************************************/
+
+#pragma once
+
+#include <pdal/Filter.hpp>
+
+namespace pdal
+{
+
+struct CSArgs;
+
+class PDAL_DLL CSFilter : public Filter
+{
+public:
+    CSFilter();
+
+    CSFilter& operator=(const CSFilter&) = delete;
+    CSFilter(const CSFilter&) = delete;
+
+    std::string getName() const;
+
+private:
+    std::unique_ptr<CSArgs> m_args;
+
+    virtual void addArgs(ProgramArgs& args);
+    virtual void addDimensions(PointLayoutPtr layout);
+    virtual void prepared(PointTableRef table);
+    virtual void ready(PointTableRef table);
+    virtual PointViewSet run(PointViewPtr view);
+
+    void classifyGround(PointViewPtr, std::vector<double>&);
+};
+
+} // namespace pdal