VoxelComplexFunctions.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 VoxelComplexFunctions.ih
 * @author Pablo Hernandez-Cerdan (\c pablo.hernandez.cerdan@outlook.com)
 * Institute of Fundamental Sciences. Massey University.
 * Palmerston North, New Zealand
 *
 * @date 2018/01/01
 *
 * Implementation of inline methods defined in VoxelComplexFunctions.h
 *
 * This file is part of the DGtal library.
 */
 
 
//////////////////////////////////////////////////////////////////////////////
#include <cstdlib>
#include <DGtal/topology/DigitalTopology.h>
#include <random>
//////////////////////////////////////////////////////////////////////////////
 
///////////////////////////////////////////////////////////////////////////////
// IMPLEMENTATION of inline methods.
///////////////////////////////////////////////////////////////////////////////
 
//-----------------------------------------------------------------------------
template < typename TComplex >
TComplex
DGtal::functions::
asymetricThinningScheme(
    TComplex & vc ,
    std::function<
    std::pair<typename TComplex::Cell, typename TComplex::Data>(
      const typename TComplex::Clique &)
    > Select ,
    std::function<
    bool(
      const TComplex & ,
      const typename TComplex::Cell & )
    > Skel,
    bool verbose )
{
  if(verbose) trace.beginBlock("Asymetric Thinning Scheme");
 
  using Cell = typename TComplex::Cell;
  using Data = typename TComplex::Data;
 
  // Initialize K set, and VoxelComplex containers;
  struct FirstComparator
  {
    bool operator() (const std::pair<Cell, Data>& lhs, const std::pair<Cell, Data>& rhs) const
    {
        return lhs.first < rhs.first;
    }
  };
  std::set<std::pair<Cell, Data>, FirstComparator> selected_voxels_per_dim;
  TComplex already_selected_voxels(vc.space());
  TComplex constraint_set(vc.space());
 
  auto & Z = selected_voxels_per_dim;
  auto & Y = already_selected_voxels;
  auto & K = constraint_set;
  auto X = vc;
  Y.copySimplicityTable(vc);
  K.copySimplicityTable(vc);
  typename TComplex::Parent x_y(vc.space());
  typename TComplex::Parent x_k(vc.space());
 
  bool stability{false};
  uint64_t generation{0};
  auto xsize =  X.nbCells(3);
  auto xsize_old =  xsize;
  typename TComplex::CliqueContainer critical_cliques;
 
  if(verbose){
      trace.info() << "generation: " << generation <<
        " ; X.nbCells(3): " << xsize << std::endl;
  }
  do {
    ++generation;
 
    Y = K ;
    /* Select voxels from critical d-cliques,
     * having priority d=3 down to 0 using Select(clique)
     * which returns the selected voxel. */
    // X.close(); // Close K-Space from current voxels.
    for (int d = 3 ; d >= 0 ; --d) {
      Z.clear();
      /* Y is the already_selected voxels (in this generation)
       * Search critical cliques only in the cells of X-Y,
       * but which are critical for X. */
      x_y = X - Y;
      // x_y.close();
      critical_cliques = X.criticalCliquesForD(d, x_y);
      for(auto & clique : critical_cliques){
        // if (d!=3)
        //   if (! (clique  <= (x_y)) )  continue ;
        auto yinclude = std::find_if(
            std::begin(clique),
            std::end(clique),
            [&Y]( const Cell& c) {
              return Y.belongs(c);
            });
        if (yinclude != std::end(clique)) continue ;
        Z.insert(Select(clique));
      } // critical_cliques
      // Y = Y union Z. From set to VoxelComplex
      for( const auto & selected_celldata_pair : Z)
        Y.insertVoxelCell(selected_celldata_pair);
    } // d-loop
    X = Y;
    // X - K is equal to X-Y, which is equal to a Ynew - Yold
    x_k = X  - K;
    for (auto it = x_k.begin(3), itE = x_k.end(3) ; it != itE ; ++it ){
      auto new_voxel = it->first ;
      if( Skel(X, new_voxel) == true){
        K.insertVoxelCell(new_voxel);
        // K.insertCell(3, new_voxel);
        // K.objectSet().insert(K.space().uCoords(new_voxel));
      }
    }
 
    // Stability Update:
    xsize = X.nbCells(3);
    if(xsize == xsize_old) stability = true;
    xsize_old = xsize;
 
    if(verbose){
      trace.info() << "generation: " << generation <<
        " ; X.nbCells(3): " << xsize <<
        " ; K (constrain set): " << K.nbCells(3) <<
        " ; Y (selected in this generation): " << Y.nbCells(3) <<
        " ; X - K ie. Ynew - Yold : " << x_k.nbCells(3) <<
        " ; X - Y  : " << x_y.nbCells(3) << std::endl;
    }
  } while( !stability );
 
  if(verbose)
      trace.endBlock();
 
  return X;
}
 
 
template < typename TComplex >
TComplex
DGtal::functions::
persistenceAsymetricThinningScheme(
    TComplex & vc ,
    std::function<
    std::pair<typename TComplex::Cell, typename TComplex::Data> (
      const typename TComplex::Clique &)
    > Select ,
    std::function<
    bool(
      const TComplex & ,
      const typename TComplex::Cell & )
    > Skel,
    uint32_t persistence,
    bool verbose )
{
  if(verbose) trace.beginBlock("Persistence asymetricThinningScheme");
 
  using Cell = typename TComplex::Cell;
  using Data = typename TComplex::Data;
 
  // Initialize Z set, and VoxelComplex containers;
  struct FirstComparator
  {
    bool operator() (const std::pair<Cell, Data>& lhs, const std::pair<Cell, Data>& rhs) const
    {
        return lhs.first < rhs.first;
    }
  };
  TComplex already_selected_voxels(vc.space());
  TComplex constraint_set(vc.space());
 
  auto & Y = already_selected_voxels;
  auto & K = constraint_set;
  TComplex X = vc;
  TComplex x_y(vc.space());
  Y.copySimplicityTable(vc);
  K.copySimplicityTable(vc);
  x_y.copySimplicityTable(vc);
 
  bool stability{false};
  unsigned int generation{0};
  auto xsize =  X.nbCells(3);
  auto xsize_old =  xsize;
  const bool close_it = true;
  typename TComplex::CliqueContainer critical_cliques;
 
  if(verbose){
      trace.info() << "Initial Voxels at generation: " << generation <<
        " ; X.nbCells(3): " << xsize << std::endl;
  }
    // cubical cell data (aka, the birth_date) is initialized to zero already.
  do {
    ++generation;
    // Update birth_date for our Skel function. (isIsthmus for example)
    for (auto it = X.begin(3), itE = X.end(3) ; it != itE ; ++it ){
      // Ignore voxels existing in K set.(ie: X-K)
      if (K.findCell(3, it->first) != K.end(3))
        continue;
      auto & voxel = it->first ;
      auto & ccdata = it->second.data ;
      if( Skel(X, voxel) == true && ccdata == 0)
        ccdata = generation;
    }
    Y = K ;
    x_y = X; //optimization instead of x_y = X-Y, use x_y -= Y;
    // d-cliques: From voxels (d=3) to pointels (d=0)
    for (int d = 3 ; d >= 0 ; --d) {
      // Search only critical cliques of X in the set X \ Y
      x_y -= Y; // Y is growing in this d-loop.
      critical_cliques = X.criticalCliquesForD(d, x_y);
 
      for(auto & clique : critical_cliques)
      {
        // Find the selected voxel in x to get associated Data.
        auto selected_voxelpair = Select(clique);
        selected_voxelpair.second = X.findCell(3,selected_voxelpair.first)->second;
        Y.insertVoxelCell(selected_voxelpair, close_it);
      }
    } // d-loop
    // Y is now  populated with some selected critical cliques for this generation,
    // plus with the fixed points in K.
    // Store result (final or for the next generation)
    X = Y;
 
    /*** Modify K (the set of voxels to conserve) ***/
    // Insert in K (inmutable set) only if the NEW critical cells of this generation (Y-K) are:
    // a) Part of the Skeleton (Skel) user wants to preserve.
    // b) persistent enough.
 
    // Update K
    Y -= K;
    for (auto it = Y.begin(3), itE = Y.end(3) ; it != itE ; ++it ){
        auto & voxel = it->first;
        auto & ccdata = it->second.data;
        bool is_skel = Skel(X, voxel);
        bool is_persistent_enough = (generation + 1 - ccdata) >= persistence;
        if (is_skel && is_persistent_enough)
          K.insertVoxelCell(voxel, close_it, ccdata);
    }
 
    if(verbose){
      trace.info() << "generation: " << generation <<
        " ; X  at start generation: " << xsize_old <<
        " ; K (constraint set): " << K.nbCells(3) <<
        " ; Y - K: " << Y.nbCells(3) <<
        " ; Y (X at the end): " << X.nbCells(3) <<
        std::endl;
    }
 
    /***** Stability Update: ****/
    xsize = X.nbCells(3);
    if(xsize == xsize_old)
      stability = true;
    xsize_old = xsize;
 
  } while( !stability );
 
  if(verbose) trace.endBlock();
 
  return X;
}
 
//////////////////////////////////////////////////////////////////////////////
// Select Functions
//////////////////////////////////////////////////////////////////////////////
template < typename TComplex >
std::pair<typename TComplex::Cell, typename TComplex::Data>
DGtal::functions::selectFirst(
      const typename TComplex::Clique & clique)
{
  return *clique.begin(3);
}
 
template < typename TComplex >
std::pair<typename TComplex::Cell, typename TComplex::Data>
DGtal::functions::selectRandom(
      const typename TComplex::Clique & clique)
{
  static std::mt19937 gen( std::random_device{}() );
  return selectRandom<TComplex>(clique, gen);
}
 
template < typename TComplex, typename TRandomGenerator >
std::pair<typename TComplex::Cell, typename TComplex::Data>
DGtal::functions::selectRandom(
      const typename TComplex::Clique & clique,
      TRandomGenerator & gen
)
{
  auto size = clique.nbCells(3);
  std::uniform_int_distribution<> dis(0, size - 1);
  auto it = clique.begin(3);
  std::advance(it, dis(gen));
  return *it;
}
 
template < typename TDistanceTransform, typename TComplex >
std::pair<typename TComplex::Cell, typename TComplex::Data>
DGtal::functions::selectMaxValue(
    const TDistanceTransform & dist_map,
    const typename TComplex::Clique & clique)
{
  typename TDistanceTransform::Value max_value{0};
  typename TDistanceTransform::Value value{0};
  using ComplexConstIterator = typename TComplex::CellMapConstIterator;
  ComplexConstIterator selected_pair;
 
  for(ComplexConstIterator it = clique.begin(3), itE = clique.end(3);
      it != itE ; ++it){
    value = dist_map(clique.space().uCoords(it->first));
    // trace.info() << "P: " << it->first << " V: " << value << std::endl;
    if(value > max_value){
      max_value = value;
      selected_pair = it;
      continue;
    }
    if(value == max_value){
      // max_value = value;
      // selected_pair = it;
      continue; // TODO choose one wisely when they have same DM value.
    }
  }
 
  return *selected_pair;
}
//////////////////////////////////////////////////////////////////////////////
// Skeleton Functions
//////////////////////////////////////////////////////////////////////////////
 
template < typename TComplex >
bool
DGtal::functions::skelUltimate( const TComplex &,
                                const typename TComplex::Cell &)
{
  return false;
}
/*--------------------------------------------------------------------------*/
 
template < typename TComplex >
bool
DGtal::functions::skelEnd( const TComplex & vc,
                           const typename TComplex::Cell & cell)
{
  const auto pnsize = vc.properNeighborhoodVoxels(cell).size();
  return (pnsize == 1);
}
 
template < typename TComplex >
bool
DGtal::functions::skelSimple( const TComplex & vc,
                              const typename TComplex::Cell & cell)
{
  return vc.isSimple(cell);
}
 
template < typename TComplex >
bool
DGtal::functions::skelIsthmus(const TComplex & vc,
    const typename TComplex::Cell & cell)
{
  if(twoIsthmus<TComplex>(vc,cell))
    return true;
  else if(oneIsthmus<TComplex>(vc,cell))
    return true;
  else
    return false;
}
 
template < typename TComplex >
bool
DGtal::functions::oneIsthmus(const TComplex & vc,
                             const typename TComplex::Cell & cell)
{
 
  auto &ks = vc.space();
  auto point_cell = ks.uCoords(cell);
  // Create an object from the set.
  using DigitalTopology = DGtal::Z3i::DT26_6;
  using DigitalSet = DGtal::DigitalSetByAssociativeContainer<
    DGtal::Z3i::Domain,
    std::unordered_set<typename DGtal::Z3i::Domain::Point>>;
  using Object = DGtal::Object<DigitalTopology, DigitalSet>;
  DigitalTopology::ForegroundAdjacency adjF;
  DigitalTopology::BackgroundAdjacency adjB;
  auto domain = DGtal::Z3i::Domain(ks.lowerBound(), ks.upperBound());
  DigitalTopology topo(
      adjF, adjB, DGtal::DigitalTopologyProperties::JORDAN_DT);
  Object object(topo, domain);
  for(auto & properNeighborhoodVoxel : vc.properNeighborhoodVoxels(cell))
  {
    object.pointSet().insertNew(ks.uCoords(properNeighborhoodVoxel));
  }
 
  // auto spN = vc.object().properNeighborhood(point_cell);
  if (isZeroSurface(object)) return true;
 
  // else thinning
  TComplex pre_thin(ks);
  pre_thin.construct(object.pointSet());
 
  auto after_thin = asymetricThinningScheme<TComplex>( pre_thin,
      selectFirst<TComplex>,
      skelUltimate<TComplex>);
 
  object.pointSet().clear();
  for (auto it = after_thin.begin(3), itE = after_thin.end(3) ; it != itE ; ++it ){
    object.pointSet().insertNew(ks.uCoords(it->first));
  }
 
  return isZeroSurface(object);
}
 
template < typename TComplex >
bool
DGtal::functions::twoIsthmus( const TComplex & vc,
                              const typename TComplex::Cell & cell)
{
  auto &ks = vc.space();
  auto point_cell = ks.uCoords(cell);
  // Create an object from the set.
  using DigitalTopology = DGtal::Z3i::DT26_6;
  using DigitalSet = DGtal::DigitalSetByAssociativeContainer<
    DGtal::Z3i::Domain,
    std::unordered_set<typename DGtal::Z3i::Domain::Point>>;
  using Object = DGtal::Object<DigitalTopology, DigitalSet>;
  DigitalTopology::ForegroundAdjacency adjF;
  DigitalTopology::BackgroundAdjacency adjB;
  auto domain = DGtal::Z3i::Domain(ks.lowerBound(), ks.upperBound());
  DigitalTopology topo(
      adjF, adjB, DGtal::DigitalTopologyProperties::JORDAN_DT);
  Object object(topo, domain);
  for(auto & properNeighborhoodVoxel : vc.properNeighborhoodVoxels(cell))
  {
    object.pointSet().insertNew(ks.uCoords(properNeighborhoodVoxel));
  }
  if (isOneSurface(object)) return true;
 
  // else thinning
  TComplex pre_thin(ks);
  pre_thin.construct(object.pointSet());
 
  auto after_thin = asymetricThinningScheme<TComplex>( pre_thin,
      selectFirst<TComplex>,
      skelUltimate<TComplex>);
 
  object.pointSet().clear();
  for (auto it = after_thin.begin(3), itE = after_thin.end(3) ; it != itE ; ++it ){
    object.pointSet().insertNew(ks.uCoords(it->first));
  }
 
  return isOneSurface(object);
}
 
template < typename TComplex >
bool
DGtal::functions::skelWithTable(
    const boost::dynamic_bitset<> & table,
    const std::unordered_map<typename TComplex::Point, unsigned int> & pointToMaskMap,
    const TComplex & vc,
    const typename TComplex::Cell & cell)
{
  auto conf = functions::getSpelNeighborhoodConfigurationOccupancy(
      vc,
      vc.space().uCoords(cell),
      pointToMaskMap);
  return table[conf];
}
///////////////////////////////////////////////////////////////////////////////
// Object Helpers
template < typename TObject >
bool
DGtal::functions::isZeroSurface(const TObject & small_obj)
{
  if (small_obj.size() != 2) return false;
  auto connectedness = small_obj.computeConnectedness();
  if( connectedness == DGtal::Connectedness::CONNECTED ) return false;
 
  return true;
}
 
template < typename TObject >
bool
DGtal::functions::isOneSurface(const TObject & small_obj)
{
    auto connectedness = small_obj.computeConnectedness();
    if( connectedness == DGtal::Connectedness::DISCONNECTED )
      return false;
 
    for (const auto & p : small_obj.pointSet()){
      auto pN = small_obj.properNeighborhood(p);
      if(!isZeroSurface(pN) ) return false;
    }
 
    return true;
}
 
template <typename TObject >
std::vector< TObject >
DGtal::functions::
connectedComponents(const TObject & input_obj, bool verbose)
{
 
  // Graph related alias
    using Graph = TObject ;
    using vertex_descriptor =
      typename boost::graph_traits<Graph>::vertex_descriptor ; // Object::Vertex
    using edge_descriptor =
      typename boost::graph_traits<Graph>::edge_descriptor ; // Object::Edge
    using vertices_size_type =
      typename boost::graph_traits<Graph>::vertices_size_type ; // Object::Size
 
  using StdColorMap = std::map< vertex_descriptor, boost::default_color_type > ;
  StdColorMap colorMap;
  boost::associative_property_map< StdColorMap > propColorMap( colorMap );
 
  using StdComponentMap = std::map< vertex_descriptor, vertices_size_type > ;
  StdComponentMap componentMap;
  boost::associative_property_map< StdComponentMap > propComponentMap( componentMap );
 
  if(verbose) trace.beginBlock( "Boost connected_components");
  vertices_size_type nbComp =
    boost::connected_components // boost graph connected components algorithm.
    ( input_obj, // the graph
      propComponentMap, // the mapping vertex -> label
      boost::color_map( propColorMap ) // this map is used internally when computing connected components.
    );
  if(verbose) trace.info() << "num_components = " << nbComp << std::endl;
 
  if(verbose) trace.beginBlock( "Filter graph and get components");
 
  using ComponentGraph =
    boost::filtered_graph<
    Graph,
    std::function<bool(edge_descriptor)>,
    std::function<bool(vertex_descriptor)>
      >;
  auto &g = input_obj;
 
  std::vector<ComponentGraph> component_graphs;
 
  for (size_t i = 0; i < nbComp; i++)
    component_graphs.emplace_back(g,
        [componentMap,i,&g](edge_descriptor e) {
        return componentMap.at(boost::source(e,g) )==i
        || componentMap.at(boost::target(e,g))==i;
        },
        [componentMap,i](vertex_descriptor v) {
        return componentMap.at(v)==i;
        });
  if(verbose) trace.endBlock();
 
  std::vector<TObject> obj_components;
  for (auto && c : component_graphs){
    // Create empty but valid object.
    obj_components.emplace_back(
        TObject( input_obj.topology(),input_obj.domain() ));
    for (auto && vp = boost::vertices(c); vp.first != vp.second ; ++vp.first){
      obj_components.back().pointSet().insertNew(*vp.first);
    }
  }
  if(verbose) trace.endBlock();
 
  return obj_components;
}
///////////////////////////////////////////////////////////////////////////////
 
///////////////////////////////////////////////////////////////////////////////