add doc about UnorientedSphereFitImpl

Ponca/src/Fitting/unorientedSphereFit.h

@@ -18,7 +18,24 @@ namespace Ponca
 /*!
     \brief Algebraic Sphere fitting procedure on point sets with non-oriented normals
 
-    Method published in \cite Chen:2013:NOMG.
+    This method published in \cite Chen:2013:NOMG maximizes the sum of squared dot product between the input normal vectors and the gradient of the algebraic sphere.
+    The maximization is done under the constraint that the norm of the gradient is unitary on average.
+    In practice, it amounts to solve the generalized eigenvalue problem
+    \f[
+        A \mathbf{u} = \lambda Q \mathbf{u}
+    \f]
+    where
+    \f[
+        \mathbf{u} = \begin{bmatrix} u_l \\ u_q \end{bmatrix}
+    \f]
+    \f[
+        A = \sum_i w_i \begin{bmatrix} \mathbf{n}_i \\ \mathbf{n}_i^T\mathbf{p}_i \end{bmatrix}\begin{bmatrix} \mathbf{n}_i^T & \mathbf{n}_i^T\mathbf{p}_i \end{bmatrix}
+    \f]
+    \f[
+        Q = \frac{1}{\sum_i w_i} \begin{bmatrix} I & \sum_i w_i \mathbf{p}_i \\ \sum_i w_i \mathbf{p}_i^T & \sum_i w_i \mathbf{p}_i^T\mathbf{p}_i \end{bmatrix}
+    \f]
+    The constant coefficient \f$u_c\f$ is computed as in \ref OrientedSphereFitImpl by minimizing the sum of squared potential evaluated at the points \f$\mathbf{p}_i\f$.
+    This fitting method corresponds to the case where the scalar field is \f$ f^* \f$ in the Section 3.3 of \cite Chen:2013:NOMG.
 
     \inherit Concept::FittingProcedureConcept
 

Ponca/src/Fitting/unorientedSphereFit.hpp

@@ -137,6 +137,25 @@ UnorientedSphereDerImpl<DataPoint, _WFunctor, DiffType, T>::finalize()
 
         const Scalar invSumW = Scalar(1) / Base::getWeightSum();
 
+        // the derivative of (A vi = li Q vi), where A and A are symmetric real, is
+        //
+        //     dA vi + A dvi = dli Q vi + li dQ vi + li Q dvi
+        //
+        // left-multiply by vj (i!=j, associated to eigenvalue lj)
+        //
+        //     vj . dA vi + vj . A dvi = dli vj . Q vi + li vj . dQ vi + li vj . Q dvi
+        //                  ^^^^^^           ^^^^^^^^^
+        //                  = lj vj . B         = 0
+        //
+        // which simplified to
+        //
+        //     (Q vj) . dvi = 1/(li - lj) vj . dQ vi
+        //
+        // therefore
+        //
+        //     dvi = sum_j (1/(li - lj) vj . dQ vi) Q vj
+        //         = Q sum_j (1/(li - lj) vj . ((dA - li dQ) vi)) vj
+        //
         for(int dim = 0; dim < Base::NbDerivatives; ++dim)
         {
             MatrixBB dQ;