
/* 
 * Useful to avoid writing DGtal:: in front of every class.
 */
namespace DGtal {

/**
   
@page moduleCubicalComplex  Cubical Complex

@writers Jacques-Olivier Lachaud

@since 0.9.1

Part of the \ref packageTopology.
  
This part of the manual describes how to represent and process
arbitrary cubical complexes. 

@note Collapse operation is a backport from \e ImaGene.
@cite ImaGene 

[TOC]


The following programs are related to this documentation:
cubical-complex-collapse.cpp, cubical-complex-illustrations.cpp,
testCubicalComplex.cpp, cubicalComplexThinning.cpp.

@section dgtal_ccomplex_sec1  Introduction to cubical complexes

We define a \b cubical \b complex \a C as a collection of cells living
in some Khalimsky space. Two cells of \a C are \b incident if and only
if they are incident in the Khalimsky space. A cubical complex
provides many services related to set of (unsigned) cells living in a
Khalimsky space: incidence operations, closure, star, link, cell set
operations, cell set relations, collapse operation. 


@note In opposition to the usual definition of simplical complex, we
do not require that all faces of a cell of \a C belong also to \a C.

@note A digital surface is not a cubical complex, because it does not
contain explicitely cells of dimension lower than n-1 and because it
requires some sort of orientation. However, it is easy to build an
ExplicitDigitalSurface or a LightExplicitDigitalSurface from a
CubicalComplex that stores a set of d-1-cells. You just need to create
a functor Surfel -> bool from the complex.

To create a cubical complex, we need to specify in which Khalimsky
space it lives and also, optionally, the type of container used for
storing cells. By default it is \c std::map but \c boost::unordered_map
or \c std::unordered_map is also possible.

\code
using namespace DGtal;
typedef KhalimskySpaceND<2,int>           KSpace; // ZxZ
typedef std::map<Cell, CubicalCellData>   Map;    // each cell must store at least a CubicalCellData
typedef CubicalComplex< KSpace, Map >     CC;     // the type of cubical complex

KSpace K;                                         // The cellular grid space
K.init( Point( 0,0 ), Point( 100,100 ), true );   // Limits of the grid space
CC complex( K );                                  // The complex
\endcode

Cells may be inserted through methods CubicalComplex::insertCell and
CubicalComplex::insertCells. Cells are easily created with a model of
CCellularGridSpaceND, for instance by specifying Khalimsky coordinates
(see \ref moduleCellularTopology). The small piece of code below creates a
ring around a pixel.

\code
typedef KSpace::Point                     Point;   // the type of each digital point.
complex.insertCell( 0, K.uCell( Point( 0, 0 ) ) ); // it is optional (but slightly faster)
complex.insertCell( 1, K.uCell( Point( 1, 0 ) ) ); // to specify the dimension of the cell
complex.insertCell( 0, K.uCell( Point( 2, 0 ) ) ); // at insertion.
complex.insertCell( 1, K.uCell( Point( 2, 1 ) ) );
complex.insertCell( 0, K.uCell( Point( 2, 2 ) ) );
complex.insertCell( 1, K.uCell( Point( 1, 2 ) ) );
complex.insertCell( 0, K.uCell( Point( 0, 2 ) ) );
complex.insertCell( 1, K.uCell( Point( 0, 1 ) ) );
trace.info() << "Euler number is " << complex.euler() << std::endl; // Should be 0
\endcode

Method CubicalComplex::insertCells accepts any range of iterators on
cells as input, hence you may for instance directly create a cubical
complex from a digital surface (see \ref moduleDigitalSurfaces).

You have also methods to remove some cells (CubicalComplex::eraseCell,
CubicalComplex::eraseCells) and methods to check if a cell belongs to
the complex (CubicalComplex::belongs).

@note You may use \c std::unordered_map or \c boost::unordered_map
instead of \c std::map. It is sometimes slightly faster. In this case
you need to define a hash function for cells. Either you provide one,
or you just include the file KhalimskyCellHashFunctions.h with:

\code
#include "DGtal/topology/KhalimskyCellHashFunctions.h"
\endcode

Last, there is a data associated with each cell of a complex. The data
type must either be CubicalCellData or a type that derives from
CubicalCellData. This data is used by the functions::collapse
operation, a function that is used to make a homotopic thinning of the
cubical complex. It may also be used by the user for other purposes,
like storing flags, coordinates or anything else associated to a
cell. Look at the documentation of CubicalCellData to see the default
stored flags and data.

@section dgtal_ccomplex_sec2  Displaying a 2D cubical complex

2D cubical complex are displayable on a Board (see \ref
moduleBoard2D). You just have to output the complex with the stream
operator of a board. You may adjust the default style before. The
following snippet creates and outputs a 2D complex.

@snippet topology/cubical-complex-illustrations.cpp cubical-complex-illustrations-X

@image html cubical-complex-illustrations-X.png "Displaying a 2D cubical complex on a Board."
@image latex cubical-complex-illustrations-X.png "Displaying a 2D cubical complex on a Board." width=5cm


@section dgtal_ccomplex_sec3  Incidence within a cubical complex

Two cells of a cubical complex are incident if and only if they are
incident in the Khalimsky space. From this relation, one can define
faces, co-faces, as well as the boundary and co-boundary of a cell.

You have methods to get the faces, the co-faces, the direct faces or
co-faces of a cell \a c, which are outputed with an output iterator:

- CubicalComplex::faces: outputs the faces of \a c, i.e. all cells of
  lesser dimension than \a c and incident to \a c. You may specify
  that the complex is closed (if it is, at least locally) to speed-up
  this function.

- CubicalComplex::directFaces: outputs the direct faces of \a c,
  i.e. all cells of dimension one less than \a c and incident to \a
  c. You may specify that the complex is closed (if it is, at least
  locally) to speed-up this function.

- CubicalComplex::coFaces: outputs the co-faces of \a c, i.e. all
  cells of greater dimension than \a c and incident to \a c. You may
  specify that the complex is open (if it is, at least locally) to
  speed-up this function.

- CubicalComplex::directCoFaces: outputs the direct co-faces of \a c,
  i.e. all cells of dimension one more than \a c and incident to \a
  c. You may specify that the complex is open (if it is, at least
  locally) to speed-up this function.

There are also versions of these methods that return the iterators on
these cells, if you need to access them directly in the complex
afterwards.

If you wish to get a vector of cells that contains all the proper
faces of some cell, then CubicalComplex::cellBoundary does the job. It
is in general slightly faster than using CubicalComplex::faces.

If you wish to get a vector of cells that contains all the proper
co-faces of some cell, then CubicalComplex::cellCoBoundary does the
job. Again it is in general slightly faster than using
CubicalComplex::coFaces.

@section dgtal_ccomplex_sec4  Closure, star and link of cells of a complex

You have three methods to compute these "subcomplexes" of a complex \a
X. In each case, you must give a subcomplex \a S of \a X which defines
the cells for which you wish to compute the closure, star or link.

@image html cubical-complex-illustrations-S.png "Cubical complex \a X in light grey and subcomplex \a S in dark grey."
@image latex cubical-complex-illustrations-S.png "Cubical complex \a X in light grey and subcomplex \a S in dark grey." width=5cm


- CubicalComplex::closure: returns the closure of the cells of \a S
  within this complex, i.e. the smallest subcomplex of \a X that
  contains each cell in \a S.

  @snippet topology/cubical-complex-illustrations.cpp cubical-complex-illustrations-closure

  @image html cubical-complex-illustrations-closure.png "Closure of \a S in \a X, displayed in red."
  @image latex cubical-complex-illustrations-closure.png "Closure of \a S in \a X, displayed in red." width=5cm

- CubicalComplex::star: returns the star of the cells of \a S
  within this complex, i.e. the union of the star (co-faces + cell) in \a X of every cell of \a S.

  @snippet topology/cubical-complex-illustrations.cpp cubical-complex-illustrations-star

  @image html cubical-complex-illustrations-star.png "Star of \a S in \a X, displayed in red."
  @image latex cubical-complex-illustrations-star.png "Star of \a S in \a X, displayed in red." width=5cm

- CubicalComplex::link: returns the link of the cells of \a S within
  this complex, i.e. the closure of the star of \a S in \a X minus the
  star of the closure of \a S in \a X.

  @snippet topology/cubical-complex-illustrations.cpp cubical-complex-illustrations-link

  @image html cubical-complex-illustrations-link.png "Link of \a S in \a X, displayed in green."
  @image latex cubical-complex-illustrations-link.png "Link of \a S in \a X, displayed in green." width=5cm

@section dgtal_ccomplex_sec5  Interior, boundary of a complex

You have three methods for computing the interior or the boundary of a
complex. Note that the boundary and the interior of a complex have an
empty intersection:

- CubicalComplex::interior: returns the (topological) interior of this complex \a X,
  i.e. the cells whose star are homeomorphic to a ball of maximal
  dimension.

  @snippet topology/cubical-complex-illustrations.cpp cubical-complex-illustrations-int

  @image html cubical-complex-illustrations-int.png "Interior of \a X, displayed in cyan."
  @image latex cubical-complex-illustrations-int.png "Interior of \a X, displayed in cyan." width=5cm
  
- CubicalComplex::boundary: returns the (topological) boundary of this
  complex \a X, i.e. the closure of \a X minus the interior of \a
  X. It is worthy to note that the boundary complex may contain cells
  that are not in \a X.

  @snippet topology/cubical-complex-illustrations.cpp cubical-complex-illustrations-bd

  @image html cubical-complex-illustrations-bd.png "Boundary of \a X, displayed in magenta."
  @image latex cubical-complex-illustrations-bd.png "Boundary of \a X, displayed in magenta." width=5cm
  
- CubicalComplex::getInteriorAndBoundary: returns both the boundary
  and the interior of X, faster then computing both separately if you
  need both.

@section dgtal_ccomplex_sec6 Set operations on complexes

A cubical complex can be seen as a generic container (it is a model of
boost::Container). Remember though that it is not a model of
boost::DefaultConstructible, because it must know a cellular grid
space to be valid. You could thus use operations defined in
DGtal::functions::setops (see \ref moduleSetFunctions) to perform set
operations and relations on complexes. However, to use the specificity
of the templated parameter CellContainer of CubicalComplex, it is
better to use operations and relations specifically defined for CubicalComplex, for instance as follows:

\code
CC X( K );            // creates the cubical complex
CC Xc  = ~X;          // closing of X
CC Xco = *Xc;         // opening of Xc (interior of X).
CC Xbd = Xc - Xco;    // boundary of X.
CC Xbd2= X.boundary();// same result for pure complex. 
\endcode

You have the following operations:
- operator|  : union of two complexes.
- operator|= : assignment union with another complex.
- operator&  : intersection of two complexes.
- operator&= : assignment intersection with another complex.
- operator-  : set difference of two complexes.
- operator-= : assignment set difference with another complex.
- operator^  : symmetric set difference of two complexes.
- operator^= : assignment symmetric set difference with another complex.
- operator~  : closing of a complex (smallest closed complex containg input complex).
- operator*  : opening of a complex (greatest open complex contained in input complex).

And the following relations:
- operator<=  : inclusion of two complexes (subset).
- operator>=  : inclusion of two complexes (supset).
- operator==  : equality of two complexes (both subset and supset).
- operator!=  : non-equality of two complexes (equal is false).

All these operations are at most \b linear in the sum of the number of
cells of input complexes, but they are generally exponential in the
dimension of complexes.

@section dgtal_ccomplex_sec7 Topological operations: closing, opening, collapsing a complex

You may \b close a complex (CubicalComplex::close), i.e. adding cells
such that the faces of every cell belongs to the complex. The complex
is then said to be closed. The updated complex corresponds to the
smallest closed complex containing it.

You may \b open a complex (CubicalComplex::open), i.e. removing cells
such that the co-faces of every cell belongs to the complex. The complex
is then said to be open. The updated complex corresponds to the
greatest open complex contained in it.

A more difficult topological operation is the \b collapse
(functions::collapse). A collapse transforms a complex by removing
cells, but at the end the complex has the same homotopy type as the
initial complex. The topology is preserved since only free pairs are
removed. The following snippet shows how to collapse complex X, while
keeping two cells \a fixed.

@snippet topology/cubical-complex-illustrations.cpp cubical-complex-illustrations-collapse

As you can see, when collapsing a complex X, you precise:
- The subcomplex of X that can be collapsed (should be closed,
  otherwise it is closed).
- A priority object that tells in which order cell are collapsed. If
  you use CubicalComplex::DefaultCellMapIteratorPriority, then the
  order follows the values stored in the data associated with each
  cell. You may use CubicalComplex::operator[] to change this value
  for each cell. The higher the value the sooner the cell is collapsed.
- Various hints to speed up the process.

Example topology/cubical-complex-collapse.cpp shows that it works also in 3D (in fact, it is nD).

@image html cubical-complex-illustrations-collapse.png "Homotopic collapse of \a X, displayed in orange, with fixed points marked in blue."
@image latex cubical-complex-illustrations-collapse.png "Homotopic collapse of \a X, displayed in orange, with fixed points marked in blue." width=5cm

@section dgtal_ccomplex_sec8 Thinning in cubical complexes

Thinning of digital objects in 2D and 3D with the guarantee that a final skeleton is thin can be achieved directly
only in \f$\mathbb{Z}^2\f$ for so-called well-composed images @cite Latecki:ER:95. Nevertheless, the problem can
 be solved for 3D objects represented by cubical complexes @cite Chaussard:IWCIA:09.

In @cite Chaussard:IWCIA:09 Chaussard and Couprie proposed a parallel directional thinning schemes in cubical complexes
 which allow to obtain a thin result. The implementation in DGtal works for both 2D and 3D cubical complexes.

The first scheme, so-called ParDirCollapse---parallel directional collapse---based on the idea such that free pairs
of faces are collapsed with respect to a fixed order of directions and order of face dimensions. In other words, in
each iteration we first collapse free pairs of a given direction starting with pairs of the highest dimension.
When there is no more free pairs to remove for a given direction then another direction is considered unless the
specified number of iterations is not reached . Note in the case of ParDirCollapse the thinness of the result
depends on the number of iterations.

The steps below present how to use ParDirCollapse in DGtal:

a) add those includes:
@snippet topology/cubicalComplexThinning.cpp include

b) create an instance of the class ParDirCollapse and initialize it with an instance of the Khalimsky space:
@snippet topology/cubicalComplexThinning.cpp create

c) attach a cubical complex \f$X\f$ to be thinned:
@snippet topology/cubicalComplexThinning.cpp init

d) run the algorithm (here two iterations):
@snippet topology/cubicalComplexThinning.cpp thinn

@image html ComplexBeforeThinning.png "The starting complex X before thinning." width=3cm
@image latex ComplexBeforeThinning.png "The starting complex X before thinning." width=3cm

@image html ParDirCollapse_2.png  "The starting complex X after two iterations of ParDirCollapse." width=5cm
@image latex ParDirCollapse_2.png  "The starting complex X after two iterations of ParDirCollapse." width=5cm

The second scheme, so-called SurfaceCollapse is similar to ParDirCollapse. The only difference is that after each
iteration of ParDirCollapse, faces of dimension \f$dim(X) - 1\f$ which are not included in faces of higher
dimension are marked as non-collapsible. The final output is guaranteed to be thin i.e. there is no faces of
\f$dim(X)\f$.

To run SurfaceCollapse we apply steps a--c for ParDirCollapse and then we apply the step below:

d) run the algorithm:
@snippet topology/cubicalComplexThinning.cpp thinnSurface

@image html ParDirCollapse_collapseSurface.png "The starting complex X after SurfaceCollapse." width=3cm
@image latex ParDirCollapse_collapseSurface.png "The starting complex X after SurfaceCollapse." width=3cm

The last scheme, so-called IsthmusCollapse is similar to SurfaceCollapse. The only difference is that after each
iteration of ParDirCollapse, we mark as non-collapsible those faces which are of dimension dim\f$(X) - 1\f$ and not
included in faces of higher dimension. Moreover, all sub-faces of dimension \f$dim(X) - 2\f$ of those faces
have to be not free -- are included in other faces of dimension \f$dim(X) - 1\f$. The final output is
guarantee to be thin i.e. there is no faces of \f$dim(X)\f$.

To run SurfaceCollapse we apply steps a--c for ParDirCollapse and then we apply the step below:

d) run the algorithm:
@snippet topology/cubicalComplexThinning.cpp thinnIsthmus

@image html ParDirCollapse_collapseIsthmus.png "The starting complex X after IsthmusCollapse." width=3cm
@image latex ParDirCollapse_collapseIsthmus.png "The starting complex X after IsthmusCollapse." width=3cm

For more details see the example topology/cubicalComplexThinning.cpp.

*/

}
