ArithmeticalDSS.ih

// IMPLEMENTATION of inline methods.

#include <cstdlib>
#include <boost/math/common_factor_rt.hpp>





// Implementation of inline methods                                          //


template <typename TIterator, typename TInteger, int connectivity>
inline
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::ArithmeticalDSS()
{
  myUf = Point(0,0);
  myLf = Point(0,0);
  myUl = Point(0,0);
  myLl = Point(0,0);

  myA = 0; 
  myB = 0; 
  myMu = 0;
  myOmega = 0;
  myNbUpPat = 0;
  myNbLowPat = 0;
}


template <typename TIterator, typename TInteger, int connectivity>
inline
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::ArithmeticalDSS(const ConstIterator& it)
{
  init(it);
}

template <typename TIterator, typename TInteger, int connectivity>
inline
void DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::init(const ConstIterator& it)
{

  myF = it;
  myL = it;

  myUf = *it;
  myLf = *it;
  myUl = *it;
  myLl = *it;

  myA = 0; 
  myB = 0; 
  myMu = 0;
  myOmega = 0;
  myNbUpPat = 0;
  myNbLowPat = 0;

  mySteps.clear();

}

template <typename TIterator, typename TInteger, int connectivity>
inline
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::ArithmeticalDSS ( 
                                                                          const ArithmeticalDSS<TIterator,TInteger,connectivity> & other ) : 
  myA(other.myA), myB(other.myB), myMu(other.myMu), myOmega(other.myOmega),
  myNbUpPat(other.myNbUpPat),myNbLowPat(other.myNbLowPat),
  myUf(other.myUf), myUl(other.myUl), myLf(other.myLf), myLl(other.myLl), 
  myF(other.myF), myL(other.myL) {

  mySteps.clear();
  mySteps = other.mySteps;
}

template <typename TIterator, typename TInteger, int connectivity>
inline
typename DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::Self&
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>
::operator=( const ArithmeticalDSS<TIterator,TInteger,connectivity> & other )
{
  if ( this != &other )
    {
      myA = other.myA;
      myB = other.myB;
      myMu = other.myMu;
      myOmega = other.myOmega;
      myNbUpPat = other.myNbUpPat;
      myNbLowPat = other.myNbLowPat;
      myUf = other.myUf;
      myLf = other.myLf;
      myUl = other.myUl;
      myLl = other.myLl;
      myF = other.myF;
      myL = other.myL;
      mySteps = other.mySteps;

    }
  return *this;
}

template <typename TIterator, typename TInteger, int connectivity>
inline
typename ArithmeticalDSS<TIterator,TInteger,connectivity>::Reverse
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>
::getReverse() const {
  return Reverse(); 
}

template <typename TIterator, typename TInteger, int connectivity>
inline
typename DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::Self
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>
::getSelf() const {
  return Self(); 
}



template <typename TIterator, typename TInteger, int connectivity>
inline
bool
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::operator==( 
                                                                    const ArithmeticalDSS<TIterator,TInteger,connectivity>& other ) const {
  return ( ( (myOmega == other.myOmega) &&
             (myUf == other.myUf) && 
             (myUl == other.myUl) &&  
             (myLf == other.myLf) &&  
             (myLl == other.myLl) && 
             (*myF == *other.myF) &&  
             (*myL == *other.myL) ) || 
           ( (myOmega == other.myOmega) &&
             (myUf == other.myLl) && 
             (myUl == other.myLf) &&  
             (myLf == other.myUl) &&  
             (myLl == other.myUf) && 
             (*myF == *other.myL) &&  
             (*myL == *other.myF) ) );
}

template <typename TIterator, typename TInteger, int connectivity>
inline
bool
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::operator!=( 
                                                                    const ArithmeticalDSS<TIterator,TInteger,connectivity> & other ) const {
  return (!(*this == other));
}


template <typename TIterator, typename TInteger, int connectivity>
inline
bool 
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>
::extendForward( const ConstIterator & it, 
                 ConstIterator & lastIt, 
                 const Vector & lastMove,
                 Point & Uf,  Point & Ul,
                 Point & Lf,  Point & Ll )
{

  Point lastPoint(*it);

  Integer deviation = Tools<Integer,connectivity>
    ::norm(lastMove[0],lastMove[1]);
  //if the two last points are not connected return false
  if (deviation > NumberTraits<Integer>::ONE) return false;
  //if the two last points are confounded return true
  else if (deviation == NumberTraits<Integer>::ZERO) 
    {
      //add the point to the DSS
      lastIt = it; 
      return true;
    } 
  //otherwise
  else if (mySteps.size()<2) {

    //initialisation stage

    if (mySteps.size()==0) {

      mySteps.push_back(lastMove);

      //update
      lastIt = it; 
      Ul = lastPoint; 
      Ll = lastPoint;
      myA += lastMove[1]; 
      myB += lastMove[0]; 
      myMu = getRemainder(myUl); 
      myOmega = Tools<Integer,connectivity>::norm(myA,myB);
      myNbUpPat = 1;
      myNbLowPat = 1;

      return true;

    } else { //mySteps.size() == 1

      if (lastMove == mySteps.at(0)) {

        //update
        lastIt = it; 
        Ul = lastPoint; 
        Ll = lastPoint;
        myNbUpPat++;
        myNbLowPat++;

        return true;

      } else {

        //difference of the two steps viewed as points
        Vector diff = ( lastMove-mySteps.at(0) );
        if (Tools<Integer,connectivity>::dualNorm(diff[0],diff[1]) == NumberTraits<Integer>::ONE) {
          //the two steps are compatible

          //update
          lastIt = it; 
          //remainder
          Integer r = getRemainder(lastPoint);
          if ( r < myMu) {  //on the left
            Ul = lastPoint; 
            Lf = Ll;
            myA = myNbUpPat*mySteps.at(0)[1] + lastMove[1];
            myB = myNbUpPat*mySteps.at(0)[0] + lastMove[0];
            myNbUpPat = 1;
            myNbLowPat = 0;
          } else {            //on the right
            Ll = lastPoint; 
            Uf = Ul;
            myA = myNbUpPat*mySteps.at(0)[1] + lastMove[1];
            myB = myNbUpPat*mySteps.at(0)[0] + lastMove[0];
            myNbUpPat = 0;
            myNbLowPat = 1;
          }
          myMu = getRemainder(myUl); 
          myOmega = Tools<Integer,connectivity>::norm(myA,myB);

          mySteps.push_back(lastMove);

          return true;

        } else {  //the two steps are not compatible
          return false;
        }
      }

    }

  } else {  //mySteps.size() == 2

    //main stage

    //if there are only two steps at most
    if (hasLessThanTwoSteps(lastMove)) {
    
      //remainder
      Integer r = getRemainder(lastPoint);
          
      if ( (r < myMu-1)||(r > myMu+myOmega) )
        return false; //strongly exterior
      else { 
        //add the point to the DSS
        lastIt = it; 
        //update
        //if weakly interior
        if (r == myMu) {
          Ul = lastPoint;  
          myNbUpPat++;
        }
        if (r == myMu+myOmega-1) {
          Ll = lastPoint;
          myNbLowPat++;
        }
        //if weakly exterior  
        if (r == myMu-1) {              //on the left
          Ul = lastPoint; 
          Lf = Ll;
          myA = myUl[1] - myUf[1]; 
          myB = myUl[0] - myUf[0]; 
          myMu = getRemainder(myUl);
          myOmega = Tools<Integer,connectivity>::norm(myA,myB);
          myNbUpPat = 1;
          myNbLowPat = 0;
        } else if (r == myMu+myOmega) { //on the right
          Ll = lastPoint; 
          Uf = Ul; 
          myA = myLl[1] - myLf[1]; 
          myB = myLl[0] - myLf[0]; 
          myMu = getRemainder(myUl);
          myOmega = Tools<Integer,connectivity>::norm(myA,myB);
          myNbUpPat = 0;
          myNbLowPat = 1;
        }
        return true;

      }

    } else { 
      //if there is more than two steps 
      return false;
    }

  } 

}

//Deprecated
//--------------------------------------------------------------------
template <typename TIterator, typename TInteger, int connectivity>
inline
bool 
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::extendForward(const ConstIterator& it)
{

  return extendForward(it, myL, (Point(*it) - Point(*myL)), myUf, myUl, myLf, myLl); 

}

//Deprecated
//--------------------------------------------------------------------
template <typename TIterator, typename TInteger, int connectivity>
inline
bool 
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::extendBackward(const ConstIterator& it)
{

  return extendForward(it, myF, (Point(*myF) - Point(*it)), myUl, myUf, myLl, myLf); 

}


//--------------------------------------------------------------------
template <typename TIterator, typename TInteger, int connectivity>
inline
bool 
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::extendForward()
{

  ConstIterator it(myL); ++it;
  return extendForward(it, myL, (Point(*it) - Point(*myL)), myUf, myUl, myLf, myLl); 

}

//--------------------------------------------------------------------
template <typename TIterator, typename TInteger, int connectivity>
inline
bool 
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::extendBackward()
{

  ConstIterator it(myF); --it;
  return extendForward(it, myF, (Point(*myF) - Point(*it)), myUl, myUf, myLl, myLf); 

}

//--------------------------------------------------------------------

template <typename TIterator, typename TInteger, int connectivity>
inline
bool 
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>
::retractForward( ConstIterator & firstIt,
                  ConstIterator & lastIt,
                  ConstIterator & nextIt,       
                  Point & Uf,  Point & Ul,
                  Point & Lf,  Point & Ll,
                  const Integer& s )
{


  Vector v = Point(*myL)-Point(*myF);
  if (Tools<Integer,connectivity>
      ::norm(v[0],v[1]) <= NumberTraits<Integer>::ONE) {//if there are only two points

    //init with nextIt
    firstIt = nextIt;
    Uf = *firstIt;
    Lf = *firstIt;

    myA = 0; 
    myB = 0; 
    myMu = 0;
    myOmega = 0;
    myNbUpPat = 0;
    myNbLowPat = 0;

    mySteps.clear();

  } else {                   //there are strictly more than two points

    //update myF
    Point ptToRemove = *firstIt;
    firstIt = nextIt;

    //if the point to remove is the first upper leaning point...
    if (ptToRemove == Uf) {
      //... and if the first and last lower leaning points are confounded
      if (Lf == Ll) {
        //update of the parameters:
        //a,b
        Vector newMainVector = ( Lf - ( Uf + vectorFrom0ToOmega() ) )*s;
        myA = newMainVector[1];
        myB = newMainVector[0];
        //leaning points
        Integer k; //repetitions of the newMainVector
        Vector toUl = (Ul - Point(*firstIt))*s;
        k = Tools<Integer,connectivity>::norm(toUl[0],toUl[1])
          /Tools<Integer,connectivity>::norm(myA,myB);
        Uf = Ul - newMainVector*k*s;
        myNbUpPat = k;
        Vector toLastIt = Point(*lastIt) - Lf;
        k = Tools<Integer,connectivity>::norm(toLastIt[0],toLastIt[1])
          /Tools<Integer,connectivity>::norm(myA,myB);
        Ll = Lf + newMainVector*k*s;
        myNbLowPat = k;
        //mu, omega
        myMu = getRemainder(myUl);
        myOmega = Tools<Integer,connectivity>::norm(myA,myB);
      } else {
        Uf = Uf + Vector(myB,myA)*s;
        myNbUpPat--;
      }
    }

    //if the point to remove is the first lower leaning point...
    if (ptToRemove == Lf) {
      //... and if the first and last upper leaning points are confounded
      if (Uf == Ul) { 
        //update of the parameters:
        //a,b
        Vector newMainVector = ( Uf - (Lf - vectorFrom0ToOmega() ) )*s;
        myA = newMainVector[1];
        myB = newMainVector[0];
        //leaning points
        Integer k; //repetitions of the newMainVector
        Vector toLl = Ll - Point(*firstIt);
        k = Tools<Integer,connectivity>::norm(toLl[0],toLl[1])
          /Tools<Integer,connectivity>::norm(myA,myB);
        Lf = Ll - newMainVector*k*s;
        myNbLowPat = k;
        Vector toLastIt = Point(*lastIt) - Uf;
        k = Tools<Integer,connectivity>::norm(toLastIt[0],toLastIt[1])
          /Tools<Integer,connectivity>::norm(myA,myB);
        Ul = Uf + newMainVector*k*s;
        myNbUpPat = k;
        //mu, omega
        myMu = getRemainder(myUl);
        myOmega = Tools<Integer,connectivity>::norm(myA,myB);
      } else {
        Lf = Lf + Vector(myB,myA)*s;
        myNbLowPat--;
      }
    }

    //steps update
    if (Vector(myB,myA) == mySteps.at(0)) {
      Vector tmp(mySteps.at(0));
      mySteps.clear();
      mySteps.push_back(tmp);
    } else if (Vector(myB,myA) == mySteps.at(1)) {
      Vector tmp(mySteps.at(1));
      mySteps.clear();
      mySteps.push_back(tmp);
    }

  }

  return true;
}


//--------------------------------------------------------------------

template <typename TIterator, typename TInteger, int connectivity>
inline
bool 
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::retractForward()
{

  Vector v = Point(*myL)-Point(*myF);
  if ( (v[0] == 0) && (v[1] == 0) ) {//if there is only one point
    return false;
  } else {                           //otherwise
    ConstIterator next(myF); ++next;
    return retractForward(myF, myL, next, myUf, myUl, myLf, myLl, 1);
  }
}

//--------------------------------------------------------------------

template <typename TIterator, typename TInteger, int connectivity>
inline
bool 
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::retractBackward()
{

  Vector v = Point(*myL)-Point(*myF);
  if ( (v[0] == 0) && (v[1] == 0) ) {//if there is only one point
    return false;
  } else {                           //otherwise
    ConstIterator previous(myL); --previous;
    return retractForward(myL, myF, previous, myUl, myUf, myLl, myLf, -1);
  }
}

//--------------------------------------------------------------------
//Deprecated
template <typename TIterator, typename TInteger, int connectivity>
inline
bool 
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::isExtendableForward(const ConstIterator& it)
{

  Point lastPoint(*it);
  Vector lastMove = lastPoint-Point(*myL);

  Integer deviation = Tools<Integer,connectivity>
    ::norm(lastMove[0],lastMove[1]);
  //if the two last points are not connected return false
  if (deviation>NumberTraits<Integer>::ONE) return false;
  //if the two last points are confounded return true
  else if (deviation==NumberTraits<Integer>::ZERO) return true; 
  //otherwise

  else if (mySteps.size()<2) {

    //initialisation stage

    if (mySteps.size()==0) {

      return true;

    } else { //mySteps.size() == 1

      if (lastMove == mySteps.at(0)) {

        return true;

      } else {

        //difference of the two steps viewed as points
        Vector diff = ( lastMove-mySteps.at(0) );
        if (Tools<Integer,connectivity>::dualNorm(diff[0],diff[1]) == NumberTraits<Integer>::ONE) {
          //the two steps are compatible
          return true;

        } else {  //the two steps are not compatible
          return false;
        }
      }

    }

  } else {  //mySteps.size() == 2

    //main stage

    //if there only two steps at most
    if (hasLessThanTwoSteps(lastMove)) {
    
      //remainder
      Integer r = getRemainder(lastPoint);
          
      if ( (r < myMu-1)||(r > myMu+myOmega) )
        return false; //strongly exterior
      
      else { 
        return true;
      }

    } else { 
      //if there is more than two steps 
      return false;
    }

  } 

}

//--------------------------------------------------------------------
template <typename TIterator, typename TInteger, int connectivity>
inline
bool 
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::isExtendableForward()
{

  ConstIterator it(myL); ++it;
  return isExtendableForward(*it, (Point(*it) - Point(*myL)) ); 

}

//--------------------------------------------------------------------
template <typename TIterator, typename TInteger, int connectivity>
inline
bool 
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::isExtendableBackward()
{

  ConstIterator it(myF); --it;
  return isExtendableForward(*it, (Point(*myF) - Point(*it)) ); 

}

template <typename TIterator, typename TInteger, int connectivity>
inline
bool 
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::isExtendableForward(
                                                                             const Point & lastPoint, 
                                                                             const Vector & lastMove )
{


  Integer deviation = Tools<Integer,connectivity>
    ::norm(lastMove[0],lastMove[1]);
  //if the two last points are not connected return false
  if (deviation > NumberTraits<Integer>::ONE) return false;
  //if the two last points are confounded return true
  else if (deviation == NumberTraits<Integer>::ZERO) return true; 
  //otherwise

  else if (mySteps.size()<2) {
    
    //initialisation stage

    if (mySteps.size()==0) {

      return true;

    } else { //mySteps.size() == 1

      if (lastMove == mySteps.at(0)) {

        return true;

      } else {

        //difference of the two steps viewed as points
        Vector diff = ( lastMove-mySteps.at(0) );
        if (Tools<Integer,connectivity>::dualNorm(diff[0],diff[1]) == 1) {
          //the two steps are compatible
          return true;

        } else {  //the two steps are not compatible
          return false;
        }
      }
      
    }
    
  } else {  //mySteps.size() == 2

    //main stage
    
    //if there only two steps at most
    if (hasLessThanTwoSteps(lastMove)) {
    
      //remainder
      Integer r = getRemainder(lastPoint);
          
      if ( (r < myMu-1)||(r > myMu+myOmega) )
        return false; //strongly exterior
      else 
        return true;
      
    } else { 
      //if there is more than two steps
      return false;
    }

  } 

}


//-------------------------------------------------------------------------
template <typename TIterator, typename TInteger, int connectivity>
inline
TInteger
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::getRemainder(const Point & aPoint) const
{
  return myA * static_cast<Integer>(aPoint[0])
    - myB * static_cast<Integer>(aPoint[1]);
}


template <typename TIterator, typename TInteger, int connectivity>
inline
TInteger
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::getRemainder(const ConstIterator & it) const
{
  return getRemainder(*it);
}


template <typename TIterator, typename TInteger, int connectivity>
inline
TInteger
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::getPosition(const Point & aPoint) const
{
  return myA * static_cast<Integer>(aPoint[0])
    + myB * static_cast<Integer>(aPoint[1]);
}

template <typename TIterator, typename TInteger, int connectivity>
inline
TInteger
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::getPosition(const ConstIterator & it) const
{
  return getPosition(*it);
}


template <typename TIterator, typename TInteger, int connectivity>
inline
bool
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::isInDSL(const Point & aPoint) const
{
  Integer r = getRemainder(aPoint);
  return ( (r >= myMu)&&(r < myMu+myOmega) );
}


template <typename TIterator, typename TInteger, int connectivity>
inline
bool
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::isInDSL(const ConstIterator & it) const
{
  return isInDSL(*it);
}

//-------------------------------------------------------------------------
template <typename TIterator, typename TInteger, int connectivity>
inline
bool
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::isInDSS(const Point & aPoint) const
{
  Integer s = getPosition(aPoint);
  Integer smin = getPosition(myF);
  Integer smax = getPosition(myL);
  return (isInDSL(aPoint) && ( (s >= smin)&&(s <= smax) ) );
}


template <typename TIterator, typename TInteger, int connectivity>
inline
bool
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::isInDSS(const ConstIterator & it) const
{
  return isInDSS(*it);
}


//------------------------- Accessors --------------------------

template <typename TIterator, typename TInteger, int connectivity>
inline
TInteger
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::getA() const {
  return myA;
}

template <typename TIterator, typename TInteger, int connectivity>
inline
TInteger
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::getB() const {
  return myB;
}

template <typename TIterator, typename TInteger, int connectivity>
inline
TInteger
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::getMu() const {
  return myMu;
}

template <typename TIterator, typename TInteger, int connectivity>
inline
TInteger
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::getOmega() const {
  return myOmega;
}

template <typename TIterator, typename TInteger, int connectivity>
inline
typename DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::Point
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::getUf() const {
  return myUf;
}
template <typename TIterator, typename TInteger, int connectivity>
inline
typename DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::Point
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::getUl() const {
  return myUl;
}

template <typename TIterator, typename TInteger, int connectivity>
inline
typename DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::Point
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::getLf() const {
  return myLf;
}

template <typename TIterator, typename TInteger, int connectivity>
inline
typename DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::Point
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::getLl() const {
  return myLl;
}

template <typename TIterator, typename TInteger, int connectivity>
inline
typename DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::Point
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::getBackPoint() const {
  return (*myF);
}

template <typename TIterator, typename TInteger, int connectivity>
inline
typename DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::Point
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::getFrontPoint() const {
  return (*myL);
}

template <typename TIterator, typename TInteger, int connectivity>
inline
typename DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::Point
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::getFirstPoint() const {
  return (*myF);
}

template <typename TIterator, typename TInteger, int connectivity>
inline
typename DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::Point
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::getLastPoint() const {
  return (*myL);
}

template <typename TIterator, typename TInteger, int connectivity>
inline
TIterator
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::getBack() const {
  return myF;
}

template <typename TIterator, typename TInteger, int connectivity>
inline
TIterator
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::getFront() const {
  return myL;
}

template <typename TIterator, typename TInteger, int connectivity>
inline
TIterator
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::begin() const {
  return myF;
}

template <typename TIterator, typename TInteger, int connectivity>
inline
TIterator
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::end() const {
  ConstIterator i(myL); ++i;
  return i;
}

//-----------------------------------------------------------------
template <typename TIterator, typename TInteger, int connectivity>
inline
bool
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::isValid() const
{

  if (Point(*myF) == Point(*myL)) { 

    //init case with one point
    if ( (myA==0)
         &&(myB==0)
         &&(myMu==0)
         &&(myOmega==0)
         &&(myNbUpPat==0)
         &&(myNbLowPat==0) ) return true;
    else return false; 

  } else {

    //remainders of leaning points
    if (getRemainder(myUf) != myMu) return false;
    else if (getRemainder(myUl) != myMu) return false;
    else if (getRemainder(myLf) != myMu+myOmega-1) return false;
    else if (getRemainder(myLl) != myMu+myOmega-1) return false;
    //number of patterns
    else if ( (Vector(myB,myA)*myNbUpPat) != (myUl - myUf) ) return false;
    else if ( (Vector(myB,myA)*myNbLowPat) != (myLl - myLf) ) return false;
    else return true;

  }

}



template <typename TIterator, typename TInteger, int connectivity>
inline
typename DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::PointD
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::project( const Point & m ) const
{
  //casts
  double a = (double) NumberTraits<TInteger>::castToInt64_t(myA);
  double b = (double) NumberTraits<TInteger>::castToInt64_t(myB);
  double mu = (double) NumberTraits<TInteger>::castToInt64_t(myMu);
  double omega = (double) NumberTraits<TInteger>::castToInt64_t(myOmega);
  double xm = (double) NumberTraits<TInteger>::castToInt64_t(m[0]);
  double ym = (double) NumberTraits<TInteger>::castToInt64_t(m[1]);

  double alpha = ( mu + (omega - 1.0)/2.0 );
  double d2 = ( a * a + b * b );
  double s = b * xm + a * ym;
  double xp = ( b * s + a * alpha ) / d2;
  double yp = ( a * s - b * alpha ) / d2;
  return PointD( xp, yp );
}

template <typename TIterator, typename TInteger, int connectivity>
inline
typename DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::PointD
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::project( const Point & m, double r ) const
{
  //casts
  double a = (double) NumberTraits<TInteger>::castToInt64_t(myA);
  double b = (double) NumberTraits<TInteger>::castToInt64_t(myB);
  double xm = (double) NumberTraits<TInteger>::castToInt64_t(m[0]);
  double ym = (double) NumberTraits<TInteger>::castToInt64_t(m[1]);

  double d2 = ( a * a + b * b );
  double s = b * xm + a * ym;
  double xp = ( b * s + a * r ) / d2;
  double yp = ( a * s - b * r ) / d2;
  
  return PointD( xp, yp );
}


template <typename TIterator, typename TInteger, int connectivity>
inline
typename DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::PointD
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::project( const Point & m, const Point & p ) const
{

  double r = (double) NumberTraits<TInteger>::castToInt64_t(getRemainder(p));
  return project(m,r);

}

template <typename TIterator, typename TInteger, int connectivity>
inline
double
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::projectedSegmentLength() const
{
  PointD v = project( *myL );
  PointD u = project( *myF );
  v -= u;
  return v.norm(PointD::L_2);
}


template <typename TIterator, typename TInteger, int connectivity>
inline
std::string
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::className() const
{
  return "ArithmeticalDSS";
}

template <typename TIterator, typename TInteger, int connectivity>
inline
void
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::selfDisplay ( std::ostream & out)
{
  
  out << "[ArithmeticalDSS]" <<  endl;
  out << "Parameters (a,b,mu,omega)=";
  out << "("<< myA << ", " << myB << ", ";
  out << myMu << ", " << myOmega << ")" <<  endl;
  out << "Number of upper patterns: " << myNbUpPat << endl;
  out << "Number of lower patterns: " << myNbLowPat << endl;
  out << "First point " << Point(*myF) << " Last point " << Point(*myL) <<  endl;
  out << "Leaning points:" <<  endl;
  out << "   Uf " << myUf << endl << "   Ul " << myUl << endl;
  out << "   Lf " << myLf << endl << "   Ll " << myLl <<  endl;
  out << "Steps:" <<  endl;
  for (unsigned int i = 0; i < mySteps.size(); i++) {
    out << "   " << mySteps.at(i) << endl;
  }
  out << "[End ArithmeticalDSS]" <<  endl;
  
}

//-------------------------------------------------------------------------
template <typename TIterator, typename TInteger, int connectivity>
inline
bool
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::hasLessThanTwoSteps(
                                                                             const Vector& aStep) const
{
  if ( (aStep == mySteps.at(0)) || 
       (aStep == mySteps.at(1)) ) {
    return true;
  } else {
    return false;
  }

}

template <typename TIterator, typename TInteger, int connectivity>
inline
typename DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::Vector
DGtal::ArithmeticalDSS<TIterator,TInteger,connectivity>::vectorFrom0ToOmega() const
{
  Vector v = mySteps.at(1) - mySteps.at(0);
  if ( getRemainder(v) == myOmega) {
    return v;
  } else {
    return ( mySteps.at(0) - mySteps.at(1) );
  }

}