VoronoiMapComplete.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
 *
 * @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 ( int i = 0; i < static_cast<int>(subRangePoints.size()); ++i )  //MSVC requires signed type for openmp
    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;
}