DGtal  0.9.4beta
Integral invariant curvature estimator 2D/3D

Table of Contents

Author
Jérémy Levallois

Overview

The algorithm implemented in the class IntegralInvariantVolumeEstimator and IntegralInvariantCovarianceEstimator are detailed in the article [21] .

In geometry processing, interesting mathematical tools have been developed to design differential estimators on smooth surfaces based on Integral Invariant [59] [60] . The principle is simple: we move a convolution kernel along the shape surface and we compute integrals on the intersection between the shape and the convolution kernel, as follow in dimension 3:

\[ V_r(x) = \int_{B_r(x)} \chi(p)dp\ \]

where \( B_r(x) \) is the Euclidean ball of radius \( r \), centered at \( x \) and \( \chi(p) \) the characteristic function of \( X \). In dimension 2, we simply denote \( A_r(x) \) such quantity (represented in orange color on the following figure).

integral2D.png
Integral invariant computation in dimension 2.
integral3D.png
Integral invariant computation in dimension 3.
notations_integral.png
Notations.

Integral Invariant for curvature computation

In [20] , we have demonstrated that some digital integral quantities provide curvature information when the kernel size tends to zero for a sufficiently smooth shape. Indeed, at \( x \) of the surface \( X \) and with a fixed radius \( r \), we obtain convergent local curvature estimators \( \tilde{\kappa}_r(X,x) \) and \( \tilde{H}_r(X,x) \) of quantities \( \kappa(X,x) \) and \( H(X,x) \) respectively:

\[ \tilde{\kappa}_r(X,x) = \frac{3\pi}{2r} - \frac{3A_r(x)}{r^3},\quad \tilde{\kappa}_r(X,x) = \kappa(X,x)+ O(r) \]

\[ \tilde{H}_r(X,x) = \frac{8}{3r} - \frac{4V_r(x)}{\pi r^4},\quad \tilde{H}_r(X,x) = H(X,x) + O(r) \]

where \( \kappa_r(X,x) \) is the 2d curvature of \( X \) at \( x \) and \( H_r(X,x) \) is the 3d mean curvature of \( X \) at \( x \).

Then we showed that we can obtain local digital curvature estimators :

\[ \forall 0 < h < r,\quad \hat{\kappa}_r(Z,x,h) = \frac{3\pi}{2r} - \frac{3\widehat{Area}(B_{r/h}(\frac{1}{h} \cdot x ) \cap Z, h)}{r^3} \]

where \( \hat{\kappa}_r \) is an integral digital curvature estimator of a digital shape \( Z \subset ℤ^2 \) at point \( x \in \rm I\! R^2 \) and step \( h \). \( B_{r/h}(\frac{1}{h} \cdot x ) \cap Z, h) \) means the intersection between \( Z \) and a Ball \( B \) of radius \( r \) digitized by \( h \) centered in \( x \).

In the same way, we have in 3d :

\[ \forall 0 < h < r,\quad \hat{H}_r(Z',x,h) = \frac{8}{3r} - \frac{4\widehat{Vol}(B_{r/h}(\frac{1}{h} \cdot x ) \cap Z', h)}{\pi r^4}. \]

where \( \hat{H}_r \) is an integral digital mean curvature estimator of a digital shape \( Z' \subset ℤ^3 \) at point \( x \in \rm I\! R^3 \) and step \( h \) .

We have demonstrated in [20] that to prove the multigrid convergence with a convergence speed of \( O(h^{\frac{1}{3}}) \), the Euclidean radius of the kernel must follow the rule \( r = k_mh^{\alpha_m} \) ( \( \alpha_m = \frac{1}{3} \) provides better worst-case errors, so we will use this value).

Experimental results can be found at https://liris.cnrs.fr/jeremy.levallois/Papers/DGCI2013/ and https://liris.cnrs.fr/jeremy.levallois/Papers/AFIG2013/

Algorithm

Overall algorithm

The user part is rather simple by using only IntegralInvariantVolumeEstimator or IntegralInvariantCovarianceEstimator.

Since 2d curvature and 3d mean curvature are related to the volume between the shape and a ball centered on the border of the shape, we need to estimate this volume and use this value to get the curvature information. Then we parametrize our volume estimator IntegralInvariantVolumeEstimator by a functor to compute the curvature (in 2d: functors::IICurvatureFunctor and in 3d for the mean curvature: functors::IIMeanCurvature3DFunctor )

In 3d we can also compute the covariance matrix of the intersection between the shape and the ball centered on the border of the shape using IntegralInvariantCovarianceEstimator. This offers us the possibility to extract many differential quantities as: 3d Gaussian curvature, first and second principal curvatures, as well first and second principal curvature directions, normal vector or tangeant vector (functors::IIGaussianCurvature3DFunctor, etc. See file DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h or namespace functors for the entire list).

Integral Invariant functors

All functors defined in IIGeometricFunctors.h have the same inititialization process: you need to use the default constructor, and initialize them with init() method. They will ask for a grid step ( \( h \)) and the Euclidean radius of the kernel ball ( \( r \)).

Integral Invariant estimators

As described below, this parameter \( r_e \) determines the level of feature detected by the estimator Since IntegralInvariantVolumeEstimator and (resp.) IntegralInvariantCovarianceEstimator are models of CSurfelLocalEstimator, they will follow rules inferred by it (init(), eval(), etc.). But some particular considerations are introduce. The following steps are the general usage of Integral Invariant estimators:

If you want to evaluate on a range of surfels, we recommend you to choose the second way. Indeed, some optimizations are available when the range of surfels is given with a 0-adjacency. Then, when you extract the digital surface of your shape, it's recommended to use a depth-first search (see DepthFirstVisitor for example). If none, no optimization are perform (it will be visible in performances for big shape).

Example code

It is important to consider a range of connected surfels when evaluating with Integral Invariant Curvature estimators in order to benefit the kernel optimization. Note that the methodology is the same with both IntegralInvariantVolumeEstimator and IntegralInvariantCovarianceEstimator. The only change is on the typedef of Functor/Estimator (see below):

double radius = std::atof(argv[3]);
typedef functors::IIMeanCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
typedef IntegralInvariantVolumeEstimator< Z3i::KSpace, ImagePredicate, MyIICurvatureFunctor > MyIICurvatureEstimator;
// For computing Gaussian curvature instead, for example, change the two typedef above by :
// typedef functors::IIGaussianCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
// typedef IntegralInvariantCovarianceEstimator< Z3i::KSpace, ImagePredicate, MyIICurvatureFunctor > MyIICurvatureEstimator;
// and it's done. The following part is exactly the same.
typedef MyIICurvatureFunctor::Value Value;
MyIICurvatureFunctor curvatureFunctor; // Functor used to convert volume -> curvature
curvatureFunctor.init( h, radius ); // Initialisation for a grid step and a given Euclidean radius of convolution kernel
MyIICurvatureEstimator curvatureEstimator( curvatureFunctor );
curvatureEstimator.attach( KSpaceShape, predicate ); // Setting a KSpace and a predicate on the object to evaluate
curvatureEstimator.setParams( radius / h ); // Setting the digital radius of the convolution kernel
curvatureEstimator.init( h, abegin, aend ); // Initialisation for a given h, and a range of surfels
std::vector< Value > results;
std::back_insert_iterator< std::vector< Value > > resultsIt( results ); // output iterator for results of Integral Invariant curvature computation
curvatureEstimator.eval( abegin, aend, resultsIt ); // Computation

Some results

Here is some results in 2d and 3d :

example-integralinvariant2D.png
Curvature mapped on a Flower2D
snapshot-mean-zero.png
Mean curvature mapped on a Goursat's surface
snapshot-K-zero.png
Gaussian curvature mapped on a Goursat's surface
Bunny_128_mean.png
Mean curvature mapped on a Stanford bunny (credit is given to the Stanford Computer Graphics Laboratory https://graphics.stanford.edu/data/3Dscanrep/)
Bunny_64_k1.png
First principal curvature direction mapped on a Stanford bunny (credit is given to the Stanford Computer Graphics Laboratory https://graphics.stanford.edu/data/3Dscanrep/)