Tutorial: Surface area estimation

Table of Contents

• Introduction
• Import volumetric image an visualization
• Border extraction
• Normal estimation
• Conclusion

Author

Jeremy Levallois

The objective of this tutorial is to perform some elementary surface analysis of 3D objects. With this tutorial, you will experiment the following DGtal features

• Build a 3D shape (parametric)
• Convert it to a Digital Shape (GaussDigitizer)
• Extract Digital border of a digital shape (LightImplicitDigitalSurface / DigitalSurface)
• Use an differential quantity estimator (IntegralInvariantCovarianceEstimator + a functor)

Please use the following template file for this tutorial: AreaSurfaceEstimation-template.cpp

Introduction

When analyzing 3D border of a shape, the surface area plays an important role for shape description. A way to estimate the area of a surface is to compute the difference between the normal of the shape and the normal of the digital point. Given a shape \( \mathcal{X}\), the surface area function can be estimated with :

\[ \hat{AreaSurface (X)} = \sum_x ( \textbf{n_e}(x) . \textbf{n_s}(\hat{x}) ) \]

where \( \textbf{n_e}(x) \) denotes the normal vector estimated at position \( x \in \mathcal{X} \), \( \textbf{n_s}(\hat{x}) \) the normal vector of surface element (surfel) \( \hat{x} \in Dig_h{\mathcal{X}} \), and \( . \) the dot product.

Import volumetric image an visualization

Let start with shape construction. First of all, create a source file with the following includes (or copy the AreaSurfaceEstimation-template.cpp source file):

 
#include "DGtal/shapes/GaussDigitizer.h"
#include "DGtal/topology/LightImplicitDigitalSurface.h"
#include "DGtal/topology/DigitalSurface.h"
#include "DGtal/graph/DepthFirstVisitor.h"
#include "DGtal/graph/GraphVisitorRange.h"
 
#include "DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h"
#include "DGtal/geometry/surfaces/estimation/IntegralInvariantCovarianceEstimator.h"
 
using namespace DGtal;

1. Define a type (typedef) for the shape you want to use (Ball3D, ImplicitBall, ....)
2. Use the GaussDigitizer to convert the Euclidean shape to a digital one, with a given gridstep \( h \).

Once you have created the digital shape, you need to get the domain and create a KhalimskySpace (see Z3i::KSpace) from it.

At this point, you have a digital object created from an Euclidean one. Now, we want to extract the border of this digital object.

Border extraction

1. First of all, have a look to the LightDigitalSurface and DigitalSurface class and check the template type you have to specified to this class to create an instance.
2. As you may have seen, we need a visitor to iterate over the border of the digital Surface. In order to use the best part of the following estimator, we will use a depth first visitor (see DepthFirstVisitor and GraphVisitorRange).
3. Thanks to the GraphVisitor, we can iterate over the digital surface.

Normal estimation

As discussed in the introduction, we will use a normal vector estimator to compute the surface area of our object. We will use IntegralInvarianceCovarianceEstimator (II) who, with the good functor (can be found in IIGeometricFunctors namespace/file ). DGtal offers also some other normal vector estimators, so feel free to use wich one you want.

The main idea is to choose the functor you want (IIGeometricFunctors::IINormalDirectionFunctor on this case) and use the II estimator on all the digital surface:

• Initialize the functor
• Initialize the covariance estimator ( constructor(functor), attach(), setParams(), init() )
• Evaluate the estimator on all surfel of the border of the shape
• Compute for all surfel the dot product of the estimated euclidean normal vector (result of the estimator) and the surfel normal vector (see below)

To extract the surfel normal vector, we need to use some method from the KSpace (see KhalimskySpaceND for all methods avalaible):

• sOrthDir( surfel ) returns the orthogonal direction of a surfel.
• sDirectIncident( surfel, orthDir ) returns the direct incident surfel of a given surfel and a direction
• sKCoords( surfel ) returns a vector of coordinates (of KhalimskySpace)

1. Is the surface area estimation multigrid convergent ? As you can see at the begining of the file, we set the gridstep \( h \) to

1. If you change this value \( h \le 1 \), how the estimation of the surface area changes ?

Conclusion

You can check AreaSurfaceEstimation-final.cpp to get some answers to the above mentioned questions.

Several optimizations can be developed. For example, instead of evaluate for all surface element separately, we can estimate all surfels in one call of eval(). This uses some optimization of II estimator and can be see in overall computation time. Some others normal vector estimators can be used, as VoronoiCovarianceMeasureLocalEstimator.