douglas_peucker.cpp

#include "internal/douglas_peucker.h"
 
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <vector>
 
// --- Signed perpendicular distance of P from the directed line A->B ---
// Positive = P lies to the "left" of the A->B travel direction (using the
// n = (-dy, dx) normal). Falls back to the point distance for a zero-length AB.
static inline float signed_left_dist(const Point& A, const Point& B, const Point& P) {
    const float dx = B.x - A.x, dy = B.y - A.y;
    const float l = std::sqrt(dx * dx + dy * dy);
    if (l < 1e-9f)
        return std::sqrt(Point::distSq(P, A));
    return ((P.y - A.y) * dx - (P.x - A.x) * dy) / l;
}
 
static void dp_reduce(
    const std::vector<Point>& pts, int lo, int hi, float eps, float retract_eps,
    float interior_sign, std::vector<uint8_t>& keep
) {
    if (hi - lo < 2)
        return; // no interior vertices to drop
 
    const Point& A = pts[lo];
    const Point& B = pts[hi];
 
    float max_abs_dev = 0.0f; // worst |deviation| -> fidelity / split choice
    float max_retract = 0.0f; // worst inward move -> watertightness guard
    int split = lo;
    for (int k = lo + 1; k < hi; ++k) {
        const float d = signed_left_dist(A, B, pts[k]);
        const float adev = std::fabs(d);
        if (adev > max_abs_dev) {
            max_abs_dev = adev;
            split = k;
        }
        // interior_sign * d > 0 => vertex is interior (concave, safe to flatten);
        // < 0 => vertex is exterior (convex), dropping it retracts the boundary.
        const float retract = -interior_sign * d;
        if (retract > max_retract)
            max_retract = retract;
    }
 
    if (max_abs_dev <= eps && max_retract <= retract_eps)
        return; // drop the run
 
    keep[split] = 1;
    dp_reduce(pts, lo, split, eps, retract_eps, interior_sign, keep);
    dp_reduce(pts, split, hi, eps, retract_eps, interior_sign, keep);
}
 
void dp_curve_reduction(
    const std::vector<std::vector<Point>>& chains, const std::vector<std::vector<uint8_t>>& fixed,
    std::vector<std::vector<QuadBezier>>& results, float eps
) {
    float retract_eps = std::min(eps, 0.5f);
 
    for (size_t i = 0; i < chains.size(); ++i) {
        const std::vector<Point>& chain = chains[i];
        const int n = static_cast<int>(chain.size());
        std::vector<QuadBezier> result;
        if (n < 2) {
            results.push_back(result);
            continue;
        }
 
        // Orientation: sign of the shoelace area decides which side of a chord is
        // the region interior (see dp_reduce). Computed over the closed chain.
        double sa = 0.0;
        for (int k = 0; k < n; ++k) {
            const Point& p = chain[k];
            const Point& q = chain[(k + 1) % n];
            sa += static_cast<double>(p.x) * q.y - static_cast<double>(q.x) * p.y;
        }
        const float interior_sign = (sa > 0.0) ? 1.0f : -1.0f;
 
        // Segment boundaries: chain ends plus every interior fixed point. DP runs
        // inside each [bounds[s], bounds[s+1]] span with its ends pinned.
        std::vector<int> bounds;
        bounds.push_back(0);
        for (int k = 1; k < n - 1; ++k)
            if (k < static_cast<int>(fixed[i].size()) && fixed[i][k])
                bounds.push_back(k);
        bounds.push_back(n - 1);
 
        std::vector<uint8_t> keep(n, 0);
        for (int b : bounds)
            keep[b] = 1;
        for (size_t s = 0; s + 1 < bounds.size(); ++s)
            dp_reduce(chain, bounds[s], bounds[s + 1], eps, retract_eps, interior_sign, keep);
 
        // Emit kept points in order as straight-line quads.
        int prev = -1;
        for (int k = 0; k < n; ++k) {
            if (!keep[k])
                continue;
            if (prev >= 0) {
                const Point& P = chain[prev];
                const Point& Q = chain[k];
                result.push_back({P, (P + Q) * 0.5f, Q});
            }
            prev = k;
        }
        results.push_back(result);
    }
}