Merge remote-tracking branch 'origin/master' into issue-2545

apps/pdal.cpp

@@ -224,7 +224,7 @@ void App::outputOptions(std::string const& stageName, std::ostream& strm)
         {
             array = NL::json::parse(ostr.str());
         }
-        catch (NL::json::parse_error)
+        catch (NL::json::parse_error&)
         {}
 
         NL::json object = { stageName, array };

cmake/policies.cmake

@@ -4,4 +4,5 @@
 if (CMAKE_MAJOR_VERSION GREATER 2)
     cmake_policy(SET CMP0022 NEW) # interface link libraries
     cmake_policy(SET CMP0042 NEW) # osx rpath
+    cmake_policy(SET CMP0075 NEW) # Mystical setting to eliminate warning
 endif()

dimbuilder/DimBuilder.cpp

@@ -129,7 +129,7 @@ bool DimBuilder::execute()
     NL::json root;
     try
     {
-        root << in;
+        in >> root;
     }
     catch (NL::json::parse_error& err)
     {

doc/embed.py

@@ -41,17 +41,17 @@ def visit_admonition(self, node):
 def visit_embed_node(self, node):
     self.body.append(self.starttag(
             node, 'div', CLASS=('admonition embed')))
-    self.set_first_last(node)
+#    self.set_first_last(node)
 
 def visit_plugin_node(self, node):
     self.body.append(self.starttag(
             node, 'div', CLASS=('admonition plugin')))
-    self.set_first_last(node)
+#    self.set_first_last(node)
 
 def visit_streamable_node(self, node):
     self.body.append(self.starttag(
             node, 'div', CLASS=('admonition streamable')))
-    self.set_first_last(node)
+#    self.set_first_last(node)
 
 def depart_node(self, node):
     self.depart_admonition(node)

doc/index.rst

@@ -99,6 +99,14 @@ Dimensions
 
    dimensions
 
+Types
+--------------------------------------------------------------------------------
+
+.. toctree::
+   :maxdepth: 2
+
+   types
+
 Python
 --------------------------------------------------------------------------------
 

doc/stages/readers.las.rst

@@ -78,9 +78,8 @@ _`filename`
 _`extra_dims`
   Extra dimensions to be read as part of each point beyond those specified by
   the LAS point format.  The format of the option is
-  <dimension_name>=<type>, ... where type is one of:
-  int8, int16, int32, int64, uint8, uint16, uint32, uint64, float, double.
-  `_t` may be added to any of the type names as well (e.g., uint32_t).
+  ``<dimension_name>=<type>[, ...]``.  Any valid PDAL :ref:`type <types>` can be
+  specified.
 
   .. note::
 

doc/stages/readers.nitf.rst

@@ -62,9 +62,8 @@ filename
 extra_dims
   Extra dimensions to be read as part of each point beyond those specified by
   the LAS point format.  The format of the option is
-  <dimension_name>=<type>, ... where type is one of:
-  int8, int16, int32, int64, uint8, uint16, uint32, uint64, float, double.
-  `_t` may be added to any of the type names as well (e.g., uint32_t).
+  ``<dimension_name>=<type>[, ...]``.  Any PDAL :ref:`type <types>` can
+  be specified.
 
   .. note::
 

doc/stages/readers.rst

@@ -19,6 +19,7 @@ like :ref:`readers.oci`, or a network service like :ref:`readers.greyhound`.
    readers.bpf
    readers.buffer
    readers.ept
+   readers.e57
    readers.faux
    readers.gdal
    readers.geowave
@@ -59,6 +60,9 @@ like :ref:`readers.oci`, or a network service like :ref:`readers.greyhound`.
 :ref:`readers.ept`
     Used for reading `Entwine Point Tile <https://entwine.io>`__ format.
 
+:ref:`readers.e57`
+    Read point clouds in the E57 format.
+
 :ref:`readers.faux`
     Used for testing pipelines. It does not read from a file or database, but
     generates synthetic data to feed into the pipeline.

doc/stages/writers.gdal.rst

@@ -133,10 +133,10 @@ gdalopts
 .. _data_type:
 
 data_type
-    The data type to use for the output raster (double, float, int32,
-    uint16, etc.).  Many GDAL drivers only
-    support a limited set of output data types. The default value depends
-    on the driver.
+    The :ref:`data type <types>` to use for the output raster.
+    Many GDAL drivers only
+    support a limited set of output data types.
+    [Default: depends on the driver]
 
 .. _nodata:
 

doc/stages/writers.las.rst

@@ -187,17 +187,21 @@ discard_high_return_numbers
 extra_dims
   Extra dimensions to be written as part of each point beyond those specified
   by the LAS point format.  The format of the option is
-  <dimension_name>=<type>, ... where type is one of:
-  int8, int16, int32, int64, uint8, uint16, uint32, uint64, float, double
-  ``_t`` may be added to any of the type names as well (e.g., uint32_t).  When
-  the version of the output file is specified as 1.4 or greater, an extra
-  bytes VLR (User ID: LASF_Spec, Record ID: 4), is created that describes the
-  extra dimensions specified by this option.
+  ``<dimension_name>=<type> [, ...]``.  Any valid PDAL :ref:`type <types>`
+  can be specified.
 
   The special value ``all`` can be used in place of a dimension/type list
   to request that all dimensions that can't be stored in the predefined
   LAS point record get added as extra data at the end of each point record.
 
+  PDAL writes an extra bytes VLR (User ID: LASF_Spec, Record ID: 4) when
+  extra dims are written.  The VLR describes the extra dimensions specified by
+  this option.  Note that reading of this VLR is only specified for LAS
+  version 1.4, though some systems will honor it for earlier file formats.
+  The :ref:`LAS reader <readers.las>` requires the option
+  use_eb_vlr in order to
+  read the extra bytes VLR for files written with 1.1 - 1.3 LAS format.
+
   Setting --verbose=Info will provide output on the names, types and order
   of dimensions being written as part of the LAS extra bytes.
 

doc/stages/writers.ply.rst

@@ -42,8 +42,9 @@ storage_mode
   of the machine. [Default: "ascii"]
 
 dims
-  List of dimensions (and sizes) in the format <dimension_name>[=<type>],...
-  to write as output.  (e.g., "Y=int32_t, X,Red=char")
+  List of dimensions (and :ref:`types`) in the format
+  ``<dimension_name>[=<type>] [,...]`` to write as output.
+  (e.g., "Y=int32_t, X,Red=char")
   [Default: All dimensions with stored types]
 
 faces

doc/stages/writers.rst

@@ -18,6 +18,7 @@ dimension type, while others only understand fixed dimension names.
 
    writers.bpf
    writers.ept_addon
+   writers.e57
    writers.gdal
    writers.geowave
    writers.greyhound
@@ -41,6 +42,9 @@ dimension type, while others only understand fixed dimension names.
 :ref:`writers.ept_addon`
     Append additional dimensions to Entwine resources.
 
+:ref:`writers.e57`
+    Write data in the E57 format.
+
 :ref:`writers.gdal`
     Create a raster from a point cloud using an interpolation algorithm.
 

doc/type-table.csv

@@ -0,0 +1,10 @@
+"Signed Integer", "8", "``int8``, ``int8_t``, ``char``"
+"Signed Integer", "16", "``int16``, ``int16_t``, ``short``"
+"Signed Integer", "32", "``int32``, ``int32_t``, ``int``"
+"Signed Integer", "64", "``int64``, ``int64_t``, ``long``"
+"Unsigned Integer", "8", "``uint8``, ``uint8_t``, ``uchar``"
+"Unsigned Integer", "16", "``uint16``, ``uint16_t``, ``ushort``"
+"Unsigned Integer", "16", "``uint32``, ``uint32_t``, ``uint``"
+"Unsigned Integer", "16", "``uint64``, ``uint64_t``, ``ulong``"
+"Floating Point", "32", "``float``, ``float32``"
+"Floating Point", "64", "``double``, ``float64``"

doc/types.rst

@@ -0,0 +1,16 @@
+.. _types:
+
+===============================================================================
+Types
+===============================================================================
+
+PDAL supports the standard integral and floating point types for
+:ref:`dimensions <dimensions>`.  This table lists the types and associated
+strings
+that can be used to describe the types in options.
+
+
+.. csv-table::
+    :file: ./type-table.csv
+    :header: "Type", "Size in Bits", "Text Representations"
+    :widths: 10, 5, 40

filters/HAGFilter.cpp

@@ -36,8 +36,11 @@
 
 #include <pdal/KDIndex.hpp>
 
+#include "private/delaunator-decl.hpp"
+
 #include <string>
 #include <vector>
+#include <cmath>
 
 namespace pdal
 {
@@ -56,6 +59,17 @@ std::string HAGFilter::getName() const
     return s_info.name;
 }
 
+void HAGFilter::addArgs(ProgramArgs& args)
+{
+    args.add("count", "The number of points to fetch to determine the ground point [default: 1].",
+            m_count, point_count_t(1));
+    args.add("max_distance", "The maximum distance to the farthest nearest neighbor before the height above ground is not calculated [default: 0 (disabled)]", m_max_distance);
+    args.add("allow_extrapolation", "If true and count > 1, allow extrapolation [default: true].",
+            m_allow_extrapolation, true);
+    args.add("delaunay", "Construct local Delaunay fans and infer heights from them [default: false].",
+            m_delaunay, false);
+}
+
 void HAGFilter::addDimensions(PointLayoutPtr layout)
 {
     layout->registerDim(Dimension::Id::HeightAboveGround);
@@ -68,6 +82,32 @@ void HAGFilter::prepared(PointTableRef table)
         throwError("Missing Classification dimension in input PointView.");
 }
 
+// https://en.wikipedia.org/wiki/Barycentric_coordinate_system
+static double distance_along_z(double x1, double y1, double z1,
+                               double x2, double y2, double z2,
+                               double x3, double y3, double z3,
+                               double x, double y) {
+    double detT = (y2-y3)*(x1-x3)+(x3-x2)*(y1-y3);
+
+    if (detT == 0.0)
+    {
+        return std::numeric_limits<double>::infinity();
+    }
+
+    double lambda1 = ((y2-y3)*(x-x3)+(x3-x2)*(y-y3))/detT;
+    double lambda2 = ((y3-y1)*(x-x3)+(x1-x3)*(y-y3))/detT;
+    double lambda3 = 1 - lambda1 - lambda2;
+
+    if ((0.0 <= lambda1 && lambda1 <= 1.0) && (0.0 <= lambda2 && lambda2 <= 1.0) && (0.0 <= lambda3 && lambda3 <= 1.0))
+    {
+        return lambda1*z1 + lambda2*z2 + lambda3*z3;
+    }
+    else
+    {
+        return std::numeric_limits<double>::infinity();
+    }
+}
+
 void HAGFilter::filter(PointView& view)
 {
     PointViewPtr gView = view.makeNew();
@@ -99,13 +139,94 @@ void HAGFilter::filter(PointView& view)
     KD2Index& kdi = gView->build2dIndex();
 
     // Second pass: Find Z difference between non-ground points and the nearest
-    // neighbor (2D) in the ground view.
+    // neighbor (2D) in the ground view or between non-ground points and the
+    // locally-computed surface (Delaunay triangultion of the neighborhood).
     for (PointId i = 0; i < ngView->size(); ++i)
     {
         PointRef point = ngView->point(i);
+        double x0 = point.getFieldAs<double>(Dimension::Id::X);
+        double y0 = point.getFieldAs<double>(Dimension::Id::Y);
         double z0 = point.getFieldAs<double>(Dimension::Id::Z);
-        auto ids = kdi.neighbors(point, 1);
-        double z1 = gView->getFieldAs<double>(Dimension::Id::Z, ids[0]);
+        auto ids = kdi.neighbors(point, m_count);
+        double z1 = std::numeric_limits<double>::quiet_NaN();
+        assert(ids.size() > 0);
+        if (ids.size() == 1) {
+            z1 = gView->getFieldAs<double>(Dimension::Id::Z, ids[0]);
+        } else if (m_delaunay == false) { // Nearest-Neighbor-based interpolation
+            auto min_x = std::numeric_limits<double>::max();
+            auto max_x = std::numeric_limits<double>::min();
+            auto min_y = std::numeric_limits<double>::max();
+            auto max_y = std::numeric_limits<double>::min();
+            double weights = 0;
+            double z_accumulator = 0;
+            bool exact_match = false;
+            for (size_t j = 0; j < ids.size(); ++j) {
+                auto x = gView->getFieldAs<double>(Dimension::Id::X, ids[j]);
+                auto y = gView->getFieldAs<double>(Dimension::Id::Y, ids[j]);
+                auto z = gView->getFieldAs<double>(Dimension::Id::Z, ids[j]);
+                auto distance = std::sqrt(std::pow(x - x0, 2) + std::pow(y - y0, 2));
+                if (distance == 0) {
+                    exact_match = true;
+                    z1 = z;
+                    break;
+                }
+                if (m_max_distance > 0 && distance > m_max_distance) {
+                    break;
+                } else {
+                    auto weight = 1 / std::pow(distance, 2);
+                    weights += weight;
+                    z_accumulator += weight * z;
+                }
+                if (x > max_x) {
+                    max_x = x;
+                }
+                if (x < min_x) {
+                    min_x = x;
+                }
+                if (y > max_y) {
+                    max_y = y;
+                }
+                if (y < min_y) {
+                    min_y = y;
+                }
+            }
+            if (exact_match || m_allow_extrapolation || (x0 > min_x && x0 < max_x && y0 > min_y && y0 < max_y)) {
+                z1 = z_accumulator / weights;
+            }
+        } else if (m_delaunay == true) { // Delaunay-based interpolation
+            auto neighbors = std::vector<double>();
+
+            for(size_t j = 0; j < ids.size(); ++j)
+            {
+                neighbors.push_back(gView->getFieldAs<double>(Dimension::Id::X, ids[j]));
+                neighbors.push_back(gView->getFieldAs<double>(Dimension::Id::Y, ids[j]));
+            }
+            auto triangulation = delaunator::Delaunator(neighbors);
+            auto triangles = triangulation.triangles;
+
+            z1 = std::numeric_limits<double>::infinity();
+            for (size_t j = 0; j < triangles.size(); j += 3)
+            {
+                auto ai = triangles[j+0];
+                auto bi = triangles[j+1];
+                auto ci = triangles[j+2];
+                double ax = gView->getFieldAs<double>(Dimension::Id::X, ids[ai]);
+                double ay = gView->getFieldAs<double>(Dimension::Id::Y, ids[ai]);
+                double az = gView->getFieldAs<double>(Dimension::Id::Z, ids[ai]);
+                double bx = gView->getFieldAs<double>(Dimension::Id::X, ids[bi]);
+                double by = gView->getFieldAs<double>(Dimension::Id::Y, ids[bi]);
+                double bz = gView->getFieldAs<double>(Dimension::Id::Z, ids[bi]);
+                double cx = gView->getFieldAs<double>(Dimension::Id::X, ids[ci]);
+                double cy = gView->getFieldAs<double>(Dimension::Id::Y, ids[ci]);
+                double cz = gView->getFieldAs<double>(Dimension::Id::Z, ids[ci]);
+
+                z1 = std::min(z1, distance_along_z(ax, ay, az, bx, by, bz, cx, cy, cz, x0, y0));
+            }
+            if (z1 == std::numeric_limits<double>::infinity())
+            {
+                z1 = gView->getFieldAs<double>(Dimension::Id::Z, ids[0]);
+            }
+        }
         view.setField(Dimension::Id::HeightAboveGround, ngIdx[i], z0 - z1);
     }
 
@@ -114,5 +235,4 @@ void HAGFilter::filter(PointView& view)
         view.setField(Dimension::Id::HeightAboveGround, i, 0.0);
 }
 
-
 } // namespace pdal

filters/HAGFilter.hpp

@@ -57,12 +57,18 @@ class PDAL_DLL HAGFilter : public Filter
 
 private:
 
+    virtual void addArgs(ProgramArgs& args);
     virtual void addDimensions(PointLayoutPtr layout);
     virtual void prepared(PointTableRef table);
     virtual void filter(PointView& view);
 
     HAGFilter& operator=(const HAGFilter&); // not implemented
     HAGFilter(const HAGFilter&); // not implemented
+
+    bool m_allow_extrapolation;
+    bool m_delaunay;
+    double m_max_distance;
+    point_count_t m_count;
 };
 
 } // namespace pdal