psimpl.h

Go to the documentation of this file.

/* ***** BEGIN LICENSE BLOCK *****
 * Version: MPL 1.1
 *
 * The contents of this file are subject to the Mozilla Public License Version
 * 1.1 (the "License"); you may not use this file except in compliance with
 * the License. You may obtain a copy of the License at
 * http://www.mozilla.org/MPL/
 *
 * Software distributed under the License is distributed on an "AS IS" basis,
 * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
 * for the specific language governing rights and limitations under the
 * License.
 *
 * The Original Code is
 * 'psimpl - generic n-dimensional polyline simplification'.
 *
 * The Initial Developer of the Original Code is
 * Elmar de Koning.
 * Portions created by the Initial Developer are Copyright (C) 2010-2011
 * the Initial Developer. All Rights Reserved.
 *
 * Contributor(s):
 *
 * ***** END LICENSE BLOCK ***** */

/*
    psimpl - generic n-dimensional polyline simplification
    Copyright (C) 2010-2011 Elmar de Koning, edekoning@gmail.com

    This file is part of psimpl, and is hosted at SourceForge:
    http://sourceforge.net/projects/psimpl/
*/

#ifndef PSIMPL_GENERIC
#define PSIMPL_GENERIC


#include <queue>
#include <stack>
#include <numeric>
#include <algorithm>
#include <cmath>


namespace psimpl
{
    namespace util
    {
        template <typename T>
        class scoped_array
        {
        public:
            scoped_array (unsigned n) {
                array = new T [n];
            }

            ~scoped_array () {
                delete [] array;
            }

            T& operator [] (int offset) {
                return array [offset];
            }

            const T& operator [] (int offset) const {
                return array [offset];
            }

            T* get () const {
                return array;
            }

            void swap (scoped_array& b) {
                T* tmp = b.array;
                b.array = array;
                array = tmp;
            }

        private:
            scoped_array (const scoped_array&);
            scoped_array& operator= (const scoped_array&);

        private:
            T* array;
        };

        template <typename T> inline void swap (scoped_array <T>& a, scoped_array <T>& b) {
            a.swap (b);
        }
    }

    namespace math
    {
        struct Statistics
        {
            Statistics () :
                max (0),
                sum (0),
                mean (0),
                std (0)
            {}

            double max;
            double sum;
            double mean;
            double std;     
        };

        template <unsigned DIM, class InputIterator>
        inline bool equal (
            InputIterator p1,
            InputIterator p2)
        {
            for (unsigned d = 0; d < DIM; ++d) {
                if (*p1 != *p2) {
                    return false;
                }
                ++p1;
                ++p2;
            }
            return true;
        }

        template <unsigned DIM, class InputIterator, class OutputIterator>
        inline OutputIterator make_vector (
            InputIterator p1,
            InputIterator p2,
            OutputIterator result)
        {
            for (unsigned d = 0; d < DIM; ++d) {
                *result = *p2 - *p1;
                ++result;
                ++p1;
                ++p2;
            }
            return result;
        }

        template <unsigned DIM, class InputIterator>
        inline typename std::iterator_traits <InputIterator>::value_type dot (
            InputIterator v1,
            InputIterator v2)
        {
            typename std::iterator_traits <InputIterator>::value_type result = 0;
            for (unsigned d = 0; d < DIM; ++d) {
                result += (*v1) * (*v2);
                ++v1;
                ++v2;
            }
            return result;
        }

        template <unsigned DIM, class InputIterator, class OutputIterator>
        inline OutputIterator interpolate (
            InputIterator p1,
            InputIterator p2,
            float fraction,
            OutputIterator result)
        {
            typedef typename std::iterator_traits <InputIterator>::value_type value_type;

            for (unsigned d = 0; d < DIM; ++d) {
                *result = *p1 + static_cast <value_type> (fraction * (*p2 - *p1));
                ++result;
                ++p1;
                ++p2;
            }
            return result;
        }

        template <unsigned DIM, class InputIterator1, class InputIterator2>
        inline typename std::iterator_traits <InputIterator1>::value_type point_distance2 (
            InputIterator1 p1,
            InputIterator2 p2)
        {
            typename std::iterator_traits <InputIterator1>::value_type result = 0;
            for (unsigned d = 0; d < DIM; ++d) {
                result += (*p1 - *p2) * (*p1 - *p2);
                ++p1;
                ++p2;
            }
            return result;
        }

        template <unsigned DIM, class InputIterator>
        inline typename std::iterator_traits <InputIterator>::value_type line_distance2 (
            InputIterator l1,
            InputIterator l2,
            InputIterator p)
        {
            typedef typename std::iterator_traits <InputIterator>::value_type value_type;

            value_type v [DIM];                 // vector l1 --> l2
            value_type w [DIM];                 // vector l1 --> p

            make_vector <DIM> (l1, l2, v);
            make_vector <DIM> (l1, p,  w);

            value_type cv = dot <DIM> (v, v);   // squared length of v
            value_type cw = dot <DIM> (w, v);   // project w onto v

            // avoid problems with divisions when value_type is an integer type
            float fraction = cv == 0 ? 0 : static_cast <float> (cw) / static_cast <float> (cv);

            value_type proj [DIM];              // p projected onto line (l1, l2)
            interpolate <DIM> (l1, l2, fraction, proj);

            return point_distance2 <DIM> (p, proj);
        }

        template <unsigned DIM, class InputIterator>
        inline typename std::iterator_traits <InputIterator>::value_type segment_distance2 (
            InputIterator s1,
            InputIterator s2,
            InputIterator p)
        {
            typedef typename std::iterator_traits <InputIterator>::value_type value_type;

            value_type v [DIM];        // vector s1 --> s2
            value_type w [DIM];        // vector s1 --> p

            make_vector <DIM> (s1, s2, v);
            make_vector <DIM> (s1, p,  w);

            value_type cw = dot <DIM> (w, v);   // project w onto v
            if (cw <= 0) {
                // projection of w lies to the left of s1
                return point_distance2 <DIM> (p, s1);
            }

            value_type cv = dot <DIM> (v, v);   // squared length of v
            if (cv <= cw) {
                // projection of w lies to the right of s2
                return point_distance2 <DIM> (p, s2);
            }

            // avoid problems with divisions when value_type is an integer type
            float fraction = cv == 0 ? 0 : static_cast <float> (cw) / static_cast <float> (cv);

            value_type proj [DIM];    // p projected onto segement (s1, s2)
            interpolate <DIM> (s1, s2, fraction, proj);

            return point_distance2 <DIM> (p, proj);
        }

        template <unsigned DIM, class InputIterator>
        inline typename std::iterator_traits <InputIterator>::value_type ray_distance2 (
            InputIterator r1,
            InputIterator r2,
            InputIterator p)
        {
            typedef typename std::iterator_traits <InputIterator>::value_type value_type;

            value_type v [DIM];        // vector r1 --> r2
            value_type w [DIM];        // vector r1 --> p

            make_vector <DIM> (r1, r2, v);
            make_vector <DIM> (r1, p,  w);

            value_type cv = dot <DIM> (v, v);    // squared length of v
            value_type cw = dot <DIM> (w, v);    // project w onto v

            if (cw <= 0) {
                // projection of w lies to the left of r1 (not on the ray)
                return point_distance2 <DIM> (p, r1);
            }

            // avoid problems with divisions when value_type is an integer type
            float fraction = cv == 0 ? 0 : static_cast <float> (cw) / static_cast <float> (cv);

            value_type proj [DIM];    // p projected onto ray (r1, r2)
            interpolate <DIM> (r1, r2, fraction, proj);

            return point_distance2 <DIM> (p, proj);
        }

        template <class InputIterator>
        inline Statistics compute_statistics (
            InputIterator first,
            InputIterator last)
        {
            typedef typename std::iterator_traits <InputIterator>::value_type value_type;
            typedef typename std::iterator_traits <InputIterator>::difference_type diff_type;

            Statistics stats;

            diff_type count = std::distance (first, last);
            if (count == 0) {
                return stats;
            }

            value_type init = 0;
            stats.max = static_cast <double> (*std::max_element (first, last));
            stats.sum = static_cast <double> (std::accumulate (first, last, init));
            stats.mean = stats.sum / count;
            std::transform (first, last, first, std::bind2nd (std::minus <value_type> (), stats.mean));
            stats.std = std::sqrt (static_cast <double> (std::inner_product (first, last, first, init)) / count);
            return stats;
        }
    }

    template <unsigned DIM, class InputIterator, class OutputIterator>
    class PolylineSimplification
    {
        typedef typename std::iterator_traits <InputIterator>::difference_type diff_type;
        typedef typename std::iterator_traits <InputIterator>::value_type value_type;
        typedef typename std::iterator_traits <const value_type*>::difference_type ptr_diff_type;

    public:
        OutputIterator NthPoint (
            InputIterator first,
            InputIterator last,
            unsigned n,
            OutputIterator result)
        {
            diff_type coordCount = std::distance (first, last);
            diff_type pointCount = DIM              // protect against zero DIM
                                   ? coordCount / DIM
                                   : 0;

            // validate input and check if simplification required
            if (coordCount % DIM || pointCount < 3 || n < 2) {
                return std::copy (first, last, result);
            }

            unsigned remaining = pointCount - 1;    // the number of points remaining after key
            InputIterator key = first;              // indicates the current key

            // the first point is always part of the simplification
            CopyKey (key, result);

            // copy each nth point
            while (Forward (key, n, remaining)) {
                CopyKey (key, result);
            }

            return result;
        }

        OutputIterator RadialDistance (
            InputIterator first,
            InputIterator last,
            value_type tol,
            OutputIterator result)
        {
            diff_type coordCount = std::distance (first, last);
            diff_type pointCount = DIM      // protect against zero DIM
                                   ? coordCount / DIM
                                   : 0;
            value_type tol2 = tol * tol;    // squared distance tolerance

            // validate input and check if simplification required
            if (coordCount % DIM || pointCount < 3 || tol2 == 0) {
                return std::copy (first, last, result);
            }

            InputIterator current = first;  // indicates the current key
            InputIterator next = first;     // used to find the next key

            // the first point is always part of the simplification
            CopyKeyAdvance (next, result);

            // Skip first and last point, because they are always part of the simplification
            for (diff_type index = 1; index < pointCount - 1; ++index) {
                if (math::point_distance2 <DIM> (current, next) < tol2) {
                    Advance (next);
                    continue;
                }
                current = next;
                CopyKeyAdvance (next, result);
            }
            // the last point is always part of the simplification
            CopyKeyAdvance (next, result);

            return result;
        }

        OutputIterator PerpendicularDistance (
            InputIterator first,
            InputIterator last,
            value_type tol,
            unsigned repeat,
            OutputIterator result)
        {
            if (repeat == 1) {
                // single pass
                return PerpendicularDistance (first, last, tol, result);
            }
            // only validate repeat; other input is validated by simplify_perpendicular_distance
            if (repeat < 1) {
                return std::copy (first, last, result);
            }
            diff_type coordCount = std::distance (first, last);

            // first pass: [first, last) --> temporary array 'tempPoly'
            util::scoped_array <value_type> tempPoly (coordCount);
            PolylineSimplification <DIM, InputIterator, value_type*> psimpl_to_array;
            diff_type tempCoordCount = std::distance (tempPoly.get (),
                psimpl_to_array.PerpendicularDistance (first, last, tol, tempPoly.get ()));

            // check if simplification did not improved
            if (coordCount == tempCoordCount) {
                return std::copy (tempPoly.get (), tempPoly.get () + coordCount, result);
            }
            std::swap (coordCount, tempCoordCount);
            --repeat;

            // intermediate passes: temporary array 'tempPoly' --> temporary array 'tempResult'
            if (1 < repeat) {
                util::scoped_array <value_type> tempResult (coordCount);
                PolylineSimplification <DIM, value_type*, value_type*> psimpl_arrays;

                while (--repeat) {
                    tempCoordCount = std::distance (tempResult.get (),
                        psimpl_arrays.PerpendicularDistance (
                            tempPoly.get (), tempPoly.get () + coordCount, tol, tempResult.get ()));

                    // check if simplification did not improved
                    if (coordCount == tempCoordCount) {
                        return std::copy (tempPoly.get (), tempPoly.get () + coordCount, result);
                    }
                    util::swap (tempPoly, tempResult);
                    std::swap (coordCount, tempCoordCount);
                }
            }

            // final pass: temporary array 'tempPoly' --> result
            PolylineSimplification <DIM, value_type*, OutputIterator> psimpl_from_array;
            return psimpl_from_array.PerpendicularDistance (
                tempPoly.get (), tempPoly.get () + coordCount, tol, result);
        }

        OutputIterator PerpendicularDistance (
            InputIterator first,
            InputIterator last,
            value_type tol,
            OutputIterator result)
        {
            diff_type coordCount = std::distance (first, last);
            diff_type pointCount = DIM      // protect against zero DIM
                                   ? coordCount / DIM
                                   : 0;
            value_type tol2 = tol * tol;    // squared distance tolerance

            // validate input and check if simplification required
            if (coordCount % DIM || pointCount < 3 || tol2 == 0) {
                return std::copy (first, last, result);
            }

            InputIterator p0 = first;
            InputIterator p1 = AdvanceCopy(p0);
            InputIterator p2 = AdvanceCopy(p1);

            // the first point is always part of the simplification
            CopyKey (p0, result);

            while (p2 != last) {
                // test p1 against line segment S(p0, p2)
                if (math::segment_distance2 <DIM> (p0, p2, p1) < tol2) {
                    CopyKey (p2, result);
                    // move up by two points
                    p0 = p2;
                    Advance (p1, 2);
                    if (p1 == last) {
                        // protect against advancing p2 beyond last
                        break;
                    }
                    Advance (p2, 2);
                }
                else {
                    CopyKey (p1, result);
                    // move up by one point
                    p0 = p1;
                    p1 = p2;
                    Advance (p2);
                }
            }
            // make sure the last point is part of the simplification
            if (p1 != last) {
                CopyKey (p1, result);
            }
            return result;
        }

        OutputIterator ReumannWitkam (
            InputIterator first,
            InputIterator last,
            value_type tol,
            OutputIterator result)
        {
            diff_type coordCount = std::distance (first, last);
            diff_type pointCount = DIM      // protect against zero DIM
                                   ? coordCount / DIM
                                   : 0;
            value_type tol2 = tol * tol;    // squared distance tolerance

            // validate input and check if simplification required
            if (coordCount % DIM || pointCount < 3 || tol2 == 0) {
                return std::copy (first, last, result);
            }

            // define the line L(p0, p1)
            InputIterator p0 = first;               // indicates the current key
            InputIterator p1 = AdvanceCopy (first); // indicates the next point after p0

            // keep track of two test points
            InputIterator pi = p1;     // the previous test point
            InputIterator pj = p1;     // the current test point (pi+1)

            // the first point is always part of the simplification
            CopyKey (p0, result);

            // check each point pj against L(p0, p1)
            for (diff_type j = 2; j < pointCount; ++j) {
                pi = pj;
                Advance (pj);

                if (math::line_distance2 <DIM> (p0, p1, pj) < tol2) {
                    continue;
                }
                // found the next key at pi
                CopyKey (pi, result);
                // define new line L(pi, pj)
                p0 = pi;
                p1 = pj;
            }
            // the last point is always part of the simplification
            CopyKey (pj, result);

            return result;
        }

        OutputIterator Opheim (
            InputIterator first,
            InputIterator last,
            value_type min_tol,
            value_type max_tol,
            OutputIterator result)
        {
            diff_type coordCount = std::distance (first, last);
            diff_type pointCount = DIM                    // protect against zero DIM
                                   ? coordCount / DIM
                                   : 0;
            value_type min_tol2 = min_tol * min_tol;    // squared minimum distance tolerance
            value_type max_tol2 = max_tol * max_tol;    // squared maximum distance tolerance

            // validate input and check if simplification required
            if (coordCount % DIM || pointCount < 3 || min_tol2 == 0 || max_tol2 == 0) {
                return std::copy (first, last, result);
            }

            // define the ray R(r0, r1)
            InputIterator r0 = first;  // indicates the current key and start of the ray
            InputIterator r1 = first;  // indicates a point on the ray
            bool rayDefined = false;

            // keep track of two test points
            InputIterator pi = r0;     // the previous test point
            InputIterator pj =         // the current test point (pi+1)
                AdvanceCopy (pi);

            // the first point is always part of the simplification
            CopyKey (r0, result);

            for (diff_type j = 2; j < pointCount; ++j) {
                pi = pj;
                Advance (pj);

                if (!rayDefined) {
                    // discard each point within minimum tolerance
                    if (math::point_distance2 <DIM> (r0, pj) < min_tol2) {
                        continue;
                    }
                    // the last point within minimum tolerance pi defines the ray R(r0, r1)
                    r1 = pi;
                    rayDefined = true;
                }

                // check each point pj against R(r0, r1)
                if (math::point_distance2 <DIM> (r0, pj) < max_tol2 &&
                    math::ray_distance2 <DIM> (r0, r1, pj) < min_tol2)
                {
                    continue;
                }
                // found the next key at pi
                CopyKey (pi, result);
                // define new ray R(pi, pj)
                r0 = pi;
                rayDefined = false;
            }
            // the last point is always part of the simplification
            CopyKey (pj, result);

            return result;
        }

        OutputIterator Lang (
            InputIterator first,
            InputIterator last,
            value_type tol,
            unsigned look_ahead,
            OutputIterator result)
        {
            diff_type coordCount = std::distance (first, last);
            diff_type pointCount = DIM      // protect against zero DIM
                                   ? coordCount / DIM
                                   : 0;
            value_type tol2 = tol * tol;    // squared minimum distance tolerance
            
            // validate input and check if simplification required
            if (coordCount % DIM || pointCount < 3 || look_ahead < 2 || tol2 == 0) {
                return std::copy (first, last, result);
            }

            InputIterator current = first;          // indicates the current key
            InputIterator next = first;             // used to find the next key

            unsigned remaining = pointCount - 1;    // the number of points remaining after current
            unsigned moved = Forward (next, look_ahead, remaining);

            // the first point is always part of the simplification
            CopyKey (current, result);

            while (moved) {
                value_type d2 = 0;
                InputIterator p = AdvanceCopy (current);

                while (p != next) {
                    d2 = std::max (d2, math::segment_distance2 <DIM> (current, next, p));
                    if (tol2 < d2) {
                        break;
                    }
                    Advance (p);
                }
                if (d2 < tol2) {
                    current = next;
                    CopyKey (current, result);
                    moved = Forward (next, look_ahead, remaining);
                }
                else {
                    Backward (next, remaining);
                }
            }
            return result;
        }

        OutputIterator DouglasPeucker (
            InputIterator first,
            InputIterator last,
            value_type tol,
            OutputIterator result)
        {
            diff_type coordCount = std::distance (first, last);
            diff_type pointCount = DIM      // protect against zero DIM
                                   ? coordCount / DIM
                                   : 0;
            // validate input and check if simplification required
            if (coordCount % DIM || pointCount < 3 || tol == 0) {
                return std::copy (first, last, result);
            }
            // radial distance routine as preprocessing
            util::scoped_array <value_type> reduced (coordCount);   // radial distance results
            PolylineSimplification <DIM, InputIterator, value_type*> psimpl_to_array;
            ptr_diff_type reducedCoordCount = std::distance (reduced.get (),
                psimpl_to_array.RadialDistance (first, last, tol, reduced.get ()));
            ptr_diff_type reducedPointCount = reducedCoordCount / DIM;

            // douglas-peucker approximation
            util::scoped_array <unsigned char> keys (pointCount);         // douglas-peucker results
            DPHelper::Approximate (reduced.get (), reducedCoordCount, tol, keys.get ());

            // copy all keys
            const value_type* current = reduced.get ();
            for (ptr_diff_type p=0; p<reducedPointCount; ++p, current += DIM) {
                if (keys [p]) {
                    for (unsigned d = 0; d < DIM; ++d) {
                        *result = current [d];
                        ++result;
                    }
                }
            }
            return result;
        }

        OutputIterator DouglasPeuckerN (
            InputIterator first,
            InputIterator last,
            unsigned count,
            OutputIterator result)
        {
            diff_type coordCount = std::distance (first, last);
            diff_type pointCount = DIM      // protect against zero DIM
                                   ? coordCount / DIM
                                   : 0;
            // validate input and check if simplification required
            if (coordCount % DIM || pointCount <= static_cast <diff_type> (count) || count < 2) {
                return std::copy (first, last, result);
            }

            // copy coords
            util::scoped_array <value_type> coords (coordCount);
            for (ptr_diff_type c=0; c<coordCount; ++c) {
                coords [c] = *first;
                ++first;
            }

            // douglas-peucker approximation
            util::scoped_array <unsigned char> keys (pointCount);
            DPHelper::ApproximateN (coords.get (), coordCount, count, keys.get ());

            // copy keys
            const value_type* current = coords.get ();
            for (ptr_diff_type p=0; p<pointCount; ++p, current += DIM) {
                if (keys [p]) {
                    for (unsigned d = 0; d < DIM; ++d) {
                        *result = current [d];
                        ++result;
                    }
                }
            }
            return result;
        }

        OutputIterator ComputePositionalErrors2 (
            InputIterator original_first,
            InputIterator original_last,
            InputIterator simplified_first,
            InputIterator simplified_last,
            OutputIterator result,
            bool* valid=0)
        {
            diff_type original_coordCount = std::distance (original_first, original_last);
            diff_type original_pointCount = DIM     // protect against zero DIM
                                            ? original_coordCount / DIM
                                            : 0;

            diff_type simplified_coordCount = std::distance (simplified_first, simplified_last);
            diff_type simplified_pointCount = DIM   // protect against zero DIM
                                              ? simplified_coordCount / DIM
                                              : 0;

            // validate input
            if (original_coordCount % DIM || original_pointCount < 2 ||
                simplified_coordCount % DIM || simplified_pointCount < 2 ||
                original_pointCount < simplified_pointCount ||
                !math::equal <DIM> (original_first, simplified_first))
            {
                if (valid) {
                    *valid = false;
                }
                return result;
            }

            // define (simplified) line segment S(simplified_prev, simplified_first)
            InputIterator simplified_prev = simplified_first;
            std::advance (simplified_first, DIM);

            // process each simplified line segment
            while (simplified_first != simplified_last) {
                // process each original point until it equals the end of the line segment
                while (original_first != original_last &&
                       !math::equal <DIM> (original_first, simplified_first))
                {
                    *result = math::segment_distance2 <DIM> (simplified_prev, simplified_first,
                                                             original_first);
                    ++result;
                    std::advance (original_first, DIM);
                }
                // update line segment S
                simplified_prev = simplified_first;
                std::advance (simplified_first, DIM);
            }
            // check if last original point matched
            if (original_first != original_last) {
                *result = 0;
                ++result;
            }

            if (valid) {
                *valid = original_first != original_last;
            }
            return result;
        }

        math::Statistics ComputePositionalErrorStatistics (
            InputIterator original_first,
            InputIterator original_last,
            InputIterator simplified_first,
            InputIterator simplified_last,
            bool* valid=0)
        {
            diff_type pointCount = std::distance (original_first, original_last) / DIM;
            util::scoped_array <double> errors (pointCount);
            PolylineSimplification <DIM, InputIterator, double*> ps;

            diff_type errorCount = 
                std::distance (
                    errors.get (), 
                    ps.ComputePositionalErrors2 (original_first, original_last,
                                                 simplified_first, simplified_last,
                                                 errors.get (), valid));

            std::transform (errors.get (), errors.get () + errorCount,
                            errors.get (),
                            std::ptr_fun <double, double> (std::sqrt));

            return math::compute_statistics (errors.get (), errors.get () + errorCount);
        }

    private:
        inline void CopyKeyAdvance (
            InputIterator& key,
            OutputIterator& result)
        {
            for (unsigned d = 0; d < DIM; ++d) {
                *result = *key;
                ++result;
                ++key;
            }
        }

        inline void CopyKey (
            InputIterator key,
            OutputIterator& result)
        {
            CopyKeyAdvance (key, result);
        }

        inline void Advance (
            InputIterator& it,
            diff_type n = 1)
        {
            std::advance (it, n * static_cast <diff_type> (DIM));
        }
        
        inline InputIterator AdvanceCopy (
            InputIterator it,
            diff_type n = 1)
        {
            Advance (it, n);
            return it;
        }

        inline unsigned Forward (
            InputIterator& it,
            unsigned n,
            unsigned& remaining)
        {
            n = std::min (n, remaining);
            Advance (it, n);
            remaining -= n;
            return n;
        }

        inline void Backward (
            InputIterator& it,
            unsigned& remaining)
        {
            Advance (it, -1);
            ++remaining;
        }

    private:
        class DPHelper
        {
            struct SubPoly {
                SubPoly (ptr_diff_type first=0, ptr_diff_type last=0) :
                    first (first), last (last) {}

                ptr_diff_type first;    
                ptr_diff_type last;     
            };

            struct KeyInfo {
                KeyInfo (ptr_diff_type index=0, value_type dist2=0) :
                    index (index), dist2 (dist2) {}

                ptr_diff_type index;    
                value_type dist2;       
            };

            struct SubPolyAlt {
                SubPolyAlt (ptr_diff_type first=0, ptr_diff_type last=0) :
                    first (first), last (last) {}

                ptr_diff_type first;    
                ptr_diff_type last;     
                KeyInfo keyInfo;        

                bool operator< (const SubPolyAlt& other) const {
                    return keyInfo.dist2 < other.keyInfo.dist2;
                }
            };

        public:
            static void Approximate (
                const value_type* coords,
                ptr_diff_type coordCount,
                value_type tol,
                unsigned char* keys)
            {
                value_type tol2 = tol * tol;    // squared distance tolerance
                ptr_diff_type pointCount = coordCount / DIM;
                // zero out keys
                std::fill_n (keys, pointCount, 0);
                keys [0] = 1;                   // the first point is always a key
                keys [pointCount - 1] = 1;      // the last point is always a key

                typedef std::stack <SubPoly> Stack;
                Stack stack;                    // LIFO job-queue containing sub-polylines

                SubPoly subPoly (0, coordCount-DIM);
                stack.push (subPoly);           // add complete poly

                while (!stack.empty ()) {
                    subPoly = stack.top ();     // take a sub poly
                    stack.pop ();               // and find its key
                    KeyInfo keyInfo = FindKey (coords, subPoly.first, subPoly.last);
                    if (keyInfo.index && tol2 < keyInfo.dist2) {
                        // store the key if valid
                        keys [keyInfo.index / DIM] = 1;
                        // split the polyline at the key and recurse
                        stack.push (SubPoly (keyInfo.index, subPoly.last));
                        stack.push (SubPoly (subPoly.first, keyInfo.index));
                    }
                }
            }

            static void ApproximateN (
                const value_type* coords,
                ptr_diff_type coordCount,
                unsigned countTol,
                unsigned char* keys)
            {
                ptr_diff_type pointCount = coordCount / DIM;
                // zero out keys
                std::fill_n (keys, pointCount, 0);
                keys [0] = 1;                   // the first point is always a key
                keys [pointCount - 1] = 1;      // the last point is always a key
                unsigned keyCount = 2;

                if (countTol == 2) {
                    return;
                }

                typedef std::priority_queue <SubPolyAlt> PriorityQueue;
                PriorityQueue queue;    // sorted (max dist2) job queue containing sub-polylines

                SubPolyAlt subPoly (0, coordCount-DIM);
                subPoly.keyInfo = FindKey (coords, subPoly.first, subPoly.last);
                queue.push (subPoly);           // add complete poly

                while (!queue.empty ()) {
                    subPoly = queue.top ();     // take a sub poly
                    queue.pop ();
                    // store the key
                    keys [subPoly.keyInfo.index / DIM] = 1;
                    // check point count tolerance
                    keyCount++;
                    if (keyCount == countTol) {
                        break;
                    }
                    // split the polyline at the key and recurse
                    SubPolyAlt left (subPoly.first, subPoly.keyInfo.index);
                    left.keyInfo = FindKey (coords, left.first, left.last);
                    if (left.keyInfo.index) {
                        queue.push (left);
                    }
                    SubPolyAlt right (subPoly.keyInfo.index, subPoly.last);
                    right.keyInfo = FindKey (coords, right.first, right.last);
                    if (right.keyInfo.index) {
                        queue.push (right);
                    }
                }
            }

        private:
            static KeyInfo FindKey (
                const value_type* coords,
                ptr_diff_type first,
                ptr_diff_type last)
            {
                KeyInfo keyInfo;

                for (ptr_diff_type current = first + DIM; current < last; current += DIM) {
                    value_type d2 = math::segment_distance2 <DIM> (coords + first, coords + last,
                                                                   coords + current);
                    if (d2 < keyInfo.dist2) {
                        continue;
                    }
                    // update maximum squared distance and the point it belongs to
                    keyInfo.index = current;
                    keyInfo.dist2 = d2;
                }
                return keyInfo;
            }
        };
    };

    template <unsigned DIM, class ForwardIterator, class OutputIterator>
    OutputIterator simplify_nth_point (
        ForwardIterator first,
        ForwardIterator last,
        unsigned n,
        OutputIterator result)
    {
        PolylineSimplification <DIM, ForwardIterator, OutputIterator> ps;
        return ps.NthPoint (first, last, n, result);
    }

    template <unsigned DIM, class ForwardIterator, class OutputIterator>
    OutputIterator simplify_radial_distance (
        ForwardIterator first,
        ForwardIterator last,
        typename std::iterator_traits <ForwardIterator>::value_type tol,
        OutputIterator result)
    {
        PolylineSimplification <DIM, ForwardIterator, OutputIterator> ps;
        return ps.RadialDistance (first, last, tol, result);
    }

    template <unsigned DIM, class ForwardIterator, class OutputIterator>
    OutputIterator simplify_perpendicular_distance (
        ForwardIterator first,
        ForwardIterator last,
        typename std::iterator_traits <ForwardIterator>::value_type tol,
        unsigned repeat,
        OutputIterator result)
    {
        PolylineSimplification <DIM, ForwardIterator, OutputIterator> ps;
        return ps.PerpendicularDistance (first, last, tol, repeat, result);
    }

    template <unsigned DIM, class ForwardIterator, class OutputIterator>
    OutputIterator simplify_perpendicular_distance (
        ForwardIterator first,
        ForwardIterator last,
        typename std::iterator_traits <ForwardIterator>::value_type tol,
        OutputIterator result)
    {
        PolylineSimplification <DIM, ForwardIterator, OutputIterator> ps;
        return ps.PerpendicularDistance (first, last, tol, result);
    }

    template <unsigned DIM, class ForwardIterator, class OutputIterator>
    OutputIterator simplify_reumann_witkam (
        ForwardIterator first,
        ForwardIterator last,
        typename std::iterator_traits <ForwardIterator>::value_type tol,
        OutputIterator result)
    {
        PolylineSimplification <DIM, ForwardIterator, OutputIterator> ps;
        return ps.ReumannWitkam (first, last, tol, result);
    }

    template <unsigned DIM, class ForwardIterator, class OutputIterator>
    OutputIterator simplify_opheim (
        ForwardIterator first,
        ForwardIterator last,
        typename std::iterator_traits <ForwardIterator>::value_type min_tol,
        typename std::iterator_traits <ForwardIterator>::value_type max_tol,
        OutputIterator result)
    {
        PolylineSimplification <DIM, ForwardIterator, OutputIterator> ps;
        return ps.Opheim (first, last, min_tol, max_tol, result);
    }

    template <unsigned DIM, class BidirectionalIterator, class OutputIterator>
    OutputIterator simplify_lang (
        BidirectionalIterator first,
        BidirectionalIterator last,
        typename std::iterator_traits <BidirectionalIterator>::value_type tol,
        unsigned look_ahead,
        OutputIterator result)
    {
        PolylineSimplification <DIM, BidirectionalIterator, OutputIterator> ps;
        return ps.Lang (first, last, tol, look_ahead, result);
    }

    template <unsigned DIM, class ForwardIterator, class OutputIterator>
    OutputIterator simplify_douglas_peucker (
        ForwardIterator first,
        ForwardIterator last,
        typename std::iterator_traits <ForwardIterator>::value_type tol,
        OutputIterator result)
    {
        PolylineSimplification <DIM, ForwardIterator, OutputIterator> ps;
        return ps.DouglasPeucker (first, last, tol, result);
    }

    template <unsigned DIM, class ForwardIterator, class OutputIterator>
    OutputIterator simplify_douglas_peucker_n (
        ForwardIterator first,
        ForwardIterator last,
        unsigned count,
        OutputIterator result)
    {
        PolylineSimplification <DIM, ForwardIterator, OutputIterator> ps;
        return ps.DouglasPeuckerN (first, last, count, result);
    }

    template <unsigned DIM, class ForwardIterator, class OutputIterator>
    OutputIterator compute_positional_errors2 (
        ForwardIterator original_first,
        ForwardIterator original_last,
        ForwardIterator simplified_first,
        ForwardIterator simplified_last,
        OutputIterator result,
        bool* valid=0)
    {
        PolylineSimplification <DIM, ForwardIterator, OutputIterator> ps;
        return ps.ComputePositionalErrors2 (original_first, original_last, simplified_first, simplified_last, result, valid);
    }

    template <unsigned DIM, class ForwardIterator>
    math::Statistics compute_positional_error_statistics (
        ForwardIterator original_first,
        ForwardIterator original_last,
        ForwardIterator simplified_first,
        ForwardIterator simplified_last,
        bool* valid=0)
    {
        PolylineSimplification <DIM, ForwardIterator, ForwardIterator> ps;
        return ps.ComputePositionalErrorStatistics (original_first, original_last, simplified_first, simplified_last, valid);
    }
}

#endif // PSIMPL_GENERIC