DGtal 1.3.0

Part of the Geometry package.
This part of the manual describes the implementation of several curvature measures on quite arbitrary meshes (any instance of SurfaceMesh). The classical curvature measures related to the Normal Cycle are provided as well as the more general curvature measures based on corrected normal currents. This module implements the research on corrected normal currents presented in papers [71] and [72], and also gives a standalone implementation of the normal cycle method, which was presented in papers [35] and [36].
The following programs are related to this documentation: geometry/meshes/curvaturemeasuresicnc3d.cpp, geometry/meshes/curvaturemeasuresicncXY3d.cpp, geometry/meshes/curvaturemeasuresnc3d.cpp, geometry/meshes/curvaturemeasuresncXY3d.cpp, geometry/meshes/objcurvaturemeasuresicnc3d.cpp, geometry/meshes/objcurvaturemeasuresicncXY3d.cpp, geometry/meshes/volcurvaturemeasuresicnc3d.cpp, geometry/meshes/volcurvaturemeasuresicncXY3d.cpp, geometry/meshes/digpolycurvaturemeasurescnc3d.cpp, geometry/meshes/digpolycurvaturemeasurescncXY3d.cpp.
The following tests are related to this documentation: testCorrectedNormalCurrentComputer.cpp, testNormalCycleComputer.cpp.
Curvature is an important geometric piece of information which has the drawback of being well defined only on smooth sets. Unfortunately, most geometric data are nonsmooth (e.g. triangulated or polygonal surfaces, point clouds, polyhedra). A lot of works have thus aimed at generalizing curvatures to such sets, while guaranteeing that such definitions tend (in some sense) towards the classical smooth curvatures when working on finer and finer non smooth approximations of smooth data.
The seminal paper of Federer [50] first defined curvature measures on sets with positive reach (which includes smooth and convex polyhedra, but not triangulated surfaces in general). Using the notion of Normal Cycle introduced by Wintgen [115], the curvature measure principle was then extended to a wider class of objects including triangulations, digitized objects and subanalytic sets by Fu [56].
The main idea of the Normal Cycle is to embed the shape into the Grassmann bundle, which encodes both positions and the normal cones. This embedding is itself piecewise smooth, and curvatures can be defined onto it by integration with the invariant differential forms (LipschitzKilling forms). This approach was used for instance by CohenSteiner and Morvan [35] [36] for triangulated surfaces, where they show stability of curvature measures for surface sampling with convergent normals. Extensions of such works to point clouds using offset surfaces were proposed by Chazal, CohenSteiner, Lieutier and Thibert in [24].
However these methods do not provide consistent curvature estimates when the normal vector field of the discretization does not tend towards the normal vector field of the reference smooth surface. A famous example is the Schwarz lantern. This is also the case of digital surfaces that are boundaries of pixel/voxel data in images. For these surfaces, even high resolution data implies only 4 (in 2D) or 6 (in 3D) possible normal vectors.
This module therefore proposes an implementation of corrected normal currents to compute curvature measures onto quite general discrete surface data. The key idea is to replace the normal vector field of a surface S by another vector field u which we assume to be geometrically more meaningful. For instance if S is a digitization of a smooth surface X, one may take for u a local average of the naive normals of X.
The general theory was proposed by Lachaud, Romon and Thibert in [71]. It deals with piecewise \( C^{1,1} \) surfaces, with piecewise \( C^1 \) corrected normal vector field. A specialized version restricted to polygonal meshes equipped with a normal vector field defined by linear interpolation was presented in [72], and leads to quite simple formulae.
For instance, the interpolated corrected curvature measures take the following values on a triangle \( \tau_{ijk} \), with vertices i, j, k:
\[ \begin{align*} \mu^{(0)}(\tau_{ijk}) = &\frac{1}{2} \langle \bar{\mathbf{u}} \mid (\mathbf{x}_j  \mathbf{x}_i) \times (\mathbf{x}_k  \mathbf{x}_i) \rangle, \\ \mu^{(1)}(\tau_{ijk}) = &\frac{1}{2} \langle \bar{\mathbf{u}} \mid (\mathbf{u}_k  \mathbf{u}_j) \times \mathbf{x}_i + (\mathbf{u}_i  \mathbf{u}_k) \times \mathbf{x}_j + (\mathbf{u}_j  \mathbf{u}_i) \times \mathbf{x}_k \rangle, \\ \mu^{(2)}(\tau_{ijk}) = &\frac{1}{2} \langle \mathbf{u}_i \mid \mathbf{u}_j \times \mathbf{u}_k \rangle, \\ \mu^{\mathbf{X},\mathbf{Y}}(\tau_{ijk}) = & \frac{1}{2} \big\langle \bar{\mathbf{u}} \big \langle \mathbf{Y}  \mathbf{u}_k \mathbf{u}_i \rangle \mathbf{X} \times (\mathbf{x}_j  \mathbf{x}_i) \big\rangle \frac{1}{2} \big\langle \bar{\mathbf{u}} \big \langle \mathbf{Y}  \mathbf{u}_j \mathbf{u}_i \rangle \mathbf{X} \times (\mathbf{x}_k  \mathbf{x}_i) \big\rangle, \end{align*} \]
where \( \langle \cdot \mid \cdot \rangle \) denotes the usual scalar product, \( \bar{\mathbf{u}}=\frac{1}{3}( \mathbf{u}_i + \mathbf{u}_j + \mathbf{u}_k )\).
The measure \( \mu^{(0)} \) is the corrected area density of the given triangle, \( \mu^{(1)} \) is twice its corrected mean curvature density, \( \mu^{(2)} \) is its corrected Gaussian curvature density. The (anisotropic) measure \( \mu^{\mathbf{X},\mathbf{Y}} \) is the trace of the corrected second fundamental form along directions \( \mathbf{X} \) and \( \mathbf{Y} \). While the smooth second fundamental form is naturally a symmetric 2tensor, there is no easy way to define tangent directions at a vertex, so the anisotropic measure depends on two 3D vectors; when \( \mathbf{X} \) and \( \mathbf{Y} \) are tangent, \( \mu^{\mathbf{X},\mathbf{Y}} \) is close to the second fundamental form applied to these vectors, while its value along normal direction tends to zero asymptotically.
These measures are naturally extended to an arbitrary subset B of \( \mathbb{R}^3 \) by measuring its intersection ratio with each triangle \( \tau \) and by suming over all triangles of the mesh.
\[ \mu^{(k)}( B ) = \sum_{\tau : \text{triangle} } \mu^{(k)}( \tau ) \frac{\mathrm{Area}( \tau \cap B )}{\mathrm{Area}(\tau)}. \]
Curvature measures formula per triangle/edge/vertex are provided in class NormalCycleFormula (for the curvatures measures induced by the Normal Cycle) and in class CorrectedNormalCurrentFormula (for the curvatures measures induced by corrected normal currents with interpolated corrected normal vector field).
However one generally wishes to compute the measures over a mesh in a ball of given center and radius. It is then more convenient to use classes NormalCycleComputer and CorrectedNormalCurrentComputer, which are dedicated to compute curvature measures over a SurfaceMesh object.
You may proceed as follows (see examples geometry/meshes/curvaturemeasuresnc3d.cpp and geometry/meshes/curvaturemeasuresncXY3d.cpp).
You need to include header file of class NormalCycleComputer:
then use some typedefs adapted to your problem.
We will compute curvature measures onto a SurfaceMesh object. Here, we build a torus of big radius 3 and small radius 1, with a discretization of 20x20.
It is not necessary here to equip the mesh with normals, since the Normal Cycle infers the normals from the positions of the mesh. We can now compute three measures onto the mesh, mu0
the area measure, mu1
twice the mean curvature measure, and mu2
the Gaussian curvature measure. You have static methods NormalCycleComputer::meanCurvature and NormalCycleComputer::GaussianCurvature that estimate respectively mean and Gaussian curvature from measures.
Normal cycle mean and Gaussian curvature measures have meaning only around edges and vertices. So it is necessary to measure them in a big enough ball that includes at least one edge or vertex, and then to normalize the result by dividing by the area measure. Therefore, if we wish to estimate mean and Gaussian curvatures at every face centroid (say), the following snippet shows how to do it, assuming a measuring ball radius of R
:
This is is the result for a measuring ball radius of 0.5 onto the torus shape.
Expected mean curvatures: min=0.25 max=0.625 Computed mean curvatures: min=0.189446 max=0.772277 Expected Gaussian curvatures: min=0.5 max=0.25 Computed Gaussian curvatures: min=0.682996 max=0.547296
A quite similar code allows you to compute anisotropic curvature measures (i.e. a kind of second fundamental form). You may extract from the returned tensor measure (using eigenvalues/eigenvectors) principal curvatures and principal directions, by using method NormalCycleComputer::principalCurvatures.
Expected k1 curvatures: min=0.5 max=0.25 Computed k1 curvatures: min=0.581281 max=0.441977 Expected k2 curvatures: min=1 max=1 Computed k2 curvatures: min=0.904081 max=1.06404
Curvature estimates are displayed with a colormap [0.625 ... 0 ... 0.625] > [Blue ... White ... Red].
You may proceed as follows (see examples geometry/meshes/curvaturemeasuresicnc3d.cpp and geometry/meshes/curvaturemeasuresicncXY3d.cpp).
You need to include header file of class CorrectedNormalCurrentComputer:
then use some typedefs adapted to your problem.
We will compute curvature measures onto a SurfaceMesh object. Here, we build a torus of big radius 3 and small radius 1, with a discretization of 20x20.
Here we just provide the corrected normal vector field \( \mathbf{u}
\) as a normal vector per vertex, the field is then linearly interpolated. We can now compute three measures onto the mesh, mu0
the area measure, mu1
twice the mean curvature measure, and mu2
the Gaussian curvature measure. You have static methods CorrectedNormalCurrentComputer::meanCurvature and CorrectedNormalCurrentComputer::GaussianCurvature that estimate respectively mean and Gaussian curvature from measures.
A big advantage of (interpolated) Corrected Normal Current over Normal Cycle is that the induced mean and Gaussian curvature measures have meaning for arbitrary measuring set (even when the area of the set tends to zero. Therefore, if we wish to estimate mean and Gaussian curvatures at every face centroid (say), the following snippet shows how to do it, assuming a measuring ball radius of R
:
This is is the result for a measuring ball radius of 0.0 onto the torus shape (perfect results !).
Expected mean curvatures: min=0.25 max=0.625 Computed mean curvatures: min=0.25 max=0.625 Expected Gaussian curvatures: min=0.5 max=0.25 Computed Gaussian curvatures: min=0.5 max=0.25
This is is the result for a measuring ball radius of 0.5 onto the torus shape (better than Normal cycle).
Expected mean curvatures: min=0.25 max=0.625 Computed mean curvatures: min=0.264763 max=0.622318 Expected Gaussian curvatures: min=0.5 max=0.25 Computed Gaussian curvatures: min=0.470473 max=0.244636
A quite similar code allows you to compute anisotropic curvature measures (i.e. a kind of second fundamental form). You may extract from the returned tensor measure (using eigenvalues/eigenvectors) principal curvatures and principal directions, by using method CorrectedNormalCurrentComputer::principalCurvatures.
This gives you such results for principal curvatures estimations, first for a radius 0:
Expected k1 curvatures: min=0.5 max=0.25 Computed k1 curvatures: min=0.500225 max=0.249888 Expected k2 curvatures: min=1 max=1 Computed k2 curvatures: min=1.00011 max=1.00678
then for a radius of 0.5:
Expected k1 curvatures: min=0.5 max=0.25 Computed k1 curvatures: min=0.454026 max=0.242436 Expected k2 curvatures: min=1 max=1 Computed k2 curvatures: min=0.924283 max=0.95338
Curvature estimates are displayed with a colormap [0.625 ... 0 ... 0.625] > [Blue ... White ... Red].
./examples/geometry/meshes/curvaturemeasuresnc3d lantern 20 20 0.5outputs Expected mean curvatures: min=0.25 max=0.25 Computed mean curvatures: min=0.795695 max=1.41211 Expected Gaussian curvatures: min=0 max=0 Computed Gaussian curvatures: min=6.79045e14 max=15.0937  ./examples/geometry/meshes/curvaturemeasuresicnc3d lantern 20 20 0.5outputs Expected mean curvatures: min=0.25 max=0.25 Computed mean curvatures: min=0.25 max=0.25 Expected Gaussian curvatures: min=0 max=0 Computed Gaussian curvatures: min=0 max=0 
The theory of corrected normal currents works for arbitrary piecewise smooth corrected normal vector fields. Class CorrectedNormalCurrentComputer allows you to choose between a constant per face corrected normal vector field (called Constant Corrected Normal Current, CCNC), or a smooth corrected normal vector field obtained by linear interpolation of vertex normals per face (called Interpolated Corrected Normal Current, ICNC).
The choice is made automatically by the CorrectedNormalCurrentComputer object according to the userdefined normal vector field of its associated SurfaceMesh object:
myMesh.vertexNormals()
is not empty, then the CorrectedNormalCurrentComputer computes ICNC measures,myMesh.faceNormals()
is not empty, then the CorrectedNormalCurrentComputer computes CCNC measures,For instance, you may force ICNC on a mesh as follows:
Otherwise, assume you have not set any vertex normals, the following code uses CCNC:
NormalCycleComputer and CorrectedNormalCurrentComputer provides methods that return measures, more precisely as SurfaceMeshMeasure objects. The values of measure may be arbitrary scalars, vectors or tensors. Here curvature measures are either scalars (like area, mean and Gaussian curvature measures) or tensors (anisotropic measures).
A SurfaceMeshMeasure provides several handy methods to evaluate measures wrt simple subsets of \( \mathbb{R}^3 \).
Examples geometry/meshes/objcurvaturemeasuresicnc3d.cpp and geometry/meshes/objcurvaturemeasuresicncXY3d.cpp show how to compute all curvature information on an arbitrary mesh.
The only difference in the above examples is the way the SurfaceMesh object is created. Here we read an OBJ file to create it.
Remember that a SurfaceMesh object may represent a nonmanifold polygonal surface.
Once the mesh is created, you just need to ensure that the SurfaceMesh object has a normal vector attached to each of its vertex. This can be done as follows (or just provide an OBJ file with vertex normals):
The rest of the code is identical. The two examples gives the following curvature estimates on "spot.obj" file.
It is difficult to estimate curvatures on digital surfaces, since they naturally have only very few normal directions. You may have a look at module Integral invariant curvature estimator 2D/3D to see a method that estimates such curvatures when the digital surface is the boundary of a volume of voxels.
Corrected Normal Currents are well adapted to digital surface geometric analysis. Indeed, we can have good estimates of the normal vector field of a digital surface (see ShortcutsGeometry), which will be used as the corrected normal vector field associated to the current.
As examples of digital surfaces, we show how to estimate curvatures on a digitization of polynomial surfaces, then on the boundary of a voxel object defined in VOL file.
Examples geometry/meshes/digpolycurvaturemeasurescnc3d.cpp and geometry/meshes/digpolycurvaturemeasurescncXY3d.cpp show how to extract all curvature information from digitized polynomial surfaces.
We first read the polynomial definition and extracts the digital surface approximating the implicit surface at the given resolution.
We build a SurfaceMesh object that represents exactly the digital surface.
We build the Corrected Normal Current on this mesh, and use a digital estimator of normal vector (here convolved trivial normals) to set the corrected normal vector at each face or vertex, depending on the choice of the user.
The remaining of the code to estimate curvatures from measures is identical to above. You can obtain the following results:
Examples geometry/meshes/volcurvaturemeasuresicnc3d.cpp and geometry/meshes/volcurvaturemeasuresicncXY3d.cpp show how to estimate curvatures on the boundary of a digital object.
We first read the VOL file and extracts the digital surface bounding the digital object.
We build a SurfaceMesh object that represents exactly the digital surface.
We build the Corrected Normal Current on this mesh, and use a digital estimator of the normal vectors (here convolved trivial normals) to set the corrected normal vector at each vertex.
The remaining of the code to estimate curvatures from measures is identical to above. We run the examples on "Al" vol file, which gives the following results.
We give below a few comparisons (in terms of accuracy and speed) between CNC curvature estimations and Integral Invariant (II) method (see Integral invariant curvature estimator 2D/3D). II method is a curvature estimation method specialized (and restricted) to boundaries of digital volumes (i.e. set of voxels). Therefore we compare their respective performance on digitizations of polynomial surfaces.
We run the tests on goursat polynomial surface, with digitization gridsteps \( h \in \{ 1, 0.5, 0.25, 0.125, 0.0625 \} \). Parameters are set up as follows:
We plot below accuracy as a function of the digitization gridstep.
We plot below accuracy as a function of the computation time (in ms).
curvature estimator  number of surfels  error loo  error l2  total time (ms)  normal estimation time (ms) 

II mean curvature  9510  0.212016  0.0835415  148 (mean+Gaussian)  
ICNC mean curvature  9510  0.0640707  0.0199467  450 (mean+Gaussian)  83 
II Gaussian curvature  9510  0.0501252  0.0150247  148 (mean+Gaussian)  
ICNC Gaussian curvature  9510  0.0263134  0.00686261  450 (mean+Gaussian)  83 
II mean curvature  151374  0.124552  0.0546997  31635 (mean+Gaussian)  
ICNC mean curvature  151374  0.0295432  0.00730423  38582 (mean+Gaussian)  19322 
II Gaussian curvature  151374  0.0433595  0.00983477  31635 (mean+Gaussian)  
ICNC Gaussian curvature  151374  0.0194384  0.00262745  38582 (mean+Gaussian)  19322 