codex/geometry/rdp_simplification.cat
# RDP simplification

# Distance carrée point-segment (2D)
dist2_point_segment = (p, a, b) => {
    ax = a[0]
    ay = a[1]
    bx = b[0]
    by = b[1]
    px = p[0]
    py = p[1]

    abx = bx - ax
    aby = by - ay
    apx = px - ax
    apy = py - ay

    ab2 = abx * abx + aby * aby
    if ab2 == 0 {
        dx = px - ax
        dy = py - ay
        dx * dx + dy * dy
    } else {
        t = (apx * abx + apy * aby) / ab2
        if t < 0 {
            t = 0
        }
        if t > 1 {
            t = 1
        }
        qx = ax + t * abx
        qy = ay + t * aby
        dx = px - qx
        dy = py - qy
        dx * dx + dy * dy
    }
}

# Ramer-Douglas-Peucker (retourne une nouvelle polyligne)
rdp = (points, epsilon) => {
    n = len(points)
    if n <= 2 {
        points
    } else {
        eps2 = epsilon * epsilon
        a = points[0]
        b = points[n - 1]

        best_i = -1
        best_d2 = -1

        i = 1
        while i < n - 1 {
            d2 = dist2_point_segment(points[i], a, b)
            if d2 > best_d2 {
                best_d2 = d2
                best_i = i
            }
            i = i + 1
        }

        if best_d2 > eps2 {
            left = rdp(points[:best_i + 1], epsilon)
            right = rdp(points[best_i:], epsilon)
            left[:-1] + right
        } else {
            list(a, b)
        }
    }
}

# Démo
polyline = list(
    list(0.0, 0.0),
    list(1.0, 0.2),
    list(2.0, -0.1),
    list(3.0, 5.0),
    list(4.0, 6.0),
    list(5.0, 7.0),
    list(6.0, 8.1),
    list(7.0, 9.0),
    list(8.0, 9.2),
    list(9.0, 9.0)
)

epsilon = 0.75
simplified = rdp(polyline, epsilon)

print("Points initiaux :", len(polyline))
print("Points simplifiés :", len(simplified))
print("Résultat :", simplified)