Update EigenUtils

* Add dilate and erode functions to EigenUtils
* Make padMatrix a templated function

pdal/EigenUtils.cpp

@@ -428,24 +428,6 @@ Eigen::MatrixXd matrixOpen(Eigen::MatrixXd data, int radius)
     return maxZ.block(radius, radius, data.rows(), data.cols());
 }
 
-Eigen::MatrixXd padMatrix(Eigen::MatrixXd d, int r)
-{
-    using namespace Eigen;
-
-    MatrixXd out = MatrixXd::Zero(d.rows()+2*r, d.cols()+2*r);
-    out.block(r, r, d.rows(), d.cols()) = d;
-    out.block(r, 0, d.rows(), r) =
-        d.block(0, 0, d.rows(), r).rowwise().reverse();
-    out.block(r, d.cols()+r, d.rows(), r) =
-        d.block(0, d.cols()-r, d.rows(), r).rowwise().reverse();
-    out.block(0, 0, r, out.cols()) =
-        out.block(r, 0, r, out.cols()).colwise().reverse();
-    out.block(d.rows()+r, 0, r, out.cols()) =
-        out.block(out.rows()-r-1, 0, r, out.cols()).colwise().reverse();
-
-    return out;
-}
-
 Eigen::MatrixXd pointViewToEigen(const PointView& view)
 {
     Eigen::MatrixXd matrix(view.size(), 3);

pdal/EigenUtils.hpp

@@ -322,7 +322,24 @@ PDAL_DLL Eigen::MatrixXd matrixOpen(Eigen::MatrixXd data, int radius);
   \param r the radius of the padding.
   \return the padded matrix.
 */
-PDAL_DLL Eigen::MatrixXd padMatrix(Eigen::MatrixXd d, int r);
+template <typename Derived>
+PDAL_DLL Derived padMatrix(const Eigen::MatrixBase<Derived>& d, int r)
+{
+    using namespace Eigen;
+
+    Derived out = Derived::Zero(d.rows()+2*r, d.cols()+2*r);
+    out.block(r, r, d.rows(), d.cols()) = d;
+    out.block(r, 0, d.rows(), r) =
+        d.block(0, 0, d.rows(), r).rowwise().reverse();
+    out.block(r, d.cols()+r, d.rows(), r) =
+        d.block(0, d.cols()-r, d.rows(), r).rowwise().reverse();
+    out.block(0, 0, r, out.cols()) =
+        out.block(r, 0, r, out.cols()).colwise().reverse();
+    out.block(d.rows()+r, 0, r, out.cols()) =
+        out.block(out.rows()-r-1, 0, r, out.cols()).colwise().reverse();
+
+    return out;
+}
 
 /**
   Converts a PointView into an Eigen::MatrixXd.
@@ -760,6 +777,118 @@ PDAL_DLL double computeSlopeFD(const Eigen::MatrixBase<Derived>& data,
     return 100.0 * std::sqrt(p);
 }
 
+/**
+  Perform a morphological dilation of the input matrix.
+
+  Performs a morphological dilation of the input matrix using a circular
+  structuring element of given radius.
+
+  \param data the input matrix.
+  \param radius the radius of the circular structuring element.
+  \return the morphological dilation of the input matrix.
+*/
+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)
+    {
+        for (int r = 0; r < A.rows(); ++r)
+        {
+            match_flag = false;
+            for (int k = 0; k < length; ++k)
+            {
+                if (match_flag)
+                    break;
+                int cdiff = k-radius;
+                int cpos = c+cdiff;
+                if (cpos < 0 || cpos >= A.cols())
+                    continue;
+                for (int l = 0; l < length; ++l)
+                {
+                    int rdiff = l-radius;
+                    int rpos = r+rdiff;
+                    if (rpos < 0 || rpos >= A.rows())
+                        continue;
+                    if ((cdiff*cdiff+rdiff*rdiff) > radius*radius)
+                        continue;
+                    if (A(rpos, cpos) == 1)
+                    {
+                        match_flag = true;
+                        break;
+                    }
+                }
+            }
+            // Assign value according to match flag
+            B(r, c) = (match_flag) ? 1 : 0;
+        }
+    }
+
+    return B;
+}
+
+/**
+  Perform a morphological erosion of the input matrix.
+
+  Performs a morphological erosion of the input matrix using a circular
+  structuring element of given radius.
+
+  \param data the input matrix.
+  \param radius the radius of the circular structuring element.
+  \return the morphological erosion of the input matrix.
+*/
+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)
+    {
+        for (int r = 0; r < A.rows(); ++r)
+        {
+            if (A(r, c) == 0)
+            {
+                B(r, c) = 0;
+                continue;
+            }
+            mismatch_flag = false;
+            for (int k = 0; k < length; k++)
+            {
+                if (mismatch_flag)
+                    break;
+                int cdiff = k-radius;
+                int cpos = c+cdiff;
+                if (cpos < 0 || cpos >= A.cols())
+                    continue;
+                for (int l = 0; l < length; l++)
+                {
+                    int rdiff = l-radius;
+                    int rpos = r+rdiff;
+                    if (rpos < 0 || rpos >= A.rows())
+                        continue;
+                    if ((cdiff*cdiff+rdiff*rdiff) > radius*radius)
+                        continue;
+                    if (A(rpos, cpos) == 0)
+                    {
+                        B(r, c) = 0;
+                        mismatch_flag = true;
+                        break;
+                    }
+                }
+            }
+            // Assign value according to mismatch flag
+            B(r, c) = (mismatch_flag) ? 0 : 1;
+        }
+    }
+
+    return B;
+}
+
 /**
   Thin Plate Spline interpolation.
 

test/unit/EigenTest.cpp

@@ -66,9 +66,9 @@ TEST(EigenTest, ComputeValues)
     using namespace Eigen;
 
     Matrix3d A;
-    A << 1.8339, 0.3188, 0.3426,
-    -2.2588, -1.3077, 3.5784,
-    0.8622, -0.4336, 2.7694;
+    A << 1.8339,  0.3188, 0.3426,
+        -2.2588, -1.3077, 3.5784,
+         0.8622, -0.4336, 2.7694;
 
     double spacing(1.4);
 
@@ -132,18 +132,18 @@ TEST(EigenTest, ComputeValues)
 
     Matrix3d gx;
     gx << -1.5151, -0.7457, 0.0238,
-    0.9511, 2.9186, 4.8861,
-    -1.2958, 0.9536, 3.2030;
+           0.9511,  2.9186, 4.8861,
+          -1.2958,  0.9536, 3.2030;
 
     for (size_t i = 0; i < 9; ++i)
         EXPECT_NEAR(gx(i), out(i), 0.0001);
 
     MatrixXd out2 = eigen::gradY(A);
 
     Matrix3d gy;
-    gy << -4.0927, -1.6265, 3.2358,
-    -0.4859, -0.3762, 1.2134,
-    3.1210, 0.8741, -0.8090;
+    gy << -4.0927, -1.6265,  3.2358,
+          -0.4859, -0.3762,  1.2134,
+           3.1210,  0.8741, -0.8090;
 
     for (size_t i = 0; i < 9; ++i)
         EXPECT_NEAR(gy(i), out2(i), 0.0001);
@@ -200,10 +200,10 @@ TEST(EigenTest, Padding)
 {
     using namespace Eigen;
 
-    Matrix3d A;
-    A << 1.8339, 0.3188, 0.3426,
-    -2.2588, -1.3077, 3.5784,
-    0.8622, -0.4336, 2.7694;
+    MatrixXd A(3, 3);
+    A << 1.8339,  0.3188, 0.3426,
+        -2.2588, -1.3077, 3.5784,
+         0.8622, -0.4336, 2.7694;
 
     MatrixXd B = eigen::padMatrix(A, 1);
 
@@ -215,3 +215,33 @@ TEST(EigenTest, Padding)
     EXPECT_EQ(3.5784, B(2, 4));
     EXPECT_EQ(-0.4336, B(4, 2));
 }
+
+TEST(EigenTest, Dilate)
+{
+    using namespace Eigen;
+
+    MatrixXd C(5, 5);
+    C << 0, 0, 0, 0, 0,
+         0, 1, 0, 0, 0,
+         0, 1, 1, 0, 0,
+         0, 0, 1, 1, 0,
+         0, 0, 0, 0, 0;
+
+    MatrixXd D = eigen::dilate(C, 1);
+
+    EXPECT_EQ(0, D(0, 0));
+    EXPECT_EQ(1, D(1, 0));
+    EXPECT_EQ(1, D(0, 1));
+
+    MatrixXd E(5, 5);
+    E << 0, 0, 0, 0, 0,
+         0, 1, 1, 1, 0,
+         0, 1, 1, 1, 0,
+         0, 1, 1, 1, 0,
+         0, 0, 0, 0, 0;
+
+    MatrixXd F = eigen::erode(E, 1);
+
+    EXPECT_EQ(0, F(1, 3));
+    EXPECT_EQ(1, F(2, 2));
+}