DGtal  0.9.3beta
Tutorial: Local morphological opening of a 3d volume by the Fast Marching Method

Table of Contents

Author(s) of this documentation:
Tristan Roussillon

Goal of the tutorial

In this tutorial, we will see how to use the FMM class in order to perform a basic processing on a 3d volume. The FMM class is part of the Geometry package. In the same time, we will learn how to read and write vol files (see IO package) and we will get familiar with digital surfaces (see Topology package).

In order to smooth the digital surface of a 3d volume, we want to perform a mophological opening (i.e. an erosion followed by a dilation) with a circular structural element of small radius. Since the radius is small, we do not want to scan the whole image but we want to work only in narrow band around the digital surface. That is why, we think that using the FMM class is a good idea.

Some steps

We will detail the implementation of some crucial steps:

The whole example is available in FMMErosion.cpp.

Reading a vol

First of all, we choose a type for the binary image representing the 3d volume. Then, we create the image and fill it by reading a vol file with VolReader.

string imageFileName = examplesPath + "samples/cat10.vol";
trace.info() << imageFileName <<std::endl;
typedef ImageContainerBySTLVector<Domain, bool> BinaryImage;
BinaryImage binaryImage = VolReader<BinaryImage>::importVol( imageFileName);

See also Image and digital object import/export.

Tracking a surface

We now track the digital surface of the volume in order to work in a narrow band around it.

This can be done in two steps:

KSpace ks;
Domain domain = binaryImage.domain();
ks.init( domain.lowerBound(), domain.upperBound(), true );
KSpace::SCell bel;
try {
//getting a bel
bel = Surfaces<KSpace>::findABel( ks, binaryImage, domain.size() );
trace.info() << "starting bel: " << bel << std::endl;
} catch (DGtal::InputException i) {
trace.emphase() << "starting bel not found" << std::endl;
return 0;
}
typedef functors::FrontierPredicate<KSpace, BinaryImage> SurfelPredicate;
typedef LightExplicitDigitalSurface<KSpace, SurfelPredicate> Frontier;
functors::SCellToIncidentPoints<KSpace> toIncidentPoints( ks );
std::pair<Point,Point> bpair = toIncidentPoints( bel );
SurfelPredicate surfelPredicate( ks, binaryImage,
binaryImage( bpair.first ),
binaryImage( bpair.second ) );
Frontier frontier( ks, surfelPredicate,
SurfelAdjacency<KSpace::dimension>( true ), bel );
Note
In this snippet, we do not explicitly track the digital surface. The digital surface is light, i.e. it does not contain any surfel, but computes them on-line.

See also Digital surfaces and Helpers for digital surfaces.

Local erosion

In order to perform a local morphological erosion, the idea is to incrementally compute the band of points to remove by marching out from the digital surface, but only inside the initial volume. We stop as soon as we detect a point whose (Euclidean) distance from the digital surface is greater than a given threshold (called maximalDistance below), which corresponds to the radius of the structuring element.

As explained in nD Fast Marching Methods, we must define

typedef ImageContainerBySTLMap<Domain,double> DistanceImage;
typedef DigitalSetFromMap<DistanceImage> AcceptedPointSet;
typedef FMM<DistanceImage, AcceptedPointSet, BinaryImage > FMM;

Then, we construct and initalize the two external data structures, i.e. the set of points and the image of distance values. We arbitrarily set the distance values of all the inner points to \( 0.5 \).

DistanceImage imageDistance( domain, 0.0 );
AcceptedPointSet pointSet( imageDistance );
functors::SCellToInnerPoint<KSpace> toInnerPoint( ks );
for (Frontier::SurfelConstIterator it = frontier.begin(), itEnd = frontier.end();
it != itEnd; ++it)
{
pointSet.insert( toInnerPoint(*it) );
imageDistance.setValue( toInnerPoint(*it), 0.5 );
}

Finally, we instanciate the FMM class and perform the computation as follows:

FMM fmm( imageDistance, pointSet, binaryImage,
domain.size(), maximalDistance );
fmm.compute();
trace.info() << fmm << std::endl;

Note that we set the last (optional) argument to maximalDistance in order to bound the computation. Consequently, it remains to flip from the volume to the background all the points belonging to pointSet. Indeed, they are all located inside the volume, at a distance from the digital surface that is less than maximalDistance by definition.

for (AcceptedPointSet::ConstIterator it = pointSet.begin(), itEnd = pointSet.end();
it != itEnd; ++it)
{
binaryImage.setValue(*it, 0);
}

Writing a vol

In order to export the updated binary image in a vol file, we use VolWriter as follows:

string outputFileName = "eroded.vol";
trace.info() << "to " << outputFileName << std::endl;
VolWriter<BinaryImage,functors::Cast<unsigned char> >::exportVol(outputFileName, binaryImage);
Note
In order to visualize a vol file with QGLViewer, you can use the tool 3dVolViewer of the DGtalTools project (https://github.com/DGtal-team/DGtalTools).

Your challenge

You know how to perform a local morphological erosion. All the previous snippets are taken from FMMErosion.cpp. The local morphological dilation as well as the final morphological opening, which is an erosion followed by a dilation, are left as an exercise.

Excercise

Your challenge is to code:

  1. a local morphological dilation,
  2. a local morphological opening,
  3. a tool that performs a local morphological opening on all connected components. The vol file and the radius of the structural element must be program options.

Hints

The dilation operation will look like the erosion operation. There are however three small differences:

using namespace DGtal::functors;
typedef NotPointPredicate<BinaryImage> BackgroundPredicate;
BackgroundPredicate,

Once these types are defined, you can construct your new point predicate as follows:

BackgroundPredicate backgroundPredicate( binaryImage );
DomainPredicate<Domain> domainPredicate( domain );
Predicate pred( backgroundPredicate, domainPredicate, andBF2 );

Note that most point predicate adapters are defined in BasicPointPredicates.h. Moreover, functors::andBF2 is a binary functor defined in BasicBoolFunctors.h