Merge pull request #1639 from PDAL/issue/1638-normal-flip

Provide for normal flipping

doc/stages/filters.normal.rst

@@ -23,6 +23,19 @@ The eigenvalue decomposition is performed using Eigen's
 ``SelfAdjointEigenSolver``. For more information see
 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.
+
+.. note::
+
+  By default, the Normal filter will invert normals such that they are always
+  pointed "up" (positive Z). If the user provides a ``viewpoint``, normals will
+  instead be inverted such that they are oriented towards the viewpoint,
+  regardless of the ``always_up`` flag. To disable all normal flipping, do not
+  provide a ``viewpoint`` and set ``always_up=false``.
+
 .. embed::
 
 Example
@@ -54,3 +67,11 @@ Options
 
 knn
   The number of k-nearest neighbors. [Default: **8**]
+
+viewpoint
+  A single WKT or GeoJSON 3D point. Normals will be inverted such that they are
+  all oriented towards the viewpoint.
+
+always_up
+  A flag indicating whether or not normals should be inverted only when the Z
+  component is negative. [Default: **true**]

filters/CropFilter.cpp

@@ -40,7 +40,8 @@
 #include <pdal/Polygon.hpp>
 #include <pdal/pdal_macros.hpp>
 #include <pdal/util/ProgramArgs.hpp>
-#include <filters/private/crop/Point.hpp>
+
+#include "private/Point.hpp"
 
 #include <sstream>
 #include <cstdarg>
@@ -251,7 +252,7 @@ void CropFilter::crop(const Polygon& g, PointView& input, PointView& output)
 }
 
 
-bool CropFilter::crop(const PointRef& point, const cropfilter::Point& center)
+bool CropFilter::crop(const PointRef& point, const filter::Point& center)
 {
     double x = point.getFieldAs<double>(Dimension::Id::X);
     double y = point.getFieldAs<double>(Dimension::Id::Y);
@@ -275,7 +276,7 @@ bool CropFilter::crop(const PointRef& point, const cropfilter::Point& center)
 }
 
 
-void CropFilter::crop(const cropfilter::Point& center, PointView& input,
+void CropFilter::crop(const filter::Point& center, PointView& input,
     PointView& output)
 {
     PointRef point = input.point(0);

filters/CropFilter.hpp

@@ -38,16 +38,13 @@
 #include <pdal/Polygon.hpp>
 #include <pdal/plugin.hpp>
 
+#include "private/Point.hpp"
+
 extern "C" int32_t CropFilter_ExitFunc();
 extern "C" PF_ExitFunc CropFilter_InitPlugin();
 
 namespace pdal
 {
-namespace cropfilter
-{
-    class Point;
-};
-
 
 class ProgramArgs;
 
@@ -69,7 +66,7 @@ class PDAL_DLL CropFilter : public Filter
     SpatialReference m_assignedSrs;
     double m_distance;
     double m_distance2;
-    std::vector<cropfilter::Point> m_centers;
+    std::vector<filter::Point> m_centers;
     std::vector<Polygon> m_geoms;
     std::vector<BOX2D> m_boxes;
 
@@ -83,8 +80,8 @@ class PDAL_DLL CropFilter : public Filter
     void crop(const BOX2D& box, PointView& input, PointView& output);
     bool crop(const PointRef& point, const Polygon& g);
     void crop(const Polygon& g, PointView& input, PointView& output);
-    bool crop(const PointRef& point, const cropfilter::Point& center);
-    void crop(const cropfilter::Point& center, PointView& input,
+    bool crop(const PointRef& point, const filter::Point& center);
+    void crop(const filter::Point& center, PointView& input,
         PointView& output);
     void transform(const SpatialReference& srs);
 
@@ -93,4 +90,3 @@ class PDAL_DLL CropFilter : public Filter
 };
 
 } // namespace pdal
-

filters/NormalFilter.cpp

@@ -1,36 +1,36 @@
 /******************************************************************************
-* Copyright (c) 2016, Bradley J Chambers (brad.chambers@gmail.com)
-*
-* 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.
-****************************************************************************/
+ * Copyright (c) 2016-2017, Bradley J Chambers (brad.chambers@gmail.com)
+ *
+ * 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 "NormalFilter.hpp"
 
@@ -58,28 +58,39 @@ std::string NormalFilter::getName() const
     return s_info.name;
 }
 
-
 void NormalFilter::addArgs(ProgramArgs& args)
 {
     args.add("knn", "k-Nearest Neighbors", m_knn, 8);
+    m_viewpointArg =
+        &args.add("viewpoint", "Viewpoint as WKT or GeoJSON", m_viewpoint);
+    args.add("always_up", "Normals always oriented with positive Z?", m_up,
+             true);
 }
 
-
 void NormalFilter::addDimensions(PointLayoutPtr layout)
 {
     using namespace Dimension;
 
-    layout->registerDims({Id::NormalX, Id::NormalY, Id::NormalZ,
-        Id::Curvature});
+    layout->registerDims(
+        {Id::NormalX, Id::NormalY, Id::NormalZ, Id::Curvature});
 }
 
-
+// public method to access filter, used by GreedyProjection and Poisson filters
 void NormalFilter::doFilter(PointView& view)
 {
     m_knn = 8;
     filter(view);
 }
 
+void NormalFilter::prepared(PointTableRef table)
+{
+    if (m_up && m_viewpointArg->set())
+    {
+        log()->get(LogLevel::Warning)
+            << "Viewpoint provided. Ignoring always_up = TRUE." << std::endl;
+        m_up = false;
+    }
+}
 
 void NormalFilter::filter(PointView& view)
 {
@@ -98,15 +109,31 @@ void NormalFilter::filter(PointView& view)
         if (solver.info() != Eigen::Success)
             throwError("Cannot perform eigen decomposition.");
         auto eval = solver.eigenvalues();
-        auto evec = solver.eigenvectors().col(0);
-
-        view.setField(Dimension::Id::NormalX, i, evec[0]);
-        view.setField(Dimension::Id::NormalY, i, evec[1]);
-        view.setField(Dimension::Id::NormalZ, i, evec[2]);
+        Eigen::Vector3f normal = solver.eigenvectors().col(0);
+
+        if (m_viewpointArg->set())
+        {
+            PointRef p = view.point(i);
+            Eigen::Vector3f vp(
+                m_viewpoint.x - p.getFieldAs<double>(Dimension::Id::X),
+                m_viewpoint.y - p.getFieldAs<double>(Dimension::Id::Y),
+                m_viewpoint.z - p.getFieldAs<double>(Dimension::Id::Z));
+            if (vp.dot(normal) < 0)
+                normal *= -1.0;
+        }
+        else if (m_up)
+        {
+            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]);
 
         double sum = eval[0] + eval[1] + eval[2];
         view.setField(Dimension::Id::Curvature, i,
-            sum ? std::fabs(eval[0] / sum) : 0);
+                      sum ? std::fabs(eval[0] / sum) : 0);
     }
 }
 

filters/NormalFilter.hpp

@@ -1,41 +1,44 @@
 /******************************************************************************
-* Copyright (c) 2016, Bradley J Chambers (brad.chambers@gmail.com)
-*
-* 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.
-****************************************************************************/
+ * Copyright (c) 2016-2017, Bradley J Chambers (brad.chambers@gmail.com)
+ *
+ * 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.
+ ****************************************************************************/
 
 #pragma once
 
 #include <pdal/Filter.hpp>
 #include <pdal/plugin.hpp>
+#include <pdal/util/ProgramArgs.hpp>
+
+#include "private/Point.hpp"
 
 #include <cstdint>
 #include <memory>
@@ -55,23 +58,27 @@ class PDAL_DLL NormalFilter : public Filter
 {
 public:
     NormalFilter() : Filter()
-    {}
+    {
+    }
     NormalFilter& operator=(const NormalFilter&) = delete;
     NormalFilter(const NormalFilter&) = delete;
 
     void doFilter(PointView& view);
 
-    static void * create();
-    static int32_t destroy(void *);
+    static void* create();
+    static int32_t destroy(void*);
     std::string getName() const;
 
 private:
     int m_knn;
+    filter::Point m_viewpoint;
+    Arg* m_viewpointArg;
+    bool m_up;
 
     virtual void addArgs(ProgramArgs& args);
     virtual void addDimensions(PointLayoutPtr layout);
+    virtual void prepared(PointTableRef table);
     virtual void filter(PointView& view);
-
 };
 
 } // namespace pdal

filters/private/crop/Point.cpp → filters/private/Point.cpp

@@ -37,6 +37,9 @@
 namespace pdal
 {
 
+namespace filter
+{
+
 namespace
 {
 
@@ -45,9 +48,6 @@ const double HIGHEST = (std::numeric_limits<double>::max)();
 
 }
 
-namespace cropfilter
-{
-
 Point::Point()
     : Geometry()
     , x(LOWEST)
@@ -114,7 +114,6 @@ bool Point::is3d() const
     return (z != LOWEST);
 }
 
-} //namespace cropfilter
-
-} //namespace pdal
+} // namespace filter
 
+} // namespace pdal