ReverseDistanceTransformation.ih

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

// IMPLEMENTATION of inline methods.

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

template <typename I, DGtal::uint32_t p, typename IntShort>
inline
DGtal::ReverseDistanceTransformation<I, p, IntShort>::ReverseDistanceTransformation(const IntShort defaultForeground, 
                        const IntShort defaultBackground):
  myForegroundValue(defaultForeground),
  myBackgroundValue(defaultBackground)
{
}

template <typename I, DGtal::uint32_t p, typename IntShort>
inline
DGtal::ReverseDistanceTransformation<I, p, IntShort>::~ReverseDistanceTransformation()
{
}


template <typename I, DGtal::uint32_t p, typename IntShort>
inline
typename DGtal::ReverseDistanceTransformation<I, p, IntShort>::OutputImage
DGtal::ReverseDistanceTransformation<I, p, IntShort>::reconstruction(const I & aImage)
{
  //We copy the image extent.
  myLowerBoundCopy = Point(); //(O,O,...O)
  myUpperBoundCopy = aImage.domain().upperBound() - aImage.domain().lowerBound();
  myDisplacementVector = aImage.domain().lowerBound();
  myExtent = myUpperBoundCopy - myLowerBoundCopy;
  
  I output ( typename I::Domain(myLowerBoundCopy, myUpperBoundCopy));
  I swap ( typename I::Domain(myLowerBoundCopy, myUpperBoundCopy));

  bool isSwap = reconstructionInternal(aImage,output,swap);

  if ( !isSwap )
    {
      //swap.translateDomain(myDisplacementVector);
      return castValues(swap);
    }
  else
    {
      //output.translateDomain(myDisplacementVector);
      return castValues(output);
    }
}

template <typename I, DGtal::uint32_t p, typename IntShort>
template <typename Set>
inline
void
DGtal::ReverseDistanceTransformation<I, p, IntShort>::reconstructionAsSet(Set &aSet,const I & aImage)
{
  //BOOST_CONCEPT_ASSERT(( CDigitalSet<Set> ));
  
  //We copy the image extent and translate the image domains to (0,..0)x(Upper-Lower)
  myLowerBoundCopy = Point(); //(O,O,...O)
  myUpperBoundCopy = aImage.domain().upperBound() - aImage.domain().lowerBound();
  myDisplacementVector = aImage.domain().lowerBound();
  myExtent = myUpperBoundCopy - myLowerBoundCopy;
  
  typename I::Domain domain(myLowerBoundCopy,myUpperBoundCopy);

  I output ( typename I::Domain(myLowerBoundCopy, myUpperBoundCopy ));
  I swap ( typename I::Domain( myLowerBoundCopy, myUpperBoundCopy ));

  bool isSwap = reconstructionInternal(aImage,output,swap);
  
  for(typename I::Domain::ConstIterator it=domain.begin(),
  itend = domain.end();
      it != itend;
      ++it)
    if (isSwap)
      {
  if ( output(*it)  > 0)
    aSet.insertNew((*it) + myDisplacementVector);
      }
    else
      {
  if ( swap(*it)  > 0)
    aSet.insertNew((*it)  + myDisplacementVector);
      }
}


template <typename I, DGtal::uint32_t p, typename IntShort>
inline
bool
DGtal::ReverseDistanceTransformation<I, p, IntShort>::reconstructionInternal(const I & aImage, I &output, I &swap)
{
 
  //We copy input image (implicit translation of values)
  typename I::Iterator ito=output.begin();
  
  for(typename I::ConstIterator it=aImage.begin(), itend = aImage.end();
      it != itend;
      ++it,++ito)
    (*ito) = (*it);

  bool isSwap = true;

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

      isSwap = !isSwap;
    }

  return isSwap;
}




template <typename I, DGtal::uint32_t p, typename IntShort>
inline
typename DGtal::ReverseDistanceTransformation<I, p, IntShort>::OutputImage
DGtal::ReverseDistanceTransformation<I, p, IntShort>::castValues(const I &input) const
{
  //We threshold input values
  OutputImage output(typename OutputImage::Domain(myLowerBoundCopy, myUpperBoundCopy));
  typename OutputImage::Iterator ito = output.begin();
  for(typename I::ConstIterator it=input.begin(),
  itend = input.end();
      it != itend;
      ++it,++ito)
    if ( (*it)  > 0)
      (*ito) =  myForegroundValue;
    else
      (*ito) =  myBackgroundValue;
  return output;
}



template <typename I, DGtal::uint32_t p, typename IntShort>
inline
void
DGtal::ReverseDistanceTransformation<I, p, IntShort>::computeSteps ( const I & input, 
                     I &output,
                     const Dimension dim) const
{
  std::string title = "RDT 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(I::dimension - 1);
  for (unsigned int k = 0; k < I::dimension ; k++)
    if ( (I::dimension - 1 - k) != dim)
      subdomain.push_back( I::dimension - 1 - k );

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

  //Stacks used in the envelope computation
  Integer *s = new Integer[maxSize+1];
  Integer *t = new Integer[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)
    {
      computeSteps1D ( input, output, (*it), dim, s, t );
    }

  delete[] s;
  delete[] t;

  trace.endBlock();

}

template <typename I, DGtal::uint32_t p, typename IntShort>
void
DGtal::ReverseDistanceTransformation<I, p, IntShort>::computeSteps1D ( const I & input,
                       I & output,
                       const Point &startingPoint,
                       const Size dim,
                       Integer s[],
                       Integer t[] ) const
{
  Coordinate w;
  Coordinate q = 0;

  Point sQ = startingPoint;
  Point pU = startingPoint;

  //trace.info()<<endl<<"Enter startingPoint="<<startingPoint<<endl;

  //init of the stack structure
  pU[dim] = myLowerBoundCopy[dim];
  while ((pU[dim] <= myUpperBoundCopy[dim]) 
   && (input ( pU ) == 0))
    pU[dim] ++;

  if ( pU[dim] > myUpperBoundCopy[dim] )
    {
      for(typename I::SpanIterator it= output.spanBegin(startingPoint,dim),
      itend=output.spanEnd(startingPoint,dim); 
    it != itend; ++it)
  output.setValue(it, 0);
      return;
    }
  
  q = 0;
  s[q] = pU[dim]; 
  sQ[dim] = s[q];
  t[q] = myLowerBoundCopy[dim];

  //Forward Scan

  for ( int u = pU[dim] + 1; u <= myUpperBoundCopy[dim] ; u++ )
    {
      pU[ dim ] = u;
      if ( input( pU ) == 0 )
  {
    //trace.info() << "Ctd " << pU ;
    continue;
  }
    
      while ( ( q >= 0 ) &&
        ( myMetric.reversedF ( t[q], s[q], (const int)input ( sQ ) ) <
    myMetric.reversedF ( t[q], u, (const int)input ( pU ) ) ) )
  {
    q--;
    if (q>=0)
      sQ[dim] = s[q];
  }

      if ( q < 0 )
  {
    q = 0;
    s[0] = u;
    t[0] = myLowerBoundCopy[dim];
    sQ[dim] = u;
  }
      else
  {
    sQ[dim] = s[q];
    w = 1 + myMetric.reversedSep ( s[q],
           (const int)input ( sQ ),
           u,
           (const int)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);

  //Backward Scan
  for (last[dim] = myUpperBoundCopy[dim];
       last[dim] >= myLowerBoundCopy[dim] ;
       last[dim]-- )
    {    
      if (myMetric.reversedF ( last[dim] , s[q], (const int)input ( sQ )) > 0 )
      output.setValue ( last, myMetric.reversedF ( last[dim] , s[q], (const int)input ( sQ ) ) );
      else
  output.setValue ( last, 0 );
    
      if (( last[dim] == t[q] ) && (q > 0))
  {
    q--;
    sQ[dim] = s[q];
  }
    }
}


//                                                                           //