Merge remote-tracking branch 'origin/master' into readers.numpy-script

filters/SkewnessBalancingFilter.cpp

@@ -56,18 +56,30 @@ void SkewnessBalancingFilter::addDimensions(PointLayoutPtr layout)
     layout->registerDim(Dimension::Id::Classification);
 }
 
-std::set<PointId> SkewnessBalancingFilter::processGround(PointViewPtr view)
+
+void SkewnessBalancingFilter::processGround(PointViewPtr view)
 {
     auto cmp = [](const PointIdxRef& p1, const PointIdxRef& p2) {
         return p1.compare(Dimension::Id::Z, p2);
     };
     std::sort(view->begin(), view->end(), cmp);
 
+    auto setClass = [&view](PointId first, PointId last, int cl)
+    {
+        for (PointId idx = first; idx <= last; ++idx)
+            view->setField(Dimension::Id::Classification, idx, cl);
+    };
+
+    const int Ground(2);
+    const int NonGround(1);
     point_count_t n(0);
     point_count_t n1(0);
     double delta, delta_n, term1, M1, M2, M3;
     M1 = M2 = M3 = 0.0;
-    std::vector<double> skewness(view->size());
+
+    PointId lastPositive = 0;
+    double skewness = 0;
+    double lastSkewness = std::numeric_limits<double>::quiet_NaN();
     for (PointId i = 0; i < view->size(); ++i)
     {
         double z = view->getFieldAs<double>(Dimension::Id::Z, i);
@@ -79,61 +91,36 @@ std::set<PointId> SkewnessBalancingFilter::processGround(PointViewPtr view)
         M1 += delta_n;
         M3 += term1 * delta_n * (n - 2) - 3 * delta_n * M2;
         M2 += term1;
-        skewness[i] = std::sqrt(n) * M3 / std::pow(M2, 1.5);
-    }
-
-    PointId j(0);
-    PointId i(view->size());
-    do
-    {
-        if (skewness[i] <= 0)
+        skewness = std::sqrt(n) * M3 / std::pow(M2, 1.5);
+        if (skewness > 0 && lastSkewness <= 0)
         {
-            j = i;
-            break;
+            setClass(lastPositive, i - 1, Ground);
+            lastPositive = i;
         }
-    } while (--i != 0);
-
-    std::set<PointId> groundIdx;
-    for (PointId i = 0; i <= j; ++i)
-        groundIdx.insert(i);
-
-    log()->get(LogLevel::Debug)
-        << "Stopped with " << groundIdx.size()
-        << " ground returns and skewness of " << skewness[j] << std::endl;
-
-    return groundIdx;
+        lastSkewness = skewness;
+    }
+    // It's possible that all our points have skewness <= 0, in which case
+    // we've never had an opportunity to set the ground state.  Do so now.
+    // Otherwise, set the remaining points to non-ground.
+    if (lastPositive == 0 && skewness <= 0)
+        setClass(lastPositive, view->size() - 1, Ground);
+    else
+        setClass(lastPositive, view->size() - 1, NonGround);
 }
 
+
 PointViewSet SkewnessBalancingFilter::run(PointViewPtr input)
 {
     PointViewSet viewSet;
     if (!input->size())
         return viewSet;
+    viewSet.insert(input);
 
     bool logOutput = log()->getLevel() > LogLevel::Debug1;
     if (logOutput)
         log()->floatPrecision(8);
 
-    auto idx = processGround(input);
-
-    if (!idx.empty())
-    {
-        // set the classification label of ground returns as 2
-        // (corresponding to ASPRS LAS specification)
-        for (const auto& i : idx)
-            input->setField(Dimension::Id::Classification, i, 2);
-
-        viewSet.insert(input);
-    }
-    else
-    {
-        if (idx.empty())
-            log()->get(LogLevel::Debug2)
-                << "Filtered cloud has no ground returns!\n";
-
-        // return the input buffer unchanged
-        viewSet.insert(input);
-    }
+    processGround(input);
 
     return viewSet;
 }

filters/SkewnessBalancingFilter.hpp

@@ -50,7 +50,7 @@ class PDAL_DLL SkewnessBalancingFilter : public Filter
 
 private:
     virtual void addDimensions(PointLayoutPtr layout);
-    std::set<PointId> processGround(PointViewPtr view);
+    void processGround(PointViewPtr view);
     virtual PointViewSet run(PointViewPtr view);
 
     SkewnessBalancingFilter&

test/unit/CMakeLists.txt

@@ -290,6 +290,7 @@ PDAL_ADD_TEST(pdal_filters_randomize_test FILES filters/RandomizeFilterTest.cpp)
 PDAL_ADD_TEST(pdal_filters_reciprocity_test FILES filters/ReciprocityFilterTest.cpp)
 PDAL_ADD_TEST(pdal_filters_returns_test FILES filters/ReturnsFilterTest.cpp)
 PDAL_ADD_TEST(pdal_filters_shell_test FILES filters/ShellFilterTest.cpp)
+PDAL_ADD_TEST(pdal_filters_skewness_test FILES filters/SkewnessFilterTest.cpp)
 
 PDAL_ADD_TEST(pdal_filters_smrf_test FILES filters/SMRFilterTest.cpp)
 PDAL_ADD_TEST(pdal_filters_sort_test

test/unit/apps/InfoTest.cpp

@@ -148,6 +148,8 @@ std::string r = R"foo(
 
     test("", r);
 
+// 10-Jan-20 - Broken by a change to proj which converts meters to ft, I think.
+/**
 std::string s = R"foo(
       "EPSG:4326":
       {
@@ -162,6 +164,7 @@ std::string s = R"foo(
         },
 )foo";
     test("", s);
+**/
 }
 
 TEST(Info, schema)
@@ -178,7 +181,6 @@ std::string r = R"foo(
         "type": "unsigned"
       },
 )foo";
-
     test("--schema", r);
 }
 
@@ -196,6 +198,5 @@ std::string r = R"foo(
         "type": "unsigned"
       },
 )foo";
-
     test("--all", r);
 }

test/unit/filters/SkewnessFilterTest.cpp

@@ -0,0 +1,138 @@
+/******************************************************************************
+* Copyright (c) 2020, Hobu Inc. (info@hobu.co)
+*
+* 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.
+****************************************************************************/
+
+#include <pdal/pdal_test_main.hpp>
+#include <pdal/StageFactory.hpp>
+#include <io/BufferReader.hpp>
+
+#include "Support.hpp"
+
+using namespace pdal;
+
+TEST(SkewnessTest, t1)
+{
+    StageFactory f;
+
+    Stage *reader(f.createStage("readers.las"));
+    Options rOpts;
+    rOpts.add("filename", Support::datapath("las/autzen_trim.las"));
+    reader->setOptions(rOpts);
+
+    Stage* filter(f.createStage("filters.skewnessbalancing"));
+    filter->setInput(*reader);
+
+    PointTable t;
+    filter->prepare(t);
+    PointViewSet s = filter->execute(t);
+
+    EXPECT_EQ(s.size(), 1U);
+    PointViewPtr v = *s.begin();
+
+    size_t ground {0};
+    for (PointId id = 0; id < v->size(); ++id)
+    {
+        uint8_t cl = v->getFieldAs<uint8_t>(Dimension::Id::Classification, id);
+        if (cl == 2)
+            ground++;
+    }
+    EXPECT_EQ(ground, 102234u);
+}
+
+TEST(SkewnessTest, t2)
+{
+    StageFactory f;
+
+    Stage *reader(f.createStage("readers.faux"));
+    Options rOpts;
+    rOpts.add("mode", "constant");
+    rOpts.add("count", 100);
+    rOpts.add("bounds", "([1, 2],[2, 3],[3, 4])");
+    reader->setOptions(rOpts);
+
+    Stage* filter(f.createStage("filters.skewnessbalancing"));
+    filter->setInput(*reader);
+
+    PointTable t;
+    filter->prepare(t);
+    PointViewSet s = filter->execute(t);
+
+    EXPECT_EQ(s.size(), 1U);
+    PointViewPtr v = *s.begin();
+
+    size_t ground {0};
+    for (PointId id = 0; id < v->size(); ++id)
+    {
+        uint8_t cl = v->getFieldAs<uint8_t>(Dimension::Id::Classification, id);
+        if (cl == 2)
+            ground++;
+    }
+    EXPECT_EQ(ground, 0);
+}
+
+
+TEST(SkewnessTest, t3)
+{
+    StageFactory f;
+    PointTable t;
+    t.layout()->registerDims(
+        {Dimension::Id::Z, Dimension::Id::Classification} );
+
+    PointViewPtr bv(new PointView(t));
+    for (int i = 0; i < 10; ++i)
+        bv->setField(Dimension::Id::Z, i, i - 5);
+
+    BufferReader reader;
+    reader.addView(bv);
+
+    Stage* filter(f.createStage("filters.skewnessbalancing"));
+    filter->setInput(reader);
+
+    filter->prepare(t);
+
+    PointViewSet s = filter->execute(t);
+
+    EXPECT_EQ(s.size(), 1U);
+    PointViewPtr v = *s.begin();
+
+    size_t ground {0};
+    for (PointId id = 0; id < v->size(); ++id)
+    {
+        uint8_t cl = v->getFieldAs<uint8_t>(Dimension::Id::Classification, id);
+        if (cl == 2)
+            ground++;
+    }
+    EXPECT_EQ(ground, 10);
+}
+
+