Repository: DGtal-team/DGtal (language: python)

Issue:
Export in PDF 3D: format U3D
A great feature for the 3D visualisation should be to have an export with the format U3D which allows to display interactive 3D view in pdf format. 
See for instance this example:  http://www.loria.fr/~kerautre/lobster.pdf
You just need to read the pdf file with acroread for instance. 
It can be generated simply from latex file and/or in a beamer presentation!

It could be a great feature for the next potential demo/tutorial session...
Perhaps Martial it could be included with the cairo related framework  ?



Discussion / hints:
I agree, it's nice but a bit limited to acroread only...

Do you really want this for milestone 0.4 ?

yes but perhaps extended in the future.
Not necessary for 0.4 but perhaps if Martial has just worked on it, it could easier to do it for 0.4 than later.
But it is perhaps not the priority ;)

I assign Martial to this issue just to keep him in touch..

I just add a very simple example to test pdflatex for creating pdf with 3d image...
http://www.loria.fr/~kerautre/ExampleLatexPDF3D.tgz
From the archive you can check by simply do a:

> pdflatex examplePdf3d.tex
> and 
> acroread examplePdf3d.pdf

That all ;)

I've got a workaround from a DGtal OFF export: meshlab which seems to be able to export to u3d.

Done in the previous pull request...  ;)


Fix patch:
diff --git a/AUTHORS b/AUTHORS
index 5e1fc33220..f1db1527e2 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -28,4 +28,5 @@ Boris Mansencal <boris.mansencal@labri.fr>
 Raphael Lenain <raphael.lenain@univ-poitiers.fr>
 Johanna Delanoy <johanna.delanoy@inria.fr>
 Phuc Ngo <hoai-diem-phuc.ngo@loria.fr>
-J. Miguel Salazar <msalazar@centrogeo.edu.mx>
\ No newline at end of file
+J. Miguel Salazar <msalazar@centrogeo.edu.mx>
+Robin Lamy <robin.lamy1@grenoble-inp.org>
\ No newline at end of file
diff --git a/ChangeLog.md b/ChangeLog.md
index 732395755d..002e0945be 100644
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -40,6 +40,12 @@
     (Bertrand Kerautret, Phuc Ngo and David Coeurjolly
     [#1593](https://github.com/DGtal-team/DGtal/pull/1593))
 
+- *Geometry*
+  - New VoronoiMapComplete class to store the full Voronoi map (with
+    all co-cycling sites (Robin Lamy, David Coeurjolly, Isabelle
+    Sivignon [#1605](https://github.com/DGtal-team/DGtal/pull/1605))
+    
+
 - *Build*
   - New cmake option (DGTAL_RANDOMIZED_TESTING_THRESHOLD) to set the
     (approximated) % of unit-tests to build and run for randomized
diff --git a/src/DGtal/doc/global.bib b/src/DGtal/doc/global.bib
index e37d921bac..7af375b1d4 100644
--- a/src/DGtal/doc/global.bib
+++ b/src/DGtal/doc/global.bib
@@ -598,6 +598,17 @@ @InProceedings{dcoeurjo_metric14
   url           = {http://liris.cnrs.fr/publis/?id=6698}
 }
 
+@article{Couprie2007,
+  title={Discrete bisector function and Euclidean skeleton in 2D and 3D},
+  author={Couprie, Michel and Coeurjolly, David and Zrour, Rita},
+  journal={Image and Vision Computing},
+  volume={25},
+  number={10},
+  pages={1543--1556},
+  year={2007},
+  publisher={Elsevier}
+}
+
 
 @article{dcoeurjo_pami_RDMA,
 	Month = {mar},
diff --git a/src/DGtal/geometry/doc/moduleVolumetric.dox b/src/DGtal/geometry/doc/moduleVolumetric.dox
index ee867621a8..eeef93a71e 100644
--- a/src/DGtal/geometry/doc/moduleVolumetric.dox
+++ b/src/DGtal/geometry/doc/moduleVolumetric.dox
@@ -21,7 +21,7 @@ namespace DGtal {
 
 [TOC]
 
-@author David Coeurjolly, Roland Denis
+@author David Coeurjolly, Roland Denis, Robin Lamy (VoronoiComplete)
 
 This part of the manual describes the DGtal volumetric module. We
 focus here on separable process based volumetric analysis such as
@@ -102,7 +102,9 @@ similarly.
 output of the result is the intersection between the Voronoi diagram
 of exterior points with @f$ \mathbb{Z}^d @f$. Furthermore, along
 Voronoi faces (@e e.i. when two sites are equidistant), only one sites
-is considered when intersection with  @f$ \mathbb{Z}^d @f$.
+is considered when intersection with  @f$ \mathbb{Z}^d @f$. If you
+want the complete Voronoi map (i.e. h @f$ \mathcal{V}\cap\mathbb{Z}^d
+@f$, have a look to @ref vorocomplete.
 
 
 In DGtal, the class VoronoiMap implements such digital Vornoi map
@@ -210,6 +212,50 @@ new Voronoi map:
 @image html voronoimap-vorol8.png "Voronoi map for the l_8 metric."
 @image latex voronoimap-vorol8.png "Voronoi map for the l_8 metric."
 
+@section vorocomplete Complete Discrete Voronoi Map
+
+In this section, we describe the VoronoiMapComplete class that implements
+the complete Voronoi map with all co-cyclic points. More precisely if
+a grid point is equidistant to two sites, the two sites are reported
+in the VoronoiComplete map .
+
+This extension of the separable approach has been described in @cite Couprie2007 
+
+The overall algorithm is roughly the same
+(separable decomposition, each grid point is processed a constant
+amount of time), the only difference is that per grid point, we need
+to iterate over all co-cyclic sites obtained from the partial Voronoi
+map of the preceding dimensions. If the domain has no co-cycling
+sites, the asymptotic computational cost is exactly the same as
+VoronoiMap, if not, there is an extra factor that is a function of the
+number of co-cyclic points.
+
+The interface is exactly the same as the VoronoiMap, except that the
+output image is an image of `std::set<Point>` instead of an image of
+`Point`. As we rely on `std::set` as a per pixel container, and if
+@f$f@f$ denotes the number of co-cyclic sites at a point, the extra
+factor is in @f$ O(f \log f) @f$ (as insertion is expected to be done
+in logarithmic time for `std::set`). 
+
+
+@note Beside the extra computational cost due to the processing of
+co-cyclic sites, the VoronoiMapComplete class as a bigger constant
+factor due to the construction of the container. Indeed, as
+illustrated in testVoronoiMapComplete, creating a VoronoiMap on a
+256x256 image without any site takes 4ms, whereas the same size
+VoronoiMapComplete (still no sites) takes 118ms (on a MacPro). To be
+complete, when adding 5000 random sites, VoronoiMap takes 6.2ms in
+total, and 162ms for the VoronoiMapComplete (maximum number of
+co-cyclic sites = 4). The construction time should be amortized on
+larger domains.
+
+@note Associative containers with O(1) amortized insertion like
+`std::unordered_set` lead to higher construction time without a clear
+improvment on the overall performances.
+
+The VoronoiMapComplete class can also be applied on toroidal domains
+(see @ref toricVol).
+
 @section DTsec Distance Transformation
 
 As discussed earlier, the distance transformation is given by
@@ -301,13 +347,13 @@ ExactPredicateLpPowerSeparableMetric, see @ref moduleMetrics
 
 In some specific applications, toric domains and volumetric analysis
 of shapes in toric domains are crucial. Thanks to the separability
-property of VoronoiMap, PowerMap (and their associated wrappers
+property of VoronoiMap, VoronoiMapComplete, PowerMap (and their associated wrappers
 DistanceTransformation and ReverseDistanceTransformation), one can
 easily consider volumetric transformation in arbitrary dimension on
 toric domains @cite Coeurjo2008. Note that changing the periodicity
 property of the domain has no impact on the computational cost.
 
-More precisely, VoronoiMap and DistanceTransformation classes have
+More precisely, VoronoiMap, VoronoiMapComplete and DistanceTransformation classes have
 constructors allowing to specify periodicity information. In dimension
 @e d, such periodicity information is an array with @e d boolean flags
 where the i-th value is \c true if the i-th dimension of the space is
diff --git a/src/DGtal/geometry/doc/packageGeometry.dox b/src/DGtal/geometry/doc/packageGeometry.dox
index f0b8a196c1..526bd2d339 100644
--- a/src/DGtal/geometry/doc/packageGeometry.dox
+++ b/src/DGtal/geometry/doc/packageGeometry.dox
@@ -74,7 +74,7 @@ of arbitrary dimension, by the means of separable and incremental distance trans
 
 - Volumetric analysis
   - \ref moduleMetrics (David Coeurjolly)
-  - \ref moduleVolumetric (David Coeurjolly, Roland Denis)
+  - \ref moduleVolumetric (David Coeurjolly, Roland Denis, Robin Lamy (VoronoiComplete))
   - \ref moduleFMM (Tristan Roussillon) 
   - \ref moduleDigitalConvexity (Jacques-Olivier Lachaud)
 
diff --git a/src/DGtal/geometry/volumes/distance/VoronoiMap.h b/src/DGtal/geometry/volumes/distance/VoronoiMap.h
index 945316b0a4..d872882b14 100644
--- a/src/DGtal/geometry/volumes/distance/VoronoiMap.h
+++ b/src/DGtal/geometry/volumes/distance/VoronoiMap.h
@@ -73,8 +73,8 @@ namespace DGtal
    *
    * Given a domain and a point predicate, an instance returns, for
    * each point in the domain, the closest point for which the
-   * predicate if false. Following Computational Geometry
-   * terminoliogy, points for which the predicate is false are "sites"
+   * predicate is false. Following Computational Geometry
+   * terminology, points for which the predicate is false are "sites"
    * for the Voronoi map construction. If a point is equi-distant to
    * two sites (e.g. if the digital point belong to a Voronoi cell
    * boundary in the Euclidean space), this Voronoi map construction
diff --git a/src/DGtal/geometry/volumes/distance/VoronoiMap.ih b/src/DGtal/geometry/volumes/distance/VoronoiMap.ih
index 179dc6d537..f648d71c38 100644
--- a/src/DGtal/geometry/volumes/distance/VoronoiMap.ih
+++ b/src/DGtal/geometry/volumes/distance/VoronoiMap.ih
@@ -29,11 +29,6 @@
 
 //////////////////////////////////////////////////////////////////////////////
 #include <cstdlib>
-
-#ifdef VERBOSE
-#include <boost/lexical_cast.hpp>
-#endif
-
 #include "DGtal/kernel/NumberTraits.h"
 
 //////////////////////////////////////////////////////////////////////////////
@@ -76,7 +71,7 @@ void
 DGtal::VoronoiMap<S,P, TSep, TImage>::computeOtherSteps ( const Dimension dim ) const
 {
 #ifdef VERBOSE
-  std::string title = "VoronoiMap dimension " +  boost::lexical_cast<std::string>( dim ) ;
+  std::string title = "VoronoiMap dimension " +  std::to_string( dim ) ;
   trace.beginBlock ( title );
 #endif
 
diff --git a/src/DGtal/geometry/volumes/distance/VoronoiMapComplete.h b/src/DGtal/geometry/volumes/distance/VoronoiMapComplete.h
new file mode 100644
index 0000000000..d192e30c43
--- /dev/null
+++ b/src/DGtal/geometry/volumes/distance/VoronoiMapComplete.h
@@ -0,0 +1,449 @@
+/**
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU Lesser General Public License as
+ *  published by the Free Software Foundation, either version 3 of the
+ *  License, or  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ **/
+
+#pragma once
+
+/**
+ * @file VoronoiMapComplete.h
+ * @brief Linear in time distance transformation
+ * @author David Coeurjolly (\c david.coeurjolly@liris.cnrs.fr )
+ * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
+ * @author Robin Lamy
+ *
+ * @date 2021/09/30
+ *
+ * Header file for module VoronoiMapComplete.cpp
+ *
+ * This file is part of the DGtal library.
+ *
+ * @see testVoronoiMapComplete.cpp
+ */
+#if defined(VoronoiMapComplete_RECURSES)
+#error Recursive header files inclusion detected in VoronoiMapComplete.h
+#else // defined(VoronoiMapComplete_RECURSES)
+/** Prevents recursive inclusion of headers. */
+#define VoronoiMapComplete_RECURSES
+
+#if !defined VoronoiMapComplete_h
+/** Prevents repeated inclusion of headers. */
+#define VoronoiMapComplete_h
+
+//////////////////////////////////////////////////////////////////////////////
+// Inclusions
+#include <iostream>
+#include <vector>
+#include <array>
+#include <set>
+#include "DGtal/base/Common.h"
+#include "DGtal/base/CountedPtr.h"
+#include "DGtal/images/ImageContainerBySTLVector.h"
+#include "DGtal/images/CImage.h"
+#include "DGtal/kernel/CPointPredicate.h"
+#include "DGtal/kernel/domains/HyperRectDomain.h"
+#include "DGtal/geometry/volumes/distance/CSeparableMetric.h"
+#include "DGtal/kernel/domains/HyperRectDomain.h"
+#include "DGtal/base/ConstAlias.h"
+//////////////////////////////////////////////////////////////////////////////
+
+namespace DGtal
+{
+
+  /////////////////////////////////////////////////////////////////////////////
+  // template class VoronoiMapComplete
+  /**
+   * Description of template class 'VoronoiMapComplete' <p>
+   * \brief Aim: Implementation of the linear in time Voronoi map
+   * construction.
+
+   * The algorithm uses a seperable process to construct Voronoi maps
+   * which has been described in @cite Maurer2003PAMI @cite dcoeurjo_these.
+   * Along periodic dimensions, the algorithm is adapted following @cite Coeurjo2008,
+   * and uses the co-cyclic site approach from @cite Couprie2007.
+   *
+   * Given a domain and a point predicate, an instance returns, for
+   * each point in the domain, the collection of closest points for which the
+   * predicate is false. Following Computational Geometry
+   * terminology, points for which the predicate is false are "sites"
+   * for the Voronoi map construction.
+   *
+   * If a point is equi-distant to several sites, all sites will be attached to
+   * that point (contrary to the class VoronoiMap that will only keep one of them).
+   * This implies a computational overhead:
+   *    - if @f$ f @f$ is the max number of co-cyclic points per grid point, we
+   * have an extra @f$ O(f log f) @f$ factor to the computational cost (cf below), plus the cost
+   * of insertion/removal in the container (O(log f) with std::set).
+   *    - as we have to use a specific container to store the per pixel equi-distant points
+   * we obviously have a memory overhead when there are equi-distant points, but a slight computational
+   * one too to access sites from the container (std::set<Point> per grid point).
+   *
+   * See @ref moduleVolumetric documnetation page.
+   *
+   * By default, the domain is considered non-periodic but per-dimension
+   * periodicity can be specified in the constructor.
+   * When the domain has periodic dimensions, the closest point
+   * coordinates \c B to a given point \c A may not be between the lower
+   * and upper bounds of the domain, in such a way that the non-periodic
+   * distance between \c A and \c B is equal to their distance considering
+   * the periodicity.
+   *
+   * The metric is specified by a model of concepts::CSeparableMetric (for
+   * instance, any instance of ExactPredicateLpSeparableMetric or
+   * InexactPredicateLpSeparableMetric).  If the separable metric has
+   * a complexity of O(h) for its "hiddenBy" predicate, the overall
+   * Voronoi construction algorithm is in @f$ O(f.h.d.n^d)@f$ for @f$
+   * n^d@f$ domains (see class constructor). For Euclidean the @f$
+   * l_2@f$ metric, the overall computation is in @f$ O(f.d.n^d)@f$,
+   * which is optimal.
+   *
+   * If DGtal has been built with OpenMP support (WITH_OPENMP flag set
+   * to "true"), the computation is done in parallel (multithreaded)
+   * in an optimal way: on @a p processors, expected runtime is in
+   * @f$ O(f.h.d.n^d / p)@f$.
+   *
+   * This class is a model of concepts::CConstImage.
+   *
+   * @tparam TSpace type of Digital Space (model of concepts::CSpace).
+   * @tparam TPointPredicate point predicate returning true for points
+   * from which we compute the distance (model of concepts::CPointPredicate)
+   * @tparam TSeparableMetric a model of concepts::CSeparableMetric
+   * @tparam TImageContainer any model of concepts::CImage to store the
+   * VoronoiMap (default: ImageContainerBySTLVector). The space of the
+   * image container and the TSpace should match. Furthermore the
+   * container value type must be a random access container on TSpace::Vector.
+   * Lastly, the domain of the container must be HyperRectDomain.
+   */
+  template < typename TSpace,
+             typename TPointPredicate,
+             typename TSeparableMetric,
+             typename TImageContainer =
+             ImageContainerBySTLVector<HyperRectDomain<TSpace>,
+                                       std::set<typename TSpace::Vector> >
+             >
+  class VoronoiMapComplete
+  {
+
+  public:
+    BOOST_CONCEPT_ASSERT(( concepts::CSpace< TSpace > ));
+    BOOST_CONCEPT_ASSERT(( concepts::CPointPredicate<TPointPredicate> ));
+    BOOST_CONCEPT_ASSERT(( concepts::CSeparableMetric<TSeparableMetric> ));
+    BOOST_CONCEPT_ASSERT(( concepts::CImage< TImageContainer > ));
+
+    ///Both Space points and PointPredicate points must be the same.
+    BOOST_STATIC_ASSERT ((boost::is_same< typename TSpace::Point,
+                          typename TPointPredicate::Point >::value ));
+
+    //ImageContainer::Domain::Space must match with TSpace
+    BOOST_STATIC_ASSERT ((boost::is_same< TSpace,
+                          typename TImageContainer::Domain::Space >::value ));
+    
+    //ImageContainer value type must be  std::vector<TSpace::Vector>
+    BOOST_STATIC_ASSERT ((boost::is_same< std::set<typename TSpace::Vector>,
+                          typename TImageContainer::Value >::value ));
+
+    //ImageContainer domain type must be  HyperRectangular
+    BOOST_STATIC_ASSERT ((boost::is_same< HyperRectDomain<TSpace>,
+                          typename TImageContainer::Domain >::value ));
+
+    ///Copy of the space type.
+    typedef TSpace Space;
+
+    ///Copy of the point predicate type.
+    typedef TPointPredicate PointPredicate;
+
+    ///Definition of the underlying domain type.
+    typedef typename TImageContainer::Domain Domain;
+
+    ///Definition of the separable metric type
+    typedef TSeparableMetric SeparableMetric;
+
+    ///Large integer type for SeparableMetricHelper construction.
+    typedef DGtal::int64_t IntegerLong;
+
+    typedef typename Space::Vector Vector;
+    typedef typename Space::Point Point;
+    typedef typename Space::Dimension Dimension;
+    typedef typename Space::Size Size;
+    typedef typename Space::Point::Coordinate Abscissa;
+
+    ///Type of resulting image
+    typedef TImageContainer OutputImage;
+
+    ///Definition of the image value type.
+    typedef typename OutputImage::Value Value;
+
+    ///Definition of the image value type.
+    typedef typename OutputImage::ConstRange  ConstRange;
+
+    ///Self type
+    typedef VoronoiMapComplete< TSpace, TPointPredicate,
+                        TSeparableMetric,TImageContainer > Self;
+
+
+    /// Periodicity specification type.
+    typedef std::array< bool, Space::dimension > PeriodicitySpec;
+
+    /**
+     * Constructor in the non-periodic case.
+     *
+     * This constructor computes the Voronoi Map of a set of point
+     * sites using a SeparableMetric metric, on a non-periodic domain.
+     *
+     * The method associates to each point satisfying the foreground predicate,
+     * the closest site for which the predicate is false. This algorithm is
+     * @f$ O(h.d.|domain size|)@f$ if the separable metric "hiddenBy"
+     * predicate is in @f$ O(h)@f$.
+     *
+     * @param aDomain a pointer to the (hyper-rectangular) domain on
+     * which the computation is performed.
+     *
+     * @param predicate a pointer to the point predicate to define the
+     * Voronoi sites (false points).
+     *
+     * @param aMetric a pointer to the separable metric instance.
+     */
+    VoronoiMapComplete(ConstAlias<Domain> aDomain,
+                       ConstAlias<PointPredicate> predicate,
+                       ConstAlias<SeparableMetric> aMetric);
+
+    /**
+     * Constructor with periodicity specification.
+     *
+     * This constructor computes the Voronoi Map of a set of point
+     * sites using a SeparableMetric metric, on a domain with specified
+     * periodicity.
+     *
+     * The method associates to each point satisfying the foreground predicate,
+     * the closest site for which the predicate is false. This algorithm is
+     * @f$ O(h.d.|domain size|)@f$ if the separable metric "hiddenBy"
+     * predicate is in @f$ O(h)@f$.
+     *
+     * @param aDomain a pointer to the (hyper-rectangular) domain on
+     * which the computation is performed.
+     *
+     * @param predicate a pointer to the point predicate to define the
+     * Voronoi sites (false points).
+     *
+     * @param aMetric a pointer to the separable metric instance.
+     *
+     * @param aPeriodicitySpec an array of size equal to the space dimension
+     *        where the i-th value is \c true if the i-th dimension of the
+     *        space is periodic, \c false otherwise.
+     */
+    VoronoiMapComplete(ConstAlias<Domain> aDomain,
+                       ConstAlias<PointPredicate> predicate,
+                       ConstAlias<SeparableMetric> aMetric,
+                       PeriodicitySpec const & aPeriodicitySpec);
+    /**
+     * Default destructor
+     */
+    ~VoronoiMapComplete() = default;
+
+    /**
+     * Disabling default constructor.
+     */
+     VoronoiMapComplete() = delete;
+
+  public:
+    // ------------------- ConstImage model ------------------------
+
+    /**
+     * Assignment operator from another Voronoi map.
+     *
+     *  @param aOtherVoronoiMap another instance of Self
+     *  @return a reference to Self
+     */
+    Self &  operator=(const Self &aOtherVoronoiMap ) = default;
+
+    /**
+     * Returns a reference (const) to the Voronoi map domain.
+     * @return a domain
+     */
+    const Domain &  domain() const
+    {
+      return *myDomainPtr;
+    }
+
+    /**
+     * Returns a const range on the Voronoi map values.
+     *  @return a const range
+     */
+    ConstRange constRange() const
+    {
+      return myImagePtr->constRange();
+    }
+
+    /**
+     * Access to a Voronoi value (a.k.a. vector to the closest site)
+     * at a point.
+     *
+     * @param aPoint the point to probe.
+     */
+    Value operator()(const Point &aPoint) const
+    {
+      return myImagePtr->operator()(aPoint);
+    }
+
+    /**
+     * @return Returns an alias to the underlying metric.
+     */
+    const SeparableMetric* metric() const
+    {
+      return myMetricPtr;
+    }
+
+    /** Periodicity specification.
+     *
+     * @returns the periodicity specification array.
+     */
+    PeriodicitySpec const & getPeriodicitySpec() const
+      {
+        return myPeriodicitySpec;
+      }
+
+    /** Periodicity specification along one dimensions.
+     *
+     * @param [in] n the dimension index.
+     * @return \c true if the n-th dimension is periodic, \c false otherwise.
+     */
+    bool isPeriodic( const Dimension n ) const
+      {
+        return myPeriodicitySpec[ n ];
+      }
+
+    /**
+     * Project point coordinates into the domain, taking into account
+     * the periodicity.
+     *
+     * @pre The given point must come from operator()(const Point &) const (for performance reasons).
+     *
+     * @param aPoint the point to project
+     * @return the coordinates projected into the domain bounds accordingly
+     *         to the periodicity specification.
+     */
+    Point projectPoint( Point aPoint ) const;
+
+    /**
+     * Self Display method.
+     *
+     * @param out output stream
+     */
+    void selfDisplay ( std::ostream & out ) const;
+
+    // ------------------- Private functions ------------------------
+  private:
+
+    /**
+     * Compute the Voronoi Map of a set of point sites using a
+     * SeparableMetric metric.  The method associates to each point
+     * satisfying the foreground predicate, the closest site for which
+     * the predicate is false. This algorithm is O(h.d.|domain size|).
+     */
+    void compute ( ) ;
+
+    /**
+     *  Compute the other steps of the separable Voronoi map.
+     *
+     * @param [in] dim the dimension to process
+     */
+    void computeOtherSteps(const Dimension dim) const;
+    /**
+     * Given  a voronoi map valid at dimension @a dim-1, this method
+     * updates the map to make it consistent at dimension @a dim along
+     * the 1D span starting at @a row along the dimension @a
+     * dim.
+     *
+     * @param [in] row starting point of the 1D process.
+     * @param [in] dim dimension of the update.
+     */
+    void computeOtherStep1D (const Point &row,
+                             const Dimension dim) const;
+
+    /**
+     * Project a coordinate into the domain, taking into account
+     * the periodicity.
+     *
+     * @pre The given coordinate must come from operator()(const Point &) const (for performance reasons).
+     *
+     * @param aCoordinate the coordinate.
+     * @param aDim  dimension of the coordinate.
+     * @return the coordinates projected into the domain bounds accordingly
+     *         to the periodicity specification.
+     */
+    typename Point::Coordinate projectCoordinate( typename Point::Coordinate aCoordinate, const Dimension aDim ) const;
+
+    // ------------------- Private members ------------------------
+  private:
+
+    ///Pointer to the computation domain
+    const Domain * myDomainPtr;
+
+    ///Pointer to the point predicate
+    const PointPredicate * myPointPredicatePtr;
+
+    ///Copy of the image lower bound
+    Point myLowerBoundCopy;
+
+    ///Copy of the image lower bound
+    Point myUpperBoundCopy;
+
+    ///Value to act as a +infinity value
+    Point myInfinity;
+
+    /// Index of the periodic dimensions
+    std::vector< Dimension > myPeriodicityIndex; // Could be boost::static_vector but it needs Boost >= 1.54.
+
+    /// Domain extent.
+    Point myDomainExtent;
+
+  protected:
+
+    ///Pointer to the separable metric instance
+    const SeparableMetric * myMetricPtr;
+
+    ///Voronoi map image
+    CountedPtr<OutputImage> myImagePtr;
+
+    /// Periodicity along each dimension.
+    PeriodicitySpec myPeriodicitySpec;
+
+  }; // end of class VoronoiMapComplete
+
+  /**
+   * Overloads 'operator<<' for displaying objects of class 'ExactPredicateLpSeparableMetric'.
+   * @param out the output stream where the object is written.
+   * @param object the object of class 'ExactPredicateLpSeparableMetric' to write.
+   * @return the output stream after the writing.
+   */
+  template <typename S, typename P,
+            typename Sep, typename TI>
+  std::ostream&
+  operator<< ( std::ostream & out, const VoronoiMapComplete<S,P,Sep,TI> & object );
+
+
+} // namespace DGtal
+
+
+///////////////////////////////////////////////////////////////////////////////
+// Includes inline functions.
+#include "DGtal/geometry/volumes/distance/VoronoiMapComplete.ih"
+
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#endif // !defined VoronoiMap_h
+
+#undef VoronoiMapComplete_RECURSES
+#endif // else defined(VoronoiMapComplete_RECURSES)
diff --git a/src/DGtal/geometry/volumes/distance/VoronoiMapComplete.ih b/src/DGtal/geometry/volumes/distance/VoronoiMapComplete.ih
new file mode 100644
index 0000000000..1a35b63288
--- /dev/null
+++ b/src/DGtal/geometry/volumes/distance/VoronoiMapComplete.ih
@@ -0,0 +1,409 @@
+/**
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU Lesser General Public License as
+ *  published by the Free Software Foundation, either version 3 of the
+ *  License, or  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ **/
+
+/**
+ * @file
+ *
+ * @author Robin Lamy
+ * @author David Coeurjolly (\c david.coeurjolly@liris.cnrs.fr )
+ * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS,
+ UMR 5205), CNRS, France
+ *
+ * initial 2012/08/14
+ * @date 2021/07/20
+
+ * Implementation of inline methods defined in VoronoiMap.h
+ *
+ * This file is part of the DGtal library.
+ */
+
+//////////////////////////////////////////////////////////////////////////////
+#include <cstdlib>
+#include "DGtal/kernel/NumberTraits.h"
+
+///////////////////////////////////////////////////////////////////////////////
+// IMPLEMENTATION of inline methods.
+///////////////////////////////////////////////////////////////////////////////
+// ----------------------- Standard services ------------------------------
+
+template <typename S, typename P, typename TSep, typename TImage>
+inline void DGtal::VoronoiMapComplete<S, P, TSep, TImage>::compute()
+{
+  // We copy the image extent
+  myLowerBoundCopy = myDomainPtr->lowerBound();
+  myUpperBoundCopy = myDomainPtr->upperBound();
+
+  // Point outside the domain
+  for ( auto & coord : myInfinity )
+    coord = DGtal::NumberTraits<typename Point::Coordinate>::max();
+
+  // Init
+  for ( auto const & pt : *myDomainPtr )
+    if ( ( *myPointPredicatePtr )( pt ) )
+    {
+      Value myVectorInfinity;
+      myVectorInfinity.insert( myInfinity );
+      myImagePtr->setValue( pt, myVectorInfinity );
+    }
+    else
+    {
+      Value external_point;
+      external_point.insert( pt );
+      myImagePtr->setValue( pt, external_point );
+    }
+
+  // We process the remaining dimensions
+  for ( Dimension dim = 0; dim < S::dimension; dim++ )
+  {
+    computeOtherSteps( dim );
+  }
+}
+
+template <typename S, typename P, typename TSep, typename TImage>
+inline void DGtal::VoronoiMapComplete<S, P, TSep, TImage>::computeOtherSteps(
+const Dimension dim ) const
+{
+#ifdef VERBOSE
+  std::string title =  "VoronoiMapComplete dimension " + std::to_string( dim );
+  trace.beginBlock( title );
+#endif
+
+  // We setup the subdomain iterator
+  // the iterator will scan dimension using the order:
+  //  {n-1, n-2, ... 1} (we skip the '0' dimension).
+  std::vector<Dimension> subdomain;
+  subdomain.reserve( S::dimension - 1 );
+  for ( int k = 0; k < (int)S::dimension; k++ )
+    if ( static_cast<Dimension>( ( (int)S::dimension - 1 - k ) ) != dim )
+      subdomain.push_back( (int)S::dimension - 1 - k );
+
+  Domain localDomain( myLowerBoundCopy, myUpperBoundCopy );
+
+#ifdef WITH_OPENMP
+  // Parallel loop
+  std::vector<Point> subRangePoints;
+  // Starting point precomputation
+  for ( auto const & pt : localDomain.subRange( subdomain ) )
+    subRangePoints.push_back( pt );
+
+    // We run the 1D problems in //
+#pragma omp parallel for schedule( dynamic )
+  for ( size_t i = 0; i < subRangePoints.size(); ++i )
+    computeOtherStep1D( subRangePoints[ i ], dim );
+
+#else
+  // We solve the 1D problems sequentially
+  for ( auto const & pt : localDomain.subRange( subdomain ) )
+    computeOtherStep1D( pt, dim );
+#endif
+
+#ifdef VERBOSE
+  trace.endBlock();
+#endif
+}
+
+// //////////////////////////////////////////////////////////////////////:
+// ////////////////////////// Other Phases
+template <typename S, typename P, typename TSep, typename TImage>
+void DGtal::VoronoiMapComplete<S, P, TSep, TImage>::computeOtherStep1D(
+const Point & startingPoint, const Dimension dim ) const
+{
+  ASSERT( dim < S::dimension );
+
+  // Default starting and ending point for a cycle
+  Point startPoint  = startingPoint;
+  Point endPoint    = startingPoint;
+  startPoint[ dim ] = myLowerBoundCopy[ dim ];
+  endPoint[ dim ]   = myUpperBoundCopy[ dim ];
+
+  Value myVectorInfinity;
+  myVectorInfinity.emplace( myInfinity );
+
+  // Extent along current dimension.
+  const auto extent = myUpperBoundCopy[ dim ] - myLowerBoundCopy[ dim ] + 1;
+
+  // Site storage.
+  std::vector<Point> Sites;
+
+  // Reserve sites storage.
+  // +1 along periodic dimension in order to store two times the site that is on
+  // break index.
+  Sites.reserve( extent + ( isPeriodic( dim ) ? 1 : 0 ) );
+
+  // Pruning the list of sites and defining cycle bounds.
+  // In the periodic case, the cycle bounds depend on the so-called break index
+  // that defines the start point.
+  if ( dim == 0 )
+  {
+    // For dim = 0, no sites are hidden.
+    for ( auto point = startPoint; point[ dim ] <= myUpperBoundCopy[ dim ];
+          ++point[ dim ] )
+    {
+      const Value psite = myImagePtr->operator()( point );
+      if ( psite != myVectorInfinity )
+      {
+        for ( Point voronoiPoint : psite )
+        {
+          Sites.push_back( voronoiPoint );
+        }
+      }
+    }
+
+    // If no sites are found, then there is nothing to do.
+    if ( Sites.size() == 0 )
+      return;
+
+    // In the periodic case and along the first dimension, the break index
+    // is at the first site found.
+    if ( isPeriodic( dim ) )
+    {
+      startPoint[ dim ] = Sites[ 0 ][ dim ];
+      endPoint[ dim ]   = startPoint[ dim ] + extent - 1;
+
+      // The first site is also the last site (with appropriate shift).
+      Sites.push_back( Sites[ 0 ] + Point::base( dim, extent ) );
+    }
+  }
+  else
+  {
+    // In the periodic case, the cycle depends on break index
+    if ( isPeriodic( dim ) )
+    {
+      // Along other than the first dimension, the break index is at the lowest
+      // site found.
+      auto minRawDist =
+      DGtal::NumberTraits<typename SeparableMetric::RawValue>::max();
+
+      for ( auto point = startPoint; point[ dim ] <= myUpperBoundCopy[ dim ];
+            ++point[ dim ] )
+      {
+        const Value psite = myImagePtr->operator()( point );
+
+        if ( psite != myVectorInfinity )
+        {
+          const auto rawDist = myMetricPtr->rawDistance(
+          point, *( psite.begin() ) ); // All distances to Voronoi sites are
+                                       // equal, we take the first one
+          if ( rawDist < minRawDist )
+          {
+            minRawDist        = rawDist;
+            startPoint[ dim ] = point[ dim ];
+          }
+        }
+      }
+
+      // If no sites are found, then there is nothing to do.
+      if ( minRawDist ==
+           DGtal::NumberTraits<typename SeparableMetric::RawValue>::max() )
+        return;
+
+      endPoint[ dim ] = startPoint[ dim ] + extent - 1;
+    }
+
+    // Pruning the list of sites for both periodic and non-periodic cases.
+    for ( auto point = startPoint; point[ dim ] <= myUpperBoundCopy[ dim ];
+          ++point[ dim ] )
+    {
+      const Value psite = myImagePtr->operator()( point );
+
+      if ( psite != myVectorInfinity )
+      {
+        for ( Point site : psite )
+        {
+          while ( ( Sites.size() >= 2 ) &&
+                  ( myMetricPtr->hiddenBy( Sites[ Sites.size() - 2 ],
+                                           Sites[ Sites.size() - 1 ], site,
+                                           startingPoint, endPoint, dim ) ) )
+          {
+
+            Sites.pop_back();
+          }
+
+          Sites.push_back( site );
+        }
+      }
+    }
+
+    // Pruning the remaining list of sites in the periodic case.
+    if ( isPeriodic( dim ) )
+    {
+      auto point   = startPoint;
+      point[ dim ] = myLowerBoundCopy[ dim ];
+      for ( ; point[ dim ] <= endPoint[ dim ] - extent + 1;
+            ++point[ dim ] ) // +1 in order to add the break-index site at the
+                             // cycle's end.
+      {
+        const Value psite = myImagePtr->operator()( point );
+
+        if ( psite != myVectorInfinity )
+        {
+          for ( Point site : psite )
+          {
+            // Site coordinates must be between startPoint and endPoint.
+            site[ dim ] += extent;
+
+            while ( ( Sites.size() >= 2 ) &&
+                    ( myMetricPtr->hiddenBy( Sites[ Sites.size() - 2 ],
+                                             Sites[ Sites.size() - 1 ], site,
+                                             startingPoint, endPoint, dim ) ) )
+              Sites.pop_back();
+
+            Sites.push_back( site );
+          }
+        }
+      }
+    }
+  }
+
+  // No sites found
+  if ( Sites.size() == 0 )
+    return;
+
+  // Rewriting for both periodic and non-periodic cases.
+  size_t siteId = 0;
+  auto point    = startPoint;
+  IntegerLong minimal_distance;
+
+  for ( ; point[ dim ] <= myUpperBoundCopy[ dim ]; ++point[ dim ] )
+  {
+    Value voronoiSites = myImagePtr->operator()( point );
+    minimal_distance =
+    myMetricPtr->rawDistance( *( voronoiSites.begin() ), point );
+
+    for ( Point site : Sites )
+    {
+      if ( myMetricPtr->rawDistance( site, point ) < minimal_distance )
+      {
+        voronoiSites.clear();
+        voronoiSites.emplace( site );
+        minimal_distance = myMetricPtr->rawDistance( site, point );
+      }
+      else if ( ( myMetricPtr->rawDistance( site, point ) ==
+                  minimal_distance ) &&
+                ( std::find( voronoiSites.begin(), voronoiSites.end(), site ) ==
+                  voronoiSites.end() ) )
+      {
+        voronoiSites.emplace( site );
+      }
+    }
+    myImagePtr->setValue( point, voronoiSites );
+  }
+
+  // Continuing rewriting in the periodic case.
+  if ( isPeriodic( dim ) )
+  {
+    for ( ; point[ dim ] <= endPoint[ dim ]; ++point[ dim ] )
+    {
+      Value voronoiSites;
+      siteId = 0;
+
+      while ( siteId < Sites.size() - 1 )
+      { // Idem ici en prenant en compte la périodicité
+        if ( myMetricPtr->closest( point, Sites[ siteId ],
+                                   Sites[ siteId + 1 ] ) ==
+             DGtal::ClosestFIRST )
+        {
+          voronoiSites.insert( Sites[ siteId ] - Point::base( dim, extent ) );
+        }
+        siteId++;
+      }
+      myImagePtr->setValue( point - Point::base( dim, extent ), voronoiSites );
+    }
+  }
+}
+
+/**
+ * Constructor.
+ */
+template <typename S, typename P, typename TSep, typename TImage>
+inline DGtal::VoronoiMapComplete<S, P, TSep, TImage>::VoronoiMapComplete(ConstAlias<Domain> aDomain,
+                                                                         ConstAlias<PointPredicate> aPredicate,
+                                                                         ConstAlias<SeparableMetric> aMetric )
+  : myDomainPtr( &aDomain )
+  , myPointPredicatePtr( &aPredicate )
+  , myDomainExtent( aDomain->upperBound() - aDomain->lowerBound() +
+                    Point::diagonal( 1 ) )
+  , myMetricPtr( &aMetric )
+{
+  myPeriodicitySpec.fill( false );
+  myImagePtr = CountedPtr<OutputImage>( new OutputImage( aDomain ) );
+  compute();
+}
+
+template <typename S, typename P, typename TSep, typename TImage>
+inline DGtal::VoronoiMapComplete<S, P, TSep, TImage>::VoronoiMapComplete(ConstAlias<Domain> aDomain,
+                                                                         ConstAlias<PointPredicate> aPredicate,
+                                                                         ConstAlias<SeparableMetric> aMetric,
+                                                                         PeriodicitySpec const & aPeriodicitySpec )
+  : myDomainPtr( &aDomain )
+  , myPointPredicatePtr( &aPredicate )
+  , myDomainExtent( aDomain->upperBound() - aDomain->lowerBound() +
+                    Point::diagonal( 1 ) )
+  , myMetricPtr( &aMetric )
+  , myPeriodicitySpec( aPeriodicitySpec )
+{
+  // Finding periodic dimension index.
+  for ( Dimension i = 0; i < Space::dimension; ++i )
+    if ( isPeriodic( i ) )
+      myPeriodicityIndex.push_back( i );
+
+  myImagePtr = CountedPtr<OutputImage>( new OutputImage( aDomain ) );
+  compute();
+}
+
+template <typename S, typename P, typename TSep, typename TImage>
+inline typename DGtal::VoronoiMapComplete<S, P, TSep, TImage>::Point
+DGtal::VoronoiMapComplete<S, P, TSep, TImage>::projectPoint(
+Point aPoint ) const
+{
+  for ( auto const & dim : myPeriodicityIndex )
+    aPoint[ dim ] = projectCoordinate( aPoint[ dim ], dim );
+
+  return aPoint;
+}
+
+template <typename S, typename P, typename TSep, typename TImage>
+inline typename DGtal::VoronoiMapComplete<S, P, TSep, TImage>::Point::Coordinate
+DGtal::VoronoiMapComplete<S, P, TSep, TImage>::projectCoordinate(typename Point::Coordinate aCoordinate,
+                                                                 const Dimension aDim ) const
+{
+  ASSERT( aCoordinate - myDomainPtr->lowerBound()[ aDim ] +
+          myDomainExtent[ aDim ] >=
+          0 );
+  return ( aCoordinate - myDomainPtr->lowerBound()[ aDim ] +
+           myDomainExtent[ aDim ] ) %
+         myDomainExtent[ aDim ] +
+         myDomainPtr->lowerBound()[ aDim ];
+}
+
+template <typename S, typename P, typename TSep, typename TImage>
+inline void DGtal::VoronoiMapComplete<S, P, TSep, TImage>::selfDisplay(
+std::ostream & out ) const
+{
+  out << "[VoronoiMapComplete] separable metric=" << *myMetricPtr;
+}
+
+// // //
+// ///////////////////////////////////////////////////////////////////////////////
+
+template <typename S, typename P, typename TSep, typename TImage>
+inline std::ostream &
+DGtal::operator<<( std::ostream & out,
+                   const VoronoiMapComplete<S, P, TSep, TImage> & object )
+{
+  object.selfDisplay( out );
+  return out;
+}


Test patch:
diff --git a/tests/geometry/volumes/distance/CMakeLists.txt b/tests/geometry/volumes/distance/CMakeLists.txt
index 5a60c12417..eacd2116e1 100644
--- a/tests/geometry/volumes/distance/CMakeLists.txt
+++ b/tests/geometry/volumes/distance/CMakeLists.txt
@@ -14,6 +14,7 @@ set(DGTAL_TESTS_SRC
   testChamferVoro
   testDigitalMetricAdapter
   testLpMetric
+  testVoronoiMapComplete
   )
 
 
diff --git a/tests/geometry/volumes/distance/testVoronoiMapComplete.cpp b/tests/geometry/volumes/distance/testVoronoiMapComplete.cpp
new file mode 100644
index 0000000000..5f775f3baa
--- /dev/null
+++ b/tests/geometry/volumes/distance/testVoronoiMapComplete.cpp
@@ -0,0 +1,240 @@
+/**
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU Lesser General Public License as
+ *  published by the Free Software Foundation, either version 3 of the
+ *  License, or  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ **/
+
+/**
+ * @file
+ * @ingroup Tests
+ * @author Robin Lamy
+ *
+ * @date 2021/07/20
+ *
+ * Functions for testing class VoronoiMapComplete.
+ *
+ * This file is part of the DGtal library.
+ */
+
+///////////////////////////////////////////////////////////////////////////////
+#include <limits>
+#include <unordered_set>
+#include "DGtal/base/Common.h"
+#include "DGtal/helpers/StdDefs.h"
+#include "DGtalCatch.h"
+
+#include "DGtal/geometry/volumes/distance/VoronoiMapComplete.h"
+#include "DGtal/geometry/volumes/distance/VoronoiMap.h"
+///////////////////////////////////////////////////////////////////////////////
+using namespace DGtal;
+
+typedef Z2i::DigitalSet DigitalSet;
+typedef Z2i::Point Point;
+typedef Z2i::Domain Domain;
+
+typedef VoronoiMapComplete<Z2i::Space, DigitalSet, Z2i::L2Metric> CompleteVMap;
+typedef VoronoiMap<Z2i::Space, DigitalSet, Z2i::L2Metric> VMap;
+typedef ImageContainerBySTLVector<HyperRectDomain<Z2i::Space>,
+                                  std::unordered_set<Z2i::Point>>
+TImageContainer;
+
+TImageContainer * brut_force_voronoi_map_complete( DigitalSet & set )
+{
+  TImageContainer * voronoi_map = new TImageContainer( set.domain() );
+  std::vector<Point> Sites;
+  // Getting all voronoi sites in one vector
+  for ( Point point : set.domain() )
+  {
+    if ( !set( point ) )
+    {
+      Sites.push_back( point );
+    }
+  }
+  std::cout << std::endl;
+
+  for ( Point point : set.domain() )
+  {
+
+    std::unordered_set<Point> voronoi_points;
+    double voronoi_distance =
+    std::numeric_limits<int>::max(); // maximum int value
+    for ( Point site : Sites )
+    {
+      if ( ( point - site ).norm() < voronoi_distance )
+      {
+        voronoi_points.clear();
+        voronoi_points.insert( site );
+        voronoi_distance = ( point - site ).norm();
+      }
+      else if ( ( point - site ).norm() == voronoi_distance )
+      {
+        voronoi_points.insert( site );
+      }
+    }
+    voronoi_map->setValue( point, voronoi_points );
+  }
+  return voronoi_map;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Functions for testing class VoronoiMapComplete
+///////////////////////////////////////////////////////////////////////////////
+TEST_CASE( "Testing VoronoiMapComplete 2D" )
+{
+  /**
+   * Random set generation
+   */
+  Point lowerBound( 0, 0 ), upperBound( 31, 31 );
+  Domain domain( lowerBound, upperBound );
+
+  DigitalSet set( domain );
+  int point_setup_index = 0;
+  int x;
+  int y;
+
+  while ( point_setup_index < 400 )
+  {
+    x = rand() % 32;
+    y = rand() % 32;
+    set.insert( Point( x, y ) );
+    point_setup_index++;
+  }
+
+  /**
+   * Voronoi Map Complete from class generation
+   */
+  CompleteVMap vmap( set.domain(), set, Z2i::l2Metric );
+
+  /**
+   * Voronoi Map Complete from brut force generation
+   */
+  TImageContainer * brut_force_vmap = brut_force_voronoi_map_complete( set );
+
+  SECTION(
+  "Testing that VoronoiMapComplete class gives the right voronoi sites sets" )
+  {
+    for ( Point point : set )
+    {
+      std::unordered_set<Point> brut_force_set =
+      brut_force_vmap->operator()( point );
+      CompleteVMap::Value class_set = vmap( point );
+
+      // Cheking that all voronoi sites from the brut force set are in the
+      // algorithm set
+      for ( Point voronoi_point : brut_force_set )
+      {
+        REQUIRE( std::find( class_set.begin(), class_set.end(),
+                            voronoi_point ) != class_set.end() );
+      }
+
+      // Cheking that all voronoi sites from the algorithm set are in the brut
+      // force set
+      for ( Point voronoi_point : class_set )
+      {
+        REQUIRE( std::find( brut_force_set.begin(), brut_force_set.end(),
+                            voronoi_point ) != brut_force_set.end() );
+      }
+    }
+  }
+
+  SECTION(
+  "Testing Complete Voronoi Map from Discrete Bisector Function paper" )
+  {
+    Point lowerBound( 0, 0 ), upperBound( 6, 7 );
+    Domain domain( lowerBound, upperBound );
+    DigitalSet set( domain );
+    for ( Point point : set.domain() )
+    {
+      if ( point != Point( 1, 0 ) && point != Point( 5, 0 ) &&
+           point != Point( 2, 2 ) && point != Point( 4, 4 ) &&
+           point != Point( 0, 6 ) && point != Point( 6, 6 ) &&
+           point != Point( 3, 7 ) )
+        set.insert( point );
+    }
+
+    TImageContainer * brutForceVoronoiMap =
+    brut_force_voronoi_map_complete( set );
+    CompleteVMap vmap( set.domain(), set, Z2i::l2Metric );
+
+    for ( Point point : set )
+    {
+      std::unordered_set<Point> brut_force_set =
+      brutForceVoronoiMap->operator()( point );
+      CompleteVMap::Value class_set = vmap( point );
+
+      // Cheking that all voronoi sites from the brut force set are in the
+      // algorithm set
+      for ( Point voronoi_point : brut_force_set )
+        REQUIRE( std::find( class_set.begin(), class_set.end(),
+                            voronoi_point ) != class_set.end() );
+
+      // Cheking that all voronoi sites from the algorithm set are in the brut
+      // force set
+      for ( Point voronoi_point : class_set )
+        REQUIRE( std::find( brut_force_set.begin(), brut_force_set.end(),
+                            voronoi_point ) != brut_force_set.end() );
+    }
+  }
+}
+
+TEST_CASE( "No sites" )
+{
+  /**
+   * Random set generation
+   */
+  Point lowerBound( 0, 0 ), upperBound( 255, 255 );
+  Domain domain( lowerBound, upperBound );
+  DigitalSet set( domain );
+
+  trace.beginBlock( "Complete Map (no site)" );
+  CompleteVMap vmap( set.domain(), set, Z2i::l2Metric );
+  trace.endBlock();
+  trace.beginBlock( "Partial Map (no site)" );
+  VMap partialmap( set.domain(), set, Z2i::l2Metric );
+  trace.endBlock();
+}
+
+TEST_CASE( "Testing Timings" )
+{
+  /**
+   * Random set generation
+   */
+  Point lowerBound( 0, 0 ), upperBound( 255, 255 );
+  Domain domain( lowerBound, upperBound );
+
+  DigitalSet set( domain );
+  int point_setup_index = 0;
+  int x;
+  int y;
+
+  while ( point_setup_index < 5000 )
+  {
+    x = rand() % 256;
+    y = rand() % 256;
+    set.insert( Point( x, y ) );
+    point_setup_index++;
+  }
+  std::cout << std::endl;
+  trace.beginBlock( "Complete Map" );
+  CompleteVMap vmap( set.domain(), set, Z2i::l2Metric );
+  trace.endBlock();
+
+  size_t maxsize = 0;
+  for ( auto & v : vmap.constRange() )
+    maxsize = std::max( maxsize, v.size() );
+  trace.info() << "Max number of co-cyclic points = " << maxsize << std::endl;
+
+  trace.beginBlock( "Partial Map" );
+  VMap partialmap( set.domain(), set, Z2i::l2Metric );
+  trace.endBlock();
+}