How to Iterate Over the Parametric Curve

• Download source code - 5.1 KB

Introduction

In working with the 3D application, we widely use curves as on plane so in 3D space. These curves often have parametric nature. Developers write cycles for walking along these curves repeatedly, again and again. I could pick out some of the tasks:

1. Walk along the curve with constant step
2. Walk along the curve with dynamically changing step
3. Stop walking by some reason

To solve these tasks, I will use some kind of iterator and this iterator will change its behaviour according to the iteration strategy.

Background

Well known mathematical abstraction Parametric Function maps interval of real numbers [start, stop] to some values in the range of function:

v = f(t)

The Intervals of Real Numbers ValueRange was described in this article.

In the common case, the type of values in the range of function and the functions itself may be different. But in 3D applications, they are Cartesian points on a plane (2D) or in space (3D) and curves could be represented as:

Point p = curve(t)

One of the most common tasks is a point retrieving while moving along the given curve, thus the developers have to write such code again and again. Such code is not trivial. It must take into account the closeness of the curve as one of the aspects.

Moreover, for a lot of purposes, we have to change the iteration strategy. For example, to visualize a curve, which curvature changes for different values of the parameter, we could use dynamically changing iteration step and there are a lot of such useful iteration strategies.

Implementation

We will start with IteratorBase to afford to write other iterators. It will contain a reference to some data which is implied to be a parametric function of any kind.

// iterator base
// 1. defines Function type
// 2. stores ref to function
// 3. declares equality operator
//
// Function - Parametric Function v = f(t)
//
template<typename Function>
class IteratorBase
{
public:
    using FunctionType = Function;

    const FunctionType& getFunction() const; 

protected:
    explicit IteratorBase(const FunctionType& function);

    template <typename Iterator>
    bool operator == (const Iterator& other) const;

protected:
    const FunctionType& m_function;
};

The next step is a generic iterator which will afford us to use it in ranged for cycles. This evolution:

1. contains the interval of real numbers [start, stop] in function definition domain
2. defines begin/end pair
3. defines operator of the increment
4. defines the equality operator

// base iterator with range
//
template<typename Function, typename Derived>
class RangeIterator : public IteratorBase<Function>
{
    using Inherited = IteratorBase<Function>;
public:
    using DerivedType = Derived;
    using ValueRange = typename Function::ValueRange;
    using ValueType = typename Function::ValueType;

protected:

    explicit RangeIterator(const typename Inherited::FunctionType& function)
        : Inherited(function)
        , m_range(function.getParamValueRange())
    {}

    RangeIterator(const Function& function, const ValueRange& range)
        : Inherited(function)
        , m_range(range)
    {}

public:
    DerivedType& begin()
    {
        m_u = m_range.start;
        return static_cast<DerivedType&>(*this);
    }

    const TheEndIterator& end() const
    {
        static const TheEndIterator TheEnd;
        return TheEnd;
    }

    bool operator == (const DerivedType& other) const;
    bool operator == (const TheEndIterator&) const;

    template <typename Iterator>
    bool operator != (const Iterator& other) const;

    DerivedType& operator++()
    {
        if (isNear(m_u, m_range.stop))
        {
            // last step always done this iterator enters beyond range.stop
            m_u = VeryLarge;
            // now it equals to end();
        }
        else
        {
            // static polymorphism: children knows better how to step
            double du = VeryLarge;
            if (static_cast<DerivedType&>(*this).getNextStep(du))
            {
                m_u += du;
                if (m_u >= m_range.stop)
                {
                    // next step is too long, let's cut m_u sharp to range.stop
                    // this can be caused by rounding errors
                    // or error in next step approximation
                    m_u = m_range.stop;
                }
            }
            else
            {
                m_u = VeryLarge;
            }
        }
        return static_cast<DerivedType&>(*this);
    }

    double getU() const
    {
        return m_u;
    }

    ValueType operator*() const
    {
        return Inherited::m_function.valueAt(getU());
    }

protected:
    const ValueRange m_range;
    double m_u = 0; // current position in u-domain
    static constexpr double VeryLarge = std::numeric_limits<double>::max();
};

The template parameters of this iterator are:

1. Function - Parametric Function (v = f(t)). The Function must provide:
   1. The type ValueRange and getter for function domain ValueRange
   2. The type ValueType and getter for value at parameter t
      C++

      struct DummyCurve
      {
          using ValueType = ...; // this is the type of value 
                                 // which function valueAt returns
      
          using ValueRange = ValueRangeImpl<
              double,
              typename std::less_equal<double>,
              typename std::greater_equal<double>>;
      
          ValueRange getParamValueRange() const;
          ValueType valueAt(double t) const;
      };

2. The next step strategy. Here, I use curiously recurring template pattern which implies less overhead. The strategy must define the inheritor.

The next step strategy must be defined by function:

bool getNextStep(double& du) const {...}

This function returns true if the next step is possible and the value of the next step. If the next step is impossible (or some condition was satisfied), it returns false and iteration stops.

The begin() returns an exemplar of Derived class as Derived knows better how to compare with the end and iterate will factually the exemplar of Derived class.

But the end() may return any class and exemplar of Derived class can conclude when the work is done by returning true in comparison with the end, thus the special TheEndIterator was introduced and some optimization was undertaken: the exemplar of TheEndIterator always the same.

Simplest Working Iterator

Now let's write the simplest working iterator. It will define walking along the curve (parametric function) with const step in the function definition domain:

// Iterator with constant step
//
template<typename Function>
class ConsecutiveIterator : public RangeIterator<Function, ConsecutiveIterator<Function>>
{
    using Inherited = RangeIterator<Function, ConsecutiveIterator<Function>>;
public:

    explicit ConsecutiveIterator(const Function& function, double du)
        : Inherited(function)
        , m_du(du)
    {}
    explicit ConsecutiveIterator(const Function& function, 
                                 const typename Inherited::ValueRange& range, 
                                 const double du)
        : Inherited(function, range)
        , m_du(du)
    {}

    ConsecutiveIterator(const ConsecutiveIterator&) = default;
    ConsecutiveIterator(ConsecutiveIterator&&) = default;

    bool getNextStep(double& du) const
    {
        du = m_du;
        return true;
    }

protected:
    const double m_du = 0; // step along u-domain
};

The function:

bool getNextStep(double& du) const  
{    
     du = m_du;
     return true;     
}

The next step is always possible and its value is constant and defined in the constructor.

Iteration with Tangent Angle Preservation

A more useful example is the iteration with dynamically changing step. Imagine that our purpose is to output curve to view. Often, it is performed by collecting points spread out along the curve, transmitting these points to some drawing engine which will present these points on the screen connected one by one. Every one of us knows it well, do we? It is obvious that it would be sufficient to generate a few points on areas of this low curvature and a lot of it on the high curvature areas of the curve. Iterator with tangent angle preservation will help us in solving such a problem.

This iterator must have a little bit more knowledge about curve to iterate:

1. The curve must provide a function for retrieving the tangent or the first derivative at given t.
2. The curve must provide a function for retrieving the normal or the second derivative at given t.

So the curve interface may look like:

struct DummyCurve
{
    using Point = ...;
    using Vector = ...;
    using ValueType = Point;     // this is the type of value
                                 // which function valueAt returns
    using CoValueType = Vector;  // this is the CoType which operation  
                                 // ValueType - ValueType returns

    using ValueRange = ValueRangeImpl<
        double,
        typename std::less_equal<double>,
        typename std::greater_equal<double>>;

    ValueRange getParamValueRange() const;
    ValueType valueAt(double t) const;

    CoValueType tangentAt(double t) const;
    CoValueType normalAt(double t) const;
};

The ValueRange was described in this tip.

CoValueType was introduced as the type of differential may vary from type of function to differentiate, for example, path and velocity.

Let's imagine that we can present the short part of the curve as part of a circle which radius is reciprocal to curvature of the curve investigated at some point with parameter value t. Basing upon such imagination, we could calculate the next step as:

ds = dAlpha * |(dc/dt)|² / |(dc/dt) ×(d²c/d²t))|

This magic formula was taken from differential geometry and its proof is beyond the scope of this article. In this magic:

1. ds - length of next step
2. dAlpha - the desired angle deviation
3. (dc/dt) - tangent to the curve or it first derivation by t
4. (d²c/d²t)) - normal to the curve or second derivative
5. × - vector cross product

Thus the const angle Iterator may look like:

template<typename Function>
class ConstAngleIterator : public RangeIterator<Function, ConstAngleIterator<Function>>
{
    using Inherited = RangeIterator<Function, ConstAngleIterator<Function>>;
public:

    ConstAngleIterator(const typename Inherited::FunctionType& function, double alpha)
        : Inherited(function)
        , m_alpha(alpha)
    {}
    ConstAngleIterator(const typename Inherited::FunctionType& function, 
                       const typename Inherited::ValueRange& range, double alpha)
        : Inherited(function, range)
        , m_alpha(alpha)
    {}

    ConstAngleIterator(const ConstAngleIterator&) = default;
    ConstAngleIterator(ConstAngleIterator&&) = default;

    bool getNextStep(double& du) const
    {
        const auto tangent = Inherited::m_function.tangentAt(Inherited::m_u);
        const auto normal = Inherited::m_function.normalAt(Inherited::m_u);
        const auto biNormal = tangent.cross(normal);

        const auto biNormalLength = biNormal.length();

        if (isNear(biNormalLength, 0.0))
        {
            // if curve does not bend, let's chose arbitrary small step
            du = Inherited::m_range.length() * 0.05;
        }
        else
        {
            const auto tangentLength = tangent.length();
            du = m_alpha * tangentLength * tangentLength / biNormalLength;
        }

        return true;
    }

protected:
    const double m_alpha = 0.015; // Pi * 0.05
};

This iterator acquires the angle deviation in the constructor and calculates the next step in accordance with it.

Using the Code

Use these iterators whenever you need to iterate over some parametric function:

void consecutiveIterateAlong(const DummyCurve& curve)
{
    // iteration along curve 
    for (DummyPoint pointAtT : ConsecutiveIterator<DummyCurve>(curve, 0.1/*step*/))
    {
        std::cout << pointAtT << std::endl;
    }
}

Compare it with the naive way of iteration:

void consecutiveIterateAlong(const DummyCurve& curve)
{
    // iteration along curve 
    const auto valueRange = curve.getParamValueRange();
    for (double t = valueRange.start, step = 0.1;
         t <= valueRange.stop + valueRange.length() * std::numeric_limits<double>::epsilon();
         t += step)
    {
        DummyPoint pointAtT = curve.valueAt(t);
        std::cout << pointAtT << std::endl;
    }
}

And what about the more difficult way of next step calculation?

History

• 26^th December, 2020: Initial version