diff --git a/pdal/EigenUtils.cpp b/pdal/EigenUtils.cpp index 2dae5d40a3..b69704debe 100644 --- a/pdal/EigenUtils.cpp +++ b/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); diff --git a/pdal/EigenUtils.hpp b/pdal/EigenUtils.hpp index 0120b957ca..e40142c2be 100644 --- a/pdal/EigenUtils.hpp +++ b/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 +PDAL_DLL Derived padMatrix(const Eigen::MatrixBase& 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& 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 +PDAL_DLL Derived dilate(const Eigen::MatrixBase& 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 +PDAL_DLL Derived erode(const Eigen::MatrixBase& 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. diff --git a/test/unit/EigenTest.cpp b/test/unit/EigenTest.cpp index cedc174974..c10154450b 100644 --- a/test/unit/EigenTest.cpp +++ b/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,8 +132,8 @@ 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); @@ -141,9 +141,9 @@ TEST(EigenTest, ComputeValues) 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)); +}