Merge branch 'master' into RimlsPlaneFit

poncateam · Jun 26, 2023 · 41e9221 · 41e9221
2 parents 84f543e + 493df42
commit 41e9221
Show file tree

Hide file tree

Showing 23 changed files with 287 additions and 104 deletions.
diff --git a/CHANGELOG b/CHANGELOG
@@ -10,9 +10,12 @@ improved doc.
 - API
     - [spatialPartitioning] Add kd-tree traits type (#80)
     - [fitting] Add Primitive::getNumNeighbors (#86)
+    - [fitting] Change naming convention and rational for principal curvatures in CurvatureEstimatorBase (#94)
+    - [fitting] Add class MeanNormalDer to compute derivatives of the mean normal vector (#78)
 
 - Bug-fixes and code improvements
     - [spatialPartitioning] Fix unwanted function hiding with DryFit::setWeightFunc (#86)
+    - [spatialPartitioning] Fix missing include directive in kdTreeTraits.h (#92)
     - [fitting] Fix a potential bug when using multi-pass fitting (#89)
     - [fitting] Fix a bug in CovarianceFit (#93)
     - [fitting] Remove deadcode in Basket (#86)
@@ -22,7 +25,7 @@ improved doc.
     - Fix WeightKernel test (failing on windows due to finite differences) (#91)
 
 -Docs
-    - [spatialPartitioning] Update module page with a minimal doc and examples (#86)
+    - [spatialPartitioning] Update module page with a minimal doc and examples (#86, #92)
     - [spatialPartitioning] Add NanoFlann example (#86)
 
 --------------------------------------------------------------------------------

diff --git a/Ponca/src/Fitting/curvature.h b/Ponca/src/Fitting/curvature.h
@@ -16,7 +16,8 @@ namespace Ponca
     template < class DataPoint, class _WFunctor, int DiffType, typename T>
 /**
  *
- * \brief Base class for any 3d curvature estimator: holds \f$k_1\f$, \f$k_2\f$ and associated vectors
+ * \brief Base class for any 3d curvature estimator: holds \f$k_{\min}\f$, \f$k_{\max}\f$ and associated vectors,
+ * such that \f$ k_{\min} <= k_{\max} \f$
  */
     class CurvatureEstimatorBase : public T
     {
@@ -30,14 +31,14 @@ namespace Ponca
         };
 
     private:
-        /// \brief Principal curvature with highest absolute magnitude
-        Scalar m_k1 {0},
-        /// \brief Principal curvature with smallest absolute magnitude
-        m_k2 {0};
-        /// \brief Direction associated to the principal curvature with highest absolute magnitude
-        VectorType m_v1 {VectorType::Zero()},
-        /// \brief Direction associated to the principal curvature with highest smallest magnitude
-        m_v2 {VectorType::Zero()};
+        /// \brief Minimal principal curvature
+        Scalar m_kmin {0},
+        /// \brief Maximal principal curvature
+        m_kmax {0};
+        /// \brief Direction associated to the minimal principal curvature
+        VectorType m_vmin {VectorType::Zero()},
+        /// \brief Direction associated to the maximal principal curvature
+        m_vmax {VectorType::Zero()};
 
         /// \brief Internal state indicating if the curvature are set to default (false), or have been computed (true)
         bool m_isValid {false};
@@ -54,38 +55,32 @@ namespace Ponca
         /**************************************************************************/
         /* Use results                                                            */
         /**************************************************************************/
-        //! \brief Returns an estimate of the first principal curvature value
-        //!
-        //! It is the greatest curvature in <b>absolute value</b>.
-        PONCA_MULTIARCH inline Scalar k1() const { return m_k1; }
-
-        //! \brief Returns an estimate of the second principal curvature value
-        //!
-        //! It is the smallest curvature in <b>absolute value</b>.
-        PONCA_MULTIARCH inline Scalar k2() const { return m_k2; }
-
-        //! \brief Returns an estimate of the first principal curvature direction
-        //!
-        //! It is the greatest curvature in <b>absolute value</b>.
-        PONCA_MULTIARCH inline VectorType k1Direction() const { return m_v1; }
-
-        //! \brief Returns an estimate of the second principal curvature direction
-        //!
-        //! It is the smallest curvature in <b>absolute value</b>.
-        PONCA_MULTIARCH inline VectorType k2Direction() const { return m_v2; }
+        //! \brief Returns an estimate of the minimal principal curvature value
+        PONCA_MULTIARCH inline Scalar kmin() const { return m_kmin; }
+
+        //! \brief Returns an estimate of the maximal principal curvature value
+        PONCA_MULTIARCH inline Scalar kmax() const { return m_kmax; }
+
+        //! \brief Returns an estimate of the minimal principal curvature direction
+        PONCA_MULTIARCH inline VectorType kminDirection() const { return m_vmin; }
+
+        //! \brief Returns an estimate of the maximal principal curvature direction
+        PONCA_MULTIARCH inline VectorType kmaxDirection() const { return m_vmax; }
 
         //! \brief Returns an estimate of the mean curvature
-        PONCA_MULTIARCH inline Scalar kMean() const { return (m_k1 + m_k2)/Scalar(2);}
+        PONCA_MULTIARCH inline Scalar kMean() const { return (m_kmin + m_kmax)/Scalar(2);}
 
         //! \brief Returns an estimate of the Gaussian curvature
-        PONCA_MULTIARCH inline Scalar GaussianCurvature() const { return m_k1 * m_k2;}
+        PONCA_MULTIARCH inline Scalar GaussianCurvature() const { return m_kmin * m_kmax;}
 
     protected:
         /// \brief Set curvature values. To be called in finalize() by child classes
         ///
-        /// Principal curvatures are re-ordered depending on their magnitude,
-        /// such that \f$ |k_1| > |k_2| \f$
-        PONCA_MULTIARCH inline void setCurvatureValues(Scalar k1, Scalar k2, const VectorType& v1, const VectorType& v2);
+        /// If the given parameters are such that \f$ k_{\min} > k_{\max} \f$, then this
+        /// method swap the two curvature values and directions and store them such that
+        /// \f$ k_{\min} <= k_{\max} \f$.
+        ///
+        PONCA_MULTIARCH inline void setCurvatureValues(Scalar kmin, Scalar kmax, const VectorType& vmin, const VectorType& vmax);
     };
 
 } //namespace Ponca

diff --git a/Ponca/src/Fitting/curvature.hpp b/Ponca/src/Fitting/curvature.hpp
@@ -1,36 +1,32 @@
-#include PONCA_MULTIARCH_INCLUDE_STD(cmath)   //abs
 
 namespace Ponca {
     template<class DataPoint, class _WFunctor, int DiffType, typename T>
     void
     CurvatureEstimatorBase<DataPoint, _WFunctor, DiffType, T>::init(const VectorType &_evalPos) {
         Base::init(_evalPos);
-        m_k1 = 0;
-        m_k2 = 0;
-        m_v1 = VectorType::Zero();
-        m_v2 = VectorType::Zero();
+        m_kmin = 0;
+        m_kmax = 0;
+        m_vmin = VectorType::Zero();
+        m_vmax = VectorType::Zero();
         m_isValid = false;
     }
 
 
     template<class DataPoint, class _WFunctor, int DiffType, typename T>
     void
     CurvatureEstimatorBase<DataPoint, _WFunctor, DiffType, T>::setCurvatureValues(
-            Scalar k1, Scalar k2, const VectorType &v1, const VectorType &v2) {
-        PONCA_MULTIARCH_STD_MATH(abs)
-
-        if (abs(k1) < abs(k2)) {
-            m_k1 = k2;
-            m_k2 = k1;
-            m_v1 = v2;
-            m_v2 = v1;
+            Scalar kmin, Scalar kmax, const VectorType &vmin, const VectorType &vmax) {
+        if(kmin <= kmax) {
+            m_kmin = kmin;
+            m_kmax = kmax;
+            m_vmin = vmin;
+            m_vmax = vmax;
         } else {
-            m_k1 = k1;
-            m_k2 = k2;
-            m_v1 = v1;
-            m_v2 = v2;
+            m_kmin = kmax;
+            m_kmax = kmin;
+            m_vmin = vmax;
+            m_vmax = vmin;
         }
-
         m_isValid = true;
     }
 }
diff --git a/Ponca/src/Fitting/curvatureEstimation.hpp b/Ponca/src/Fitting/curvatureEstimation.hpp
@@ -148,35 +148,35 @@ NormalCovarianceCurvatureEstimator<DataPoint, _WFunctor, DiffType, T>::finalize
 
         m_solver.computeDirect(m_cov);
 
-        Scalar k1 = m_solver.eigenvalues()(1);
-        Scalar k2 = m_solver.eigenvalues()(2);
+        Scalar kmin = m_solver.eigenvalues()(1);
+        Scalar kmax = m_solver.eigenvalues()(2);
 
-        VectorType v1 = m_solver.eigenvectors().col(1);
-        VectorType v2 = m_solver.eigenvectors().col(2);
+        VectorType vmin = m_solver.eigenvectors().col(1);
+        VectorType vmax = m_solver.eigenvectors().col(2);
 
         // fallback
         // TODO(thib) which epsilon value should be chosen ?
         Scalar epsilon = Scalar(1e-3);
-        if(k1<epsilon && k2<epsilon)
+        if(kmin<epsilon && kmax<epsilon)
         {
-            k1 = Scalar(0);
-            k2 = Scalar(0);
+            kmin = Scalar(0);
+            kmax = Scalar(0);
 
             // set principal directions from normals center of gravity
             VectorType n = m_cog.normalized();
             Index i0 = -1, i1 = -1, i2 = -1;
             n.rowwise().squaredNorm().minCoeff(&i0);
             i1 = (i0+1)%3;
             i2 = (i0+2)%3;
-            v1[i0] = 0;
-            v1[i1] = n[i2];
-            v1[i2] = -n[i1];
-            v2[i0] = n[i1]*n[i1] + n[i2]*n[i2];
-            v2[i1] = -n[i1]*n[i0];
-            v2[i2] = -n[i2]*n[i0];
+            vmin[i0] = 0;
+            vmin[i1] = n[i2];
+            vmin[i2] = -n[i1];
+            vmax[i0] = n[i1]*n[i1] + n[i2]*n[i2];
+            vmax[i1] = -n[i1]*n[i0];
+            vmax[i2] = -n[i2]*n[i0];
         }
 
-        Base::setCurvatureValues(k1, k2, v1, v2);
+        Base::setCurvatureValues(kmin, kmax, vmin, vmax);
     }
     return res;
 }
@@ -259,8 +259,8 @@ ProjectedNormalCovarianceCurvatureEstimator<DataPoint, _WFunctor, DiffType, T>::
 
         m_solver.computeDirect(m_cov);
 
-        Base::m_k1 = m_solver.eigenvalues()(0);
-        Base::m_k2 = m_solver.eigenvalues()(1);
+        Base::m_kmin = m_solver.eigenvalues()(0);
+        Base::m_kmax = m_solver.eigenvalues()(1);
 
         // transform from local plane coordinates to world coordinates
         Base::m_v1 = m_tframe * m_solver.eigenvectors().col(0);
@@ -269,10 +269,10 @@ ProjectedNormalCovarianceCurvatureEstimator<DataPoint, _WFunctor, DiffType, T>::
         //TODO(thib) which epsilon value should be chosen ?
 //        Scalar epsilon = Eigen::NumTraits<Scalar>::dummy_precision();
         Scalar epsilon = Scalar(1e-3);
-        if(Base::m_k1<epsilon && Base::m_k2<epsilon)
+        if(Base::m_kmin<epsilon && Base::m_kmax<epsilon)
         {
-            Base::m_k1 = Scalar(0);
-            Base::m_k2 = Scalar(0);
+            Base::m_kmin = Scalar(0);
+            Base::m_kmax = Scalar(0);
 
             // set principal directions from fitted plane
             Base::m_v2 = m_tframe.col(0);

diff --git a/Ponca/src/Fitting/mean.h b/Ponca/src/Fitting/mean.h
@@ -87,7 +87,7 @@ namespace Ponca {
             PROVIDES_MEAN_POSITION_DERIVATIVE,    /*!< \brief Provides derivative of the mean position*/
         };
 
-        /*! \brief Derivatives of the of the input points vectors */
+        /*! \brief Derivatives of the input points vectors */
         VectorArray m_dSumP {VectorArray::Zero()};
 
     public:
@@ -122,6 +122,56 @@ namespace Ponca {
 
     }; //class MeanPositionDer
 
+    template<class DataPoint, class _WFunctor, int DiffType, typename T>
+    class MeanNormalDer : public T {
+        PONCA_FITTING_DECLARE_DEFAULT_TYPES
+        PONCA_FITTING_DECLARE_DEFAULT_DER_TYPES
+
+    protected:
+        enum {
+            Check = Base::PROVIDES_PRIMITIVE_DERIVATIVE && 
+                    Base::PROVIDES_MEAN_NORMAL,
+            PROVIDES_MEAN_NORMAL_DERIVATIVE,  /*!< \brief Provides derivative of the mean normal*/
+        };
+
+        private:
+        /*! \brief Derivatives of the input normals of the input points vectors*/
+        VectorArray m_dSumN {VectorArray::Zero()};
+
+    public:
+
+        PONCA_EXPLICIT_CAST_OPERATORS_DER(MeanNormalDer,meanNormalDer)
+        PONCA_FITTING_DECLARE_INIT
+        PONCA_FITTING_DECLARE_ADDNEIGHBOR_DER
+
+    /// \brief Compute the derivative of the mean normal vector of the input points. 
+    /// 
+    /// ### Step-by-step derivation of the mean normal : 
+    ///  Given the definition of the mean normal \f$ n(\mathbf{x}) = \frac{\sum_i w_\mathbf{x}(\mathbf{n_i}) \mathbf{n_i}}{\sum_i w_\mathbf{x}(\mathbf{n_i})} \f$, 
+    /// where \f$\left[\mathbf{n_i} \in \text{neighborhood}(\mathbf{x})\right]\f$ are all the normal samples in \f$\mathbf{x}\f$'s neighborhood.
+    ///
+    /// We denote \f$ t(\mathbf{x}) = \sum_i w_\mathbf{x}(\mathbf{n_i}) \mathbf{n_i} \f$ and \f$ s(\mathbf{x}) = \sum_i w_\mathbf{x}(\mathbf{n_i})\f$, 
+    /// such that \f$ n(\mathbf{x}) = \frac{t(\mathbf{x})}{s(\mathbf{x})}\f$.
+    ///
+    /// By definition, \f$ n'(\mathbf{x}) = \frac{s(\mathbf{x})t'(\mathbf{x}) - t(\mathbf{x})s'(\mathbf{x})}{s(\mathbf{x})^2}\f$.
+    ///
+    /// Assuming the weight of each normal is dependent on the position, we have \f$ s'(\mathbf{x}) = \sum_i w'_\mathbf{x}(\mathbf{n_i}) \f$.
+    ///
+    /// We rewrite \f$ t(\mathbf{x}) = \sum u(\mathbf{x})v(\mathbf{x}) \f$, with \f$ u(\mathbf{x}) = w_\mathbf{x}(\mathbf{n_i}) \f$ and \f$ v(\mathbf{x}) = \mathbf{n_i} \f$.
+    ///
+    /// Assuming the normal vectors themselves do not change with position, \f$v(\mathbf{x})\f$ is constant, its derivative is null, 
+    /// and so \f$ t'(\mathbf{x}) = \sum_i u'(\mathbf{x}) v(\mathbf{x}) = \sum_i w'_\mathbf{x}(\mathbf{n_i}) \mathbf{n_i} \f$.
+    ///
+    /// Which leads to \f$ n'(\mathbf{x}) = \frac{\sum_i w'\mathbf{x}(\mathbf{n_i}) \mathbf{n_i} - n(\mathbf{x})\sum w'(\mathbf{x})}{\sum_i w\mathbf{x}(\mathbf{n_i})} \f$.
+    /// \note This code is not directly tested. 
+
+    PONCA_MULTIARCH inline const VectorArray& dMeanNormal() const
+    { 
+        return m_dSumN; 
+    }
+
+    }; //class MeanNormalDer
+
 #include "mean.hpp"
 
 } //namespace Ponca
diff --git a/Ponca/src/Fitting/mean.hpp b/Ponca/src/Fitting/mean.hpp
@@ -68,3 +68,25 @@ MeanPositionDer<DataPoint, _WFunctor, DiffType, T>::addLocalNeighbor(Scalar w,
 
     return false;
 }
+
+
+template<class DataPoint, class _WFunctor, int DiffType, typename T>
+void
+MeanNormalDer<DataPoint, _WFunctor, DiffType, T>::init(const VectorType &_evalPos) {
+    Base::init(_evalPos);
+    m_dSumN.setZero();
+}
+
+template<class DataPoint, class _WFunctor, int DiffType, typename T>
+bool
+MeanNormalDer<DataPoint, _WFunctor, DiffType, T>::addLocalNeighbor(Scalar w,
+                                                                 const VectorType &localQ,
+                                                                 const DataPoint &attributes,
+                                                                 ScalarArray &dw) {
+    if (Base::addLocalNeighbor(w, localQ, attributes, dw)) {
+        m_dSumN += attributes.normal() * dw;
+        return true;
+    }
+
+    return false;
+}
diff --git a/Ponca/src/Fitting/orientedSphereFit.h b/Ponca/src/Fitting/orientedSphereFit.h
@@ -135,7 +135,7 @@ class OrientedSphereDerImpl : public T
 
 }; //class OrientedSphereDerImpl
 
-/// \brief Helper alias for Plane fitting on 3D points using CovariancePlaneFitImpl
+/// \brief Helper alias for Oriented Sphere fitting on 3D points using OrientedSphereDerImpl
     template < class DataPoint, class _WFunctor, int DiffType, typename T>
     using OrientedSphereDer =
         OrientedSphereDerImpl<DataPoint, _WFunctor, DiffType,

diff --git a/Ponca/src/SpatialPartitioning/KdTree/kdTree.h b/Ponca/src/SpatialPartitioning/KdTree/kdTree.h
@@ -300,7 +300,7 @@ public :
     NodeContainer m_nodes;
     IndexContainer m_indices;
 
-    LeafSizeType m_min_cell_size;
+    LeafSizeType m_min_cell_size; ///< Minimal number of points per leaf
     NodeCountType m_leaf_count; ///< Number of leaves in the Kdtree (computed during construction)
 };
 

diff --git a/Ponca/src/SpatialPartitioning/KdTree/kdTreeTraits.h b/Ponca/src/SpatialPartitioning/KdTree/kdTreeTraits.h
@@ -4,8 +4,11 @@
  file, You can obtain one at http://mozilla.org/MPL/2.0/.
 */
 
+
 #pragma once
 
+#include <Eigen/Geometry>
+
 namespace Ponca {
 #ifndef PARSED_WITH_DOXYGEN
 namespace internal

diff --git a/doc/src/example.mdoc b/doc/src/example.mdoc
@@ -18,6 +18,7 @@
     - \subpage example_cu_ssc_page : Calculate Screen Space Curvature using CUDA/C++.
     - \subpage example_python_ssc_page : Calculate Screen Space Curvature using CUDA/Python.
   - @ref spatialpartitioning "Spatial Partitioning Module"
+    - \subpage example_cxx_neighbor_search : Nearest, k-nearest and range neighbor searches using Ponca::KdTree.
     - \subpage example_cxx_nanoflann_page : Comparison between Nanoflann and Ponca KdTree APIs.
 
  */
diff --git a/doc/src/example_cxx_neighbor_search.mdoc b/doc/src/example_cxx_neighbor_search.mdoc
@@ -0,0 +1,14 @@
+/*
+ This Source Code Form is subject to the terms of the Mozilla Public
+ License, v. 2.0. If a copy of the MPL was not distributed with this
+ file, You can obtain one at http://mozilla.org/MPL/2.0/.
+*/
+
+
+/*!
+  \page example_cxx_neighbor_search Ponca::KdTree neighbor searches
+
+  This is an example of how to use Ponca::KdTree to perform nearest, k-nearest and range neighbor search.
+
+  \include cpp/ponca_neighbor_search.cpp
+*/