v0.3.3

include/pca.hpp

@@ -79,26 +79,21 @@ static inline PCAResult<real_t> pca_from_pointcloud(const PointCloud<real_t>& cl
     Eigen::SelfAdjointEigenSolver<Eigen::Matrix<real_t, 3, 3>> es(cov);
 
     // Sort the values and vectors in order of increasing eigenvalue
-    const auto ev = es.eigenvalues().real();
+    const auto ev = es.eigenvalues();
 
-    std::array<Eigen::Index, 3> indices = {0, 1, 2};
+    // Clamp eigen values
+    Vec3<real_t> val = {(std::max(ev(0), real_t(0.))), (std::max(ev(1), real_t(0.))), (std::max(ev(2), real_t(0.)))};
 
-    std::sort(
-        std::begin(indices), std::end(indices), [&](Eigen::Index i1, Eigen::Index i2) { return ev(i1) > ev(i2); });
-
-    Vec3<real_t> val = {
-        (std::max(ev(indices[0]), real_t(0.))), (std::max(ev(indices[1]), real_t(0.))),
-        (std::max(ev(indices[2]), real_t(0.)))};
-    Vec3<real_t> v0 = es.eigenvectors().col(indices[0]).real();
-    Vec3<real_t> v1 = es.eigenvectors().col(indices[1]).real();
-    Vec3<real_t> v2 = es.eigenvectors().col(indices[2]).real();
+    Vec3<real_t> v0 = es.eigenvectors().col(0);
+    Vec3<real_t> v1 = es.eigenvectors().col(1);
+    Vec3<real_t> v2 = es.eigenvectors().col(2);
 
     // To standardize the orientation of eigenvectors, we choose to enforce all eigenvectors
     // to be expressed in the Z+ half-space.
-    // Only the third eigenvector (v2) needs to be reoriented because it is the
+    // Only the third eigenvector (v2 the surface normal) needs to be reoriented because it is the
     // only one used in further computations.
     // TODO: In case we want to orient normal, this should be improved
-    if (v2(2) < real_t(0.)) { v2 = real_t(-1.) * v2; }
+    if (v2(2) < real_t(0.)) { v2 *= real_t(-1.); };
     return {val, v0, v1, v2};
 };
 
@@ -280,14 +275,12 @@ void compute_selected_features(
                             pca.val(2) * std::abs(pca.v2(2))};
 
                     feature_results[output_id] = unary_vector(2) / unary_vector.norm();
-                    // TODO: Jakteristics compute this as feature_results[output_id] = real_t(1.0) -
-                    // std::abs(pca.v2(2));
-                    // It seems to be the most common formula for the verticality in the literature
                 }
                 break;
             case EFeatureID::Verticality:
                 // The verticality as defined in most of the papers
                 // http://lareg.ensg.eu/labos/matis/pdf/articles_revues/2015/isprs_wjhm_15.pdf
+                // It's the way it's computed in jakteristics
                 feature_results[output_id] = real_t(1.0) - std::abs(pca.v2(2));
                 break;
             case EFeatureID::Eigentropy:

pyproject.toml

@@ -8,7 +8,7 @@ build-backend = "scikit_build_core.build"
 
 [project]
 name = "pgeof"
-version = "0.3.2"
+version = "0.3.3"
 readme = "README.md"
 description = "Compute the geometric features associated with each point's neighborhood:"
 requires-python = ">=3.8,<3.14"