#!/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: float; y: float;

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

    # Distance carrée à un segment [a, b]
    dist2_to_segment(self, a: Point, b: Point): float => {
        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
        }
    }
}

furthest_index = (points, start, end, eps2): int => {
    best_i = -1
    best_d2 = -1
    a = points[start]
    b = points[end]
    i = start + 1
    while i < end {
        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 { best_i } else { -1 }
}

# Ramer-Douglas-Peucker (retourne une nouvelle polyligne)
rdp = (points: list[Point], epsilon: float): list[Point] => {
    n = len(points)
    if n <= 2 {
        points
    } else {
        eps2 = epsilon * epsilon

        keep = list()
        i = 0
        while i < n {
            keep.append(True)
            i = i + 1
        }

        start_stack = list(0)
        end_stack = list(n - 1)
        while len(start_stack) > 0 {
            start = start_stack[len(start_stack) - 1]
            start_stack = start_stack[:len(start_stack) - 1]
            end = end_stack[len(end_stack) - 1]
            end_stack = end_stack[:len(end_stack) - 1]

            if end <= start + 1 { continue }

            best_i = furthest_index(points, start, end, eps2)

            if best_i >= 0 {
                start_stack = start_stack + list(start, best_i)
                end_stack = end_stack + list(best_i, end)
            } else {
                j = start + 1
                while j < end {
                    keep[j] = False
                    j = j + 1
                }
            }
        }

        result = list()
        i = 0
        while i < n {
            if keep[i] { result.append(points[i]) }
            i = i + 1
        }
        result
    }
}

# 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)}")
pretty = " → ".join(simplified.[(p) => { p.display() }])
print(f"Résultat : {pretty}")