Merge branch 'master' into e57plugin

doc/references.rst

@@ -38,8 +38,12 @@ Reference
 
 .. [Cook1986] Cook, Robert L. "Stochastic sampling in computer graphics." *ACM Transactions on Graphics (TOG)* 5.1 (1986): 51-72.
 
+.. [Demantke2011] Demantké J., Mallet C., David N., Vallet, B. "Dimensionality Based Scale Selection in 3d LIDAR Point Clouds."  Int. Arch. Photogramm. Remote Sens. Spatial Inf. Sci, XXXVIII-5/W12, 97-102, 2011
+
 .. [Dippe1985] Dippé, Mark AZ, and Erling Henry Wold. "Antialiasing through stochastic sampling." *ACM Siggraph Computer Graphics* 19.3 (1985): 69-78.
 
+.. [Guinard2017] Guinard S., Landrieu L. "Weakly Supervised Segmented-Aided Classification of Urban Scenes From 3D LIDAR Point Clouds." Int. Arch. Photogramm. Remote Sens. Spatial Inf. Sci., XLII-1/W1, 151-157, 2017
+
 .. [Kazhdan2006] Kazhdan, Michael, Matthew Bolitho, and Hugues Hoppe. "Poisson surface reconstruction." Proceedings of the fourth Eurographics symposium on Geometry processing. Vol. 7. 2006.
 
 .. [Li2010] Li, Ruosi, et al. "Polygonizing extremal surfaces with manifold guarantees." Proceedings of the 14th ACM Symposium on Solid and Physical Modeling. ACM, 2010.

doc/stages/filters.covariancefeatures.rst

@@ -0,0 +1,65 @@
+.. _filters.covariancefeatures:
+
+===============================================================================
+filters.covariancefeatures
+===============================================================================
+
+This filter implements various local feature descriptors introduced that are based on the covariance matrix of a
+point's neighborhood. The user can pick a set of feature descriptors by
+setting the ``feature_set`` option. Currently, the only supported feature is the dimensionality_
+set of feature descriptors introduced below.
+
+Example
+-------------------------------------------------------------------------------
+
+.. code-block:: json
+
+  [
+      "input.las",
+      {
+          "type":"filters.covariancefeatures",
+          "knn":8,
+          "threads": 2,
+          "feature_set": "Dimensionality"
+      },
+      {
+          "type":"writers.bpf",
+          "filename":"output.las",
+          "output_dims":"X,Y,Z,Linearity,Planarity,Scattering,Verticality"
+      }
+  ]
+
+Options
+-------------------------------------------------------------------------------
+
+knn
+  The number of k nearest neighbors used for calculating the covariance matrix. [Default: 10]
+
+threads
+  The number of threads used for computing the feature descriptors. [Default: 1]
+
+feature_set
+  The features to be computed. Currently only supports ``Dimensionality``. [Default: "Dimensionality"]
+
+.. _dimensionality:
+
+Dimensionality feature set
+................................................................................
+The features introduced in [Demantke2011]_ describe the shape
+of the neighborhood, indicating whether
+the local geometry is more linear (1D), planar (2D) or volumetric (3D) while the one introduced in
+[Guinard2017]_ adds the idea of a structure being vertical.
+
+
+The dimensionality filter introduces the following four descriptors that are computed from the covariance matrix of the ``knn`` neighbors:
+
+* linearity - higher for long thin strips
+* planarity - higher for planar surfaces
+* scattering - higher for complex 3d neighbourhoods
+* verticality - higher for vertical structures, highest for thin vertical strips
+
+It introduces four new dimensions that hold each one of these values: ``Linearity``  ``Planarity``  ``Scattering``
+and  ``Verticality``.
+
+
+

doc/stages/filters.hexbin.rst

@@ -115,5 +115,5 @@ threshold
   is considered "in" the data set. [Default: 15]
 
 precision
-  Coordinate precision to use in writing out the well-known text
-  of the boundary polygon. [Default: 8]
+  Minimum number of significant digits to use in writing out the
+  well-known text of the boundary polygon. [Default: 8]

doc/stages/filters.reprojection.rst

@@ -24,7 +24,7 @@ reprojecting.
 Example 1
 --------------------------------------------------------------------------------
 
-This pipeline reprojecs terrain points with Z-values between 0 and 100 by first
+This pipeline reprojects terrain points with Z-values between 0 and 100 by first
 applying a range filter and then specifing both the input and output spatial
 reference as EPSG-codes. The X and Y dimensions are scaled to allow enough
 precision in the output coordinates.

doc/stages/writers.ply.rst

@@ -42,12 +42,21 @@ storage_mode
   of the machine. [Default: "ascii"]
 
 dims
-  List of dimensions to write as elements. [Default: all dimensions]
+  List of dimensions (and sizes) 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
   Write a mesh as faces in addition to writing points as vertices.
   [Default: false]
 
+sized_types
+  PLY has variously been written with explicitly sized type strings
+  ('int8', 'float32", 'uint32', etc.) and implied sized type strings
+  ('char', 'float', 'int', etc.).  If true, explicitly sized type strings
+  are used.  If false, implicitly sized type strings are used.
+  [Default: true]
+
 precision
   If specified, the number of digits to the right of the decimal place
   using f-style formatting.  Only permitted when 'storage_mode' is 'ascii'.

doc/tutorial/iowa-entwine.rst

@@ -10,7 +10,7 @@ This tutorial describes how to use `Conda`_, `Entwine`_, `PDAL`_, and `GDAL`_ to
 read data from the `USGS 3DEP AWS Public Dataset`_. We will be using PDAL's
 `readers.ept`_ to fetch data, we will filter it for noise using `filters.outlier`_,
 we will classify the data as ground/not-ground using `filters.smrf`_, and we will
-write out a digital terrain model with `writers.gdal`. Once our elevation model
+write out a digital terrain model with :ref:`writers.gdal`. Once our elevation model
 is constructed, we will use GDAL `gdaldem`_ operations to create hillshade, slope,
 and color relief.
 
@@ -150,7 +150,7 @@ each point to ``0``.
 .. figure:: ../images/pipeline-example-filters.assign.png
    :scale: 50%
 
-   :ref`filters.assign` can also take in an option to apply assignments
+   :ref:`filters.assign` can also take in an option to apply assignments
    based on a conditional. If you want to assign values based on a
    bounding geometry, use :ref:`filters.overlay`.
 

filters/CovarianceFeaturesFilter.cpp

@@ -0,0 +1,155 @@
+/******************************************************************************
+* Copyright (c) 2019, Helix Re Inc. nicolas@helix.re
+*
+* 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 Helix Re Inc. 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.
+****************************************************************************/
+
+// This is an implementation of the local feature descriptors introduced in
+// WEAKLY SUPERVISED SEGMENTATION-AIDED CLASSIFICATION OF URBANSCENES FROM 3D LIDAR POINT CLOUDS
+// Stéphane Guinard, Loïc Landrieu, 2017
+
+#include "CovarianceFeaturesFilter.hpp"
+
+#include <pdal/EigenUtils.hpp>
+#include <pdal/KDIndex.hpp>
+#include <pdal/util/ProgramArgs.hpp>
+
+#include <Eigen/Dense>
+
+#include <string>
+#include <vector>
+#include <cmath>
+
+namespace pdal
+{
+
+static StaticPluginInfo const s_info
+{
+    "filters.covariancefeatures",
+    "Filter that calculates local features based on the covariance matrix of a point's neighborhood.",
+    "http://pdal.io/stages/filters.covariancefeatures.html"
+};
+
+CREATE_STATIC_STAGE(CovarianceFeaturesFilter, s_info)
+
+std::string CovarianceFeaturesFilter::getName() const
+{
+    return s_info.name;
+}
+
+void CovarianceFeaturesFilter::addArgs(ProgramArgs& args)
+{
+    args.add("knn", "k-Nearest neighbors", m_knn, 10);
+    args.add("threads", "Number of threads used to run this filter", m_threads, 1);
+    args.add("feature_set", "Set of features to be computed", m_featureSet, "Dimensionality");
+}
+
+void CovarianceFeaturesFilter::addDimensions(PointLayoutPtr layout)
+{
+    if (m_featureSet == "Dimensionality")
+    {
+        for (auto dim: {"Linearity", "Planarity", "Scattering", "Verticality"})
+            m_extraDims[dim] = layout->registerOrAssignDim(dim, Dimension::Type::Float);
+    }
+}
+
+void CovarianceFeaturesFilter::filter(PointView& view)
+{
+
+    KD3Index& kdi = view.build3dIndex();
+
+    point_count_t nloops = view.size();
+    std::vector<std::thread> threadPool(m_threads);
+    for(int t = 0;t<m_threads;t++)
+    {
+        threadPool[t] = std::thread(std::bind(
+                [&](const PointId start, const PointId end, const PointId t)
+                {
+                    for(PointId i = start;i<end;i++)
+                        setDimensionality(view, i, kdi);
+                },
+                t*nloops/m_threads,(t+1)==m_threads?nloops:(t+1)*nloops/m_threads,t));
+    }
+    for (auto &t: threadPool)
+        t.join();
+}
+
+void CovarianceFeaturesFilter::setDimensionality(PointView &view, const PointId &id, const KD3Index &kid)
+{
+    using namespace Eigen;
+
+    // find the k-nearest neighbors
+    auto ids = kid.neighbors(id, m_knn + 1);
+
+    // compute covariance of the neighborhood
+    auto B = eigen::computeCovariance(view, ids);
+
+    // perform the eigen decomposition
+    SelfAdjointEigenSolver<Matrix3f> solver(B);
+    if (solver.info() != Success)
+        throwError("Cannot perform eigen decomposition.");
+
+    // Extract eigenvalues and eigenvectors in decreasing order (largest eigenvalue first)
+    auto ev = solver.eigenvalues();
+    std::vector<float> lambda = {(std::max(ev[2],0.f)),
+                                 (std::max(ev[1],0.f)),
+                                 (std::max(ev[0],0.f))};
+
+    if (lambda[0] == 0)
+        throwError("Eigenvalues are all 0. Can't compute local features.");
+
+    auto eigenVectors = solver.eigenvectors();
+    std::vector<float> v1(3), v2(3), v3(3);
+    for (int i=0; i < 3; i++)
+    {
+        v1[i] = eigenVectors.col(2)(i);
+        v2[i] = eigenVectors.col(1)(i);
+        v3[i] = eigenVectors.col(0)(i);
+    }
+
+    float linearity  = (sqrtf(lambda[0]) - sqrtf(lambda[1])) / sqrtf(lambda[0]);
+    float planarity  = (sqrtf(lambda[1]) - sqrtf(lambda[2])) / sqrtf(lambda[0]);
+    float scattering =  sqrtf(lambda[2]) / sqrtf(lambda[0]);
+    view.setField(m_extraDims["Linearity"], id, linearity);
+    view.setField(m_extraDims["Planarity"], id, planarity);
+    view.setField(m_extraDims["Scattering"], id, scattering);
+
+    std::vector<float> unary_vector(3);
+    float norm = 0;
+    for (int i=0; i <3 ; i++)
+    {
+        unary_vector[i] = lambda[0] * fabsf(v1[i]) + lambda[1] * fabsf(v2[i]) + lambda[2] * fabsf(v3[i]);
+        norm += unary_vector[i] * unary_vector[i];
+    }
+    norm = sqrtf(norm);
+    view.setField(m_extraDims["Verticality"], id, unary_vector[2] / norm);
+}
+}

filters/CovarianceFeaturesFilter.hpp

@@ -0,0 +1,66 @@
+/******************************************************************************
+* Copyright (c) 2019, Helix Re Inc. nicolas@helix.re
+*
+* 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 Helix Re Inc. 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 <thread>
+
+#include <pdal/Filter.hpp>
+
+namespace pdal {
+
+class PDAL_DLL CovarianceFeaturesFilter: public Filter
+{
+public:
+    CovarianceFeaturesFilter() : Filter() {}
+    CovarianceFeaturesFilter &operator=(const CovarianceFeaturesFilter &) = delete;
+    CovarianceFeaturesFilter(const CovarianceFeaturesFilter &) = delete;
+
+    std::string getName() const;
+
+private:
+
+    int m_knn;
+    int m_threads;
+    std::string m_featureSet;
+    std::map<std::string,Dimension::Id> m_extraDims;
+
+    virtual void addDimensions(PointLayoutPtr layout);
+    virtual void addArgs(ProgramArgs &args);
+    virtual void filter(PointView &view);
+
+    void setDimensionality(PointView &view, const PointId &id, const KD3Index &kid);
+};
+}
+