DGtal  1.3.beta
Digital surface regularization using the Shrouds algorithm
Author(s) of this documentation:
Jacques-Olivier Lachaud

Part of the Geometry package.

This part of the manual describes classes and functions related to the regularization of digital surfaces using the "Shrouds" algorithm, by Nielson et al [88].

On binary 3D images, the "Shrouds" algorithm builds a smooth surface that separates the interior voxels from the exterior voxels. It starts from the (closed) digital surface that separates interior voxels from exterior voxels. The 3D smooth surface reconstruction problem is transformed into a series of 2D smooth contour reconstruction problems, by slicing the 3D volume (and thus the digital surface) with axis-align planes. Each contour is regularized to minimize its squared curvature and the separation property is obtained by forcing each vertex to lie on its unit edge joining an interior voxel to an exterior voxel (i.e. the edge that is dual to the surfel). The 2D problems are intertwined since each vertex is shared by two slices.

Note
This method is limited to closed surface reconstruction, in opposition to the method presented in Digital surface regularization by normal vector alignment.

Disclaimer. This implementation is our interpretation of the work presented in [88], and tries to follow at best the ideas presented therein. However, a major part of the numerical resolution is not detailed. We have inferred a reasonable way of solving the fourth-order equations but results may differ from the original implementation. For instance, visually on the bunny example of the paper, our results are very close (and look better). This implementation comes thus without any guarantee.

Related example: geometry/surfaces/testShroudsRegularization.cpp

Introduction

Let I be a binary image and S the digital surface that is the boundary of the object in I. Slicing S by a plane aligned orthogonal to an axis k defines a ring of consecutive surfels, that can be seen as a simple contour joining points \( X[s] \) associated to each surfel \( s \in S \). If we define a meaningful adjacency on the digital surface (see Digital surfaces), when joining points \( X[\cdot] \) corresponding adjacent surfels within the slice, we build closed 2D contours. Since each surfel in 3D has only two tangent directions, each surfel belongs to exactly two orthogonal contours. During the regularization, the "Shrouds" algorithm principle is to move the point \( X[s] \) associated to each surfel s in a constrained way : it can only move on the unit segment that is orthogonal to the surfel, otherwise said the point is constrained to lie on the edge dual to the primal surfel.

The "Shrouds" algorithm minimizes a global energy that is determined by the positions of each point on its two contours. The user can choose between minimizing the area of the output surface (ShroudsRegularization::Regularization::AREA), the snake energy of the output surface (ShroudsRegularization::Regularization::SNAKE) or the squared curvature of the output surface (ShroudsRegularization::Regularization::SQUARED_CURVATURE). The latter one gives the best results.

The snake energy \( E^{snk} \) and the squared curvature energy \( E^{\kappa^2} \) are defined per slice on each 2D contour.

For a slice contour \( C=(x(s),y(s)) \) and boundary constraints, we define:

\[ E^{snk}(C) = \int_C \alpha (x'(s)^2 + y'(s)^2) + \beta (x''(s)^2 + y''(s)^2) ds, \]

and

\[ E^{\kappa^2}(C) = \int_C (x'(s) y''(s) + x''(s) y'(s))^2 / (x'(s)^2 + y'(s)^2)^3 ds, \]

Note
This process does not optimize the true 3D squared mean curvature of the surface or its 3D elastic and thin plate energy.

Apart from area energy, energies \( E^{snk}(C) \) and \( E^{\kappa^2}(C) \) are not convex and movement constraints on \( X \) make things worse. We propose here dedicated finite difference algorithms giving a reasonnable gradient descent for the regularization process, together with randomization at the beginning.

Main usage with squared curvature energy

Starting from an implicit digital surface (with a gridstep set to 0.3 instead of 1 on images):

typedef Shortcuts<Z3i::KSpace> SH3;
typedef SH3::ExplicitSurfaceContainer Container;
auto params = SH3::defaultParameters();
params( "polynomial", "goursat" )( "gridstep", 1)("verbose", 0);
auto implicit_shape = SH3::makeImplicitShape3D ( params );
auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
auto K = SH3::getKSpace( params );
auto surface = SH3::makeDigitalSurface( digitized_shape, K, params );

we obtain the following slice geometry, best seen on the polygonal surface dual to the digital surface:

Initialization for the Shrouds algorithm (dual surface to goursat digital surface).

The regularization class instance can be set up with default parameters using the following syntax:

auto idxsurface = SH3::makeIdxDigitalSurface( surface, params );
ShroudsRegularization< Container > shrouds_reg( idxsurface );
auto originalPos = shrouds_reg.positions();

You may use ShroudsRegularization::setParams to change some parameters:

  • \( \epsilon \) the bounds for varying the positions of points on dual edge in \( [\epsilon,1-\epsilon] \), where 0 is the center of the bordering inner voxel and 1 is the center of the bordering outer voxel.
  • \( \alpha \) parameter for Snake first order regularization (~ area).
  • \( \beta \) parameter for Snake second order regularization (~ curvature).

We can now minimize the squared curvature energy as follows:

double loo = 0.0;
double l2 = 0.0;
double energyInitK2 = shrouds_reg.energy ( RegType::SQUARED_CURVATURE );
std::tie( loo, l2 ) = shrouds_reg.regularize( RegType::SQUARED_CURVATURE,
0.5, 0.0001, 100 );
double energyRegK2 = shrouds_reg.energy ( RegType::SQUARED_CURVATURE );

The user can specify the number of steps of the gradient descent as well as the randomization. If the ShroudsRegularization::regularize method is called another time, the descent starts from the previous results (aka warm restart).

Using the default settings, we obtain the following reconstruction, after respectively 100 and 1000 iterations:

Optimized surface according to squared curvature (100 iterations).
Optimized surface according to squared curvature (1000 iterations).

Clearly, the result after 1000 iterations is more satisfactory. You may output the result for instance as:

auto regularizedPos = shrouds_reg.positions();
auto polySurf = SH3::makeDualPolygonalSurface( idxsurface );
auto polySurfPos = polySurf->positions();
for ( auto i = 0; i < regularizedPos.size(); i++ )
polySurfPos[ i ] = regularizedPos[ i ];
SH3::saveOBJ( polySurf, "goursat-shrouds-reg-k2.obj" );
Note
The regularized position are retrieved using the ShroudsRegularization::positions. The current energy can be retrieved with ShroudsRegularization::energy.

Shrouds regularization with area and snake energy

You may optimize according the shape according to area minimization as follows:

shrouds_reg.init();
double energyInitArea= shrouds_reg.energy ( RegType::AREA );
std::tie( loo, l2 ) = shrouds_reg.regularize( RegType::AREA,
0.5, 0.0001, 100 );
double energyRegArea = shrouds_reg.energy ( RegType::AREA );

and according to snake energy as follows (use ShroudsRegularization::setParams to change \( \alpha,\beta \) parameters):

shrouds_reg.init();
double energyInitSnk= shrouds_reg.energy ( RegType::SNAKE );
std::tie( loo, l2 ) = shrouds_reg.regularize( RegType::SNAKE,
0.5, 0.0001, 100 );
double energyRegSnk = shrouds_reg.energy ( RegType::SNAKE );

We display below the three possible optimizations (after 1000 iterations), which shows that squared curvature energy \( E^{\kappa^2}(C) \) produces the best output.

Optimized surface according to squared curvature.
Optimized surface according to snake energy.
Optimized surface according to area.
Note
Results are visually very good, notably due to the nice slice / isocontouring effect. However, the output polygonal surface contains many cells with bad geometric properties (very small or elongated faces).
SH3
Shortcuts< KSpace > SH3
Definition: testArithmeticalDSSComputerOnSurfels.cpp:49
K
KSpace K
Definition: testCubicalComplex.cpp:62
DGtal::Shortcuts::makeDigitalSurface
static CountedPtr< DigitalSurface > makeDigitalSurface(CountedPtr< TPointPredicate > bimage, const KSpace &K, const Parameters &params=parametersDigitalSurface())
Definition: Shortcuts.h:1209
DGtal::Shortcuts::defaultParameters
static Parameters defaultParameters()
Definition: Shortcuts.h:203
DGtal::Shortcuts::makeImplicitShape3D
static CountedPtr< ImplicitShape3D > makeImplicitShape3D(const Parameters &params=parametersImplicitShape3D())
Definition: Shortcuts.h:282
DGtal::Shortcuts::makeIdxDigitalSurface
static CountedPtr< IdxDigitalSurface > makeIdxDigitalSurface(CountedPtr< BinaryImage > bimage, const KSpace &K, const Parameters &params=parametersDigitalSurface())
Definition: Shortcuts.h:1275
DGtal::Shortcuts::makeDigitizedImplicitShape3D
static CountedPtr< DigitizedImplicitShape3D > makeDigitizedImplicitShape3D(CountedPtr< ImplicitShape3D > shape, Parameters params=parametersDigitizedImplicitShape3D())
Definition: Shortcuts.h:523
DGtal::Shortcuts::getKSpace
static KSpace getKSpace(const Point &low, const Point &up, Parameters params=parametersKSpace())
Definition: Shortcuts.h:332
DGtal::SetOfSurfels
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as connected surfels....
Definition: SetOfSurfels.h:73
DGtal::ShroudsRegularization::Regularization
Regularization
The enum class specifying the possible shrouds regularization.
Definition: ShroudsRegularization.h:96
DGtal::Shortcuts::saveOBJ
static bool saveOBJ(CountedPtr< ::DGtal::DigitalSurface< TDigitalSurfaceContainer > > digsurf, const TCellEmbedder &embedder, const RealVectors &normals, const Colors &diffuse_colors, std::string objfile, const Color &ambient_color=Color(32, 32, 32), const Color &diffuse_color=Color(200, 200, 255), const Color &specular_color=Color::White)
Definition: Shortcuts.h:1739
DGtal::Shortcuts::makeDualPolygonalSurface
static CountedPtr< PolygonalSurface > makeDualPolygonalSurface(Surfel2Index &s2i, CountedPtr< ::DGtal::DigitalSurface< TContainer > > aSurface)
Definition: Shortcuts.h:2268