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.a());      
        Value b = (Value) NumberTraits<typename DSS::Integer>
          ::castToDouble(aDSS.b());      
 
              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.b() ); 
        double y = NumberTraits<typename DSS::Integer>
          ::castToDouble( aDSS.a() );
        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.b(), aDSS.a());
      }
    }; 
    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::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.getStabbingLineComputerPtr()->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::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.getStabbingLineComputerPtr()->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.b(), aDSS.a() );
        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(( concepts::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::L_1) <= NumberTraits<Integer>::ONE);  
        bool isConnectedAtFront = isNotEmpty(front, myEnd)
          &&((*boost::prior(front)-*front).norm(Vector::L_1) <= NumberTraits<Integer>::ONE);  
  
        Integer mu = mySCPtr->mu();
        Integer omega = mySCPtr->omega();
        if (isConnectedAtBack)
        {
          --back;
          //parameters
                if (isConnectedAtFront)
          {
            if ((mySCPtr->remainder(*back)<=mu-1)&&
                (mySCPtr->remainder(*front)<=mu-1) )
            {
              //convex
              k = myFunctor(*mySCPtr) / myH;
            }
            else if ( (mySCPtr->remainder(*back)>=mu+omega)&&
                      (mySCPtr->remainder(*front)>=mu+omega) )
            {
              //concave
              k = -myFunctor(*mySCPtr) / myH;
            } //else inflection
              }
          else
          {
            if ( (mySCPtr->remainder(*back)<=mu-1) )
            {
              //convex
              k = myFunctor(*mySCPtr) / myH; 
            }
            else if ( (mySCPtr->remainder(*back)>=mu+omega) )
            {
              //concave
              k = -myFunctor(*mySCPtr) / myH; 
            } //else inflection
                }
              }
        else if (isConnectedAtFront)
        {
          if ( (mySCPtr->remainder(*front)<=mu-1) )
          {
            //convex
            k = myFunctor(*mySCPtr) / myH; 
          }
          else if ( (mySCPtr->remainder(*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)