kurbo/

mindist.rs

// Copyright 2021 the Kurbo Authors
// SPDX-License-Identifier: Apache-2.0 OR MIT

//! Minimum distance between two Bézier curves
//!
//! This implements the algorithm in "Computing the minimum distance between
//! two Bézier curves", Chen et al., *Journal of Computational and Applied
//! Mathematics* 229(2009), 294-301

use crate::Vec2;
use core::cmp::Ordering;

#[cfg(not(feature = "std"))]
use crate::common::FloatFuncs;

/// Given the control points of n- and m-dimensional Bézier curves in `bez1` and `bez2`
/// respectively, calculate the minimum distance between the two curves.
///
/// This should work for arbitrary Bézier curves up to and including degree 255, but only lines,
/// quadratics and cubics are tested.
///
/// Returns a tuple of the distance, the path time `t1` of the closest point on the first Bézier,
/// and the path time `t2` of the closest point on the second Bézier.
///
/// # Panics
///
/// This panics if more than 256 control points are given in `bez1` or `bez2`, i.e., if either
/// Bézier is of degree greater than 255.
pub(crate) fn min_dist_param(
    bez1: &[Vec2],
    bez2: &[Vec2],
    u: (f64, f64),
    v: (f64, f64),
    epsilon: f64,
    best_alpha: Option<f64>,
) -> (f64, f64, f64) {
    assert!(!bez1.is_empty() && !bez2.is_empty());
    let n = u8::try_from(bez1.len() - 1).expect("At most 256 control points");
    let m = u8::try_from(bez2.len() - 1).expect("At most 256 control points");
    let (umin, umax) = u;
    let (vmin, vmax) = v;
    let umid = (umin + umax) / 2.0;
    let vmid = (vmin + vmax) / 2.0;
    let svalues: [(f64, f64, f64); 4] = [
        (S(umin, vmin, bez1, bez2), umin, vmin),
        (S(umin, vmax, bez1, bez2), umin, vmax),
        (S(umax, vmin, bez1, bez2), umax, vmin),
        (S(umax, vmax, bez1, bez2), umax, vmax),
    ];
    let alpha: f64 = svalues.iter().map(|(a, _, _)| *a).reduce(f64::min).unwrap();
    if let Some(best) = best_alpha {
        if alpha > best {
            return (alpha, umid, vmid);
        }
    }

    if (umax - umin).abs() < epsilon || (vmax - vmin).abs() < epsilon {
        return (alpha, umid, vmid);
    }

    // Property one: D(r>k) > alpha
    let mut is_outside = true;
    let mut min_drk = None;
    let mut min_ij = None;
    for r in 0..(2 * n) {
        for k in 0..(2 * m) {
            let d_rk = D_rk(r, k, bez1, bez2);
            if d_rk < alpha {
                is_outside = false;
            }
            if min_drk.is_none() || Some(d_rk) < min_drk {
                min_drk = Some(d_rk);
                min_ij = Some((r, k));
            }
        }
    }
    if is_outside {
        return (alpha, umid, vmid);
    }

    // Property two: boundary check
    let mut at_boundary0on_bez1 = true;
    let mut at_boundary1on_bez1 = true;
    let mut at_boundary0on_bez2 = true;
    let mut at_boundary1on_bez2 = true;
    for i in 0..2 * n {
        for j in 0..2 * m {
            let dij = D_rk(i, j, bez1, bez2);
            let dkj = D_rk(0, j, bez1, bez2);
            if dij < dkj {
                at_boundary0on_bez1 = false;
            }
            let dkj = D_rk(2 * n, j, bez1, bez2);
            if dij < dkj {
                at_boundary1on_bez1 = false;
            }
            let dkj = D_rk(i, 0, bez1, bez2);
            if dij < dkj {
                at_boundary0on_bez2 = false;
            }
            let dkj = D_rk(i, 2 * m, bez1, bez2);
            if dij < dkj {
                at_boundary1on_bez2 = false;
            }
        }
    }
    if at_boundary0on_bez1 && at_boundary0on_bez2 {
        return svalues[0];
    }
    if at_boundary0on_bez1 && at_boundary1on_bez2 {
        return svalues[1];
    }
    if at_boundary1on_bez1 && at_boundary0on_bez2 {
        return svalues[2];
    }
    if at_boundary1on_bez1 && at_boundary1on_bez2 {
        return svalues[3];
    }

    let (min_i, min_j) = min_ij.unwrap();
    let new_umid = umin + (umax - umin) * (min_i as f64 / (2 * n) as f64);
    let new_vmid = vmin + (vmax - vmin) * (min_j as f64 / (2 * m) as f64);

    // Subdivide
    let results = [
        min_dist_param(
            bez1,
            bez2,
            (umin, new_umid),
            (vmin, new_vmid),
            epsilon,
            Some(alpha),
        ),
        min_dist_param(
            bez1,
            bez2,
            (umin, new_umid),
            (new_vmid, vmax),
            epsilon,
            Some(alpha),
        ),
        min_dist_param(
            bez1,
            bez2,
            (new_umid, umax),
            (vmin, new_vmid),
            epsilon,
            Some(alpha),
        ),
        min_dist_param(
            bez1,
            bez2,
            (new_umid, umax),
            (new_vmid, vmax),
            epsilon,
            Some(alpha),
        ),
    ];

    *results
        .iter()
        .min_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(Ordering::Less))
        .unwrap()
}

// $ S(u,v)=\sum_{r=0}^{2n} \sum _{k=0}^{2m} D_{r,k} B_{2n}^r(u) B_{2m}^k(v) $

#[allow(non_snake_case)]
fn S(u: f64, v: f64, bez1: &[Vec2], bez2: &[Vec2]) -> f64 {
    #[expect(
        clippy::cast_possible_truncation,
        reason = "By the checks in `min_dist_param` we know `bez1.len()` and `bez2.len()` are less than or equal to 256"
    )]
    let (n, m) = { ((bez1.len() - 1) as u8, (bez2.len() - 1) as u8) };
    let mut summand = 0.0;
    for r in 0..=2 * n {
        for k in 0..=2 * m {
            summand +=
                D_rk(r, k, bez1, bez2) * basis_function(2 * n, r, u) * basis_function(2 * m, k, v);
        }
    }
    summand
}

// $ C_{r,k} = ( \sum_{i=\theta}^\upsilon P_i C_n^i C_n^{r-i} / C_{2n}^r ) \cdot (\sum_{j=\sigma}^\varsigma Q_j C_m^j C_m^{k-j} / C_{2m}^k ) $
#[allow(non_snake_case)]
fn C_rk(r: u8, k: u8, bez1: &[Vec2], bez2: &[Vec2]) -> f64 {
    #[expect(
        clippy::cast_possible_truncation,
        reason = "By the checks in `min_dist_param` we know `bez1.len()` and `bez2.len()` are less than or equal to 256"
    )]
    let (n, m) = { ((bez1.len() - 1) as u8, (bez2.len() - 1) as u8) };
    let upsilon = r.min(n);
    let theta = r - n.min(r);
    let mut left: Vec2 = Vec2::ZERO;
    for i in theta..upsilon + 1 {
        let item = bez1[usize::from(i)];
        left += item * (choose(n, i) * choose(n, r - i)) as f64 / (choose(2 * n, r) as f64);
    }

    let varsigma = k.min(m);
    let sigma = k - m.min(k);
    let mut right: Vec2 = Vec2::ZERO;
    for j in sigma..varsigma + 1 {
        let item = bez2[usize::from(j)];
        right += item * (choose(m, j) * choose(m, k - j)) as f64 / (choose(2 * m, k) as f64);
    }
    left.dot(right)
}

// $ A_r=\sum _{i=\theta} ^\upsilon (P_i \cdot P_{r-i}) C_n^i C_n^{r-i} / C_{2n}^r $
// ($ B_k $ is just the same as $ A_r $ but for the other curve.)

#[allow(non_snake_case)]
fn A_r(r: u8, p: &[Vec2]) -> f64 {
    #[expect(
        clippy::cast_possible_truncation,
        reason = "By the checks in `min_dist_param` we know `p.len()` is less than or equal to 256"
    )]
    let n = (p.len() - 1) as u8;
    let upsilon = r.min(n);
    let theta = r - n.min(r);
    (theta..=upsilon)
        .map(|i| {
            let dot = p[usize::from(i)].dot(p[usize::from(r - i)]); // These are bounds checked by the sum limits
            let factor = (choose(n, i) * choose(n, r - i)) as f64 / (choose(2 * n, r) as f64);
            dot * factor
        })
        .sum()
}

#[allow(non_snake_case)]
fn D_rk(r: u8, k: u8, bez1: &[Vec2], bez2: &[Vec2]) -> f64 {
    // In the paper, B_k is used for the second factor, but it's the same thing
    A_r(r, bez1) + A_r(k, bez2) - 2.0 * C_rk(r, k, bez1, bez2)
}

// Bézier basis function
fn basis_function(n: u8, i: u8, u: f64) -> f64 {
    choose(n, i) as f64 * (1.0 - u).powi((n - i) as i32) * u.powi(i as i32)
}

// Binomial co-efficient, but returning zeros for values outside of domain
fn choose(n: u8, k: u8) -> u32 {
    let mut n = n as u32;
    let k = k as u32;
    if k > n {
        return 0;
    }
    let mut p = 1;
    for i in 1..=(n - k) {
        p *= n;
        p /= i;
        n -= 1;
    }
    p
}

#[cfg(test)]
mod tests {
    use crate::mindist::A_r;
    use crate::mindist::{D_rk, choose};
    use crate::param_curve::ParamCurve;
    use crate::{CubicBez, Line, PathSeg, Vec2};

    #[test]
    fn test_choose() {
        assert_eq!(choose(6, 0), 1);
        assert_eq!(choose(6, 1), 6);
        assert_eq!(choose(6, 2), 15);
    }

    #[test]
    fn test_d_rk() {
        let bez1 = vec![
            Vec2::new(129.0, 139.0),
            Vec2::new(190.0, 139.0),
            Vec2::new(201.0, 364.0),
            Vec2::new(90.0, 364.0),
        ];
        let bez2 = vec![
            Vec2::new(309.0, 159.0),
            Vec2::new(178.0, 159.0),
            Vec2::new(215.0, 408.0),
            Vec2::new(309.0, 408.0),
        ];
        let b = A_r(1, &bez2);
        assert!((b - 80283.0).abs() < 0.005, "B_1(Q)={b:?}");
        let d = D_rk(0, 1, &bez1, &bez2);
        assert!((d - 9220.0).abs() < 0.005, "D={d:?}");
    }

    #[test]
    fn test_mindist() {
        let bez1 = PathSeg::Cubic(CubicBez::new(
            (129.0, 139.0),
            (190.0, 139.0),
            (201.0, 364.0),
            (90.0, 364.0),
        ));
        let bez2 = PathSeg::Cubic(CubicBez::new(
            (309.0, 159.0),
            (178.0, 159.0),
            (215.0, 408.0),
            (309.0, 408.0),
        ));
        let mindist = bez1.min_dist(bez2, 0.001);
        assert!((mindist.distance - 50.9966).abs() < 0.5);
    }

    #[test]
    fn test_overflow() {
        let bez1 = PathSeg::Cubic(CubicBez::new(
            (232.0, 126.0),
            (134.0, 126.0),
            (139.0, 232.0),
            (141.0, 301.0),
        ));
        let bez2 = PathSeg::Line(Line::new((359.0, 416.0), (367.0, 755.0)));
        let mindist = bez1.min_dist(bez2, 0.001);
        assert!((mindist.distance - 246.4731222669117).abs() < 0.5);
    }

    #[test]
    fn test_out_of_order() {
        let bez1 = PathSeg::Cubic(CubicBez::new(
            (287.0, 182.0),
            (346.0, 277.0),
            (356.0, 299.0),
            (359.0, 416.0),
        ));
        let bez2 = PathSeg::Line(Line::new((141.0, 301.0), (152.0, 709.0)));
        let mindist1 = bez1.min_dist(bez2, 0.5);
        let mindist2 = bez2.min_dist(bez1, 0.5);
        assert!((mindist1.distance - mindist2.distance).abs() < 0.5);
    }

    #[test]
    fn test_line_curve() {
        let line = PathSeg::Line(Line::new((929.0, 335.0), (911.0, 340.0)));

        let line_as_bez = PathSeg::Cubic(CubicBez::new(
            line.eval(0.0),
            line.eval(1.0 / 3.0),
            line.eval(2.0 / 3.0),
            line.eval(1.0),
        ));

        let bez2 = PathSeg::Cubic(CubicBez::new(
            (1052.0, 401.0),
            (1048.0, 305.0),
            (1046.0, 216.0),
            (1054.0, 146.0),
        ));
        let mindist_as_bez = line_as_bez.min_dist(bez2, 0.5);
        let mindist_as_line = line.min_dist(bez2, 0.5);
        assert!((mindist_as_line.distance - mindist_as_bez.distance).abs() < 0.5);
    }
}