#!/usr/bin/env catnip
# RDP simplification
#
# Ramer-Douglas-Peucker : simplification de polyligne en 2D.
# Élimine les points dont la distance au segment est inférieure à epsilon.
#
# Ref : https://en.wikipedia.org/wiki/Ramer%E2%80%93Douglas%E2%80%93Peucker_algorithm

struct Point {
    x; y;

    display(self) => { f"({self.x}, {self.y})" }

    # Distance carrée à un segment [a, b]
    dist2_to_segment(self, a, b) => {
        abx = b.x - a.x
        aby = b.y - a.y
        apx = self.x - a.x
        apy = self.y - a.y

        ab2 = abx * abx + aby * aby
        if ab2 == 0 {
            # Segment dégénéré (a == b)
            dx = self.x - a.x
            dy = self.y - a.y
            dx * dx + dy * dy
        } else {
            t = (apx * abx + apy * aby) / ab2
            if t < 0 { t = 0 }
            if t > 1 { t = 1 }
            qx = a.x + t * abx
            qy = a.y + t * aby
            dx = self.x - qx
            dy = self.y - 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 = points[i].dist2_to_segment(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(
    Point(0.0, 0.0),
    Point(1.0, 0.2),
    Point(2.0, -0.1),
    Point(3.0, 5.0),
    Point(4.0, 6.0),
    Point(5.0, 7.0),
    Point(6.0, 8.1),
    Point(7.0, 9.0),
    Point(8.0, 9.2),
    Point(9.0, 9.0),
)

epsilon = 0.75
simplified = rdp(polyline, epsilon)

print(f"Points initiaux : {len(polyline)}")
print(f"Points simplifiés : {len(simplified)}")
print(f"Résultat : {simplified}")