diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml index 60feb80e..b561d637 100644 --- a/.github/workflows/rust.yml +++ b/.github/workflows/rust.yml @@ -28,4 +28,4 @@ jobs: - name: Build workspace, std run: cargo build --workspace --verbose --all-targets --all-features - name: Run tests, std - run: cargo test --workspace --verbose --features "std" --exclude 'retrofire-demos*' + run: cargo test --workspace --verbose --features "std" --exclude 'retrofire-*demo*' diff --git a/.gitignore b/.gitignore index b2f54076..a4fdb2fc 100644 --- a/.gitignore +++ b/.gitignore @@ -1,8 +1,8 @@ -/target -/Cargo.lock -/.idea -/.vscode -/tmp +target/ +Cargo.lock +.idea/ +.vscode/ +tmp/ *.iml *.DS_Store diff --git a/Cargo.toml b/Cargo.toml index b6dfea4b..24a552a5 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -19,6 +19,7 @@ license.workspace = true keywords.workspace = true categories.workspace = true repository.workspace = true +documentation.workspace = true [workspace] members = [ @@ -27,7 +28,7 @@ members = [ "geom", "front", "demos", - "demos/wasm" + "demos/wasm", ] resolver = "2" @@ -36,9 +37,10 @@ edition = "2024" version = "0.4.0" authors = ["Johannes 'Sharlin' Dahlström "] license = "MIT OR Apache-2.0" -repository = "https://github.com/jdahlstrom/retrofire" keywords = ["graphics", "gamedev", "demoscene", "retrocomputing", "rendering"] categories = ["graphics", "game-development", "no-std"] +repository = "https://github.com/jdahlstrom/retrofire" +documentation = "https://crates.io/crates/retrofire" [workspace.lints] clippy.manual_range_contains = "allow" @@ -53,17 +55,33 @@ retrofire-core = { version = "0.4.0", path = "core" } retrofire-front = { version = "0.4.0", path = "front" } retrofire-geom = { version = "0.4.0", path = "geom" } +[dev-dependencies] +divan = "0.1.21" + [profile.release] opt-level = 2 -codegen-units = 1 +codegen-units = 4 lto = "thin" debug = 1 [profile.bench] opt-level = 3 codegen-units = 1 -lto = "fat" +lto = "thin" [profile.dev] opt-level = 1 split-debuginfo = "unpacked" + + +[[bench]] +name = "fill" +harness = false + +[[bench]] +name = "clip" +harness = false + +[[bench]] +name = "isect" +harness = false diff --git a/README.md b/README.md index 20aee6c3..2e264210 100644 --- a/README.md +++ b/README.md @@ -54,13 +54,13 @@ for custom allocators is planned in order to make `alloc` optional as well. * Type-tagged affine and linear transforms and projections * Perspective-correct texture mapping * Triangle mesh data structure and a library of shapes +* Cubic Bézier, Hermite, Catmull–Rom, and B-splines +* Simple random number generation and distributions * Simple text rendering with bitmap fonts * Fully customizable rasterization stage * Collecting rendering performance data * Reading and writing pnm image files * Reading and writing Wavefront .obj files -* Cubic Bezier curves and splines -* Simple random number generation and distributions * Minifb, SDL2, and Wasm frontends * Forever emoji-free README and docs * Forever LLM-free code diff --git a/benches/clip.rs b/benches/clip.rs new file mode 100644 index 00000000..c2450ff3 --- /dev/null +++ b/benches/clip.rs @@ -0,0 +1,91 @@ +//! Triangle clipping benchmarks. + +use core::{array, iter::repeat_with}; + +use divan::{Bencher, counter::ItemsCount}; + +use retrofire_core::{ + geom::{Tri, vertex}, + math::rand::{DEFAULT_RNG, DefaultRng, Distrib}, + math::{orthographic, pt3}, + render::clip::{ClipVert, view_frustum}, +}; + +//#[global_allocator] +//static ALLOC: AllocProfiler = AllocProfiler::system(); + +#[divan::bench(args = [1, 10, 100, 1000, 10_000])] +fn clip_mixed(b: Bencher, n: usize) { + let rng = &mut DefaultRng::default(); + let pts = pt3(-10.0, -10.0, -10.0)..pt3(10.0, 10.0, 10.0); + let proj = orthographic(pt3(-1.0, -1.0, -1.0), pt3(1.0, 1.0, 1.0)); + + b.with_inputs(|| { + repeat_with(|| { + let vs = array::from_fn(|_| { + ClipVert::new(vertex(proj.apply(&pts.sample(rng)), ())) + }); + Tri(vs) + }) + .take(n) + .collect::>() + }) + .input_counter(|tris| ItemsCount::of_iter(tris)) + .bench_local_values(|tris| { + let mut out = Vec::new(); + view_frustum::clip(tris.as_slice(), &mut out); + out + }) +} + +#[divan::bench(args = [1, 10, 100, 1000, 10_000])] +fn clip_all_inside(b: Bencher, n: usize) { + let rng = &mut DefaultRng::default(); + let pts = pt3(-1.0, -1.0, -1.0)..pt3(1.0, 1.0, 1.0); + let proj = orthographic(pt3(-1.0, -1.0, -1.0), pt3(1.0, 1.0, 1.0)); + + b.with_inputs(|| { + repeat_with(|| { + let vs = array::from_fn(|_| { + ClipVert::new(vertex(proj.apply(&pts.sample(rng)), ())) + }); + Tri(vs) + }) + .take(n) + .collect::>() + }) + .input_counter(|tris| ItemsCount::of_iter(tris)) + .bench_local_values(|tris| { + let mut out = Vec::new(); + view_frustum::clip(tris.as_slice(), &mut out); + out + }) +} + +#[divan::bench(args = [1, 10, 100, 1000, 10_000])] +fn clip_all_outside(b: Bencher, n: usize) { + let mut rng = DEFAULT_RNG; + let pts = pt3(2.0, -10.0, -10.0)..pt3(10.0, 10.0, 10.0); + let proj = orthographic(pt3(-1.0, -1.0, -1.0), pt3(1.0, 1.0, 1.0)); + + b.with_inputs(|| { + repeat_with(|| { + let vs = ([pts.start; 3]..[pts.end; 3]) + .sample(&mut rng) + .map(|pt| ClipVert::new(vertex(proj.apply(&pt), ()))); + Tri(vs) + }) + .take(n) + .collect::>() + }) + .input_counter(|tris| ItemsCount::of_iter(tris)) + .bench_local_values(|tris| { + let mut out = Vec::with_capacity(tris.len()); + view_frustum::clip(tris.as_slice(), &mut out); + out + }) +} + +fn main() { + divan::main() +} diff --git a/benches/fill.rs b/benches/fill.rs new file mode 100644 index 00000000..a6f38133 --- /dev/null +++ b/benches/fill.rs @@ -0,0 +1,101 @@ +//! Fillrate benchmarks. + +use core::iter::zip; + +use divan::{Bencher, counter::ItemsCount}; + +use retrofire_core::{ + geom::{Tri, vertex}, + math::{Color3, Color3f, color::gray, pt3, rgb}, + render::{ + Texture, raster::ScreenPt, raster::tri_fill, tex::SamplerRepeatPot, uv, + }, + util::{buf::Buf2, pnm::save_ppm}, +}; + +const SIZES: [f32; 5] = [4.0, 16.0, 64.0, 256.0, 1024.0]; + +const VERTS: [ScreenPt; 3] = + [pt3(0.1, 0.1, 0.0), pt3(0.9, 0.3, 0.5), pt3(0.4, 0.9, 1.0)]; + +#[divan::bench(args = SIZES)] +fn flat(b: Bencher, sz: f32) { + let mut buf: Buf2 = Buf2::new((1024, 1024)); + + b.with_inputs(|| VERTS.map(|p| vertex(p * sz, ()))) + .input_counter(move |vs| ItemsCount::new(Tri(*vs).area() as usize)) + .bench_local_values(|vs| { + tri_fill(vs, |sl| { + buf[sl.y][sl.xs].fill(gray(0xCC)); + }); + }); + + save_ppm("benches_fill_flat.ppm", buf).unwrap(); +} +#[divan::bench(args = SIZES)] +fn gouraud(b: Bencher, sz: f32) { + let mut buf: Buf2 = Buf2::new((1024, 1024)); + + b.with_inputs(|| { + [ + vertex(VERTS[0] * sz, rgb(0.9, 0.1, 0.0)), + vertex(VERTS[1] * sz, rgb(0.1, 0.8, 0.1)), + vertex(VERTS[2] * sz, rgb(0.2, 0.3, 1.0)), + ] + }) + .input_counter(move |vs| ItemsCount::new(Tri(*vs).area() as usize)) + .bench_local_values(|vs| { + tri_fill(vs, |sl| { + let y = sl.y; + let xs = sl.xs.clone(); + let span = &mut buf[y][xs]; + + for ((_, col), pix) in zip(sl.vs, span) { + *pix = col; + } + }); + }); + + let buf = Buf2::new_from( + (1024, 1024), + buf.data().into_iter().map(|c| c.to_color3()), + ); + save_ppm("benches_fill_color.ppm", buf).unwrap(); +} + +#[divan::bench(args = SIZES)] +fn texture(b: Bencher, sz: f32) { + let mut buf: Buf2 = Buf2::new((1024, 1024)); + + let tex = Texture::from(Buf2::::new_from( + (2, 2), + [gray(0xFF), gray(0x33), gray(0x33), gray(0xFF)], + )); + let sampler = SamplerRepeatPot::new(&tex); + + b.with_inputs(|| { + [ + vertex(VERTS[0] * sz, uv(0.0, 0.0)), + vertex(VERTS[1] * sz, uv(4.0, 0.0)), + vertex(VERTS[2] * sz, uv(0.0, 4.0)), + ] + }) + .input_counter(move |vs| ItemsCount::new(Tri(*vs).area() as usize)) + .bench_local_values(|vs| { + tri_fill(vs, |sl| { + let y = sl.y; + let xs = sl.xs.clone(); + let span = &mut buf[y][xs]; + + for ((_, uv), pix) in zip(sl.vs, span) { + *pix = sampler.sample(&tex, uv); + } + }); + }); + + save_ppm("benches_fill_tex.ppm", buf).unwrap(); +} + +fn main() { + divan::main() +} diff --git a/benches/isect.rs b/benches/isect.rs new file mode 100644 index 00000000..2775dc42 --- /dev/null +++ b/benches/isect.rs @@ -0,0 +1,151 @@ +//! Intersection testing benchmarks. + +use core::hint::black_box; + +use divan::{Bencher, counter::ItemsCount}; + +use retrofire::core::{ + geom::{Ray, Sphere}, + math::rand::*, + math::{Point3, degs, pt3, spherical}, + render::scene::BBox, +}; +use retrofire::geom::Intersect; + +#[divan::bench] +fn ray_bbox_hit(b: Bencher) { + let mut rng = DefaultRng::default(); + let bbox = BBox::<()>(pt3(-1.0, -1.0, -1.0), pt3(1.0, 1.0, 1.0)); + + b.with_inputs(|| { + let v = 100.0 * UnitSphere.sample(&mut rng); + Ray(v.to_pt(), -v * (0.0..10.0).sample(&mut rng)) + }) + .counter(ItemsCount::new(1usize)) + .bench_local_values(|ray| { + assert!(ray.intersect(&black_box(bbox)).is_some()) + }); +} +#[divan::bench] +fn ray_bbox_hit_2(b: Bencher) { + let mut rng = DefaultRng::default(); + let bbox = BBox::<()>(pt3(-1.0, -1.0, -1.0), pt3(1.0, 1.0, 1.0)); + + let min = spherical(0.0, degs(-180.0), degs(-90.0)); + let max = spherical(10.0, degs(180.0), degs(-45.0)); + b.with_inputs(|| { + let v = (min..max).sample(&mut rng); + Ray(pt3(0.0, 2.0, 0.0), v.to_cart()) + }) + .counter(ItemsCount::new(1usize)) + .bench_local_values(|ray| { + assert!(ray.intersect(&black_box(bbox)).is_some()) + }); +} +#[divan::bench] +fn ray_bbox_inside(b: Bencher) { + let mut rng = DefaultRng::default(); + let bbox = BBox::<()>(pt3(-1.0, -1.0, -1.0), pt3(1.0, 1.0, 1.0)); + + b.with_inputs(|| { + let pt = PointsInUnitBall.sample(&mut rng); + let dir = VectorsInUnitBall.sample(&mut rng); + Ray(pt, dir) + }) + .counter(ItemsCount::new(1usize)) + .bench_local_values(|ray| { + assert!(ray.intersect(&black_box(bbox)).is_some()) + }); +} + +#[divan::bench] +fn ray_bbox_miss(b: Bencher) { + let mut rng = DefaultRng::default(); + let bbox = BBox::<()>(pt3(-1.0, -1.0, -1.0), pt3(1.0, 1.0, 1.0)); + + b.with_inputs(|| { + let v = (spherical(0.0, degs(-180.0), degs(-45.0)) + ..spherical(10.0, degs(180.0), degs(90.0))) + .sample(&mut rng); + + Ray(pt3(0.0, 3.0, 0.0), v.to_cart()) + }) + .counter(ItemsCount::new(1usize)) + .bench_local_values(|ray| { + assert!(ray.intersect(&black_box(bbox)).is_none()) + }); +} + +#[divan::bench] +fn ray_bbox_mixed(b: Bencher) { + let mut rng = DefaultRng::default(); + let (p, q) = (pt3(-1.0, -1.0, -1.0), pt3(1.0, 1.0, 1.0)); + let bbox = BBox::<()>(p, q); + + b.with_inputs(|| { + // Approximately one third of the rays hits the box + let orig = (p..q).sample(&mut rng); + let dir = (p..q).sample(&mut rng); + Ray(2.0 * orig, 100.0 * dir.to_vec()) + }) + .counter(ItemsCount::new(1usize)) + .bench_local_values(|ray| ray.intersect(&black_box(bbox))); +} + +#[divan::bench] +fn ray_sphere_miss(b: Bencher) { + let mut rng = DefaultRng::default(); + + let sphere = Sphere(::origin(), 1.0); + + let min = spherical(0.0, degs(-180.0), degs(-45.0)); + let max = spherical(10.0, degs(180.0), degs(90.0)); + b.with_inputs(|| { + let v = (min..max).sample(&mut rng); + Ray(pt3(0.0, 3.0, 0.0), v.to_cart()) + }) + .counter(ItemsCount::new(1usize)) + .bench_local_values(|ray| { + let ip = ray.intersect(&black_box(sphere)); + assert!(ip.is_none()); + ip + }); +} + +#[divan::bench] +fn ray_sphere_hit(b: Bencher) { + let mut rng = DefaultRng::default(); + + let sphere = Sphere(::origin(), 1.0); + + b.with_inputs(|| { + let v = (spherical(0.0, degs(-180.0), degs(-90.0)) + ..spherical(10.0, degs(180.0), degs(-45.0))) + .sample(&mut rng); + Ray(pt3(0.0, 2.0f32.sqrt(), 0.0), v.to_cart()) + }) + .counter(ItemsCount::new(1usize)) + .bench_local_values(|ray| { + let ip = ray.intersect(&black_box(sphere)); + assert!(ip.is_some()); + ip + }); +} + +#[divan::bench] +fn ray_sphere_mixed(b: Bencher) { + let mut rng = DefaultRng::default(); + + let sphere = Sphere(::origin(), 1.0); + + b.with_inputs(|| { + let v = VectorsInUnitBall.sample(&mut rng); + Ray(pt3(0.0, 2.0, 0.0), v) + }) + .counter(ItemsCount::new(1usize)) + .bench_local_values(|ray| ray.intersect(&black_box(sphere))); +} + +fn main() { + divan::main(); +} diff --git a/core/Cargo.toml b/core/Cargo.toml index 8925a066..e06566cf 100644 --- a/core/Cargo.toml +++ b/core/Cargo.toml @@ -19,6 +19,7 @@ license.workspace = true keywords.workspace = true categories.workspace = true repository.workspace = true +documentation.workspace = true [features] default = ["std"] diff --git a/core/examples/hello_tri.rs b/core/examples/hello_tri.rs index ec17a6be..3e33b97e 100644 --- a/core/examples/hello_tri.rs +++ b/core/examples/hello_tri.rs @@ -1,3 +1,4 @@ +use retrofire_core::render::{Model, render, shader}; use retrofire_core::{prelude::*, util::*}; fn main() { @@ -48,7 +49,7 @@ fn main() { if cfg!(feature = "fp") { assert_eq!(center_pixel, rgba(151, 128, 187, 255)); } else { - assert_eq!(center_pixel, rgba(114, 102, 127, 255)); + assert_eq!(center_pixel, rgba(114, 102, 128, 255)); } #[cfg(feature = "std")] { diff --git a/core/src/geom.rs b/core/src/geom.rs index fb4af3be..3a6aebac 100644 --- a/core/src/geom.rs +++ b/core/src/geom.rs @@ -1,920 +1,9 @@ -//! Basic geometric primitives. +//! TODO doc +//! -use alloc::vec::Vec; -use core::fmt::Debug; - -use crate::math::{ - Affine, Lerp, Linear, Mat4, Parametric, Point, Point2, Point3, Vec2, Vec3, - Vector, space::Real, vec2, vec3, -}; - -use crate::render::Model; - -pub use mesh::Mesh; +pub use mesh::{Builder, Mesh}; +pub use prim::*; pub mod mesh; -/// Vertex with a position and arbitrary other attributes. -#[derive(Copy, Clone, Debug, Eq, PartialEq)] -pub struct Vertex { - pub pos: P, - pub attrib: A, -} - -/// Two-dimensional vertex type. -pub type Vertex2 = Vertex, A>; - -/// Three-dimensional vertex type. -pub type Vertex3 = Vertex, A>; - -/// Triangle, defined by three vertices. -#[derive(Copy, Clone, Debug, Eq, PartialEq)] -#[repr(transparent)] -pub struct Tri(pub [V; 3]); - -/// Plane, defined by the four parameters of the plane equation. -#[derive(Copy, Clone, Debug, Eq, PartialEq)] -#[repr(transparent)] -pub struct Plane(pub(crate) V); - -/// Plane embedded in 3D space, splitting the space into two half-spaces. -pub type Plane3 = Plane>>; - -/// A ray, or a half line, composed of an initial point and a direction vector. -#[derive(Copy, Clone, Debug, Eq, PartialEq)] -pub struct Ray(pub T, pub T::Diff); - -/// A curve composed of a chain of line segments. -/// -/// The polyline is represented as a list of points, or vertices, with each -/// pair of consecutive vertices sharing an edge. -#[derive(Clone, Debug, Eq, PartialEq)] -pub struct Polyline(pub Vec); - -/// A closed curve composed of a chain of line segments. -/// -/// The polygon is represented as a list of points, or vertices, with each pair -/// of consecutive vertices, as well as the first and last vertex, sharing an edge. -#[derive(Clone, Debug, Eq, PartialEq)] -pub struct Polygon(pub Vec); - -/// A line segment between two vertices. -#[derive(Copy, Clone, Debug, Eq, PartialEq)] -pub struct Edge(pub T, pub T); - -/// A surface normal in 3D. -// TODO Use distinct type rather than alias -pub type Normal3 = Vec3; -/// A surface normal in 2D. -pub type Normal2 = Vec2; - -/// Polygon winding order. -/// -/// The triangle *ABC* below has clockwise winding, while -/// the triangle *DEF* has counter-clockwise winding. -/// -/// ```text -/// B F -/// / \ / \ -/// / \ / \ -/// / \ / \ -/// A-------C D-------E -/// Cw Ccw -/// ``` -#[derive(Copy, Clone, Debug, Default, Eq, PartialEq)] -pub enum Winding { - /// Clockwise winding. - Cw, - /// Counter-clockwise winding. - #[default] - Ccw, -} - -/// Creates a `Vertex` with the give position and attribute values. -#[inline] -pub const fn vertex(pos: P, attrib: A) -> Vertex { - Vertex { pos, attrib } -} - -/// Creates a `Tri` with the given vertices. -#[inline] -pub const fn tri(a: V, b: V, c: V) -> Tri { - Tri([a, b, c]) -} - -// -// Inherent impls -// - -impl Tri { - /// Given a triangle ABC, returns the edges [AB, BC, CA]. - /// - /// # Examples - /// ``` - /// use retrofire_core::geom::{Tri, Edge}; - /// use retrofire_core::math::{Point2, pt2}; - /// - /// let pts: [Point2; _] = [pt2(-1.0, 0.0), pt2(2.0, 0.0), pt2(1.0, 2.0)]; - /// let tri = Tri(pts); - /// - /// let [e0, e1, e2] = tri.edges(); - /// assert_eq!(e0, Edge(&pts[0], &pts[1])); - /// assert_eq!(e1, Edge(&pts[1], &pts[2])); - /// assert_eq!(e2, Edge(&pts[2], &pts[0])); - /// - /// ``` - #[inline] - pub fn edges(&self) -> [Edge<&V>; 3] { - let [a, b, c] = &self.0; - [Edge(a, b), Edge(b, c), Edge(c, a)] - } -} - -impl Tri> { - /// Given a triangle ABC, returns the vectors [AB, AC]. - #[inline] - pub fn tangents(&self) -> [P::Diff; 2] { - let [a, b, c] = &self.0; - [b.pos.sub(&a.pos), c.pos.sub(&a.pos)] - } -} - -impl Tri> { - /// Returns the winding order of `self`. - /// - /// # Examples - /// ``` - /// use retrofire_core::geom::{Tri, vertex, Winding}; - /// use retrofire_core::math::pt2; - /// - /// let mut tri = Tri([ - /// vertex(pt2::<_, ()>(0.0, 0.0), ()), - /// vertex(pt2(0.0, 3.0), ()), - /// vertex(pt2(4.0, 0.0), ()), - /// ]); - /// assert_eq!(tri.winding(), Winding::Cw); - /// - /// tri.0.swap(1, 2); - /// assert_eq!(tri.winding(), Winding::Ccw); - /// ``` - pub fn winding(&self) -> Winding { - let [t, u] = self.tangents(); - if t.perp_dot(u) < 0.0 { - Winding::Cw - } else { - Winding::Ccw - } - } - - /// Returns the signed area of `self`. - /// - /// The area is positive *iff* `self` is wound counter-clockwise. - /// - /// # Examples - /// ``` - /// use retrofire_core::geom::{Tri, vertex}; - /// use retrofire_core::math::pt2; - /// - /// let tri = Tri([ - /// vertex(pt2::<_, ()>(0.0, 0.0), ()), - /// vertex(pt2(0.0, 3.0), ()), - /// vertex(pt2(4.0, 0.0), ()), - /// ]); - /// assert_eq!(tri.signed_area(), -6.0); - /// ``` - pub fn signed_area(&self) -> f32 { - let [t, u] = self.tangents(); - t.perp_dot(u) / 2.0 - } - - /// Returns the (positive) area of `self`. - /// - /// # Examples - /// ``` - /// use retrofire_core::geom::{vertex, Tri}; - /// use retrofire_core::math::pt2; - /// - /// let tri = Tri([ - /// vertex(pt2::<_, ()>(0.0, 0.0), ()), - /// vertex(pt2(0.0, 3.0), ()), - /// vertex(pt2(4.0, 0.0), ()), - /// ]); - /// assert_eq!(tri.area(), 6.0); - /// ``` - pub fn area(&self) -> f32 { - self.signed_area().abs() - } -} - -impl Tri> { - /// Returns the normal vector of `self`. - /// - /// The result is normalized to unit length. - /// - /// # Examples - /// ``` - /// use core::f32::consts::SQRT_2; - /// use retrofire_core::geom::{Tri, vertex}; - /// use retrofire_core::math::{pt3, vec3}; - /// - /// // Triangle lying in a 45° angle - /// let tri = Tri([ - /// vertex(pt3::<_, ()>(0.0, 0.0, 0.0), ()), - /// vertex(pt3(0.0, 3.0, 3.0), ()), - /// vertex(pt3(4.0, 0.0,0.0), ()), - /// ]); - /// assert_eq!(tri.normal(), vec3(0.0, SQRT_2 / 2.0, -SQRT_2 / 2.0)); - /// ``` - pub fn normal(&self) -> Normal3 { - let [t, u] = self.tangents(); - // TODO normal with basis - t.cross(&u).normalize().to() - } - - /// Returns the plane that `self` lies on. - /// - /// # Examples - /// ``` - /// use retrofire_core::geom::{Tri, Plane3, vertex}; - /// use retrofire_core::math::{pt3, Vec3}; - /// - /// let tri = Tri([ - /// vertex(pt3::(0.0, 0.0, 2.0), ()), - /// vertex(pt3(1.0, 0.0, 2.0), ()), - /// vertex(pt3(0.0, 1.0, 2.0), ()) - /// ]); - /// assert_eq!(tri.plane().normal(), Vec3::Z); - /// assert_eq!(tri.plane().offset(), 2.0); - /// ``` - pub fn plane(&self) -> Plane3 { - let [a, b, c] = &self.0; - let [p, q, r] = [a.pos, b.pos, c.pos]; - Plane::from_points(p, q, r) - } - - /// Returns the winding order of `self`, as projected to the XY plane. - // TODO is this 3D version meaningful/useful enough? - pub fn winding(&self) -> Winding { - // TODO better way to xyz->xy... - let [u, v] = self.tangents(); - let ([ux, uy, _], [vx, vy, _]) = (u.0, v.0); - let z = vec2::<_, ()>(ux, uy).perp_dot(vec2(vx, vy)); - if z < 0.0 { Winding::Cw } else { Winding::Ccw } - } - - /// Returns the area of `self`. - /// - /// # Examples - /// ``` - /// use retrofire_core::geom::{tri, vertex}; - /// use retrofire_core::math::pt3; - /// - /// let tri = tri( - /// vertex(pt3::<_, ()>(0.0, 0.0, 0.0), ()), - /// vertex(pt3(4.0, 0.0, 0.0), ()), - /// vertex(pt3(0.0, 3.0, 0.0), ()), - /// ); - /// assert_eq!(tri.area(), 6.0); - /// ``` - #[cfg(feature = "fp")] - pub fn area(&self) -> f32 { - let [t, u] = self.tangents(); - t.cross(&u).len() / 2.0 - } -} - -impl Plane3 { - /// The x = 0 coordinate plane. - pub const YZ: Self = Self::new(1.0, 0.0, 0.0, 0.0); - - /// The y = 0 coordinate plane. - pub const XZ: Self = Self::new(0.0, 1.0, 0.0, 0.0); - - /// The z = 0 coordinate plane. - pub const XY: Self = Self::new(0.0, 0.0, 1.0, 0.0); - - /// Creates a new plane with the given coefficients. - /// - // TODO not normalized because const - // The coefficients are normalized to - // - // (a', b', c', d') = (a, b, c, d) / |(a, b, c)|. - /// - /// The returned plane satisfies the plane equation - /// - /// *ax* + *by* + *cz* = *d*, - /// - /// or equivalently - /// - /// *ax* + *by* + *cz* - *d* = 0. - /// - /// Note the sign of the *d* coefficient. - /// - /// The coefficients (a, b, c) make up a vector normal to the plane, - /// and d is proportional to the plane's distance to the origin. - /// If (a, b, c) is a unit vector, then d is exactly the offset of the - /// plane from the origin in the direction of the normal. - /// - /// # Examples - /// ``` - /// use retrofire_core::{geom::Plane3, math::Vec3}; - /// - /// let p = ::new(1.0, 0.0, 0.0, -2.0); - /// assert_eq!(p.normal(), Vec3::X); - /// assert_eq!(p.offset(), -2.0); - /// - /// ``` - #[inline] - pub const fn new(a: f32, b: f32, c: f32, d: f32) -> Self { - Self(Vector::new([a, b, c, -d])) - } - - /// Creates a plane given three points on the plane. - /// - /// # Panics - /// If the points are collinear or nearly so. - /// - /// # Examples - /// ``` - /// use retrofire_core::{geom::Plane3, math::{pt3, vec3}}; - /// - /// let p = ::from_points( - /// pt3(0.0, 0.0, 2.0), - /// pt3(1.0, 0.0, 2.0), - /// pt3(0.0, 1.0, 2.0), - /// ); - /// assert_eq!(p.normal(), vec3(0.0, 0.0, 1.0)); - /// assert_eq!(p.offset(), 2.0); - /// - /// ``` - pub fn from_points(a: Point3, b: Point3, c: Point3) -> Self { - let n = (b - a).cross(&(c - a)).to(); - Self::from_point_and_normal(a, n) - } - - /// Creates a plane given a point on the plane and a normal. - /// - /// `n` does not have to be normalized. - /// - /// # Panics - /// If `n` is non-finite or nearly zero-length. - /// - /// # Examples - /// ``` - /// use retrofire_core::{geom::Plane3, math::{Vec3, pt3, vec3}}; - /// - /// let p = ::from_point_and_normal(pt3(1.0, 2.0, 3.0), Vec3::Z); - /// assert_eq!(p.normal(), Vec3::Z); - /// assert_eq!(p.offset(), 3.0); - /// - /// ``` - pub fn from_point_and_normal(pt: Point3, n: Normal3) -> Self { - let n = n.normalize(); - // For example, if pt = (0, 1, 0) and n = (0, 1, 0), d has to be 1 - // to satisfy the plane equation n_x + n_y + n_z = d - let d = pt.to_vec().dot(&n.to()); - Plane::new(n.x(), n.y(), n.z(), d) - } - - /// Returns the normal vector of `self`. - /// - /// The normal returned is unit length. - /// - /// # Examples - /// ``` - /// use retrofire_core::{geom::Plane3, math::Vec3}; - /// - /// assert_eq!(::XY.normal(), Vec3::Z); - /// assert_eq!(::YZ.normal(), Vec3::X); - #[inline] - pub fn normal(&self) -> Normal3 { - self.abc().normalize().to() - } - - /// Returns the signed distance of `self` from the origin. - /// - /// This distance is negative if the origin is [*outside*][Self::is_inside] - /// the plane and positive if the origin is *inside* the plane. - /// - /// # Examples - /// ``` - /// use retrofire_core::{geom::Plane3, math::{Vec3, pt3}}; - /// - /// assert_eq!(::new(0.0, 1.0, 0.0, 3.0).offset(), 3.0); - /// assert_eq!(::new(0.0, 2.0, 0.0, 6.0).offset(), 3.0); - /// assert_eq!(::new(0.0, -1.0, 0.0, -3.0).offset(), -3.0); - /// ``` - #[inline] - pub fn offset(&self) -> f32 { - // plane dist from origin is origin dist from plane, negated - -self.signed_dist(Point3::origin()) - } - - /// Returns the perpendicular projection of a point on `self`. - /// - /// In other words, returns *P'*, the point on the plane closest to *P*. - /// - /// ```text - /// ^ P - /// / · - /// / · - /// / · · · · · P' - /// / · - /// / · - /// O------------------> - /// ``` - /// - /// # Examples - /// ``` - /// use retrofire_core::geom::Plane3; - /// use retrofire_core::math::{Point3, pt3}; - /// - /// let pt: Point3 = pt3(1.0, 2.0, -3.0); - /// - /// assert_eq!(::XZ.project(pt), pt3(1.0, 0.0, -3.0)); - /// assert_eq!(::XY.project(pt), pt3(1.0, 2.0, 0.0)); - /// - /// assert_eq!(::new(0.0, 0.0, 1.0, 2.0).project(pt), pt3(1.0, 2.0, 2.0)); - /// assert_eq!(::new(0.0, 0.0, 2.0, 2.0).project(pt), pt3(1.0, 2.0, 1.0)); - /// ``` - pub fn project(&self, pt: Point3) -> Point3 { - // t = -(plane dot orig) / (plane dot dir) - // In this case dir is parallel to plane normal - - let dir = self.abc().to(); - - // TODO add to_homog()/to_real() methods - let pt_hom = [pt.x(), pt.y(), pt.z(), 1.0].into(); - - // Use homogeneous pt to get self · pt = ax + by + cz + d - // Could also just add d manually to ax + by + cz - let plane_dot_orig = self.0.dot(&pt_hom); - - // Vector, so w = 0, so dir_hom · dir_hom = dir · dir - let plane_dot_dir = dir.len_sqr(); // = dir · dir - - let t = -plane_dot_orig / plane_dot_dir; - - pt + t * dir - } - - /// Returns the signed distance of a point to `self`. - /// - /// # Examples - /// ``` - /// use retrofire_core::geom::Plane3; - /// use retrofire_core::math::{Point3, pt3, Vec3}; - /// - /// let pt: Point3 = pt3(1.0, 2.0, -3.0); - /// - /// assert_eq!(::XZ.signed_dist(pt), 2.0); - /// assert_eq!(::XY.signed_dist(pt), -3.0); - /// - /// let p = ::new(-1.0, 0.0, 0.0, 2.0); - /// assert_eq!(p.signed_dist(pt), -3.0); - /// ``` - #[inline] - pub fn signed_dist(&self, pt: Point3) -> f32 { - use crate::math::float::*; - let len_sqr = self.abc().len_sqr(); - // TODO use to_homog once committed - let pt = [pt.x(), pt.y(), pt.z(), 1.0].into(); - self.0.dot(&pt) * f32::recip_sqrt(len_sqr) - } - - /// Returns whether a point is in the half-space that the normal of `self` - /// points away from. - /// - /// # Examples - /// ``` - /// use retrofire_core::geom::Plane3; - /// use retrofire_core::math::{Point3, pt3}; - /// - /// let pt: Point3 = pt3(1.0, 2.0, -3.0); - /// - /// assert!(!::XZ.is_inside(pt)); - /// assert!(::XY.is_inside(pt)); - /// ``` - // TODO "plane.is_inside(point)" reads wrong - #[cfg(feature = "fp")] - #[inline] - pub fn is_inside(&self, pt: Point3) -> bool { - self.signed_dist(pt) <= 0.0 - } - - /// Returns an orthonormal affine basis on `self`. - /// - /// The y-axis of the basis is the normal vector; the x- and z-axes are - /// two arbitrary orthogonal unit vectors tangent to the plane. The origin - /// point is the point on the plane closest to the origin. - /// - /// # Examples - /// ``` - /// use retrofire_core::assert_approx_eq; - /// use retrofire_core::geom::Plane3; - /// use retrofire_core::math::{Point3, pt3, vec3, Apply}; - /// - /// let p = ::from_point_and_normal(pt3(0.0,1.0,0.0), vec3(0.0,1.0,1.0)); - /// let m = p.basis::<()>(); - /// - /// assert_approx_eq!(m.apply(&Point3::origin()), pt3(0.0, 0.5, 0.5)); - /// ``` - pub fn basis(&self) -> Mat4 { - let up = self.abc(); - - let right: Vec3 = - if up.x().abs() < up.y().abs() && up.x().abs() < up.z().abs() { - Vec3::X - } else { - Vec3::Z - }; - let fwd = right.cross(&up).normalize(); - let right = up.normalize().cross(&fwd); - - let origin = self.offset() * up; - - Mat4::from_affine(right, up, fwd, origin.to_pt()) - } - - /// Helper that returns the plane normal non-normalized. - fn abc(&self) -> Vec3 { - let [a, b, c, _] = self.0.0; - vec3(a, b, c) - } -} - -impl Polyline { - /// Creates a new polyline from an iterator of vertex points. - pub fn new(verts: impl IntoIterator) -> Self { - Self(verts.into_iter().collect()) - } - - /// Returns an iterator over the line segments of `self`. - /// - /// # Examples - /// ``` - /// use retrofire_core::geom::{Polyline, Edge}; - /// use retrofire_core::math::{pt2, Point2}; - /// - /// let pts: [Point2; _] = [pt2(0.0, 0.0), pt2(1.0, 1.0), pt2(2.0, 1.0)]; - /// - /// let pline = Polyline::new(pts); - /// let mut edges = pline.edges(); - /// - /// assert_eq!(edges.next(), Some(Edge(&pts[0], &pts[1]))); - /// assert_eq!(edges.next(), Some(Edge(&pts[1], &pts[2]))); - /// assert_eq!(edges.next(), None); - /// ``` - pub fn edges(&self) -> impl Iterator> + '_ { - self.0.windows(2).map(|e| Edge(&e[0], &e[1])) - } -} - -impl Polyline { - /// Returns the sum of the lengths of the edges using a custom metric. - /// - /// The function passed can be arbitrary; this method does not assume any - /// actual metric properties such as positivity or triangle inequality. - /// - /// # Examples - /// ``` - /// use retrofire_core::geom::{Polyline, Edge}; - /// use retrofire_core::math::{pt2, Point2}; - /// - /// let pts: [Point2; _] = [pt2(0.0, 0.0), pt2(1.0, 2.0), pt2(2.0, -3.0)]; - /// let pline = Polyline::new(pts); - /// - /// // The taxicab, or Manhattan, distance. - /// fn taxicab(a: &Point2, b: &Point2) -> f32 { - /// let d = *b - *a; - /// d.x().abs() + d.y().abs() - /// } - /// - /// assert_eq!(pline.len_by(taxicab), 9.0); - /// ``` - pub fn len_by(&self, mut m: impl FnMut(&T, &T) -> f32) -> f32 { - self.edges().map(|Edge(a, b)| m(a, b)).sum() - } -} - -impl Polyline>> { - /// Returns the sum of the lengths of the edges of `self`. - /// - /// # Examples - /// ``` - /// use retrofire_core::geom::{Polyline, Edge}; - /// use retrofire_core::math::{pt2, Point2}; - /// - /// let pts: [Point2; _] = [pt2(0.0, 0.0), pt2(1.0, 0.0), pt2(1.0, -3.0)]; - /// let pline = Polyline::new(pts); - /// - /// assert_eq!(pline.len(), 4.0); - /// ``` - #[cfg(feature = "fp")] - pub fn len(&self) -> f32 { - self.len_by(Point::distance) - } -} - -impl Polygon { - /// Creates a new polygon from an iterator of vertex points. - pub fn new(verts: impl IntoIterator) -> Self { - Self(verts.into_iter().collect()) - } - - /// Returns an iterator over the edges of `self`. - /// - /// Given a polygon ABC...XYZ, returns the edges AB, BC, ..., XY, YZ, ZA. - /// If `self` has zero or one vertices, returns an empty iterator. - /// - /// # Examples - /// ``` - /// use retrofire_core::geom::{Polygon, Edge}; - /// use retrofire_core::math::{Point2, pt2}; - /// - /// let pts: [Point2; _] = [pt2(0.0, 0.0), pt2(1.0, 1.0), pt2(2.0, 1.0)]; - /// - /// let poly = Polygon::new(pts); - /// let mut edges = poly.edges(); - /// - /// assert_eq!(edges.next(), Some(Edge(&pts[0], &pts[1]))); - /// assert_eq!(edges.next(), Some(Edge(&pts[1], &pts[2]))); - /// assert_eq!(edges.next(), Some(Edge(&pts[2], &pts[0]))); - /// assert_eq!(edges.next(), None); - /// ``` - pub fn edges(&self) -> impl Iterator> + '_ { - let last_first = if let [f, .., l] = &self.0[..] { - Some(Edge(l, f)) - } else { - None - }; - self.0 - .windows(2) - .map(|e| Edge(&e[0], &e[1])) - .chain(last_first) - } -} - -// -// Local trait impls -// - -impl Parametric for Ray -where - T: Affine>, -{ - fn eval(&self, t: f32) -> T { - self.0.add(&self.1.mul(t)) - } -} - -impl Parametric for Polyline { - /// Returns the point on `self` at *t*. - /// - /// If the number of vertices in `self` is *n* > 1, the vertex at index - /// *k* < *n* corresponds to `t` = *k* / (*n* - 1). Intermediate values - /// of *t* are linearly interpolated between the two closest vertices. - /// Values *t* < 0 and *t* > 1 are clamped to 0 and 1 respectively. - /// A polyline with a single vertex returns the value of that vertex - /// for any value of *t*. - /// - /// # Panics - /// If `self` has no vertices. - /// - /// # Examples - /// ``` - /// use retrofire_core::geom::{Polyline, Edge}; - /// use retrofire_core::math::{pt2, Point2, Parametric}; - /// - /// let pl = Polyline::( - /// vec![pt2(0.0, 0.0), pt2(1.0, 2.0), pt2(2.0, 1.0)]); - /// - /// assert_eq!(pl.eval(0.0), pl.0[0]); - /// assert_eq!(pl.eval(0.5), pl.0[1]); - /// assert_eq!(pl.eval(1.0), pl.0[2]); - /// - /// // Values not corresponding to a vertex are interpolated: - /// assert_eq!(pl.eval(0.25), pt2(0.5, 1.0)); - /// assert_eq!(pl.eval(0.75), pt2(1.5, 1.5)); - /// - /// // Values of t outside 0.0..=1.0 are clamped: - /// assert_eq!(pl.eval(-1.23), pl.eval(0.0)); - /// assert_eq!(pl.eval(7.68), pl.eval(1.0)); - /// ``` - fn eval(&self, t: f32) -> T { - let pts = &self.0; - assert!(!pts.is_empty(), "cannot eval an empty polyline"); - - let max = pts.len() - 1; - let i = t.clamp(0.0, 1.0) * max as f32; - let t_rem = i % 1.0; - let i = i as usize; - - if i == max { - pts[i].clone() - } else { - pts[i].lerp(&pts[i + 1], t_rem) - } - } -} - -impl Lerp for Vertex { - fn lerp(&self, other: &Self, t: f32) -> Self { - vertex( - self.pos.lerp(&other.pos, t), - // TODO Normals shouldn't be lerped - self.attrib.lerp(&other.attrib, t), - ) - } -} - -#[cfg(test)] -mod tests { - use crate::assert_approx_eq; - use crate::math::*; - use alloc::vec; - - use super::*; - - type Pt = Point<[f32; N], Real>; - - fn tri( - a: Pt, - b: Pt, - c: Pt, - ) -> Tri, ()>> { - Tri([a, b, c].map(|p| vertex(p, ()))) - } - - #[test] - fn triangle_winding_2_cw() { - let tri = tri(pt2(-1.0, 0.0), pt2(0.0, 1.0), pt2(1.0, -1.0)); - assert_eq!(tri.winding(), Winding::Cw); - } - #[test] - fn triangle_winding_2_ccw() { - let tri = tri(pt2(-2.0, 0.0), pt2(1.0, 0.0), pt2(0.0, 1.0)); - assert_eq!(tri.winding(), Winding::Ccw); - } - #[test] - fn triangle_winding_3_cw() { - let tri = - tri(pt3(-1.0, 0.0, 0.0), pt3(0.0, 1.0, 1.0), pt3(1.0, -1.0, 0.0)); - assert_eq!(tri.winding(), Winding::Cw); - } - #[test] - fn triangle_winding_3_ccw() { - let tri = - tri(pt3(-1.0, 0.0, 0.0), pt3(1.0, 0.0, 0.0), pt3(0.0, 1.0, -1.0)); - assert_eq!(tri.winding(), Winding::Ccw); - } - - #[test] - fn triangle_area_2() { - let tri = tri(pt2(-1.0, 0.0), pt2(2.0, 0.0), pt2(2.0, 1.0)); - assert_eq!(tri.area(), 1.5); - } - #[cfg(feature = "fp")] - #[test] - fn triangle_area_3() { - // base = 3, height = 2 - let tri = tri( - pt3(-1.0, 0.0, -1.0), - pt3(2.0, 0.0, -1.0), - pt3(0.0, 0.0, 1.0), - ); - assert_approx_eq!(tri.area(), 3.0); - } - - #[test] - fn triangle_plane() { - let tri = tri( - pt3(-1.0, -2.0, -1.0), - pt3(2.0, -2.0, -1.0), - pt3(0.0, -2.0, 1.0), - ); - assert_approx_eq!(tri.plane().0, Plane3::new(0.0, -1.0, 0.0, 2.0).0); - } - - #[test] - fn plane_from_points() { - let p = ::from_points( - pt3(1.0, 0.0, 0.0), - pt3(0.0, 1.0, 0.0), - pt3(0.0, 0.0, 1.0), - ); - - assert_approx_eq!(p.normal(), vec3(1.0, 1.0, 1.0).normalize()); - assert_approx_eq!(p.offset(), f32::sqrt(1.0 / 3.0)); - } - #[test] - #[should_panic] - fn plane_from_collinear_points_panics() { - ::from_points( - pt3(1.0, 2.0, 3.0), - pt3(-2.0, -4.0, -6.0), - pt3(0.5, 1.0, 1.5), - ); - } - #[test] - #[should_panic] - fn plane_from_zero_normal_panics() { - ::from_point_and_normal( - pt3(1.0, 2.0, 3.0), - vec3(0.0, 0.0, 0.0), - ); - } - #[test] - fn plane_from_point_and_normal() { - let p = ::from_point_and_normal( - pt3(1.0, 2.0, -3.0), - vec3(0.0, 0.0, 12.3), - ); - assert_approx_eq!(p.normal(), vec3(0.0, 0.0, 1.0)); - assert_approx_eq!(p.offset(), -3.0); - } - #[cfg(feature = "fp")] - #[test] - fn plane_is_point_inside_xz() { - let p = ::from_point_and_normal(pt3(1.0, 2.0, 3.0), Vec3::Y); - - // Inside - assert!(p.is_inside(pt3(0.0, 0.0, 0.0))); - // Coincident=inside - assert!(p.is_inside(pt3(0.0, 2.0, 0.0))); - assert!(p.is_inside(pt3(1.0, 2.0, 3.0))); - // Outside - assert!(!p.is_inside(pt3(0.0, 3.0, 0.0))); - assert!(!p.is_inside(pt3(1.0, 3.0, 3.0))); - } - #[cfg(feature = "fp")] - #[test] - fn plane_is_point_inside_neg_xz() { - let p = ::from_point_and_normal(pt3(1.0, 2.0, 3.0), -Vec3::Y); - - // Outside - assert!(!p.is_inside(pt3(0.0, 0.0, 0.0))); - // Coincident=inside - assert!(p.is_inside(pt3(0.0, 2.0, 0.0))); - assert!(p.is_inside(pt3(1.0, 2.0, 3.0))); - // Inside - assert!(p.is_inside(pt3(0.0, 3.0, 0.0))); - assert!(p.is_inside(pt3(1.0, 3.0, 3.0))); - } - #[cfg(feature = "fp")] - #[test] - fn plane_is_point_inside_diagonal() { - let p = ::from_point_and_normal(pt3(0.0, 1.0, 0.0), splat(1.0)); - - // Inside - assert!(p.is_inside(pt3(0.0, 0.0, 0.0))); - assert!(p.is_inside(pt3(-1.0, 1.0, -1.0))); - // Coincident=inside - assert!(p.is_inside(pt3(0.0, 1.0, 0.0))); - // Outside - assert!(!p.is_inside(pt3(0.0, 2.0, 0.0))); - assert!(!p.is_inside(pt3(1.0, 1.0, 1.0))); - assert!(!p.is_inside(pt3(1.0, 0.0, 1.0))); - } - - #[test] - fn plane_project_point() { - let p = ::from_point_and_normal(pt3(0.0, 2.0, 0.0), Vec3::Y); - - // Outside - assert_approx_eq!(p.project(pt3(5.0, 10.0, -3.0)), pt3(5.0, 2.0, -3.0)); - // Coincident - assert_approx_eq!(p.project(pt3(5.0, 2.0, -3.0)), pt3(5.0, 2.0, -3.0)); - // Inside - assert_approx_eq!( - p.project(pt3(5.0, -10.0, -3.0)), - pt3(5.0, 2.0, -3.0) - ); - } - - #[test] - fn polyline_eval_f32() { - let pl = Polyline(vec![0.0, 1.0, -0.5]); - - assert_eq!(pl.eval(-5.0), 0.0); - assert_eq!(pl.eval(0.00), 0.0); - assert_eq!(pl.eval(0.25), 0.5); - assert_eq!(pl.eval(0.50), 1.0); - assert_eq!(pl.eval(0.75), 0.25); - assert_eq!(pl.eval(1.00), -0.5); - assert_eq!(pl.eval(5.00), -0.5); - } - - #[test] - #[should_panic] - fn empty_polyline_eval() { - Polyline::(vec![]).eval(0.5); - } - - #[test] - fn singleton_polyline_eval() { - let pl = Polyline(vec![1.23]); - assert_eq!(pl.eval(0.0), 1.23); - assert_eq!(pl.eval(1.0), 1.23); - } -} +mod prim; diff --git a/core/src/geom/mesh.rs b/core/src/geom/mesh.rs index c6d91d98..068e8faf 100644 --- a/core/src/geom/mesh.rs +++ b/core/src/geom/mesh.rs @@ -88,25 +88,22 @@ impl Mesh { pub fn faces(&self) -> impl Iterator>> { self.faces .iter() - .map(|Tri(vs)| Tri(vs.map(|i| &self.verts[i]))) + .map(|tri| tri.map(|i| &self.verts[i])) } /// Returns a mesh with the faces and vertices of both `self` and `other`. pub fn merge(mut self, Self { faces, verts }: Self) -> Self { let n = self.verts.len(); self.verts.extend(verts); - self.faces.extend( - faces - .into_iter() - .map(|Tri([i, j, k])| Tri([n + i, n + j, n + k])), - ); + self.faces + .extend(faces.into_iter().map(|tri| tri.map(|i| i + n))); self } } #[inline(never)] fn assert_indices_in_bounds(faces: &[Tri], len: usize) { - for (Tri(vs), i) in zip(faces.iter(), 0..) { + for (Tri(vs), i) in zip(faces, 0..) { assert!( vs.iter().all(|&j| j < len), "vertex index out of bounds at faces[{i}]: {vs:?}" @@ -209,10 +206,10 @@ impl Builder { let Mesh { verts, faces } = self.mesh; // Compute weighted face normals... - let face_normals = faces.iter().map(|Tri(vs)| { + let face_normals = faces.iter().map(|tri| { // TODO If n-gonal faces are supported some day, the cross // product is not proportional to area anymore - let [a, b, c] = vs.map(|i| verts[i].pos); + let [a, b, c] = tri.map(|i| verts[i].pos).0; (b - a).cross(&(c - a)).to() }); // ...initialize vertex normals to zero... diff --git a/core/src/geom/prim.rs b/core/src/geom/prim.rs new file mode 100644 index 00000000..5cd48129 --- /dev/null +++ b/core/src/geom/prim.rs @@ -0,0 +1,982 @@ +//! Basic geometric primitives. +//! +//! Includes vertices, polygons, planes, rays, and more. + +use alloc::vec::Vec; +use core::fmt::{self, Debug, Formatter}; + +use crate::math::{ + Affine, Lerp, Linear, Mat4, Parametric, Point, Point2, Point3, Vec2, Vec3, + Vector, pt3, space::Real, vec2, vec3, +}; +use crate::render::Model; + +/// Vertex with a position and arbitrary other attributes. +#[derive(Copy, Clone, Debug, Default, Eq, PartialEq)] +pub struct Vertex { + pub pos: P, + pub attrib: A, +} + +/// Two-dimensional vertex type. +pub type Vertex2 = Vertex, A>; + +/// Three-dimensional vertex type. +pub type Vertex3 = Vertex, A>; + +/// Triangle, defined by three vertices. +#[derive(Copy, Clone, Debug, Default, Eq, PartialEq)] +#[repr(transparent)] +pub struct Tri(pub [V; 3]); + +/// Plane, defined by the four parameters of the plane equation. +#[derive(Copy, Clone, Debug, Eq, PartialEq)] +#[repr(transparent)] +pub struct Plane(pub(crate) V); + +/// Plane embedded in 3D space, splitting the space into two half-spaces. +pub type Plane3 = Plane>>; + +/// A ray, or a half line, composed of an initial point and a direction vector. +#[derive(Copy, Clone, Debug, Eq, PartialEq)] +pub struct Ray(pub T, pub T::Diff); + +pub type Ray3 = Ray>; + +/// A curve composed of a chain of line segments. +/// +/// The polyline is represented as a list of points, or vertices, with each +/// pair of consecutive vertices sharing an edge. +#[derive(Clone, Debug, Default, Eq, PartialEq)] +pub struct Polyline(pub Vec); + +/// A closed curve composed of a chain of line segments. +/// +/// The polygon is represented as a list of points, or vertices, with each pair +/// of consecutive vertices, as well as the first and last vertex, sharing an edge. +#[derive(Clone, Debug, Default, Eq, PartialEq)] +pub struct Polygon(pub Vec); + +/// A line segment between two vertices. +#[derive(Copy, Clone, Debug, Default, Eq, PartialEq)] +pub struct Edge(pub T, pub T); + +#[derive(Copy, Clone, PartialEq)] +pub struct Sphere(pub Point3, pub f32); + +/// A surface normal in 3D. +// TODO Use distinct type rather than alias +pub type Normal3 = Vec3; +/// A surface normal in 2D. +pub type Normal2 = Vec2; + +/// Polygon winding order. +/// +/// The triangle *ABC* below has clockwise winding, while +/// the triangle *DEF* has counter-clockwise winding. +/// +/// ```text +/// B F +/// / \ / \ +/// / \ / \ +/// / \ / \ +/// A-------C D-------E +/// Cw Ccw +/// ``` +#[derive(Copy, Clone, Debug, Default, Eq, PartialEq)] +pub enum Winding { + /// Clockwise winding. + Cw, + /// Counter-clockwise winding. + #[default] + Ccw, +} + +/// Creates a [`Vertex`] with the give position and attribute values. +#[inline] +pub const fn vertex(pos: P, attrib: A) -> Vertex { + Vertex { pos, attrib } +} + +/// Creates a [`Tri`] with the given vertices. +#[inline] +pub const fn tri(a: V, b: V, c: V) -> Tri { + Tri([a, b, c]) +} + +// +// Inherent impls +// + +impl Tri { + /// Given a triangle ABC, returns the edges [AB, BC, CA]. + /// + /// # Examples + /// ``` + /// use retrofire_core::geom::{Tri, Edge}; + /// use retrofire_core::math::{Point2, pt2}; + /// + /// let pts: [Point2; _] = [pt2(-1.0, 0.0), pt2(2.0, 0.0), pt2(1.0, 2.0)]; + /// let tri = Tri(pts); + /// + /// let [e0, e1, e2] = tri.edges(); + /// assert_eq!(e0, Edge(&pts[0], &pts[1])); + /// assert_eq!(e1, Edge(&pts[1], &pts[2])); + /// assert_eq!(e2, Edge(&pts[2], &pts[0])); + /// + /// ``` + #[inline] + pub fn edges(&self) -> [Edge<&V>; 3] { + let [a, b, c] = &self.0; + [Edge(a, b), Edge(b, c), Edge(c, a)] + } + + /// Returns `self` with each vertex mapped with a function. + #[inline] + pub fn map(self, mut f: impl FnMut(V) -> U) -> Tri { + let [a, b, c] = self.0; + Tri([f(a), f(b), f(c)]) + } +} + +impl Tri> { + /// Given a triangle ABC, returns the vectors [AB, AC]. + #[inline] + pub fn tangents(&self) -> [P::Diff; 2] { + let [a, b, c] = &self.0; + [b.pos.sub(&a.pos), c.pos.sub(&a.pos)] + } + + /// Returns the geometric center, or "balance point", of `self`. + /// + /// The centroid is simply the average of the three vertex positions. + pub fn centroid(&self) -> P + where + P::Diff: Linear, + { + let [ab, ac] = self.tangents(); + self.0[0].pos.add(&ab.add(&ac).mul(1.0 / 3.0)) + } +} + +impl Tri> { + /// Returns the winding order of `self`. + /// + /// # Examples + /// ``` + /// use retrofire_core::geom::{Tri, vertex, Winding}; + /// use retrofire_core::math::pt2; + /// + /// let mut tri = Tri([ + /// vertex(pt2::<_, ()>(0.0, 0.0), ()), + /// vertex(pt2(0.0, 3.0), ()), + /// vertex(pt2(4.0, 0.0), ()), + /// ]); + /// assert_eq!(tri.winding(), Winding::Cw); + /// + /// tri.0.swap(1, 2); + /// assert_eq!(tri.winding(), Winding::Ccw); + /// ``` + pub fn winding(&self) -> Winding { + let [t, u] = self.tangents(); + if t.perp_dot(u) < 0.0 { + Winding::Cw + } else { + Winding::Ccw + } + } + + /// Returns the signed area of `self`. + /// + /// The area is positive *iff* `self` is wound counter-clockwise. + /// + /// # Examples + /// ``` + /// use retrofire_core::geom::{Tri, vertex}; + /// use retrofire_core::math::pt2; + /// + /// let tri = Tri([ + /// vertex(pt2::<_, ()>(0.0, 0.0), ()), + /// vertex(pt2(0.0, 3.0), ()), + /// vertex(pt2(4.0, 0.0), ()), + /// ]); + /// assert_eq!(tri.signed_area(), -6.0); + /// ``` + pub fn signed_area(&self) -> f32 { + let [t, u] = self.tangents(); + t.perp_dot(u) / 2.0 + } + + /// Returns the (positive) area of `self`. + /// + /// # Examples + /// ``` + /// use retrofire_core::geom::{vertex, Tri}; + /// use retrofire_core::math::pt2; + /// + /// let tri = Tri([ + /// vertex(pt2::<_, ()>(0.0, 0.0), ()), + /// vertex(pt2(0.0, 3.0), ()), + /// vertex(pt2(4.0, 0.0), ()), + /// ]); + /// assert_eq!(tri.area(), 6.0); + /// ``` + pub fn area(&self) -> f32 { + self.signed_area().abs() + } +} + +impl Tri> { + /// Returns the normal vector of `self`. + /// + /// The result is normalized to unit length. If self is degenerate and + /// has no normal, returns a zero vector. + /// + /// # Examples + /// ``` + /// use core::f32::consts::SQRT_2; + /// use retrofire_core::geom::{Tri, vertex}; + /// use retrofire_core::math::{pt3, vec3}; + /// + /// // Triangle lying in a 45° angle + /// let tri = Tri([ + /// vertex(pt3::<_, ()>(0.0, 0.0, 0.0), ()), + /// vertex(pt3(0.0, 3.0, 3.0), ()), + /// vertex(pt3(4.0, 0.0,0.0), ()), + /// ]); + /// assert_eq!(tri.normal(), vec3(0.0, SQRT_2 / 2.0, -SQRT_2 / 2.0)); + /// ``` + pub fn normal(&self) -> Normal3 { + let [t, u] = self.tangents(); + // TODO normal with basis + t.cross(&u).normalize_or_zero().to() + } + + /// Returns the plane that `self` lies on. + /// + /// # Examples + /// ``` + /// use retrofire_core::geom::{Tri, Plane3, vertex}; + /// use retrofire_core::math::{pt3, Vec3}; + /// + /// let tri = Tri([ + /// vertex(pt3::(0.0, 0.0, 2.0), ()), + /// vertex(pt3(1.0, 0.0, 2.0), ()), + /// vertex(pt3(0.0, 1.0, 2.0), ()) + /// ]); + /// assert_eq!(tri.plane().normal(), Vec3::Z); + /// assert_eq!(tri.plane().offset(), 2.0); + /// ``` + pub fn plane(&self) -> Plane3 { + let [a, b, c] = &self.0; + let [p, q, r] = [a.pos, b.pos, c.pos]; + Plane::from_points(p, q, r) + } + + /// Returns the winding order of `self`, as projected to the XY plane. + // TODO is this 3D version meaningful/useful enough? + pub fn winding(&self) -> Winding { + // TODO better way to xyz->xy... + let [u, v] = self.tangents(); + let ([ux, uy, _], [vx, vy, _]) = (u.0, v.0); + let z = vec2::<_, ()>(ux, uy).perp_dot(vec2(vx, vy)); + if z < 0.0 { Winding::Cw } else { Winding::Ccw } + } + + /// Returns the area of `self`. + /// + /// # Examples + /// ``` + /// use retrofire_core::geom::{tri, vertex}; + /// use retrofire_core::math::pt3; + /// + /// let tri = tri( + /// vertex(pt3::<_, ()>(0.0, 0.0, 0.0), ()), + /// vertex(pt3(4.0, 0.0, 0.0), ()), + /// vertex(pt3(0.0, 3.0, 0.0), ()), + /// ); + /// assert_eq!(tri.area(), 6.0); + /// ``` + pub fn area(&self) -> f32 { + let [t, u] = self.tangents(); + t.cross(&u).len() / 2.0 + } +} + +impl Plane3 { + /// The x = 0 coordinate plane. + pub const YZ: Self = Self::new(1.0, 0.0, 0.0, 0.0); + + /// The y = 0 coordinate plane. + pub const XZ: Self = Self::new(0.0, 1.0, 0.0, 0.0); + + /// The z = 0 coordinate plane. + pub const XY: Self = Self::new(0.0, 0.0, 1.0, 0.0); + + /// Creates a new plane with the given coefficients. + /// + // TODO not normalized because const + // The coefficients are normalized to + // + // (a', b', c', d') = (a, b, c, d) / |(a, b, c)|. + /// + /// The returned plane satisfies the plane equation + /// + /// *ax* + *by* + *cz* = *d*, + /// + /// or equivalently + /// + /// *ax* + *by* + *cz* - *d* = 0. + /// + /// Note the sign of the *d* coefficient. + /// + /// The coefficients (a, b, c) make up a vector normal to the plane, + /// and d is proportional to the plane's distance to the origin. + /// If (a, b, c) is a unit vector, then d is exactly the offset of the + /// plane from the origin in the direction of the normal. + /// + /// # Examples + /// ``` + /// use retrofire_core::{geom::Plane3, math::Vec3}; + /// + /// let p = ::new(1.0, 0.0, 0.0, -2.0); + /// assert_eq!(p.normal(), Vec3::X); + /// assert_eq!(p.offset(), -2.0); + /// + /// ``` + #[inline] + pub const fn new(a: f32, b: f32, c: f32, d: f32) -> Self { + assert!(a != 0.0 || b != 0.0 || c != 0.0, "degenerate plane"); + Self(Vector::new([a, b, c, -d])) + } + + /// Creates a plane given three points on the plane. + /// + /// # Panics + /// If the points are collinear or nearly so. + /// + /// # Examples + /// ``` + /// use retrofire_core::{geom::Plane3, math::{pt3, vec3}}; + /// + /// let p = ::from_points( + /// pt3(0.0, 0.0, 2.0), + /// pt3(1.0, 0.0, 2.0), + /// pt3(0.0, 1.0, 2.0), + /// ); + /// assert_eq!(p.normal(), vec3(0.0, 0.0, 1.0)); + /// assert_eq!(p.offset(), 2.0); + /// + /// ``` + pub fn from_points(a: Point3, b: Point3, c: Point3) -> Self { + let n = (b - a).cross(&(c - a)).to(); + Self::from_point_and_normal(a, n) + } + + /// Creates a plane given a point on the plane and a normal. + /// + /// `n` does not have to be normalized. + /// + /// # Panics + /// If `n` is non-finite or nearly zero-length. + /// + /// # Examples + /// ``` + /// use retrofire_core::{geom::Plane3, math::{Vec3, pt3, vec3}}; + /// + /// let p = ::from_point_and_normal(pt3(1.0, 2.0, 3.0), Vec3::Z); + /// assert_eq!(p.normal(), Vec3::Z); + /// assert_eq!(p.offset(), 3.0); + /// + /// ``` + pub fn from_point_and_normal(pt: Point3, n: Normal3) -> Self { + let n = n.normalize(); + let d = pt.to_vec().dot(&n.to()); + Plane::new(n.x(), n.y(), n.z(), d) + } + + /// Returns the normal vector of `self`. + /// + /// The normal returned is unit length. + /// + /// # Examples + /// ``` + /// use retrofire_core::{geom::Plane3, math::Vec3}; + /// + /// assert_eq!(::XY.normal(), Vec3::Z); + /// assert_eq!(::YZ.normal(), Vec3::X); + #[inline] + pub fn normal(&self) -> Normal3 { + self.abc().normalize().to() + } + + /// Returns the signed distance of `self` from the origin. + /// + /// This distance is negative if the origin is [*outside*][Self::is_inside] + /// the plane and positive if the origin is *inside* the plane. + /// + /// # Examples + /// ``` + /// use retrofire_core::{geom::Plane3, math::{Vec3, pt3}}; + /// + /// assert_eq!(::new(0.0, 1.0, 0.0, 3.0).offset(), 3.0); + /// assert_eq!(::new(0.0, 2.0, 0.0, 6.0).offset(), 3.0); + /// assert_eq!(::new(0.0, -1.0, 0.0, -3.0).offset(), -3.0); + /// ``` + #[inline] + pub fn offset(&self) -> f32 { + // plane dist from origin is origin dist from plane, negated + -self.signed_dist(Point3::origin()) + } + + /// Returns the perpendicular projection of a point on `self`. + /// + /// In other words, returns *P'*, the point on the plane closest to *P*. + /// + /// ```text + /// ^ P + /// / · + /// / · + /// / · · · · · P' + /// / · + /// / · + /// O------------------> + /// ``` + /// + /// # Examples + /// ``` + /// use retrofire_core::geom::Plane3; + /// use retrofire_core::math::{Point3, pt3}; + /// + /// let pt: Point3 = pt3(1.0, 2.0, -3.0); + /// + /// assert_eq!(::XZ.project(pt), pt3(1.0, 0.0, -3.0)); + /// assert_eq!(::XY.project(pt), pt3(1.0, 2.0, 0.0)); + /// + /// assert_eq!(::new(0.0, 0.0, 1.0, 2.0).project(pt), pt3(1.0, 2.0, 2.0)); + /// assert_eq!(::new(0.0, 0.0, 2.0, 2.0).project(pt), pt3(1.0, 2.0, 1.0)); + /// ``` + pub fn project(&self, pt: Point3) -> Point3 { + // t = -(plane dot orig) / (plane dot dir) + // In this case dir is parallel to plane normal + + let dir = self.abc(); + + // TODO add to_homog()/to_real() methods + let pt_hom = [pt.x(), pt.y(), pt.z(), 1.0].into(); + + // Use homogeneous pt to get self · pt = ax + by + cz + d + // Could also just add d manually to ax + by + cz + let plane_dot_orig = self.0.dot(&pt_hom); + + // Vector, so w = 0, so dir_hom · dir_hom = dir · dir = |dir|² + let plane_dot_dir = dir.len_sqr(); + + let t = -plane_dot_orig / plane_dot_dir; + + pt + t * dir + } + + /// Returns the signed distance of a point to `self`. + /// + /// # Examples + /// ``` + /// use retrofire_core::geom::Plane3; + /// use retrofire_core::math::{Point3, pt3, Vec3}; + /// + /// let pt: Point3 = pt3(1.0, 2.0, -3.0); + /// + /// assert_eq!(::XZ.signed_dist(pt), 2.0); + /// assert_eq!(::XY.signed_dist(pt), -3.0); + /// + /// let p = ::new(-1.0, 0.0, 0.0, 2.0); + /// assert_eq!(p.signed_dist(pt), -3.0); + /// ``` + #[inline] + pub fn signed_dist(&self, pt: Point3) -> f32 { + self.dot(pt) / self.abc().len() + } + + /// Returns whether a point is in the half-space that the normal of `self` + /// points away from. + /// + /// # Examples + /// ``` + /// use retrofire_core::geom::Plane3; + /// use retrofire_core::math::{Point3, pt3}; + /// + /// let pt: Point3 = pt3(1.0, 2.0, -3.0); + /// + /// assert!(!::XZ.is_inside(pt)); + /// assert!(::XY.is_inside(pt)); + /// ``` + // TODO "plane.is_inside(point)" reads wrong + #[inline] + pub fn is_inside(&self, pt: Point3) -> bool { + self.dot(pt) <= 0.0 + } + + fn dot(&self, pt: Point3) -> f32 { + // TODO add to_homog method + let [x, y, z] = pt.0; + self.0.dot(&[x, y, z, 1.0].into()) + } + + /// Returns an orthonormal affine basis on `self`. + /// + /// The y-axis of the basis is the normal vector; the x- and z-axes are + /// two arbitrary orthogonal unit vectors tangent to the plane. The origin + /// point is the point on the plane closest to the origin. + /// + /// # Examples + /// ``` + /// use retrofire_core::assert_approx_eq; + /// use retrofire_core::geom::Plane3; + /// use retrofire_core::math::{Point3, pt3, vec3}; + /// + /// let p = ::from_point_and_normal(pt3(0.0,1.0,0.0), vec3(0.0,1.0,1.0)); + /// let m = p.basis::<()>(); + /// + /// assert_approx_eq!(m.apply(&Point3::origin()), pt3(0.0, 0.5, 0.5)); + /// ``` + pub fn basis(&self) -> Mat4 { + let up = self.abc(); + + let right: Vec3 = + if up.x().abs() <= up.y().abs() && up.x().abs() <= up.z().abs() { + Vec3::X + } else { + Vec3::Z + }; + let fwd = right.cross(&up).normalize(); + let up = up.normalize(); + let right = up.cross(&fwd); + + let origin = self.offset() * up; + + Mat4::from_affine(right, up, fwd, origin.to_pt()) + } + + /// Helper that returns the a, b, and c coefficients non-normalized. + fn abc(&self) -> Vec3 { + let [a, b, c, _] = self.0.0; + vec3(a, b, c) + } +} + +impl Polyline { + /// Creates a new polyline from an iterator of vertex points. + pub fn new(verts: impl IntoIterator) -> Self { + Self(verts.into_iter().collect()) + } + + /// Returns an iterator over the line segments of `self`. + /// + /// # Examples + /// ``` + /// use retrofire_core::geom::{Polyline, Edge}; + /// use retrofire_core::math::{pt2, Point2}; + /// + /// let pts: [Point2; _] = [pt2(0.0, 0.0), pt2(1.0, 1.0), pt2(2.0, 1.0)]; + /// + /// let pline = Polyline::new(pts); + /// let mut edges = pline.edges(); + /// + /// assert_eq!(edges.next(), Some(Edge(&pts[0], &pts[1]))); + /// assert_eq!(edges.next(), Some(Edge(&pts[1], &pts[2]))); + /// assert_eq!(edges.next(), None); + /// ``` + pub fn edges(&self) -> impl Iterator> + '_ { + self.0.windows(2).map(|e| Edge(&e[0], &e[1])) + } +} + +impl Polyline { + /// Returns the sum of the lengths of the edges using a custom metric. + /// + /// The function passed can be arbitrary; this method does not assume any + /// actual metric properties such as positivity or triangle inequality. + /// + /// # Examples + /// ``` + /// use retrofire_core::geom::{Polyline, Edge}; + /// use retrofire_core::math::{pt2, Point2}; + /// + /// let pts: [Point2; _] = [pt2(0.0, 0.0), pt2(1.0, 2.0), pt2(2.0, -3.0)]; + /// let pline = Polyline::new(pts); + /// + /// // The taxicab, or Manhattan, distance. + /// fn taxicab(a: &Point2, b: &Point2) -> f32 { + /// let d = *b - *a; + /// d.x().abs() + d.y().abs() + /// } + /// + /// assert_eq!(pline.len_by(taxicab), 9.0); + /// ``` + pub fn len_by(&self, mut m: impl FnMut(&T, &T) -> f32) -> f32 { + self.edges().map(|Edge(a, b)| m(a, b)).sum() + } +} + +impl Polyline>> { + /// Returns the sum of the lengths of the edges of `self`. + /// + /// # Examples + /// ``` + /// use retrofire_core::geom::{Polyline, Edge}; + /// use retrofire_core::math::{pt2, Point2}; + /// + /// let pts: [Point2; _] = [pt2(0.0, 0.0), pt2(1.0, 0.0), pt2(1.0, -3.0)]; + /// let pline = Polyline::new(pts); + /// + /// assert_eq!(pline.len(), 4.0); + /// ``` + pub fn len(&self) -> f32 { + self.len_by(Point::distance) + } +} + +impl Polygon { + /// Creates a new polygon from an iterator of vertex points. + pub fn new(verts: impl IntoIterator) -> Self { + Self(verts.into_iter().collect()) + } + + /// Returns an iterator over the edges of `self`. + /// + /// Given a polygon ABC...XYZ, returns the edges AB, BC, ..., XY, YZ, ZA. + /// If `self` has zero or one vertices, returns an empty iterator. + /// + /// # Examples + /// ``` + /// use retrofire_core::geom::{Polygon, Edge}; + /// use retrofire_core::math::{Point2, pt2}; + /// + /// let pts: [Point2; _] = [pt2(0.0, 0.0), pt2(1.0, 1.0), pt2(2.0, 1.0)]; + /// + /// let poly = Polygon::new(pts); + /// let mut edges = poly.edges(); + /// + /// assert_eq!(edges.next(), Some(Edge(&pts[0], &pts[1]))); + /// assert_eq!(edges.next(), Some(Edge(&pts[1], &pts[2]))); + /// assert_eq!(edges.next(), Some(Edge(&pts[2], &pts[0]))); + /// assert_eq!(edges.next(), None); + /// ``` + pub fn edges(&self) -> impl Iterator> + '_ { + let last_first = if let [f, .., l] = &self.0[..] { + Some(Edge(l, f)) + } else { + None + }; + self.0 + .windows(2) + .map(|e| Edge(&e[0], &e[1])) + .chain(last_first) + } +} + +// +// Local trait impls +// + +impl Parametric for Ray +where + T: Affine>, +{ + fn eval(&self, t: f32) -> T { + self.0.add(&self.1.mul(t)) + } +} + +impl Parametric for Polyline { + /// Returns the point on `self` at *t*. + /// + /// If the number of vertices in `self` is *n* > 1, the vertex at index + /// *k* < *n* corresponds to `t` = *k* / (*n* - 1). Intermediate values + /// of *t* are linearly interpolated between the two closest vertices. + /// Values *t* < 0 and *t* > 1 are clamped to 0 and 1 respectively. + /// A polyline with a single vertex returns the value of that vertex + /// for any value of *t*. + /// + /// # Panics + /// If `self` has no vertices. + /// + /// # Examples + /// ``` + /// use retrofire_core::geom::{Polyline, Edge}; + /// use retrofire_core::math::{pt2, Point2, Parametric}; + /// + /// let pl = Polyline::( + /// vec![pt2(0.0, 0.0), pt2(1.0, 2.0), pt2(2.0, 1.0)]); + /// + /// assert_eq!(pl.eval(0.0), pl.0[0]); + /// assert_eq!(pl.eval(0.5), pl.0[1]); + /// assert_eq!(pl.eval(1.0), pl.0[2]); + /// + /// // Values not corresponding to a vertex are interpolated: + /// assert_eq!(pl.eval(0.25), pt2(0.5, 1.0)); + /// assert_eq!(pl.eval(0.75), pt2(1.5, 1.5)); + /// + /// // Values of t outside 0.0..=1.0 are clamped: + /// assert_eq!(pl.eval(-1.23), pl.eval(0.0)); + /// assert_eq!(pl.eval(7.68), pl.eval(1.0)); + /// ``` + fn eval(&self, t: f32) -> T { + let pts = &self.0; + assert!(!pts.is_empty(), "cannot eval an empty polyline"); + + let max = pts.len() - 1; + let i = t.clamp(0.0, 1.0) * max as f32; + let t_rem = i % 1.0; + let i = i as usize; + + if i == max { + pts[i].clone() + } else { + pts[i].lerp(&pts[i + 1], t_rem) + } + } +} + +impl Lerp for Vertex { + fn lerp(&self, other: &Self, t: f32) -> Self { + vertex( + self.pos.lerp(&other.pos, t), + // TODO Normals shouldn't be lerped + self.attrib.lerp(&other.attrib, t), + ) + } +} + +// +// Foreign trait impls +// + +impl Default for Plane3 { + /// Returns the XZ coordinate plane. + fn default() -> Self { + Plane3::XZ + } +} + +impl Debug for Sphere { + fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result { + f.debug_tuple("Sphere") + .field(&self.0) + .field(&self.1) + .finish() + } +} + +impl Default for Sphere { + /// Returns a unit sphere, with the center at the origin and radius 1. + fn default() -> Self { + Self(pt3(0.0, 0.0, 0.0), 1.0) + } +} + +#[cfg(test)] +mod tests { + use alloc::vec; + use core::f32::consts::FRAC_1_SQRT_2; + + use crate::math::*; + use crate::{assert_approx_eq, mat}; + + use super::*; + + type Pt = Point<[f32; N], Real>; + + fn tri( + a: Pt, + b: Pt, + c: Pt, + ) -> Tri, ()>> { + Tri([a, b, c]).map(|p| vertex(p, ())) + } + + #[test] + fn triangle_winding_2_cw() { + let tri = tri(pt2(-1.0, 0.0), pt2(0.0, 1.0), pt2(1.0, -1.0)); + assert_eq!(tri.winding(), Winding::Cw); + } + #[test] + fn triangle_winding_2_ccw() { + let tri = tri(pt2(-2.0, 0.0), pt2(1.0, 0.0), pt2(0.0, 1.0)); + assert_eq!(tri.winding(), Winding::Ccw); + } + #[test] + fn triangle_winding_3_cw() { + let tri = + tri(pt3(-1.0, 0.0, 0.0), pt3(0.0, 1.0, 1.0), pt3(1.0, -1.0, 0.0)); + assert_eq!(tri.winding(), Winding::Cw); + } + #[test] + fn triangle_winding_3_ccw() { + let tri = + tri(pt3(-1.0, 0.0, 0.0), pt3(1.0, 0.0, 0.0), pt3(0.0, 1.0, -1.0)); + assert_eq!(tri.winding(), Winding::Ccw); + } + + #[test] + fn triangle_area_2() { + let tri = tri(pt2(-1.0, 0.0), pt2(2.0, 0.0), pt2(2.0, 1.0)); + assert_eq!(tri.area(), 1.5); + } + #[test] + fn triangle_area_3() { + // base = 3, height = 2 + let tri = tri( + pt3(-1.0, 0.0, -1.0), + pt3(2.0, 0.0, -1.0), + pt3(0.0, 0.0, 1.0), + ); + assert_approx_eq!(tri.area(), 3.0); + } + + #[test] + fn triangle_plane() { + let tri = tri( + pt3(-1.0, -2.0, -1.0), + pt3(2.0, -2.0, -1.0), + pt3(0.0, -2.0, 1.0), + ); + assert_approx_eq!(tri.plane().0, Plane3::new(0.0, -1.0, 0.0, 2.0).0); + } + + #[test] + fn plane_from_points() { + let p = ::from_points( + pt3(1.0, 0.0, 0.0), + pt3(0.0, 1.0, 0.0), + pt3(0.0, 0.0, 1.0), + ); + + assert_approx_eq!(p.normal(), vec3(1.0, 1.0, 1.0).normalize()); + assert_approx_eq!(p.offset(), f32::sqrt(1.0 / 3.0)); + } + #[test] + #[should_panic] + fn plane_from_collinear_points_panics() { + ::from_points( + pt3(1.0, 2.0, 3.0), + pt3(-2.0, -4.0, -6.0), + pt3(0.5, 1.0, 1.5), + ); + } + #[test] + #[should_panic] + fn plane_from_zero_normal_panics() { + ::from_point_and_normal( + pt3(1.0, 2.0, 3.0), + vec3(0.0, 0.0, 0.0), + ); + } + #[test] + fn plane_from_point_and_normal() { + let p = ::from_point_and_normal( + pt3(1.0, 2.0, -3.0), + vec3(0.0, 0.0, 12.3), + ); + assert_approx_eq!(p.normal(), vec3(0.0, 0.0, 1.0)); + assert_approx_eq!(p.offset(), -3.0); + } + #[test] + fn plane_is_point_inside_xz() { + let p = ::from_point_and_normal(pt3(1.0, 2.0, 3.0), Vec3::Y); + + // Inside + assert!(p.is_inside(pt3(0.0, 0.0, 0.0))); + // Coincident=inside + assert!(p.is_inside(pt3(0.0, 2.0, 0.0))); + assert!(p.is_inside(pt3(1.0, 2.0, 3.0))); + // Outside + assert!(!p.is_inside(pt3(0.0, 3.0, 0.0))); + assert!(!p.is_inside(pt3(1.0, 3.0, 3.0))); + } + #[test] + fn plane_is_point_inside_neg_xz() { + let p = ::from_point_and_normal(pt3(1.0, 2.0, 3.0), -Vec3::Y); + + // Outside + assert!(!p.is_inside(pt3(0.0, 0.0, 0.0))); + // Coincident=inside + assert!(p.is_inside(pt3(0.0, 2.0, 0.0))); + assert!(p.is_inside(pt3(1.0, 2.0, 3.0))); + // Inside + assert!(p.is_inside(pt3(0.0, 3.0, 0.0))); + assert!(p.is_inside(pt3(1.0, 3.0, 3.0))); + } + #[test] + fn plane_is_point_inside_diagonal() { + let p = ::from_point_and_normal(pt3(0.0, 1.0, 0.0), splat(1.0)); + + // Inside + assert!(p.is_inside(pt3(0.0, 0.0, 0.0))); + assert!(p.is_inside(pt3(-1.0, 1.0, -1.0))); + // Coincident=inside + assert!(p.is_inside(pt3(0.0, 1.0, 0.0))); + // Outside + assert!(!p.is_inside(pt3(0.0, 2.0, 0.0))); + assert!(!p.is_inside(pt3(1.0, 1.0, 1.0))); + assert!(!p.is_inside(pt3(1.0, 0.0, 1.0))); + } + + #[test] + fn plane_project_point() { + let p = ::from_point_and_normal(pt3(0.0, 2.0, 0.0), Vec3::Y); + + // Outside + assert_approx_eq!(p.project(pt3(5.0, 10.0, -3.0)), pt3(5.0, 2.0, -3.0)); + // Coincident + assert_approx_eq!(p.project(pt3(5.0, 2.0, -3.0)), pt3(5.0, 2.0, -3.0)); + // Inside + assert_approx_eq!( + p.project(pt3(5.0, -10.0, -3.0)), + pt3(5.0, 2.0, -3.0) + ); + } + + #[test] + fn plane_basis() { + let p = ::new(0.0, 2.0, 2.0, 4.0); + + let m = p.basis::(); + + assert_approx_eq!( + m, + mat![ + 1.0, 0.0, 0.0, 0.0; + 0.0, FRAC_1_SQRT_2, -FRAC_1_SQRT_2, 1.0; + 0.0, FRAC_1_SQRT_2, FRAC_1_SQRT_2, 1.0; + 0.0, 0.0, 0.0, 1.0; + ] + ); + } + + #[test] + fn polyline_eval_f32() { + let pl = Polyline(vec![0.0, 1.0, -0.5]); + + assert_eq!(pl.eval(-5.0), 0.0); + assert_eq!(pl.eval(0.00), 0.0); + assert_eq!(pl.eval(0.25), 0.5); + assert_eq!(pl.eval(0.50), 1.0); + assert_eq!(pl.eval(0.75), 0.25); + assert_eq!(pl.eval(1.00), -0.5); + assert_eq!(pl.eval(5.00), -0.5); + } + + #[test] + #[should_panic] + fn empty_polyline_eval() { + Polyline::(vec![]).eval(0.5); + } + + #[test] + fn singleton_polyline_eval() { + let pl = Polyline(vec![1.23]); + assert_eq!(pl.eval(0.0), 1.23); + assert_eq!(pl.eval(1.0), 1.23); + } +} diff --git a/core/src/lib.rs b/core/src/lib.rs index d3cc55fe..11d2a072 100644 --- a/core/src/lib.rs +++ b/core/src/lib.rs @@ -19,9 +19,9 @@ //! # Crate features //! //! * `std`: -//! Makes available items requiring I/O, timekeeping, or any floating-point -//! functions not included in `core`. In particular this means trigonometric -//! and transcendental functions. +//! Enabled by default. Makes available items requiring I/O, timekeeping, or +//! any floating-point functions not included in `core`. In particular, this +//! means trigonometric and transcendental functions. //! //! If this feature is disabled, the crate only depends on `alloc`. //! @@ -33,7 +33,10 @@ //! Provides fast approximate implementations of floating-point functions //! via the [micromath](https://crates.io/crates/micromath) crate. //! -//! All features are disabled by default. +//! If none of the above features is enabled, fallback implementations of +//! a critical subset of floating-point functions are used. However, all APIs +//! whose implementation relies on trigonometric or transcendental functions +//! are disabled. //! //! # Example //! @@ -61,8 +64,8 @@ pub mod prelude { geom::{ Mesh, Normal2, Normal3, Tri, Vertex, Vertex2, Vertex3, tri, vertex, }, - math::*, - render::*, + math::re_exports::*, + render::re_exports::*, util::buf::{AsMutSlice2, AsSlice2, Buf2, MutSlice2, Slice2}, }; } diff --git a/core/src/math.rs b/core/src/math.rs index 0e5c6105..d90a9aae 100644 --- a/core/src/math.rs +++ b/core/src/math.rs @@ -17,28 +17,39 @@ //! to matching vectors. Angles are strongly typed as well, to allow working //! with different angular units without confusion. -pub use { - angle::{ - Angle, PolarVec, SphericalVec, degs, polar, rads, spherical, turns, - }, - approx::ApproxEq, - color::{Color, Color3, Color3f, Color4, Color4f, rgb, rgba}, - mat::{ - Apply, Mat2, Mat3, Mat4, Matrix, ProjMat3, orthographic, perspective, - scale, scale3, translate, translate3, viewport, - }, - param::Parametric, - point::{Point, Point2, Point2u, Point3, pt2, pt3}, - space::{Affine, Linear}, - spline::{BezierSpline, CubicBezier, smootherstep, smoothstep}, - vary::Vary, - vec::{ProjVec3, Vec2, Vec2i, Vec3, Vec3i, Vector, splat, vec2, vec3}, -}; -#[cfg(feature = "fp")] -pub use { - angle::{acos, asin, atan2}, - mat::{orient_y, orient_z, rotate, rotate_x, rotate_y, rotate_z, rotate2}, -}; +use core::fmt::Debug; + +pub(super) mod re_exports { + pub use super::{ + Lerp, + angle::{ + Angle, PolarVec, SphericalVec, degs, polar, rads, spherical, turns, + }, + approx::ApproxEq, + color::{Color, Color3, Color3f, Color4, Color4f, rgb, rgba}, + lerp, + mat::{ + Apply, Mat2, Mat3, Mat4, Matrix, ProjMat3, orthographic, + perspective, scale, scale3, translate, translate3, viewport, + }, + param::Parametric, + point::{Point, Point2, Point2u, Point3, pt2, pt3}, + space::{Affine, Linear}, + spline::{BezierSpline, CubicBezier, smootherstep, smoothstep}, + vary::Vary, + vec::{ProjVec3, Vec2, Vec2i, Vec3, Vec3i, Vector, splat, vec2, vec3}, + }; + #[cfg(feature = "fp")] + pub use super::{ + angle::{acos, asin, atan2}, + mat::{ + orient_y, orient_z, rotate, rotate_pyr, rotate_x, rotate_y, + rotate_z, rotate2, + }, + }; +} + +pub use re_exports::*; /// Implements an operator trait in terms of an op-assign trait. macro_rules! impl_op { @@ -64,6 +75,7 @@ pub mod angle; pub mod approx; pub mod color; pub mod float; +pub mod grad; pub mod mat; pub mod param; pub mod point; @@ -74,7 +86,7 @@ pub mod vary; pub mod vec; /// Trait for linear interpolation between two values. -pub trait Lerp: Sized { +pub trait Lerp: Clone + Debug + Sized { /// Linearly interpolates between `self` and `other`. /// /// if `t` = 0, returns `self`; if `t` = 1, returns `other`. @@ -146,7 +158,7 @@ pub const SQRT_3: f32 = 1.7320508; impl Lerp for T where - T: Affine>, + T: Affine> + Clone + Debug, { /// Linearly interpolates between `self` and `other`. /// diff --git a/core/src/math/angle.rs b/core/src/math/angle.rs index 1fcdbb71..83c30016 100644 --- a/core/src/math/angle.rs +++ b/core/src/math/angle.rs @@ -2,7 +2,7 @@ use core::{ f32::consts::{PI, TAU}, - fmt::{self, Debug, Display}, + fmt::{self, Debug, Display, Formatter}, marker::PhantomData, ops::{Add, Div, Mul, Neg, Rem, Sub}, ops::{AddAssign, DivAssign, MulAssign, SubAssign}, @@ -26,11 +26,11 @@ use crate::math::{Vec2, Vec3, float::f32, vec2, vec3}; pub struct Angle(f32); /// Tag type for a polar coordinate space -#[derive(Copy, Clone, Debug, Default, Eq, PartialEq)] +#[derive(Copy, Clone, Default, Eq, PartialEq)] pub struct Polar(PhantomData); /// Tag type for a spherical coordinate space. -#[derive(Copy, Clone, Debug, Default, Eq, PartialEq)] +#[derive(Copy, Clone, Default, Eq, PartialEq)] pub struct Spherical(PhantomData); /// A polar coordinate vector, with radius and azimuth components. @@ -107,7 +107,7 @@ pub fn acos(x: f32) -> Angle { /// /// assert_eq!(atan2(0.0, 1.0), degs(0.0)); /// assert_eq!(atan2(2.0, 2.0), degs(45.0)); -/// assert_eq!(atan2(3.0, 0.0), degs(90.0)); +/// assert_eq!(atan2(-3.0, 0.0), degs(-90.0)); /// ``` #[cfg(feature = "fp")] pub fn atan2(y: f32, x: f32) -> Angle { @@ -199,6 +199,22 @@ impl Angle { pub fn clamp(self, min: Self, max: Self) -> Self { Self(self.0.clamp(min.0, max.0)) } + + /// Returns `self` "wrapped around" to the range `min..max`. + /// + /// # Examples + /// ``` + /// use retrofire_core::assert_approx_eq; + /// use retrofire_core::math::{degs, turns}; + /// + /// // 400 (mod 360) = 40 + /// assert_approx_eq!(degs(400.0).wrap(turns(0.0), turns(1.0)), degs(40.0)) + /// ``` + #[must_use] + pub fn wrap(self, min: Self, max: Self) -> Self { + use super::float::f32; + Self(min.0 + f32::rem_euclid(self.0 - min.0, max.0 - min.0)) + } } #[cfg(feature = "fp")] @@ -247,21 +263,6 @@ impl Angle { pub fn tan(self) -> f32 { f32::tan(self.0) } - - /// Returns `self` "wrapped around" to the range `min..max`. - /// - /// # Examples - /// ``` - /// use retrofire_core::assert_approx_eq; - /// use retrofire_core::math::{degs, turns}; - /// - /// // 400 (mod 360) = 40 - /// assert_approx_eq!(degs(400.0).wrap(turns(0.0), turns(1.0)), degs(40.0)) - /// ``` - #[must_use] - pub fn wrap(self, min: Self, max: Self) -> Self { - Self(min.0 + f32::rem_euclid(self.0 - min.0, max.0 - min.0)) - } } impl PolarVec { @@ -330,17 +331,34 @@ impl SphericalVec { /// Returns `self` converted to the equivalent Cartesian 3-vector. /// /// # Examples - /// TODO examples + /// ``` + /// use retrofire_core::assert_approx_eq; + /// use retrofire_core::math::{degs, spherical, vec3, SphericalVec}; + /// + /// let mut v = spherical::<()>(1.0, degs(0.0), degs(0.0)); + /// assert_approx_eq!(v.to_cart(), vec3(1.0, 0.0, 0.0)); + /// + /// v = spherical(2.0, degs(90.0), degs(0.0)); + /// assert_approx_eq!(v.to_cart(), vec3(0.0, 0.0, -2.0)); + /// + /// v = spherical(3.0, degs(0.0), degs(90.0)); + /// assert_approx_eq!(v.to_cart(), vec3(0.0, 3.0, 0.0)); + /// ``` #[cfg(feature = "fp")] pub fn to_cart(&self) -> Vec3 { + // First about z by alt, then about y by az: + // + // ( caz 0 saz ) ( calt -salt 0 ) ( r ) + // ( 0 1 0 ) · ( salt calt 0 ) · ( 0 ) + // (-saz 0 caz ) ( 0 0 1 ) ( 0 ) + // + // ( caz 0 saz ) ( calt·r ) ( caz·calt·r ) + // = ( 0 1 0 ) · ( salt·r ) = ( salt·r ) + // (-saz 0 caz ) ( 0 ) (-saz·calt·r ) + let (sin_alt, cos_alt) = self.alt().sin_cos(); let (sin_az, cos_az) = self.az().sin_cos(); - - let x = cos_az * cos_alt; - let z = sin_az * cos_alt; - let y = sin_alt; - - self.r() * vec3(x, y, z) + self.r() * vec3(cos_az * cos_alt, sin_alt, -sin_az * cos_alt) } } @@ -392,12 +410,11 @@ impl Vec2 { impl Vec3 { /// Converts `self` into the equivalent spherical coordinate vector. /// - /// The `r` component of the result equals `self.len()`. - /// - /// The `az` component is the angle between `self` and the xy-plane in the - /// range (-180°, 180°] such that positive `z` maps to positive `az`. - /// - /// The `alt` component is the angle between `self` and the xz-plane in the + /// Returns a vector (r, az, alt) such that: + /// * `r` equals `self.len()` + /// * `az`is the angle between `self` and the xy-plane in the range + /// (-180°, 180°] such that positive `z` maps to *negative* `az`, and + /// * `alt` is the angle between `self` and the xz-plane in the /// range [-90°, 90°] such that positive `y` maps to positive `alt`. /// /// # Examples @@ -406,23 +423,23 @@ impl Vec3 { /// /// // The positive x-axis lies at zero azimuth and altitude /// assert_eq!( - /// vec3(2.0, 0.0, 0.0).to_spherical(), - /// spherical::<()>(2.0, degs(0.0), degs(0.0)) + /// vec3(1.0, 0.0, 0.0).to_spherical(), + /// spherical::<()>(1.0, degs(0.0), degs(0.0)) /// ); /// // The positive y-axis lies at 90° altitude /// assert_eq!( /// vec3(0.0, 2.0, 0.0).to_spherical(), /// spherical::<()>(2.0, degs(0.0), degs(90.0)) /// ); - /// // The positive z axis lies at 90° azimuth + /// // The positive z-axis lies at *-90°* azimuth /// assert_eq!( - /// vec3(0.0, 0.0, 2.0).to_spherical(), - /// spherical::<()>(2.0, degs(90.0), degs(0.0)) + /// vec3(0.0, 0.0, 3.0).to_spherical(), + /// spherical::<()>(3.0, degs(-90.0), degs(0.0)) /// ); /// ``` pub fn to_spherical(&self) -> SphericalVec { let [x, y, z] = self.0; - let az = atan2(z, x); + let az = atan2(-z, x); let alt = atan2(y, f32::sqrt(x * x + z * z)); let r = self.len(); spherical(r, az, alt) @@ -477,6 +494,17 @@ impl ZDiv for Angle {} // Foreign trait impls // +impl Debug for Polar { + fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result { + f.write_str("Pol") + } +} +impl Debug for Spherical { + fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result { + f.write_str("Sph") + } +} + impl Default for SphericalVec { fn default() -> Self { Self::new([1.0, 0.0, 0.0]) @@ -491,7 +519,7 @@ impl Default for PolarVec { impl Display for Angle { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { let (val, unit) = if f.alternate() { - (self.to_rads() / PI, "𝜋 rad") + (self.to_rads() / PI, "·𝜋 rad") } else { (self.to_degs(), "°") }; @@ -503,8 +531,8 @@ impl Display for Angle { impl Debug for Angle { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { f.write_str("Angle(")?; - Display::fmt(self, f)?; - f.write_str(")") + Debug::fmt(&(self.to_rads() / PI), f)?; + f.write_str("·𝜋 rad)") } } @@ -618,6 +646,7 @@ impl From> for SphericalVec { #[allow(unused, nonstandard_style)] mod tests { use core::f32::consts::{PI, TAU}; + use std::eprintln; use crate::{ assert_approx_eq, @@ -626,8 +655,12 @@ mod tests { use super::*; - const vec2: fn(f32, f32) -> Vec2 = math::vec2; - const vec3: fn(f32, f32, f32) -> Vec3 = math::vec3; + const fn vec2(x: f32, y: f32) -> Vec2 { + math::vec2(x, y) + } + const fn vec3(x: f32, y: f32, z: f32) -> Vec3 { + math::vec3(x, y, z) + } #[test] fn rads_to_degs() { @@ -710,7 +743,6 @@ mod tests { assert_approx_eq!(atan2(-1.0, 1.0), degs(-45.0)); } - #[cfg(feature = "fp")] #[test] fn wrapping() { use crate::assert_approx_eq; @@ -788,68 +820,46 @@ mod tests { assert_eq!(vec2(0.0, -4.0).to_polar(), polar(4.0, degs(-90.0))); } - #[cfg(feature = "fp")] - #[test] - fn spherical_to_cartesian() { - let spherical = spherical::<()>; - assert_eq!( - spherical(0.0, degs(0.0), degs(0.0)).to_cart(), - vec3(0.0, 0.0, 0.0) - ); - assert_eq!( - spherical(1.0, degs(0.0), degs(0.0)).to_cart(), - vec3(1.0, 0.0, 0.0) - ); - assert_approx_eq!( - spherical(2.0, degs(60.0), degs(0.0)).to_cart(), - vec3(1.0, 0.0, SQRT_3) - ); - assert_approx_eq!( - spherical(2.0, degs(90.0), degs(0.0)).to_cart(), - vec3(0.0, 0.0, 2.0) - ); - assert_approx_eq!( - spherical(3.0, degs(123.0), degs(90.0)).to_cart(), - vec3(0.0, 3.0, 0.0) - ); + const fn sph(r: f32, az: f32, alt: f32) -> SphericalVec { + spherical(r, degs(az), degs(alt)) } + #[rustfmt::skip] + const CART_SPH: [(Vec3, SphericalVec); 10] = [ + (vec3( 0.0, 0.0, 0.0), sph(0.0, 0.0, 0.0)), + + (vec3( 1.0, 0.0, 0.0), sph(1.0, 0.0, 0.0)), + (vec3( SQRT_3, 0.0, -1.0), sph(2.0, 30.0, 0.0)), + (vec3( 1.0, 0.0, -SQRT_3), sph(2.0, 60.0, 0.0)), + (vec3( 0.0, 0.0, -2.0), sph(2.0, 90.0, 0.0)), + (vec3(-SQRT_3, 0.0, -1.0), sph(2.0, 150.0, 0.0)), + + // Doesn't roundtrip due to imprecision and + // the discontinuity from 180° to -180° :( + (vec3( -3.0, 0.0, 0.0), sph(3.0, 180.0, 0.0)), + + (vec3(SQRT_3, 1.0, 0.0), sph(2.0, 0.0, 30.0)), + (vec3( 1.0, SQRT_3, 0.0), sph(2.0, 0.0, 60.0)), + (vec3( 0.0, 2.0, 0.0), sph(2.0, 0.0, 90.0)), + (vec3( 0.0, -3.0, 0.0), sph(3.0, 0.0, -90.0)), + ]; #[cfg(feature = "fp")] #[test] - fn cartesian_to_spherical_zero_alt() { - assert_approx_eq!( - vec3(0.0, 0.0, 0.0).to_spherical(), - spherical(0.0, degs(0.0), degs(0.0)) - ); - assert_eq!( - vec3(1.0, 0.0, 0.0).to_spherical(), - spherical(1.0, degs(0.0), degs(0.0)) - ); - assert_approx_eq!( - vec3(1.0, SQRT_3, 0.0).to_spherical(), - spherical(2.0, degs(0.0), degs(60.0)) - ); - assert_eq!( - vec3(0.0, 2.0, 0.0).to_spherical(), - spherical(2.0, degs(0.0), degs(90.0)) - ); + fn spherical_to_cartesian() { + for (cart, sp) in CART_SPH { + let actual = sp.to_cart(); + eprintln!("Testing {sp:?} -> {cart:?}"); + assert_approx_eq!(actual, cart); + } } #[cfg(feature = "fp")] #[test] fn cartesian_to_spherical() { - use core::f32::consts::SQRT_2; - assert_approx_eq!( - vec3(SQRT_3, 0.0, 1.0).to_spherical(), - spherical(2.0, degs(30.0), degs(0.0)) - ); - assert_approx_eq!( - vec3(1.0, SQRT_2, 1.0).to_spherical(), - spherical(2.0, degs(45.0), degs(45.0)) - ); - assert_approx_eq!( - vec3(0.0, 0.0, 3.0).to_spherical(), - spherical(3.0, degs(90.0), degs(0.0)) - ); + for (cart, sp) in CART_SPH { + let actual = cart.to_spherical(); + eprintln!("Testing {cart:?} -> {sp:?}"); + assert_approx_eq!(actual, sp); + } } } diff --git a/core/src/math/color.rs b/core/src/math/color.rs index 51fc2c44..eab9acab 100644 --- a/core/src/math/color.rs +++ b/core/src/math/color.rs @@ -24,7 +24,6 @@ use super::{Affine, Linear, Vector, vary::ZDiv}; /// Color components are also called *channels*. /// * `Space`: the color space that `Self` is an element of. #[repr(transparent)] -#[derive(Copy, Clone, Default, Eq, PartialEq)] pub struct Color(pub Repr, PhantomData); /// The (S)RGB (red, green, blue) color space. @@ -123,6 +122,28 @@ impl Color<[Ch; N], Sp> { Color::new(self.0.map(f)) } } +impl Color<[f32; N], Sp> { + /// Returns `self` clamped channel-wise to the given range. + /// + /// In other words, for each channel `self[i]`, the result `r` has + /// `r[i]` equal to `self[i].clamp(min[i], max[i])`. + /// + /// # Examples + /// ``` + /// use retrofire_core::math::color::{Color3f, gray, rgb}; + /// let c: Color3f = rgb(0.0, 0.5, 1.0); + /// + /// let clamped = c.clamp(&gray(0.1), &gray(0.5)); + /// assert_eq!(clamped, rgb(0.1, 0.5, 0.5)); + // TODO f32 and f64 have inherent clamp methods because they're not Ord. + // A generic clamp for Sc: Ord would conflict with this one. There is + // currently no clean way to support both floats and impl Ord types. + // However, ColorX should have its own inherent impls. + #[must_use] + pub fn clamp(&self, min: &Self, max: &Self) -> Self { + array::from_fn(|i| self[i].clamp(min[i], max[i])).into() + } +} impl Color3 { /// Returns `self` as RGBA, with alpha set to 0xFF (fully opaque). @@ -550,7 +571,7 @@ where // Local trait impls // -impl Affine for Color<[u8; DIM], Sp> { +impl Affine for Color<[u8; DIM], Sp> { type Space = Sp; // Color is currently not Linear, so use Vector for now type Diff = Vector<[i32; DIM], Sp>; @@ -569,7 +590,10 @@ impl Affine for Color<[u8; DIM], Sp> { } } -impl Affine for Color<[f32; DIM], Sp> { +impl Affine for Color<[f32; DIM], Sp> +where + Self: Debug, +{ type Space = Sp; type Diff = Self; @@ -585,7 +609,10 @@ impl Affine for Color<[f32; DIM], Sp> { } } -impl Linear for Color<[f32; DIM], Sp> { +impl Linear for Color<[f32; DIM], Sp> +where + Self: Debug, +{ type Scalar = f32; /// Returns the all-zeroes color (black). @@ -604,12 +631,36 @@ impl ZDiv for Color<[Sc; N], Sp> where Sc: ZDiv + Copy { // Foreign trait impls // +// Implemented manually to avoid bound on space + +impl Copy for Color {} + +impl Clone for Color { + fn clone(&self) -> Self { + self.0.clone().into() + } +} + impl Debug for Color { fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result { write!(f, "Color<{:?}>{:?}", Space::default(), self.0) } } +impl Default for Color { + fn default() -> Self { + R::default().into() + } +} + +impl Eq for Color {} + +impl PartialEq for Color { + fn eq(&self, other: &Self) -> bool { + self.0 == other.0 + } +} + impl From for Color { #[inline] fn from(els: R) -> Self { diff --git a/core/src/math/float.rs b/core/src/math/float.rs index 73a56ed7..956dd65b 100644 --- a/core/src/math/float.rs +++ b/core/src/math/float.rs @@ -26,8 +26,9 @@ pub mod libm { pub use super::fallback::rem_euclid; + #[inline] pub fn recip_sqrt(x: f32) -> f32 { - powf(x, -0.5) + 1.0 / sqrt(x) } } @@ -47,7 +48,8 @@ pub mod mm { #[inline] pub fn sqrt(x: f32) -> f32 { let y = mm::sqrt(x); - // One round of Newton's method + // Two rounds of Newton's method + let y = 0.5 * (y + (x / y)); 0.5 * (y + (x / y)) } /// Returns the approximate reciprocal of the square root of `x`. @@ -61,6 +63,7 @@ pub mod mm { pub fn powf(x: f32, y: f32) -> f32 { mm::powf(x, y) } + #[inline] pub fn sin(x: f32) -> f32 { mm::sin(x) @@ -90,27 +93,45 @@ pub mod mm { } mm::atan2(y, x) } + + #[inline] + pub fn exp(x: f32) -> f32 { + mm::exp(x) + } + #[inline] + pub fn log2(x: f32) -> f32 { + mm::log2(x) + } } pub mod fallback { /// Returns the largest integer less than or equal to `x`. #[inline] pub fn floor(x: f32) -> f32 { - (x as i64 - x.is_sign_negative() as i64) as f32 + (x as i64 - (x < 0.0) as i64) as f32 } - // Returns the least non-negative remainder of `x` (mod `m`). + /// Returns the least non-negative remainder of `x` (mod `m`). #[inline] pub fn rem_euclid(x: f32, m: f32) -> f32 { - x % m + (x.is_sign_negative() as u32 as f32) * m + let r = x % m; + r + if r < 0.0 { m.abs() } else { 0.0 } } /// Returns the approximate reciprocal of the square root of `x`. #[inline] pub fn recip_sqrt(x: f32) -> f32 { + if x < 0.0 { + return f32::NAN; + } // https://en.wikipedia.org/wiki/Fast_inverse_square_root let y = f32::from_bits(0x5f37_5a86 - (x.to_bits() >> 1)); - // A round of Newton's method + // Two rounds of Newton's method + let y = y * (1.5 - 0.5 * x * y * y); y * (1.5 - 0.5 * x * y * y) } + #[inline] + pub fn sqrt(x: f32) -> f32 { + 1.0 / recip_sqrt(x) + } } #[cfg(feature = "std")] @@ -141,59 +162,134 @@ pub use fallback as f32; #[cfg(test)] #[allow(unused_imports)] mod tests { + use core::f32::consts::*; + use super::{RecipSqrt, f32, *}; use crate::assert_approx_eq; #[cfg(feature = "libm")] #[test] fn libm_functions() { - use super::libm; - use core::f32::consts::PI; - assert_eq!(libm::cos(PI), -1.0); + assert_eq!(libm::floor(1.5), 1.0); + assert_eq!(libm::floor(0.99), 0.0); + assert_eq!(libm::floor(-0.0), 0.0); + assert_eq!(libm::floor(-1.1), -2.0); + + assert_approx_eq!(libm::rem_euclid(1.6, 0.5), 0.1); + assert_approx_eq!(libm::rem_euclid(-1.6, 0.5), 0.4); + assert_approx_eq!(libm::rem_euclid(1.6, -0.5), 0.1); + assert_approx_eq!(libm::rem_euclid(-1.6, -0.5), 0.4); + assert_eq!(libm::sqrt(9.0), 3.0); + assert_eq!(libm::sqrt(16.0), 4.0); + assert!(libm::sqrt(-1.0).is_nan()); + assert_eq!(libm::recip_sqrt(9.0), 1.0 / 3.0); + assert_eq!(libm::recip_sqrt(0.0), f32::INFINITY); + assert!(libm::recip_sqrt(-1.0).is_nan()); + + assert_eq!(libm::powf(3.0, 2.0), 9.0); + assert_eq!(libm::powf(-3.0, 2.0), 9.0); + assert_eq!(libm::powf(3.0, -2.0), 1.0 / 9.0); + assert_eq!(libm::powf(-3.0, 3.0), -27.0); + + assert_approx_eq!(libm::sin(FRAC_PI_6), 0.5); + assert_eq!(libm::cos(PI), -1.0); + + assert_eq!(libm::exp(1.0), E); + assert_approx_eq!(libm::exp(2.0), E * E); + assert_eq!(libm::log2(8.0), 3.0); + assert!(libm::log2(-1.0).is_nan()); } #[cfg(feature = "mm")] #[test] fn mm_functions() { - use core::f32::consts::*; + assert_eq!(mm::floor(1.5), 1.0); + assert_eq!(mm::floor(0.99), 0.0); + assert_eq!(mm::floor(-0.0), 0.0); + assert_eq!(mm::floor(-1.1), -2.0); - use super::f32; + assert_approx_eq!(mm::rem_euclid(1.6, 0.5), 0.1); + assert_approx_eq!(mm::rem_euclid(-1.6, 0.5), 0.4); + assert_approx_eq!(mm::rem_euclid(1.6, -0.5), 0.1); + assert_approx_eq!(mm::rem_euclid(-1.6, -0.5), 0.4); - assert_approx_eq!(f32::sin(FRAC_PI_6), 0.5); - assert_eq!(f32::cos(PI), -1.0); - assert_eq!(f32::sqrt(16.0), 4.0); - assert_approx_eq!(f32::sqrt(9.0), 3.0, eps = 1e-3); + assert_approx_eq!(mm::sqrt(9.0), 3.0); + assert_eq!(mm::sqrt(16.0), 4.0); + assert!(mm::sqrt(-1.0).is_nan()); + assert_approx_eq!(mm::recip_sqrt(9.0), 1.0 / 3.0); + // mm doesn't check for zero, just gives a big number + assert_approx_eq!(mm::recip_sqrt(0.0), 1.9818e19); + // mm doesn't check for negative, panics due to sub overflow + //assert!(mm::recip_sqrt(-1.0).is_nan()); + + assert_approx_eq!(mm::powf(3.0, 2.0), 9.0); + assert_approx_eq!(mm::powf(-3.0, 2.0), 9.0); + assert_approx_eq!(mm::powf(3.0, -2.0), 1.0 / 9.0); + assert_approx_eq!(mm::powf(-3.0, 3.0), -27.0); + + assert_approx_eq!(mm::sin(FRAC_PI_6), 0.5); + assert_eq!(mm::cos(PI), -1.0); + + assert_eq!(mm::exp(1.0), E); + assert_approx_eq!(mm::exp(2.0), E * E); + assert_approx_eq!(mm::log2(8.0), 3.0); + // mm doesn't check for negative, panics due to sub overflow + //assert!(mm::log2(-1.0).is_nan()); } #[cfg(feature = "std")] #[test] fn std_functions() { - use super::f32; - use core::f32::consts::PI; - assert_eq!(f32::cos(PI), -1.0); + assert_eq!(f32::floor(-0.0), 0.0); + + assert_approx_eq!(f32::rem_euclid(1.6, 0.5), 0.1); + assert_approx_eq!(f32::rem_euclid(-1.6, 0.5), 0.4); + assert_approx_eq!(f32::rem_euclid(1.6, -0.5), 0.1); + assert_approx_eq!(f32::rem_euclid(-1.6, -0.5), 0.4); + assert_eq!(f32::sqrt(9.0), 3.0); + assert!(f32::sqrt(-1.0).is_nan()); + assert_eq!(f32::recip_sqrt(9.0), 1.0 / 3.0); + assert_eq!(f32::recip_sqrt(0.0), f32::INFINITY); + assert!(f32::recip_sqrt(-1.0).is_nan()); + + assert_eq!(f32::cos(PI), -1.0); } #[cfg(not(feature = "fp"))] #[test] fn fallback_functions() { - use super::{RecipSqrt, f32}; + use fallback as fb; + assert_eq!(fb::floor(1.5), 1.0); + assert_eq!(fb::floor(0.99), 0.0); + assert_eq!(fb::floor(-0.0), 0.0); + assert_eq!(fb::floor(-1.1), -2.0); - assert_eq!(f32::floor(1.23), 1.0); - assert_eq!(f32::floor(0.0), 0.0); - assert_eq!(f32::floor(-1.23), -2.0); + assert_approx_eq!(fb::rem_euclid(1.6, 0.5), 0.1); + assert_approx_eq!(fb::rem_euclid(-1.6, 0.5), 0.4); + assert_approx_eq!(fb::rem_euclid(1.6, -0.5), 0.1); + assert_approx_eq!(fb::rem_euclid(-1.6, -0.5), 0.4); - assert_approx_eq!(f32::rem_euclid(1.23, 4.0), 1.23); - assert_approx_eq!(f32::rem_euclid(4.0, 4.0), 0.0); - assert_approx_eq!(f32::rem_euclid(5.67, 4.0), 1.67); - assert_approx_eq!(f32::rem_euclid(-1.23, 4.0), 2.77); - } + assert_approx_eq!(fb::sqrt(9.0), 3.0); + assert_approx_eq!(fb::sqrt(16.0), 4.0); + assert!(fb::sqrt(-1.0).is_nan()); + assert_approx_eq!(fb::recip_sqrt(9.0), 1.0 / 3.0); + // doesn't check for infinity, just returns a big number + assert_approx_eq!(fb::recip_sqrt(0.0), 2.9727e19); + assert!(fb::recip_sqrt(-1.0).is_nan()); - #[test] - fn recip_sqrt() { - use super::{RecipSqrt, f32}; - assert_approx_eq!(f32::recip_sqrt(2.0), f32::sqrt(0.5), eps = 1e-2); - assert_approx_eq!(f32::recip_sqrt(9.0), 1.0 / 3.0, eps = 1e-3); + // assert_eq!(fb::powf(3.0, 2.0), 9.0); + // assert_eq!(fb::powf(-3.0, 2.0), 9.0); + // assert_eq!(fb::powf(3.0, -2.0), 1.0 / 9.0); + // assert_eq!(fb::powf(-3.0, 3.0), -27.0); + // + // assert_approx_eq!(libm::sin(FRAC_PI_6), 0.5); + // assert_eq!(libm::cos(PI), -1.0); + // + // assert_eq!(libm::exp(1.0), E); + // assert_approx_eq!(libm::exp(2.0), E * E); + // assert_eq!(libm::log2(8.0), 3.0); + // assert!(libm::log2(-1.0).is_nan()); } } diff --git a/core/src/math/grad.rs b/core/src/math/grad.rs new file mode 100644 index 00000000..2df68f17 --- /dev/null +++ b/core/src/math/grad.rs @@ -0,0 +1,266 @@ +use alloc::vec::Vec; +use core::fmt::Debug; + +use super::{Lerp, Parametric, Point2, inv_lerp}; + +/// A position-based color progression that can be used to fill a 2D surface. +#[derive(Clone, Debug, PartialEq)] +pub struct Gradient2 { + /// The shape of the gradient. + pub kind: Kind, + /// The sequence of colors to interpolate between. + pub map: ColorMap, +} + +/// The shape of a gradient. +/// +/// Maps a point to a *t* value used to look up the respective color +/// in the color sequence of a gradient. +#[derive(Copy, Clone, Debug, PartialEq)] +pub enum Kind { + /// A linear, or axial, gradient between two points. + /// + /// Given two points Q and R and an input point P, perpendicularly projects + /// P onto the line crossing Q and R and returns the corresponding *t* + /// value clamped to [0, 1] such that t = 0 at Q and t = 1 at R. + Linear(Pt, Pt), + /// A circularly symmetric gradient of some radius around a point. + /// + /// Given a center point Q, radius *r*, and input point P, maps the distance + /// |P - Q| to *t* values such that: + /// * *t* = 0 when P = Q, and + /// * *t* = 1 when |P - Q| >= *r*. + Radial(Pt, f32), + /// TODO + #[cfg(feature = "fp")] + Conical(Pt), +} + +/// A sequence of (number, color) pairs, mapping t values to colors. +/// The numbers must be in a *nondecreasing* order. +#[derive(Clone, Debug, PartialEq)] +pub struct ColorMap(Vec<(f32, T)>); + +impl Gradient2 { + /// Creates a new gradient. + /// + /// # Panics + /// If there are no stops, or not all the stop values are nondecreasing + pub fn new( + kind: Kind, + stops: impl IntoIterator, + ) -> Self { + Self { + kind, + map: ColorMap::new(stops), + } + } + + /// Returns the value of `self` at the given point. + pub fn eval(&self, p: Point2) -> T { + let t = match self.kind { + Kind::Linear(p0, p1) => (p - p0).scalar_project(&(p1 - p0)), + Kind::Radial(p0, r) => p.distance(&p0) / r, + #[cfg(feature = "fp")] + Kind::Conical(p0) => { + let angle = (p - p0).atan(); + // map negative angles to positive + use super::float::f32; + f32::rem_euclid(angle.to_turns(), 1.0) + } + }; + self.map.eval(t) + } +} + +impl ColorMap { + /// Creates a new color map. + /// + /// # Panics + /// If there are no stops, or not all the stop values are nondecreasing + pub fn new(it: impl IntoIterator) -> Self { + let mut t0 = f32::MIN; + let stops: Vec<_> = it + .into_iter() + .inspect(|&(t, _)| { + assert!(t >= t0, "t values must be nondecreasing"); + t0 = t; + }) + .collect(); + assert!(!stops.is_empty(), "at least one stop is required"); + Self(stops) + } +} + +impl Parametric for ColorMap { + fn eval(&self, t: f32) -> T { + let v = &self.0[..]; + debug_assert!(!v.is_empty(), "failed invariant"); + let res = v.binary_search_by(|(u, _)| u.total_cmp(&t)); + match res { + // t == t_i + Ok(i) => v[i].1.clone(), + // t < t_0 + Err(0) => v[0].1.clone(), + // t > t_n + Err(i) if i == v.len() => v[i - 1].1.clone(), + // 0 < i < len + Err(i) => { + let (t1, v1) = &v[i - 1]; // ok: 0 < i + let (t2, v2) = &v[i]; // ok: i < len + // Remap t such that t=0 -> v1 and t=1 -> v2 + v1.lerp(v2, inv_lerp(t, *t1, *t2)) + } + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::assert_approx_eq; + use crate::math::pt2; + use alloc::vec; + + #[test] + fn linear_gradient() { + let p = pt2(-1.0, 0.0); + let q = pt2(2.0, 1.0); + let g = Gradient2::new(Kind::Linear(p, q), [(0.25, 0.9), (0.6, 0.1)]); + let v = q - p; + + // start point, t=0 + assert_eq!(g.eval(p), 0.9); + + // perpendicular to start point + assert_eq!(g.eval(p + v.perp()), 0.9); + assert_eq!(g.eval(p + 100.0 * v.perp()), 0.9); + + // t < 0 + assert_eq!(g.eval(pt2(-2.0, 0.0)), 0.9); + assert_eq!(g.eval(pt2(-100.0, 0.0)), 0.9); + + // t = 0.25 + assert_eq!(g.eval(p + 0.25 * v), 0.9); + + // t = (0.6 + 0.25)/2 + assert_approx_eq!(g.eval(p + 0.425 * v), 0.5); + assert_approx_eq!(g.eval(p + 0.425 * v + v.perp()), 0.5); + + // t = 0.6 + assert_eq!(g.eval(p + 0.6 * v), 0.1); + + // end point, t = 1 + assert_eq!(g.eval(q), 0.1); + + // t > 1 + assert_eq!(g.eval(pt2(3.0, 1.0)), 0.1); + assert_eq!(g.eval(pt2(100.0, 1.0)), 0.1); + } + + #[cfg(feature = "fp")] + #[test] + fn conical_gradient() { + let g = Gradient2::new( + Kind::Conical(pt2(2.0, -1.0)), + [(0.25, 0.9), (0.6, 0.1)], + ); + + // center point and points to its right with same y should be t=0 + assert_eq!(g.eval(pt2(2.0, -1.0)), 0.9); + assert_eq!(g.eval(pt2(3.0, -1.0)), 0.9); + assert_eq!(g.eval(pt2(100.0, -1.0)), 0.9); + + // t=0.25 corresponds to directly up, should be 0.9 until that + assert_eq!(g.eval(pt2(2.0, 0.0)), 0.9); + assert_eq!(g.eval(pt2(2.0, 100.0)), 0.9); + + // left + assert_approx_eq!(g.eval(pt2(-1.0, -1.0)), 0.328572); + assert_approx_eq!(g.eval(pt2(-100.0, -1.0)), 0.328572); + + // down + assert_eq!(g.eval(pt2(2.0, -2.0)), 0.1); + assert_eq!(g.eval(pt2(2.0, -100.0)), 0.1); + } + + #[test] + fn radial_gradient() { + let g = Gradient2::new( + Kind::Radial(pt2(2.0, -1.0), 2.0), + [(0.25, 0.9), (0.5, 0.1)], + ); + + assert_eq!(g.eval(pt2(2.0, -1.0)), 0.9); // t=0.0 + + assert_approx_eq!(g.eval(pt2(2.0, -1.5)), 0.9); // t=0.25 + assert_approx_eq!(g.eval(pt2(1.5, -1.0)), 0.9); // t=0.25 + + assert_approx_eq!(g.eval(pt2(2.0, -1.75)), 0.5); // t=0.375 + assert_approx_eq!(g.eval(pt2(1.25, -1.0)), 0.5); // t=0.375 + + assert_eq!(g.eval(pt2(2.0, -3.0)), 0.1); // t=0.5 + assert_eq!(g.eval(pt2(0.0, -1.0)), 0.1); // t=0.5 + + assert_eq!(g.eval(pt2(2.0, 2.0)), 0.1); // t=1.0 + assert_eq!(g.eval(pt2(0.0, -1.0)), 0.1); // t=1.0 + } + + #[test] + fn t_less_than_or_eq_to_min_gives_min() { + let g = ColorMap(vec![(0.2, -1.23f32), (0.8, 1.23f32)]); + assert_eq!(g.eval(-10.0), -1.23); + assert_eq!(g.eval(0.0), -1.23); + assert_eq!(g.eval(0.2), -1.23); + } + #[test] + fn t_greater_than_or_eq_to_max_gives_min() { + let g = ColorMap(vec![(0.2, -1.23f32), (0.8, 1.23f32)]); + assert_eq!(g.eval(0.8), 1.23); + assert_eq!(g.eval(1.0), 1.23); + assert_eq!(g.eval(10.0), 1.23); + } + + #[test] + fn t_between_min_max_interpolates() { + let g = ColorMap(vec![(0.2, 1.23f32), (0.6, 2.23f32)]); + + assert_eq!(g.eval(0.2), 1.23); + assert_eq!(g.eval(0.4), 1.73); + assert_eq!(g.eval(0.6), 2.23); + } + + #[test] + fn nonincreasing_t_gives_sharp_change() { + let g = ColorMap(vec![(0.4, -1.23f32), (0.4, 1.23f32)]); + + assert_eq!(g.eval(0.2), -1.23); + assert_eq!(g.eval(0.4f32.next_down()), -1.23); + assert_eq!(g.eval(0.4), 1.23); + assert_eq!(g.eval(0.6), 1.23); + } + + #[test] + fn single_stop_gradient_has_constant_value() { + let g = ColorMap(vec![(0.5, 1.23f32)]); + assert_eq!(g.eval(10.0), 1.23); + assert_eq!(g.eval(-1.0), 1.23); + assert_eq!(g.eval(0.0), 1.23); + assert_eq!(g.eval(0.2), 1.23); + assert_eq!(g.eval(0.8), 1.23); + assert_eq!(g.eval(1.0), 1.23); + assert_eq!(g.eval(10.0), 1.23); + } + + #[test] + #[should_panic(expected = "at least one stop is required")] + fn stops_with_zero_entries_panics() { + _ = ColorMap::<()>::new([]); + } + + #[test] + #[should_panic(expected = "t values must be nondecreasing")] + fn stops_with_nondecreasing_t_panics() { + _ = ColorMap::<()>::new([(0.8, ()), (0.2, ())]); + } +} diff --git a/core/src/math/mat.rs b/core/src/math/mat.rs index 99d07594..a0d472c0 100644 --- a/core/src/math/mat.rs +++ b/core/src/math/mat.rs @@ -550,15 +550,17 @@ impl Mat4 { /// use retrofire_core::assert_approx_eq; /// use retrofire_core::math::*; /// - /// let m = rotate_y(degs(90.0)).then(&translate3(1.0, 2.0, 3.0)); - /// let lin = m.linear(); - /// assert_approx_eq!(lin.apply(&pt3(1.0, 0.0, 0.0)), pt3(0.0, 0.0, -1.0)); + /// let m = scale(splat(5.0)).then(&translate3(1.0, 2.0, 3.0)); + /// let pt = pt3(1.0, -1.0, 0.5); + /// + /// // Only the scale is applied because the translate is not linear + /// assert_approx_eq!(m.linear().apply(&pt), pt3(5.0, -5.0, 2.5)); pub const fn linear(&self) -> Mat3 { let [r, s, t, _] = self.0; mat![ r[0], r[1], r[2]; - s[0], r[1], r[2]; - t[0], r[1], r[2]; + s[0], s[1], s[2]; + t[0], t[1], t[2]; ] } @@ -569,12 +571,25 @@ impl Mat4 { /// use retrofire_core::math::*; /// /// let trans = vec3(1.0, 2.0, 3.0); - /// let m = rotate_y(degs(45.0)).then(&translate(trans)); + /// let m = scale(splat(5.0)).then(&translate(trans)); /// assert_eq!(m.translation(), trans); pub const fn translation(&self) -> Vec3 { vec3(self.0[0][3], self.0[1][3], self.0[2][3]) } + /// Returns the translation column vector of `self` as a point. + /// + /// # Example + /// ``` + /// use retrofire_core::math::*; + /// + /// let trans = vec3(1.0, 2.0, 3.0); + /// let m = scale(splat(5.0)).then(&translate(trans)); + /// assert_eq!(m.origin(), pt3(1.0, 2.0, 3.0)); + pub const fn origin(&self) -> Point3 { + self.translation().to_pt() + } + /// Returns the determinant of `self`. /// /// Given a matrix M, @@ -1154,6 +1169,23 @@ pub fn rotate_z(a: Angle) -> Mat4 { ] } +/// Returns a rotation matrix based on three local (intrinsic) angles: +/// pitch, yaw, and roll. +/// +/// The pitch (elevation) angle controls rotation about the local lateral (x) +/// axis, the yaw (bearing) angle about the local vertical (y) axis, and the +/// roll (bank) angle about the local longitudinal (z) axis. These angles are +/// also known to as Tait–Bryan angles. +/// +/// See also: [`cam::PitchYawRoll`][crate::render::cam::PitchYawRoll]. +#[cfg(feature = "fp")] +pub fn rotate_pyr(pitch: Angle, yaw: Angle, roll: Angle) -> Mat4 { + let p = rotate_x(pitch); + let y = rotate_y(yaw); + let r = rotate_z(roll); + p.then(&y).then(&r) +} + /// Returns a matrix applying a 2D rotation by an angle. #[cfg(feature = "fp")] pub fn rotate2(a: Angle) -> Mat3 { @@ -1283,7 +1315,7 @@ mod tests { #[allow(unused)] const O: Vec3 = Vec3::new([0.0; 3]); - mod mat2x2 { + mod mat2 { use super::*; #[test] @@ -1317,7 +1349,7 @@ mod tests { } } - mod mat3x3 { + mod mat3 { use super::*; const MAT: Mat3 = mat![ @@ -1454,6 +1486,24 @@ mod tests { assert_eq!(MAT.col_vec(3), [3.0, 13.0, 23.0, 33.0].into()); } + #[test] + fn linear_part() { + assert_eq!( + MAT.linear(), + mat![ + 0.0, 1.0, 2.0; + 10.0, 11.0, 12.0; + 20.0, 21.0, 22.0; + ] + ); + } + + #[test] + fn translation_part() { + assert_eq!(MAT.translation(), vec3(3.0, 13.0, 23.0)); + assert_eq!(MAT.origin(), pt3(3.0, 13.0, 23.0)); + } + #[test] fn composition() { let tr = translate3(1.0, 2.0, 3.0).to::(); @@ -1465,9 +1515,14 @@ mod tests { assert_eq!(tr_sc, sc.compose(&tr)); assert_eq!(sc_tr, tr.compose(&sc)); - let o = ::origin(); - assert_eq!(tr_sc.apply(&o.to()), pt3::<_, B1>(3.0, 4.0, 3.0)); - assert_eq!(sc_tr.apply(&o.to()), pt3::<_, B2>(1.0, 2.0, 3.0)); + assert_eq!( + tr_sc.apply(&Point3::origin()), + pt3::<_, B1>(3.0, 4.0, 3.0) + ); + assert_eq!( + sc_tr.apply(&Point3::origin()), + pt3::<_, B2>(1.0, 2.0, 3.0) + ); } #[test] @@ -1606,7 +1661,7 @@ mod tests { } #[test] - fn from_basis() { + fn from_linear_basis() { let m = Mat4::from_linear(Y, 2.0 * Z, -3.0 * X); assert_eq!(m.apply(&X), Y); @@ -1636,7 +1691,7 @@ mod tests { fn orientation_no_op() { let m = orient_y(Y, X); - assert_eq!(m.apply(&X), X); + assert_approx_eq!(m.apply(&X), X); assert_eq!(m.apply(&X.to_pt()), X.to_pt()); assert_eq!(m.apply(&Y), Y); @@ -1650,7 +1705,7 @@ mod tests { fn orientation_y_to_z() { let m = orient_y(Z, X); - assert_eq!(m.apply(&X), X); + assert_approx_eq!(m.apply(&X), X); assert_eq!(m.apply(&X.to_pt()), X.to_pt()); assert_eq!(m.apply(&Y), Z); @@ -1664,7 +1719,7 @@ mod tests { fn orientation_z_to_y() { let m = orient_z(Y, X); - assert_eq!(m.apply(&X), X); + assert_approx_eq!(m.apply(&X), X); assert_eq!(m.apply(&X.to_pt()), X.to_pt()); assert_eq!(m.apply(&Y), -Z); diff --git a/core/src/math/param.rs b/core/src/math/param.rs index 8941df56..4fd6d4e8 100644 --- a/core/src/math/param.rs +++ b/core/src/math/param.rs @@ -4,6 +4,7 @@ use super::Lerp; /// Represents a single-variable parametric curve. // TODO More documentation +// TODO Associated type instead of parameter? pub trait Parametric { /// Returns the value of `self` at `t`. /// diff --git a/core/src/math/point.rs b/core/src/math/point.rs index cd642498..3dc9a513 100644 --- a/core/src/math/point.rs +++ b/core/src/math/point.rs @@ -109,7 +109,6 @@ impl Point<[f32; N], Real> { /// let y4 = pt2(0.0, 4.0); /// assert_eq!(x3.distance(&y4), 5.0); /// ``` - #[cfg(feature = "fp")] #[inline] pub fn distance(&self, other: &Self) -> f32 { self.sub(other).len() @@ -157,6 +156,42 @@ impl Point<[f32; N], Real> { pub fn clamp(&self, min: &Self, max: &Self) -> Self { array::from_fn(|i| self.0[i].clamp(min.0[i], max.0[i])).into() } + + /// Returns `self`, moved towards another point by an offset. + /// + /// If the distance to `other` is less than `d`, returns `other`. That is, + /// this method never "overshoots" the target. Negative values of `d` result + /// in moving away from `other`. If `self` equals `other`, returns `other` + /// independent of the value of `d`. + /// + /// To translate by a *relative* offset instead, use [`lerp`][super::Lerp::lerp]. + /// + /// # Examples + /// ``` + /// use retrofire_core::math::{Point2, pt2}; + /// + /// let a: Point2 = pt2(0.0, 0.0); + /// let b: Point2 = pt2(4.0, 3.0); + /// // Move two units along the line y = 3x/4: + /// assert_eq!(a.approach(&b, 2.0), pt2(1.6, 1.2)); + /// // Movement is clamped to b: + /// assert_eq!(a.approach(&b, 10.0), b); + /// // Negative values of `d` move away from b: + /// assert_eq!(a.approach(&b, -1.0), pt2(-0.8, -0.6)); + /// // Approaching the point itself does nothing: + /// assert_eq!(a.approach(&a, 1.0), a); + /// assert_eq!(a.approach(&a, -1.0), a); + /// ``` + #[must_use] + pub fn approach(&self, other: &Self, d: f32) -> Self { + let v = *other - *self; + let l = v.len(); + if d < l && l != 0.0 { + *self + d / l * v + } else { + *other + } + } } impl Point<[Sc; 2], Real<2, B>> { @@ -436,6 +471,7 @@ mod tests { mod f32 { use super::*; + use crate::assert_approx_eq; const pt2: fn(f32, f32) -> Point2 = super::pt2; const pt3: fn(f32, f32, f32) -> Point3 = super::pt3; @@ -485,10 +521,12 @@ mod tests { ); } #[test] - #[cfg(feature = "fp")] fn point_point_distance() { - assert_eq!(pt2(1.0, -1.0).distance(&pt2(-2.0, 3.0)), 5.0); - assert_eq!(pt3(1.0, -3.0, 2.0).distance(&pt3(-2.0, 3.0, 4.0)), 7.0); + assert_approx_eq!(pt2(1.0, -1.0).distance(&pt2(-2.0, 3.0)), 5.0); + assert_approx_eq!( + pt3(1.0, -3.0, 2.0).distance(&pt3(-2.0, 3.0, 4.0)), + 7.0 + ); } #[test] fn point2_clamp() { @@ -549,7 +587,7 @@ mod tests { mod u32 { use super::*; - const pt2: fn(u32, u32) -> Point2u = super::super::pt2; + const pt2: fn(u32, u32) -> Point2u = super::pt2; #[test] fn vector_addition() { diff --git a/core/src/math/rand.rs b/core/src/math/rand.rs index 5b2fbf9a..cb37df7c 100644 --- a/core/src/math/rand.rs +++ b/core/src/math/rand.rs @@ -10,6 +10,8 @@ use super::{Angle, Color, Point, Point2, Point3, Vec2, Vec3, Vector, rads}; pub type DefaultRng = Xorshift64; +pub const DEFAULT_RNG: DefaultRng = Xorshift64(Xorshift64::DEFAULT_SEED); + /// Trait for generating values sampled from a probability distribution. pub trait Distrib { /// The type of the elements of the sample space of `Self`, also called @@ -138,8 +140,8 @@ impl Xorshift64 { /// # Panics /// /// If `seed` equals 0. - pub fn from_seed(seed: u64) -> Self { - assert_ne!(seed, 0, "xorshift seed cannot be zero"); + pub const fn from_seed(seed: u64) -> Self { + assert!(seed != 0, "xorshift seed cannot be zero"); Self(seed) } @@ -172,7 +174,7 @@ impl Xorshift64 { /// Successive calls to this function (with the same `self`) will yield /// every value in the interval [1, 264) exactly once before /// starting to repeat the sequence. - pub fn next_bits(&mut self) -> u64 { + pub const fn next_bits(&mut self) -> u64 { let Self(x) = self; *x ^= *x << 13; *x ^= *x >> 7; @@ -207,8 +209,7 @@ impl Default for Xorshift64 { /// assert_eq!(g.next_bits(), 11039719294064252060); /// ``` fn default() -> Self { - // Random 64-bit prime - Self::from_seed(Self::DEFAULT_SEED) + DEFAULT_RNG } } @@ -461,7 +462,6 @@ where } } -#[cfg(feature = "fp")] impl Distrib for UnitCircle { type Sample = Vec2; @@ -493,7 +493,7 @@ impl Distrib for VectorsOnUnitDisk { /// let rng = &mut DefaultRng::default(); /// /// let vec = VectorsOnUnitDisk.sample(rng); - /// assert!(vec.len_sqr() <= 1.0); + /// assert!(vec.len() <= 1.0); /// ``` fn sample(&self, rng: &mut DefaultRng) -> Vec2 { let d = Uniform([-1.0f32; 2]..[1.0; 2]); @@ -507,7 +507,6 @@ impl Distrib for VectorsOnUnitDisk { } } -#[cfg(feature = "fp")] impl Distrib for UnitSphere { type Sample = Vec3; @@ -520,7 +519,7 @@ impl Distrib for UnitSphere { /// let rng = &mut DefaultRng::default(); /// /// let vec = UnitSphere.sample(rng); - /// assert_approx_eq!(vec.len_sqr(), 1.0); + /// assert_approx_eq!(vec.len(), 1.0); /// ``` fn sample(&self, rng: &mut DefaultRng) -> Vec3 { let d = Uniform([-1.0; 3]..[1.0; 3]); @@ -679,7 +678,6 @@ mod tests { assert_eq!(approx_100, 82); } - #[cfg(feature = "fp")] #[test] fn unit_circle() { use crate::assert_approx_eq; @@ -695,7 +693,6 @@ mod tests { } } - #[cfg(feature = "fp")] #[test] fn unit_sphere() { use crate::assert_approx_eq; diff --git a/core/src/math/space.rs b/core/src/math/space.rs index ccec9ac0..93396b4f 100644 --- a/core/src/math/space.rs +++ b/core/src/math/space.rs @@ -2,11 +2,13 @@ //! //! TODO -use core::fmt::{Debug, Formatter}; -use core::iter::zip; -use core::marker::PhantomData; +use core::{ + fmt::{Debug, Formatter}, + iter::zip, + marker::PhantomData, +}; -use crate::math::vary::{Iter, Vary, ZDiv}; +use super::vary::{Iter, Vary, ZDiv}; /// Trait for types representing elements of an affine space. /// @@ -187,7 +189,7 @@ impl Affine for u32 { } } -impl Vary for V +impl Vary for V where Self: Affine + Clone> + ZDiv, { diff --git a/core/src/math/spline.rs b/core/src/math/spline.rs index af2f55e5..7dfb5e80 100644 --- a/core/src/math/spline.rs +++ b/core/src/math/spline.rs @@ -1,18 +1,20 @@ //! Bézier curves and splines. - use alloc::vec::Vec; -use core::array; +use core::{array::from_fn, fmt::Debug, marker::PhantomData}; use crate::geom::{Polyline, Ray}; +use crate::mat; -use super::{Affine, Lerp, Linear, Parametric, Vector}; +use super::{ + Affine, Lerp, Linear, Mat4, Parametric, Point, Vary, Vector, inv_lerp, + space::Real, +}; /// A cubic Bézier curve, defined by four control points. /// /// TODO More info about Béziers /// /// ```text -/// /// p1 /// \ ____ /// \ _-´ `--_ p3 @@ -23,9 +25,68 @@ use super::{Affine, Lerp, Linear, Parametric, Vector}; /// p0 \ /// p2 /// ``` -#[derive(Debug, Clone, Eq, PartialEq)] +#[derive(Debug, Copy, Clone, Eq, PartialEq)] pub struct CubicBezier(pub [T; 4]); +/// A cubic Hermite curve, defined by two control points and the velocity +/// vectors at the control points. +/// +/// Hermite curves are closely related to Bézier curves, via the identity +/// +/// H(p0, d0, p1, d1) = B(p0, p0 + d0/3, p1 - d1/3, p1), +/// +/// or, equivalently, +/// +/// B(p0, p1, p2, p3) = H(p0, 3(p1 - p0), p3, 3(p2 - p3)). +/// +/// ```text +/// d0 +/// ^ +/// \ +/// \ ____ +/// \ _-´ `--_ p1 +/// \ / `-_ \ +/// \| `-_ |\ +/// \ `-__ / \ +/// p0 `---_____--´ \ +/// \ +/// \ +/// v +/// d1 +/// ``` +#[derive(Debug, Copy, Clone, Eq, PartialEq)] +pub struct CubicHermite(pub [P; 2], pub [D; 2]); + +/// A piecewise curve composed of concatenated [cubic Bézier curves][CubicBezier]. +#[derive(Debug, Clone, Eq, PartialEq)] +pub struct BezierSpline(Vec); + +/// A piecewise curve composed of concatenated [cubic Hermite curves][CubicHermite]. +#[derive(Debug, Clone, Eq, PartialEq)] +// HACK: The PhantomData field only exists to force the derive impls +// to include the correct `T::Diff: Trait` bounds +pub struct HermiteSpline(Vec>, PhantomData); + +/// A piecewise curve composed of concatenated cubic curves. +#[derive(Debug, Clone, Eq, PartialEq)] +pub struct CatmullRomSpline(Vec); + +/// A piecewise curve composed of concatenated cubic curves. +#[derive(Debug, Clone, Eq, PartialEq)] +pub struct BSpline(Vec); + +/// Euclidean (arc-length) parameterization for splines. +/// +/// Instead of *t* ∈ [0, 1], a `Euclidean` spline is parameterized by *s*, +/// the actual Euclidean distance travelled along the curve. This makes it easy, +/// for example, to position objects along the spline at regular intervals, or +/// to travel along the spline at a desired speed. +/// +/// Because arc-length parameterization is not generally possible in closed form, +/// this implementation uses a look-up table to map *s* values to *t* values, +/// interpolating linearly between entries. +pub struct Euclidean(Spl, Vec<(f32, f32)>); + /// Interpolates smoothly from 0.0 to 1.0 as `t` goes from 0.0 to 1.0. /// /// Returns 0 for all `t` <= 0 and 1 for all `t` >= 1. Has a continuous @@ -58,49 +119,163 @@ where } } -impl CubicBezier +/// Approximates a curve as a chain of line segments. +/// +/// Recursively subdivides the curve into two half-curves, stopping once +/// the approximation error is less than `error`. +/// +/// # Examples +/// ``` +/// use retrofire_core::math::{BezierSpline, Point2, pt2}; +/// use retrofire_core::math::spline::approximate; +/// +/// let curve = BezierSpline::::new( +/// [pt2(0.0, 0.0), pt2(0.0, 1.0), pt2(1.0, 1.0), pt2(1.0, 0.0)] +/// ); +/// // Find an approximation with error less than 0.01 +/// let approx = approximate(&curve, 0.01); +/// +/// // Euclidean length of the polyline approximation +/// assert_eq!(approx.len(), 1.9969313); +/// +/// // Number of line segments used by the approximation +/// assert_eq!(approx.0.len(), 17); +/// ``` +/// +/// # Panics +/// If `err` ≤ 0. +pub fn approximate( + curve: &impl Parametric, + error: f32, +) -> Polyline where - T: Affine> + Clone, + T: Affine>, { - /// Evaluates the value of `self` at `t`. + assert!(error > 0.0); + approximate_with(curve, &|e: &T::Diff| e.len_sqr() < error * error) +} + +/// Approximates a curve as a chain of line segments. +/// +/// Recursively subdivides the curve into two half-curves, stopping once +/// the approximation error is small enough, as determined by the `halt` +/// function. +/// +/// Given a curve segment between some points `p` and `r`, the parameter +/// passed to `halt` is the distance to the real midpoint `q` from its +/// linear approximation `q'`. If `halt` returns `true`, the line segment +/// `pr` is returned as the approximation of this curve segment, otherwise +/// the bisection continues. +/// +/// Note that this heuristic does not work well in certain edge cases +/// (consider, for example, an S-shaped curve where `q'` is very close +/// to `q`, yet a straight line would be a poor approximation). However, +/// in practice it tends to give reasonable results. +/// +/// ```text +/// ___--- q ---___ +/// _--´ | `--_ +/// _--´ | `--_ +/// _-p ------------- q' ------------ r-_ +/// _-´ `-_ +/// ``` +/// +/// # Examples +/// ``` +/// use retrofire_core::math::{BezierSpline, Point2, pt2}; +/// use retrofire_core::math::spline::approximate_with; +/// +/// let curve = BezierSpline::::new( +/// [pt2(0.0, 0.0), pt2(0.0, 1.0), pt2(1.0, 1.0), pt2(1.0, 0.0)] +/// ); +/// // Find an approximation with error less than 0.01 +/// let approx = approximate_with(&curve, |err| err.len_sqr() < 0.01 * 0.01); +/// +/// // Euclidean length of the polyline approximation +/// assert_eq!(approx.len(), 1.9969313); +/// +/// // Number of line segments used by the approximation +/// assert_eq!(approx.0.len(), 17); +/// ``` +pub fn approximate_with>>( + curve: &impl Parametric, + halt: impl Fn(&T::Diff) -> bool, +) -> Polyline { + let mut res = Vec::new(); + do_approx(curve, 0.0, 1.0, 10, &halt, &mut res); + res.push(curve.eval(1.0)); + Polyline(res) +} + +fn do_approx>>( + c: &impl Parametric, + a: f32, + b: f32, + max_dep: u32, + halt: &impl Fn(&T::Diff) -> bool, + accum: &mut Vec, +) { + let mid = a.midpoint(b); + + let ap = c.eval(a); + let bp = c.eval(b); + + let real = c.eval(mid); + let approx = ap.add(&bp.sub(&ap).mul(0.5)); + + if max_dep == 0 || halt(&real.sub(&approx)) { + accum.push(ap); + } else { + do_approx(c, a, mid, max_dep - 1, halt, accum); + do_approx(c, mid, b, max_dep - 1, halt, accum); + } +} + +// +// Inherent impls +// + +impl CubicBezier { + /// Returns the point of `self` at the given *t* value. Uses + /// [De Casteljau's algorithm][1]. /// - /// For t < 0, returns the first control point. For t > 1, returns the last - /// control point. Uses [De Casteljau's algorithm][1]. + /// Values of *t* outside the interval [0, 1] are accepted and extrapolate + /// the curve beyond the control points. /// /// [1]: https://en.wikipedia.org/wiki/De_Casteljau%27s_algorithm pub fn eval(&self, t: f32) -> T { let [p0, p1, p2, p3] = &self.0; - step(t, p0, p3, |t| { - let p01 = p0.lerp(p1, t); - let p12 = p1.lerp(p2, t); - let p23 = p2.lerp(p3, t); - p01.lerp(&p12, t).lerp(&p12.lerp(&p23, t), t) - }) + let p01 = p0.lerp(p1, t); + let p12 = p1.lerp(p2, t); + let p23 = p2.lerp(p3, t); + p01.lerp(&p12, t).lerp(&p12.lerp(&p23, t), t) } +} - /// Evaluates the value of `self` at `t`. - /// - /// For t < 0, returns the first control point. For t > 1, returns the last - /// control point. +impl CubicBezier +where + T: Affine> + Clone, +{ + /// Returns the point of `self` at the given *t* value. /// - /// Directly evaluates the cubic. Faster but possibly less numerically - /// stable than [`Self::eval`]. + /// Directly evaluates the cubic polynomial. Faster but possibly less + /// numerically stable than [`Self::eval`]. Values of *t* outside the + /// interval [0, 1] are accepted and extrapolate the curve beyond the + /// control points. pub fn fast_eval(&self, t: f32) -> T { - let [p0, .., p3] = &self.0; - step(t, p0, p3, |t| { - // Add a linear combination of the three coefficients - // to `p0` to get the result - let [co3, co2, co1] = self.coefficients(); - p0.add(&co3.mul(t).add(&co2).mul(t).add(&co1).mul(t)) - }) + // Add a linear combination of the three coefficients + // to `p0` to get the result + let p0 = &self.0[0]; + let [co3, co2, co1] = self.coefficients(); + p0.add(&co3.mul(t).add(&co2).mul(t).add(&co1).mul(t)) } - /// Returns the tangent, or direction vector, of `self` at `t`. + /// Returns the velocity vector of `self` at the given *t* value. /// - /// Clamps `t` to the range [0, 1]. - pub fn tangent(&self, t: f32) -> T::Diff { + /// Values of *t* outside the interval [0, 1] are accepted and extrapolate + /// the curve beyond the control points. + pub fn velocity(&self, t: f32) -> T::Diff { let [p0, p1, p2, p3] = &self.0; - let t = t.clamp(0.0, 1.0); // 3 (3 (p1 - p2) + (p3 - p0)) * t^2 // + 6 ((p0 - p1 + p2 - p1) * t @@ -157,181 +332,450 @@ where } } -/// A curve composed of one or more concatenated -/// [cubic Bézier curves][CubicBezier]. -#[derive(Debug, Clone, Eq, PartialEq)] -pub struct BezierSpline(Vec); +impl

CubicHermite +where + P: Affine> + Clone, +{ + // Characteristic matrix M + // 1.0, 0.0, 0.0, 0.0; + // 0.0, 1.0, 0.0, 0.0; + // -3.0, -2.0, 3.0, -1.0; + // 2.0, 1.0, -2.0, 1.0; + + /// Returns the point of `self` at the given *t* value. + /// + /// Values of *t* outside the interval [0, 1] are accepted and extrapolate + /// the curve beyond the control points. + pub fn eval(&self, t: f32) -> P { + let Self([p0, p1], [d0, d1]) = self; + let [_0, t1, t2, t3] = [1.0, t, t * t, t * t * t]; + + // b = ts * M + let _0 = 1.0 - 3.0 * t2 + 2.0 * t3; // = 1 - b2 + let b1 = t1 - 2.0 * t2 + t3; + let b2 = 3.0 * t2 - 2.0 * t3; + let b3 = -t2 + t3; + + // H(t) = b * P + + // b0 * p0 + b1 * d0 + b2 * p1 + b3 * d1 + // = b0 * p0 + b2 * p1 // Affine: b0 + b2 = 1: lerp + // + b1 * d0 + b3 * d1 // Linear + + // b0 * p0 + b2 * p1 + // = (1 - b2) * p0 + b2 * p1 + // = p0 + b2 * (p1 - p0) + + p0.add(&p1.sub(&p0).mul(b2)) // Affine part + .add(&d0.mul(b1).add(&d1.mul(b3))) // Linear part + } + + /// Returns the velocity vector of `self` at the given *t* value. + /// + /// Values of *t* outside the interval [0, 1] are accepted and extrapolate + /// the curve beyond the control points. + pub fn velocity(&self, t: f32) -> P::Diff { + let Self([p0, p1], [d0, d1]) = self; + // Derivatives of the powers of t + let [_0, _1, t2, t3] = [0.0, 1.0, 2.0 * t, 3.0 * t * t]; + + // b = ts * M + let _0 = 0.0 - 3.0 * t2 + 2.0 * t3; // = -b2 + let b1 = 1.0 - 2.0 * t2 + t3; + let b2 = 3.0 * t2 - 2.0 * t3; + let b3 = -t2 + t3; + + // H(t) = b * P + + // b0 * p0 + b1 * d0 + b2 * p1 + b3 * d1 + // = b0 * p0 + b2 * p1 // b0 = -b2 + // + b1 * d0 + b3 * d1 // Linear + + // b0 * p0 + b2 * p1 + // = -b2 * p0 + b2 * p1 + // = b2 * p1 - b2 * p0 + // = b2 * (p1 - p0) + + // Only vectors as expected: + // b2·(p1 - p0) + b1·d0 + b3·d1 + p1.sub(&p0) + .mul(b2) + .add(&d0.mul(b1).add(&d1.mul(b3))) + } +} impl BezierSpline where T: Affine + Clone> + Clone, { - /// Creates a Bézier spline from the given control points. The number of - /// elements in `pts` must be 3n + 1 for some positive integer n. + /// Creates a Bézier spline from the given control points. /// - /// Consecutive points in `pts` make up Bézier curves such that: - /// * `pts[0..=3]` define the first curve, - /// * `pts[3..=6]` define the second curve, + /// The number of elements in `pts` must be 3*n* + 1 for some positive integer + /// *n*. Consecutive points in `pts` make up Bézier curves such that: + /// * (p0, p1, p2, p3) define the first curve, + /// * (p3, p4, p5, p6) define the second curve, /// /// and so on. /// /// # Panics - /// If `pts.len() < 4` or if `pts.len() % 3 != 1`. - pub fn new(pts: &[T]) -> Self { + /// If the number of points *n* < 4 or if *n* ≠ 1 (mod 3). + pub fn new(pts: impl IntoIterator) -> Self { + let pts = pts.into_iter().collect::>(); + let len = pts.len(); assert!( - pts.len() >= 4 && pts.len() % 3 == 1, - "length must be 3n+1 for some integer n > 0, was {}", - pts.len() + len >= 4 && len % 3 == 1, + "length must be 3n+1 for some integer n > 0, was {len}", ); - Self(pts.to_vec()) + Self(pts) } - /// Constructs a Bézier spline - pub fn from_rays(rays: I) -> Self - where - I: IntoIterator>, - { - let mut rays = rays.into_iter().peekable(); - let mut first = true; - let mut pts = Vec::with_capacity(2 * rays.size_hint().0); - while let Some(ray) = rays.next() { - if !first { - pts.push(ray.eval(-1.0)); - } - first = false; - pts.push(ray.0.clone()); - if rays.peek().is_some() { - pts.push(ray.eval(1.0)); - } - } - Self::new(&pts) + /// Constructs a Bézier spline from (position, tangent) pairs. + /// + /// For each pair of consecutive rays (P, **v**) and (Q, **u**), the result + /// contains one cubic Bézier curve segment with control points (P, P + + /// **v**, Q - **u**, Q). + /// + /// # Panics + /// If the number of rays < 2. + pub fn from_rays(rays: impl IntoIterator>) -> Self { + let pts: Vec<_> = rays + .into_iter() + .flat_map(|Ray(p, d)| [p.add(&d.neg()), p.clone(), p.add(&d)]) + .collect(); + Self::new(pts[1..pts.len() - 1].into_iter().cloned()) } - /// Evaluates `self` at position `t`. + /// Returns the point of `self` at the given *t* value. /// - /// Returns the first point if `t < 0` and the last point if `t > 1`. + /// Values of *t* outside the interval [0, 1] are accepted and extrapolate + /// the curve beyond the control points. pub fn eval(&self, t: f32) -> T { - // invariant self.0.len() != 0 -> last always exists - step(t, &self.0[0], self.0.last().unwrap(), |t| { - let (t, seg) = self.segment(t); - CubicBezier(seg).fast_eval(t) - }) + let (u, seg) = self.segment(t); + seg.fast_eval(u) } - /// Returns the tangent of `self` at `t`. + /// Returns the velocity vector of `self` at the given *t* value. /// - /// Clamps `t` to the range [0, 1]. - pub fn tangent(&self, t: f32) -> T::Diff { - let (t, seg) = self.segment(t); - CubicBezier(seg).tangent(t) + /// Values of *t* outside the interval [0, 1] are accepted and extrapolate + /// the curve beyond the control points. + pub fn velocity(&self, t: f32) -> T::Diff { + let (u, seg) = self.segment(t); + seg.velocity(u) } - fn segment(&self, t: f32) -> (f32, [T; 4]) { - let segs = ((self.0.len() - 1) / 3) as f32; - // TODO use floor and make the code cleaner - let seg = ((t * segs) as u32 as f32).min(segs - 1.0); - let t2 = t * segs - seg; - let idx = 3 * (seg as usize); - (t2, array::from_fn(|k| self.0[idx + k].clone())) + /// Returns the list of control points of `self`. + pub fn control_points(&self) -> &[T] { + &self.0 } - /// Approximates `self` as a chain of line segments. + /// Returns the spline segment and local *t* value corresponding to + /// the given global *t* value. /// - /// Recursively subdivides the curve into two half-curves, stopping once - /// the approximation error is less than `error`. - /// - /// # Examples - /// ``` - /// use retrofire_core::math::{BezierSpline, vec2, Vec2}; + fn segment(&self, t: f32) -> (f32, CubicBezier) { + // Consecutive segments share an endpoint: + // [B0 B1 B2 B3] + // [B3 B4 B5 B6] + // [B6 B7 B8 B9] + // [B9 ... + // If the number of segs is n, the number of control points is 3n + 1, + // thus if the number of points is l, the number of segs is (l - 1) / 3. + let num_segs = (self.0.len() - 1) / 3; + // Rescale from [0, 1] to [0, num_segs] + let t = t * num_segs as f32; + use super::float::f32; + // Calculate the segment index. + let seg_i = (t as usize).min(num_segs - 1); + // The leftover part is the local t value. This is the fractional part + // for 0 <= t < segs. t = segs maps to u = 1 of the last subsegment. + // Values of t < 0 or t > segs result in u < 0 or u > 1 and extrapolate + // beyond the first or last subsegment, respectively. + let u = t - seg_i as f32; + // Index of the first control point of the segment + let i = 3 * seg_i; + let seg = from_fn(|j| self.0[i + j].clone()); + (u, CubicBezier(seg)) + } +} + +impl HermiteSpline +where + T: Affine + Clone> + Clone, +{ + /// Creates a new Hermite spline from a sequence of rays. /// - /// let curve = BezierSpline::::new( - /// &[vec2(0.0, 0.0), vec2(0.0, 1.0), vec2(1.0, 1.0), vec2(1.0, 0.0)] - /// ); - /// let approx = curve.approximate(0.01); - /// assert_eq!(approx.0.len(), 17); - /// ``` + /// Each ray (Pi, **v**i) makes up a point + /// Pi on the curve and the velocity (velocity) vector + /// **v**i of the curve at that point. Thus, + /// the ray lies tangent to the curve at point Pi. /// /// # Panics - /// If `err` ≤ 0. - pub fn approximate(&self, error: f32) -> Polyline - where - T: Affine>, - { - assert!(error > 0.0); - self.approximate_with(&|e: &T::Diff| e.len_sqr() < error * error) + /// If `rays` has fewer than two items. + pub fn new(rays: impl IntoIterator>) -> Self { + let rays: Vec<_> = rays.into_iter().collect(); + assert!( + rays.len() >= 2, + "a Hermite spline requires at least two points and two vectors" + ); + Self(rays, PhantomData) } - /// Approximates `self` as a chain of line segments. + /// Returns the subsegment and local *t* value corresponding to the given + /// global *t* value. + fn segment(&self, t: f32) -> (f32, CubicHermite) { + // Scale from [0, 1] to [0, len-1] + let t = t * (self.0.len() - 1) as f32; + // Calculate the index of the subsegment. There are len-1 subsegments: + // (0, 1), (1,2), ..., (len-2, len-1). + let i = (t as usize).min(self.0.len() - 2); + // The leftover part is the local t value. This is the fractional part + // for 0 <= t < len-1. t = len-1 maps to u = 1 of the last subsegment. + // Values of t < 0 or t > len-1 result in u < 0 or u > 1 and extrapolate + // beyond the first or last subsegment, respectively. + let u = t - i as f32; + // Ok: i <= self.0.len() - 2 + let Ray(p0, d0) = self.0[i].clone(); + let Ray(p1, d1) = self.0[i + 1].clone(); + (u, CubicHermite([p0, p1], [d0, d1])) + } + + /// Returns the point of `self` at the given *t* value. /// - /// Recursively subdivides the curve into two half-curves, stopping once - /// the approximation error is small enough, as determined by the `halt` - /// function. + /// Values of *t* outside the interval [0, 1] are accepted and extrapolate + /// the curve beyond the control points. + pub fn eval(&self, t: f32) -> T { + let (u, seg) = self.segment(t); + seg.eval(u) + } + + /// Returns the velocity vector of `self` + /// at the given *t* value. /// - /// Given a curve segment between some points `p` and `r`, the parameter - /// passed to `halt` is the distance to the real midpoint `q` from its - /// linear approximation `q'`. If `halt` returns `true`, the line segment - /// `pr` is returned as the approximation of this curve segment, otherwise - /// the bisection continues. + /// Values of *t* outside the interval [0, 1] are accepted and extrapolate + /// the curve beyond the control points. + pub fn velocity(&self, t: f32) -> T::Diff { + let (u, seg) = self.segment(t); + seg.velocity(u) + } +} + +impl CatmullRomSpline +where + T: Affine> + Clone, +{ + const _CHAR_MAT: Mat4 = mat![ + 0.0, 1.0, 0.0, 0.0; + -0.5, 0.0, 0.5, 0.0; + 1.0, -2.5, 2.0, -0.5; + -0.5, 1.5, -1.5, 0.5; + ]; + + pub fn new(pts: impl IntoIterator) -> Self { + let pts: Vec<_> = pts.into_iter().collect(); + assert!( + pts.len() >= 4, + "a Catmull–Rom spline requires at least four points" + ); + Self(pts) + } + + /// Returns the point of `self` at the given *t* value. /// - /// Note that this heuristic does not work well in certain edge cases - /// (consider, for example, an S-shaped curve where `q'` is very close - /// to `q`, yet a straight line would be a poor approximation). However, - /// in practice it tends to give reasonable results. + /// Values of *t* outside the interval [0, 1] are accepted and extrapolate + /// the curve beyond the control points. + pub fn eval(&self, t: f32) -> T { + let (t, [p0, p1, p2, p3]) = crb_segment(&self.0, t); + let [_0, t1, t2, t3] = [1.0, t, t * t, t * t * t]; + + let _0 = (-t1 + 2.0 * t2 - t3) / 2.0; + let b1 = (2.0 - 5.0 * t2 + 3.0 * t3) / 2.0; + let b2 = (t1 + 4.0 * t2 - 3.0 * t3) / 2.0; + let b3 = (-t2 + t3) / 2.0; + + // b0 + b1 + b2 + b3 = 1 + // b0 = 1 - b1 - b2 - b3 + + // b0·P0 + b1·P1 + b2·P2 + b3·P3 + // = (1 - b1 - b2 - b3)·P0 + b1·P1 + b2·P2 + b3·P3 + // = P0 - b1·P0 - b2·P0 - b3·P0 + b1·P1 + b2·P2 + b3·P3 + // = P0 + b1·(P1 - P0) + b2·(P2 - P0) + b3·(P3 - P0) + + p0.add(&p1.sub(&p0).mul(b1)) + .add(&p2.sub(&p0).mul(b2)) + .add(&p3.sub(&p0).mul(b3)) + } + + /// Returns the gradient, or velocity, vector of `self` at the given *t* + /// value. /// - /// ```text - /// ___--- q ---___ - /// _--´ | `--_ - /// _--´ | `--_ - /// _-p ------------- q' ------------ r-_ - /// _-´ `-_ - /// ``` + /// Values of *t* outside the interval [0, 1] are accepted and extrapolate + /// the curve beyond the control points. + pub fn gradient(&self, t: f32) -> T::Diff { + let (t, [p0, p1, p2, p3]) = crb_segment(&self.0, t); + let [_0, _1, t2, t3] = [0.0, 1.0, 2.0 * t, 3.0 * t * t]; + + // ⎛ b0 ⎞ ⎛ -1 + 2·t2 - t3 ⎞ + // ⎜ b1 ⎟ = 1/2 ⎜ - 5·t2 + 3·t3 ⎟ + // ⎜ b2 ⎟ ⎜ 1 + 4·t2 - 3·t3 ⎟ + // ⎝ b3 ⎠ ⎝ - t2 + t3 ⎠ + + let b1 = -2.5 * t2 + 1.5 * t3; + let b2 = 0.5 + 2.0 * t2 - 1.5 * t3; + let b3 = -0.5 * t2 + 0.5 * t3; + + // b0 + b1 + b2 + b3 = 0 <=> b0 = -(b1 + b2 + b3) + // + // b0·P0 + b1·P1 + b2·P2 + b3·P3 + // = -(b1 + b2 + b3)·P0 + b1·P1 + b2·P2 + b3·P3 + // = b1·P1 + b2·P2 + b3·P3 - b1·P0 - b2·P0 - b3·P0 + // = b1·(P1 - P0) + b2·(P2 - P0) + b3·(P3 - P0) + + p1.sub(&p0) + .mul(b1) + .add(&p2.sub(&p0).mul(b2)) + .add(&p3.sub(&p0).mul(b3)) + } +} + +impl BSpline +where + T: Affine> + Clone, +{ + const _CHAR_MAT: Mat4 = { + const _1_6: f32 = 1.0 / 6.0; + const _2_3: f32 = 2.0 / 3.0; + mat![ + _1_6, _2_3, _1_6, 0.0; + -0.5, 0.0, 0.5, 0.0; + 0.5, -1.0, 0.5, 0.0; + -_1_6, 0.5, -0.5, _1_6; + ] + }; + + pub fn new(pts: impl IntoIterator) -> Self { + let pts: Vec<_> = pts.into_iter().collect(); + assert!(pts.len() >= 4, "a B-spline requires at least four points"); + Self(pts) + } + + /// Returns the point of `self` at the given *t* value. /// - /// # Examples - /// ``` - /// use retrofire_core::math::{BezierSpline, vec2, Vec2}; + /// Values of *t* outside the interval [0, 1] are accepted and extrapolate + /// the curve beyond the control points. + pub fn eval(&self, t: f32) -> T { + let (t, [p0, p1, p2, p3]) = crb_segment(&self.0, t); + let [_0, t1, t2, t3] = [1.0, t, t * t, t * t * t]; + + let _0 = (1.0 - 3.0 * t1 + 3.0 * t2 + t3) / 6.0; + let b1 = (4.0 - 6.0 * t2 + 3.0 * t3) / 6.0; + let b2 = (1.0 + 3.0 * t1 + 3.0 * t2 - 3.0 * t3) / 6.0; + let b3 = t3 / 6.0; + + p0.add(&p1.sub(&p0).mul(b1)) + .add(&p2.sub(&p0).mul(b2)) + .add(&p3.sub(&p0).mul(b3)) + } + + /// Returns the gradient, or velocity, vector of `self` at the given *t* + /// value. /// - /// let curve = BezierSpline::::new( - /// &[vec2(0.0, 0.0), vec2(0.0, 1.0), vec2(1.0, 1.0), vec2(1.0, 0.0)] - /// ); - /// let approx = curve.approximate_with(|err| err.len_sqr() < 0.01 * 0.01); - /// assert_eq!(approx.0.len(), 17); - /// ``` - pub fn approximate_with( - &self, - halt: impl Fn(&T::Diff) -> bool, - ) -> Polyline { - let len = self.0.len(); - let mut res = Vec::with_capacity(3 * len); - self.do_approx(0.0, 1.0, 10 + len.ilog2(), &halt, &mut res); - res.push(self.0[len - 1].clone()); - Polyline(res) - } - - fn do_approx( - &self, - a: f32, - b: f32, - max_dep: u32, - halt: &impl Fn(&T::Diff) -> bool, - accum: &mut Vec, - ) { - let mid = a.midpoint(b); - - let ap = self.eval(a); - let bp = self.eval(b); - - let real = self.eval(mid); - let approx = ap.midpoint(&bp); - - if max_dep == 0 || halt(&real.sub(&approx)) { - accum.push(ap); - } else { - self.do_approx(a, mid, max_dep - 1, halt, accum); - self.do_approx(mid, b, max_dep - 1, halt, accum); + /// Values of *t* outside the interval [0, 1] are accepted and extrapolate + /// the curve beyond the control points. + pub fn gradient(&self, t: f32) -> T::Diff { + let (t, [p0, p1, p2, p3]) = crb_segment(&self.0, t); + let [_0, _1, t2, t3] = [0.0, 1.0, 2.0 * t, 3.0 * t * t]; + + // ⎛ b0 ⎞ ⎛ -3 + 3·t2 - t3 ⎞ + // ⎜ b1 ⎟ = 1/6 ⎜ - 6·t2 + 3·t3 ⎟ + // ⎜ b2 ⎟ ⎜ 3 + 3·t2 - 3·t3 ⎟ + // ⎝ b3 ⎠ ⎝ t3 ⎠ + + let b1 = -t2 + 0.5 * t3; + let b2 = 0.5 + 0.5 * t2 - 0.5 * t3; + let b3 = t3 / 6.0; + + // b0 + b1 + b2 + b3 = 0 <=> b0 = -(b1 + b2 + b3) + // + // b0·P0 + b1·P1 + b2·P2 + b3·P3 + // = -(b1 + b2 + b3)·P0 + b1·P1 + b2·P2 + b3·P3 + // = b1·P1 + b2·P2 + b3·P3 - b1·P0 - b2·P0 - b3·P0 + // = b1·(P1 - P0) + b2·(P2 - P0) + b3·(P3 - P0) + + p1.sub(&p0) + .mul(b1) + .add(&p2.sub(&p0).mul(b2)) + .add(&p3.sub(&p0).mul(b3)) + } +} + +/// Returns the curve segment and local *t* value corresponding to +/// the global *t* value of a Catmull–Rom or B-spline. +fn crb_segment(pts: &[T], t: f32) -> (f32, &[T; 4]) { + let t = 1.0 + t * (pts.len() as f32 - 3.0); + + let i = (t as usize).clamp(1, pts.len() - 3); + let u = t - i as f32; + let pts = pts[i - 1..i + 3] // OK: 1 <= i < len - 3 + .try_into() + .expect("3 - (-1) = 4"); + (u, pts) +} + +impl Euclidean { + /// Creates a new `Euclidean` wrapper of the given spline. + pub fn new(spline: Spl) -> Self + where + Spl: Parametric>>, + { + let mut lut = Vec::new(); + let mut s = 0.0; + let mut p0 = spline.eval(0.0); + + // TODO smarter sampling + for t in 0.0.vary_to(1.0, 256) { + lut.push((s, t)); + let p1 = spline.eval(t); + s += p1.distance(&p0); + p0 = p1; + } + Self(spline, lut) + } + + /// Returns the point of `self` at distance *s* from the start, + /// as measured along the curve. + pub fn eval(&self, s: f32) -> T + where + Spl: Parametric, + { + self.0.eval(self.t(s)) + } + + /// Returns the approximate arc length of the spline. + pub fn len(&self) -> f32 { + self.1.last().map_or(0.0, |m| m.0) + } + + /// Returns the approximate *t* value for the given *s*. + fn t(&self, s: f32) -> f32 { + let lut = &self.1; + let i = lut.binary_search_by(|x| x.0.total_cmp(&s)); + match i { + Ok(i) => lut[i].1, + Err(0) => lut[0].1, + Err(i) if i == lut.len() => lut[lut.len() - 1].1, + Err(i) => { + let (s0, t0) = lut[i - 1]; // Ok: 0 < i + let (s1, t1) = lut[i]; // Ok: i < lut.len() + // Interpolate t in [t0, t1] given s in [s0, s1] + t0.lerp(&t1, inv_lerp(s, s0, s1)) + } } } } +// +// Local trait impls +// + impl Parametric for CubicBezier where T: Affine + Clone> + Clone, @@ -341,6 +785,15 @@ where } } +impl Parametric for CubicHermite +where + T: Affine + Clone> + Clone, +{ + fn eval(&self, t: f32) -> T { + self.eval(t) + } +} + impl Parametric for BezierSpline where T: Affine + Clone> + Clone, @@ -350,15 +803,51 @@ where } } +impl Parametric for HermiteSpline +where + T: Affine + Clone> + Clone, +{ + fn eval(&self, t: f32) -> T { + self.eval(t) + } +} + +impl Parametric for CatmullRomSpline +where + T: Affine + Clone> + Clone, +{ + fn eval(&self, t: f32) -> T { + self.eval(t) + } +} + +impl Parametric for BSpline +where + T: Affine + Clone> + Clone, +{ + fn eval(&self, t: f32) -> T { + self.eval(t) + } +} + +impl Parametric for Euclidean +where + Spl: Parametric, +{ + fn eval(&self, s: f32) -> T { + self.eval(s) + } +} + #[cfg(test)] mod tests { - use alloc::vec; - use crate::assert_approx_eq; use crate::math::{Parametric, Point2, Vec2, pt2, vec2}; use super::*; + const TEST_T_VALS: [f32; 7] = [-1.0, 0.0, 0.25, 0.5, 0.75, 1.0, 2.0]; + #[test] fn smoothstep_test() { assert_eq!(0.0, smoothstep(-10.0)); @@ -390,9 +879,9 @@ mod tests { } #[test] - fn bezier_spline_eval_eq_fast_eval() { - let b: CubicBezier = CubicBezier( - [[0.0, 0.0], [0.0, 2.0], [1.0, -1.0], [1.0, 1.0]].map(Vec2::from), + fn cubic_bezier_eval_eq_fast_eval() { + let b = CubicBezier( + [[0.0, 0.0], [0.0, 2.0], [1.0, -1.0], [1.0, 1.0]].map(::from), ); for i in 0..11 { let t = i as f32 / 10.0; @@ -403,88 +892,290 @@ mod tests { } #[test] - fn bezier_spline_eval_1d() { + fn cubic_bezier_f32_eval() { let b = CubicBezier([0.0, 2.0, -1.0, 1.0]); - assert_eq!(b.eval(-1.0), 0.0); - assert_eq!(b.eval(0.00), 0.0); - assert_eq!(b.eval(0.25), 0.71875); - assert_eq!(b.eval(0.50), 0.5); - assert_eq!(b.eval(0.75), 0.28125); - assert_eq!(b.eval(1.00), 1.0); - assert_eq!(b.eval(2.00), 1.0); + let expected = [-31.0, 0.0, 0.71875, 0.5, 0.28125, 1.0, 32.0]; + let actual = TEST_T_VALS.map(|t| b.eval(t)); + + assert_eq!(expected, actual); } #[test] - fn bezier_spline_eval_2d_vec() { - let b = CubicBezier( - [[0.0, 0.0], [0.0, 2.0], [1.0, -1.0], [1.0, 1.0]] - .map(Vec2::<()>::from), - ); + fn cubic_bezier_vec2_eval() { + let b = CubicBezier::([ + vec2(0.0, 0.0), + vec2(0.0, 2.0), + vec2(1.0, -1.0), + vec2(1.0, 1.0), + ]); - assert_eq!(b.eval(-1.0), vec2(0.0, 0.0)); - assert_eq!(b.eval(0.00), vec2(0.0, 0.0)); - assert_eq!(b.eval(0.25), vec2(0.15625, 0.71875)); - assert_eq!(b.eval(0.50), vec2(0.5, 0.5)); - assert_eq!(b.eval(0.75), vec2(0.84375, 0.281250)); - assert_eq!(b.eval(1.00), vec2(1.0, 1.0)); - assert_eq!(b.eval(2.00), vec2(1.0, 1.0)); + #[rustfmt::skip] + let expected = [ + [5.0, -31.0], [0.0, 0.0], [0.15625, 0.71875], [0.5, 0.5], + [0.84375, 0.281250], [1.0, 1.0], [-4.0, 32.0], + ]; + let actual = TEST_T_VALS.map(|t| b.eval(t).0); + + assert_eq!(expected, actual); } #[test] - fn bezier_spline_eval_2d_point() { - let b = CubicBezier( - [[0.0, 0.0], [0.0, 2.0], [1.0, -1.0], [1.0, 1.0]] - .map(Point2::<()>::from), - ); + fn cubic_bezier_point2_eval() { + let b = CubicBezier::([ + pt2(0.0, 0.0), + pt2(0.0, 2.0), + pt2(1.0, -1.0), + pt2(1.0, 1.0), + ]); - assert_eq!(b.eval(-1.0), pt2(0.0, 0.0)); - assert_eq!(b.eval(0.00), pt2(0.0, 0.0)); - assert_eq!(b.eval(0.25), pt2(0.15625, 0.71875)); - assert_eq!(b.eval(0.50), pt2(0.5, 0.5)); - assert_eq!(b.eval(0.75), pt2(0.84375, 0.281250)); - assert_eq!(b.eval(1.00), pt2(1.0, 1.0)); - assert_eq!(b.eval(2.00), pt2(1.0, 1.0)); + #[rustfmt::skip] + let expected = [ + [5.0, -31.0], [0.0, 0.0], [0.15625, 0.71875], [0.5, 0.5], + [0.84375, 0.281250], [1.0, 1.0], [-4.0, 32.0], + ]; + let actual = TEST_T_VALS.map(|t| b.eval(t).0); + + assert_eq!(expected, actual); } #[test] - fn bezier_spline_tangent_1d() { + fn cubic_bezier_f32_velocity() { let b = CubicBezier([0.0, 2.0, -1.0, 1.0]); - assert_eq!(b.tangent(-1.0), 6.0); - assert_eq!(b.tangent(0.00), 6.0); - assert_eq!(b.tangent(0.25), 0.375); - assert_eq!(b.tangent(0.50), -1.5); - assert_eq!(b.tangent(0.75), 0.375); - assert_eq!(b.tangent(1.00), 6.0); - assert_eq!(b.tangent(2.00), 6.0); + let expected = [66.0, 6.0, 0.375, -1.5, 0.375, 6.0, 66.0]; + let actual = TEST_T_VALS.map(|t| b.velocity(t)); + + assert_eq!(expected, actual); } #[test] - fn bezier_spline_tangent_2d() { - let b = CubicBezier( - [[0.0, 0.0], [0.0, 1.0], [1.0, 0.0], [1.0, 1.0]] - .map(Point2::<()>::from), + fn cubic_bezier_point2_velocity() { + #[rustfmt::skip] + let b = CubicBezier::([ + pt2(0.0, 0.0), pt2(0.0, 1.0), pt2(1.0, 0.0), pt2(1.0, 1.0), + ]); + + #[rustfmt::skip] + let expected = [ + [-12.0, 27.0], [0.0, 3.0], [1.125, 0.75], [1.5, 0.0], + [1.125, 0.75], [0.0, 3.0], [-12.0, 27.0], + ]; + let actual = TEST_T_VALS.map(|t| b.velocity(t).0); + + assert_eq!(expected, actual); + } + + #[test] + fn cubic_hermite_f32_eval() { + let h = CubicHermite([0.0, 0.0], [2.0, -2.0]); + + let expected = [-4.0, 0.0, 0.375, 0.5, 0.375, 0.0, -4.0]; + let actual = TEST_T_VALS.map(|t| h.eval(t)); + + assert_eq!(expected, actual); + } + + #[test] + fn cubic_hermite_f32_velocity() { + let h = CubicHermite([0.0, 0.0], [1.0, -1.0]); + + let expected = [3.0, 1.0, 0.5, 0.0, -0.5, -1.0, -3.0]; + let actual = TEST_T_VALS.map(|t| h.velocity(t)); + + assert_eq!(expected, actual); + } + + #[test] + fn bezier_hermite_equivalence() { + let [p0, p1, p2, p3] = + [pt2(0.0, 0.0), pt2(0.0, 1.0), pt2(1.0, 0.0), pt2(1.0, 1.0)]; + let b = CubicBezier::([p0, p1, p2, p3]); + let h = CubicHermite::( + [p0, p3], + [3.0 * (p1 - p0), 3.0 * (p3 - p2)], + ); + assert_eq!( + TEST_T_VALS.map(|t| b.eval(t).0), + TEST_T_VALS.map(|t| h.eval(t).0) + ); + } + + #[test] + fn bezier_spline_segment() { + let b = BezierSpline::new([0.0, 0.2, 0.4, 0.5, 0.6, 0.8, 1.0]); + + let first = CubicBezier([0.0, 0.2, 0.4, 0.5]); + let second = CubicBezier([0.5, 0.6, 0.8, 1.0]); + + assert_eq!(b.segment(-1.0), (-2.0, first)); + assert_eq!(b.segment(0.0), (0.0, first)); + assert_eq!(b.segment(0.5), (0.0, second)); + assert_eq!(b.segment(1.0), (1.0, second)); + assert_eq!(b.segment(2.0), (3.0, second)); + } + + #[test] + fn bezier_spline_f32_eval() { + let b = BezierSpline::new([0.0, 0.8, 0.9, 1.0, 0.6, 0.5, 0.5]); + + let expected = [-18.8, 0.0, 0.7625, 1.0, 0.6, 0.5, 0.1]; + let actual = TEST_T_VALS.map(|t| b.eval(t)); + + assert_approx_eq!(expected, actual); + } + + #[test] + fn bezier_spline_point2_from_rays() { + #[rustfmt::skip] + let expected = BezierSpline::::new([ + pt2(0.0, 0.0),pt2(0.0, 2.0),pt2(1.0, -1.0),pt2(1.0, 1.0) + ]); + let actual = BezierSpline::::from_rays([ + Ray(pt2(0.0, 0.0), vec2(0.0, 2.0)), + Ray(pt2(1.0, 1.0), vec2(0.0, 2.0)), + ]); + + assert_eq!( + TEST_T_VALS.map(|t| expected.eval(t)), + TEST_T_VALS.map(|t| actual.eval(t)) ); + } + + #[test] + fn bezier_spline_point2_velocity() { + #[rustfmt::skip] + let b = BezierSpline::::new([ + pt2(0.0, 0.0), pt2(0.0, 1.0), pt2(1.0, 0.0), pt2(1.0, 1.0), + ]); + #[rustfmt::skip] + let expected = [ + vec2(-12.0, 27.0), vec2(0.0, 3.0), vec2(1.125, 0.75), vec2(1.5, 0.0), + vec2(1.125, 0.75), vec2(0.0, 3.0), vec2(-12.0, 27.0), + ]; + let actual = TEST_T_VALS.map(|t| b.velocity(t)); + + assert_eq!(expected, actual); + } + + #[test] + fn hermite_spline_segment() { + let b = + HermiteSpline::new([Ray(0.0, 0.4), Ray(0.5, 0.1), Ray(1.0, 0.6)]); + + let first = CubicHermite([0.0, 0.5], [0.4, 0.1]); + let second = CubicHermite([0.5, 1.0], [0.1, 0.6]); + + assert_eq!(b.segment(-1.0), (-2.0, first)); + assert_eq!(b.segment(0.0), (0.0, first)); + assert_eq!(b.segment(0.5), (0.0, second)); + assert_eq!(b.segment(1.0), (1.0, second)); + assert_eq!(b.segment(2.0), (3.0, second)); + } + + #[test] + fn hermite_spline_point2_eval() { + let h = HermiteSpline::::new([ + Ray(pt2(0.0, 0.0), vec2(1.0, 0.0)), + Ray(pt2(1.0, 1.0), vec2(1.0, 0.0)), + ]); + + #[rustfmt::skip] + let expected = [ + [-1.0, 5.0], [0.0, 0.0], [0.25, 0.15625], [0.5, 0.5], + [0.75, 0.84375], [1.0, 1.0], [2.0, -4.0] + ]; + let actual = TEST_T_VALS.map(|t| h.eval(t).0); + + assert_eq!(expected, actual); + } + + #[test] + fn hermite_spline_point2_velocity() { + let h = HermiteSpline::::new([ + Ray(pt2(0.0, 0.0), vec2(1.0, 0.0)), + Ray(pt2(1.0, 1.0), vec2(1.0, 0.0)), + ]); - assert_eq!(b.tangent(-1.0), vec2(0.0, 3.0),); - assert_eq!(b.tangent(0.0), vec2(0.0, 3.0),); - assert_eq!(b.tangent(0.25), vec2(1.125, 0.75),); - assert_eq!(b.tangent(0.5), vec2(1.5, 0.0),); - assert_eq!(b.tangent(0.75), vec2(1.125, 0.75),); - assert_eq!(b.tangent(1.0), vec2(0.0, 3.0),); - assert_eq!(b.tangent(2.0), vec2(0.0, 3.0),); + #[rustfmt::skip] + let expected = [ + [1.0, -12.0], [1.0, 0.0], [1.0, 1.125], [1.0, 1.5], + [1.0, 1.125], [1.0, 0.0], [1.0, -12.0] + ]; + let actual = TEST_T_VALS.map(|t| h.velocity(t).0); + + assert_eq!(expected, actual); + } + + #[test] + fn catmull_rom_spline_point2_eval() { + #[rustfmt::skip] + let c = CatmullRomSpline::::new([ + pt2(-1.0, 0.0), pt2(0.0, 0.0), pt2(1.0, 0.0), + pt2(0.0, 1.0), pt2(1.0, 1.0), pt2(2.0, 1.0), + ]); + + #[rustfmt::skip] + let expected = [ + [33.0, -18.0], [0.0, 0.0], [0.890625, -0.0703125], [0.5, 0.5], + [0.109375, 1.0703125], [1.0, 1.0], [-32.0, 19.0] + ]; + let actual = TEST_T_VALS.map(|t| c.eval(t).0); + + assert_eq!(expected, actual); + } + + #[test] + fn catmull_rom_spline_point2_gradient() { + #[rustfmt::skip] + let c = CatmullRomSpline::::new([ + pt2(-1.0, 0.0), pt2(0.0, 0.0), pt2(1.0, 0.0), + pt2(0.0, 1.0), pt2(1.0, 1.0), pt2(2.0, 1.0), + ]); + + #[rustfmt::skip] + let expected = [ + [-32.0, 16.5], [1.0, 0.0], [0.8125, 0.09375], [-1.5, 1.25], + [0.8125, 0.09375], [1.0, 0.0], [-32.0, 16.5] + ]; + let actual = TEST_T_VALS.map(|t| c.gradient(t).0); + + assert_eq!(expected, actual); } #[test] - fn bezier_spline_eval() { - let c = BezierSpline(vec![0.0, 0.8, 0.9, 1.0, 0.6, 0.5, 0.5]); - assert_eq!(c.eval(-1.0), 0.0); - assert_eq!(c.eval(0.0), 0.0); - assert_approx_eq!(c.eval(0.25), 0.7625); - assert_eq!(c.eval(0.5), 1.0); - assert_eq!(c.eval(0.75), 0.6); - assert_eq!(c.eval(1.0), 0.5); - assert_eq!(c.eval(2.0), 0.5); + fn b_spline_point2_eval() { + #[rustfmt::skip] + let b = BSpline::::new([ + pt2(-1.0, 0.0), pt2(0.0, 0.0), pt2(1.0, 0.0), + pt2(0.0, 1.0), pt2(1.0, 1.0), pt2(2.0, 1.0), + ]); + + #[rustfmt::skip] + let expected = [ + [6.0, -4.5], [0.0, 0.0], [0.609375, 0.0703125], [0.5, 0.5], + [0.39062497, 0.92968756], [1.0, 1.0], [-5.0, 5.5] + ]; + let actual = TEST_T_VALS.map(|t| b.eval(t).0); + + assert_approx_eq!(expected, actual); + } + + #[test] + fn b_spline_point2_gradient() { + #[rustfmt::skip] + let b = BSpline::::new([ + pt2(-1.0, 0.0), pt2(0.0, 0.0), pt2(1.0, 0.0), + pt2(0.0, 1.0), pt2(1.0, 1.0), pt2(2.0, 1.0), + ]); + + #[rustfmt::skip] + let expected = [ + [-8.0, 4.5], [1.0, 0.0], [0.4375, 0.28125], [-0.5, 0.75], + [0.4375, 0.28125], [1.0, 0.0], [-8.0, 4.5] + ]; + let actual = TEST_T_VALS.map(|t| b.gradient(t).0); + + assert_eq!(expected, actual); } } diff --git a/core/src/math/vary.rs b/core/src/math/vary.rs index f23b5712..ec9fe73a 100644 --- a/core/src/math/vary.rs +++ b/core/src/math/vary.rs @@ -20,7 +20,7 @@ pub trait ZDiv: Sized { /// This trait is designed particularly for *varyings:* types that are /// meant to be interpolated across the face of a polygon when rendering, /// but the methods are useful for various purposes. -pub trait Vary: Lerp + ZDiv + Sized + Clone { +pub trait Vary: Lerp + ZDiv { /// The iterator returned by the [vary][Self::vary] method. type Iter: Iterator; /// The difference type of `Self`. diff --git a/core/src/math/vec.rs b/core/src/math/vec.rs index e5ab5717..dcc91f03 100644 --- a/core/src/math/vec.rs +++ b/core/src/math/vec.rs @@ -11,7 +11,7 @@ use core::{ ops::{AddAssign, DivAssign, MulAssign, SubAssign}, }; -use crate::math::{ +use super::{ Affine, ApproxEq, Linear, Point, space::{Proj3, Real}, vary::ZDiv, @@ -42,6 +42,7 @@ pub struct Vector(pub Repr, Pd); pub type Vec2 = Vector<[f32; 2], Real<2, Basis>>; /// A 3-vector with `f32` components. pub type Vec3 = Vector<[f32; 3], Real<3, Basis>>; + /// A `f32` 4-vector in the projective 3-space over ℝ, aka P3(ℝ). pub type ProjVec3 = Vector<[f32; 4], Proj3>; @@ -58,11 +59,13 @@ pub type Vec3i = Vector<[i32; 3], Real<3, Basis>>; // /// Returns a real 2-vector with components `x` and `y`. +#[inline] pub const fn vec2(x: Sc, y: Sc) -> Vector<[Sc; 2], Real<2, B>> { Vector([x, y], Pd) } /// Returns a real 3-vector with components `x`, `y`, and `z`. +#[inline] pub const fn vec3(x: Sc, y: Sc, z: Sc) -> Vector<[Sc; 3], Real<3, B>> { Vector([x, y, z], Pd) } @@ -127,7 +130,6 @@ impl Vector { // TODO Many of these functions could be more generic impl Vector<[f32; N], Sp> { /// Returns the length (magnitude) of `self`. - #[cfg(feature = "fp")] #[inline] pub fn len(&self) -> f32 { super::float::f32::sqrt(self.dot(self)) @@ -146,7 +148,7 @@ impl Vector<[f32; N], Sp> { /// ``` /// /// # Panics - /// Panics in dev mode if `self` is a zero vector. + /// Panics if the length of `self` is approximately zero. #[inline] #[must_use] pub fn normalize(&self) -> Self { @@ -162,6 +164,34 @@ impl Vector<[f32; N], Sp> { *self * f32::recip_sqrt(len_sqr) } + /// Returns `self` normalized to unit length, or a zero vector if the + /// length of `self` is approximately zero. + /// + /// # Examples + /// ``` + /// use retrofire_core::assert_approx_eq; + /// use retrofire_core::math::{vec2, Vec2}; + /// + /// let normalized: Vec2 = vec2(3.0, 4.0).normalize_or_zero(); + /// assert_approx_eq!(normalized, vec2(0.6, 0.8), eps=1e-2); + /// + /// let zero: Vec2 = vec2(0.0, 0.0).normalize_or_zero(); + /// assert_eq!(zero, vec2(0.0, 0.0)); + /// ``` + #[inline] + #[must_use] + pub fn normalize_or_zero(&self) -> Self { + #[cfg(feature = "std")] + use super::float::RecipSqrt; + use super::float::f32; + let len_sqr = self.len_sqr(); + if len_sqr.approx_eq_eps(&0.0, &1e-12) { + Vector::zero() + } else { + *self * f32::recip_sqrt(len_sqr) + } + } + /// Returns `self` clamped component-wise to the given range. /// /// In other words, for each component `self[i]`, the result `r` has @@ -184,8 +214,7 @@ impl Vector<[f32; N], Sp> { array::from_fn(|i| self[i].clamp(min[i], max[i])).into() } - /// Returns `true` if every component of `self` is finite, - /// `false` otherwise. + /// Returns `true` if every component of `self` is finite, `false` otherwise. /// /// See [`f32::is_finite()`]. pub fn is_finite(&self) -> bool { @@ -200,7 +229,7 @@ where { /// Returns the length of `self`, squared. /// - /// This avoids taking the square root in cases it's not needed, + /// This avoids taking the square root in cases where it's not needed, /// and works with scalars for which a square root is not defined. #[inline] pub fn len_sqr(&self) -> Sc { @@ -208,6 +237,8 @@ where } /// Returns the dot product of `self` and `other`. + /// + /// TODO docs #[inline] pub fn dot(&self, other: &Self) -> Sc { zip(&self.0, &other.0) @@ -331,8 +362,8 @@ impl Vector<[Sc; N], Sp> { /// ``` #[inline] #[must_use] - pub fn map(self, mut f: impl FnMut(Sc) -> T) -> Vector<[T; N], Sp> { - array::from_fn(|i| f(self.0[i])).into() + pub fn map(self, f: impl FnMut(Sc) -> T) -> Vector<[T; N], Sp> { + self.0.map(f).into() } /// Returns a vector of the same dimension as `self` by applying `f` /// component-wise to `self` and `other`. @@ -355,6 +386,70 @@ impl Vector<[Sc; N], Sp> { ) -> Vector<[U; N], Sp> { array::from_fn(|i| f(self.0[i], other.0[i])).into() } + + #[inline] + pub fn max(&self) -> Sc + where + Sc: PartialOrd, + { + const { assert!(N > 0, "0D vectors have no maximum") } + let mut max = self.0[0]; + for c in self.0 { + if c > max { + max = c; + } + } + max + } + + #[inline] + pub fn min(&self) -> Sc + where + Sc: PartialOrd, + { + const { assert!(N > 0, "0D vectors have no minimum") } + let mut min = self.0[0]; + for c in self.0 { + if c < min { + min = c; + } + } + min + } + + #[inline] + pub fn argmax(&self) -> usize + where + Sc: PartialOrd, + { + const { assert!(N > 0, "0D vectors have no maximum") } + let mut max = self.0[0]; + let mut max_i = 0; + for (c, i) in zip(self.0, 0..) { + if c > max { + max = c; + max_i = i; + } + } + max_i + } + + #[inline] + pub fn argmin(&self) -> usize + where + Sc: PartialOrd, + { + const { assert!(N > 0, "0D vectors have no minimum") } + let mut min = self.0[0]; + let mut min_i = 0; + for (c, i) in zip(self.0, 0..) { + if c < min { + min = c; + min_i = i; + } + } + min_i + } } impl Vector<[Sc; 2], Real<2, B>> { @@ -406,17 +501,17 @@ impl Vec2 { /// /// where *θ* is the (signed) angle between **a** and **b**. In particular, /// the result is zero if **a** and **b** are parallel (or either is zero), - /// positive if the angle from **a** to **b** is positive, and negative if - /// the angle is negative: + /// negative if the angle from **a** to **b** is negative (clockwise), and + /// positive if the angle is positive (counter-clockwise) /// /// ```text - /// ^ b ^ a - /// / ^ b / ^ a - /// ^ a \ / \ - /// / \ / \ - /// O O O-----> b + /// ^ b + /// / a ^ ^ a + /// ^ a / \ + /// / / \ + /// O b <-----O O-----> b /// - /// a⟂·b = 0 a⟂·b > 0 a⟂·b < 0 + /// a⟂·b = 0 a⟂·b > 0 a⟂·b < 0 /// ``` /// /// # Examples @@ -433,6 +528,24 @@ impl Vec2 { pub fn perp_dot(self, other: Self) -> f32 { self.perp().dot(&other) } + + /// Returns the angle between `self` and the positive x-axis. + /// + /// Equivalent to `atan2(self.y(), self.x())` or `self.to_polar().az()`. + /// + /// # Examples + /// ``` + /// use retrofire_core::{assert_approx_eq, math::{vec2, degs}}; + /// let vec2 = vec2::; + /// + /// assert_approx_eq!(vec2(1.0, 1.0).atan(), degs(45.0)); + /// assert_approx_eq!(vec2(-2.0, 0.0).atan(), degs(180.0)); + /// assert_approx_eq!(vec2(0.0, -3.0).atan(), degs(-90.0)); + /// ``` + #[cfg(feature = "fp")] + pub fn atan(self) -> super::Angle { + super::atan2(self.y(), self.x()) + } } impl Vector<[Sc; 3], Real<3, B>> @@ -461,15 +574,13 @@ where /// proportional to the area of the parallelogram formed by the vectors. /// Specifically, the length is given by the identity: /// - /// ```text - /// |𝗮 × 𝗯| = |𝗮| |𝗯| sin 𝜽 - /// ``` + /// |**a** × **b**| = |**a**| |**b**| sin *θ*, /// /// where |·| denotes the length of a vector and 𝜽 equals the angle - /// between 𝗮 and 𝗯. Specifically, the result has unit length if 𝗮 and 𝗯 - /// are orthogonal and |𝗮| = |𝗯| = 1. The cross product can be used to - /// produce an *orthonormal basis* from any two non-parallel non-zero - /// 3-vectors. + /// between **a** and **b**. Specifically, the result has unit length if + /// **a** and **b** are orthogonal and |**a**| = |**b**| = 1. The cross + /// product can be used to produce an *orthonormal basis* from any two + /// non-parallel non-zero 3-vectors. /// /// ```text /// ^ @@ -609,10 +720,11 @@ impl Copy for Vector {} impl Clone for Vector { fn clone(&self) -> Self { - Self(self.0.clone(), Pd) + Self::new(self.0.clone()) } } +// Limited to Cartesian vectors because Spherical/PolarVec have own impl impl Default for Vector> { fn default() -> Self { Self::new(R::default()) @@ -678,7 +790,7 @@ where /// Note that `Self::DIM` can be less than the number of elements in `R`. #[inline] fn index(&self, i: usize) -> &Self::Output { - assert!(i < Self::DIM, "index {i} out of bounds ({})", Self::DIM); + debug_assert!(i < Self::DIM, "index {i} out of bounds ({})", Self::DIM); &self.0[i] } } @@ -695,7 +807,7 @@ where /// Note that `Self::DIM` can be less than the number of elements in `R`. #[inline] fn index_mut(&mut self, i: usize) -> &mut Self::Output { - assert!(i < Self::DIM, "index {i} out of bounds ({})", Self::DIM); + debug_assert!(i < Self::DIM, "index {i} out of bounds ({})", Self::DIM); &mut self.0[i] } } @@ -865,7 +977,6 @@ mod tests { mod f32 { use super::*; - #[cfg(feature = "fp")] #[test] fn length() { assert_approx_eq!(vec2(1.0, 1.0).len(), SQRT_2); @@ -992,8 +1103,8 @@ mod tests { #[test] fn perp() { - assert_eq!(Vec2::<()>::zero().perp(), Vec2::zero()); - assert_eq!(Vec2::<()>::X.perp(), Vec2::Y); + assert_eq!(::zero().perp(), ::zero()); + assert_eq!(::X.perp(), Vec2::Y); assert_eq!(vec2(-0.2, -1.5).perp(), vec2(1.5, -0.2)); } diff --git a/core/src/render.rs b/core/src/render.rs index b61b61ce..f6d8c6f0 100644 --- a/core/src/render.rs +++ b/core/src/render.rs @@ -20,23 +20,27 @@ use self::{ raster::Scanline, }; -pub use self::{ - batch::Batch, - cam::Camera, - clip::Clip, - ctx::Context, - raster::Frag, - shader::{FragmentShader, VertexShader}, - stats::Stats, - target::{Colorbuf, Framebuf, Target}, - tex::{TexCoord, Texture, uv}, - text::Text, -}; +pub(super) mod re_exports { + pub use super::{ + batch::Batch, + cam::Camera, + clip::Clip, + ctx::Context, + raster::Frag, + shader::{FragmentShader, VertexShader}, + stats::Stats, + target::{Colorbuf, Framebuf, Target}, + tex::{TexCoord, Texture, uv}, + text::Text, + }; +} +pub use re_exports::*; pub mod batch; pub mod cam; pub mod clip; pub mod ctx; +pub mod debug; pub mod prim; pub mod raster; pub mod scene; diff --git a/core/src/render/batch.rs b/core/src/render/batch.rs index 287e3abb..bf917455 100644 --- a/core/src/render/batch.rs +++ b/core/src/render/batch.rs @@ -3,10 +3,8 @@ use alloc::vec::Vec; use core::borrow::Borrow; -use crate::{ - geom::{Mesh, Tri, Vertex3}, - math::{Mat4, Vary}, -}; +use crate::geom::{Edge, Mesh, Tri, Vertex3}; +use crate::math::{Mat4, Vary}; use super::{Clip, Context, Ndc, Render, Screen, Shader, Target}; @@ -31,13 +29,13 @@ use super::{Clip, Context, Ndc, Render, Screen, Shader, Target}; // [instances]: https://en.wikipedia.org/wiki/Geometry_instancing #[derive(Clone, Debug, Default)] pub struct Batch { - prims: Vec, - verts: Vec, - uniform: Uni, - shader: Shd, - viewport: Mat4, - target: Tgt, - ctx: Ctx, + pub prims: Vec, + pub verts: Vec, + pub uniform: Uni, + pub shader: Shd, + pub viewport: Mat4, + pub target: Tgt, + pub ctx: Ctx, } macro_rules! update { @@ -147,3 +145,25 @@ impl Batch { ); } } + +impl Batch, Vtx, Uni, Shd, Tgt, Ctx> { + pub fn append(&mut self, other: Self) { + let Batch { prims, verts, .. } = other; + let n = self.verts.len(); + let prims = prims.into_iter().map(|e| Edge(e.0 + n, e.1 + n)); + + self.verts.extend(verts); + self.prims.extend(prims) + } +} + +impl Batch, Vtx, Uni, Shd, Tgt, Ctx> { + pub fn append(&mut self, other: Self) { + let Batch { prims, verts, .. } = other; + let n = self.verts.len(); + let prims = prims.into_iter().map(|tri| tri.map(|i| i + n)); + + self.verts.extend(verts); + self.prims.extend(prims); + } +} diff --git a/core/src/render/cam.rs b/core/src/render/cam.rs index 61ae5e14..0b9e38d5 100644 --- a/core/src/render/cam.rs +++ b/core/src/render/cam.rs @@ -4,11 +4,11 @@ use core::ops::Range; #[cfg(feature = "fp")] use crate::math::{ - Angle, Vec3, orient_z, rotate_x, rotate_y, spherical, translate, turns, + Angle, Vec3, orient_z, rotate_pyr, rotate_x, rotate_y, spherical, turns, }; use crate::math::{ - Lerp, Mat4, Point3, ProjMat3, SphericalVec, Vary, mat::RealToReal, - orthographic, perspective, pt2, viewport, + Mat4, Point3, ProjMat3, SphericalVec, Vary, orthographic, perspective, pt2, + translate, viewport, }; use crate::util::{Dims, rect::Rect}; @@ -67,7 +67,7 @@ pub struct Camera { /// /// This is the familiar "FPS" movement mode, based on camera /// position and heading (look-at vector). -#[derive(Copy, Clone, Debug)] +#[derive(Copy, Clone, Debug, Default)] pub struct FirstPerson { /// Current position of the camera in **world** space. pub pos: Point3, @@ -75,19 +75,12 @@ pub struct FirstPerson { pub heading: SphericalVec, } -pub type ViewToWorld = RealToReal<3, View, World>; - -/// Creates a unit `SphericalVec` from azimuth and altitude. -#[cfg(feature = "fp")] -fn az_alt(az: Angle, alt: Angle) -> SphericalVec { - spherical(1.0, az, alt) -} /// Orbiting camera transform. /// /// Keeps the camera centered on a **world-space** point, and allows free /// 360°/180° azimuth/altitude rotation around that point as well as setting /// the distance from the point. -#[derive(Copy, Clone, Debug)] +#[derive(Copy, Clone, Debug, Default)] pub struct Orbit { /// The camera's target point in **world** space. pub target: Point3, @@ -95,6 +88,24 @@ pub struct Orbit { pub dir: SphericalVec, } +/// Camera transform implementing airplane-like controls based on three angles. +/// +/// The pitch (elevation) angle controls rotation about the local lateral (x) +/// axis, the yaw (bearing) angle about the local vertical (y) axis, and the +/// roll (bank) angle about the local longitudinal (z) axis. These angles are +/// also known to as Tait–Bryan angles. +#[derive(Copy, Clone, Debug, Default)] +pub struct PitchYawRoll { + pub rot: Mat4, + pub pos: Point3, +} + +/// Creates a unit `SphericalVec` from azimuth and altitude. +#[cfg(feature = "fp")] +fn az_alt(az: Angle, alt: Angle) -> SphericalVec { + spherical(1.0, az, alt) +} + // // Inherent impls // @@ -195,13 +206,22 @@ impl Camera { } impl Camera { + /// Returns the camera matrix. + pub fn world_to_view(&self) -> Mat4 { + self.transform.world_to_view() + } + /// Returns the inverse camera matrix. + pub fn view_to_world(&self) -> Mat4 { + self.world_to_view().inverse() + } + /// Returns the composed camera and projection matrix. pub fn world_to_project(&self) -> ProjMat3 { - self.transform.world_to_view().then(&self.project) + self.world_to_view().then(&self.project) } /// Renders the given geometry from the viewpoint of this camera. - pub fn render( + pub fn render( &self, prims: impl AsRef<[Prim]>, verts: impl AsRef<[Vtx]>, @@ -234,10 +254,7 @@ impl FirstPerson { /// Creates a first-person transform with position in the origin /// and heading in the direction of the positive x-axis. pub fn new() -> Self { - Self { - pos: Point3::origin(), - heading: az_alt(turns(0.0), turns(0.0)), - } + Self::default() } /// Rotates the camera to center the view on a **world-space** point. @@ -276,6 +293,11 @@ impl FirstPerson { #[cfg(feature = "fp")] impl Orbit { + /// TODO + pub fn new() -> Self { + Self::default() + } + /// Adds the azimuth and altitude to the camera's current direction. /// /// Wraps the resulting azimuth to [-180°, 180°) and clamps the altitude to [-90°, 90°]. @@ -323,10 +345,55 @@ impl Orbit { } } +#[cfg(feature = "fp")] +impl PitchYawRoll { + /// Creates a new pitch/yaw/roll transform, with the initial position + /// at the world space origin and facing along the negative z-axis. + pub fn new() -> Self { + Self::default() + } + + /// Adjusts the camera position in view space (relative to the current + /// position, along the current orientation axes). + pub fn translate(&mut self, v: Vec3) { + self.pos += self.view_to_world().apply(&v); + } + + /// Moves the camera to the given position in world space. + pub fn translate_to(&mut self, pt: Point3) { + self.pos = pt; + } + + /// Adjusts the orientation of the camera by the given delta angles. + pub fn rotate(&mut self, pitch: Angle, yaw: Angle, roll: Angle) { + self.rot = rotate_pyr(pitch, yaw, roll).to().then(&self.rot); + } + + /// Sets the orientation of the camera to the given **world-space** angles. + /// + /// To adjust the camera in view space (that is, relative to the current + /// orientation), use [`rotate()`][Self::rotate]. + pub fn rotate_to(&mut self, pitch: Angle, yaw: Angle, roll: Angle) { + self.rot = rotate_pyr(pitch, yaw, roll).to(); + } + + /// Returns the matrix from view to world space. + pub fn view_to_world(&self) -> Mat4 { + let trans: Mat4 = translate(self.pos.to_vec().to()).to(); + self.rot.then(&trans) + } +} + // // Local trait impls // +impl Transform for Mat4 { + fn world_to_view(&self) -> Mat4 { + *self + } +} + #[cfg(feature = "fp")] impl Transform for FirstPerson { fn world_to_view(&self) -> Mat4 { @@ -360,40 +427,25 @@ impl Transform for Orbit { } } -impl Transform for Mat4 { +impl Transform for PitchYawRoll { fn world_to_view(&self) -> Mat4 { - *self - } -} + let trans = translate(-self.pos.to_vec().to()).to(); + let rot = self.rot.transpose(); -// -// Foreign trait impls -// - -#[cfg(feature = "fp")] -impl Default for FirstPerson { - /// Returns [`FirstPerson::new`]. - fn default() -> Self { - Self::new() - } -} - -#[cfg(feature = "fp")] -impl Default for Orbit { - fn default() -> Self { - Self { - target: Point3::default(), - dir: az_alt(turns(0.0), turns(0.0)), - } + trans.then(&rot) } } #[cfg(test)] mod tests { - use crate::assert_approx_eq; - use super::*; + #[cfg(feature = "fp")] + use crate::{ + assert_approx_eq, + math::{SQRT_3, degs}, + }; + use Fov::*; #[test] @@ -413,7 +465,6 @@ mod tests { #[cfg(feature = "fp")] #[test] fn angle_of_view_focal_ratio_with_unit_aspect_ratio() { - use crate::math::{SQRT_3, degs}; use core::f32::consts::SQRT_2; assert_approx_eq!(Horizontal(degs(60.0)).focal_ratio(1.0), SQRT_3); @@ -427,8 +478,6 @@ mod tests { #[cfg(feature = "fp")] #[test] fn angle_of_view_focal_ratio_with_other_aspect_ratio() { - use crate::math::{SQRT_3, degs}; - assert_approx_eq!(Horizontal(degs(60.0)).focal_ratio(SQRT_3), SQRT_3); assert_approx_eq!(Vertical(degs(60.0)).focal_ratio(SQRT_3), 1.0); assert_approx_eq!(Diagonal(degs(60.0)).focal_ratio(SQRT_3), 2.0); diff --git a/core/src/render/clip.rs b/core/src/render/clip.rs index 040c0008..2b64f3ff 100644 --- a/core/src/render/clip.rs +++ b/core/src/render/clip.rs @@ -16,7 +16,6 @@ use alloc::vec::Vec; use core::{iter::zip, mem::swap}; - use crate::geom::{Edge, Tri, Vertex, vertex}; use crate::math::{Lerp, ProjVec3}; @@ -165,7 +164,7 @@ impl ClipPlane { /// `---__ \ /// `---B /// ``` - pub fn clip_simple_polygon( + pub fn clip_simple_polygon( &self, verts_in: &[ClipVert], verts_out: &mut Vec>, @@ -281,7 +280,7 @@ pub mod view_frustum { /// /// [^1]: Ivan Sutherland, Gary W. Hodgman: Reentrant Polygon Clipping. /// Communications of the ACM, vol. 17, pp. 32–42, 1974 -pub fn clip_simple_polygon<'a, A: Lerp + Clone>( +pub fn clip_simple_polygon<'a, A: Lerp>( planes: &[ClipPlane], verts_in: &'a mut Vec>, verts_out: &'a mut Vec>, @@ -309,7 +308,7 @@ impl ClipVert { } } -impl Clip for [Edge>] { +impl Clip for [Edge>] { type Item = Edge>; fn clip(&self, planes: &[ClipPlane], out: &mut Vec) { @@ -348,7 +347,7 @@ impl Clip for [Edge>] { } } -impl Clip for [Tri>] { +impl Clip for [Tri>] { type Item = Tri>; fn clip(&self, planes: &[ClipPlane], out: &mut Vec) { @@ -419,7 +418,7 @@ mod tests { } fn tri(a: ClipVec, b: ClipVec, c: ClipVec) -> Tri> { - Tri([a, b, c].map(vtx)) + Tri([a, b, c]).map(vtx) } #[test] diff --git a/core/src/render/debug.rs b/core/src/render/debug.rs new file mode 100644 index 00000000..b1b6646a --- /dev/null +++ b/core/src/render/debug.rs @@ -0,0 +1,195 @@ +//! Routines for drawing wireframe visualizations of geometric objects +//! for debugging purposes. Includes normals, bounding boxes, and more. + +#[cfg(feature = "fp")] +use alloc::vec::Vec; +use core::fmt::Debug; + +use crate::geom::{Edge, Tri, Vertex, Vertex3, vertex}; +use crate::math::{ + Color, Color4, Color4f, Mat4, Point3, Vec3, color::gray, mat::ProjMat3, + pt3, vec::ProjVec3, +}; +#[cfg(feature = "fp")] +use crate::math::{Vary, polar, turns, vec3}; + +use super::{Context, Frag, FragmentShader, VertexShader, scene::BBox}; + +#[derive(Default)] +pub struct Shader; + +impl<'a, B> VertexShader, &'a ProjMat3> for Shader { + type Output = Vertex; + + fn shade_vertex( + &self, + v: Vertex3, + m: &'a ProjMat3, + ) -> Self::Output { + vertex(m.apply(&v.pos), v.attrib) + } +} + +impl FragmentShader for Shader { + fn shade_fragment(&self, f: Frag) -> Option { + Some(f.var.to_color4()) + } +} + +pub type DbgBatch = + super::Batch, Vertex3, (), Shader, (), Context>; + +/// Returns a color visualizing the direction of a vector. +/// +/// # Examples +/// ``` +/// use retrofire_core::math::{Vec3, rgba, vec3}; +/// use retrofire_core::render::debug::dir_to_rgb; +/// +/// let right: Vec3 = vec3(1.0, 0.0, 0.0); +/// assert_eq!(dir_to_rgb(right), rgba(1.0, 0.5, 0.5, 1.0)); +/// +/// let down: Vec3 = vec3(0.0, -1.0, 0.0); +/// assert_eq!(dir_to_rgb(down), rgba(0.5, 0.0, 0.5, 1.0)); +/// +/// ``` +pub fn dir_to_rgb(v: Vec3) -> Color4f { + (0.5 * Color::new(v.0) + gray(0.5)) + .clamp(&gray(0.0), &gray(1.0)) + .to_rgba() +} + +/// Draws an illustration of a ray. +pub fn ray<'a, B>(o: Point3, dir: Vec3) -> DbgBatch { + let mut b = dir.cross(&Vec3::Y); + if b.len_sqr() < 1e-6 { + b = dir.cross(&Vec3::X); + } + let b = b.normalize_or_zero(); + let c = dir.cross(&b).normalize_or_zero(); + + let (head_w, head_h) = (0.04, 0.1); + let a = o + dir - head_h * dir.normalize_or_zero(); + let b = head_w * b; + let c = head_w * c; + + let verts = [o, o + dir, a + b, a - b, a + c, a - c] + .map(|p| vertex(p, dir_to_rgb(dir))); + #[rustfmt::skip] + let edges = [ + [0, 1], [1, 2], [1, 3], [1, 4], [1, 5], + [2, 4], [2, 5], [3, 4], [3, 5], + ].map(|[i, j]| Edge(i, j)); + + DbgBatch::new(&edges, &verts) +} + +/// Draws a unit-length ray denoting the normal vector of a triangle. +/// +/// The ray originates from the triangle's centroid. +pub fn face_normal( + tri: Tri>, +) -> DbgBatch { + ray(tri.centroid(), tri.normal().to()) +} + +/// Draws a visualization of an affine basis. +/// +/// Draws three rays representing the coordinate axes. The rays originate +/// from the origin point of the basis. +pub fn basis(m: Mat4) -> DbgBatch { + let xyz = m.linear(); + let x = xyz.col_vec(0); + let y = xyz.col_vec(1); + let z = xyz.col_vec(2); + let o = m.origin(); + + let mut b = ray(o, x); + b.append(ray(o, y)); + b.append(ray(o, z)); + b +} + +/// Draws an axis-aligned box with the given opposite vertices. +pub fn cuboid(v0: Point3, v1: Point3) -> DbgBatch { + let [x0, y0, z0] = v0.0; + let [x1, y1, z1] = v1.0; + #[rustfmt::skip] + let verts = [ + v0, pt3(x0, y0, z1), + pt3(x0, y1, z1), pt3(x0, y1, z0), + pt3(x1, y0, z0), pt3(x1, y0, z1), + v1, pt3(x1, y1, z0), + ].map(|p| vertex(p, dir_to_rgb(p.to_vec()))); + #[rustfmt::skip] + let edges = [ + [0, 1], [1, 2], [2, 3], [3, 0], + [4, 5], [5, 6], [6, 7], [7, 4], + [0, 4], [1, 5], [2, 6], [3, 7], + ].map(|[i, j]| Edge(i, j)); + + DbgBatch::new(&edges, &verts) +} + +/// Draws the smallest axis-aligned box that contains a set of vertices. +pub fn bbox(vs: &[Vertex3]) -> DbgBatch { + let BBox(min, max) = vs.iter().map(|v| &v.pos).collect(); + cuboid(min, max) +} + +/// Draws a circle on the XY plane with the given center and radius. +#[cfg(feature = "fp")] +pub fn circle(o: Point3, r: f32) -> DbgBatch { + const RES: usize = 64; // TODO constant, use array rather than Vec + + let verts: Vec<_> = 0.0 + .vary_to(1.0, RES as u32 + 1) + .map(|a| { + let v = polar(r, turns(a)).to_cart().to_vec3(); + vertex(o + v, dir_to_rgb(v)) + }) + .collect(); + + let edges: Vec<_> = (0..RES).map(|i| Edge(i, i + 1)).collect(); + + DbgBatch::new(&edges, &verts) +} + +/// Draws a wireframe sphere with the given center and radius. +/// +/// The sphere is represented by three circles lying on the XY, XZ, and +/// YZ planes. +#[cfg(feature = "fp")] +pub fn sphere(o: Point3, r: f32) -> DbgBatch { + const RES: usize = 64; + + let verts: Vec<_> = 0.0 + .vary_to(1.0, RES as u32 + 1) + .flat_map(|a| { + let [x, y] = polar::<()>(r, turns(a)).to_cart().0; + [vec3(x, y, 0.0), vec3(x, 0.0, y), vec3(0.0, x, y)] + .map(|v| vertex(o + v, dir_to_rgb(v))) + }) + .collect(); + + let edges: Vec<_> = (0..RES) + .flat_map(|i| { + [ + Edge(3 * i, 3 * i + 3), + Edge(3 * i + 1, 3 * i + 4), + Edge(3 * i + 2, 3 * i + 5), + ] + }) + .collect(); + + DbgBatch::new(&edges, &verts) +} + +impl DbgBatch { + fn new(prims: &[Edge], verts: &[Vertex3]) -> Self { + DbgBatch::::default() + .primitives(prims) + .vertices(verts) + .shader(Shader) + } +} diff --git a/core/src/render/prim.rs b/core/src/render/prim.rs index 609cb85d..a510be9a 100644 --- a/core/src/render/prim.rs +++ b/core/src/render/prim.rs @@ -14,11 +14,8 @@ impl Render for Tri { type Clips = [Tri>]; type Screen = Tri>; - fn inline( - Tri([i, j, k]): Tri, - vs: &[ClipVert], - ) -> Tri> { - Tri([vs[i].clone(), vs[j].clone(), vs[k].clone()]) + fn inline(tri: Self, vs: &[ClipVert]) -> Tri> { + tri.map(|i| vs[i].clone()) } fn depth(Tri([a, b, c]): &Self::Clip) -> f32 { diff --git a/core/tests/rendering.rs b/core/tests/rendering.rs index 953ac5ad..314537d7 100644 --- a/core/tests/rendering.rs +++ b/core/tests/rendering.rs @@ -3,7 +3,7 @@ use retrofire_core::prelude::*; use retrofire_core::{ - render::tex::SamplerClamp, + render::{Model, render, shader, tex::SamplerClamp}, util::{self, pixfmt::Xrgb8888, pnm::parse_pnm}, }; @@ -44,8 +44,8 @@ fn textured_quad() { assert_eq!(framebuf[255][0], rgb(0x7F, 0, 0)); assert_eq!(framebuf[0][255], rgb(0x7F, 0, 0)); - let comp = *include_bytes!("textured_quad.ppm"); - let comp = parse_pnm(comp).expect("should be a valid ppm"); + static COMP: &[u8] = include_bytes!("textured_quad.ppm"); + let comp = parse_pnm(COMP.iter().copied()).expect("should be a valid ppm"); assert_eq!(framebuf, comp); diff --git a/core/triangle.ppm b/core/triangle.ppm new file mode 100644 index 00000000..a0887752 Binary files /dev/null and b/core/triangle.ppm differ diff --git a/demos/Cargo.toml b/demos/Cargo.toml index 4960e7c7..4e6e3ac5 100644 --- a/demos/Cargo.toml +++ b/demos/Cargo.toml @@ -19,6 +19,7 @@ license.workspace = true keywords.workspace = true categories.workspace = true repository.workspace = true +documentation.workspace = true [dependencies] re = { version = "0.4.0", path = "..", package = "retrofire" } diff --git a/demos/nostd/Cargo.toml b/demos/nostd/Cargo.toml new file mode 100644 index 00000000..67bebf0f --- /dev/null +++ b/demos/nostd/Cargo.toml @@ -0,0 +1,27 @@ +[package] +name = "retrofire-no-std-demo" +edition = "2024" +version = "0.1.0" +license = "MIT" +authors = ["Johannes 'Sharlin' Dahlström"] +repository = "https://github.com/jdahlstrom/retrofire" + +[workspace] + +[dependencies.libc] +version = "0.2.177" +default-features = false + +[dependencies.re] +package = "retrofire-core" +version = "0.4.0-pre4" +path = "../../core" +default-features = false + +[profile.release] +opt-level = 3 +codegen-units = 1 +lto = "fat" +debug = 0 +panic = "abort" +strip = true diff --git a/demos/nostd/LICENSE-MIT b/demos/nostd/LICENSE-MIT new file mode 100644 index 00000000..5e178dc2 --- /dev/null +++ b/demos/nostd/LICENSE-MIT @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2025 Johannes Dahlström + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/demos/nostd/README.md b/demos/nostd/README.md new file mode 100644 index 00000000..99e0a6b7 --- /dev/null +++ b/demos/nostd/README.md @@ -0,0 +1,60 @@ +``` + ______ + ___ /´ ___/\ + __ ______ _____ / /\_ _ ______ _____ ____/ /_/___/\ __ _____ ______ + ==/ ´ ____/ __ \ ____/ ´ ____/ __ ` __ ___, /==/ ´ ___/ __ \ + ==/ /´=/ ______/ /==/ /´=/ /==/ /=/ /=/ /==/ /´=/ ______/\ + ==/ /==/ /____/ /__/ /==/ /__/ /=/ /=/ /__/ /==/ /______\/ +==/___/ ==\_______/\______/__/ ==\________,´_/ /==\______/__/ ==\________/\ +==\___\/ ==\______\/\_____\__\/ ==\______/_____,´ /==\_____\___\/==\_______\/ + \_____\,´ +``` + +## Retrofire no_std demo + +A small demo program aimed at minimizing the executable size. It uses +[`retrofire`][rf] to transform, clip, project, and rasterize a colorful +triangle and prints the result in PPM format to stdout. Pipe or redirect +the output if you don't want to fill your terminal with binary nonsense: + +``` + $ cargo run -r > output.ppm +``` + +Unless `std` is compiled without panic support, the build is somewhat finicky +and requires the use of the release profile specified in Cargo.toml to make +sure that LLVM rips out everything referring to `rust_eh_personality`, lest +the linker complain about an undefined symbol. + +[rf]: https://github.com/jdahlstrom/retrofire + +### Dependencies + +The demo is `#[no_std]` and the only thing it needs from `alloc` is `Vec`. +It provides a global allocator that is a minimal wrapper of libc `malloc` +and `free`, and uses `libc` `puts` and `putchar` to output the image. The +program has neither direct nor transitive dependencies beyond `alloc`, +`libc`, and`retrofire-core`. + +### Executable size + +On macOS `x86_64-apple-darwin`, the size of the compiled executable is around +20 kB when built in release mode. Of that, the `text` section takes up 10 kB +or so. On Windows 10 `x86_64-pc-windows-msvc` the executable size is roughly +16 kB. A 4k intro this is not, but not a bad result given that apart from +ripping out`std`, no specific manual size optimization has been done. + +### The output + +![test](output.png) + +### Acknowledgements + +Big thanks to user @jonthagen for their awesome [`min-sized-rust`][msr] +guide to minimizing the size of Rust binaries. + +[msr]: https://github.com/johnthagen/min-sized-rust + +### License + +Licensed under the MIT license, found in [LICENSE-MIT](LICENSE-MIT). diff --git a/demos/nostd/output.png b/demos/nostd/output.png new file mode 100644 index 00000000..3a0d4e65 Binary files /dev/null and b/demos/nostd/output.png differ diff --git a/demos/nostd/src/main.rs b/demos/nostd/src/main.rs new file mode 100644 index 00000000..40abb6f4 --- /dev/null +++ b/demos/nostd/src/main.rs @@ -0,0 +1,77 @@ +#![no_std] +#![no_main] + +extern crate alloc; + +use alloc::alloc::*; +use core::{ffi::c_void, panic::PanicInfo}; + +use libc::{abort, c_char, c_int, free, malloc, putchar, puts}; + +use re::prelude::*; + +use re::math::mat::ProjMat3; + +#[global_allocator] +static ALLOC: Malloc = Malloc; + +struct Malloc; + +unsafe impl GlobalAlloc for Malloc { + unsafe fn alloc(&self, layout: Layout) -> *mut u8 { + unsafe { malloc(layout.size()) as *mut u8 } + } + unsafe fn dealloc(&self, ptr: *mut u8, _: Layout) { + unsafe { free(ptr as *mut c_void) } + } +} + +#[panic_handler] +unsafe fn panic(_info: &PanicInfo) -> ! { + unsafe { abort() } +} + +#[unsafe(no_mangle)] +fn main() -> i32 { + let verts = [ + vertex(pt3(-1.0, -1.0, 0.0), rgb(1.0, 0.0, 0.0)), + vertex(pt3(0.0, 1.0, 0.0), rgb(0.4, 0.4, 1.0)), + vertex(pt3(1.0, -1.0, 0.0), rgb(0.0, 0.8, 0.0)), + ]; + + let shader = shader::new( + |v: Vertex3, mvp: &ProjMat3| { + vertex(mvp.apply(&v.pos), v.attrib) + }, + |frag: Frag>| frag.var.to_color4(), + ); + + let dims @ (w, h) = (640, 480); + let modelview = translate3(0.0, 0.0, 2.0).to(); + let project = perspective(1.0, w as f32 / h as f32, 0.1..1000.0); + let viewport = viewport(pt2(0, h)..pt2(w, 0)); + + let mut framebuf = Buf2::::new(dims); + + render( + [tri(0, 1, 2)], + verts, + &shader, + &modelview.then(&project), + viewport, + &mut framebuf, + &Context::default(), + ); + + unsafe { + puts("P6\n640 480 255\n\0".as_ptr() as *const c_char); + } + for &col in framebuf.data() { + unsafe { + putchar(col.r() as c_int); + putchar(col.g() as c_int); + putchar(col.b() as c_int); + }; + } + 0 +} diff --git a/demos/src/bin/bezier.rs b/demos/src/bin/bezier.rs index 6fe1fc27..26823715 100644 --- a/demos/src/bin/bezier.rs +++ b/demos/src/bin/bezier.rs @@ -2,9 +2,12 @@ use core::ops::ControlFlow::Continue; use re::prelude::*; -use re::core::geom::{Edge, Ray}; -use re::core::math::rand::{Distrib, Uniform, VectorsOnUnitDisk, Xorshift64}; -use re::core::render::raster::line; +use re::core::{ + geom::{Edge, Ray}, + math::rand::{Distrib, Uniform, VectorsOnUnitDisk, Xorshift64}, + math::spline::{Euclidean, approximate}, + render::raster::line, +}; use re::front::{Frame, dims, minifb::Window}; fn main() { @@ -24,18 +27,18 @@ fn main() { let vel = VectorsOnUnitDisk; let mut pos_vels: Vec<(Point2, Vec2)> = - (pos, vel).samples(rng).take(32).collect(); + (pos, vel).samples(rng).take(16).collect(); // Disable some unneeded things win.ctx.color_clear = None; win.ctx.depth_clear = None; - win.run(|Frame { dt, buf, .. }| { + win.run(|Frame { t, dt, buf, .. }| { let buf = &mut buf.borrow_mut().color_buf.buf; // Fade out previous frame a bit buf.iter_mut() - .for_each(|c| *c = c.saturating_sub(0x08_08_02)); + .for_each(|c| *c = c.saturating_sub(0xFF_FF_FF)); let rays: Vec> = pos_vels .chunks(2) @@ -44,15 +47,26 @@ fn main() { let b = BezierSpline::from_rays(rays); // Stop once error is less than one pixel - let approx = b.approximate(1.0); + let approx = approximate(&b, 1.0); for Edge(p0, p1) in approx.edges() { let vs = [p0, p1].map(|p| vertex(p.to_pt3().to(), ())); line(vs, |sl| { - buf[sl.y][sl.xs].fill(0xFF_FF_FF); + buf[sl.y][sl.xs].fill(0x33_00_00); }) } + let euc = Euclidean::new(b); + for s in 0.0.vary_to(euc.len(), 32) { + let p = euc.eval((s + 40.0 * t.as_secs_f32()) % euc.len()); + let x = p.x() as usize; + let y = p.y() as usize; + buf[y - 1][x] = 0xFF_FF_FF; + buf[y][x - 1] = 0xFF_FF_FF; + buf[y + 1][x] = 0xFF_FF_FF; + buf[y][x + 1] = 0xFF_FF_FF; + } + let dt = dt.as_secs_f32(); for (pos, vel) in &mut pos_vels { *pos = (*pos + 80.0 * *vel * dt).clamp(&min, &max); diff --git a/demos/src/bin/crates.rs b/demos/src/bin/crates.rs index e5590da6..d05d8a7e 100644 --- a/demos/src/bin/crates.rs +++ b/demos/src/bin/crates.rs @@ -4,9 +4,11 @@ use re::prelude::*; use re::core::math::color::gray; use re::core::render::{ + Model, cam::{FirstPerson, Fov}, clip::Status::*, scene::Obj, + shader, tex::SamplerClamp, }; // Try also Rgb565 or Rgba4444 @@ -124,13 +126,8 @@ fn main() { // TODO Try to get rid of clone .clone() .mesh(geom) - .uniform(&model_to_project) - // TODO Allow setting shader before uniform .shader(crate_shader) - // TODO storing &mut target makes Batch not Clone, maybe - // pass to render() instead. OTOH then a Frame::batch - // helper wouldn't be as useful. Maybe just wrap the - // target in a RefCell? + .uniform(&model_to_project) .render(); frame.ctx.stats.borrow_mut().objs.o += 1; @@ -138,7 +135,7 @@ fn main() { Continue(()) }) - .expect("should run") + .expect("should run"); } fn crates() -> Vec> { @@ -150,6 +147,7 @@ fn crates() -> Vec> { for j in (-n..=n).step_by(5) { res.push(Obj { tf: translate3(i as f32, 0.0, j as f32).to(), + // TODO Same geometry cloned many times ..obj.clone() }); } diff --git a/demos/src/bin/curses.rs b/demos/src/bin/curses.rs index fc102235..4661ad08 100644 --- a/demos/src/bin/curses.rs +++ b/demos/src/bin/curses.rs @@ -5,7 +5,8 @@ use pancurses::*; use re::prelude::*; use re::core::render::{ - ctx::DepthSort::BackToFront, raster::Scanline, stats::Throughput, + Model, ctx::DepthSort::BackToFront, raster::Scanline, render, shader, + stats::Throughput, }; use re::geom::solids::{Build, Torus}; diff --git a/demos/src/bin/hello.rs b/demos/src/bin/hello.rs index 37790783..55c8227b 100644 --- a/demos/src/bin/hello.rs +++ b/demos/src/bin/hello.rs @@ -3,7 +3,7 @@ use std::{env, fmt::Write, ops::ControlFlow::Continue}; use re::prelude::*; use re::core::{ - render::{Text, tex::Atlas, tex::Layout}, + render::{Model, Text, World, render, shader, tex::Atlas, tex::Layout}, util::pnm::parse_pnm, }; diff --git a/demos/src/bin/solids.rs b/demos/src/bin/solids.rs index 777a1540..59d63504 100644 --- a/demos/src/bin/solids.rs +++ b/demos/src/bin/solids.rs @@ -4,12 +4,14 @@ use minifb::{Key, KeyRepeat}; use re::prelude::*; -use re::core::geom::Polyline; -use re::core::math::{ProjMat3, ProjVec3, color::gray}; -use re::core::render::cam::Fov; - +use re::core::{ + geom::Polyline, + math::{ProjMat3, ProjVec3, color::gray}, + render::cam::Fov, + render::{Model, ModelToWorld, shader}, +}; use re::front::{Frame, minifb::Window}; -use re::geom::{io::parse_obj, solids::*}; +use re::geom::{io::read_obj, solids::*}; // Carousel animation for switching between objects. #[derive(Default)] @@ -67,7 +69,7 @@ fn main() { fn vtx_shader(v: VertexIn, (mvp, spin): Uniform) -> VertexOut { // Transform vertex normal - let norm = spin.apply(&v.attrib.to()); + let norm = spin.apply(&v.attrib); // Calculate diffuse shading let diffuse = (norm.z() + 0.2).max(0.2) * 0.8; // Visualize normal by mapping to RGB values @@ -108,14 +110,16 @@ fn main() { let object = &objects[carousel.idx % objects.len()]; - Batch::new() - .mesh(object) - .uniform((&model_view_project, &spin)) - .shader(shader) - .viewport(cam.viewport) - .target(frame.buf) - .context(&*frame.ctx) - .render(); + Batch { + prims: object.faces.clone(), + verts: object.verts.clone(), + uniform: (&model_view_project, &spin), + shader: shader, + viewport: cam.viewport, + target: frame.buf, + ctx: &*frame.ctx, + } + .render(); Continue(()) }); @@ -173,20 +177,21 @@ fn lathe(secs: u32) -> Mesh { // Loads the Utah teapot model. fn teapot() -> Mesh { - parse_obj(*include_bytes!("../../assets/teapot.obj")) + static TEAPOT: &[u8] = include_bytes!("../../assets/teapot.obj"); + read_obj(TEAPOT) .unwrap() .transform( &scale(splat(0.4)) .then(&translate(-0.5 * Vec3::Y)) .to(), ) - //.with_vertex_normals() .build() } // Loads the Stanford bunny model. fn bunny() -> Mesh { - parse_obj::<()>(*include_bytes!("../../assets/bunny.obj")) + static BUNNY: &[u8] = include_bytes!("../../assets/bunny.obj"); + read_obj::<()>(BUNNY) .unwrap() .transform(&scale(splat(0.12)).then(&translate(-Vec3::Y)).to()) .with_vertex_normals() @@ -196,7 +201,7 @@ fn bunny() -> Mesh { // Loads the Stanford dragon model. fn dragon() -> Mesh { static DRAGON: &[u8] = include_bytes!("../../assets/dragon.obj"); - parse_obj::<()>(DRAGON.iter().copied()) + read_obj::<()>(DRAGON) .unwrap() .with_vertex_normals() .transform( diff --git a/demos/src/bin/sprites.rs b/demos/src/bin/sprites.rs index 290ad058..2892bcc3 100644 --- a/demos/src/bin/sprites.rs +++ b/demos/src/bin/sprites.rs @@ -6,7 +6,7 @@ use re::core::math::{ color::gray, rand::{Distrib, PointsInUnitBall, Xorshift64}, }; -use re::core::render::{Model, cam::*, render}; +use re::core::render::{Model, View, cam::*, render, shader}; use re_front::minifb::Window; @@ -20,7 +20,7 @@ fn main() { let count = 10000; let rng = &mut Xorshift64::default(); - let verts: Vec>> = PointsInUnitBall + let verts: Vec>> = PointsInUnitBall .samples(rng) .take(count) .flat_map(|pos| verts.map(|v| vertex(pos.to(), v))) @@ -39,7 +39,7 @@ fn main() { let shader = shader::new( |v: Vertex3>, (mv, proj): (&Mat4, &ProjMat3)| { - let vertex_pos = 0.008 * v.attrib.to_vec3().to(); + let vertex_pos = 0.008 * v.attrib.to_vec3().to(); // Model->View let view_pos = mv.apply(&v.pos) + vertex_pos; vertex(proj.apply(&view_pos), v.attrib) }, @@ -64,7 +64,7 @@ fn main() { let modelview = rotate_x(theta * 0.2) .then(&rotate_z(theta * 0.14)) .to() - .then(&cam.transform.world_to_view()); + .then(&cam.world_to_view()); render( &tris, diff --git a/demos/src/bin/square.rs b/demos/src/bin/square.rs index 09a920aa..aa2b66cf 100644 --- a/demos/src/bin/square.rs +++ b/demos/src/bin/square.rs @@ -2,7 +2,7 @@ use core::ops::ControlFlow::*; use re::prelude::*; -use re::core::render::tex::SamplerClamp; +use re::core::render::{render, shader, tex::SamplerClamp}; use re::front::minifb::Window; fn main() { diff --git a/front/Cargo.toml b/front/Cargo.toml index 21da7a15..2d22f058 100644 --- a/front/Cargo.toml +++ b/front/Cargo.toml @@ -19,6 +19,7 @@ license.workspace = true keywords.workspace = true categories.workspace = true repository.workspace = true +documentation.workspace = true [features] wasm = ["dep:wasm-bindgen", "dep:web-sys"] diff --git a/front/src/minifb.rs b/front/src/minifb.rs index 7013025f..9c8184b8 100644 --- a/front/src/minifb.rs +++ b/front/src/minifb.rs @@ -10,7 +10,7 @@ use std::time::Instant; use minifb::{Key, WindowOptions}; use retrofire_core::{ - render::{Colorbuf, Context, target}, + render::{Colorbuf, Context, Stats, target}, util::{Dims, buf::Buf2, buf::MutSlice2, pixfmt::Xrgb8888}, }; @@ -111,9 +111,9 @@ impl Window { /// * the user closes the window via the GUI (e.g. titlebar close button); /// * the Esc key is pressed; or /// * the callback returns `ControlFlow::Break`. - pub fn run(&mut self, mut frame_fn: F) + pub fn run(&mut self, mut frame_fn: F) -> Stats where - F: FnMut(&mut Frame>) -> ControlFlow<()>, + F: FnMut(&mut Frame>) -> ControlFlow<()>, { let (w, h) = self.dims; let mut cbuf = Buf2::new((w, h)); @@ -146,7 +146,9 @@ impl Window { ctx.stats.borrow_mut().frames += 1.0; } - println!("{}", ctx.stats.borrow()); + let stats = ctx.stats.into_inner(); + println!("{stats}"); + stats } fn should_quit(&self) -> bool { diff --git a/front/src/sdl2.rs b/front/src/sdl2.rs index 01dc8a7e..d340c3bc 100644 --- a/front/src/sdl2.rs +++ b/front/src/sdl2.rs @@ -13,7 +13,7 @@ use sdl2::{ use retrofire_core::math::{Color4, Vary}; use retrofire_core::render::{ - Colorbuf, Context, FragmentShader, Target, raster::Scanline, + Colorbuf, Context, FragmentShader, Stats, Target, raster::Scanline, stats::Throughput, target::rasterize_fb, }; use retrofire_core::util::{ @@ -181,7 +181,7 @@ impl, const N: usize> Window { /// * the user closes the window via the GUI (e.g. a title bar button); /// * the Esc key is pressed; or /// * the callback returns [`ControlFlow::Break`][ControlFlow]. - pub fn run(&mut self, mut frame_fn: F) -> Result<(), Error> + pub fn run(&mut self, mut frame_fn: F) -> Result where F: FnMut(&mut Frame>>) -> ControlFlow<()>, Color4: IntoPixel, @@ -244,8 +244,9 @@ impl, const N: usize> Window { break; } } - println!("{}", ctx.stats.borrow()); - Ok(()) + let stats = ctx.stats.into_inner(); + println!("{stats}"); + Ok(stats) } } diff --git a/geom/Cargo.toml b/geom/Cargo.toml index e9446b4c..5084220b 100644 --- a/geom/Cargo.toml +++ b/geom/Cargo.toml @@ -19,6 +19,7 @@ license.workspace = true keywords.workspace = true categories.workspace = true repository.workspace = true +documentation.workspace = true [dependencies] retrofire-core = { version = "0.4.0", path = "../core", default-features = false } diff --git a/geom/src/io.rs b/geom/src/io.rs index 86076d31..885f9e6d 100644 --- a/geom/src/io.rs +++ b/geom/src/io.rs @@ -134,7 +134,7 @@ where /// Parses an OBJ format mesh from an iterator. /// /// # Errors -/// Returns [`self::Error`] if OBJ parsing fails. +/// Returns [`self::Error`][Error] if OBJ parsing fails. pub fn parse_obj(src: impl IntoIterator) -> Result> where Builder: TryFrom, diff --git a/geom/src/isect.rs b/geom/src/isect.rs new file mode 100644 index 00000000..e6e579db --- /dev/null +++ b/geom/src/isect.rs @@ -0,0 +1,436 @@ +use core::fmt::Debug; + +#[cfg(feature = "std")] // TODO separate fp feature for geom +use retrofire_core::geom::Sphere; +use retrofire_core::{ + geom::{Plane3, Ray, Ray3}, + math::{ApproxEq, Point3, vec3}, + render::scene::BBox, +}; + +/// Trait for calculating whether and at which points two objects intersect. +pub trait Intersect { + /// The result of an intersection test. + type Result; + + /// Finds the point(s) where `self` and another object intersect, if any. + /// + /// It is implementation defined whether this method returns all the + /// intersection points or, for instance, only the closest one. + fn intersect(&self, other: &T) -> Self::Result; +} + +type RayIntersect3 = Option<(f32, Point3)>; + +impl Intersect> for Ray3 { + type Result = RayIntersect3; + + /// Returns the unique intersection point of `self` and a plane, + /// or `None` if they do not intersect. + /// + /// If an intersection point exists, returns `Some((t, point))`, where + /// `point` is the intersection point and `t` is the ray parameter value + /// such that `self.orig + t * self.dir == point`. + fn intersect(&self, p: &Plane3) -> Self::Result { + let Self(orig, dir) = self; + + // TODO checking two very unlikely conditions + + let num = p.signed_dist(*orig); + if num.approx_eq(&0.0) { + // Origin point coincident with the plane + return Some((0.0, *orig)); + } + let denom = dir.dot(&p.normal().to()); + if denom.approx_eq(&0.0) { + // Ray parallel with but not coincident with the plane + // (or dir is a zero vector) -> no intersection + return None; + } + let t = -num / denom; + if t.approx_le(&0.0) { + // Ray points away from the plane, intersection "behind" it + return None; + } + Some((t, *orig + t * *dir)) + } +} + +impl Intersect> for Ray3 { + type Result = RayIntersect3; // Only closest for now + + /// Returns the nearest intersection point of `self` and a box, + /// or `None` if they do not intersect. + /// + /// If an intersection point exists, returns `Some((t, point))`, where + /// `point` is the intersection point and `t` is the ray parameter value + /// such that `self.orig + t * self.dir == point`. + fn intersect(&self, bbox: &BBox) -> Self::Result { + let &BBox(low, upp) = bbox; + let Ray(orig, dir) = *self; + + // Ray equation: + // p(t) = O + d·t + // x(t) = Ox + dx·t + // y(t) = Oy + dy·t + // z(t) = Oz + dz·t + // + // Plane equations: + // x = lx, x = ux + // y = ly, x = uy + // z = lz, x = uz + // + // For each slab, ie. pair of parallel planes: + // Substitute eg. + // x_l = Ox + x_d·t0 + // x_u = Ox + x_d·t1 + // + // t0 = (x_l - x_O) / x_d | x_d=0 iff ray parallel with planes + // t1 = (x_u - x_O) / x_d + // + // Same for y and z slabs + + if bbox.is_empty() { + return None; + } + + let r_d = vec3(1.0 / dir.x(), 1.0 / dir.y(), 1.0 / dir.z()); + let low = (low - orig) * r_d; + let upp = (upp - orig) * r_d; + + let near = low.zip_map(upp, |l, u| l.min(u)); + let far = low.zip_map(upp, |l, u| l.max(u)); + + let near_t = near[0].max(near[1]).max(near[2]); + let far_t = far[0].min(far[1]).min(far[2]); + + if far_t.is_infinite() { + return None; + } + + if near_t > far_t || far_t < 0.0 { + // ---max---min--- (misses the box) or + // ---min---max---0---> (box behind ray) + return None; + } + let t = if near_t >= 0.0 { + // ---0---min---max---> (hits box) + near_t + } else { + // ---min---0---max---> (inside the box, unlikely) + far_t + }; + Some((t, orig + t * dir)) + } +} + +#[cfg(feature = "std")] +impl Intersect> for Ray3 { + type Result = RayIntersect3; // Only closest for now + + /// Returns the intersection point of `self` and a sphere closest to the + /// origin of `self`, or `None` if they do not intersect. + /// + /// # Examples + /// ``` + /// ``` + fn intersect(&self, &Sphere(center, r): &Sphere) -> Self::Result { + let &Ray(orig, dir) = self; + + // > r, no intersection + // If |C - C'| = r, one -"- + // < r, two -"- + // _______ + // / \ + // / \ + // | C--r--| + // \ | / + // \___|___/ + // _| + // O--------+-C'----> d + // + // Find point P = (x, y, z) given sphere (C, r) and ray (O, d) + // + // Sphere equation: + // (x - c_x)² + (y - c_y)² + (z - c_z)² = r² + // or in vector form + // (P - C) · (P - C) = r² + // + // Ray equation: + // P = O + t·d + // + // Intersection: + // + // Substitute ray equation to sphere equation: + // (o + t·d - c) · (o + t·d - c) = r² + // + // Multiply out + // o (o + td - c) + t·d (o + td - c) - c (o + td - c) = r² + // + // Distribute + // o·o + o·td - o·c + o·td + td·td - td·c - c·o - c·td + c·c = r² + // + // Reorder + // td·td + o·td + o·td - td·c + o·o - o·c - c·o + c·c - r² = 0 + // + // Factor out t's + // (d·d) t² + (o·d + o·d - c·d) t + o·o - 2o·c + c·c - r² = 0 + // + // Solve quadratic equation: + // (d·d) t² + 2(o - c)·d t + (o - c)² - r² = 0 + // + // t = (-b ± √(b² - 4ac)) / 2a + + let c_to_o = orig - center; + let a = dir.len_sqr(); // >= 0 + let b = 2.0 * c_to_o.dot(&dir); + let c = c_to_o.len_sqr() - r * r; + + let discriminant = b * b - 4.0 * a * c; + if discriminant < 0.0 { + // the line of the ray does not hit the sphere + return None; + } + + use retrofire_core::math::float::f32; + let sqrt = f32::sqrt(discriminant); + // sqrt >= 0.0, thus t0 <= t1 always + let (t0, t1) = (-b - sqrt, -b + sqrt); + let t = if t0 >= 0.0 { + // ray hits both points + t0 + } else if t1 >= 0.0 { + // ray origin is inside sphere + t1 + } else { + // sphere is behind ray + return None; + }; + let t = t / (2.0 * a); + Some((t, orig + t * dir)) + } +} + +#[cfg(test)] +mod tests { + use retrofire_core::math::{Linear, Vec3, pt3}; + + use super::*; + + mod ray_plane { + use super::*; + + const PLANE: Plane3<()> = Plane3::new(0.0, 1.0, 0.0, 2.0); + + #[test] + #[ignore] + fn ray_plane_xxx() { + let r = Ray::( + pt3(-3.308549, 6.2584567, -3.351655), + vec3(3.308549, -6.2584567, 3.351655), + ); + let bbox = BBox(pt3(-1.0, -1.0, -1.0), pt3(1.0, 1.0, 1.0)); + + assert_eq!(r.intersect(&bbox), Some((0.0, pt3(0.0, 0.0, 0.0)))); + } + + #[test] + fn ray_towards_plane_has_intersection() { + // Outside + let r = Ray(pt3(0.0, 3.0, 0.0), vec3(1.0, -1.0, 1.0)); + assert_eq!(r.intersect(&PLANE), Some((1.0, pt3(1.0, 2.0, 1.0)))); + + // Inside + let r = Ray(pt3(0.0, 1.0, 0.0), vec3(1.0, 1.0, 1.0)); + assert_eq!(r.intersect(&PLANE), Some((1.0, pt3(1.0, 2.0, 1.0)))); + } + #[test] + fn ray_origin_on_plane_has_intersection() { + let r = Ray(pt3(0.0, 2.0, 0.0), vec3(1.0, -1.0, 1.0)); + assert_eq!(r.intersect(&PLANE), Some((0.0, pt3(0.0, 2.0, 0.0)))); + } + #[test] + fn ray_coincident_with_plane_has_intersection() { + let r = Ray(pt3(0.0, 2.0, 0.0), vec3(1.0, 0.0, 1.0)); + assert_eq!(r.intersect(&PLANE), Some((0.0, pt3(0.0, 2.0, 0.0)))); + } + #[test] + fn ray_parallel_with_plane_no_intersection() { + let r = Ray(pt3(0.0, 3.0, 0.0), vec3(1.0, 0.0, 1.0)); + assert_eq!(r.intersect(&PLANE), None); + } + #[test] + fn ray_points_away_from_plane_no_intersection() { + // Outside + let r = Ray(pt3(0.0, 3.0, 0.0), vec3(1.0, 1.0, 1.0)); + assert_eq!(r.intersect(&PLANE), None); + + // Inside + let r = Ray(pt3(0.0, 1.0, 0.0), vec3(1.0, -1.0, 1.0)); + assert_eq!(r.intersect(&PLANE), None); + } + #[test] + fn degenerate_ray_only_intersects_if_coincident() { + let r = Ray(pt3(0.0, 3.0, 0.0), Vec3::zero()); + assert_eq!(r.intersect(&PLANE), None); + + let r = Ray(pt3(0.0, 1.0, 0.0), Vec3::zero()); + assert_eq!(r.intersect(&PLANE), None); + + let r = Ray(pt3(0.0, 2.0, 0.0), Vec3::zero()); + assert_eq!(r.intersect(&PLANE), Some((0.0, pt3(0.0, 2.0, 0.0)))); + } + } + + mod ray_bbox { + use super::*; + + const BBOX: BBox<()> = BBox(pt3(-1.0, -1.0, -1.0), pt3(1.0, 1.0, 1.0)); + + #[test] + fn parallel_simple() { + // +----+ + // x--> | | + // +----+ + let ray = Ray(pt3(0.0, 0.0, -2.0), vec3(0.0, 0.0, 1.0)); + assert_eq!(ray.intersect(&BBOX), Some((1.0, pt3(0.0, 0.0, -1.0)))); + } + #[test] + fn diagonal_simple() { + // + // +----+ + // ^| | + // / +----+ + // x + let ray = Ray(pt3(-1.5, 0.0, -2.0), vec3(1.0, 0.0, 1.0)); + assert_eq!(ray.intersect(&BBOX), Some((1.0, pt3(-0.5, 0.0, -1.0)))); + } + #[test] + #[ignore] + fn parallel_intersect_at_vertex() { + // x--> ,_____. + // / /| + // /_____/ | + // | | / + // |_____|/ + let ray = Ray(pt3(1.0, 1.0, -2.0), vec3(0.0, 0.0, 1.0)); + assert_eq!(ray.intersect(&BBOX), Some((1.0, pt3(1.0, 1.0, -1.0)))); + } + #[test] + #[ignore] + fn parallel_intersect_at_edge() { + // ,_____. + // x--> / /| + // /_____/ | + // | | / + // |_____|/ + // + let ray = Ray(pt3(0.0, 1.0, -2.0), vec3(0.0, 0.0, 1.0)); + assert_eq!(ray.intersect(&BBOX), Some((1.0, pt3(0.0, 1.0, -1.0)))); + } + #[test] + fn ray_starts_inside() { + // +--^----+ + // | | | + // | x | + // +-------+ + let ray = Ray(pt3(0.0, 0.0, -0.5), vec3(0.0, 1.0, 0.0)); + assert_eq!(ray.intersect(&BBOX), Some((1.0, pt3(0.0, 1.0, -0.5)))); + } + #[test] + fn ray_starts_on_side_plane() { + // Points away + // +-----+ + // <--x | + // +-----+ + let ray = Ray(pt3(0.0, 0.0, -1.0), vec3(0.0, 0.0, -1.0)); + assert_eq!(ray.intersect(&BBOX), Some((0.0, pt3(0.0, 0.0, -1.0)))); + // Points inside + // +-----+ + // x--> | + // +-----+ + let ray = Ray(pt3(0.0, 0.0, -1.0), vec3(0.0, 0.0, 1.0)); + assert_eq!(ray.intersect(&BBOX), Some((0.0, pt3(0.0, 0.0, -1.0)))); + } + #[test] + fn no_intersection() { + // Diagonal ray + // ^ + // / +----+ + // x | | + // +----+ + let ray = Ray(pt3(0.0, 0.0, -2.5), vec3(0.0, 1.0, 1.0)); + assert_eq!(ray.intersect(&BBOX), None); + + // Parallel but offset ray + // x---> + // +----+ + // | | + // +----+ + let ray = Ray(pt3(0.0, 1.5, -2.0), vec3(0.0, 0.0, 1.0)); + assert_eq!(ray.intersect(&BBOX), None); + } + #[test] + fn opposite_direction() { + // +----+ + // | | x---> + // +----+ + let ray = Ray(pt3(0.0, 0.0, 2.0), vec3(0.0, 0.0, 1.0)); + assert_eq!(ray.intersect(&BBOX), None); + } + #[test] + fn zero_length_ray() { + let ray = Ray(pt3(0.0, 0.0, -2.0), vec3(0.0, 0.0, 0.0)); + assert_eq!(ray.intersect(&BBOX), None); + } + #[test] + fn empty_box() { + let empty = BBox::<()>(pt3(-1.0, -1.0, 1.0), pt3(1.0, 1.0, -1.0)); + let ray = Ray(pt3(0.0, 0.0, -2.0), vec3(0.0, 0.0, 1.0)); + assert_eq!(ray.intersect(&empty), None); + } + } + + // TODO until sqrt has a fallback + #[cfg(feature = "std")] + mod ray_sphere { + use super::*; + + const SPHERE: Sphere = Sphere(pt3(0.0, 0.0, 1.0), 2.0); + + #[test] + fn ray_passes_through_sphere() { + let ray: Ray3 = Ray(pt3(0.0, 0.0, -3.0), vec3(0.0, 0.0, 2.0)); + assert_eq!( + ray.intersect(&SPHERE), + Some((1.0, pt3(0.0, 0.0, -1.0))) + ); + } + #[test] + fn ray_tangent_to_sphere() { + let ray: Ray3 = Ray(pt3(0.0, 2.0, -3.0), vec3(0.0, 0.0, 2.0)); + assert_eq!(ray.intersect(&SPHERE), Some((2.0, pt3(0.0, 2.0, 1.0)))); + } + + #[test] + fn ray_origin_inside_sphere() { + let ray: Ray3 = Ray(pt3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 2.0)); + assert_eq!(ray.intersect(&SPHERE), Some((1.5, pt3(0.0, 0.0, 3.0)))); + + let ray: Ray3 = Ray(pt3(0.0, 0.0, 2.0), vec3(0.0, 0.0, 2.0)); + assert_eq!(ray.intersect(&SPHERE), Some((0.5, pt3(0.0, 0.0, 3.0)))); + } + + #[test] + fn sphere_behind_ray() { + let ray: Ray3 = Ray(pt3(0.0, 0.0, -3.0), vec3(0.0, 0.0, -1.0)); + assert_eq!(ray.intersect(&SPHERE), None); + } + + #[test] + fn ray_misses_sphere() { + let ray: Ray3 = Ray(pt3(0.0, 0.0, -3.0), vec3(0.0, 1.0, 1.0)); + assert_eq!(ray.intersect(&SPHERE), None); + } + } +} diff --git a/geom/src/lib.rs b/geom/src/lib.rs index 692111bc..35e4e144 100644 --- a/geom/src/lib.rs +++ b/geom/src/lib.rs @@ -6,4 +6,7 @@ extern crate core; extern crate std; pub mod io; +pub mod isect; pub mod solids; + +pub use isect::Intersect;