Helper functor in PolygonalCalculus to project vertices onto their tangent plane.

ChangeLog.md

@@ -43,6 +43,11 @@
   - The WindingNumberShape class can output the raw winding number values
     (David Coeurjolly,[#1719](https://github.com/DGtal-team/DGtal/issues/1719))
 
+- *DEC*
+  - New helper functor to construct an embedder to correct the PolygonalCalculs
+    (projection onto estimated tangent planes) (David Coeurjolly,
+    [#1730](https://github.com/DGtal-team/DGtal/issues/17309))
+
 - *Geometry package*
   - Add creation of polytopes from segments and triangles in
     ConvexityHelper and 3-5xfaster full subconvexity tests for triangles

examples/polyscope-examples/dgtalCalculus-geodesic.cpp

@@ -93,25 +93,11 @@ void precompute()
     auto surfels   = SH3::getSurfelRange( surface, params2 );
     iinormals = SHG3::getIINormalVectors(binary_image, surfels,params2);
     trace.info()<<iinormals.size()<<std::endl;
-    auto myProjEmbedder = [&](Face f, Vertex v)
-    {
-      const auto nn = iinormals[f];
-      RealPoint centroid(0.0,0.0,0.0); //centroid of the original face
-      auto cpt=0;
-      for(auto v: surfmesh.incidentVertices(f))
-      {
-        cpt++;
-        centroid += surfmesh.position(v);
-      }
-      centroid = centroid / (double)cpt;
-      RealPoint p = surfmesh.position(v);
-      auto cp = p-centroid;
-      RealPoint q = p - nn.dot(cp)*nn;
-      return q;
-    };
     psMesh->addFaceVectorQuantity("II normals", iinormals);
+
     calculus = new PolyCalculus(surfmesh);
-    calculus->setEmbedder( myProjEmbedder );
+    functors::EmbedderFromNormalVectors<Z3i::RealPoint, Z3i::RealVector> embedderFromNormals(iinormals,surfmesh);
+    calculus->setEmbedder( embedderFromNormals );
   }
 
   heat = new GeodesicsInHeat<PolyCalculus>(calculus);

src/DGtal/dec/PolygonalCalculus.h

@@ -44,6 +44,64 @@
 namespace DGtal
 {
 
+  namespace functors {
+    /**
+     *
+     * \brief Functor that projects a face vertex of a surface mesh onto the tangent plane
+     * given by a per-face normal vector.
+     * This functor can be used in PolygonalCalculus to correct the embedding of
+     * digital surfaces using an estimated normal vector field (see @cite coeurjolly2022simple).
+     *
+     * @note when used in PolygonalCalculus, all operators being invariant by translation, all
+     * tangent planes pass through the origin (0,0,0) (no offest).
+     *
+     * @tparam TRealPoint a model of points @f$\mathbb{R}^3@f$ (e.g. PointVector).
+     * @tparam TRealVector a model of vectors in @f$\mathbb{R}^3@f$ (e.g. PointVector).
+     */    template <typename TRealPoint, typename TRealVector>
+    struct EmbedderFromNormalVectors
+    {
+      ///Type of SurfaceMesh
+      typedef SurfaceMesh<TRealPoint, TRealVector> MySurfaceMesh;
+      ///Vertex type
+      typedef typename MySurfaceMesh::Vertex Vertex;
+      ///Face type
+      typedef typename MySurfaceMesh::Face Face;
+      ///Position type
+      typedef typename MySurfaceMesh::RealPoint Real3dPoint;
+
+      EmbedderFromNormalVectors() = delete;
+
+      /// Constructor from an array of normal vectors and a surface mesh instance.
+      /// @param normals a vector of per face normal vectors (same ordering as the SurfaceMesh face indicies).
+      /// @param surfmesh an instance of SurfaceMesh
+      EmbedderFromNormalVectors(ConstAlias<std::vector<Real3dPoint>> normals,
+                                ConstAlias<MySurfaceMesh> surfmesh)
+      {
+        myNormals     = &normals;
+        mySurfaceMesh = &surfmesh;
+      }
+
+      /// Project a face vertex onto its tangent plane (given by the per-face estimated
+      /// normal vector).
+      ///
+      /// @param f the face that contains the vertex
+      /// @param v the vertex to project
+      Real3dPoint operator()(const Face &f,const Vertex &v)
+      {
+        const auto nn = (*myNormals)[f];
+        Real3dPoint p = mySurfaceMesh->position(v);
+        return p - nn.dot(p)*nn;
+      }
+
+      ///Alias to the normal vectors
+      const std::vector<Real3dPoint> *myNormals;
+      ///Alias to the surface mesh
+      const MySurfaceMesh *mySurfaceMesh;
+    };
+  }
+
+
+
 /////////////////////////////////////////////////////////////////////////////
 // template class PolygonalCalculus
 /**
@@ -223,7 +281,7 @@ class PolygonalCalculus
   {
     myEmbedder = externalFunctor;
   }
-
+  /// @}
 
   // ----------------------- Per face operators --------------------------------------
   //MARK: Per face operator on scalars

src/DGtal/dec/doc/modulePolygonalCalculus.dox

@@ -120,7 +120,11 @@ PolygonalCalculus<SurfMesh>::setEmbedder for an example.
 can have different embeddings for all its incident faces.
 
 
+
 \subsection susub1 Basic operators
+
+
+
 
 We first describe some standard per face operators. Note that for all
 extrinsic operators that require the vertex position in @f$
@@ -253,40 +257,75 @@ used when one solves a weak problem and wishes to get a pointwise
 per-vertex solution.
 
 
-\section secLap Example: Solving a Laplace problem
-
-Let suppose we want to solve the following Laplace problem for data interpolation:
-\f{eqnarray*}{
-        \Delta_\Omega u& = 0  \\ 
-        & s.t. u = g \text{ on } \partial\Omega
-\f}
-
-We want to solve that problem on a polygonal mesh @f$\Omega@f$
-(digital surface here) with a boundary and some scalar values attached
-to boundary vertices, or sampled on the object surface.
-
-Furthermore, the discrete version of the Laplace problem boils down to
-a simple linear problem using on the discrete Laplace-Beltrami sparse
-matrix.
-
-We also use class DirichletConditions to enforce Dirichlet boundary
-conditions on the system.
-
-The overall code is:
-\snippet dgtalCalculus-poisson.cpp PolyDEC-init
-
-Leading to the following results (see \ref dgtalCalculus-poisson.cpp):
-
-Surface  | Boundary condition @f$ g@f$ | Solution @f$ u @f$
---|--|--
-@image html images/poly/poisson-surf.png "" | @image html images/poly/poisson-g.png "" | @image html images/poly/poisson-u.png "" 
-@image html images/poly/bunny-init.png "" | @image html images/poly/bunny-g.png "" | @image html images/poly/bunny-u.png "" 
-@image html images/poly/cat-init.png "" | @image html images/poly/cat-g.png "" | @image html images/poly/cat-u.png "" 
+
+\section sectCorrected Corrected Calculus using Estimated Normal Vectors
+
+On digital surfaces, solving PDE on original embedding with axis aligned quad surfaces may fail to correctly capture the surface metric.
+ As discussed in @cite coeurjolly2022simple, given an estimation of the tangent bundle of the discrete surface (for instance using
+ estimated normal vectors from @ref moduleIntegralInvariant, or @ref moduleVCM, cf @ref moduleShortcuts),
+ one can implicitly project each face to a prescribed tangent plane and perform the computations on this  new embedding of the geometry.
 
-\subsection Global Vector Calculus
+Geodesic distances without correction | Geodesic distances with correction
+ --|--
+@image html images/poly/corrected-without.png "" | @image html images/poly/corrected-with.png ""
 
-Global Vector Laplace/Poisson problems can also be solved by the same way, using instead PolygonalCalculus<SurfMesh>::globalConnectionLaplace() and PolygonalCalculus<SurfMesh>::doubledGlobalLumpedMassMatrix(). One can find examples of such use in the \ref VectorsInHeat class.
+
+The functor functors::EmbedderFromNormalVectors can be used to implicitly project Face vertices onto the prescribed tangent plane. A calssical usage is the following one:
+
+ @code
+ //A surface mesh. Eg. primal surface of a digital surface
+ SurfaceMesh< Z3i::RealPoint, Z3i::RealVector > surfmesh(...);
 
+ //Per face normal vector estimation
+ std::vector<Z3i::RealVector> ii_normals = ....
+
+ //New embedder using tangent plane projection of face vertices
+ functors::EmbedderFromNormalVectors<Z3i::RealPoint, Z3i::RealVector> embedderFromNormals(ii_normals,surfmesh);
+
+ //The calculus instance with the new embedder
+ PolygonalCalculus<Z3i::RealPoint,23i::RealVector> calculus(surfmesh);
+ calculus->setEmbedder( embedderFromNormals );
+ @endcode
+ A complete code is given in the @ref dgtalCalculus-geodesic.cpp example.
+
+ \section secLap Example: Solving a Laplace problem
+
+ Let suppose we want to solve the following Laplace problem for data interpolation:
+ \f{eqnarray*}{
+ \Delta_\Omega u& = 0  \\
+ & s.t. u = g \text{ on } \partial\Omega
+ \f}
+
+ We want to solve that problem on a polygonal mesh @f$\Omega@f$
+ (digital surface here) with a boundary and some scalar values attached
+ to boundary vertices, or sampled on the object surface.
+
+ Furthermore, the discrete version of the Laplace problem boils down to
+ a simple linear problem using on the discrete Laplace-Beltrami sparse
+ matrix.
+
+ We also use class DirichletConditions to enforce Dirichlet boundary
+ conditions on the system.
+
+ The overall code is:
+ \snippet dgtalCalculus-poisson.cpp PolyDEC-init
+
+ Leading to the following results (see \ref dgtalCalculus-poisson.cpp):
+
+ Surface  | Boundary condition @f$ g@f$ | Solution @f$ u @f$
+ --|--|--
+ @image html images/poly/poisson-surf.png "" | @image html images/poly/poisson-g.png "" | @image html images/poly/poisson-u.png ""
+ @image html images/poly/bunny-init.png "" | @image html images/poly/bunny-g.png "" | @image html images/poly/bunny-u.png ""
+ @image html images/poly/cat-init.png "" | @image html images/poly/cat-g.png "" | @image html images/poly/cat-u.png ""
+
+ \subsection Global Vector Calculus
+
+ Global Vector Laplace/Poisson problems can also be solved by the same way, using instead PolygonalCalculus<SurfMesh>::globalConnectionLaplace() and PolygonalCalculus<SurfMesh>::doubledGlobalLumpedMassMatrix(). One can find examples of such use in the \ref VectorsInHeat class.
+
+
+
+
+
 \section sectMisc Miscellaneous
 
 \subsection sectPolygonalCalculusHP Cache mechanisms and high-performance computing
@@ -325,7 +364,7 @@ Then, cached operators can be accessed and combined:
 auto Mf =   cachefaceArea[f] * cacheU[f].transpose()*cacheU[f] + lambda * cacheP[f].transpose() * cacheP[f];
 @endcode 
 
-
+ 
 
 */
 

src/DGtal/doc/global.bib

@@ -81,6 +81,14 @@ @inproceedings{nielson2003shrouds
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 %% DEC/AT
+@inproceedings{coeurjolly2022simple,
+  title        = {A Simple Discrete Calculus for Digital Surfaces},
+  author       = {Coeurjolly, David and Lachaud, Jacques-Olivier},
+  booktitle    = {International Conference on Discrete Geometry and Mathematical Morphology},
+  pages        = {341--353},
+  year         = {2022},
+  organization = {Springer}
+}
 
 @article{ambrosio1990approximation,
   title={Approximation of functional depending on jumps by elliptic functional via t-convergence},