unoptimized for loop

[kkonevets]

This code works MUCH faster. It does not use for loop, instead it uses numpy vectorization

import numpy as np


def line_dists(points, start, end):
    if np.all(start == end):
        return np.linalg.norm(points - start, axis=1)

    vec = end - start
    cross = np.cross(vec, start - points)
    return np.divide(abs(cross), np.linalg.norm(vec))


def rdp(M, epsilon=0):
    M = np.array(M)
    start, end = M[0], M[-1]
    dists = line_dists(M, start, end)

    index = np.argmax(dists)
    dmax = dists[index]

    if dmax > epsilon:
        result1 = rdp(M[:index + 1], epsilon)
        result2 = rdp(M[index:], epsilon)

        result = np.vstack((result1[:-1], result2))
    else:
        result = np.array([start, end])

    return result

[JustasB]

I can confirm that it works faster than the pure python version for me as well.

[soichih]

@kkonevets I can't figure out how to convert your code to work on 3D data.. Is it possible to do something similar with 3D data?

[ivicajan]

Dear all,

is there a way to limit number of simplified points together with max distance (epsilon)?
Or to be more precise, to pick only subset of N the most relevant points for given dmax.

Thanks in advance!

Cheers,
Ivica

[darwinharianto]

import numpy as np

def line_dists2D(points, start, end):
    if np.all(start == end):
        return np.linalg.norm(points - start, axis=1)

    vec = end - start
    cross = np.cross(vec, start - points)
    return np.divide(abs(cross), np.linalg.norm(vec))


# https://stackoverflow.com/questions/56463412/distance-from-a-point-to-a-line-segment-in-3d-python
def lineseg_dist3D(p: np.ndarray, a: np.ndarray, b: np.ndarray):

    # normalized tangent vector
    d = np.divide(b - a, np.linalg.norm(b - a))

    # signed parallel distance components
    s = np.dot(a - p, d)
    t = np.dot(p - b, d)

    # clamped parallel distance
    h = np.maximum.reduce([s, t, np.zeros(len(p))])

    # perpendicular distance component
    c = np.cross(p - a, d)

    res = np.hypot(h, np.linalg.norm(c, axis=1))
    res[res < np.finfo(np.float64).eps] = 0

    return res


def rdp(M, epsilon=0):
    M = np.array(M)
    start, end = M[0], M[-1]
    if M.shape[1] == 2:
        dists = line_dists2D(M, start, end)
    elif M.shape[1] == 3:
        dists = lineseg_dist3D(M, start, end)
    else:
        raise ValueError("point dimension must be 2d or 3d")
    index = np.argmax(dists)
    dmax = dists[index]

    if dmax == np.nan:
        raise ValueError

    if dmax > epsilon:
        result1 = rdp(M[: index + 1], epsilon)
        result2 = rdp(M[index:], epsilon)

        result = np.vstack((result1[:-1], result2))
    else:
        result = np.array([start, end])

    return result

  points = np.array([1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4]).reshape(-1, 3)
  points = np.array([1, 1, 2, 2, 3, 3, 4, 4, 3, 3, 2, 3]).reshape(-1, 2)
  print(rdp(points))

@soichih
I tried to make a sample on 3d Points, i think this is working

edit: fixed close to 0 and linalg norm over rows

[star-plu-cn-sk2]

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