//! Extrude & Bevel, Inflate and Rotate (a flat plane) as projected, shaded vector facets.
//!
//! The 2D art is the front of the solid: its outline is extruded backwards along the view's
//! z-axis, about the centre of its bounds, at the art's own position and size. Planar faces (caps
//! and the rotated plane) are drawn whole, holes included, so no triangulation is needed; curved
//! faces (bevels, side walls under twist or taper, inflated surfaces) are cut into small facets
//! and flat shaded, like Revolve.
//!
//! Visibility is a painter's algorithm, as for Revolve: faces of a closed solid that face away
//! from the viewer are left out, the rest are drawn back to front.

use vectorcraft_geom::{BezPath, ParamCurve, PathData, Point, Rect};

use crate::{MAX_FACES, MAX_PROFILE_POINTS, Revolve, V3};

/// The view shared by every 3D effect: rotation, perspective and one directional light.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct View {
    pub rotation_x: f64,
    pub rotation_y: f64,
    pub rotation_z: f64,
    /// Perspective strength, 0 for orthographic projection, up to 1.
    pub perspective: f64,
    pub light_azimuth: f64,
    pub light_elevation: f64,
    pub light_intensity: f64,
    pub ambient: f64,
    pub shade: bool,
}

impl Default for View {
    fn default() -> Self {
        let r = Revolve::default();
        Self {
            rotation_x: r.rotation_x,
            rotation_y: r.rotation_y,
            rotation_z: r.rotation_z,
            perspective: r.perspective,
            light_azimuth: r.light_azimuth,
            light_elevation: r.light_elevation,
            light_intensity: r.light_intensity,
            ambient: r.ambient,
            shade: r.shade,
        }
    }
}

impl View {
    /// The same rotation, projection and lighting maths as Revolve's.
    fn revolve(self) -> Revolve {
        Revolve {
            rotation_x: self.rotation_x,
            rotation_y: self.rotation_y,
            rotation_z: self.rotation_z,
            perspective: self.perspective,
            light_azimuth: self.light_azimuth,
            light_elevation: self.light_elevation,
            light_intensity: self.light_intensity,
            ambient: self.ambient,
            shade: self.shade,
            ..Revolve::default()
        }
        .sanitized()
    }
}

/// The shape of a bevel's cross-section, from the cap's edge (0) to the side wall (1).
#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
pub enum BevelShape {
    /// A straight 45° chamfer when width and height are equal.
    #[default]
    Classic,
    /// A convex quarter round.
    Rounded,
    /// A concave quarter round.
    Cove,
    /// Two flat steps.
    Step,
    /// A convex then concave S curve.
    Ogee,
}

impl BevelShape {
    pub const ALL: [BevelShape; 5] = [Self::Classic, Self::Rounded, Self::Cove, Self::Step, Self::Ogee];

    pub fn id(self) -> &'static str {
        match self {
            Self::Classic => "classic",
            Self::Rounded => "rounded",
            Self::Cove => "cove",
            Self::Step => "step",
            Self::Ogee => "ogee",
        }
    }

    pub fn parse(s: &str) -> Option<Self> {
        Self::ALL.into_iter().find(|b| b.id() == s)
    }

    /// The cross-section at `s` (0..1): (share of the width, share of the height), both 0..1.
    fn at(self, s: f64) -> (f64, f64) {
        let s = s.clamp(0.0, 1.0);
        let q = s * std::f64::consts::FRAC_PI_2;
        match self {
            Self::Classic => (s, s),
            Self::Rounded => (1.0 - q.cos(), q.sin()),
            Self::Cove => (q.sin(), 1.0 - q.cos()),
            Self::Step if s >= 1.0 => (1.0, 1.0),
            Self::Step => {
                // Down, across, down, across.
                let k = s * 4.0;
                let (i, f) = (k.floor(), k.fract());
                match i as u8 {
                    0 => (0.0, f * 0.5),
                    1 => (f * 0.5, 0.5),
                    2 => (0.5, 0.5 + f * 0.5),
                    _ => (0.5 + f * 0.5, 1.0),
                }
            }
            Self::Ogee => {
                let h = (1.0 - (s * std::f64::consts::PI).cos()) / 2.0;
                (s, h)
            }
        }
    }

    /// Samples along the cross-section: the step shape needs its corners exactly.
    fn samples(self) -> Vec<f64> {
        match self {
            Self::Classic => vec![0.0, 1.0],
            Self::Step => vec![0.0, 0.25, 0.5, 0.75, 1.0],
            _ => (0..=8).map(|i| f64::from(i) / 8.0).collect(),
        }
    }
}

/// A bevel along the front (and with `both_sides` the back) edge of an extrusion.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Bevel {
    pub shape: BevelShape,
    /// How far the bevel reaches across the face, in points.
    pub width: f64,
    /// How deep it reaches along the extrusion, in points.
    pub height: f64,
    /// How many times the cross-section repeats across the width.
    pub repeat: usize,
    /// The share (0..1) of the width left flat between repeats.
    pub space: f64,
    /// Carved into the shape (true) or added around it (false).
    pub inside: bool,
    pub both_sides: bool,
}

impl Default for Bevel {
    fn default() -> Self {
        Self { shape: BevelShape::Classic, width: 4.0, height: 4.0, repeat: 1, space: 0.0, inside: false, both_sides: false }
    }
}

/// Extrude & Bevel's options. Lengths are points, angles degrees.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Extrude {
    pub depth: f64,
    /// The turn of the back face against the front, about the shape's centre.
    pub twist: f64,
    /// The size of the back face, as a share (0.01..1) of the front's.
    pub taper: f64,
    /// Closed ends (a solid) or open ones (a hollow tube).
    pub cap: bool,
    pub bevel: Option<Bevel>,
}

impl Default for Extrude {
    fn default() -> Self {
        Self { depth: 50.0, twist: 0.0, taper: 1.0, cap: true, bevel: None }
    }
}

/// Inflate's options.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Inflate {
    /// The flat depth between the front and back, in points.
    pub depth: f64,
    /// How far the surface puffs out (0..100 %; 50 % rounds it like a half sphere).
    pub volume: f64,
    /// Puff the back out too; otherwise it stays flat.
    pub both_sides: bool,
}

impl Default for Inflate {
    fn default() -> Self {
        Self { depth: 0.0, volume: 50.0, both_sides: false }
    }
}

/// One projected, flat-shaded polygon. Holes are contours of their own (nonzero winding).
#[derive(Clone, Debug, PartialEq)]
pub struct Facet {
    pub contours: Vec<Vec<Point>>,
    pub depth: f64,
    pub brightness: f64,
}

impl Facet {
    /// About how thick the facet is across: twice its area over its perimeter (a w × l rectangle
    /// long against its width gives about w).
    pub fn thickness(&self) -> f64 {
        let area: f64 = self.contours.iter().map(|c| signed_area(c)).sum::<f64>().abs();
        let perimeter: f64 =
            self.contours.iter().map(|c| c.iter().zip(c.iter().cycle().skip(1)).take(c.len()).map(|(a, b)| a.distance(*b)).sum::<f64>()).sum();
        if perimeter > 0.0 { 2.0 * area / perimeter } else { 0.0 }
    }

    pub fn path(&self) -> PathData {
        let mut bp = BezPath::new();
        for c in &self.contours {
            let Some(first) = c.first() else { continue };
            bp.move_to(*first);
            for p in c.iter().skip(1) {
                bp.line_to(*p);
            }
            bp.close_path();
        }
        PathData::from_bezpath(&bp)
    }
}

fn finite(v: f64, fallback: f64, min: f64, max: f64) -> f64 {
    if v.is_finite() { v.clamp(min, max) } else { fallback }
}

impl Extrude {
    fn sanitized(self) -> Self {
        Self {
            depth: finite(self.depth, 50.0, 0.0, 2000.0),
            twist: finite(self.twist, 0.0, -360.0, 360.0),
            taper: finite(self.taper, 1.0, 0.01, 1.0),
            bevel: self.bevel.map(|b| Bevel {
                width: finite(b.width, 4.0, 0.0, 1000.0),
                height: finite(b.height, 4.0, 0.0, 1000.0),
                repeat: b.repeat.clamp(1, 32),
                space: finite(b.space, 0.0, 0.0, 0.9),
                ..b
            }),
            ..self
        }
    }
}

impl Inflate {
    fn sanitized(self) -> Self {
        Self { depth: finite(self.depth, 0.0, 0.0, 2000.0), volume: finite(self.volume, 50.0, 0.0, 100.0), ..self }
    }
}

/// The rings of a shape: closed polygons with the filled side on their left, (-dy, dx): outer
/// outlines run one way and holes the other, so nonzero winding fills them as the art is.
fn rings(path: &PathData) -> Result<Vec<Vec<Point>>, String> {
    let mut count = 0;
    let mut out = vec![];
    for sp in &path.subpaths {
        let segments = if sp.closed { sp.anchors.len() } else { sp.anchors.len().saturating_sub(1) };
        if segments > MAX_PROFILE_POINTS {
            return Err("The shape is too complex for a 3D effect".into());
        }
        let mut ring: Vec<Point> = vec![];
        if let Some(first) = sp.anchors.first() {
            ring.push(first.p);
        }
        for i in 0..segments {
            let curve = sp.segment(i);
            let samples: u32 = if sp.segment_is_line(i) { 1 } else { 16 };
            for k in 1..=samples {
                ring.push(curve.eval(f64::from(k) / f64::from(samples)));
            }
            count += samples as usize;
            if count > MAX_PROFILE_POINTS {
                return Err("The shape is too complex for a 3D effect".into());
            }
        }
        if ring.iter().any(|p| !p.x.is_finite() || !p.y.is_finite() || p.x.abs() > 1e8 || p.y.abs() > 1e8) {
            return Err("The shape contains invalid coordinates".into());
        }
        ring.dedup_by(|a, b| a.distance(*b) < 1e-9);
        while ring.len() > 1 && ring.first().zip(ring.last()).is_some_and(|(a, b)| a.distance(*b) < 1e-9) {
            ring.pop();
        }
        if ring.len() >= 3 && signed_area(&ring).abs() > 1e-9 {
            out.push(ring);
        }
    }
    // A ring inside an odd number of others is a hole.
    let holes: Vec<bool> = out
        .iter()
        .enumerate()
        .map(|(i, r)| {
            let p = r.first().copied().unwrap_or_default();
            out.iter().enumerate().filter(|(j, o)| *j != i && contains(o, p)).count() % 2 == 1
        })
        .collect();
    for (r, hole) in out.iter_mut().zip(holes) {
        if (signed_area(r) > 0.0) == hole {
            r.reverse();
        }
    }
    Ok(out)
}

fn signed_area(r: &[Point]) -> f64 {
    r.iter().zip(r.iter().cycle().skip(1)).take(r.len()).map(|(a, b)| a.x * b.y - b.x * a.y).sum::<f64>() / 2.0
}

fn contains(r: &[Point], p: Point) -> bool {
    let mut inside = false;
    for (a, b) in r.iter().zip(r.iter().cycle().skip(1)).take(r.len()) {
        if (a.y > p.y) != (b.y > p.y) && p.x < a.x + (p.y - a.y) / (b.y - a.y) * (b.x - a.x) {
            inside = !inside;
        }
    }
    inside
}

/// The unit normal on the left of `a`→`b`: towards the filled side.
fn left(a: Point, b: Point) -> (f64, f64) {
    let (dx, dy) = (b.x - a.x, b.y - a.y);
    let l = dx.hypot(dy);
    if l < 1e-12 { (0.0, 0.0) } else { (-dy / l, dx / l) }
}

/// `ring` moved `w` points towards its filled side (negative: away from it), each corner along
/// its bisector, its mitre limited to four times `w` so sharp corners don't shoot off. Each point
/// stays the partner of its source point, so rings at different insets join face to face.
fn inset(ring: &[Point], w: f64) -> Vec<Point> {
    if w == 0.0 {
        return ring.to_vec();
    }
    let n = ring.len();
    (0..n)
        .map(|i| {
            let p = ring[i];
            let prev = ring[(i + n - 1) % n];
            let next = ring[(i + 1) % n];
            let (ax, ay) = left(prev, p);
            let (bx, by) = left(p, next);
            let (mut mx, mut my) = (ax + bx, ay + by);
            let ml = mx.hypot(my);
            if ml < 1e-9 {
                (mx, my) = (ax, ay);
            } else {
                (mx, my) = (mx / ml, my / ml);
            }
            let cos = (mx * ax + my * ay).max(0.25);
            Point::new(p.x + mx * w / cos, p.y + my * w / cos)
        })
        .collect()
}

/// Turns 3D points into document points and shades faces for one view.
struct Camera {
    view: Revolve,
    origin: Point,
    distance: f64,
    light: V3,
}

impl Camera {
    fn new(view: View, bounds: Rect, depth: f64) -> Self {
        let view = view.revolve();
        let [lx, ly, lz] = view.light_direction();
        let size = (bounds.width() / 2.0).max(bounds.height() / 2.0).max(depth / 2.0).max(1.0);
        Self { view, origin: bounds.center(), distance: 4.0 * size, light: V3(lx, ly, lz) }
    }
    fn rotate(&self, v: V3) -> V3 {
        self.view.rotate(v)
    }
    /// A point of the art, `z` points in front of its plane, in view space.
    fn place(&self, p: Point, z: f64) -> V3 {
        self.rotate(V3(p.x - self.origin.x, p.y - self.origin.y, z))
    }
    fn project(&self, v: V3) -> Point {
        let f = self.distance / (self.distance - v.2 * self.view.perspective);
        Point::new(self.origin.x + v.0 * f, self.origin.y + v.1 * f)
    }
    /// The cosine between a face's outward unit normal `n` and the direction from `at` (both in
    /// view space) to the eye: positive when the face is turned towards the viewer.
    fn facing(&self, at: V3, n: V3) -> f64 {
        let to_eye = if self.view.perspective > 0.0 { V3(0.0, 0.0, self.distance / self.view.perspective).sub(at).unit() } else { V3(0.0, 0.0, 1.0) };
        n.dot(to_eye)
    }
    fn brightness(&self, n: V3) -> f64 {
        if self.view.shade { (self.view.ambient + self.view.light_intensity * n.dot(self.light).max(0.0)).clamp(0.0, 1.0) } else { 1.0 }
    }
}

/// Facets being collected for one solid.
struct Builder<'a> {
    camera: &'a Camera,
    /// A closed solid hides the faces turned away from the viewer; an open one shows their
    /// insides.
    closed: bool,
    facets: Vec<Facet>,
}

impl Builder<'_> {
    /// Add a face of view-space contours with outward normal `n` (view space).
    fn face(&mut self, contours: Vec<Vec<V3>>, n: V3) -> Result<(), String> {
        let count: usize = contours.iter().map(Vec::len).sum();
        if count < 3 {
            return Ok(());
        }
        let centre = contours.iter().flatten().fold(V3::default(), |s, v| V3(s.0 + v.0, s.1 + v.1, s.2 + v.2));
        let centre = V3(centre.0 / count as f64, centre.1 / count as f64, centre.2 / count as f64);
        let n = n.unit();
        let facing = self.camera.facing(centre, n);
        // A closed solid hides the faces turned away. One seen almost exactly edge-on is kept but
        // drawn first, under everything: on a curved surface rounding alone turns such faces
        // either way, and dropping every other one would leave gaps.
        let (n, behind) = if facing > 0.0 {
            (n, false)
        } else if self.closed && facing > -0.02 {
            (n, true)
        } else if self.closed {
            return Ok(());
        } else {
            (V3(-n.0, -n.1, -n.2), false)
        };
        let contours: Vec<Vec<Point>> = contours.iter().map(|c| c.iter().map(|v| self.camera.project(*v)).collect()).collect();
        if contours.iter().flatten().any(|p| !p.x.is_finite() || !p.y.is_finite()) {
            return Err("The 3D projection is not finite".into());
        }
        if contours.iter().map(|c| signed_area(c).abs()).sum::<f64>() < 1e-10 {
            return Ok(());
        }
        let depth = if behind { centre.2 - 4.0 * self.camera.distance } else { centre.2 };
        self.facets.push(Facet { contours, depth, brightness: self.camera.brightness(n) });
        if self.facets.len() > MAX_FACES {
            return Err("The 3D surface is too complex; simplify the shape".into());
        }
        Ok(())
    }

    /// The band between two partner rings (same length), as one quad per edge. `outward` is the
    /// direction (before rotation) the band faces away from the solid, for orienting normals.
    fn band(&mut self, a: &[V3], b: &[V3], outward: impl Fn(usize) -> V3) -> Result<(), String> {
        let n = a.len().min(b.len());
        for i in 0..n {
            let j = (i + 1) % n;
            let quad = vec![a[i], a[j], b[j], b[i]];
            let mut normal = quad[1].sub(quad[0]).cross(quad[2].sub(quad[0]));
            if normal.dot(normal) < 1e-18 {
                normal = quad[2].sub(quad[0]).cross(quad[3].sub(quad[0]));
            }
            if normal.dot(normal) < 1e-18 {
                continue;
            }
            if normal.dot(self.camera.rotate(outward(i))) < 0.0 {
                normal = V3(-normal.0, -normal.1, -normal.2);
            }
            self.face(vec![quad], normal)?;
        }
        Ok(())
    }

    fn finish(mut self) -> Vec<Facet> {
        self.facets.sort_by(|a, b| a.depth.total_cmp(&b.depth));
        self.facets
    }
}

/// The area-weighted normal of (possibly non-planar) contours, by Newell's method: for contours
/// winding with their filled side on the left as seen from in front, it points to the front.
fn newell(contours: &[Vec<V3>]) -> V3 {
    let mut n = V3::default();
    for c in contours {
        for (a, b) in c.iter().zip(c.iter().cycle().skip(1)).take(c.len()) {
            let x = a.cross(*b);
            n = V3(n.0 + x.0, n.1 + x.1, n.2 + x.2);
        }
    }
    n
}

fn check_bounds(b: Rect) -> Result<(), String> {
    if [b.x0, b.y0, b.x1, b.y1].iter().all(|v| v.is_finite() && v.abs() <= 1e8) {
        Ok(())
    } else {
        Err("3D effects require finite coordinates within 100 million points".into())
    }
}

/// The art's outward direction at edge `i` of `ring` (before rotation), tilted towards the face
/// it borders by `z` (+1 front, -1 back, 0 a side wall).
fn outward(ring: &[Point], i: usize, z: f64) -> V3 {
    let n = ring.len();
    let (lx, ly) = left(ring[i], ring[(i + 1) % n]);
    V3(-lx, -ly, z)
}

/// The cross-section of a side: each ring of the shape at depth fraction `t` (0 front, 1 back),
/// turned by the twist and scaled by the taper about the centre of `bounds`.
fn section(ring: &[Point], bounds: Rect, o: &Extrude, t: f64) -> Vec<Point> {
    let c = bounds.center();
    let (s, cos) = (o.twist * t).to_radians().sin_cos();
    let k = 1.0 + (o.taper - 1.0) * t;
    ring.iter()
        .map(|p| {
            let (x, y) = ((p.x - c.x) * k, (p.y - c.y) * k);
            Point::new(c.x + x * cos - y * s, c.y + x * s + y * cos)
        })
        .collect()
}

/// One step of an extrusion's outline: how far it is inset (negative: outset) and its depth
/// fraction (0 front, 1 back), with the face it borders (+1 front bevel, -1 back bevel, 0 wall).
#[derive(Clone, Copy)]
struct Step {
    inset: f64,
    t: f64,
    face: f64,
}

/// The outline's steps from the front cap's edge to the back cap's edge.
fn steps(o: &Extrude) -> Vec<Step> {
    let d = o.depth.max(1e-6);
    let mut out = vec![];
    let bevel = o.bevel.filter(|b| b.width > 0.0 || b.height > 0.0);
    // A bevel can't be deeper than the extrusion (half of it on each side with both sides).
    let sides = if bevel.is_some_and(|b| b.both_sides) { 2.0 } else { 1.0 };
    let height = bevel.map_or(0.0, |b| b.height.min(o.depth / sides));
    // The bevel's cross-section across the width: (share across, share down).
    let profile = |b: &Bevel| -> Vec<(f64, f64)> {
        let mut pts = vec![(0.0, 0.0)];
        let n = b.repeat.max(1) as f64;
        let run = (1.0 - b.space) / n;
        let gap = if b.repeat > 1 { b.space / (n - 1.0) } else { 0.0 };
        for k in 0..b.repeat.max(1) {
            let start = k as f64 * (run + gap);
            for s in b.shape.samples().into_iter().skip(1) {
                let (x, y) = b.shape.at(s);
                pts.push((start + x * run, (k as f64 + y) / n));
            }
            if b.repeat > 1 && k + 1 < b.repeat {
                pts.push((start + run + gap, (k as f64 + 1.0) / n));
            }
        }
        if let Some(last) = pts.last_mut() {
            last.0 = 1.0;
        }
        pts
    };
    // With the bevel inside, the cap is inset by its width and the wall is the outline; outside,
    // the cap is the outline and the wall is outset by the width.
    let wall = match bevel {
        Some(b) if !b.inside => -b.width,
        _ => 0.0,
    };
    if let Some(b) = bevel {
        for (x, y) in profile(&b) {
            let inset = if b.inside { b.width * (1.0 - x) } else { -b.width * x };
            out.push(Step { inset, t: y * height / d, face: 1.0 });
        }
    } else {
        out.push(Step { inset: 0.0, t: 0.0, face: 1.0 });
    }
    // Twist and taper bend the walls: cut them into slices.
    let slices = if o.twist != 0.0 || o.taper != 1.0 { ((o.twist.abs() / 6.0).ceil() as usize).clamp(8, 64) } else { 1 };
    let back_height = if bevel.is_some_and(|b| b.both_sides) { height } else { 0.0 };
    let (t0, t1) = (height / d, 1.0 - back_height / d);
    for k in 0..=slices {
        let t = t0 + (t1 - t0) * k as f64 / slices as f64;
        out.push(Step { inset: wall, t, face: 0.0 });
    }
    if let Some(b) = bevel.filter(|b| b.both_sides) {
        for (x, y) in profile(&b).into_iter().rev() {
            let inset = if b.inside { b.width * (1.0 - x) } else { -b.width * x };
            out.push(Step { inset, t: 1.0 - y * height / d, face: -1.0 });
        }
    }
    out.dedup_by(|a, b| (a.inset - b.inset).abs() < 1e-9 && (a.t - b.t).abs() < 1e-12);
    out
}

/// Extrude & Bevel: the shape's outline extruded backwards, with optional twist, taper, bevels
/// and caps.
pub fn extrude(path: &PathData, bounds: Rect, options: Extrude, view: View) -> Result<Vec<Facet>, String> {
    check_bounds(bounds)?;
    let o = options.sanitized();
    let rings = rings(path)?;
    // No depth: the flat shape, seen from either side.
    if o.depth <= 0.0 {
        return plane(path, bounds, view);
    }
    let camera = Camera::new(view, bounds, o.depth);
    let mut b = Builder { camera: &camera, closed: o.cap, facets: vec![] };
    let steps = steps(&o);
    let z = |t: f64| o.depth / 2.0 - o.depth * t;
    let ring_at =
        |ring: &[Point], s: &Step| -> Vec<V3> { section(&inset(ring, s.inset), bounds, &o, s.t).iter().map(|p| camera.place(*p, z(s.t))).collect() };
    for ring in &rings {
        let placed: Vec<Vec<V3>> = steps.iter().map(|s| ring_at(ring, s)).collect();
        for (k, pair) in placed.windows(2).enumerate() {
            let ([a, c], Some(s0), Some(s1)) = (pair, steps.get(k), steps.get(k + 1)) else { continue };
            // A band between two front bevel steps borders the front, one reaching a back bevel
            // step the back; the rest are walls.
            let face = if s1.face < 0.0 {
                -1.0
            } else if s0.face > 0.0 && s1.face > 0.0 {
                1.0
            } else {
                0.0
            };
            b.band(a, c, |i| outward(ring, i, face))?;
        }
    }
    if o.cap {
        for (s, n) in [(steps.first(), 1.0), (steps.last(), -1.0)] {
            let Some(s) = s else { continue };
            b.face(rings.iter().map(|r| ring_at(r, s)).collect(), camera.rotate(V3(0.0, 0.0, n)))?;
        }
    }
    Ok(b.finish())
}

/// Rotate: the flat art turned in 3D about the centre of `bounds` (which may be larger than the
/// shape's own, so that every part of one object turns about the same point).
pub fn plane(path: &PathData, bounds: Rect, view: View) -> Result<Vec<Facet>, String> {
    check_bounds(bounds)?;
    let rings = rings(path)?;
    let camera = Camera::new(view, bounds, 0.0);
    let mut b = Builder { camera: &camera, closed: false, facets: vec![] };
    let contours: Vec<Vec<V3>> = rings.iter().map(|r| r.iter().map(|p| camera.place(*p, 0.0)).collect()).collect();
    b.face(contours, camera.rotate(V3(0.0, 0.0, 1.0)))?;
    Ok(b.finish())
}

/// The distance from `p` to the nearest edge of `rings`.
fn edge_distance(rings: &[Vec<Point>], p: Point) -> f64 {
    let mut best = f64::INFINITY;
    for r in rings {
        for (a, b) in r.iter().zip(r.iter().cycle().skip(1)).take(r.len()) {
            let (dx, dy) = (b.x - a.x, b.y - a.y);
            let l2 = dx * dx + dy * dy;
            let t = if l2 < 1e-18 { 0.0 } else { (((p.x - a.x) * dx + (p.y - a.y) * dy) / l2).clamp(0.0, 1.0) };
            best = best.min(p.distance(Point::new(a.x + dx * t, a.y + dy * t)));
        }
    }
    best
}

/// Is `p` filled by `rings` (nonzero winding, as they are oriented)?
fn filled(rings: &[Vec<Point>], p: Point) -> bool {
    rings.iter().filter(|r| contains(r, p)).count() % 2 == 1
}

/// `ring` cut to the rectangle `cell` (Sutherland–Hodgman: the rectangle is convex, the ring
/// needn't be).
fn clip(ring: &[Point], cell: Rect) -> Vec<Point> {
    let mut pts = ring.to_vec();
    let sides: [fn(Point, Rect) -> f64; 4] = [|p, c| p.x - c.x0, |p, c| c.x1 - p.x, |p, c| p.y - c.y0, |p, c| c.y1 - p.y];
    for side in sides {
        if pts.is_empty() {
            break;
        }
        let mut out = vec![];
        for (a, b) in pts.iter().zip(pts.iter().cycle().skip(1)).take(pts.len()) {
            let (da, db) = (side(*a, cell), side(*b, cell));
            if da >= 0.0 {
                out.push(*a);
            }
            if (da >= 0.0) != (db >= 0.0) {
                out.push(a.lerp(*b, da / (da - db)));
            }
        }
        pts = out;
    }
    pts.dedup_by(|a, b| a.distance(*b) < 1e-9);
    if pts.len() >= 3 && signed_area(&pts).abs() > 1e-12 { pts } else { vec![] }
}

/// The shape's surface as a membrane under pressure, held at zero along its outline (and its
/// holes): the solution of ∇²φ = −1 on a grid, with each node's distance to the outline along the
/// grid lines taken into account where it is less than a grid step (Shortley–Weller), so the
/// outline needn't fall on the grid. √φ rises vertically from the outline and rounds smoothly
/// inside, without the creases a distance-to-the-edge height has along a shape's middle.
struct Membrane {
    origin: Point,
    step: f64,
    nx: usize,
    ny: usize,
    /// φ at each node; outside the shape, extrapolated (negative) from the inside, so
    /// interpolation crosses zero near the outline.
    phi: Vec<f64>,
    peak: f64,
}

impl Membrane {
    fn solve(rings: &[Vec<Point>], origin: Point, step: f64, nx: usize, ny: usize) -> Self {
        let node = |i: usize, j: usize| Point::new(origin.x + i as f64 * step, origin.y + j as f64 * step);
        let idx = |i: usize, j: usize| j * nx + i;
        // A node on the outline itself is held at zero with the outside.
        let inside: Vec<bool> = (0..nx * ny)
            .map(|k| {
                let p = node(k % nx, k / nx);
                filled(rings, p) && edge_distance(rings, p) > 1e-6 * step
            })
            .collect();
        // The nearest crossing of the outline from `p` along (dx, dy), within one step.
        let crossing = |p: Point, dx: f64, dy: f64| -> f64 {
            let mut best = step;
            for r in rings {
                for (a, b) in r.iter().zip(r.iter().cycle().skip(1)).take(r.len()) {
                    let t = if dy == 0.0 {
                        if (a.y > p.y) == (b.y > p.y) {
                            continue;
                        }
                        (a.x + (p.y - a.y) / (b.y - a.y) * (b.x - a.x) - p.x) * dx
                    } else {
                        if (a.x > p.x) == (b.x > p.x) {
                            continue;
                        }
                        (a.y + (p.x - a.x) / (b.x - a.x) * (b.y - a.y) - p.y) * dy
                    };
                    if t > 1e-9 * step && t < best {
                        best = t;
                    }
                }
            }
            best
        };
        // Each inside node's four neighbours: (weight, neighbour inside or None at the outline).
        let mut stencil: Vec<[(f64, Option<usize>); 4]> = vec![[(0.0, None); 4]; nx * ny];
        for j in 0..ny {
            for i in 0..nx {
                let k = idx(i, j);
                if !inside[k] {
                    continue;
                }
                let p = node(i, j);
                let dirs = [(1.0, 0.0), (-1.0, 0.0), (0.0, 1.0), (0.0, -1.0)];
                let h: Vec<f64> = dirs.iter().map(|(dx, dy)| crossing(p, *dx, *dy)).collect();
                for (d, (dx, dy)) in dirs.iter().enumerate() {
                    let opp = h[d ^ 1];
                    let w = 2.0 / (h[d] * (h[d] + opp));
                    let (ni, nj) = (i as isize + *dx as isize, j as isize + *dy as isize);
                    let neighbour = (h[d] >= step && ni >= 0 && nj >= 0 && (ni as usize) < nx && (nj as usize) < ny)
                        .then(|| idx(ni as usize, nj as usize))
                        .filter(|n| inside[*n]);
                    stencil[k][d] = (w, neighbour);
                }
            }
        }
        // Successive over-relaxation.
        let mut phi = vec![0.0; nx * ny];
        let tolerance = 1e-9 * step * step * (nx * ny) as f64;
        for _ in 0..2000 {
            let mut change: f64 = 0.0;
            for k in 0..nx * ny {
                if !inside[k] {
                    continue;
                }
                let (mut sum, mut weight) = (1.0, 0.0);
                for (w, n) in stencil[k] {
                    weight += w;
                    sum += w * n.map_or(0.0, |n| phi[n]);
                }
                let next = phi[k] + 1.8 * (sum / weight - phi[k]);
                change = change.max((next - phi[k]).abs());
                phi[k] = next;
            }
            if change < tolerance {
                break;
            }
        }
        let peak = phi.iter().copied().fold(0.0, f64::max);
        // Outside nodes next to the shape: φ continued past the outline at the slope of their
        // inside neighbours, so values interpolated between them reach zero at the outline.
        let mut ext = phi.clone();
        for j in 0..ny {
            for i in 0..nx {
                let k = idx(i, j);
                if inside[k] {
                    continue;
                }
                let d_out = edge_distance(rings, node(i, j));
                let mut slopes = vec![];
                for (di, dj) in [(1isize, 0isize), (-1, 0), (0, 1), (0, -1)] {
                    let (ni, nj) = (i as isize + di, j as isize + dj);
                    if ni < 0 || nj < 0 || ni as usize >= nx || nj as usize >= ny || !inside[idx(ni as usize, nj as usize)] {
                        continue;
                    }
                    let d_in = edge_distance(rings, node(ni as usize, nj as usize)).max(1e-9);
                    slopes.push(phi[idx(ni as usize, nj as usize)] / d_in);
                }
                if !slopes.is_empty() {
                    ext[k] = -d_out * slopes.iter().sum::<f64>() / slopes.len() as f64;
                }
            }
        }
        Self { origin, step, nx, ny, phi: ext, peak }
    }

    /// φ at `p`, interpolated between the four nodes around it.
    fn at(&self, p: Point) -> f64 {
        let fx = ((p.x - self.origin.x) / self.step).clamp(0.0, (self.nx - 1) as f64);
        let fy = ((p.y - self.origin.y) / self.step).clamp(0.0, (self.ny - 1) as f64);
        let (i, j) = ((fx.floor() as usize).min(self.nx.saturating_sub(2)), (fy.floor() as usize).min(self.ny.saturating_sub(2)));
        let (tx, ty) = (fx - i as f64, fy - j as f64);
        let v = |i: usize, j: usize| self.phi.get(j * self.nx + i).copied().unwrap_or(0.0);
        if self.nx < 2 || self.ny < 2 {
            return v(0, 0);
        }
        let top = v(i, j) * (1.0 - tx) + v(i + 1, j) * tx;
        let bottom = v(i, j + 1) * (1.0 - tx) + v(i + 1, j + 1) * tx;
        top * (1.0 - ty) + bottom * ty
    }
}

/// The part of `outline` inside `cell`, as rings (holes running the other way).
fn material(outline: &PathData, cell: Rect) -> Option<Vec<Vec<Point>>> {
    let r = PathData::single(vectorcraft_geom::SubPath::polyline(
        &[Point::new(cell.x0, cell.y0), Point::new(cell.x1, cell.y0), Point::new(cell.x1, cell.y1), Point::new(cell.x0, cell.y1)],
        true,
    ));
    let rule = vectorcraft_geom::FillRule::NonZero;
    let part =
        vectorcraft_pathops::try_boolean(&r, rule, outline, rule, vectorcraft_pathops::BoolOp::Intersect, vectorcraft_pathops::DEFAULT_PRECISION)
            .ok()?;
    rings(&part).ok()
}

/// Inflate: the shape puffed out like a cushion, its edge kept where the art's outline is.
pub fn inflate(path: &PathData, bounds: Rect, options: Inflate, view: View) -> Result<Vec<Facet>, String> {
    check_bounds(bounds)?;
    let o = options.sanitized();
    let rings = rings(path)?;
    if rings.is_empty() {
        return Ok(vec![]);
    }
    let size = bounds.width().max(bounds.height());
    if size <= 0.0 {
        return Ok(vec![]);
    }
    const CELLS: usize = 48;
    let cell = size / CELLS as f64;
    let (nx, ny) = (((bounds.width() / cell).ceil() as usize).max(1), ((bounds.height() / cell).ceil() as usize).max(1));
    let m = Membrane::solve(&rings, bounds.origin(), cell, nx + 1, ny + 1);
    let outline = PathData::new(rings.iter().map(|r| vectorcraft_geom::SubPath::polyline(r, true)).collect());
    // A circle of radius r puffs into a half sphere of radius r at 50 % (its membrane peaks at r²/4).
    let height = 2.0 * m.peak.sqrt() * o.volume / 50.0;
    let lift = |p: Point| {
        if edge_distance(&rings, p) < 1e-9 || m.peak <= 0.0 {
            return 0.0;
        }
        height * (m.at(p).max(0.0) / m.peak).sqrt()
    };
    let camera = Camera::new(view, bounds, o.depth + 2.0 * height);
    let mut b = Builder { camera: &camera, closed: true, facets: vec![] };
    let half = o.depth / 2.0;
    // The grid's cells, those the surface rises steeply across (near the outline) cut into
    // halves or quarters each way, so the steep rim is drawn in facets as fine as its curve.
    // Every corner sits on one fine lattice (a quarter cell), computed the same way wherever it
    // is used, and a cell's border takes its finer neighbour's lattice points too: where cells
    // of different fineness meet, their edges follow the surface through the same points, with
    // no cracks.
    const FINE: usize = 4;
    let unit = cell / FINE as f64;
    let at = |qx: usize, qy: usize| Point::new(bounds.x0 + qx as f64 * unit, bounds.y0 + qy as f64 * unit);
    let base = |i: usize, j: usize| Rect::from_points(at(i * FINE, j * FINE), at((i + 1) * FINE, (j + 1) * FINE));
    // The shape's part of a rectangle, as rings.
    let cut = |r: Rect| -> Vec<Vec<Point>> {
        let pieces: Vec<Vec<Point>> = rings.iter().map(|ring| clip(ring, r)).filter(|c| !c.is_empty()).collect();
        if pieces.iter().map(|p| signed_area(p)).sum::<f64>() <= 1e-9 * r.area() {
            return vec![];
        }
        // Where a hole's piece cancels part of an outline's, that only holds flat: lifted onto
        // the surface the two no longer meet. Such a cell takes the shape's own part of it.
        let mixed = pieces.iter().any(|p| signed_area(p) < 0.0) && pieces.iter().any(|p| signed_area(p) > 0.0);
        if mixed { material(&outline, r).unwrap_or(pieces) } else { pieces }
    };
    // How many times each cell is halved each way (0, 1 or 2), from how steeply it rises.
    let levels: Vec<u32> = (0..nx * ny)
        .map(|k| {
            let pieces = cut(base(k % nx, k / nx));
            let (lo, hi) = pieces.iter().flatten().map(|q| lift(*q)).fold((f64::INFINITY, f64::NEG_INFINITY), |(lo, hi), h| (lo.min(h), hi.max(h)));
            let rise = if pieces.is_empty() { 0.0 } else { (hi - lo) / cell };
            if rise > 6.0 {
                2
            } else if rise > 2.0 {
                1
            } else {
                0
            }
        })
        .collect();
    let level =
        |i: isize, j: isize| -> u32 { if i < 0 || j < 0 || i as usize >= nx || j as usize >= ny { 0 } else { levels[j as usize * nx + i as usize] } };
    for j in 0..ny {
        for i in 0..nx {
            let l = level(i as isize, j as isize);
            let step = FINE >> l;
            // The lattice spacing along each side of this cell: its own or its neighbour's.
            let (ii, jj) = (i as isize, j as isize);
            let side_step = |dl: u32| FINE >> l.max(dl);
            let (left, right, top, bottom) =
                (side_step(level(ii - 1, jj)), side_step(level(ii + 1, jj)), side_step(level(ii, jj - 1)), side_step(level(ii, jj + 1)));
            let (bx0, by0, bx1, by1) = (i * FINE, j * FINE, (i + 1) * FINE, (j + 1) * FINE);
            for sy in (by0..by1).step_by(step) {
                for sx in (bx0..bx1).step_by(step) {
                    let r = Rect::from_points(at(sx, sy), at(sx + step, sy + step));
                    let pieces = cut(r);
                    if pieces.is_empty() {
                        continue;
                    }
                    // Lattice points along the cell's own border, between each edge's ends.
                    let stitch = |a: Point, b: Point| -> Vec<Point> {
                        let tol = 1e-7 * unit;
                        let along = |fixed_x: Option<usize>, fixed_y: Option<usize>, every: usize| -> Vec<Point> {
                            let mut out = vec![];
                            match (fixed_x, fixed_y) {
                                (Some(qx), None) => {
                                    let (lo, hi) = (a.y.min(b.y), a.y.max(b.y));
                                    for qy in (by0..=by1).step_by(every) {
                                        let p = at(qx, qy);
                                        if p.y > lo + tol && p.y < hi - tol {
                                            out.push(p);
                                        }
                                    }
                                }
                                (None, Some(qy)) => {
                                    let (lo, hi) = (a.x.min(b.x), a.x.max(b.x));
                                    for qx in (bx0..=bx1).step_by(every) {
                                        let p = at(qx, qy);
                                        if p.x > lo + tol && p.x < hi - tol {
                                            out.push(p);
                                        }
                                    }
                                }
                                _ => {}
                            }
                            if (a.x > b.x) || (a.y > b.y) {
                                out.reverse();
                            }
                            out
                        };
                        let on = |v: f64, w: f64| (v - w).abs() <= tol;
                        let (x0, x1, y0, y1) = (at(bx0, by0).x, at(bx1, by1).x, at(bx0, by0).y, at(bx1, by1).y);
                        if on(a.x, x0) && on(b.x, x0) {
                            along(Some(bx0), None, left)
                        } else if on(a.x, x1) && on(b.x, x1) {
                            along(Some(bx1), None, right)
                        } else if on(a.y, y0) && on(b.y, y0) {
                            along(None, Some(by0), top)
                        } else if on(a.y, y1) && on(b.y, y1) {
                            along(None, Some(by1), bottom)
                        } else {
                            vec![]
                        }
                    };
                    let pieces: Vec<Vec<Point>> = pieces
                        .iter()
                        .map(|p| {
                            let mut out = Vec::with_capacity(p.len());
                            for (a, b) in p.iter().zip(p.iter().cycle().skip(1)).take(p.len()) {
                                out.push(*a);
                                out.extend(stitch(*a, *b));
                            }
                            out
                        })
                        .collect();
                    let lifts: Vec<Vec<f64>> = pieces.iter().map(|p| p.iter().map(|q| lift(*q)).collect()).collect();
                    // Each piece lifted onto the surface. Its normal is the lifted piece's own average
                    // (Newell's), not the slope at one point: near the edge the surface turns
                    // vertical, and a point sample there misjudges which way a whole piece faces.
                    let lifted = |sign: f64| -> Vec<Vec<V3>> {
                        pieces.iter().zip(&lifts).map(|(p, h)| p.iter().zip(h).map(|(q, h)| camera.place(*q, sign * (half + h))).collect()).collect()
                    };
                    let front = lifted(1.0);
                    let n = newell(&front);
                    b.face(front, n)?;
                    if o.both_sides {
                        // The back winds the same way round as the front, so its outward normal is
                        // reversed.
                        let back = lifted(-1.0);
                        let n = newell(&back);
                        b.face(back, V3(-n.0, -n.1, -n.2))?;
                    }
                }
            }
        }
    }
    if !o.both_sides {
        let back: Vec<Vec<V3>> = rings.iter().map(|r| r.iter().map(|p| camera.place(*p, -half)).collect()).collect();
        b.face(back, camera.rotate(V3(0.0, 0.0, -1.0)))?;
    }
    if o.depth > 0.0 {
        for ring in &rings {
            let front: Vec<V3> = ring.iter().map(|p| camera.place(*p, half)).collect();
            let back: Vec<V3> = ring.iter().map(|p| camera.place(*p, -half)).collect();
            b.band(&front, &back, |i| outward(ring, i, 0.0))?;
        }
    }
    Ok(b.finish())
}

#[cfg(test)]
#[path = "tests_solid.rs"]
mod tests;