DistanceTransformation.ih

#include <cstdlib>
#include <boost/lexical_cast.hpp>

// IMPLEMENTATION of inline methods.

// ----------------------- Standard services ------------------------------

template <typename S,typename P,DGtal::uint32_t p, typename IntLong>
inline
DGtal::DistanceTransformation<S,P, p, IntLong>::DistanceTransformation(const Domain &aDomain,
                                                                       const PointPredicate & aPredicate):
  myDomain(aDomain), myPointPredicate(aPredicate)
{
}


template <typename S, typename P,DGtal::uint32_t p, typename IntLong>
inline
DGtal::DistanceTransformation<S, P,p, IntLong>::~DistanceTransformation()
{
}




template <typename S, typename P,DGtal::uint32_t p, typename IntLong>
inline
bool
DGtal::DistanceTransformation<S,P,p, IntLong>::checkTypesValidity (  ) const
{
  typename Point::UnsignedComponent maxExtent = ( myDomain.upperBound() - 
                                                  myDomain.lowerBound() ).normInfinity();
  
  //Estimate worst-case bit size (special case for p=0 == Linfinity
  double bitSize;
  
  if (SeparableMetric::p != 0)
    bitSize = SeparableMetric::p * (log ((double) S::dimension * 2 * maxExtent ) / log(2.0));
  else
    bitSize =  sizeof ( typename OutputImage::Value ) * 8;
  
  if ( bitSize > sizeof ( typename OutputImage::Value ) * 8 )
    {
      trace.warning() << "(DistanceTransformation) The output image value range may"
                      <<" not be sufficient to store the exact values  according to"
                      <<" the input domain extent." << endl;
      return false;
    }
  return true;
}


template <typename S, typename P, DGtal::uint32_t p, typename IntLong>
inline
typename DGtal::DistanceTransformation<S,P, p, IntLong>::OutputImage
DGtal::DistanceTransformation<S,P, p, IntLong>::compute ( )
{

  //We trace type validdity check result;
  checkTypesValidity ( );

  //We copy the image extent and translate the image domains to (0,..0)x(Upper-Lower)
  myLowerBoundCopy = Point(); //(O,O,...O)
  myUpperBoundCopy = myDomain.upperBound() - myDomain.lowerBound();
  myDisplacementVector = myDomain.lowerBound();

  myExtent = myUpperBoundCopy - myLowerBoundCopy;
  myInfinity  = myMetric.power(static_cast<typename S::Integer>(S::dimension) * 
                               myExtent.normInfinity() + 1);
  
  Domain workingDomain(myLowerBoundCopy, myUpperBoundCopy);
  OutputImage output ( workingDomain );
  OutputImage swap ( workingDomain );
  bool isSwap = true;

  //First step
  computeFirstStep ( output );

  //We process the dimensions swaping the temporary buffers
  for ( Dimension dim = 1; dim < S::dimension ; dim++ )
    {
      if ( isSwap )
        computeOtherSteps ( output, swap, dim );
      else
        computeOtherSteps ( swap, output, dim );

      isSwap = !isSwap;
    }

  //We translate the output image to the correct position and return.
  if ( !isSwap )
    {
      swap.translateDomain(myDisplacementVector);
      return swap;
    }

  else
    {
      output.translateDomain(myDisplacementVector);
      return output;
    }
}


template <typename S, typename P,DGtal::uint32_t p, typename IntLong>
inline
void
DGtal::DistanceTransformation<S, P, p, IntLong>::computeFirstStep ( OutputImage &output ) const
{
  trace.beginBlock ( "DT dimension 0" );
  Point startingPoint = myLowerBoundCopy;

  typedef typename Domain::ConstSubRange::ConstIterator ConstDomIt;
  
  //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<Size> subdomain;
  subdomain.reserve(S::dimension - 1);
  for (Dimension k = S::dimension-1; k>0 ; --k)
    subdomain.push_back( k );

  Domain localDomain(myLowerBoundCopy, myUpperBoundCopy);

  //We process the dimensions to construct a Point
  for (ConstDomIt it = localDomain.subRange( subdomain ).begin(),
         itend = localDomain.subRange( subdomain ).end(); it != itend; ++it)
    {
      // trace.info()  << "Processing 1D slice starting at =" << (*it) << endl;
      computeFirstStep1D (output, (*it) );
    }

  trace.endBlock();
}


template <typename S, typename P,DGtal::uint32_t p, typename IntLong>
inline
void
DGtal::DistanceTransformation<S,P, p, IntLong>::computeOtherSteps ( const OutputImage &input, 
                                                                    OutputImage &output,
                                                                    const Dimension dim ) const
{
  std::string title = "DT dimension " +  boost::lexical_cast<string>( dim ) ;
  trace.beginBlock ( title );

  typedef typename Domain::ConstSubRange::ConstIterator ConstDomIt;
  
  //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<Size> subdomain;
  subdomain.reserve(S::dimension - 1);
  for (unsigned int k = 0; k < S::dimension ; k++)
    if ( (S::dimension - 1 - k) != dim)
      subdomain.push_back( S::dimension - 1 - k );

  Domain localDomain(myLowerBoundCopy, myUpperBoundCopy);
  Size maxSize = myExtent.normInfinity();

  //Stacks used in the envelope computation
  Abscissa *s = new Abscissa[maxSize+1];
  Abscissa *t = new Abscissa[maxSize+1];
  
  ASSERT( s != NULL);
  ASSERT( t != NULL);
  
  //We process the dimensions to construct a Point
  for (ConstDomIt it = localDomain.subRange( subdomain ).begin(),
         itend = localDomain.subRange( subdomain ).end();
       it != itend; ++it)
    {
      computeOtherStep1D ( input, output, (*it), dim, s, t );
    }
  
  delete[] s;
  delete[] t;
  trace.endBlock();

}

template <typename S, typename P,DGtal::uint32_t p, typename IntLong>
void
DGtal::DistanceTransformation<S,P,p, IntLong>::computeFirstStep1D ( OutputImage & output,
                                                                    const Point &startingPoint ) const
{
  Point prec = startingPoint;
  Point point = startingPoint;

  //PRECOND : output can store 2*(myUpperBoundCopy[0] -  myLowerBoundCopy[0]) in its valuetype
  //INFTY = something > myUpperBoundCopy[0] -  myLowerBoundCopy[0]
  if ( myPointPredicate ( point + myDisplacementVector ) )
    output.setValue ( point, myInfinity );
  else
    output.setValue ( point, 0 );

  //Forward scan 
  for ( point[0] = 1; point[0] <= myUpperBoundCopy[0]; point[0]++ )
    {
      // trace.warning() << point << " ";
      if ( myPointPredicate ( point  + myDisplacementVector ))
        output.setValue ( point, 1 + output ( prec ) );
      else
        output.setValue ( point, 0 );

      prec[0] = point[0];
    }


  //prec is the the rightmost point of "point'
  prec[0] =  myUpperBoundCopy[0];

  //Backward scan
  for ( point[0] = myUpperBoundCopy[0] - 1; point[0] >= 0 ; point[0]-- )
    {
      if ( output ( prec ) < output ( point ) )
        output.setValue ( point, 1 + output ( prec ) );
      prec[0] = point[0];
    }

  //final computation
  for ( point[0] = 0; point[0] <= myUpperBoundCopy[0]; point[0]++ )
    if (output( point ) < myInfinity)
      output.setValue ( point, myMetric.power( (const int)output ( point ) ));
    else
      output.setValue ( point, myInfinity);
}


template <typename S,typename P, DGtal::uint32_t p, typename IntLong>
void
DGtal::DistanceTransformation<S,P,p, IntLong>::computeOtherStep1D ( const OutputImage & input,
                                                                    OutputImage & output,
                                                                    const Point &startingPoint,
                                                                    const Size dim,
                                                                    Abscissa s[],
                                                                    Abscissa t[] ) const
{
  Abscissa w;
  Abscissa q = 0;  //index for the stack "head"
  
  Point sQ = startingPoint;
  Point pU = startingPoint;

  // We look for the first point in the 1D column with distance different from
  // myInfinity
  pU[dim] = myLowerBoundCopy[dim];
  while ((pU[dim] <= myUpperBoundCopy[dim]) && (input ( pU ) == myInfinity))
    pU[dim] ++;


  // All points are set to +infinity, we just copy the infinity value and return
  if ( pU[dim] > myUpperBoundCopy[dim] )
    {
      for(typename OutputImage::SpanIterator it= output.spanBegin(startingPoint,dim), 
            itend=output.spanEnd(startingPoint,dim); 
          it != itend; ++it)
        output.setValue(it, myInfinity);
      return;
    }
  
  //Stack structure
  q = 0;
  s[q] = pU[dim]; // first point with DT!=infinity
  sQ[dim] = s[q]; // nD Point corresponding to the head of the stack
  t[q] = myLowerBoundCopy[dim]; 

  //Forward Scan 
  //We scan all the pixels
  for ( Abscissa u = pU[dim] + 1; u <= myUpperBoundCopy[dim] ; u++ )
    {
      pU[ dim ] = u;
      if ( input( pU ) == myInfinity )
        {
          continue;
        }
      
      while ( ( q >= 0 ) &&
              ( myMetric.F ( t[q], s[q], input ( sQ ) ) >
                myMetric.F ( t[q], u, input ( pU ) ) ) )
        {
          q--;
          if (q>=0)
            sQ[dim] = s[q];
        }
      
      if ( q < 0 )
        {
          q = 0;
          s[0] = u;
          sQ[dim] = u;
          t[0] = myLowerBoundCopy[dim];
        }
      else
        {
          sQ[dim] = s[q];
          w = 1 + myMetric.Sep ( s[q],
                                 input ( sQ ),
                                 u,
                                 input ( pU ) );
    
          if (( w <= myUpperBoundCopy[dim] ) 
              && (w>= myLowerBoundCopy[dim]))
            {
              q++;
              s[q] = u;
              sQ[dim] = u;
              t[q] = w;
            }
        }
    }
  
  Point last = startingPoint;
  
  ASSERT(q>=0);
  
  sQ[dim] = s[q];
  
  //Backward Scan
  for (last[dim] = myUpperBoundCopy[dim];
       last[dim] >= myLowerBoundCopy[dim] ;
       last[dim]-- )
    {
      output.setValue ( last, myMetric.F ( last[dim] , s[q], input ( sQ ) ) );
      if (( last[dim] == t[q] ) && (q > 0))
        {
          q--;
          sQ[dim] = s[q];
        }
    }
}


//                                                                           //