CubicalComplexFunctions.ih

/**
 *  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 CubicalComplexFunctions.ih
 * @author Jacques-Olivier Lachaud (\c jacques-olivier.lachaud@univ-savoie.fr )
 * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
 *
 * @date 2015/11/24
 *
 * Implementation of inline methods defined in CubicalComplexFunctions.h
 *
 * This file is part of the DGtal library.
 */
 
 
//////////////////////////////////////////////////////////////////////////////
#include <cstdlib>
#include "DGtal/kernel/domains/HyperRectDomain.h"
#include "DGtal/topology/DigitalTopology.h"
#include "DGtal/topology/helpers/NeighborhoodConfigurationsHelper.h"
//////////////////////////////////////////////////////////////////////////////
 
///////////////////////////////////////////////////////////////////////////////
// IMPLEMENTATION of inline methods.
///////////////////////////////////////////////////////////////////////////////
 
//-----------------------------------------------------------------------------
template <typename TKSpace, typename TCellContainer,
          typename CellConstIterator,
          typename CellMapIteratorPriority >
DGtal::uint64_t
DGtal::functions::
collapse( CubicalComplex< TKSpace, TCellContainer > & K,
          CellConstIterator S_itB, CellConstIterator S_itE,
          const CellMapIteratorPriority& priority,
          bool hintIsSClosed, bool hintIsKClosed,
          bool verbose )
{
  using namespace std;
  typedef CubicalComplex< TKSpace, TCellContainer > CC;
  typedef typename CC::Cell                         Cell;
  typedef typename CC::CellType                     CellType;
  typedef typename CC::CellMapIterator              CellMapIterator;
  typedef vector< CellMapIterator >                 CMIVector;
  typedef typename CMIVector::const_iterator        CMIVectorConstIterator;
  // NB : a maximal k-cell is collapsible if it has a free incident k-1-cell.
  Dimension n = K.dim();
  CMIVector S;            // stores the cells to process
  CMIVector Q_low;        // stores the iterators on direct faces of the maximal cell.
  CMIVector Q_collapsible;// stores collapsible cells in order to clean them at the end.
  CellMapIterator        it_cell;   // generic iterator on a cell.
  CellMapIterator        it_cell_c; // points on c in the free pair (c,d)
  CellMapIterator        it_cell_d; // points on d in the free pair (c,d)
  CMIVectorConstIterator itlow;
  CMIVectorConstIterator itqup;
 
  if ( verbose ) trace.info() << "[CC::collapse]-+ tag collapsible elements... " << flush;
  // Restricts the set of elements that are collapsible.
  if ( hintIsSClosed )
    for ( CellConstIterator S_it = S_itB; S_it != S_itE; ++S_it )
      {
        Cell c           = *S_it;
        Dimension k      = K.dim( c );
        it_cell          = K.findCell( k, c );
        uint32_t& ccdata = it_cell->second.data;
        ASSERT( it_cell != K.end( k ) );
        S.push_back( it_cell );
        if ( ! ( ccdata & (CC::FIXED | CC::COLLAPSIBLE ) ) )
          {
            ccdata |= CC::COLLAPSIBLE;
            Q_collapsible.push_back( it_cell );
          }
      }
  else // not ( hintIsSClosed )
    for ( CellConstIterator S_it = S_itB; S_it != S_itE; ++S_it )
      {
        Cell c           = *S_it;
        Dimension k      = K.dim( c );
        it_cell          = K.findCell( k, c );
        uint32_t& ccdata = it_cell->second.data;
        ASSERT( it_cell != K.end( k ) );
        S.push_back( it_cell );
        if ( ! ( ccdata & (CC::FIXED | CC::COLLAPSIBLE ) ) )
          {
            ccdata |= CC::COLLAPSIBLE;
            Q_collapsible.push_back( it_cell );
          }
        vector<Cell> cells;
        back_insert_iterator< vector<Cell> > back_it( cells );
        K.faces( back_it, c, hintIsKClosed );
        for ( typename vector<Cell>::const_iterator
                it = cells.begin(), itE = cells.end(); it != itE; ++it )
          {
            it_cell           = K.findCell( *it );
            uint32_t& ccdata2 = it_cell->second.data;
            if ( ! ( ccdata2 & (CC::FIXED | CC::COLLAPSIBLE ) ) )
              {
                ccdata2 |= CC::COLLAPSIBLE;
                Q_collapsible.push_back( it_cell );
              }
          }
      }
  if ( verbose ) trace.info() << " " << Q_collapsible.size() << " found." << endl;
 
  // Fill queue
  priority_queue<CellMapIterator, CMIVector, CellMapIteratorPriority> PQ( priority );
 
  if ( verbose ) trace.info() << "[CC::collapse]-+ entering collapsing loop. " << endl;
  uint64_t nb_pass     = 0;
  uint64_t nb_examined = 0;
  uint64_t nb_removed  = 0;
  while ( ! S.empty() )
    {
      for ( CMIVectorConstIterator it = S.begin(), itE = S.end();
            it != itE; ++it )
        {
          PQ.push( *it );
          (*it)->second.data |= CC::USER1;
        }
      S.clear();
      if ( verbose ) trace.info() << "[CC::collapse]---+ Pass " << ++nb_pass
                                  << ", Card(PQ)=" << PQ.size() << " elements, "
                                  << "nb_exam=" << nb_examined << endl;
 
      // Try to collapse elements according to priority queue.
      while ( ! PQ.empty() )
        {
          // Get top element.
          CellMapIterator itcur = PQ.top();
          uint32_t& cur_data    = itcur->second.data;
          PQ.pop();
          ++nb_examined;
          // Check if the cell is removable
          if ( ( cur_data & CC::REMOVED ) || ( ! ( cur_data & CC::COLLAPSIBLE ) ) )
            continue;
          // Check if the cell was several time in the queue and is already processed.
          if ( ! ( cur_data & CC::USER1 ) )
            continue;
          ASSERT( cur_data & CC::USER1 );
          cur_data             &= ~CC::USER1;
 
          // Cell may be removable.
          // Check if it is a maximal cell
          CellMapIterator itup;
          const Cell & cur_c  = itcur->first;
          CellType cur_c_type = K.computeCellType( cur_c, itup, n );
          bool found_pair     = false;
          // trace.info() << "  - Cell " << cur_c << " Dim=" << dim( cur_c ) << " Type=" << cur_c_type << std::endl;
          if ( cur_c_type == CC::Maximal )
            { // maximal cell... must find a free face
              // check faces to find a free face.
              back_insert_iterator< CMIVector > back_it( Q_low );
              K.directFacesIterators( back_it, cur_c );
              bool best_free_face_found = false;
              CellMapIterator best_free_face_it;
              for ( CMIVectorConstIterator it = Q_low.begin(), itE = Q_low.end();
                    it != itE; ++it )
                {
                  CellMapIterator low_ic = *it;
                  uint32_t& data         = low_ic->second.data;
                  // trace.info() << "    + Cell " << low_ic->first << " data=" << data << std::endl;
                  if ( ( data & CC::REMOVED ) || ! ( data & CC::COLLAPSIBLE ) ) continue;
                  const Cell& cur_d      = low_ic->first;
                  CellType cur_d_type    = K.computeCellType( cur_d, itup, n );
                  // trace.info() << "      + Type=" << cur_d_type << std::endl;
                  if ( cur_d_type == CC::Free )
                    { // found a free n-1-face ic
                      if ( ( ! best_free_face_found )
                           || ( ! priority( low_ic, best_free_face_it ) ) )
                        {
                          best_free_face_it    = low_ic;
                          best_free_face_found = true;
                        }
                    }
                }
              if ( best_free_face_found )
                {
                  // delete c and ic.
                  found_pair = true;
                  it_cell_c  = itcur;
                  it_cell_d  = best_free_face_it;
                  // Q_low already contains cells that should be
                  // checked again
                }
            }
          else if ( cur_c_type == CC::Free )
            { // free face... check that its 1-up-incident face is maximal.
              CellMapIterator it_up_up;
              const Cell& cur_d   = itup->first;
              CellType cur_d_type = K.computeCellType( cur_d, it_up_up, n );
              if ( cur_d_type == CC::Maximal )
                { // found a maximal face.
                  found_pair = true;
                  it_cell_c  = itup;
                  it_cell_d  = itcur;
                  // Q_low will contain cells that should be checked
                  // again
                  back_insert_iterator< CMIVector > back_it( Q_low );
                  K.directFacesIterators( back_it, it_cell_c->first );
                }
            }
          if ( found_pair )
            { // If found, remove pair from complex (logical removal).
              it_cell_c->second.data |= CC::REMOVED;
              it_cell_d->second.data |= CC::REMOVED;
              nb_removed             += 2;
              // Incident cells have to be checked again.
              for ( CMIVectorConstIterator it = Q_low.begin(), itE = Q_low.end();
                    it != itE; ++it )
                {
                  it_cell             = *it;
                  uint32_t& data_qlow = it_cell->second.data;
                  if ( ( ! ( data_qlow & CC::REMOVED ) )
                       && ( data_qlow & CC::COLLAPSIBLE )
                       && ( ! ( data_qlow & CC::USER1 ) ) )
                    {
                      S.push_back( it_cell );
                    }
                }
            }
          Q_low.clear();
        } // while ( ! PQ.empty() )
    } // while ( ! S.empty() )
 
  if ( verbose ) trace.info() << "[CC::collapse]-+ cleaning complex." << std::endl;
 
  // Now clean the complex so that removed cells are effectively
  // removed and no more cell is tagged as collapsible.
  for ( CMIVectorConstIterator it = Q_collapsible.begin(), itE = Q_collapsible.end();
        it != itE; ++it )
    {
      CellMapIterator cmIt  = *it;
      uint32_t& cur_data    = cmIt->second.data;
      if ( cur_data & CC::REMOVED ) K.eraseCell( cmIt );
      else                          cur_data &= ~CC::COLLAPSIBLE;
      // if ( (*it)->second.data & CC::REMOVED )
      //   K.eraseCell( *it );
      // else
      //   (*it)->second.data &= ~CC::COLLAPSIBLE;
    }
  return nb_removed;
}
 
 
//-----------------------------------------------------------------------------
template <typename TKSpace, typename TCellContainer,
          typename BdryCellOutputIterator,
          typename InnerCellOutputIterator>
void
DGtal::functions::
filterCellsWithinBounds( const CubicalComplex< TKSpace, TCellContainer > & K,
                         const typename TKSpace::Point& kLow,
                         const typename TKSpace::Point& kUp,
                         BdryCellOutputIterator itBdry,
                         InnerCellOutputIterator itInner )
{
  typedef CubicalComplex< TKSpace, TCellContainer > CC;
  typedef typename CC::Cell                         Cell;
  typedef typename CC::Point                        Point;
  typedef typename CC::CellMapConstIterator         CellMapConstIterator;
  Dimension d = K.dim();
  for ( Dimension i = 0; i <= d; ++i )
    {
      for ( CellMapConstIterator it = K.begin( i ), itE = K.end( i ); it != itE; ++it )
        {
          Cell cell = it->first;
          Point kCell = K.space().uKCoords( cell );
          if ( ( kLow.inf( kCell ) == kLow ) && ( kUp.sup( kCell ) == kUp ) )
            { // Inside or on boundary.
              bool bdry = false;
              for ( Dimension j = 0; j < Point::dimension; ++j )
                {
                  if ( ( kCell[ j ] == kLow[ j ] ) || ( kCell[ j ] == kUp[ j ] ) )
                    {
                      bdry = true;
                      break;
                    }
                }
              if ( bdry ) *itBdry++  = cell;
              else        *itInner++ = cell;
            }
        }
    }
}
 
//-----------------------------------------------------------------------------
template <typename TObject, typename TKSpace, typename TCellContainer>
std::unique_ptr<TObject>
DGtal::functions::
objectFromSpels(
                const CubicalComplex< TKSpace, TCellContainer > & C)
{
 
  const auto & cs = C.space();
  // Get domain of C KSpace.
  typename TObject::Domain domain(cs.lowerBound(), cs.upperBound());
  typename TObject::DigitalSet dset(domain);
  for(auto it = C.begin(TKSpace::DIM);
      it!= C.end(TKSpace::DIM) ; ++it)
    dset.insertNew(cs.uCoords(it->first));
 
  typename TObject::ForegroundAdjacency fAdj;
  typename TObject::BackgroundAdjacency bAdj;
  typename TObject::DigitalTopology topo(fAdj,bAdj, JORDAN_DT);
  std::unique_ptr<TObject> obj( new TObject(topo, dset) );
  return obj;
}
 
template<typename TComplex>
DGtal::NeighborhoodConfiguration
DGtal::functions::
getSpelNeighborhoodConfigurationOccupancy
( const          TComplex        & input_complex,
  const typename TComplex::Point & center,
  const std::unordered_map<
  typename TComplex::Point, NeighborhoodConfiguration> & mapPointToMask )
{
  using Point  = typename TComplex::Point;
  using Space  = typename TComplex::Space;
  using Domain = HyperRectDomain< Space >;
 
  Point p1 = Point::diagonal( -1 );
  Point p2 = Point::diagonal(  1 );
  Point c = Point::diagonal( 0 );
  Domain domain( p1, p2 );
 
  using KPreSpace = typename TComplex::KSpace::PreCellularGridSpace;
  const auto & kspace = input_complex.space();
  const auto & not_found( input_complex.end(3) );
  NeighborhoodConfiguration cfg{0};
  for ( auto it = domain.begin(); it != domain.end(); ++it ) {
    if( *it == c) continue;
    const auto pre_cell = KPreSpace::uSpel(center + *it);
    if (input_complex.belongs(pre_cell) &&
        input_complex.findCell(3, kspace.uCell(pre_cell)) != not_found )
      cfg |= mapPointToMask.at(*it) ;
  }
  return cfg;
}
//                                                                           //
///////////////////////////////////////////////////////////////////////////////