main
Caleb Publish 0.2.0 cbe04a0 11h ago
//! Spatial index: one bucket grid per copper layer, plus one for holes and one for the edge.
//! Items can be added and removed as routes are ripped up.

use crate::board::NetId;
use crate::geom::*;
use std::cell::{Cell, RefCell};
use std::rc::Rc;

#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum Owner {
    Pad(u32),
    /// A committed route (connection id) on this run.
    Conn(u32),
    /// Copper that was on the board before this run; it never moves.
    Fixed,
    Zone(u32),
    Edge,
    Keepout,
    /// A fine-pitch pin's way out, held until its route is laid (4.7). Other nets keep clear.
    Lane(u32),
}

#[derive(Clone, Debug)]
pub struct Obs {
    pub shape: Shape,
    pub net: NetId,
    pub owner: Owner,
    pub bb: (f64, f64, f64, f64),
    /// A big pour's edges by place, so a query near it doesn't test them all.
    pub fast: Option<Rc<PolyIndex>>,
}

impl Obs {
    /// Edge distance from a segment to this item, exact up to `limit`; past that it may answer
    /// anything at or above `limit`. Pours of 700 to 1100 points took a search on benchmark B to
    /// 25 microseconds a node, against about 1 elsewhere.
    #[inline]
    pub fn dist_within(&self, a: V, b: V, limit: f64) -> f64 {
        match &self.fast {
            Some(f) => f.dist_within(a, b, limit),
            None => self.shape.dist_seg(a, b),
        }
    }
}

/// A polygon's edges in buckets: by cell for distance, by row for inside or out.
#[derive(Debug)]
pub struct PolyIndex {
    pts: Vec<V>,
    r: f64,
    x0: f64,
    y0: f64,
    cell: f64,
    w: usize,
    h: usize,
    cells: Vec<Vec<u32>>,
    rows: Vec<Vec<u32>>,
}

impl PolyIndex {
    /// Worth it only past this many points.
    const MIN_POINTS: usize = 32;

    fn new(pts: &[V], r: f64) -> PolyIndex {
        let mut bb = (f64::MAX, f64::MAX, f64::MIN, f64::MIN);
        for p in pts {
            bb = (bb.0.min(p.0), bb.1.min(p.1), bb.2.max(p.0), bb.3.max(p.1));
        }
        let cell = 1.0 * NM;
        let w = ((bb.2 - bb.0) / cell) as usize + 1;
        let h = ((bb.3 - bb.1) / cell) as usize + 1;
        let mut ix = PolyIndex { pts: pts.to_vec(), r, x0: bb.0, y0: bb.1, cell, w, h, cells: vec![Vec::new(); w * h], rows: vec![Vec::new(); h] };
        let n = pts.len();
        for e in 0..n {
            let (p, q) = (pts[e], pts[(e + 1) % n]);
            let (j0, j1) = (ix.row(p.1.min(q.1)), ix.row(p.1.max(q.1)));
            for j in j0..=j1 {
                ix.rows[j].push(e as u32);
                // Only the columns this edge reaches inside this row: long diagonals are common.
                let (ya, yb) = (ix.y0 + j as f64 * cell, ix.y0 + (j + 1) as f64 * cell);
                let at = |y: f64| if (q.1 - p.1).abs() < 1.0 { p.0 } else { p.0 + (q.0 - p.0) * ((y - p.1) / (q.1 - p.1)).clamp(0.0, 1.0) };
                let (xa, xb) = (at(ya), at(yb));
                let (i0, i1) = (ix.col(xa.min(xb).min(if j == j0 { p.0.min(q.0) } else { f64::MAX })), ix.col(xa.max(xb).max(if j == j0 { p.0.min(q.0) } else { f64::MIN })));
                let (i0, i1) = if (q.1 - p.1).abs() < 1.0 { (ix.col(p.0.min(q.0)), ix.col(p.0.max(q.0))) } else { (i0, i1) };
                for i in i0..=i1 {
                    ix.cells[j * w + i].push(e as u32);
                }
            }
        }
        ix
    }

    #[inline]
    fn col(&self, x: f64) -> usize {
        (((x - self.x0) / self.cell).max(0.0) as usize).min(self.w - 1)
    }
    #[inline]
    fn row(&self, y: f64) -> usize {
        (((y - self.y0) / self.cell).max(0.0) as usize).min(self.h - 1)
    }

    fn inside(&self, p: V) -> bool {
        if p.1 < self.y0 || p.1 > self.y0 + self.h as f64 * self.cell {
            return false;
        }
        let n = self.pts.len();
        let mut inside = false;
        for &e in &self.rows[self.row(p.1)] {
            let (a, b) = (self.pts[e as usize], self.pts[(e as usize + 1) % n]);
            if (a.1 > p.1) != (b.1 > p.1) && p.0 < (b.0 - a.0) * (p.1 - a.1) / (b.1 - a.1) + a.0 {
                inside = !inside;
            }
        }
        inside
    }

    fn dist_within(&self, a: V, b: V, limit: f64) -> f64 {
        if self.inside(a) || self.inside(b) {
            return 0.0;
        }
        let reach = limit + self.r;
        let (i0, i1) = (self.col(a.0.min(b.0) - reach), self.col(a.0.max(b.0) + reach));
        let (j0, j1) = (self.row(a.1.min(b.1) - reach), self.row(a.1.max(b.1) + reach));
        let n = self.pts.len();
        let mut d = limit + self.r;
        for j in j0..=j1 {
            for i in i0..=i1 {
                for &e in &self.cells[j * self.w + i] {
                    let (p, q) = (self.pts[e as usize], self.pts[(e as usize + 1) % n]);
                    d = d.min(seg_seg(a, b, p, q));
                    if d <= 0.0 {
                        return 0.0;
                    }
                }
            }
        }
        (d - self.r).max(0.0)
    }
}

pub struct Grid {
    x0: f64,
    y0: f64,
    cell: f64,
    w: usize,
    h: usize,
    cells: Vec<Vec<u32>>,
    items: Vec<Option<Obs>>,
    stamp: Vec<Cell<u32>>,
    tick: Cell<u32>,
    free: Vec<u32>,
    /// What changed since the coarse map last looked: (added, item).
    log: RefCell<Vec<(bool, Obs)>>,
}

impl Grid {
    pub fn new(bb: (f64, f64, f64, f64), cell: f64) -> Grid {
        let pad = 5.0 * NM;
        let x0 = bb.0 - pad;
        let y0 = bb.1 - pad;
        let w = (((bb.2 - bb.0 + 2.0 * pad) / cell).ceil() as usize + 1).max(1);
        let h = (((bb.3 - bb.1 + 2.0 * pad) / cell).ceil() as usize + 1).max(1);
        Grid { x0, y0, cell, w, h, cells: vec![Vec::new(); w * h], items: Vec::new(), stamp: Vec::new(), tick: Cell::new(0), free: Vec::new(), log: RefCell::new(Vec::new()) }
    }

    fn range(&self, bb: (f64, f64, f64, f64)) -> (usize, usize, usize, usize) {
        let cl = |v: f64, n: usize| -> usize { (v.max(0.0) as usize).min(n - 1) };
        (
            cl((bb.0 - self.x0) / self.cell, self.w),
            cl((bb.1 - self.y0) / self.cell, self.h),
            cl((bb.2 - self.x0) / self.cell, self.w),
            cl((bb.3 - self.y0) / self.cell, self.h),
        )
    }

    pub fn add(&mut self, shape: Shape, net: NetId, owner: Owner) -> u32 {
        let bb = shape.bbox();
        let fast = match &shape {
            Shape::Poly { pts, r } if pts.len() >= PolyIndex::MIN_POINTS => Some(Rc::new(PolyIndex::new(pts, *r))),
            _ => None,
        };
        self.log.borrow_mut().push((true, Obs { shape: shape.clone(), net, owner, bb, fast: fast.clone() }));
        let id = if let Some(i) = self.free.pop() {
            self.items[i as usize] = Some(Obs { shape, net, owner, bb, fast: fast.clone() });
            i
        } else {
            self.items.push(Some(Obs { shape, net, owner, bb, fast: fast.clone() }));
            self.stamp.push(Cell::new(0));
            (self.items.len() - 1) as u32
        };
        let (a, b, c, d) = self.range(bb);
        for j in b..=d {
            for i in a..=c {
                self.cells[j * self.w + i].push(id);
            }
        }
        id
    }

    pub fn remove(&mut self, id: u32) {
        let Some(o) = self.items[id as usize].take() else { return };
        self.log.borrow_mut().push((false, o.clone()));
        let (a, b, c, d) = self.range(o.bb);
        for j in b..=d {
            for i in a..=c {
                self.cells[j * self.w + i].retain(|&x| x != id);
            }
        }
        self.free.push(id);
    }

    /// Visit every item whose box meets `bb`. The visitor returns false to stop early.
    #[inline]
    pub fn visit(&self, bb: (f64, f64, f64, f64), mut f: impl FnMut(&Obs) -> bool) {
        let t = self.tick.get().wrapping_add(1);
        self.tick.set(t);
        let (a, b, c, d) = self.range(bb);
        for j in b..=d {
            for i in a..=c {
                for &id in &self.cells[j * self.w + i] {
                    let s = &self.stamp[id as usize];
                    if s.get() == t {
                        continue;
                    }
                    s.set(t);
                    let o = self.items[id as usize].as_ref().unwrap();
                    if o.bb.0 > bb.2 || o.bb.2 < bb.0 || o.bb.1 > bb.3 || o.bb.3 < bb.1 {
                        continue;
                    }
                    if !f(o) {
                        return;
                    }
                }
            }
        }
    }

    pub fn drain_log(&self) -> Vec<(bool, Obs)> {
        std::mem::take(&mut *self.log.borrow_mut())
    }

    pub fn iter(&self) -> impl Iterator<Item = &Obs> {
        self.items.iter().flatten()
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    /// The index must agree with the exact distance wherever it matters, on a polygon with long
    /// edges of every slope and a notch, the shape a real pour has.
    #[test]
    fn indexed_distance_agrees_with_exact() {
        let mut pts = vec![V(0.0, 0.0), V(20.0, 3.0), V(23.0, 18.0), V(12.0, 9.0), V(9.0, 21.0), V(-2.0, 13.0)];
        // Enough points to be worth an index: subdivide every edge.
        let mut fine = Vec::new();
        for i in 0..pts.len() {
            let (a, b) = (pts[i], pts[(i + 1) % pts.len()]);
            for k in 0..8 {
                let t = k as f64 / 8.0;
                fine.push(V(a.0 + (b.0 - a.0) * t, a.1 + (b.1 - a.1) * t));
            }
        }
        pts = fine.iter().map(|p| V(p.0 * NM, p.1 * NM)).collect();
        let shape = Shape::Poly { pts: pts.clone(), r: 0.0 };
        let ix = PolyIndex::new(&pts, 0.0);
        let limit = 0.8 * NM;
        let mut x = 12345u64;
        let mut rnd = |m: f64| {
            x = x.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
            (x >> 33) as f64 / (1u64 << 31) as f64 * m
        };
        for _ in 0..20000 {
            let a = V(rnd(28.0 * NM) - 3.0 * NM, rnd(28.0 * NM) - 3.0 * NM);
            let b = V(a.0 + rnd(6.0 * NM) - 3.0 * NM, a.1 + rnd(6.0 * NM) - 3.0 * NM);
            let exact = shape.dist_seg(a, b);
            let fast = ix.dist_within(a, b, limit);
            if exact < limit {
                assert!((exact - fast).abs() < 1.0, "{a:?} {b:?}: exact {exact}, indexed {fast}");
            } else {
                assert!(fast >= limit - 1.0, "{a:?} {b:?}: exact {exact}, indexed {fast}");
            }
        }
    }
}