Merge remote-tracking branch 'origin/master' into updateazp

cmake/gtest.cmake

@@ -3,4 +3,20 @@ if (MSVC)
     #link dynamically too (default is /MT[d])
     option(gtest_force_shared_crt "Always use shared Visual C++ run-time DLL" ON)
 endif()
+
+set(GOOGLETEST_VERSION 1.10.0)
 add_subdirectory(vendor/gtest)
+
+find_package(absl)
+if (absl_FOUND)
+    cmake_policy(SET CMP0079 NEW)
+    target_compile_definitions(gtest PUBLIC GTEST_HAS_ABSL=1)
+    target_compile_definitions(gtest_main PUBLIC GTEST_HAS_ABSL=1)
+    target_link_libraries(gtest PRIVATE absl::algorithm
+                                        absl::base
+                                        absl::debugging
+                                        absl::numeric
+                                        absl::strings
+                                        absl::utility
+                                        absl::failure_signal_handler)
+endif()

doc/stages/filters.normal.rst

@@ -23,10 +23,10 @@ stages for more advanced filtering.
 The eigenvalue decomposition is performed using Eigen's
 `SelfAdjointEigenSolver <https://eigen.tuxfamily.org/dox/classEigen_1_1SelfAdjointEigenSolver.html>`_.
 
-Normals will be automatically flipped towards the viewpoint to be consistent. By
-default the viewpoint is located at the midpoint of the X and Y extents, and
-1000 units above the max Z value. Users can override any of the XYZ coordinates,
-or set them all to zero to effectively disable the normal flipping.
+Normals will be automatically flipped towards positive Z, unless the always_up_
+flag is set to `false`. Users can optionally set any of the XYZ coordinates to
+specify a custom viewpoint_ or set them all to zero to effectively disable the
+normal flipping.
 
 .. note::
 
@@ -36,6 +36,14 @@ or set them all to zero to effectively disable the normal flipping.
   regardless of the always_up_ flag. To disable all normal flipping, do not
   provide a viewpoint_ and set `always_up`_ to false.
 
+In addition to always_up_ and viewpoint_, users can run a refinement step (on
+by default) that propagates normals using a minimum spanning tree. The
+propagated normals can lead to much more consistent results across the dataset.
+
+.. note::
+
+  To disable normal propagation, users can set refine_ to `false`.
+
 .. embed::
 
 Example
@@ -73,3 +81,7 @@ _`viewpoint`
 _`always_up`
   A flag indicating whether or not normals should be inverted only when the Z
   component is negative. [Default: true]
+
+_`refine`
+  A flag indicating whether or not to reorient normals using minimum spanning
+  tree propagation. [Default: true]

filters/NormalFilter.cpp

@@ -1,5 +1,5 @@
 /******************************************************************************
- * Copyright (c) 2016-2017, Bradley J Chambers (brad.chambers@gmail.com)
+ * Copyright (c) 2016, 2017, 2019, 2020 Bradley J Chambers (brad.chambers@gmail.com)
  *
  * All rights reserved.
  *
@@ -32,6 +32,14 @@
  * OF SUCH DAMAGE.
  ****************************************************************************/
 
+// Normal refinement algorithm is presented in [1] and adapted to PDAL based on
+// the implementation provided in [2].
+//
+// [1] H. Hoppe, T.  DeRose, T. Duchamp, J. McDonald, and W. Stuetzle, "Surface
+// reconstruction from unorganized points," Computer Graphics, vol. 26. no. 2,
+// pp. 71-78, 1992.
+// [2] https://github.com/CloudCompare/CloudCompare.
+
 #include "NormalFilter.hpp"
 #include "private/Point.hpp"
 
@@ -47,12 +55,12 @@
 namespace pdal
 {
 
-static StaticPluginInfo const s_info
-{
-    "filters.normal",
-    "Normal Filter",
-    "http://pdal.io/stages/filters.normal.html"
-};
+using namespace Eigen;
+using namespace Dimension;
+
+static StaticPluginInfo const s_info{
+    "filters.normal", "Normal Filter",
+    "http://pdal.io/stages/filters.normal.html"};
 
 CREATE_STATIC_STAGE(NormalFilter, s_info)
 
@@ -61,15 +69,12 @@ struct NormalArgs
     int m_knn;
     filter::Point m_viewpoint;
     bool m_up;
+    bool m_refine;
 };
 
-NormalFilter::NormalFilter() : m_args(new NormalArgs)
-{}
-
-
-NormalFilter::~NormalFilter()
-{}
+NormalFilter::NormalFilter() : m_args(new NormalArgs), m_count(0) {}
 
+NormalFilter::~NormalFilter() {}
 
 std::string NormalFilter::getName() const
 {
@@ -79,16 +84,17 @@ std::string NormalFilter::getName() const
 void NormalFilter::addArgs(ProgramArgs& args)
 {
     args.add("knn", "k-Nearest Neighbors", m_args->m_knn, 8);
-    m_viewpointArg = &args.add("viewpoint",
-        "Viewpoint as WKT or GeoJSON", m_args->m_viewpoint);
+    m_viewpointArg = &args.add("viewpoint", "Viewpoint as WKT or GeoJSON",
+                               m_args->m_viewpoint);
     args.add("always_up", "Normals always oriented with positive Z?",
-        m_args->m_up, true);
+             m_args->m_up, true);
+    args.add("refine",
+             "Refine normals using minimum spanning tree propagation?",
+             m_args->m_refine, true);
 }
 
 void NormalFilter::addDimensions(PointLayoutPtr layout)
 {
-    using namespace Dimension;
-
     layout->registerDims(
         {Id::NormalX, Id::NormalY, Id::NormalZ, Id::Curvature});
 }
@@ -112,53 +118,174 @@ void NormalFilter::prepared(PointTableRef table)
             << "Viewpoint provided. Ignoring always_up = TRUE." << std::endl;
         m_args->m_up = false;
     }
+
+    // The query point is returned as a neighbor of itself, so we must increase
+    // k by one to get the desired number of neighbors.
+    ++m_args->m_knn;
 }
 
-void NormalFilter::filter(PointView& view)
+void NormalFilter::compute(PointView& view, KD3Index& kdi)
 {
-    KD3Index& kdi = view.build3dIndex();
-
-    for (PointId i = 0; i < view.size(); ++i)
+    log()->get(LogLevel::Debug) << "Computing normal vectors\n";
+    for (auto&& p : view)
     {
-        // find the k-nearest neighbors
-        auto ids = kdi.neighbors(i, m_args->m_knn);
-
-        // compute covariance of the neighborhood
-        auto B = computeCovariance(view, ids);
-
-        // perform the eigen decomposition
-        Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> solver(B);
-        if (solver.info() != Eigen::Success)
+        // Perform eigen decomposition of covariance matrix computed from
+        // neighborhood composed of k-nearest neighbors.
+        PointIdList neighbors = kdi.neighbors(p.pointId(), m_args->m_knn);
+        auto B = computeCovariance(view, neighbors);
+        SelfAdjointEigenSolver<Matrix3d> solver(B);
+        if (solver.info() != Success)
             throwError("Cannot perform eigen decomposition.");
+
+        // The curvature is computed as the ratio of the first (smallest)
+        // eigenvalue to the sum of all eigenvalues.
         auto eval = solver.eigenvalues();
-        Eigen::Vector3d normal = solver.eigenvectors().col(0);
+        double sum = eval[0] + eval[1] + eval[2];
+        double curvature = sum ? std::fabs(eval[0] / sum) : 0;
+
+        // The normal is defined by the eigenvector corresponding to the
+        // smallest eigenvalue.
+        Vector3d normal = solver.eigenvectors().col(0);
 
         if (m_viewpointArg->set())
         {
-            using namespace Dimension;
-
-            PointRef p = view.point(i);
-            Eigen::Vector3d vp(
-                (float)(m_args->m_viewpoint.x() - p.getFieldAs<double>(Id::X)),
-                (float)(m_args->m_viewpoint.y() - p.getFieldAs<double>(Id::Y)),
-                (float)(m_args->m_viewpoint.z() - p.getFieldAs<double>(Id::Z)));
+            // If a viewpoint has been specified, orient the normals to face the
+            // viewpoint by taking the dot product of the vector connecting the
+            // point with the viewpoint and the normal. Flip the normal, where
+            // the dot product is negative.
+            double dx = m_args->m_viewpoint.x() - p.getFieldAs<double>(Id::X);
+            double dy = m_args->m_viewpoint.y() - p.getFieldAs<double>(Id::Y);
+            double dz = m_args->m_viewpoint.z() - p.getFieldAs<double>(Id::Z);
+            Vector3d vp(dx, dy, dz);
             if (vp.dot(normal) < 0)
                 normal *= -1.0;
         }
         else if (m_args->m_up)
         {
+            // If normals are expected to be upward facing, invert them when the
+            // Z component is negative.
             if (normal[2] < 0)
                 normal *= -1.0;
         }
 
-        view.setField(Dimension::Id::NormalX, i, normal[0]);
-        view.setField(Dimension::Id::NormalY, i, normal[1]);
-        view.setField(Dimension::Id::NormalZ, i, normal[2]);
+        // Set the computed normal and curvature dimensions.
+        p.setField(Id::NormalX, normal[0]);
+        p.setField(Id::NormalY, normal[1]);
+        p.setField(Id::NormalZ, normal[2]);
+        p.setField(Id::Curvature, curvature);
+    }
+}
 
-        double sum = eval[0] + eval[1] + eval[2];
-        view.setField(Dimension::Id::Curvature, i,
-                      sum ? std::fabs(eval[0] / sum) : 0);
+void NormalFilter::update(
+    PointView& view, KD3Index& kdi, std::vector<bool> inMST,
+    std::priority_queue<Edge, EdgeList, CompareEdgeWeight> edge_queue,
+    PointId updateIdx)
+{
+    // Add the current PointId to the minimum spanning tree.
+    inMST[updateIdx] = true;
+    ++m_count;
+
+    // Consider neighbors of the newly added PointId, adding them to
+    // the edge queue if they are not already part of the minimum
+    // spanning tree. The first neighbor is the query point which is
+    // already part of the minimum spanning tree and can safely be
+    // skipped.
+    PointIdList neighbors = kdi.neighbors(updateIdx, m_args->m_knn);
+    neighbors.erase(neighbors.begin());
+    PointRef p = view.point(updateIdx);
+    Vector3d N0(p.getFieldAs<double>(Id::NormalX),
+                p.getFieldAs<double>(Id::NormalY),
+                p.getFieldAs<double>(Id::NormalZ));
+    for (PointId const& neighborIdx : neighbors)
+    {
+        if (!inMST[neighborIdx])
+        {
+            PointRef q = view.point(neighborIdx);
+            Vector3d N1(q.getFieldAs<double>(Id::NormalX),
+                        q.getFieldAs<double>(Id::NormalY),
+                        q.getFieldAs<double>(Id::NormalZ));
+            double weight = 1.0 - std::fabs(N0.dot(N1));
+            edge_queue.emplace(updateIdx, neighborIdx, weight);
+        }
+    }
+}
+
+void NormalFilter::refine(PointView& view, KD3Index& kdi)
+{
+    log()->get(LogLevel::Debug)
+        << "Refining normals using minimum spanning tree\n";
+
+    std::priority_queue<Edge, EdgeList, CompareEdgeWeight> edge_queue;
+    std::vector<bool> inMST(view.size(), false);
+    PointId nextIdx(0);
+    while (m_count < view.size())
+    {
+        // Find the PointId of the next point not currently part of the minimum
+        // spanning tree.
+        while (inMST[nextIdx])
+            ++nextIdx;
+
+        update(view, kdi, inMST, edge_queue, nextIdx);
+
+        // Iterate on the edge queue until empty (or all points have been added
+        // to the minimum spanning tree).
+        while (!edge_queue.empty() && (m_count < view.size()))
+        {
+            // Retrieve the edge with the smallest weight.
+            Edge edge(edge_queue.top());
+            edge_queue.pop();
+
+            // Record the PointId and normal of the PointId (if one exists)
+            // that is not already in the minimum spanning tree.
+            PointId newIdx(0);
+            Vector3d normal;
+            PointRef p = view.point(edge.m_v0);
+            Vector3d N0(p.getFieldAs<double>(Id::NormalX),
+                        p.getFieldAs<double>(Id::NormalY),
+                        p.getFieldAs<double>(Id::NormalZ));
+            PointRef q = view.point(edge.m_v1);
+            Vector3d N1(q.getFieldAs<double>(Id::NormalX),
+                        q.getFieldAs<double>(Id::NormalY),
+                        q.getFieldAs<double>(Id::NormalZ));
+            if (!inMST[edge.m_v0])
+            {
+                newIdx = edge.m_v0;
+                normal = N0;
+            }
+            else if (!inMST[edge.m_v1])
+            {
+                newIdx = edge.m_v1;
+                normal = N1;
+            }
+            else
+                continue;
+
+            // Where the dot product of the normals is less than 0, invert the
+            // normal of the selected PointId.
+            if (N0.dot(N1) < 0)
+            {
+                normal *= -1;
+                view.setField(Id::NormalX, newIdx, normal(0));
+                view.setField(Id::NormalY, newIdx, normal(1));
+                view.setField(Id::NormalZ, newIdx, normal(2));
+            }
+
+            update(view, kdi, inMST, edge_queue, newIdx);
+        }
     }
 }
 
+void NormalFilter::filter(PointView& view)
+{
+    KD3Index& kdi = view.build3dIndex();
+
+    // Compute the normal/curvature and optionally orient toward viewpoint or
+    // positive Z.
+    compute(view, kdi);
+
+    // If requested, refine normals through minimum spanning tree propagation.
+    if (m_args->m_refine)
+        refine(view, kdi);
+}
+
 } // namespace pdal

filters/NormalFilter.hpp

@@ -1,5 +1,5 @@
 /******************************************************************************
- * Copyright (c) 2016-2017, Bradley J Chambers (brad.chambers@gmail.com)
+ * Copyright (c) 2016, 2017, 2020 Bradley J Chambers (brad.chambers@gmail.com)
  *
  * All rights reserved.
  *
@@ -48,6 +48,28 @@ class PointLayout;
 class PointView;
 struct NormalArgs;
 
+struct Edge
+{
+    PointId m_v0;
+    PointId m_v1;
+    double m_weight;
+
+    Edge(PointId i, PointId j, double weight)
+        : m_v0(i), m_v1(j), m_weight(weight)
+    {
+    }
+};
+
+struct CompareEdgeWeight
+{
+    bool operator()(Edge const& lhs, Edge const& rhs)
+    {
+        return lhs.m_weight > rhs.m_weight;
+    }
+};
+
+typedef std::vector<Edge> EdgeList;
+
 class PDAL_DLL NormalFilter : public Filter
 {
 public:
@@ -63,8 +85,16 @@ class PDAL_DLL NormalFilter : public Filter
 
 private:
     std::unique_ptr<NormalArgs> m_args;
+    point_count_t m_count;
     Arg* m_viewpointArg;
 
+    void compute(PointView& view, KD3Index& kdi);
+    void refine(PointView& view, KD3Index& kdi);
+    void
+    update(PointView& view, KD3Index& kdi, std::vector<bool> inMST,
+           std::priority_queue<Edge, EdgeList, CompareEdgeWeight> edge_queue,
+           PointId updateIdx);
+
     virtual void addArgs(ProgramArgs& args);
     virtual void addDimensions(PointLayoutPtr layout);
     virtual void prepared(PointTableRef table);