Add some logging and error handling for covariancefeature edge cases (PDAL#4707)

* Add some logging and error handling for covariancefeature edge cases

* Log empty covariance in several filters

* bump fixed builder, add geotiff package to current

---------

Co-authored-by: Howard Butler <[email protected]>

Author: chambbj

filters/ApproximateCoplanarFilter.cpp

@@ -89,6 +89,17 @@ void ApproximateCoplanarFilter::filter(PointView& view)
         // compute covariance of the neighborhood
         Matrix3d B = math::computeCovariance(view, ids);
 
+        // Check if the covariance matrix is all zeros
+        if (B.isZero())
+        {
+            log()->get(LogLevel::Info)
+                << "Skipping point " << p.pointId()
+                << ". Covariance matrix is all zeros. This suggests a large "
+                   "number of redundant points. Consider using filters.sample "
+                   "with a small radius to remove redundant points.\n";
+            continue;
+        }
+
         // perform the eigen decomposition
         Eigen::SelfAdjointEigenSolver<Matrix3d> solver(B);
         if (solver.info() != Eigen::Success)

filters/CovarianceFeaturesFilter.cpp

@@ -181,6 +181,9 @@ void CovarianceFeaturesFilter::filter(PointView& view)
     if (npoints % m_threads) chunk_size++;
     std::vector<std::thread> threadList(m_threads);
 
+    log()->get(LogLevel::Debug) << "Processing " << npoints << " points in "
+                                << m_threads << " threads.\n";
+
     for(int t = 0; t < m_threads; t++)
     {
         threadList[t] = std::thread(
@@ -213,8 +216,12 @@ void CovarianceFeaturesFilter::setDimensionality(PointView &view, const PointId
     else if (m_radiusArg->set())
     {
         ids = kdi.radius(p, m_radius);
-        if (ids.size() < (size_t)m_minK)
+        if (ids.size() < (size_t)m_minK) {
+            log()->get(LogLevel::Info)
+                << "Skipping point " << id << ". Found " << ids.size()
+                << " neighbors but required " << m_minK << ".\n";
             return;
+        }
     }
     else
     {
@@ -224,6 +231,17 @@ void CovarianceFeaturesFilter::setDimensionality(PointView &view, const PointId
     // compute covariance of the neighborhood
     auto B = math::computeCovariance(view, ids);
 
+    // Check if the covariance matrix is all zeros
+    if (B.isZero())
+    {
+        log()->get(LogLevel::Info)
+            << "Skipping point " << id
+            << ". Covariance matrix is all zeros. This suggests a large number "
+               "of redundant points. Consider using filters.sample with a "
+               "small radius to remove redundant points.\n";
+        return;
+    }
+
     // perform the eigen decomposition
     SelfAdjointEigenSolver<Matrix3d> solver(B);
     if (solver.info() != Success)

filters/EigenvaluesFilter.cpp

@@ -127,8 +127,12 @@ void EigenvaluesFilter::filter(PointView& view)
             // if insufficient number of neighbors, eigen solver will fail
             // anyway, it may be okay to silently return without setting any of
             // the computed features?
-            if (ids.size() < (size_t)m_args->m_minK)
+            if (ids.size() < (size_t)m_args->m_minK) {
+                log()->get(LogLevel::Info)
+                    << "Skipping point " << p.pointId() << ". Found " << ids.size()
+                    << " neighbors but required " << m_args->m_minK << ".\n";
                 continue;
+            }
         }
         else
         {
@@ -138,6 +142,17 @@ void EigenvaluesFilter::filter(PointView& view)
         // compute covariance of the neighborhood
         Matrix3d B = math::computeCovariance(view, ids);
 
+        // Check if the covariance matrix is all zeros
+        if (B.isZero())
+        {
+            log()->get(LogLevel::Info)
+                << "Skipping point " << p.pointId()
+                << ". Covariance matrix is all zeros. This suggests a large "
+                   "number of redundant points. Consider using filters.sample "
+                   "with a small radius to remove redundant points.\n";
+            continue;
+        }
+
         // perform the eigen decomposition
         Eigen::SelfAdjointEigenSolver<Matrix3d> solver(B);
         if (solver.info() != Eigen::Success)

filters/NormalFilter.cpp

@@ -133,6 +133,18 @@ void NormalFilter::compute(PointView& view, KD3Index& kdi)
         // neighborhood composed of k-nearest neighbors.
         PointIdList neighbors = kdi.neighbors(p.pointId(), m_args->m_knn);
         auto B = math::computeCovariance(view, neighbors);
+
+        // Check if the covariance matrix is all zeros
+        if (B.isZero())
+        {
+            log()->get(LogLevel::Info)
+                << "Skipping point " << p.pointId()
+                << ". Covariance matrix is all zeros. This suggests a large "
+                   "number of redundant points. Consider using filters.sample "
+                   "with a small radius to remove redundant points.\n";
+            continue;
+        }
+
         SelfAdjointEigenSolver<Matrix3d> solver(B);
         if (solver.info() != Success)
             throwError("Cannot perform eigen decomposition.");

filters/PlaneFitFilter.cpp

@@ -134,6 +134,17 @@ void PlaneFitFilter::setPlaneFit(PointView& view, const PointId& i)
     // Compute covariance of the neighbors.
     Matrix3d B = math::computeCovariance(view, neighbors);
 
+    // Check if the covariance matrix is all zeros
+    if (B.isZero())
+    {
+        log()->get(LogLevel::Info)
+            << "Skipping point " << i
+            << ". Covariance matrix is all zeros. This suggests a large "
+               "number of redundant points. Consider using filters.sample "
+               "with a small radius to remove redundant points.\n";
+        return;
+    }
+
     // Perform the eigen decomposition, using the eigenvector of the smallest
     // eigenvalue as the normal.
     Eigen::SelfAdjointEigenSolver<Matrix3d> solver(B);

scripts/ci/current-environment.yml

@@ -2,7 +2,8 @@ name: pdal-build
 channels:
   - conda-forge
 dependencies:
-  - libgdal
+  - libgdal-core
+  - geotiff
   - libpq
   - curl
   - hdf5

scripts/ci/fixed-environment.yml

@@ -2,8 +2,7 @@ name: pdal-build
 channels:
   - conda-forge
 dependencies:
-  - libgdal=3.8
-  - boa
+  - libgdal=3.9
   - libpq
   - curl
   - hdf5
@@ -13,7 +12,7 @@ dependencies:
   - openssl
   - nitro
   - ceres-solver
-  - python=3.10
+  - python=3.12
   - conda
   - mamba
   - conda-build