SegmentComputerEstimators.h

#pragma once

#if defined(SegmentComputerEstimators_RECURSES)
#error Recursive header files inclusion detected in SegmentComputerEstimators.h
#else // defined(SegmentComputerEstimators_RECURSES)

#define SegmentComputerEstimators_RECURSES

#if !defined SegmentComputerEstimators_h

#define SegmentComputerEstimators_h

// Inclusions
#include <cmath>
#include "DGtal/base/Common.h"
#include "DGtal/helpers/StdDefs.h"
#include "boost/utility.hpp"



namespace DGtal
{

  namespace detail
  {
    // class PosIndepScaleIndepSCEstimator

    template <typename TSegmentComputer, typename Functor, 
              typename ReturnType = typename Functor::Value>
    class PosIndepScaleIndepSCEstimator
    {

    public: 

      // ----------------------- inner type ------------------------------
      typedef TSegmentComputer SegmentComputer;
      typedef typename SegmentComputer::ConstIterator ConstIterator; 
      typedef ReturnType Quantity;

      // ----------------------- Internal data  ------------------------------
    public:
      ConstIterator myBegin; 
      ConstIterator myEnd;
      const SegmentComputer* mySCPtr; 
      Functor myFunctor; 

      // ----------------------- Standard services ------------------------------
    public:

      bool isValid() const
      {
        return (mySCPtr != 0); 
      };

      // ----------------------- Interface --------------------------------------
    public:

      PosIndepScaleIndepSCEstimator(): myBegin(), myEnd(), mySCPtr(0), myFunctor() {}

      void init(const double /*h*/, const ConstIterator& itb, const ConstIterator& ite)
      {
        myBegin = itb;
        myEnd = ite; 
      }

      void attach(const SegmentComputer& aSC) 
      {
        mySCPtr = &aSC; 
      };

      Quantity eval(const ConstIterator& /*it*/) const
      {
        ASSERT( mySCPtr ); 
        return myFunctor( *mySCPtr ); 
      }

      template <typename OutputIterator>
      OutputIterator eval(const ConstIterator& itb, const ConstIterator& ite, 
                          OutputIterator result) const
      {
        ASSERT( mySCPtr ); 

        // do-while loop to deal with the case of a whole circular range
        if (isNotEmpty(itb, ite))
          {
            ConstIterator it = itb; 
            do
              {
                *result++ = myFunctor( *mySCPtr ); 
                ++it; 
              } while (it != ite);
          }

        return result; 
      }

    }; // end of class PosIndepScaleIndepSCEstimator

    // class PosIndepScaleDepSCEstimator

    template <typename TSegmentComputer, typename Functor, 
              typename ReturnType = typename Functor::Value>
    class PosIndepScaleDepSCEstimator
    {

    public: 

      // ----------------------- inner type ------------------------------
      typedef TSegmentComputer SegmentComputer;
      typedef typename SegmentComputer::ConstIterator ConstIterator; 
      typedef ReturnType Quantity;

      // ----------------------- internal data ------------------------------
    public:
      double myH; 
      ConstIterator myBegin; 
      ConstIterator myEnd;
      const SegmentComputer* mySCPtr; 
      Functor myFunctor; 

      // ----------------------- Standard services ------------------------------
    public:

      PosIndepScaleDepSCEstimator()
        : myH( 0.0 ), myBegin(), myEnd(), mySCPtr(0), myFunctor()
      {
      }
      PosIndepScaleDepSCEstimator( const PosIndepScaleDepSCEstimator& other )
        : myH( other.myH ), myBegin( other.myBegin ), myEnd( other.myEnd ), 
          mySCPtr( other.mySCPtr ), myFunctor( other.myFunctor )
      {
      }
      PosIndepScaleDepSCEstimator& operator=( const PosIndepScaleDepSCEstimator& other )
      {
        if (this != &other)
          {
            myH = other.myH; 
            myBegin = other.myBegin; 
            myEnd = other.myEnd; 
            mySCPtr = other.mySCPtr; 
            myFunctor = other.myFunctor; 
          }
        return *this; 
      }
      ~PosIndepScaleDepSCEstimator() {}

      bool isValid() const 
      {
        return (myH > 0)&&(mySCPtr != 0); 
      };



      // ----------------------- Interface --------------------------------------
    public:

      void init(const double h, const ConstIterator& itb, const ConstIterator& ite)
      {
        myH = h; 
        myBegin = itb;
        myEnd = ite;
        ASSERT( myH > 0 );  
      }

      void attach(const SegmentComputer& aSC) 
      {
        mySCPtr = &aSC; 
        ASSERT( mySCPtr ); 
      };

      Quantity eval(const ConstIterator& /*it*/) const
      {
        ASSERT( isValid() ); 
        return myFunctor( *mySCPtr, myH ); 
      }

      template <typename OutputIterator>
      OutputIterator eval(const ConstIterator& itb, const ConstIterator& ite, 
                          OutputIterator result) const
      {
        ASSERT( isValid() ); 

        // do-while loop to deal with the case of a whole circular range
        if (isNotEmpty(itb, ite))
          {
            ConstIterator it = itb; 
            do
              {
                *result++ = myFunctor( *mySCPtr, myH ); 
                ++it; 
              } while (it != ite);
          }

        return result; 
      }

    }; // end of class PosIndepScaleDepSCEstimator

    // class PosDepScaleIndepSCEstimator

    template <typename TSegmentComputer, typename Functor, 
              typename ReturnType = typename Functor::Value>
    class PosDepScaleIndepSCEstimator
    {

    public: 

      // ----------------------- inner type ------------------------------
      typedef TSegmentComputer SegmentComputer;
      typedef typename SegmentComputer::ConstIterator ConstIterator; 
      typedef ReturnType Quantity;

      // ----------------------- Internal data  ------------------------------
    public:
      ConstIterator myBegin; 
      ConstIterator myEnd;
      const SegmentComputer* mySCPtr; 
      Functor myFunctor; 

      // ----------------------- Standard services ------------------------------
    public:

      bool isValid() const 
      {
        return (mySCPtr != 0); 
      };

      // ----------------------- Interface --------------------------------------
    public:

      PosDepScaleIndepSCEstimator(): myBegin(), myEnd(), mySCPtr(0), myFunctor() {}

      void init(const double /*h*/, const ConstIterator& itb, const ConstIterator& ite)
      {
        myBegin = itb;
        myEnd = ite; 
      }

      void attach(const SegmentComputer& aSC) 
      {
        mySCPtr = &aSC; 
      };

      Quantity eval(const ConstIterator& it) const
      {
        ASSERT( mySCPtr ); 
        return myFunctor( it, *mySCPtr ); 
      }

      template <typename OutputIterator>
      OutputIterator eval(const ConstIterator& itb, const ConstIterator& ite, 
                          OutputIterator result) const
      {
        ASSERT( mySCPtr ); 

        // do-while loop to deal with the case of a whole circular range
        if (isNotEmpty(itb, ite))
          {
            ConstIterator it = itb; 
            do
              {
                *result++ = myFunctor( it, *mySCPtr ); 
                ++it; 
              } while (it != ite);
          }

        return result; 
      }

    }; // end of class PosDepScaleIndepSCEstimator

    // class PosDepScaleDepSCEstimator

    template <typename TSegmentComputer, typename Functor, 
              typename ReturnType = typename Functor::Value>
    class PosDepScaleDepSCEstimator
    {

    public: 

      // ----------------------- inner type ------------------------------
      typedef TSegmentComputer SegmentComputer;
      typedef typename SegmentComputer::ConstIterator ConstIterator; 
      typedef ReturnType Quantity;

      // ----------------------- internal data ------------------------------
    public:
      double myH; 
      ConstIterator myBegin; 
      ConstIterator myEnd;
      const SegmentComputer* mySCPtr; 
      Functor myFunctor; 

      // ----------------------- Standard services ------------------------------
    public:

      PosDepScaleDepSCEstimator()
        : myH( 0.0 ), myBegin(), myEnd(), mySCPtr(0), myFunctor()
      {
      }
      PosDepScaleDepSCEstimator( const PosDepScaleDepSCEstimator& other )
        : myH( other.myH ), myBegin( other.myBegin ), myEnd( other.myEnd ), 
          mySCPtr( other.mySCPtr ), myFunctor( other.myFunctor )
      {
      }
      PosDepScaleDepSCEstimator& operator=( const PosDepScaleDepSCEstimator& other )
      {
        if (this != &other)
          {
            myH = other.myH; 
            myBegin = other.myBegin; 
            myEnd = other.myEnd; 
            mySCPtr = other.mySCPtr; 
            myFunctor = other.myFunctor; 
          }
        return *this; 
      }
      ~PosDepScaleDepSCEstimator() {}

      bool isValid() const 
      {
        return (myH > 0)&&(mySCPtr != 0); 
      };



      // ----------------------- Interface --------------------------------------
    public:

      void init(const double h, const ConstIterator& itb, const ConstIterator& ite)
      {
        myH = h; 
        myBegin = itb;
        myEnd = ite;
        ASSERT( myH > 0 );  
      }

      void attach(const SegmentComputer& aSC) 
      {
        mySCPtr = &aSC; 
        ASSERT( mySCPtr ); 
      };

      Quantity eval(const ConstIterator& it) const
      {
        ASSERT( isValid() ); 
        return myFunctor( it, *mySCPtr, myH ); 
      }

      template <typename OutputIterator>
      OutputIterator eval(const ConstIterator& itb, const ConstIterator& ite, 
                          OutputIterator result) const 
      {
        ASSERT( isValid() ); 

        // do-while loop to deal with the case of a whole circular range
        if (isNotEmpty(itb, ite))
          {
            ConstIterator it = itb; 
            do
              {
                *result++ = myFunctor( it, *mySCPtr, myH ); 
                ++it; 
              } while (it != ite);
          }

        return result; 
      }

    }; // end of class PosDepScaleDepSCEstimator


    struct TangentAngleFromDSS
    {      
    public:
      typedef double Value; 


      template<typename DSS>
      Value operator() (const DSS& aDSS) const 
      {
        Value a = (Value) NumberTraits<typename DSS::Integer>
          ::castToDouble(aDSS.getA());      
        Value b = (Value) NumberTraits<typename DSS::Integer>
          ::castToDouble(aDSS.getB());      

        return std::atan2(a,b);
      }
    }; 
    struct NormalizedTangentVectorFromDSS
    {      
    public:
      typedef DGtal::PointVector<2,double> RealVector; 
      typedef RealVector Value;

      template<typename DSS>
      Value operator() (const DSS& aDSS) const 
      {
        double x = NumberTraits<typename DSS::Integer>
	  ::castToDouble( aDSS.getB() ); 
        double y = NumberTraits<typename DSS::Integer>
	  ::castToDouble( aDSS.getA() );
        RealVector v(x,y); 
        double norm = v.norm(RealVector::L_2);
        v /= norm; 
        return v;
      }
    }; 
    template<typename DSS>
    struct TangentVectorFromDSS
    {      
    public:
      typedef typename DSS::Vector Value;

      Value operator() (const DSS& aDSS) const 
      {
        return Value(aDSS.getB(), aDSS.getA());
      }
    }; 
    template<bool isCCW = true>
    struct CurvatureFromDCA
    {      
    public:
      typedef double Value; 

      template<typename DCA>
      Value operator() (const DCA& aDCA, const double& aH = 1.0) const 
      {
        if ( aDCA.isStraight() )
          return 0.0; 
        else
          return ( aDCA.getSeparatingCircle().getCurvature() / aH );
      }
    }; 
    template<>
    struct CurvatureFromDCA<false>
    {      
    public:
      typedef double Value; 

      template<typename DCA>
      Value operator() (const DCA& aDCA, const Value& aH = 1.0) const 
      {
        if ( aDCA.isStraight() )
          return 0.0; 
        else
          return - ( aDCA.getSeparatingCircle().getCurvature() / aH );
      }
    }; 
    struct NormalVectorFromDCA
    {      
    public:
      typedef PointVector<2,double> Value; 


      template<typename DCA>
      Value operator() (const typename DCA::ConstIterator& it, 
                        const DCA& aDCA) const 
      {
        typedef typename DCA::ConstIterator ConstIterator; 
        typedef typename DCA::Pair Pair; 
        typedef typename DCA::Point Point;
        typedef typename Point::Coordinate Coordinate; 
        
        if ( !aDCA.isStraight() )
          {
            //separating circle center
            double c0, c1, r; 
            aDCA.getSeparatingCircle().getParameters(c0, c1, r);
            //point
            Pair pair = *it; 
            Point i = pair.first; 
            Point o = pair.second;
            double m0 = NumberTraits<Coordinate>::castToDouble(i[0]+o[0]) / 2.0; 
            double m1 = NumberTraits<Coordinate>::castToDouble(i[1]+o[1]) / 2.0;
            //normal vector 
            double v0 = m0 - c0; 
            double v1 = m1 - c1; 
            //norm
            double n = std::sqrt(v0*v0 + v1*v1); 
            return Value( v0/n, v1/n );
          }
        else
          {
            //separating straight line and normal vector
            double a, b, c; 
            aDCA.getGeometricalDSSPtr()->getParameters(a, b, c);
            //norm
            double n = std::sqrt(a*a + b*b); 
            return Value( a/n, b/n ); 
          }
      }
    }; 

    struct TangentVectorFromDCA
    {      
    public:
      typedef PointVector<2,double> Value; 

      template<typename DCA>
      Value operator() (const typename DCA::ConstIterator& it, 
                        const DCA& aDCA) const 
      {
        NormalVectorFromDCA f; 
        Value normal = f(it, aDCA); 
        return Value( normal[1], normal[0] ); 
      }
    }; 

    struct DistanceFromDCA
    {      
    public:
      typedef std::pair<double,double> Value; 

      template<typename DCA>
      Value operator() (const typename DCA::ConstIterator& it, 
                        const DCA& aDCA, const double& aH) const 
      {
        typedef typename DCA::ConstIterator ConstIterator; 
        typedef typename DCA::Pair Pair; 
        typedef typename DCA::Point Point;
        typedef typename Point::Coordinate Coordinate; 
        
        if ( !aDCA.isStraight() )
          {
            //separating circle center
            double c0, c1, r; 
            aDCA.getSeparatingCircle().getParameters(c0, c1, r);
            //points
            Pair pair = *it; 
            Point i = pair.first; 
            Point o = pair.second;
            //distances
            double distI0 = NumberTraits<Coordinate>::castToDouble(i[0]) - c0; 
            double distI1 = NumberTraits<Coordinate>::castToDouble(i[1]) - c1;
            double distI = std::sqrt( distI0*distI0 + distI1*distI1 ) - r; 
            double distO0 = NumberTraits<Coordinate>::castToDouble(o[0]) - c0; 
            double distO1 = NumberTraits<Coordinate>::castToDouble(o[1]) - c1;
            double distO = std::sqrt( distO0*distO0 + distO1*distO1 ) - r; 
            return Value( distI*aH, distO*aH );
          }
        else
          {
            //separating straight line
            double a, b, c; 
            aDCA.getGeometricalDSSPtr()->getParameters(a, b, c); 
            //norm
            double n = std::sqrt(a*a + b*b); 
            //points
            Pair pair = *it; 
            Point i = pair.first; 
            Point o = pair.second;
            //distances
            double rI = NumberTraits<Coordinate>::castToDouble(i[0])*a + 
              NumberTraits<Coordinate>::castToDouble(i[1])*b + c;
            double distI = rI / n; 
            double rO = NumberTraits<Coordinate>::castToDouble(o[0])*a + 
              NumberTraits<Coordinate>::castToDouble(o[1])*b + c;
            double distO = rO / n; 
            return Value( distI*aH, distO*aH ); 
          }
      }
    }; 

  }//namespace detail



  //-------------------------------------------------------------------------------------------
  template <typename DSSComputer>
  class TangentFromDSSEstimator: 
    public detail::PosIndepScaleIndepSCEstimator<DSSComputer, detail::NormalizedTangentVectorFromDSS>
  {
    typedef 
    detail::PosIndepScaleIndepSCEstimator<DSSComputer, detail::NormalizedTangentVectorFromDSS> 
    Super; 

  public: 
    TangentFromDSSEstimator(): Super() {};
    TangentFromDSSEstimator( const TangentFromDSSEstimator & other ): Super(other) {};
  }; 

  //-------------------------------------------------------------------------------------------
  template <typename DSSComputer>
  class TangentVectorFromDSSEstimator: 
    public detail::PosIndepScaleIndepSCEstimator<DSSComputer, detail::TangentVectorFromDSS<DSSComputer> >
  {
    typedef 
    detail::PosIndepScaleIndepSCEstimator<DSSComputer, detail::TangentVectorFromDSS<DSSComputer> > 
    Super; 

  public: 
    TangentVectorFromDSSEstimator(): Super() {};
    TangentVectorFromDSSEstimator( const TangentVectorFromDSSEstimator & other ): Super(other) {};
  }; 

  //-------------------------------------------------------------------------------------------
  template <typename DSSComputer>
  class TangentAngleFromDSSEstimator: 
    public detail::PosIndepScaleIndepSCEstimator<DSSComputer, detail::TangentAngleFromDSS>
  {
    typedef 
    detail::PosIndepScaleIndepSCEstimator<DSSComputer, detail::TangentAngleFromDSS> 
    Super; 

  public: 
    TangentAngleFromDSSEstimator(): Super() {};
    TangentAngleFromDSSEstimator( const TangentAngleFromDSSEstimator & other ): Super(other) {};
  }; 

  //-------------------------------------------------------------------------------------------
  template <typename DCAComputer, bool isCCW = true>
  class CurvatureFromDCAEstimator: 
    public detail::PosIndepScaleDepSCEstimator<DCAComputer, 
                                               detail::CurvatureFromDCA<isCCW> >
  {
    typedef 
    detail::PosIndepScaleDepSCEstimator<DCAComputer, detail::CurvatureFromDCA<isCCW> > 
    Super; 

  public: 
    CurvatureFromDCAEstimator(): Super() {};
    CurvatureFromDCAEstimator( const CurvatureFromDCAEstimator & other ): Super(other) {};
  }; 

  //-------------------------------------------------------------------------------------------
  template <typename DCAComputer>
  class NormalFromDCAEstimator: 
    public detail::PosDepScaleIndepSCEstimator<DCAComputer, 
                                               detail::NormalVectorFromDCA>
  {
    typedef 
    detail::PosDepScaleIndepSCEstimator<DCAComputer, detail::NormalVectorFromDCA> 
    Super; 

  public: 
    NormalFromDCAEstimator(): Super() {};
    NormalFromDCAEstimator( const NormalFromDCAEstimator & other ): Super(other) {};
  }; 

  //-------------------------------------------------------------------------------------------
  template <typename DCAComputer>
  class TangentFromDCAEstimator: 
    public detail::PosDepScaleIndepSCEstimator<DCAComputer, 
                                               detail::TangentVectorFromDCA>
  {
    typedef 
    detail::PosDepScaleIndepSCEstimator<DCAComputer, detail::TangentVectorFromDCA> 
    Super; 

  public: 
    TangentFromDCAEstimator(): Super() {};
    TangentFromDCAEstimator( const TangentFromDCAEstimator & other ): Super(other) {};
  }; 

  //-------------------------------------------------------------------------------------------
  template <typename DCAComputer>
  class DistanceFromDCAEstimator: 
    public detail::PosDepScaleDepSCEstimator<DCAComputer, 
                                             detail::DistanceFromDCA>
  {
    typedef 
    detail::PosDepScaleDepSCEstimator<DCAComputer, detail::DistanceFromDCA> 
    Super; 

  public: 
    DistanceFromDCAEstimator(): Super() {};
    DistanceFromDCAEstimator( const DistanceFromDCAEstimator & other ): Super(other) {};
  }; 

  namespace detail
  {
    struct CurvatureFromDSSLength
    {      
    public:
      typedef double Value; 

      template<typename DSS>
      Value operator() (const DSS& aDSS) const 
      {
        typedef typename DSS::Vector Vector; 
        //length
        Vector v = ( *aDSS.begin() - *boost::prior(aDSS.end()) ); 
        Value l = v.norm(Vector::L_2);
        //result
        return 1/( (l*l)/8 + 0.5 );  
      }
    }; 

    struct CurvatureFromDSSLengthAndWidth
    {      
    public:
      typedef double Value; 

      template<typename DSS>
      Value operator() (const DSS& aDSS) const 
      {
        typedef typename DSS::Vector Vector; 
        //length
        Vector v = ( *aDSS.begin() - *boost::prior(aDSS.end()) ); 
        Value l = v.norm(Vector::L_2);
        //width
        Vector t( aDSS.getB(), aDSS.getA() );
        Value w = 1.0 / v.norm(Vector::L_2); 
        //result
        return 1.0/( (l*l)/(8*w) + w/2 ); 
      }
    }; 

    // class CurvatureFromDSSBaseEstimator

    template <typename DSSComputer, typename Functor = detail::CurvatureFromDSSLength >
    class CurvatureFromDSSBaseEstimator
    {

    public: 

      // ----------------------- inner type ------------------------------------
      typedef DSSComputer SegmentComputer; 
      typedef typename DSSComputer::ConstIterator ConstIterator; 
      typedef double Quantity;

      BOOST_CONCEPT_ASSERT(( CUnaryFunctor< Functor, SegmentComputer, Quantity > ));  

      // ----------------------- internal data ------------------------------
    public:
      double myH; 
      ConstIterator myBegin; 
      ConstIterator myEnd;
      const SegmentComputer* mySCPtr; 
      Functor myFunctor; 

      // ----------------------- Standard services ------------------------------
    public:

      CurvatureFromDSSBaseEstimator()
        : myH( 0.0 ), myBegin(), myEnd(), mySCPtr(0), myFunctor()
      {
      }
      CurvatureFromDSSBaseEstimator( const CurvatureFromDSSBaseEstimator& other )
        : myH( other.myH ), myBegin( other.myBegin ), myEnd( other.myEnd ), 
          mySCPtr( other.mySCPtr ), myFunctor( other.myFunctor )
      {
      }
      CurvatureFromDSSBaseEstimator& operator=( const CurvatureFromDSSBaseEstimator& other )
      {
        if (this != &other)
          {
            myH = other.myH; 
            myBegin = other.myBegin; 
            myEnd = other.myEnd; 
            mySCPtr = other.mySCPtr; 
            myFunctor = other.myFunctor; 
          }
        return *this; 
      }
      ~CurvatureFromDSSBaseEstimator() {}

      bool isValid() const
      {
        return (myH > 0)&&(mySCPtr != 0); 
      };

      // ----------------------- Interface --------------------------------------
    public:

      void init(const double h, const ConstIterator& itb, const ConstIterator& ite)
      {
        myH = h; 
        myBegin = itb;
        myEnd = ite;
        ASSERT( myH > 0 );  
      }

      Quantity eval(const ConstIterator& it)
      {
        ASSERT( isValid() ); 

        //types
        typedef typename DSSComputer::Integer Integer; 
        typedef typename DSSComputer::Vector Vector; 

        //curvature value
        Quantity k = 0;  


        //begin and end iterators
        //(back point on the first point)
        //(front point on the last point)
        ConstIterator back = mySCPtr->begin();  
        ConstIterator front = mySCPtr->end();
        bool isConnectedAtBack = isNotEmpty(myBegin, back)
          &&((*boost::prior(back)-*back).norm(Vector::L1) <= NumberTraits<Integer>::ONE);  
        bool isConnectedAtFront = isNotEmpty(front, myEnd)
          &&((*boost::prior(front)-*front).norm(Vector::L1) <= NumberTraits<Integer>::ONE);  
  

        if (isConnectedAtBack) {
          if (isConnectedAtFront) {

            --back;

            //parameters
            Integer mu = mySCPtr->getMu();
            Integer omega = mySCPtr->getOmega();

            //cases
            if ( (mySCPtr->getRemainder(*back)<=mu-1)&&
                 (mySCPtr->getRemainder(*front)<=mu-1) ) {                //convex
              k = myFunctor(*mySCPtr) / myH; 
            } else if ( (mySCPtr->getRemainder(*back)>=mu+omega)&&
                        (mySCPtr->getRemainder(*front)>=mu+omega) ) {           //concave
              k = -myFunctor(*mySCPtr) / myH; 
            } //else                                                  //inflection

          } else {

            --back;

            //parameters
            Integer mu = mySCPtr->getMu();
            Integer omega = mySCPtr->getOmega();

            //cases
            if ( (mySCPtr->getRemainder(*back)<=mu-1) ) {                //convex
              k = myFunctor(*mySCPtr) / myH; 
            } else if ( (mySCPtr->getRemainder(*back)>=mu+omega) ) {     //concave
              k = -myFunctor(*mySCPtr) / myH; 
            } //else                                                 //inflection

          }
        } else if (isConnectedAtFront) {

          //parameters
          Integer mu = mySCPtr->getMu();
          Integer omega = mySCPtr->getOmega();

          //cases
          if ( (mySCPtr->getRemainder(*front)<=mu-1) ) {                //convex
            k = myFunctor(*mySCPtr) / myH; 
          } else if ( (mySCPtr->getRemainder(*front)>=mu+omega) ) {     //concave
            k = -myFunctor(*mySCPtr) / myH; 
          } //else                                                  //inflection

        } //else cannot be extended: k is set to 0

        return k;

      }

      template <typename OutputIterator>
      OutputIterator eval(const ConstIterator& itb, const ConstIterator& ite, 
                          OutputIterator result)
      {
        ASSERT( isValid() ); 

        // do-while loop to deal with the case of a whole circular range
        if (isNotEmpty(itb, ite))
          {
            ConstIterator it = itb; 
            do
              {
                *result++ = eval( it ); 
                ++it; 
              } while (it != ite);
          }

        return result; 
      }

      void attach(const SegmentComputer& aSC) 
      {
        mySCPtr = &aSC; 
        ASSERT( mySCPtr ); 
      };


    }; // end of class CurvatureFromDSSBaseEstimator

  }//namespace detail



  //-------------------------------------------------------------------------------------------
  template <typename DSSComputer>
  class CurvatureFromDSSLengthEstimator: 
    public detail::CurvatureFromDSSBaseEstimator<DSSComputer, detail::CurvatureFromDSSLength >
  {

    typedef detail::CurvatureFromDSSBaseEstimator<DSSComputer, detail::CurvatureFromDSSLength > Super; 
 
  public: 
    CurvatureFromDSSLengthEstimator(): Super() {};
    CurvatureFromDSSLengthEstimator( const CurvatureFromDSSLengthEstimator & other ): Super(other) {};
  }; 

  //-------------------------------------------------------------------------------------------
  template <typename DSSComputer>
  class CurvatureFromDSSEstimator: 
    public detail::CurvatureFromDSSBaseEstimator<DSSComputer, detail::CurvatureFromDSSLengthAndWidth >
  {

    typedef detail::CurvatureFromDSSBaseEstimator<DSSComputer, detail::CurvatureFromDSSLengthAndWidth > Super; 
 
  public: 
    CurvatureFromDSSEstimator(): Super() {};
    CurvatureFromDSSEstimator( const CurvatureFromDSSEstimator & other ): Super(other) {};
  }; 


} // namespace DGtal

//                                                                           //

#endif // !defined SegmentComputerEstimators_h

#undef SegmentComputerEstimators_RECURSES
#endif // else defined(SegmentComputerEstimators_RECURSES)