Fix sphere fitting eigen solver

CHANGELOG

@@ -14,6 +14,7 @@ Current head (v.1.2 RC)
     - [fitting] Use variadic template for basket extensions (#85)
     - [fitting] Fix current status issue (#108)
     - [spatialPartitioning] Fix potential compilation issues in KnnGraph (#111) 
+    - [fitting] Fix sphere fit eigen solver (#112)
 
 -Docs
     - [fitting] Clarify documentation on FIT_RESULT (#108)

Ponca/src/Fitting/sphereFit.h

@@ -13,9 +13,40 @@ namespace Ponca
 {
 
 /*!
-    \brief Algebraic Sphere fitting procedure on point sets without normals
-
-    Method published in \cite Guennebaud:2007:APSS.
+    \brief Algebraic Sphere fitting procedure on point set without normals
+
+    This method published in \cite Guennebaud:2007:APSS minimizes
+    \f[
+        \mathcal{L}(\mathbf{u}) = \frac{1}{2} \sum_i w_i f_{\mathbf{u}}(\mathbf{x}_i)^2 = \frac{1}{2} \mathbf{u}^T A \mathbf{u}
+    \f]
+    with \f$ A = \sum_i w_i \tilde{\mathbf{x}_i}  \tilde{\mathbf{x}_i}^T\f$ and \f$ f_{\mathbf{u}} \f$ the algebraic sphere defined by
+    \f[
+        f_{\mathbf{u}}(\mathbf{x}) =
+        u_c + \mathbf{u}_l.\mathbf{x} + u_q \mathbf{x}.\mathbf{x} =
+        \begin{bmatrix}
+        1 & \mathbf{x}^T & \mathbf{x}.\mathbf{x}
+        \end{bmatrix}
+        \begin{bmatrix}
+        u_c \\ u_l \\ u_q
+        \end{bmatrix}
+        = \tilde{\mathbf{x}}^T \mathbf{u},
+    \f]
+    under the constraint (unitary gradient onto the surface)
+    \f[
+        \|\mathbf{u}_l\|^2 - 4 u_c u_q = \mathbf{u}^T C \mathbf{u} = 1
+    \f]
+    where
+    \f[
+        C =
+        \begin{bmatrix}
+            0 &   &         &   & -2 \\
+              & 1 &         &   &    \\
+              &   &  \ddots &   &    \\
+              &   &         & 1 &    \\
+           -2 &   &         &   &  0
+        \end{bmatrix}
+    \f]
+    which amounts to solve the generalized eigenvalue problem \f$ A\mathbf{u} = \lambda C \mathbf{u} \f$.
 
     \inherit Concept::FittingProcedureConcept
 
@@ -35,12 +66,21 @@ PONCA_FITTING_DECLARE_DEFAULT_TYPES
     typedef Eigen::Matrix<Scalar, DataPoint::Dim+2, 1>      VectorA;
     typedef Eigen::Matrix<Scalar, DataPoint::Dim+2, DataPoint::Dim+2>  MatrixA;
 
+public:
+    using Solver = Eigen::EigenSolver<MatrixA>;
+
+protected:
     // computation data
     MatrixA  m_matA {MatrixA::Zero()};  /*!< \brief Covariance matrix of [1, p, p^2] */
 
+    Solver m_solver;
+
 public:
     PONCA_EXPLICIT_CAST_OPERATORS(SphereFitImpl,sphereFit)
     PONCA_FITTING_DECLARE_INIT_ADD_FINALIZE
+
+    PONCA_MULTIARCH inline const Solver& solver() const { return m_solver; }
+
 }; //class SphereFit
 
 /// \brief Helper alias for Sphere fitting on 3D points using SphereFitImpl

Ponca/src/Fitting/sphereFit.hpp

@@ -64,17 +64,18 @@ SphereFitImpl<DataPoint, _WFunctor, T>::finalize ()
     invCpratt.template topLeftCorner<1,1>()     << 0;
     invCpratt.template bottomRightCorner<1,1>() << 0;
 
-    MatrixA M = invCpratt * m_matA;
-    // go to positive semi-definite matrix to be compatible with
-    // SelfAdjointEigenSolver requirements
-    // Note: This does not affect the eigen vectors order
-    Eigen::SelfAdjointEigenSolver<MatrixA> solver;
+    // Remarks:
+    //   A and C are symmetric so all eigenvalues and eigenvectors are real
+    //   we look for the minimal positive eigenvalue (eigenvalues may be negative)
+    //   C^{-1}A is not symmetric
+    //   calling Eigen::GeneralizedEigenSolver on (A,C) and Eigen::EigenSolver on C^{-1}A is equivalent
+    //   C is not positive definite so Eigen::GeneralizedSelfAdjointEigenSolver cannot be used
 #ifdef __CUDACC__
-    solver.computeDirect(M.transpose() * M);
+    m_solver.computeDirect(invCpratt * m_matA);
 #else
-    solver.compute(M.transpose() * M);
+    m_solver.compute(invCpratt * m_matA);
 #endif
-    VectorA eivals = solver.eigenvalues().real();
+    VectorA eivals = m_solver.eigenvalues().real();
     int minId = -1;
     for(int i=0 ; i<DataPoint::Dim+2 ; ++i)
     {
@@ -84,7 +85,7 @@ SphereFitImpl<DataPoint, _WFunctor, T>::finalize ()
     }
 
     //mLambda = eivals(minId);
-    VectorA vecU = solver.eigenvectors().col(minId).real();
+    VectorA vecU = m_solver.eigenvectors().col(minId).real();
     Base::m_uq = vecU[1+DataPoint::Dim];
     Base::m_ul = vecU.template segment<DataPoint::Dim>(1);
     Base::m_uc = vecU[0];

examples/cpp/ponca_basic_cpu.cpp

@@ -19,6 +19,7 @@ file, You can obtain one at http://mozilla.org/MPL/2.0/.
 #include <Ponca/src/Fitting/gls.h>
 #include <Ponca/src/Fitting/orientedSphereFit.h>
 #include <Ponca/src/Fitting/unorientedSphereFit.h>
+#include <Ponca/src/Fitting/sphereFit.h>
 #include <Ponca/src/Fitting/weightFunc.h>
 #include <Ponca/src/Fitting/weightKernel.h>
 #include <Ponca/src/Fitting/curvatureEstimation.h>
@@ -72,7 +73,7 @@ typedef DistWeightFunc<MyPoint,SmoothWeightKernel<Scalar> > WeightFunc;
 using Fit1 = Basket<MyPoint,WeightFunc,OrientedSphereFit,   GLSParam>;
 using Fit2 = Basket<MyPoint,WeightFunc,UnorientedSphereFit, GLSParam>;
 using Fit3 = BasketDiff< Fit1, FitSpaceDer, OrientedSphereDer, GLSDer, CurvatureEstimatorBase, NormalDerivativesCurvatureEstimator>;
-
+using Fit4 = Basket<MyPoint,WeightFunc,SphereFit, GLSParam>;
 
 template<typename Fit>
 void test_fit(Fit& _fit, const KdTree<MyPoint>& tree, const VectorType& _p)
@@ -158,4 +159,10 @@ int main()
     cout << fit3.kminDirection() << endl << endl;
     cout << fit3.kmaxDirection() << endl;
   }
+
+  std::cout << "\n\n====================\nSphereFit:\n";
+  Fit4 fit4;
+  test_fit(fit4, tree, p);
+
+  return 0;
 }