fhirschmann / rdp

Python/Numpy implementation of the Ramer-Douglas-Peucker algorithm
https://pypi.python.org/pypi/rdp
MIT License
244 stars 117 forks source link

Wrong calculation of pldist #13

Open flbbronnec opened 4 years ago

flbbronnec commented 4 years ago

The current pldist function returns the distance between a point and the infinte line generated by two other points (i.e the distance to the orthogonal projection), which is incorrect. It should return the distance between a point and a line-segment.

Here is a fix coming from this answer https://stackoverflow.com/a/54442561:

def pldist(point, start, end):
    """
    Calculates the distance from ``point`` to the line given
    by the points ``start`` and ``end``.

    :param point: a point
    :type point: numpy array
    :param start: a point of the line
    :type start: numpy array
    :param end: another point of the line
    :type end: numpy array
    """
    if np.all(start == end)):
        return np.linalg.norm(point - start)

    # normalized tangent vector
    d = np.divide(end - start, np.linalg.norm(end - start))

    # signed parallel distance components
    s = np.dot(start - point, d)
    t = np.dot(point - end, d)

    # clamped parallel distance
    h = np.max([s, t, 0])

    # perpendicular distance component, as before
    c = np.cross(point - start, d)

    # use hypot for Pythagoras to improve accuracy
    return np.hypot(h, c)
district10 commented 2 years ago

pybind11-rdp handles this pretty well:

https://github.com/cubao/pybind11-rdp/blob/master/src/main.cpp#L26-L37

struct LineSegment
{
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
    const Eigen::Vector3d A, B, AB;
    const double len2, inv_len2;
    LineSegment(const Eigen::Vector3d &a, const Eigen::Vector3d &b)
        : A(a), B(b), AB(b - a), //
          len2((b - a).squaredNorm()), inv_len2(1.0 / len2)
    {
    }
    double distance2(const Eigen::Vector3d &P) const
    {
        double dot = (P - A).dot(AB);
        if (dot <= 0) {
            return (P - A).squaredNorm();
        } else if (dot >= len2) {
            return (P - B).squaredNorm();
        }
        // P' = A + dot/length * normed(AB)
        //    = A + dot * AB / (length^2)
        return (A + (dot * inv_len2 * AB) - P).squaredNorm();
    }
    double distance(const Eigen::Vector3d &P) const
    {
        return std::sqrt(distance2(P));
    }
};

Some tests:

TEST(RdpTest, Distance)
{
    auto seg = LineSegment({0.0, 0.0, 0.0}, {10.0, 0.0, 0.0});
    ASSERT_EQ(4.0, dbg(seg.distance({5.0, 4.0, 0.0})));
    ASSERT_EQ(5.0, dbg(seg.distance({-4.0, 3.0, 0.0})));
    ASSERT_EQ(5.0, dbg(seg.distance({14.0, 3.0, 0.0})));
}

TEST(RdpTest, Distance_dup)
{
    auto seg = LineSegment({0.0, 0.0, 0.0}, {0.0, 0.0, 0.0});
    ASSERT_EQ(5.0, dbg(seg.distance({3.0, 4.0, 0.0})));
    ASSERT_EQ(5.0, dbg(seg.distance({-4.0, 3.0, 0.0})));
    ASSERT_EQ(13.0, dbg(seg.distance({5.0, 12.0, 0.0})));
}
star-plu-cn-sk2 commented 2 years ago

这是来自QQ邮箱的假期自动回复邮件。你好,我最近正在休假中,无法亲自回复你的邮件。我将在假期结束后,尽快给你回复。

district10 commented 8 months ago

As https://github.com/fhirschmann/rdp/issues/13 points out, pdist in rdp is WRONGLY Point-to-Line distance. We (pybind11_rdp) use Point-to-LineSegment distance.

from rdp import rdp
print(rdp([[0, 0], [10, 0.1], [1, 0]], epsilon=1.0)) # wrong
# [[0.0, 0.0],
#  [1.0, 0.0]]

from pybind11_rdp import rdp
print(rdp([[0, 0], [10, 0.1], [1, 0]], epsilon=1.0)) # correct
# [[ 0.   0. ]
#  [10.   0.1]
#  [ 1.   0. ]]