KiCad - the KiCad Bridge

//! 3D bodies (the EDA Bridge Standard's `bodies` companion): each part's 3D model as a box, its
//! footprint rectangle on the board and its height range above the surface it sits on, in mm.
//!
//! KiCad keeps no heights in the board: they are in the model files a footprint names. This reads a
//! STEP file's topological vertices and its circles (each circle's real extent in its own plane, so a
//! radial capacitor is not made a radius taller) or a VRML file's points (KiCad's 0.1 inch unit), then
//! places the box the way KiCad's 3D viewer does (`render_3d_opengl.cpp`):
//! model = T(offset) * Rz(-rz) * Ry(-ry) * Rx(-rx) * S(scale), in a y-up space;
//! board = T(position) * Rz(orientation) * [bottom: (x, y, z) -> (x, -y, -z)] * model.
//! Nothing here talks to KiCad; the caller resolves the model paths (`${KICAD10_3DMODEL_DIR}` etc.).

use std::collections::HashMap;
use std::path::{Path, PathBuf};

use serde_json::{json, Value};

use crate::model_check::{parse_sexpr, Sx};

/// An axis-aligned box: min and max corners, mm.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Bounds {
    pub min: [f64; 3],
    pub max: [f64; 3],
}

impl Bounds {
    fn empty() -> Self {
        Bounds { min: [f64::MAX; 3], max: [f64::MIN; 3] }
    }
    fn add(&mut self, p: [f64; 3]) {
        for i in 0..3 {
            self.min[i] = self.min[i].min(p[i]);
            self.max[i] = self.max[i].max(p[i]);
        }
    }
    fn is_empty(&self) -> bool {
        self.min[0] > self.max[0]
    }
    fn corners(&self) -> Vec<[f64; 3]> {
        let mut out = Vec::with_capacity(8);
        for &x in &[self.min[0], self.max[0]] {
            for &y in &[self.min[1], self.max[1]] {
                for &z in &[self.min[2], self.max[2]] {
                    out.push([x, y, z]);
                }
            }
        }
        out
    }
}

// -- STEP ---------------------------------------------------------------------------------------

/// The numbers in a STEP argument list such as `'',(0.8,0.4,0.)`.
fn numbers(s: &str) -> Vec<f64> {
    s.split(|c: char| c == ',' || c == '(' || c == ')' || c.is_whitespace()).filter_map(|t| t.parse::<f64>().ok()).collect()
}

/// The `#n` references in a STEP argument list, in order.
fn refs(s: &str) -> Vec<u64> {
    s.split(|c: char| !(c.is_ascii_digit() || c == '#')).filter_map(|t| t.strip_prefix('#')).filter_map(|t| t.parse().ok()).collect()
}

/// The bounds of a STEP model in mm: its vertices, and every circular edge swept between its two
/// vertices in its own plane (a cylinder's side bulges past its vertices; a trimmed arc of a larger
/// circle does not reach the whole circle). Falls back to every point when the file has no vertices.
/// Ignores assembly placements (KiCad's library models are one body in one frame).
pub fn step_bounds(text: &str) -> Option<Bounds> {
    let scale = if text.contains("CONVERSION_BASED_UNIT('INCH'") || text.contains("CONVERSION_BASED_UNIT ( 'INCH'") {
        25.4
    } else if text.contains("SI_UNIT($,.METRE.)") || text.contains("SI_UNIT ( $ , .METRE. )") {
        1000.0
    } else {
        1.0
    };
    let mut points: HashMap<u64, [f64; 3]> = HashMap::new();
    let mut dirs: HashMap<u64, [f64; 3]> = HashMap::new();
    let mut axes: HashMap<u64, (u64, Option<u64>, Option<u64>)> = HashMap::new();
    let mut vertices: HashMap<u64, u64> = HashMap::new();
    let mut circles: HashMap<u64, (u64, f64)> = HashMap::new();
    // (start vertex, end vertex, curve, same sense as the curve)
    let mut edges: Vec<(u64, u64, u64, bool)> = Vec::new();
    // statements end at ';' outside a quoted string
    let mut stmt = String::new();
    let mut quoted = false;
    let data = text.find("DATA;").map(|i| &text[i + 5..]).unwrap_or(text);
    for c in data.chars() {
        if c == '\'' {
            quoted = !quoted;
        }
        if c == ';' && !quoted {
            let st = stmt.trim();
            if let Some(rest) = st.strip_prefix('#') {
                if let Some((id, body)) = rest.split_once('=') {
                    if let Ok(id) = id.trim().parse::<u64>() {
                        let body = body.trim();
                        let (kind, args) = body.split_once('(').map(|(k, a)| (k.trim(), a)).unwrap_or((body, ""));
                        // skip the entity's name string before reading numbers
                        let after_name = args.split_once("',").map(|(_, r)| r).unwrap_or(args);
                        match kind {
                            "CARTESIAN_POINT" => {
                                let n = numbers(after_name);
                                if n.len() >= 3 {
                                    points.insert(id, [n[0], n[1], n[2]]);
                                }
                            }
                            "DIRECTION" => {
                                let n = numbers(after_name);
                                if n.len() >= 3 {
                                    dirs.insert(id, [n[0], n[1], n[2]]);
                                }
                            }
                            "AXIS2_PLACEMENT_3D" => {
                                let r = refs(after_name);
                                if let Some(&loc) = r.first() {
                                    axes.insert(id, (loc, r.get(1).copied(), r.get(2).copied()));
                                }
                            }
                            "VERTEX_POINT" => {
                                if let Some(&p) = refs(after_name).first() {
                                    vertices.insert(id, p);
                                }
                            }
                            "CIRCLE" => {
                                let r = refs(after_name);
                                let rad = after_name.rsplit(',').next().and_then(|t| t.trim().trim_end_matches(')').trim().parse::<f64>().ok());
                                if let (Some(&ax), Some(rad)) = (r.first(), rad) {
                                    circles.insert(id, (ax, rad));
                                }
                            }
                            "EDGE_CURVE" => {
                                let r = refs(after_name);
                                if r.len() >= 3 {
                                    edges.push((r[0], r[1], r[2], !after_name.contains(".F.")));
                                }
                            }
                            _ => {}
                        }
                    }
                }
            }
            stmt.clear();
            continue;
        }
        stmt.push(c);
    }
    let mut b = Bounds::empty();
    for p in vertices.values() {
        if let Some(p) = points.get(p) {
            b.add(*p);
        }
    }
    let sub = |a: [f64; 3], c: [f64; 3]| [a[0] - c[0], a[1] - c[1], a[2] - c[2]];
    let dot = |a: [f64; 3], c: [f64; 3]| a[0] * c[0] + a[1] * c[1] + a[2] * c[2];
    let unit = |a: [f64; 3]| { let l = dot(a, a).sqrt().max(1e-12); [a[0] / l, a[1] / l, a[2] / l] };
    for (v1, v2, curve, same) in &edges {
        let Some((ax, r)) = circles.get(curve) else { continue };
        let Some((loc, dir, refd)) = axes.get(ax) else { continue };
        let Some(c) = points.get(loc).copied() else { continue };
        // the circle's own frame: x along its reference direction, y = normal x x
        let n = unit(dir.and_then(|d| dirs.get(&d)).copied().unwrap_or([0.0, 0.0, 1.0]));
        let mut x = refd.and_then(|d| dirs.get(&d)).copied().unwrap_or(if n[2].abs() < 0.9 { [0.0, 0.0, 1.0] } else { [1.0, 0.0, 0.0] });
        let k = dot(x, n);
        x = unit([x[0] - k * n[0], x[1] - k * n[1], x[2] - k * n[2]]);
        let y = [n[1] * x[2] - n[2] * x[1], n[2] * x[0] - n[0] * x[2], n[0] * x[1] - n[1] * x[0]];
        let angle = |p: [f64; 3]| { let d = sub(p, c); dot(d, y).atan2(dot(d, x)) };
        let (p1, p2) = (vertices.get(v1).and_then(|p| points.get(p)).copied(), vertices.get(v2).and_then(|p| points.get(p)).copied());
        // the edge runs counterclockwise about the normal from its start to its end when it has the
        // curve's sense, else from its end to its start; one vertex twice is the whole circle
        let (from, sweep) = match (p1, p2) {
            (Some(a), Some(e)) if v1 != v2 => {
                let (s0, s1) = if *same { (angle(a), angle(e)) } else { (angle(e), angle(a)) };
                (s0, (s1 - s0).rem_euclid(std::f64::consts::TAU))
            }
            _ => (0.0, std::f64::consts::TAU),
        };
        let steps = ((sweep.to_degrees() / 2.0).ceil() as usize).max(1);
        for i in 0..=steps {
            let t = from + sweep * i as f64 / steps as f64;
            let (sn, cs) = t.sin_cos();
            b.add([c[0] + r * (cs * x[0] + sn * y[0]), c[1] + r * (cs * x[1] + sn * y[1]), c[2] + r * (cs * x[2] + sn * y[2])]);
        }
    }
    if b.is_empty() {
        for p in points.values() {
            b.add(*p);
        }
    }
    if b.is_empty() {
        return None;
    }
    for i in 0..3 {
        b.min[i] *= scale;
        b.max[i] *= scale;
    }
    Some(b)
}

/// The bounds of a VRML model in mm: every `point [ ... ]` list, in KiCad's VRML unit (0.1 inch).
pub fn wrl_bounds(text: &str) -> Option<Bounds> {
    let mut b = Bounds::empty();
    let mut rest = text;
    while let Some(i) = rest.find("point") {
        rest = &rest[i + 5..];
        let Some(open) = rest.find('[') else { break };
        if !rest[..open].trim().is_empty() {
            continue;
        }
        let Some(close) = rest[open..].find(']') else { break };
        let n = numbers(&rest[open + 1..open + close]);
        for t in n.chunks_exact(3) {
            b.add([t[0] * 2.54, t[1] * 2.54, t[2] * 2.54]);
        }
        rest = &rest[open + close..];
    }
    (!b.is_empty()).then_some(b)
}

// -- placement ----------------------------------------------------------------------------------

/// One `(model ...)` of a footprint.
#[derive(Clone, Debug, PartialEq)]
pub struct ModelRef {
    pub path: String,
    pub offset: [f64; 3],
    pub scale: [f64; 3],
    pub rotate: [f64; 3],
    pub hidden: bool,
}

fn xyz(m: &Sx, tag: &str, default: [f64; 3]) -> [f64; 3] {
    let Some(x) = m.find(tag).and_then(|t| t.find("xyz")) else { return default };
    let n = |i| x.atom_at(i).and_then(|a| a.parse::<f64>().ok());
    match (n(1), n(2), n(3)) {
        (Some(a), Some(b), Some(c)) => [a, b, c],
        _ => default,
    }
}

pub fn model_refs(fp: &Sx) -> Vec<ModelRef> {
    fp.find_all("model").map(|m| ModelRef {
        path: m.value().unwrap_or("").to_string(),
        offset: xyz(m, "offset", [0.0; 3]),
        scale: xyz(m, "scale", [1.0; 3]),
        rotate: xyz(m, "rotate", [0.0; 3]),
        hidden: m.find("hide").map(|h| h.value() != Some("no")).unwrap_or(false) || m.as_list().unwrap_or(&[]).iter().any(|a| a.as_atom() == Some("hide")),
    }).collect()
}

fn rot(p: [f64; 3], axis: usize, deg: f64) -> [f64; 3] {
    let (s, c) = deg.to_radians().sin_cos();
    let (i, j) = match axis { 0 => (1, 2), 1 => (2, 0), _ => (0, 1) };
    let mut q = p;
    q[i] = p[i] * c - p[j] * s;
    q[j] = p[i] * s + p[j] * c;
    q
}

/// The model's box in the footprint's own frame (y up, z up from the surface): KiCad's model transform.
pub fn in_footprint(b: &Bounds, m: &ModelRef) -> Bounds {
    let mut out = Bounds::empty();
    for p in b.corners() {
        let mut q = [p[0] * m.scale[0], p[1] * m.scale[1], p[2] * m.scale[2]];
        q = rot(q, 0, -m.rotate[0]);
        q = rot(q, 1, -m.rotate[1]);
        q = rot(q, 2, -m.rotate[2]);
        out.add([q[0] + m.offset[0], q[1] + m.offset[1], q[2] + m.offset[2]]);
    }
    out
}

/// A footprint-frame box placed on the board (y up): its footprint rectangle and its height range
/// above the surface it sits on.
pub fn on_board(b: &Bounds, at: (f64, f64), angle: f64, bottom: bool) -> ([f64; 4], [f64; 2]) {
    let mut x = (f64::MAX, f64::MIN);
    let mut y = (f64::MAX, f64::MIN);
    for p in b.corners() {
        let q = if bottom { [p[0], -p[1], -p[2]] } else { p };
        let r = rot(q, 2, angle);
        x = (x.0.min(r[0] + at.0), x.1.max(r[0] + at.0));
        y = (y.0.min(r[1] + at.1), y.1.max(r[1] + at.1));
    }
    ([x.0, y.0, x.1, y.1], [b.min[2], b.max[2]])
}

/// The footprint rectangle of a footprint-frame box placed on the board (y up), as its four turned corners: what a
/// part not at a multiple of 90° really covers (its `on_board` extent is the bounding box of these).
pub fn on_board_corners(b: &Bounds, at: (f64, f64), angle: f64, bottom: bool) -> [[f64; 2]; 4] {
    let (x0, y0, x1, y1) = (b.min[0], b.min[1], b.max[0], b.max[1]);
    [[x0, y0], [x1, y0], [x1, y1], [x0, y1]].map(|[x, y]| {
        let q = if bottom { [x, -y, 0.0] } else { [x, y, 0.0] };
        let r = rot(q, 2, angle);
        [r[0] + at.0, r[1] + at.1]
    })
}

// -- paths --------------------------------------------------------------------------------------

/// Expand `${VAR}` in a model path: the 3D model dirs of every KiCad version name the install's
/// 3dmodels folder, KIPRJMOD the board's folder, anything else the variables given (KiCad's own
/// settings, then the process environment).
pub fn expand(path: &str, models_dir: Option<&Path>, project_dir: &Path, vars: &HashMap<String, String>) -> Option<String> {
    let mut out = String::new();
    let mut rest = path;
    while let Some(i) = rest.find("${") {
        out.push_str(&rest[..i]);
        let after = &rest[i + 2..];
        let end = after.find('}')?;
        let name = &after[..end];
        let val = if name == "KIPRJMOD" {
            Some(project_dir.to_string_lossy().to_string())
        } else if let Some(v) = vars.get(name) {
            Some(v.clone())
        } else if (name.ends_with("_3DMODEL_DIR") && name.starts_with("KICAD")) || name == "KISYS3DMOD" {
            models_dir.map(|d| d.to_string_lossy().to_string())
        } else {
            std::env::var(name).ok()
        };
        out.push_str(&val?);
        rest = &after[end + 1..];
    }
    out.push_str(rest);
    Some(out)
}

/// The file KiCad would load for a model path: the path itself, or the same name as STEP when a
/// VRML is named but only the STEP ships (KiCad 10's library), or the other way round.
pub fn find_model(expanded: &str, project_dir: &Path) -> Option<PathBuf> {
    let p = PathBuf::from(expanded);
    let p = if p.is_relative() { project_dir.join(p) } else { p };
    if p.is_file() {
        return Some(p);
    }
    let alts: &[&str] = match p.extension().and_then(|e| e.to_str()).map(str::to_ascii_lowercase).as_deref() {
        Some("wrl") | Some("wrz") => &["step", "stp", "STEP", "STP"],
        Some("step") | Some("stp") => &["wrl", "WRL"],
        _ => &["step", "stp", "wrl"],
    };
    alts.iter().map(|e| p.with_extension(e)).find(|q| q.is_file())
}

/// Bounds of one model file, by its extension.
pub fn file_bounds(path: &Path) -> Result<Bounds, String> {
    let ext = path.extension().and_then(|e| e.to_str()).unwrap_or("").to_ascii_lowercase();
    let bytes = std::fs::read(path).map_err(|e| format!("could not read {}: {e}", path.display()))?;
    let text = String::from_utf8_lossy(&bytes);
    let b = match ext.as_str() {
        "step" | "stp" => step_bounds(&text),
        "wrl" => wrl_bounds(&text),
        other => return Err(format!("a .{other} model is not read (STEP and VRML are)")),
    };
    b.ok_or_else(|| format!("no geometry found in {}", path.display()))
}

/// Every part's bodies as the standard's boxes, from the board text. `load` returns a model file's
/// bounds (cached by the caller) or why it could not.
pub fn bodies(board_text: &str, y_up: bool, resolve: &mut dyn FnMut(&str) -> Result<(Bounds, String), String>) -> Result<Value, String> {
    let root = parse_sexpr(board_text).ok_or("empty board text")?;
    let mut out = Vec::new();
    let mut missing = Vec::new();
    let mut parts = 0;
    for fp in root.find_all("footprint") {
        parts += 1;
        let reference = fp.find_all("property").find(|p| p.value() == Some("Reference")).and_then(|p| p.atom_at(2))
            .or_else(|| fp.find_all("fp_text").find(|t| t.value() == Some("reference")).and_then(|t| t.atom_at(2)))
            .unwrap_or("").to_string();
        let at = fp.find("at");
        let num = |i| at.and_then(|a| a.atom_at(i)).and_then(|s| s.parse::<f64>().ok());
        let (x, y, angle) = (num(1).unwrap_or(0.0), -num(2).unwrap_or(0.0), num(3).unwrap_or(0.0));
        let bottom = fp.find("layer").and_then(Sx::value) == Some("B.Cu");
        let side = if bottom { "bottom" } else { "top" };
        let models: Vec<ModelRef> = model_refs(fp).into_iter().filter(|m| !m.hidden).collect();
        if models.is_empty() {
            missing.push(json!({"ref": reference, "why": "no 3D model"}));
            continue;
        }
        for m in models {
            if m.path.starts_with("kicad-embed://") {
                missing.push(json!({"ref": reference, "model": m.path, "why": "an embedded model (kicad-embed://) is not read yet"}));
                continue;
            }
            match resolve(&m.path) {
                Ok((b, file)) => {
                    let local = in_footprint(&b, &m);
                    let (r, z) = on_board(&local, (x, y), angle, bottom);
                    let fy = |v: f64| if y_up { v } else { -v };
                    let (y0, y1) = if y_up { (r[1], r[3]) } else { (fy(r[3]), fy(r[1])) };
                    let r6 = |v: f64| (v * 1e4).round() / 1e4 + 0.0;
                    let mut row = json!({"ref": reference, "side": side, "x0": r6(r[0]), "y0": r6(y0), "x1": r6(r[2]), "y1": r6(y1),
                        "z0": r6(z[0]), "z1": r6(z[1]), "heightMm": r6(z[1].max(0.0)), "model": m.path, "file": file});
                    // a part turned off the 90° grid: its real rectangle, so the checks don't compare the bigger bounding
                    // box (standard #26)
                    if (angle / 90.0 - (angle / 90.0).round()).abs() > 1e-6 {
                        let cs = on_board_corners(&local, (x, y), angle, bottom);
                        row["corners"] = json!(cs.iter().map(|c| [r6(c[0]), r6(fy(c[1]))]).collect::<Vec<_>>());
                    }
                    out.push(row);
                }
                Err(why) => missing.push(json!({"ref": reference, "model": m.path, "why": why})),
            }
        }
    }
    Ok(json!({"bodies": out, "noBody": missing, "parts": parts}))
}

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

    const STEP: &str = "ISO-10303-21;\nDATA;\n#1=CARTESIAN_POINT('',(-0.8,-0.4,0.));\n#2=VERTEX_POINT('',#1);\n#3=CARTESIAN_POINT('',(0.8,0.4,0.45));\n#4=VERTEX_POINT('',#3);\n\
        #5=CARTESIAN_POINT('',(9.,9.,9.));\n#6=CARTESIAN_POINT('Origin',(0.,0.,5.));\n#7=DIRECTION('',(0.,0.,1.));\n#8=DIRECTION('',(1.,0.,0.));\n#9=AXIS2_PLACEMENT_3D('',#6,#7,#8);\n#10=CIRCLE('',#9,1.5);\n\
        #11=CARTESIAN_POINT('',(1.5,0.,5.));\n#12=VERTEX_POINT('',#11);\n#13=EDGE_CURVE('',#12,#12,#10,.T.);\n\
        #14=CIRCLE('',#9,9.);\n#15=CARTESIAN_POINT('',(9.,0.,5.));\n#16=VERTEX_POINT('',#15);\n#17=CARTESIAN_POINT('',(0.,9.,5.));\n#18=VERTEX_POINT('',#17);\n#19=EDGE_CURVE('',#18,#16,#14,.F.);\nENDSEC;\n";

    #[test]
    fn step_bounds_are_vertices_and_circles_in_their_plane() {
        let b = step_bounds(STEP).unwrap();
        // the stray control point (9, 9, 9) is not a vertex; the full circle at z 5 of radius 1.5 lies
        // flat; the quarter arc of radius 9 from (9, 0) to (0, 9) (an edge against the curve's sense,
        // so the short way) reaches no further than its ends
        assert!((b.min[0] + 1.5).abs() < 1e-9 && (b.min[1] + 1.5).abs() < 1e-9 && b.min[2] == 0.0, "{b:?}");
        assert!((b.max[0] - 9.0).abs() < 1e-9 && (b.max[1] - 9.0).abs() < 1e-9 && b.max[2] == 5.0, "{b:?}");
        let inch = STEP.replace("DATA;", "CONVERSION_BASED_UNIT('INCH',#99);\nDATA;");
        assert!((step_bounds(&inch).unwrap().max[2] - 127.0).abs() < 1e-6);
    }

    #[test]
    fn wrl_points_are_in_tenths_of_an_inch() {
        let b = wrl_bounds("#VRML V2.0 utf8\nShape { geometry IndexedFaceSet { coord Coordinate { point [ 0 0 0, 1 2 0.5 ] } } }").unwrap();
        assert_eq!(b.max, [2.54, 5.08, 1.27]);
    }

    #[test]
    fn a_body_is_placed_turned_and_flipped_like_kicads_3d_viewer() {
        // a 2 x 1 x 0.5 body offset 1 mm to +x in its footprint
        let b = Bounds { min: [-1.0, -0.5, 0.0], max: [1.0, 0.5, 0.5] };
        let m = ModelRef { path: String::new(), offset: [1.0, 0.0, 0.0], scale: [1.0; 3], rotate: [0.0; 3], hidden: false };
        let local = in_footprint(&b, &m);
        assert_eq!(local.min, [0.0, -0.5, 0.0]);
        let (r, z) = on_board(&local, (10.0, -10.0), 90.0, false);
        // turned 90 counterclockwise: the body now reaches from the origin up +y
        let r: Vec<f64> = r.iter().map(|v| (v * 1e9).round() / 1e9).collect();
        assert_eq!(r, vec![9.5, -10.0, 10.5, -8.0]);
        assert_eq!(z, [0.0, 0.5]);
        // on the bottom the footprint is mirrored top to bottom in its own frame: an offset in +y flips
        let m2 = ModelRef { offset: [0.0, 1.0, 0.0], ..m };
        let (r, _) = on_board(&in_footprint(&b, &m2), (0.0, 0.0), 0.0, true);
        assert_eq!(r, [-1.0, -1.5, 1.0, -0.5]);
    }

    #[test]
    fn model_paths_expand_like_kicad_and_fall_back_to_the_step() {
        let vars = HashMap::from([("MYLIB".to_string(), "D:/libs".to_string())]);
        let models = Path::new("C:/KiCad/share/kicad/3dmodels");
        let proj = Path::new("C:/p");
        assert_eq!(expand("${KICAD8_3DMODEL_DIR}/R.3dshapes/R.wrl", Some(models), proj, &vars).unwrap(), "C:/KiCad/share/kicad/3dmodels/R.3dshapes/R.wrl");
        assert_eq!(expand("${KIPRJMOD}/m.step", None, proj, &vars).unwrap(), "C:/p/m.step");
        assert_eq!(expand("${MYLIB}/x.step", None, proj, &vars).unwrap(), "D:/libs/x.step");
        assert!(expand("${NOPE_NOT_SET_ANYWHERE}/x.step", None, proj, &vars).is_none());
        let dir = std::env::temp_dir().join(format!("kb-bodies-{}", std::process::id()));
        std::fs::create_dir_all(&dir).unwrap();
        std::fs::write(dir.join("R.step"), STEP).unwrap();
        assert_eq!(find_model(&dir.join("R.wrl").to_string_lossy(), &dir), Some(dir.join("R.step")));
        std::fs::remove_dir_all(&dir).ok();
    }

    #[test]
    fn a_part_turned_45_degrees_covers_its_rectangle_not_the_bounding_box() {
        let mut b = Bounds::empty();
        b.add([-1.0, -0.5, 0.0]);
        b.add([1.0, 0.5, 0.6]);
        let cs = on_board_corners(&b, (10.0, 20.0), 45.0, false);
        let (r, _) = on_board(&b, (10.0, 20.0), 45.0, false);
        // the corners lie inside the bounding box and touch each of its sides
        for c in cs {
            assert!(c[0] >= r[0] - 1e-9 && c[0] <= r[2] + 1e-9 && c[1] >= r[1] - 1e-9 && c[1] <= r[3] + 1e-9);
        }
        // the rectangle keeps its 2 x 1 size: two neighbouring corners are 2 mm apart, the next pair 1 mm
        let d = |a: [f64; 2], b: [f64; 2]| ((a[0] - b[0]).powi(2) + (a[1] - b[1]).powi(2)).sqrt();
        assert!((d(cs[0], cs[1]) - 2.0).abs() < 1e-9 && (d(cs[1], cs[2]) - 1.0).abs() < 1e-9);
        // the bounding box is larger than 2 x 1 at 45°
        assert!(r[2] - r[0] > 2.0 && r[3] - r[1] > 1.0);
    }
}