Merge branch 'master' into issue-3282

CMakeLists.txt

@@ -129,7 +129,6 @@ include(${PDAL_CMAKE_DIR}/gdal.cmake)
 include(${PDAL_CMAKE_DIR}/geotiff.cmake)  # Optional (not really)
 include(${PDAL_CMAKE_DIR}/lazperf.cmake)  # Optional
 include(${PDAL_CMAKE_DIR}/laszip.cmake)  # Optional
-include(${PDAL_CMAKE_DIR}/draco.cmake)  # Optional
 include(${PDAL_CMAKE_DIR}/threads.cmake)
 include(${PDAL_CMAKE_DIR}/zlib.cmake)
 include(${PDAL_CMAKE_DIR}/lzma.cmake)

cmake/draco.cmake

@@ -1,19 +1,6 @@
 #
-option(WITH_DRACO "Choose if Draco support should be built" TRUE)
-find_package(Draco EXACT 1.3.6)
+# the `pkg_check_modules` function is created with this call
+find_package(PkgConfig REQUIRED)
+pkg_check_modules(DRACO REQUIRED draco>1.4.0)
 
-if (WITH_DRACO)
-    set_package_properties(Draco PROPERTIES TYPE RECOMMENDED
-        PURPOSE "Provides Draco compression")
-    if(Draco_FOUND)
-        include_directories(${draco_INCLUDE_DIR})
-        set(CMAKE_REQUIRED_LIBRARIES ${draco_LIBRARIES})
-        set(PDAL_HAVE_DRACO 1)
-        set(BUILD_PLUGIN_DRACO 1)
-        set(DRACO_LIBRARY "draco")
-    else()
-        set(WITH_DRACO FALSE)
-    endif()
-else()
-        set(WITH_DRACO FALSE)
-endif()
+set(PDAL_HAVE_DRACO 1)

doc/stages/filters.cluster.rst

@@ -46,3 +46,5 @@ 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]
+
+.. include:: filter_opts.rst

doc/stages/filters.zsmooth.rst

@@ -0,0 +1,52 @@
+.. _filters.zsmooth:
+
+filters.zsmooth
+===============================================================================
+
+The **Zsmooth Filter** computes a new Z value as another dimension that is based
+on the Z values of neighboring points.
+
+All points within some distance in the X-Y plane from a reference point are ordered by Z value.
+The reference point's new smoothed Z value is chosen to be that of the Nth median value of
+the neighboring points, where N is specified as the _`medianpercent` option.
+
+Use :ref:`filters.assign` to assign the smoothed Z value to the actual Z dimension if
+desired.
+
+Example
+-------
+
+Compute the smoothed Z value as the median Z value of the neighbors within 2 units and
+assign the value back to the Z dimension.
+
+[
+    "input.las",
+    {
+        "type": "filters.zsmooth",
+        "radius": 2,
+        "dim": "Zadj"
+    },
+    {
+        "type": "filters.assign",
+        "value": "Z = Zadj"
+    },
+    "output.las"
+]
+
+Options
+-------------------------------------------------------------------------------
+
+radius
+  All points within `radius` units from the reference point in the X-Y plane are considered
+  to determine the smoothed Z value. [Default: 1]
+
+medianpercent
+  A value between 0 and 100 that specifies the relative position of ordered Z values of neighbors
+  to use as the new smoothed Z value. 0 specifies the minimum value. 100 specifies the
+  maximum value. 50 specifies the mathematical median of the values. [Default: 50]
+
+dim
+  The name of a dimension to use for the adjusted Z value. Cannot be 'Z'. [Required]
+
+.. include:: filter_opts.rst
+

filters/ZsmoothFilter.cpp

@@ -0,0 +1,133 @@
+/******************************************************************************
+ * Copyright (c) 2020, University Nevada, Reno
+ *
+ * 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.
+ *
+ * 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/KDIndex.hpp>
+
+#include "ZsmoothFilter.hpp"
+
+namespace pdal
+{
+
+static PluginInfo const ptstatInfo
+{
+    "filters.zsmooth",
+    "Zsmooth Filter",
+    "http://pdal.io/stages/filters.zsmooth.html"
+};
+
+struct ZsmoothFilter::Private
+{
+    double radius;
+    double pos;
+    std::string dimName;
+    Dimension::Id statDim;
+};
+
+CREATE_SHARED_STAGE(ZsmoothFilter, ptstatInfo)
+
+std::string ZsmoothFilter::getName() const
+{
+    return ptstatInfo.name;
+}
+
+ZsmoothFilter::ZsmoothFilter() : m_p(new Private)
+{}
+
+
+ZsmoothFilter::~ZsmoothFilter()
+{}
+
+
+void ZsmoothFilter::addArgs(ProgramArgs& args)
+{
+    args.add("radius", "Radius in X/Y plane in which to find neighboring points", m_p->radius, 1.0);
+    args.add("medianpercent", "Location (percent) in neighbor list at which to find "
+        "neighbor Z value (min == 0, max == 100, median == 50, etc.)", m_p->pos, 50.0);
+    args.add("dim", "Name of dimension in which to store statistic", m_p->dimName).setPositional();
+}
+
+
+void ZsmoothFilter::addDimensions(PointLayoutPtr layout)
+{
+    m_p->statDim = layout->registerOrAssignDim(m_p->dimName, Dimension::Type::Double);
+    if (m_p->statDim == Dimension::Id::Z)
+        throwError("Can't use 'Z' as output dimension.");
+}
+
+
+void ZsmoothFilter::prepared(PointTableRef)
+{
+    if (m_p->pos < 0.0 || m_p->pos > 100.0)
+        throwError("'medicanpercent' value must be in the range [0, 100]");
+    m_p->pos /= 100.0;
+}
+
+
+void ZsmoothFilter::filter(PointView& view)
+{
+    const KD2Index& kdi = view.build2dIndex();
+
+    for (PointId idx = 0; idx < view.size(); ++idx)
+    {
+        double d = view.getFieldAs<double>(Dimension::Id::Z, idx);
+
+        std::vector<double> valList;
+        PointIdList nears = kdi.radius(idx, m_p->radius);
+        for (PointId n = 1; n < nears.size(); ++n)
+        {
+            double z = view.getFieldAs<double>(Dimension::Id::Z, nears[n]);
+            valList.push_back(z);
+        }
+        std::sort(valList.begin(), valList.end());
+
+        double val;
+        if (valList.empty())
+            val = view.getFieldAs<double>(Dimension::Id::Z, idx);
+        else if (valList.size() == 1)
+            val = valList[0];
+        else if (m_p->pos == 0.0)
+            val = valList[0];
+        else if (m_p->pos == 1.0)
+            val = valList[valList.size() - 1];
+        else
+        {
+            double pos = m_p->pos * (valList.size() - 1);
+            size_t low = std::floor(pos);
+            size_t high = low + 1;
+            double highfrac = pos - low;
+            double lowfrac = 1 - highfrac;
+            val = valList[low] * lowfrac + valList[high] * highfrac;
+
+        }
+        view.setField(m_p->statDim, idx, val);
+    }
+}
+
+} // namespace pdal

filters/ZsmoothFilter.hpp

@@ -0,0 +1,57 @@
+/***************************************************************************
+ * Copyright (c) 2020, University Nevada, Reno
+ * 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.
+ *
+ * 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 <memory>
+
+#include <pdal/Filter.hpp>
+
+namespace pdal
+{
+
+class ZsmoothFilter : public Filter
+{
+    struct Private;
+
+public:
+    std::string getName() const;
+
+    ZsmoothFilter();
+    ~ZsmoothFilter();
+
+private:
+    void addArgs(ProgramArgs& args);
+    void addDimensions(PointLayoutPtr layout);
+    void prepared(PointTableRef table);
+    void filter(PointView& view);
+
+    std::unique_ptr<Private> m_p;
+};
+
+} // namespace pdal
+

io/LasHeader.cpp

@@ -202,25 +202,26 @@ Dimension::IdList LasHeader::usedDims() const
 
 void LasHeader::setSrs()
 {
-    bool useWkt = false;
-
-    if (incompatibleSrs())
+    if (has14Format() && !useWkt())
     {
         m_log->get(LogLevel::Error) << "Global encoding WKT flag not set "
             "for point format 6 - 10." << std::endl;
     }
-    else if (findVlr(TRANSFORM_USER_ID, WKT_RECORD_ID) &&
+    if (findVlr(TRANSFORM_USER_ID, WKT_RECORD_ID) &&
         findVlr(TRANSFORM_USER_ID, GEOTIFF_DIRECTORY_RECORD_ID))
     {
         m_log->get(LogLevel::Debug) << "File contains both "
             "WKT and GeoTiff VLRs which is disallowed." << std::endl;
     }
-    else
-        useWkt = (m_versionMinor >= 4);
 
+    // We always use WKT for formats 6+, regardless of the WKT global encoding bit (warning
+    // issued above.)
+    // Otherwise (formats 0-5), we only use it of the WKT bit is set and it's version 1.4 
+    // or better.  For valid files the WKT bit won't be set for files < version 1.4, but
+    // we can't be sure, so we check here.
     try
     {
-        if (useWkt)
+        if ((useWkt() && m_versionMinor >= 4) || has14Format())
             setSrsFromWkt();
         else
             setSrsFromGeotiff();

io/LasHeader.hpp

@@ -347,9 +347,6 @@ class PDAL_DLL LasHeader
     bool useWkt() const
         { return (bool)((m_globalEncoding >> 4) & 1); }
 
-    bool incompatibleSrs() const
-        { return !useWkt() && has14Format(); }
-
     /// Returns true iff the file is compressed (laszip),
     /// as determined by the high bit in the point type
     bool compressed() const

io/LasWriter.cpp

@@ -935,7 +935,6 @@ bool LasWriter::writeLasZipBuf(PointRef& point)
     if (has14Format)
     {
         p.classification = (classification & 0x1F) | (classFlags << 5);
-        p.scan_angle_rank = point.getFieldAs<int8_t>(Id::ScanAngleRank);
         p.number_of_returns = (std::min)((uint8_t)7, numberOfReturns);
         p.return_number = (std::min)((uint8_t)7, returnNumber);
 
@@ -1169,7 +1168,7 @@ void LasWriter::finishOutput()
     OLeStream out(m_ostream);
 
     // addVlr prevents any eVlrs from being added before version 1.4.
-    m_lasHeader.setEVlrOffset((uint32_t)m_ostream->tellp());
+    m_lasHeader.setEVlrOffset(m_eVlrs.size() ? (uint32_t)m_ostream->tellp() : 0);
     for (auto vi = m_eVlrs.begin(); vi != m_eVlrs.end(); ++vi)
     {
         ExtLasVLR evlr = *vi;

io/private/ept/EptInfo.cpp

@@ -151,6 +151,7 @@ void EptInfo::initialize()
         double scale = element.value("scale", 1.0);
         double offset = element.value("offset", 0);
 
+        name = Dimension::fixName(name);
         Dimension::Id id = m_remoteLayout.registerOrAssignFixedDim(name, type);
         m_dims[name] = DimType(id, type, scale, offset);
     }

io/private/ept/TileContents.hpp

@@ -63,6 +63,10 @@ class TileContents
         { return m_overlap.m_key; }
     point_count_t nodeId() const
         { return m_overlap.m_nodeId; }
+    //ABELL - This is bad. We're assuming that the actual number of points we have matches
+    // what our index information told us. This may not be the case because of some
+    // issue. Downstream handling may depend on this being the actual number of points
+    // in the tile, rather than the number that were *supposed to be* in the tile.
     point_count_t size() const
         { return m_overlap.m_count; }
     const std::string& error() const