Ramer Douglas Peucker Algorithm for Simplification, iterative c++ implementation

Hello there,

After a while, I decided to write this post not to explain RDP algorithm, because I suppose it's widely available over internet e.g.: wikipedia
but since mostly implementations provided are recursive ( that could lead to a stack over flow ), I wanted to provide a C++ iterative implementation I wrote last year in one of my projects.

I will start with a small reminder of what is Ramer-Douglas-Peuker ( Quote from Wikipedia ), then I will provide the implementation ( C++ ) that could be easily translated to Java or other languages.

Definition:

The Ramer–Douglas–Peucker algorithm (RDP) is an algorithm for reducing the number of points in a curve that is approximated by a series of points.

The starting curve is an ordered set of points or lines and the distance dimension ε > 0.
The algorithm recursively divides the line. Initially it is given all the points between the first and last point. It automatically marks the first and last point to be kept. It then finds the point that is furthest from the line segment with the first and last points as end points (this point is obviously furthest on the curve from the approximating line segment between the end points). If the point is closer than ε to the line segment then any points not currently marked to keep can be discarded without the simplified curve being worse than ε.
If the point furthest from the line segment is greater than ε from the approximation then that point must be kept. The algorithm recursively calls itself with the first point and the worst point and then with the worst point and the last point (which includes marking the worst point being marked as kept).
When the recursion is completed a new output curve can be generated consisting of all (and only) those points that have been marked as kept.

 Iterative implementation( C++ ) :

 Let us assume we have a type called point that represents the coordinates of a point in the curve.

Since we want to implement the RDP iteratively, we would need a stack data structure to store the elements.

<code>
    struct StackElement
    {
        point point;
        size_type index;
    };

    std::vector< StackElement > rdp_stack; 

<code>

Now here is the rdp function that we would apply to an input vector of points, with a given threshold i.e. ε . It would return 0 if success, -1 if not.

<code>

 int
rdp_implimentation( const std::vector& input,
                                  std::vector& output,
                                  double ε )
{
    output.resize( 0 );

    if ( input.empty( ) )
    {
        return 0;
    }

    // Clean up
    rdp_stack.resize( 0 );

    if (rdp_stack.reserve( input.size( ) ) < 0
        || output.reserve( input.size( ) ) < 0 )
    {
        // Failed to allocate memory
        return -1;
    }

    point anchor = input.front( );
    size_t index_of_anchor = 0;
    point floater = input.back( );
    size_t index_of_floater = input.size( ) - 1;

    // Add the first point in the poly to the result
    output.push_back( anchor );

    StackElement stack_element =
    {
        floater,
        index_of_floater
    };

    rdp_stack.push_back( stack_element );

    while ( !rdp_stack.empty( ) )
    {
        double max_squared_distance = 0.0;
        point farthest = anchor;
        size_t index_of_farthest = index_of_anchor;

        // Find point furthest from line defined by anchor and floater

        // function depends on your projection/dimension
        DistanceHelper distance_helper( anchor, floater );

        for ( size_t i = index_of_anchor + 1; i < index_of_floater; ++i )
        {

            const double squared_distance = distance_helper.squared_distance_to( input[ i ] );
            if ( squared_distance > max_squared_distance )
            {
                max_squared_distance = squared_distance;
                farthest = input[ i ];
                index_of_farthest = i;
            }
        }

        // Furthest point nearer than tolerance?
        if ( max_squared_distance <= ε )
        {
            output.push_back( rdp_stack.back( ).point );
            rdp_stack.pop_back( );
            anchor = floater;
            index_of_anchor = index_of_floater;
            if ( !rdp_stack.empty( ) )
            {
                floater = rdp_stack.back( ).point;
                index_of_floater = rdp_stack.back( ).index;
            }
        }
        else
        {
            floater = farthest;
            index_of_floater = index_of_farthest;
            stack_element.point = floater;
            stack_element.index = index_of_floater;
            rdp_stack.push_back( stack_element );
        }
    }

    // Success
    return 0;
}

<code>

Useful links :

- Matlab demo

- Code project

Conclusion

RDP could be useful in many use cases such us in preparing and transforming polylines before rendering them, this is so called line simplification. It can enhance the rendering performance and reduce GPU work.

I'm blogging today from Chicago, next post will come soon hopefully from  home in Berlin/Germany. I hope you enjoy the post and hope to hear back from you :)
God bless internet!

Cheers
Anis