Feature proposal: improve `approximate_polygon_dp` on closed contours.

[SamDM]

I am working on rust port of the OpenCV2 aruco detector code, for which I'm leaning on some imageproc functions.
In doing so, I encountered a problem with the approximate_polygon_dp algorithm on closed contours, which caused my detector to fail on fairly trivial cases.

Below I illustrate the problem, give a minimal reproducible example, and provide a robuster alternative, based on the OpenCV implementation of the DP algorithm. If this is deemed to be valuable, I can make a PR to apply this in the imageproc code.

The problem

The current implementation sometimes performs fairly poorly in approximating corners of closed contours. Even though the result is not in violation of the "epsilon" bound, it is suboptimal for my application, and perhaps for others too.
My application is finding aruco markers in images (the QR-code looking things below) with as precise corner localization as possible. This has many uses in robotics.

Before (current implementation):

The blue rectangles are approximated contours found by the approximate_polygon_dp function in imageproc.
Notice the poor alignment to the actual corners of each tag. The two left tags of the first row were discarded because the current implementation returned a 5-point polygon instead of a 4-point polygon.

After (opencv based implementation):

While this aruco application is not a direct concern here, I want to illustrate that the current implementation has limitations that affect real-world usage.

MRE

The script below provides an MRE, the first ~260 LOC are the code for the OpenCV based implementation of the DP polygon approximation algorithm. The rest of the code sets up some pathological examples and runs the current and proposed implementation on them.

#!/usr/bin/env -S cargo +nightly -Zscript
---
[package]
edition = "2024"

[dependencies]
image = "0.25"
imageproc = "0.26"
clap = { version = "4", features = ["derive"] }
num-traits = "0.2"
---

use clap::Parser;
use image::{imageops, DynamicImage, GrayImage, Luma, Rgb, RgbImage};
use num_traits::NumCast;
use imageproc::contours::find_contours;
use imageproc::drawing::{draw_line_segment_mut, draw_polygon_mut};
use imageproc::geometry::approximate_polygon_dp;
use imageproc::point::Point;

// ── Improved Ramer-Douglas-Peucker ───────────────────────────────────────────
//
// Two changes relative to `imageproc::approximate_polygon_dp`:
//
//   1. Distance to the chord **segment** (clamped to endpoints) instead of
//      the infinite line through the chord endpoints.
//
//   2. For closed contours: locate two approximate-diameter endpoints A and B
//      first, split the ring into arcs A→B and B→A, and run the standard
//      open-arc DP on each independently.  The existing approach runs a single
//      open-chain DP whose base chord connects the two *adjacent* endpoints of
//      the Suzuki-Abe contour, which can produce incorrect corners.
//
// A single-pass cleanup step removes output points that lie nearly on the line
// between their neighbours.

/// Cast a `Point<T>` to `Point<f64>`.
/// `imageproc` has `pub(crate)` method `Point::to_f64(&self)`, which I can't access here.
/// Hence this quick re-implementation
fn to_f64<T: NumCast + Copy>(p: Point<T>) -> Point<f64> {
    Point::new(
        NumCast::from(p.x).unwrap(),
        NumCast::from(p.y).unwrap(),
    )
}

/// Distance from `pt` to the segment `ab`, clamped to the nearest endpoint
/// when `pt` does not project onto the segment interior.
///
/// Drop-in replacement for `Line::distance_from_point` — the only change is
/// "infinite line" → "segment with endpoint clamping".
fn segment_distance<T: NumCast + Copy>(pt: Point<T>, a: Point<T>, b: Point<T>) -> f64 {
    let pt = to_f64(pt);
    let a = to_f64(a);
    let b = to_f64(b);
    let dx = b.x - a.x;
    let dy = b.y - a.y;
    let seg_len2 = dx * dx + dy * dy;
    if seg_len2 == 0.0 {
        let ex = pt.x - a.x;
        let ey = pt.y - a.y;
        return (ex * ex + ey * ey).sqrt();
    }
    let ptx = pt.x - a.x;
    let pty = pt.y - a.y;
    let t = ((ptx * dx + pty * dy) / seg_len2).clamp(0.0, 1.0);
    let ex = ptx - t * dx;
    let ey = pty - t * dy;
    (ex * ex + ey * ey).sqrt()
}

/// Open-arc Ramer-Douglas-Peucker.
///
/// Identical in structure to `imageproc::approximate_polygon_dp` but calls
/// `segment_distance` instead of `Line::distance_from_point`.  Does not
/// handle the `closed` case — that logic lives in `approximate_polygon_dp_new`.
fn dp_arc<T: NumCast + Copy>(curve: &[Point<T>], epsilon: f64) -> Vec<Point<T>> {
    if curve.len() < 2 {
        return curve.to_vec();
    }

    // Find the point with the maximum distance from the base chord.
    let mut dmax = 0.0;
    let mut index = 0;
    let end = curve.len() - 1;
    for (i, &point) in curve.iter().enumerate().skip(1) {
        let d = segment_distance(point, curve[0], curve[end]);
        if d > dmax {
            index = i;
            dmax = d;
        }
    }

    // If max distance is greater than epsilon, recursively simplify.
    if dmax > epsilon {
        let mut partial1 = dp_arc(&curve[0..=index], epsilon);
        let partial2 = dp_arc(&curve[index..=end], epsilon);
        partial1.pop();
        partial1.extend(partial2);
        partial1
    } else {
        vec![curve[0], curve[end]]
    }
}

/// Remove output points that lie nearly on the line between their neighbours.
///
/// The DP guarantees that every *removed* point was within `ε` of its arc
/// chord.  It does not guarantee that every *kept* point is far from the line
/// through its output neighbours.  Near the diameter endpoints A and B the
/// two arcs are stitched together, so a point that was legitimately far from
/// its own arc chord can still be nearly collinear with neighbours that came
/// from the other arc.  This pass catches those cases.
///
/// A point `pt` (predecessor `start`, successor `end`) is dropped when its
/// perpendicular distance to the line `start→end` is below `ε` and `pt`
/// projects onto the segment (not beyond either endpoint).
///
/// Single-pass; when a point is dropped the window immediately advances to
/// the next `end_pt`, enabling cascaded removal of collinear runs.
fn cleanup<T: NumCast + Copy>(points: Vec<Point<T>>, epsilon: f64, closed: bool) -> Vec<Point<T>> {
    let n = points.len();
    if n <= 2 {
        return points;
    }

    let src = points;
    let mut keep = vec![true; n];
    let mut new_count = n;

    // Closed: window wraps — start_pt begins as the last point, loop covers all n.
    // Open:   endpoints are never removed — start_pt begins as src[0], loop covers [1..n-2].
    let (mut start_pt, mut pt, mut pos, mut wpos, i_start, i_end) = if closed {
        (to_f64(src[n - 1]), to_f64(src[0]), 1usize, 0usize, 0usize, n)
    } else {
        (to_f64(src[0]), to_f64(src[1]), 2usize, 1usize, 1usize, n - 1)
    };

    let mut i = i_start;
    while i < i_end && new_count > 2 {
        let end_pt = to_f64(src[pos % n]);
        pos += 1;

        let dx = end_pt.x - start_pt.x;
        let dy = end_pt.y - start_pt.y;
        let seg_len = (dx * dx + dy * dy).sqrt();

        // Perpendicular distance from `pt` to the infinite line through start→end.
        let perp_dist = ((pt.y - start_pt.y) * dx - (pt.x - start_pt.x) * dy).abs() / seg_len;
        // Positive when `pt` projects onto the segment (not past either endpoint).
        let projects_between = (pt.x - start_pt.x) * (end_pt.x - pt.x)
            + (pt.y - start_pt.y) * (end_pt.y - pt.y) >= 0.0;

        if perp_dist <= epsilon && projects_between {
            keep[wpos] = false;
            new_count -= 1;
            start_pt = end_pt;
            pt = to_f64(src[pos % n]);
            pos += 1;
            i += 2;
            continue;
        }

        start_pt = pt;
        pt = end_pt;
        wpos = (wpos + 1) % n;
        i += 1;
    }

    src.into_iter()
        .zip(keep)
        .filter_map(|(p, k)| k.then_some(p))
        .collect()
}

/// Ramer-Douglas-Peucker polygon approximation with correct closed-contour
/// handling.
///
/// This implementation differs from `imageproc::approximate_polygon_dp` in
/// two ways — see the module-level comment above for details.
pub fn approximate_polygon_dp_new<T: NumCast + Copy>(
    contour: &[Point<T>],
    epsilon: f64,
    closed: bool,
) -> Vec<Point<T>> {
    let n = contour.len();
    if n == 0 {
        return vec![];
    }
    if n == 1 {
        return vec![contour[0]];
    }

    let eps2 = epsilon * epsilon;

    if !closed {
        return dp_arc(contour, epsilon);
    }

    // ── Closed contour: locate two approximate-diameter endpoints A and B ─────
    //
    // Three "farthest-from-current" iterations starting at index 0.  Each
    // iteration moves the start to the previously found farthest point,
    // converging on a stable diameter pair.

    let mut a_idx = 0usize;
    let mut b_offset = 0usize;

    for _ in 0..3 {
        a_idx = (a_idx + b_offset) % n;
        let start_pt = to_f64(contour[a_idx]);

        let mut max_d2 = 0.0f64;
        b_offset = 0;
        for j in 1..n {
            let pt = to_f64(contour[(a_idx + j) % n]);
            let dx = pt.x - start_pt.x;
            let dy = pt.y - start_pt.y;
            let d2 = dx * dx + dy * dy;
            if d2 > max_d2 {
                max_d2 = d2;
                b_offset = j;
            }
        }
        // All points within an epsilon-ball — collapse to a single point.
        if max_d2 <= eps2 {
            return vec![contour[a_idx]];
        }
    }

    let b_idx = (a_idx + b_offset) % n;

    // Materialise both arcs as contiguous slices for the recursive DP.
    let collect_arc = |start: usize, end: usize| -> Vec<Point<T>> {
        let mut v = Vec::new();
        let mut i = start;
        loop {
            v.push(contour[i]);
            if i == end {
                break;
            }
            i = (i + 1) % n;
        }
        v
    };

    // Run open-arc DP on arc A→B and arc B→A independently.
    let mut left = dp_arc(&collect_arc(a_idx, b_idx), epsilon);
    let mut right = dp_arc(&collect_arc(b_idx, a_idx), epsilon);

    // Merge: `left` ends at B (= `right[0]`), `right` ends at A (= `left[0]`).
    left.pop(); // drop duplicate B
    right.pop(); // drop duplicate A
    left.extend(right);

    // Remove any nearly-collinear output points (single-pass).
    cleanup(left, epsilon, true)
}


/// MRE for imageproc::approximate_polygon_dp closed-contour inaccuracy.
///
/// For each SIZE,ANGLE case, saves a 4-column grid PNG:
///   col 1: drawn rectangle
///   col 2: contour points (blue, green=start, red=end)
///   col 3: imageproc (current) approx polygon
///   col 4: opencv based approx polygon
#[derive(Parser)]
struct Args {
    /// SIZE,ANGLE pairs. Defaults show the failure modes:
    ///   100,5  → 4 corners but BR is ~5 px off
    ///   50,-7  → 5 corners instead of 4
    ///   60,-5  → 6 corners instead of 4
    /// And one correct result:
    ///   60,0   → 4 corners, at most ~1 px off
    #[arg(default_values = ["100,5", "50,-7", "60,-5", "60,0"])]
    cases: Vec<String>,
}

fn to_rgb(gray: &GrayImage) -> RgbImage {
    DynamicImage::from(gray.clone()).to_rgb8()
}

/// Greedy nearest-neighbour matching of `discovered` points to `gt` points.
/// Returns a `Vec` parallel to `discovered`: `Some(gt_index)` or `None`.
/// Each GT point is used at most once.
fn greedy_match(discovered: &[Point<i32>], gt: &[Point<i32>]) -> Vec<Option<usize>> {
    let mut matched = vec![None; discovered.len()];
    let mut used = vec![false; gt.len()];
    let mut candidates: Vec<(f64, usize, usize)> = (0..discovered.len())
        .flat_map(|i| {
            (0..gt.len()).map(move |j| {
                let p = discovered[i];
                let q = gt[j];
                let d = (((p.x - q.x).pow(2) + (p.y - q.y).pow(2)) as f64).sqrt();
                (d, i, j)
            })
        })
        .collect();
    candidates.sort_by(|a, b| a.0.total_cmp(&b.0));
    for (_, i, j) in candidates {
        if matched[i].is_none() && !used[j] {
            matched[i] = Some(j);
            used[j] = true;
        }
    }
    matched
}

fn run(size: u32, angle_deg: f64) {
    let tag = format!("{size}px_{angle_deg:+.0}deg");

    // ── 1. Draw rotated square ────────────────────────────────────────────────
    let sz = size as f64;
    let (cx, cy) = (sz / 2.0, sz / 2.0);
    let theta = angle_deg.to_radians();
    let (ct, st) = (theta.cos(), theta.sin());
    // Largest half-extent that fits the canvas, with 10 % margin.
    let half = (sz / 2.0) / (ct.abs() + st.abs()) * 0.90;
    let gt_poly: Vec<Point<i32>> = [(-half, -half), (half, -half), (half, half), (-half, half)]
        .iter()
        .map(|&(dx, dy)| {
            Point::new(
                (cx + dx * ct - dy * st).round() as i32,
                (cy + dx * st + dy * ct).round() as i32,
            )
        })
        .collect();

    let mut gray = GrayImage::from_pixel(size, size, Luma([255u8]));
    draw_polygon_mut(&mut gray, &gt_poly, Luma([0u8]));
    let panel1 = to_rgb(&gray);

    // ── 2. Contour ────────────────────────────────────────────────────────────
    let contour = find_contours::<i32>(&gray)
        .into_iter()
        .max_by_key(|c| c.points.len())
        .expect("no contours found");
    let cont_pts = &contour.points;
    let n = cont_pts.len(); // approximation of "actual" length, can be up to ~1/√2 off

    let mut panel2 = to_rgb(&gray);
    for p in cont_pts {
        panel2.put_pixel(p.x as u32, p.y as u32, Rgb([0, 0, 255]));
    }
    panel2.put_pixel(cont_pts[0].x as u32, cont_pts[0].y as u32, Rgb([0, 255, 0]));
    panel2.put_pixel(cont_pts[n - 1].x as u32, cont_pts[n - 1].y as u32, Rgb([255, 0, 0]));

    // ── 3. Approximate polygon — two implementations ──────────────────────────
    let epsilon = n as f64 * 0.03;

    // 3a. Current imageproc implementation.
    let poly_current = approximate_polygon_dp(cont_pts, epsilon, true);
    // 3b. OpenCV-based implementation.
    let poly_opencv = approximate_polygon_dp_new(cont_pts, epsilon, true);

    let make_panel = |pts: &Vec<Point<i32>>| {
        let m = pts.len();
        let mut panel = to_rgb(&gray);
        for i in 0..m {
            let a = pts[i];
            let b = pts[(i + 1) % m];
            draw_line_segment_mut(
                &mut panel,
                (a.x as f32, a.y as f32),
                (b.x as f32, b.y as f32),
                Rgb([0u8, 255, 255]),
            );
        }
        for p in pts {
            panel.put_pixel(p.x as u32, p.y as u32, Rgb([0, 0, 255]));
        }
        if !pts.is_empty() {
            panel.put_pixel(pts[0].x as u32, pts[0].y as u32, Rgb([255, 155, 0]));
        }
        panel
    };

    let panel3 = make_panel(&poly_current);
    let panel4 = make_panel(&poly_opencv);

    // ── Assemble 1×4 grid ─────────────────────────────────────────────────────
    let mut grid = RgbImage::new(size * 4, size);
    imageops::overlay(&mut grid, &panel1, 0, 0);
    imageops::overlay(&mut grid, &panel2, size as i64, 0);
    imageops::overlay(&mut grid, &panel3, size as i64 * 2, 0);
    imageops::overlay(&mut grid, &panel4, size as i64 * 3, 0);
    grid.save(format!("{tag}.png")).unwrap();

    // ── Summary ───────────────────────────────────────────────────────────────
    let corner_labels = ["TL", "TR", "BR", "BL"];

    for (label, poly_pts) in [("imageproc (current):", &poly_current), ("opencv based:", &poly_opencv)] {
        let m = poly_pts.len();
        let matched = greedy_match(poly_pts, &gt_poly);
        println!("  [{label}] {m} corners:");
        for (i, p) in poly_pts.iter().enumerate() {
            let is_anchor = if i == 0 { "anchor" } else { "      " };
            let disc = format!("({:4},{:4})", p.x, p.y);
            let (gt_str, delta_str) = match matched[i] {
                Some(j) => {
                    let q = gt_poly[j];
                    let lbl = corner_labels.get(j).copied().unwrap_or("?");
                    (
                        format!("{} ({:4},{:4})", lbl, q.x, q.y),
                        format!("({:+4},{:+4})", p.x - q.x, p.y - q.y),
                    )
                }
                None => (format!("{:<14}", "—"), format!("{:<11}", "—")),
            };
            println!("    {is_anchor}  {disc}  →  {gt_str}  Δ{delta_str}");
        }
    }
    println!("'{tag}.png' saved  (contour={n}pts, ε={epsilon:.2})\n");
}

fn parse_case(s: &str) -> (u32, f64) {
    let (size_s, angle_s) = s
        .split_once(',')
        .unwrap_or_else(|| panic!("case must be SIZE,ANGLE — got {s:?}"));
    (
        size_s
            .trim()
            .parse()
            .unwrap_or_else(|_| panic!("invalid size in {s:?}")),
        angle_s
            .trim()
            .parse()
            .unwrap_or_else(|_| panic!("invalid angle in {s:?}")),
    )
}

fn main() {
    for case in Args::parse().cases {
        let (size, angle) = parse_case(&case);
        run(size, angle);
    }
}

This creates the following illustrations.
The first column is the input data for contour finding (not directly relevant here), the second column is the input contour to the DP algorithm. The third is the current output, and the fourth is the proposed alternative implementation.
The green and red points are start/end of the input contour.
The orange point is where the contour is "glued" in the output.
The light blue lines are the simplified segments.

The stdout of the script:

$ ./mre.rs 
   Compiling mre v0.0.0 (/workspaces/arucotools-rs/_scratch/approximate_polygon_dp_mre/mre.rs)
    Finished `dev` profile [unoptimized + debuginfo] target(s) in 0.85s
     Running `/home/vscode/.cargo/build/be/dbeff7b667e356/target/debug/mre`
  [imageproc (current):] 4 corners:
    anchor  (  11,   5)  →  TL (  12,   5)  Δ(  -1,  +0)
            (  96,  12)  →  TR (  95,  12)  Δ(  +1,  +0)
            (  83,  96)  →  BR (  88,  95)  Δ(  -5,  +1)
            (   4,  88)  →  BL (   5,  88)  Δ(  -1,  +0)
  [opencv based:] 4 corners:
    anchor  (  12,   4)  →  TL (  12,   5)  Δ(  +0,  -1)
            (  96,  12)  →  TR (  95,  12)  Δ(  +1,  +0)
            (  88,  96)  →  BR (  88,  95)  Δ(  +0,  +1)
            (   4,  88)  →  BL (   5,  88)  Δ(  -1,  +0)
'100px_+5deg.png' saved  (contour=336pts, ε=10.08)

  [imageproc (current):] 5 corners:
    anchor  (  38,   3)  →  —               Δ—          
            (  44,   3)  →  TR (  43,   3)  Δ(  +1,  +0)
            (  49,  43)  →  BR (  48,  43)  Δ(  +1,  +0)
            (  11,  49)  →  BL (   7,  48)  Δ(  +4,  +1)
            (   2,   7)  →  TL (   3,   7)  Δ(  -1,  +0)
  [opencv based:] 4 corners:
    anchor  (  43,   2)  →  TR (  43,   3)  Δ(  +0,  -1)
            (  49,  43)  →  BR (  48,  43)  Δ(  +1,  +0)
            (   7,  49)  →  BL (   7,  48)  Δ(  +0,  +1)
            (   2,   7)  →  TL (   3,   7)  Δ(  -1,  +0)
'50px_-7deg.png' saved  (contour=166pts, ε=4.98)

  [imageproc (current):] 6 corners:
    anchor  (  46,   3)  →  —               Δ—          
            (  54,   3)  →  TR (  53,   3)  Δ(  +1,  +0)
            (  58,  53)  →  BR (  57,  53)  Δ(  +1,  +0)
            (  13,  58)  →  —               Δ—          
            (   6,  57)  →  BL (   7,  57)  Δ(  -1,  +0)
            (   2,   7)  →  TL (   3,   7)  Δ(  -1,  +0)
  [opencv based:] 4 corners:
    anchor  (  53,   2)  →  TR (  53,   3)  Δ(  +0,  -1)
            (  58,  53)  →  BR (  57,  53)  Δ(  +1,  +0)
            (   7,  58)  →  BL (   7,  57)  Δ(  +0,  +1)
            (   2,   7)  →  TL (   3,   7)  Δ(  -1,  +0)
'60px_-5deg.png' saved  (contour=204pts, ε=6.12)

  [imageproc (current):] 4 corners:
    anchor  (   2,   3)  →  TL (   3,   3)  Δ(  -1,  +0)
            (  58,   3)  →  TR (  57,   3)  Δ(  +1,  +0)
            (  57,  58)  →  BR (  57,  57)  Δ(  +0,  +1)
            (   3,  58)  →  BL (   3,  57)  Δ(  +0,  +1)
  [opencv based:] 4 corners:
    anchor  (   2,   3)  →  TL (   3,   3)  Δ(  -1,  +0)
            (  57,   2)  →  TR (  57,   3)  Δ(  +0,  -1)
            (  58,  57)  →  BR (  57,  57)  Δ(  +1,  +0)
            (   3,  58)  →  BL (   3,  57)  Δ(  +0,  +1)
'60px_+0deg.png' saved  (contour=220pts, ε=6.60)

The proposed alternative

The reason this happens is because the current implementation does not allow changing the start and end point of a closed contour input.
Therefore, the furthest point from the line given by the start/end point is picked as the first "split" point. This point is often sub-optimal, i.e. better alternative points exist that result in the same (or less) amount of points such that the output is still within "epsilon" and/or where the resulting points are all closer to the original contour.

• The opencv implementation, for a closed loop, jumps around the closed contour for 3 iterations to find two points with maximal distance, then splits on that.

• It also computes distances to a line segment, which is not the same as the distance to an infinite line.

• In the end, it also does an extra pass to remove near collinear points that can exist after gluing segments together.

The opencv port I put in the MRE intetionally matches the existing imageproc code. But I have an alternative implementation too that keeps a bunch of the speedup tricks employed by the opencv devs.
I didn't put that optimized version in this MRE because the core logic is near incomprehensible due to all the optimizations applied.
But I can put the more "optimized" version in the PR, after some benchmarking to make sure the illegibility is worth the speed tradeoff.

Final disclaimer, I used Claude Code heavily for generating the alternative implementation, and most of this MRE. But I promise I am a human :)