CellGeometry.h

 
#pragma once
 
#if defined(CellGeometry_RECURSES)
#error Recursive header files inclusion detected in CellGeometry.h
#else // defined(CellGeometry_RECURSES)
#define CellGeometry_RECURSES
 
#if !defined CellGeometry_h
#define CellGeometry_h
 
// Inclusions
#include <iostream>
#include <list>
#include <vector>
#include <string>
#include <unordered_set>
#include "DGtal/base/Common.h"
#include "DGtal/kernel/UnorderedSetByBlock.h"
#include "DGtal/kernel/PointHashFunctions.h"
#include "DGtal/topology/CCellularGridSpaceND.h"
#include "DGtal/topology/KhalimskySpaceND.h"
#include "DGtal/topology/KhalimskyCellHashFunctions.h"
#include "DGtal/geometry/volumes/BoundedLatticePolytope.h"
#include "DGtal/geometry/volumes/BoundedRationalPolytope.h"
 
namespace DGtal
{
 
  // template class CellGeometry
  template < typename TKSpace >
  class CellGeometry
  {
    BOOST_CONCEPT_ASSERT(( concepts::CCellularGridSpaceND< TKSpace > ));
 
  public:
    typedef CellGeometry<TKSpace>           Self;
    typedef TKSpace                         KSpace;
    typedef typename KSpace::Integer        Integer;
    typedef typename KSpace::Point          Point;
    typedef typename KSpace::Vector         Vector;
    typedef typename KSpace::Cell           Cell;
    typedef typename KSpace::Space          Space;
    typedef typename KSpace::Size           Size;
#ifdef WITH_BIGINTEGER
    typedef DGtal::BigInteger               BigInteger;
#else
    typedef DGtal::int64_t                  BigInteger;
#endif
    typedef DGtal::BoundedLatticePolytope< Space >  Polytope;
    typedef DGtal::BoundedLatticePolytope< Space >  LatticePolytope;
    typedef DGtal::BoundedRationalPolytope< Space > RationalPolytope;
 
    static const Dimension dimension = KSpace::dimension;
 
 
 
    ~CellGeometry() = default;
 
    CellGeometry();
 
    CellGeometry ( const Self & other ) = default;
 
    CellGeometry ( Self && other ) = default;
 
    CellGeometry ( const KSpace & K,
                   Dimension min_cell_dim = 0,
                   Dimension max_cell_dim = KSpace::dimension,
                   bool verbose = false );
 
    Self & operator= ( const Self & other ) = default;
 
    void init( const KSpace & K,
               Dimension min_cell_dim = 0,
               Dimension max_cell_dim = KSpace::dimension,
               bool verbose = false );
 
    // ----------------------- Cells services ------------------------------
  public:
 
    template <typename PointIterator>
    void addCellsTouchingPoints( PointIterator itB, PointIterator itE );
 
    template <typename PointelIterator>
    void addCellsTouchingPointels( PointelIterator itB, PointelIterator itE );
 
    void addCellsTouchingPolytopePoints( const LatticePolytope& polytope );
 
    void addCellsTouchingPolytopePoints( const RationalPolytope& polytope );
 
    void addCellsTouchingPolytope( const LatticePolytope& polytope );
 
    void addCellsTouchingPolytope( const RationalPolytope& polytope );
 
    CellGeometry& operator+=( const CellGeometry& other );
 
 
    // ----------------------- Accessor services ------------------------------
  public:
 
    Size nbCells() const;
    Size computeNbCells(const Dimension k ) const;
    Integer computeEuler() const;
    Dimension minCellDim() const;
    Dimension maxCellDim() const;
 
 
    // ----------------------- Cell services ------------------------------
  public:
 
    bool subset( const CellGeometry& other ) const;
 
    bool subset( const CellGeometry& other, const Dimension k ) const;
 
 
    // ----------------------- helper services ------------------------------
  public:
 
    std::vector< Point >
    getIntersectedKPoints( const LatticePolytope& polytope, const Dimension i ) const;
 
    std::vector< Point >
    getIntersectedKPoints( const RationalPolytope& polytope, const Dimension i ) const;
 
    std::vector< Point >
    getTouchedKPoints( const std::vector< Point >& points, const Dimension i ) const;
 
 
    std::vector< Cell >
    getIntersectedCells( const LatticePolytope& polytope, const Dimension i ) const;
 
    std::vector< Cell >
    getIntersectedCells( const RationalPolytope& polytope, const Dimension i ) const;
 
    std::vector< Cell >
    getTouchedCells( const std::vector< Point >& points, const Dimension i ) const;
 
 
    // ----------------------- Khalimsky point services -------------------------------
  public:
 
    static Dimension dim( const Point& kp );
 
 
    // ----------------------- Interface --------------------------------------
  public:
 
    void selfDisplay ( std::ostream & out ) const;
 
    bool isValid() const;
 
    std::string className() const;
 
 
    // ------------------------- Protected Datas ------------------------------
  protected:
 
    KSpace myK;
    UnorderedSetByBlock< Point,
                         Splitter< Point, uint64_t> > myKPoints;
    Dimension myMinCellDim;
    Dimension myMaxCellDim;
    bool myVerbose;
 
    // ------------------------- Private Datas --------------------------------
  private:
 
 
    // ------------------------- Internals ------------------------------------
  private:
    // Internal method for cheking if sorted range [it1,itE1] is a
    // subset of sorted range [it2, itE2].  Different from
    // std::includes since it performs exponential march and dichotomy
    // to walk faster along range [it1,itE1].
    template <typename RandomIterator>
    static
    bool includes( RandomIterator it2, RandomIterator itE2,
                   RandomIterator it1, RandomIterator itE1 );
 
  }; // end of class CellGeometry
 
 
  template <typename TKSpace>
  std::ostream&
  operator<< ( std::ostream & out,
               const CellGeometry<TKSpace> & object );
 
 
  template <typename TKSpace, int i, int N>
  struct CellGeometryFunctions
  {
    BOOST_CONCEPT_ASSERT(( concepts::CCellularGridSpaceND< TKSpace > ));
    typedef TKSpace                KSpace;
    typedef typename KSpace::Space Space;
    typedef typename KSpace::Cell  Cell;
    typedef typename KSpace::Point Point;
 
    template <typename PointelIterator>
    static
    std::unordered_set<typename KSpace::Cell>
    getIncidentCellsToPointels( const KSpace& K,
                                PointelIterator itB, PointelIterator itE )
    {
      std::unordered_set<typename KSpace::Cell> cells;
      if ( i == 0 )
        for ( auto it = itB; it != itE; ++it )
          cells.insert( *it );
      else
        for ( auto it = itB; it != itE; ++it )
          {
            auto pointel = *it;
            auto cofaces = K.uCoFaces( pointel );
            for ( auto&& f : cofaces )
              if ( K.uDim( f ) == i ) cells.insert( f );
          }
      return cells;
    }
 
    template <typename PointIterator>
    static
    std::unordered_set<typename KSpace::Cell>
    getIncidentCellsToPoints( const KSpace& K,
                              PointIterator itB, PointIterator itE )
    {
      std::unordered_set<typename KSpace::Cell> cells;
      if ( i == 0 )
        for ( auto it = itB; it != itE; ++it )
          cells.insert( K.uPointel( *it ) );
      else
        for ( auto it = itB; it != itE; ++it )
          {
            auto pointel = K.uPointel( *it );
            auto cofaces = K.uCoFaces( pointel );
            for ( auto&& f : cofaces )
              if ( K.uDim( f ) == i ) cells.insert( f );
          }
      return cells;
    }
 
    template <typename PointIterator>
    static
    UnorderedSetByBlock< typename KSpace::Point,
                         Splitter< typename KSpace::Point, uint64_t > >
    getIncidentKPointsToPoints( const KSpace& K,
                                PointIterator itB, PointIterator itE )
    {
      UnorderedSetByBlock< typename KSpace::Point,
                           Splitter< typename KSpace::Point, uint64_t > > kpoints;
      if ( i == 0 )
        for ( auto it = itB; it != itE; ++it )
          kpoints.insert( K.uKCoords( K.uPointel( *it ) ) );
      else
        for ( auto it = itB; it != itE; ++it )
          {
            auto pointel = K.uPointel( *it );
            auto cofaces = K.uCoFaces( pointel );
            for ( auto&& f : cofaces )
              if ( K.uDim( f ) == i ) kpoints.insert( K.uKCoords( f ) );
          }
      return kpoints;
    }
 
 
  }; // end struct CellGeometryFunctions
 
  template <typename TKSpace>
  struct CellGeometryFunctions< TKSpace, 1, 2 >
  {
    BOOST_CONCEPT_ASSERT(( concepts::CCellularGridSpaceND< TKSpace > ));
    BOOST_STATIC_ASSERT( TKSpace::dimension == 2 );
    typedef TKSpace                KSpace;
    typedef typename KSpace::Space Space;
    typedef typename KSpace::Cell  Cell;
 
    template <typename PointelIterator>
    static
    std::unordered_set<typename KSpace::Cell>
    getIncidentCellsToPointels( const KSpace& K,
                                PointelIterator itB, PointelIterator itE )
    {
      std::unordered_set<typename KSpace::Cell> cells;
      for ( auto it = itB; it != itE; ++it )
        {
          auto pointel = *it;
          cells.insert( K.uIncident( pointel, 0, true ) );
          cells.insert( K.uIncident( pointel, 0, false ) );
          cells.insert( K.uIncident( pointel, 1, true ) );
          cells.insert( K.uIncident( pointel, 1, false ) );
        }
      return cells;
    }
 
    template <typename PointIterator>
    static
    std::unordered_set<typename KSpace::Cell>
    getIncidentCellsToPoints( const KSpace& K,
                              PointIterator itB, PointIterator itE )
    {
      std::unordered_set<typename KSpace::Cell> cells;
      for ( auto it = itB; it != itE; ++it )
        {
          auto pointel = K.uPointel( *it );
          cells.insert( K.uIncident( pointel, 0, true ) );
          cells.insert( K.uIncident( pointel, 0, false ) );
          cells.insert( K.uIncident( pointel, 1, true ) );
          cells.insert( K.uIncident( pointel, 1, false ) );
        }
      return cells;
    }
 
    template <typename PointIterator>
    static
    UnorderedSetByBlock< typename KSpace::Point,
                         Splitter< typename KSpace::Point, uint64_t > >
    getIncidentKPointsToPoints( const KSpace& K,
                                PointIterator itB, PointIterator itE )
    {
      UnorderedSetByBlock< typename KSpace::Point,
                           Splitter< typename KSpace::Point, uint64_t > > kpoints;
      for ( auto it = itB; it != itE; ++it )
        {
          auto kp = K.uKCoords( K.uPointel( *it ) );
          kpoints.emplace( kp[ 0 ]    , kp[ 1 ] - 1 );
          kpoints.emplace( kp[ 0 ] - 1, kp[ 1 ]     );
          kpoints.emplace( kp[ 0 ] + 1, kp[ 1 ]     );
          kpoints.emplace( kp[ 0 ]    , kp[ 1 ] + 1 );
        }
      return kpoints;
    }
 
  }; // end struct CellGeometryFunctions
 
  template <typename TKSpace>
  struct CellGeometryFunctions< TKSpace, 1, 3 >
  {
    BOOST_CONCEPT_ASSERT(( concepts::CCellularGridSpaceND< TKSpace > ));
    BOOST_STATIC_ASSERT( TKSpace::dimension == 3 );
    typedef TKSpace                KSpace;
    typedef typename KSpace::Space Space;
    typedef typename KSpace::Cell  Cell;
 
    template <typename PointelIterator>
    static
    std::unordered_set<typename KSpace::Cell>
    getIncidentCellsToPointels( const KSpace& K,
                                PointelIterator itB, PointelIterator itE )
    {
      std::unordered_set<typename KSpace::Cell> cells;
      for ( auto it = itB; it != itE; ++it )
        {
          auto pointel = *it;
          cells.insert( K.uIncident( pointel, 0, true ) );
          cells.insert( K.uIncident( pointel, 0, false ) );
          cells.insert( K.uIncident( pointel, 1, true ) );
          cells.insert( K.uIncident( pointel, 1, false ) );
          cells.insert( K.uIncident( pointel, 2, true ) );
          cells.insert( K.uIncident( pointel, 2, false ) );
        }
      return cells;
    }
 
    template <typename PointIterator>
    static
    std::unordered_set<typename KSpace::Cell>
    getIncidentCellsToPoints( const KSpace& K,
                              PointIterator itB, PointIterator itE )
    {
      std::cout << "<1,3> specialization" << std::endl;
      std::unordered_set<typename KSpace::Cell> cells;
      for ( auto it = itB; it != itE; ++it )
        {
          auto pointel = K.uPointel( *it );
          cells.insert( K.uIncident( pointel, 0, true ) );
          cells.insert( K.uIncident( pointel, 0, false ) );
          cells.insert( K.uIncident( pointel, 1, true ) );
          cells.insert( K.uIncident( pointel, 1, false ) );
          cells.insert( K.uIncident( pointel, 2, true ) );
          cells.insert( K.uIncident( pointel, 2, false ) );
        }
      return cells;
    }
 
    template <typename PointIterator>
    static
    UnorderedSetByBlock< typename KSpace::Point,
                         Splitter< typename KSpace::Point, uint64_t > >
    getIncidentKPointsToPoints( const KSpace& K,
                                PointIterator itB, PointIterator itE )
    {
      UnorderedSetByBlock< typename KSpace::Point,
                           Splitter< typename KSpace::Point, uint64_t > > kpoints;
      for ( auto it = itB; it != itE; ++it )
        {
          auto kp = K.uKCoords( K.uPointel( *it ) );
          kpoints.emplace( kp[ 0 ] - 1, kp[ 1 ]    , kp[ 2 ]     );
          kpoints.emplace( kp[ 0 ] + 1, kp[ 1 ]    , kp[ 2 ]     );
          kpoints.emplace( kp[ 0 ]    , kp[ 1 ] - 1, kp[ 2 ]     );
          kpoints.emplace( kp[ 0 ]    , kp[ 1 ] + 1, kp[ 2 ]     );
          kpoints.emplace( kp[ 0 ]    , kp[ 1 ]    , kp[ 2 ] - 1 );
          kpoints.emplace( kp[ 0 ]    , kp[ 1 ]    , kp[ 2 ] + 1 );
        }
      return kpoints;
    }
 
  }; // end struct CellGeometryFunctions
 
  template <typename TKSpace>
  struct CellGeometryFunctions< TKSpace, 2, 2 >
  {
    BOOST_CONCEPT_ASSERT(( concepts::CCellularGridSpaceND< TKSpace > ));
    BOOST_STATIC_ASSERT( TKSpace::dimension == 2 );
    typedef TKSpace                KSpace;
    typedef typename KSpace::Space Space;
    typedef typename KSpace::Cell  Cell;
 
    template <typename PointelIterator>
    static
    std::unordered_set<typename KSpace::Cell>
    getIncidentCellsToPointels( const KSpace& K,
                                PointelIterator itB, PointelIterator itE )
    {
      std::unordered_set<typename KSpace::Cell> cells;
      for ( auto it = itB; it != itE; ++it )
        {
          auto pointel = *it;
          auto linelxp = K.uIncident( pointel, 0, true );
          auto linelxm = K.uIncident( pointel, 0, false );
          cells.insert( K.uIncident( linelxp, 1, true ) );
          cells.insert( K.uIncident( linelxp, 1, false ) );
          cells.insert( K.uIncident( linelxm, 1, true ) );
          cells.insert( K.uIncident( linelxm, 1, false ) );
        }
      return cells;
    }
 
    template <typename PointIterator>
    static
    std::unordered_set<typename KSpace::Cell>
    getIncidentCellsToPoints( const KSpace& K,
                              PointIterator itB, PointIterator itE )
    {
      std::cout << "<2,2> specialization" << std::endl;
      std::unordered_set<typename KSpace::Cell> cells;
      for ( auto it = itB; it != itE; ++it )
        {
          auto pointel = K.uPointel( *it );
          auto linelxp = K.uIncident( pointel, 0, true );
          auto linelxm = K.uIncident( pointel, 0, false );
          cells.insert( K.uIncident( linelxp, 1, true ) );
          cells.insert( K.uIncident( linelxp, 1, false ) );
          cells.insert( K.uIncident( linelxm, 1, true ) );
          cells.insert( K.uIncident( linelxm, 1, false ) );
        }
      return cells;
    }
 
    template <typename PointIterator>
    static
    UnorderedSetByBlock< typename KSpace::Point,
                         Splitter< typename KSpace::Point, uint64_t> >
    getIncidentKPointsToPoints( const KSpace& K,
                                PointIterator itB, PointIterator itE )
    {
      UnorderedSetByBlock< typename KSpace::Point,
                           Splitter< typename KSpace::Point, uint64_t> > kpoints;
      for ( auto it = itB; it != itE; ++it )
        {
          auto kp = K.uKCoords( K.uPointel( *it ) );
          kpoints.emplace( kp[ 0 ] - 1, kp[ 1 ] - 1 );
          kpoints.emplace( kp[ 0 ] + 1, kp[ 1 ] - 1 );
          kpoints.emplace( kp[ 0 ] - 1, kp[ 1 ] + 1 );
          kpoints.emplace( kp[ 0 ] + 1, kp[ 1 ] + 1 );
        }
      return kpoints;
    }
 
  }; // end struct CellGeometryFunctions
 
  template <typename TKSpace>
  struct CellGeometryFunctions< TKSpace, 2, 3 >
  {
    BOOST_CONCEPT_ASSERT(( concepts::CCellularGridSpaceND< TKSpace > ));
    BOOST_STATIC_ASSERT( TKSpace::dimension == 2 );
    typedef TKSpace                KSpace;
    typedef typename KSpace::Space Space;
    typedef typename KSpace::Cell  Cell;
 
    template <typename PointelIterator>
    static
    std::unordered_set<typename KSpace::Cell>
    getIncidentCellsToPointels( const KSpace& K,
                                PointelIterator itB, PointelIterator itE )
    {
      std::unordered_set<typename KSpace::Cell> cells;
      for ( auto it = itB; it != itE; ++it )
        {
          auto pointel = *it;
          auto linelxp = K.uIncident( pointel, 0, true );
          auto linelxm = K.uIncident( pointel, 0, false );
          auto linelyp = K.uIncident( pointel, 1, true );
          auto linelym = K.uIncident( pointel, 1, false );
          cells.insert( K.uIncident( linelxp, 1, true ) );
          cells.insert( K.uIncident( linelxp, 1, false ) );
          cells.insert( K.uIncident( linelxp, 2, true ) );
          cells.insert( K.uIncident( linelxp, 2, false ) );
          cells.insert( K.uIncident( linelxm, 1, true ) );
          cells.insert( K.uIncident( linelxm, 1, false ) );
          cells.insert( K.uIncident( linelxm, 2, true ) );
          cells.insert( K.uIncident( linelxm, 2, false ) );
          cells.insert( K.uIncident( linelyp, 2, true ) );
          cells.insert( K.uIncident( linelyp, 2, false ) );
          cells.insert( K.uIncident( linelym, 2, true ) );
          cells.insert( K.uIncident( linelym, 2, false ) );
        }
      return cells;
    }
 
    template <typename PointIterator>
    static
    std::unordered_set<typename KSpace::Cell>
    getIncidentCellsToPoints( const KSpace& K,
                              PointIterator itB, PointIterator itE )
    {
      std::cout << "<2,3> specialization" << std::endl;
      std::unordered_set<typename KSpace::Cell> cells;
      for ( auto it = itB; it != itE; ++it )
        {
          auto pointel = K.uPointel( *it );
          auto linelxp = K.uIncident( pointel, 0, true );
          auto linelxm = K.uIncident( pointel, 0, false );
          auto linelyp = K.uIncident( pointel, 1, true );
          auto linelym = K.uIncident( pointel, 1, false );
          cells.insert( K.uIncident( linelxp, 1, true ) );
          cells.insert( K.uIncident( linelxp, 1, false ) );
          cells.insert( K.uIncident( linelxp, 2, true ) );
          cells.insert( K.uIncident( linelxp, 2, false ) );
          cells.insert( K.uIncident( linelxm, 1, true ) );
          cells.insert( K.uIncident( linelxm, 1, false ) );
          cells.insert( K.uIncident( linelxm, 2, true ) );
          cells.insert( K.uIncident( linelxm, 2, false ) );
          cells.insert( K.uIncident( linelyp, 2, true ) );
          cells.insert( K.uIncident( linelyp, 2, false ) );
          cells.insert( K.uIncident( linelym, 2, true ) );
          cells.insert( K.uIncident( linelym, 2, false ) );
        }
      return cells;
    }
 
    template <typename PointIterator>
    static
    UnorderedSetByBlock< typename KSpace::Point,
                         Splitter< typename KSpace::Point, uint64_t> >
    getIncidentKPointsToPoints( const KSpace& K,
                                PointIterator itB, PointIterator itE )
    {
      UnorderedSetByBlock< typename KSpace::Point,
                           Splitter< typename KSpace::Point, uint64_t> > kpoints;
      for ( auto it = itB; it != itE; ++it )
        {
          auto kp = K.uKCoords( K.uPointel( *it ) );
          kpoints.emplace( kp[ 0 ] - 1, kp[ 1 ] - 1, kp[ 2 ]     );
          kpoints.emplace( kp[ 0 ] + 1, kp[ 1 ] - 1, kp[ 2 ]     );
          kpoints.emplace( kp[ 0 ] - 1, kp[ 1 ] + 1, kp[ 2 ]     );
          kpoints.emplace( kp[ 0 ] + 1, kp[ 1 ] + 1, kp[ 2 ]     );
          kpoints.emplace( kp[ 0 ] - 1, kp[ 1 ]    , kp[ 2 ] - 1 );
          kpoints.emplace( kp[ 0 ] + 1, kp[ 1 ]    , kp[ 2 ] - 1 );
          kpoints.emplace( kp[ 0 ] - 1, kp[ 1 ]    , kp[ 2 ] + 1 );
          kpoints.emplace( kp[ 0 ] + 1, kp[ 1 ]    , kp[ 2 ] + 1 );
          kpoints.emplace( kp[ 0 ]    , kp[ 1 ] - 1, kp[ 2 ] - 1 );
          kpoints.emplace( kp[ 0 ]    , kp[ 1 ] + 1, kp[ 2 ] - 1 );
          kpoints.emplace( kp[ 0 ]    , kp[ 1 ] - 1, kp[ 2 ] + 1 );
          kpoints.emplace( kp[ 0 ]    , kp[ 1 ] + 1, kp[ 2 ] + 1 );
        }
      return kpoints;
    }
 
  }; // end struct CellGeometryFunctions
 
  template <typename TKSpace>
  struct CellGeometryFunctions< TKSpace, 3, 3 >
  {
    BOOST_CONCEPT_ASSERT(( concepts::CCellularGridSpaceND< TKSpace > ));
    BOOST_STATIC_ASSERT( TKSpace::dimension == 3 );
    typedef TKSpace                KSpace;
    typedef typename KSpace::Space Space;
    typedef typename KSpace::Cell  Cell;
 
    template <typename PointelIterator>
    static
    std::unordered_set<typename KSpace::Cell>
    getIncidentCellsToPointels( const KSpace& K,
                                PointelIterator itB, PointelIterator itE )
    {
      std::unordered_set<typename KSpace::Cell> cells;
      for ( auto it = itB; it != itE; ++it )
        {
          auto pointel = *it;
          auto linelxp = K.uIncident( pointel, 0, true );
          auto linelxm = K.uIncident( pointel, 0, false );
          auto surfxpyp = K.uIncident( linelxp, 1, true );
          auto surfxpym = K.uIncident( linelxp, 1, false );
          auto surfxmyp = K.uIncident( linelxm, 1, true );
          auto surfxmym = K.uIncident( linelxm, 1, false );
          cells.insert( K.uIncident( surfxpyp, 2, true ) );
          cells.insert( K.uIncident( surfxpyp, 2, false ) );
          cells.insert( K.uIncident( surfxpym, 2, true ) );
          cells.insert( K.uIncident( surfxpym, 2, false ) );
          cells.insert( K.uIncident( surfxmyp, 2, true ) );
          cells.insert( K.uIncident( surfxmyp, 2, false ) );
          cells.insert( K.uIncident( surfxmym, 2, true ) );
          cells.insert( K.uIncident( surfxmym, 2, false ) );
        }
      return cells;
    }
 
    template <typename PointIterator>
    static
    std::unordered_set<typename KSpace::Cell>
    getIncidentCellsToPoints( const KSpace& K,
                              PointIterator itB, PointIterator itE )
    {
      std::unordered_set<typename KSpace::Cell> cells;
      for ( auto it = itB; it != itE; ++it )
        {
          auto pointel = K.uPointel( *it );
          auto linelxp = K.uIncident( pointel, 0, true );
          auto linelxm = K.uIncident( pointel, 0, false );
          auto surfxpyp = K.uIncident( linelxp, 1, true );
          auto surfxpym = K.uIncident( linelxp, 1, false );
          auto surfxmyp = K.uIncident( linelxm, 1, true );
          auto surfxmym = K.uIncident( linelxm, 1, false );
          cells.insert( K.uIncident( surfxpyp, 2, true ) );
          cells.insert( K.uIncident( surfxpyp, 2, false ) );
          cells.insert( K.uIncident( surfxpym, 2, true ) );
          cells.insert( K.uIncident( surfxpym, 2, false ) );
          cells.insert( K.uIncident( surfxmyp, 2, true ) );
          cells.insert( K.uIncident( surfxmyp, 2, false ) );
          cells.insert( K.uIncident( surfxmym, 2, true ) );
          cells.insert( K.uIncident( surfxmym, 2, false ) );
        }
      return cells;
    }
 
    template <typename PointIterator>
    static
    UnorderedSetByBlock< typename KSpace::Point,
                         Splitter< typename KSpace::Point, uint64_t> >
    getIncidentKPointsToPoints( const KSpace& K,
                                PointIterator itB, PointIterator itE )
    {
      UnorderedSetByBlock< typename KSpace::Point,
                           Splitter< typename KSpace::Point, uint64_t> > kpoints;
      for ( auto it = itB; it != itE; ++it )
        {
          auto kp = K.uKCoords( K.uPointel( *it ) );
          kpoints.emplace( kp[ 0 ] - 1, kp[ 1 ] - 1, kp[ 2 ] - 1 );
          kpoints.emplace( kp[ 0 ] + 1, kp[ 1 ] - 1, kp[ 2 ] - 1 );
          kpoints.emplace( kp[ 0 ] - 1, kp[ 1 ] + 1, kp[ 2 ] - 1 );
          kpoints.emplace( kp[ 0 ] + 1, kp[ 1 ] + 1, kp[ 2 ] - 1 );
          kpoints.emplace( kp[ 0 ] - 1, kp[ 1 ] - 1, kp[ 2 ] + 1 );
          kpoints.emplace( kp[ 0 ] + 1, kp[ 1 ] - 1, kp[ 2 ] + 1 );
          kpoints.emplace( kp[ 0 ] - 1, kp[ 1 ] + 1, kp[ 2 ] + 1 );
          kpoints.emplace( kp[ 0 ] + 1, kp[ 1 ] + 1, kp[ 2 ] + 1 );
        }
      return kpoints;
    }
 
  }; // end struct CellGeometryFunctions
 
 
} // namespace DGtal
 
 
// Includes inline functions.
#include "CellGeometry.ih"
 
//                                                                           //
 
#endif // !defined CellGeometry_h
 
#undef CellGeometry_RECURSES
#endif // else defined(CellGeometry_RECURSES)