Merge pull request #1488 from PDAL/slope-update

Improve quality of derivative writer outputs

io/DerivativeWriter.cpp

@@ -57,6 +57,8 @@ void DerivativeWriter::addArgs(ProgramArgs& args)
     args.add("filename", "Output filename", m_filename).setPositional();
     args.add("edge_length", "Edge length", m_edgeLength, 15.0);
     args.add("primitive_type", "Primitive type", m_primTypesSpec, {"slope_d8"});
+    args.add("slope_units", "Slope units (degrees/percent)", m_slopeUnit,
+             "percent");
     args.add("altitude", "Illumination altitude (degrees)", m_illumAltDeg,
              45.0);
     args.add("azimuth", "Illumination azimuth (degrees)", m_illumAzDeg, 315.0);
@@ -125,19 +127,16 @@ void DerivativeWriter::write(const PointViewPtr data)
     size_t rows = ((bounds.maxy - bounds.miny) / m_edgeLength) + 1;
 
     // Begin by creating a DSM of max elevations per XY cell.
-    MatrixXd DSM = createMaxMatrix(*data.get(), rows, cols, m_edgeLength,
-                                   bounds);
-
-    // Continue by cleaning the DSM.
-    MatrixXd cleanedDSM = cleanDSM(DSM);
+    MatrixXd DSM = createMaxMatrix2(*data.get(), rows, cols, m_edgeLength,
+                                    bounds);
 
     // We will pad the edges by 1 cell, though according to some texts we should
     // simply compute forward- or backward-difference as opposed to centered
     // difference at these points.
-    MatrixXd paddedDSM = padMatrix(cleanedDSM, 1);
+    MatrixXd paddedDSM = padMatrix(DSM, 1);
 
     // Prepare the out matrix.
-    MatrixXd out(cleanedDSM.rows(), cleanedDSM.cols());
+    MatrixXd out(DSM.rows(), DSM.cols());
     out.setConstant(std::numeric_limits<double>::quiet_NaN());
 
     for (TypeOutput& to : m_primitiveTypes)
@@ -151,9 +150,17 @@ void DerivativeWriter::write(const PointViewPtr data)
                     continue;
                 Matrix3d block = paddedDSM.block(r-1, c-1, 3, 3);
                 if (to.m_type == SLOPE_D8)
+                {
                     out(r-1, c-1) = computeSlopeD8(block, m_edgeLength);
+                    if (Utils::iequals(m_slopeUnit, "degrees"))
+                        out(r-1, c-1) = percentSlopeToDegrees(out(r-1, c-1));
+                }
                 if (to.m_type == SLOPE_FD)
+                {
                     out(r-1, c-1) = computeSlopeFD(block, m_edgeLength);
+                    if (Utils::iequals(m_slopeUnit, "degrees"))
+                        out(r-1, c-1) = percentSlopeToDegrees(out(r-1, c-1));
+                }
                 if (to.m_type == ASPECT_D8)
                     out(r-1, c-1) = computeAspectD8(block, m_edgeLength);
                 if (to.m_type == ASPECT_FD)

io/DerivativeWriter.hpp

@@ -90,12 +90,20 @@ class PDAL_DLL DerivativeWriter : public Writer
 
     std::string m_filename;
     std::string m_driver;
+    std::string m_slopeUnit;
     double m_edgeLength;
     double m_illumAltDeg;
     double m_illumAzDeg;
     StringList m_primTypesSpec;
     std::vector<TypeOutput> m_primitiveTypes;
 
+    const double c_PI = std::acos(-1.0);
+    const double c_rad2deg = 180.0 / c_PI;
+    double percentSlopeToDegrees(double slope)
+    {
+        return std::atan(slope / 100.0) * c_rad2deg;
+    }
+
     DerivativeWriter& operator=(const DerivativeWriter&); // not implemented
     DerivativeWriter(const DerivativeWriter&); // not implemented
 };

pdal/EigenUtils.cpp

@@ -34,6 +34,7 @@
 
 #include <pdal/EigenUtils.hpp>
 #include <pdal/GDALUtils.hpp>
+#include <pdal/KDIndex.hpp>
 #include <pdal/PointView.hpp>
 #include <pdal/SpatialReference.hpp>
 #include <pdal/util/Bounds.hpp>
@@ -255,6 +256,41 @@ Eigen::MatrixXd createMaxMatrix(PointView& view, int rows, int cols,
     return ZImax;
 }
 
+Eigen::MatrixXd createMaxMatrix2(PointView& view, int rows, int cols,
+                                 double cell_size, BOX2D bounds)
+{
+    using namespace Dimension;
+    using namespace Eigen;
+
+    KD2Index kdi(view);
+    kdi.build();
+
+    MatrixXd ZImax(rows, cols);
+    ZImax.setConstant(std::numeric_limits<double>::lowest());
+
+    // for each grid center, search PointView for neighbors, and find max of those
+    for (int c = 0; c < cols; ++c)
+    {
+        double x = bounds.minx + (c + 0.5) * cell_size;
+
+        for (int r = 0; r < rows; ++r)
+        {
+            double y = bounds.miny + (r + 0.5) * cell_size;
+
+            auto neighbors = kdi.radius(x, y, cell_size * std::sqrt(2.0));
+
+            for (auto const& n : neighbors)
+            {
+                PointRef point = view.point(n);
+                if (point.getFieldAs<double>(Id::Z) > ZImax(r, c))
+                    ZImax(r, c) = point.getFieldAs<double>(Id::Z);
+            }
+        }
+    }
+
+    return ZImax;
+}
+
 Eigen::MatrixXd extendedLocalMinimum(PointView& view, int rows, int cols,
                                      double cell_size, BOX2D bounds)
 {

pdal/EigenUtils.hpp

@@ -250,6 +250,9 @@ PDAL_DLL uint8_t computeRank(PointView& view, std::vector<PointId> ids,
 PDAL_DLL Eigen::MatrixXd createMaxMatrix(PointView& view, int rows, int cols,
         double cell_size, BOX2D bounds);
 
+PDAL_DLL Eigen::MatrixXd createMaxMatrix2(PointView& view, int rows, int cols,
+        double cell_size, BOX2D bounds);
+
 /**
   Create matrix of minimum Z values.
 
@@ -738,12 +741,14 @@ PDAL_DLL double computeSlopeD8(const Eigen::MatrixBase<Derived>& data,
     submatrix.setConstant(data(1, 1));
     submatrix -= data;
     submatrix /= spacing;
-    submatrix(0, 1) /= std::sqrt(2.0);
-    submatrix(1, 0) /= std::sqrt(2.0);
-    submatrix(1, 2) /= std::sqrt(2.0);
-    submatrix(2, 1) /= std::sqrt(2.0);
-
-    // find max and convert to degrees
+    submatrix(0, 0) /= std::sqrt(2.0);
+    submatrix(0, 2) /= std::sqrt(2.0);
+    submatrix(2, 0) /= std::sqrt(2.0);
+    submatrix(2, 2) /= std::sqrt(2.0);
+
+    // Why not just use Eigen's maxCoeff reduction to find the max? Well, as it
+    // turns out, if there is a chance that we will have NaN's then maxCoeff
+    // has no way to ignore the NaN.
     double maxval = std::numeric_limits<double>::lowest();
     for (int i = 0; i < submatrix.size(); ++i)
     {
@@ -791,7 +796,7 @@ template <typename Derived>
 PDAL_DLL Derived dilate(const Eigen::MatrixBase<Derived>& A, int radius)
 {
     Derived B = Derived::Constant(A.rows(), A.cols(), 0);
-    
+
     int length = 2 * radius + 1;
     bool match_flag;
     for (int c = 0; c < A.cols(); ++c)
@@ -826,7 +831,7 @@ PDAL_DLL Derived dilate(const Eigen::MatrixBase<Derived>& A, int radius)
             B(r, c) = (match_flag) ? 1 : 0;
         }
     }
-    
+
     return B;
 }
 
@@ -844,7 +849,7 @@ template <typename Derived>
 PDAL_DLL Derived erode(const Eigen::MatrixBase<Derived>& A, int radius)
 {
     Derived B = Derived::Constant(A.rows(), A.cols(), 1);
-    
+
     int length = 2 * radius + 1;
     bool mismatch_flag;
     for (int c = 0; c < A.cols(); ++c)
@@ -885,7 +890,7 @@ PDAL_DLL Derived erode(const Eigen::MatrixBase<Derived>& A, int radius)
             B(r, c) = (mismatch_flag) ? 0 : 1;
         }
     }
-    
+
     return B;
 }
 

test/unit/EigenTest.cpp

@@ -115,7 +115,7 @@ TEST(EigenTest, ComputeValues)
     EXPECT_NEAR(0.8236099869, slope, 0.0001);
 
     double sd8 = eigen::computeSlopeD8(A, spacing);
-    EXPECT_NEAR(48.0378, sd8, 0.0001);
+    EXPECT_NEAR(67.9357, sd8, 0.0001);
 
     double sfd = eigen::computeSlopeFD(A, spacing);
     EXPECT_NEAR(210.1961, sfd, 0.0001);