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:

reading a vol image.
tracking a digital surface.
performing a local erosion by FMM.
writing a vol image.

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);

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:

We search for a boundary element (bel for short), i.e. a surfel lying between the volume and its background.

  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, (unsigned int)domain.size() );
    trace.info() << "starting bel: " << bel << std::endl;
 
  } catch (DGtal::InputException& i) {
    trace.emphase() << "starting bel not found" << std::endl; 
    return 0; 
  }

We define and track the digital surface from the starting bel.

  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.

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

a type for the set of points for which we compute the distance values (called AcceptedPointSet below),
a type for the image that stores the distance values (called DistanceImage below),
a type for the point predicate that implicitly represents the volume. Here, we use the binary image as a point predicate.

  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:

a local morphological dilation,
a local morphological opening,
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:

Contrary to the erosion operation that is bounded by the volume, the dilation operation must be limited to the background, but inside the domain. On one hand, you may adapt the binary image, used as a point predicate in FMMErosion.cpp, with functors::NotPointPredicate in order to bound the computation to the background. On the other hand, you may adapt the domain with functors::DomainPredicate in order to limit the computation to the domain. Finally, you may use functors::BinaryPointPredicate to combine these two point predicates. Thus, you may define the following types:

using namespace DGtal::functors;
typedef NotPointPredicate<BinaryImage> BackgroundPredicate; 
typedef BinaryPointPredicate<
  BackgroundPredicate, 
  DomainPredicate<Domain> > Predicate;            

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

You must initialize the accepted point set from the outer points of the digital surface. You may use functors::SCellToOuterPoint defined in SCellsFunctors.h instead of functors::SCellToInnerPoint.
Obviously, after the computation, you must flip the points of the accepted point set from the background to the volume in order to complete the dilation operation.

Table of Contents