From 2a41743fe516affabb876acede97be84de46e87e Mon Sep 17 00:00:00 2001 From: numzero Date: Sat, 14 Sep 2024 15:02:03 +0300 Subject: [PATCH] fmt! --- rustfmt.toml | 1 + src/bin/flat/float_fun.rs | 8 +- src/bin/flat/fns.rs | 131 +++++++++++--- src/bin/flat/main.rs | 115 +++++++++---- src/bin/flat/riemann.rs | 111 +++++++++--- src/bin/flat/tube/coords.rs | 330 ++++++++++++++++++++++++++++++------ src/bin/flat/tube/metric.rs | 91 ++++++---- src/bin/flat/tube/mod.rs | 223 ++++++++++++++++++++---- src/bin/flat/types.rs | 12 +- src/bin/mesh/main.rs | 43 +++-- src/bin/mesh/mesh_loader.rs | 34 ++-- 11 files changed, 881 insertions(+), 218 deletions(-) create mode 100644 rustfmt.toml diff --git a/rustfmt.toml b/rustfmt.toml new file mode 100644 index 0000000..218e203 --- /dev/null +++ b/rustfmt.toml @@ -0,0 +1 @@ +hard_tabs = true diff --git a/src/bin/flat/float_fun.rs b/src/bin/flat/float_fun.rs index 82411e9..e9eb4c3 100644 --- a/src/bin/flat/float_fun.rs +++ b/src/bin/flat/float_fun.rs @@ -12,8 +12,12 @@ pub trait FloatExt2: bounds::Pair { } impl FloatExt2 for (F, F) { - fn lerp(self, t: F) -> F { F::lerp(self.0, self.1, t) } - fn inverse_lerp(self, y: F) -> F { F::inverse_lerp(self.0, self.1, y) } + fn lerp(self, t: F) -> F { + F::lerp(self.0, self.1, t) + } + fn inverse_lerp(self, y: F) -> F { + F::inverse_lerp(self.0, self.1, y) + } } #[cfg(test)] diff --git a/src/bin/flat/fns.rs b/src/bin/flat/fns.rs index c52bdaa..077820b 100644 --- a/src/bin/flat/fns.rs +++ b/src/bin/flat/fns.rs @@ -52,29 +52,73 @@ pub struct QuadraticAccelerator { /// Продолжает функцию f с [-lim, lim] линейно в предположении f(±lim) = ±val, f'(±lim) = 1. fn extend_linear(t: f32, f: impl FnOnce(f32) -> f32, lim: f32, val: f32) -> f32 { - if t.abs() <= lim { f(t) } else { t + t.signum() * (val - lim) } + if t.abs() <= lim { + f(t) + } else { + t + t.signum() * (val - lim) + } } /// Продолжает функцию f с [-lim, lim] константой в предположении f(±lim) = val, f'(±lim) = 0. fn extend_const(t: f32, f: impl FnOnce(f32) -> f32, lim: f32, val: f32) -> f32 { - if t.abs() <= lim { f(t) } else { val } + if t.abs() <= lim { + f(t) + } else { + val + } } impl QuadraticAccelerator { - fn a(&self) -> f32 { -(self.external - self.internal) / self.internal.powi(2) } - fn b(&self) -> f32 { 2.0 * self.external / self.internal - 1.0 } - fn root(&self, x: f32) -> f32 { (self.b().powi(2) + 4.0 * self.a() * x.abs()).sqrt() } + fn a(&self) -> f32 { + -(self.external - self.internal) / self.internal.powi(2) + } + fn b(&self) -> f32 { + 2.0 * self.external / self.internal - 1.0 + } + fn root(&self, x: f32) -> f32 { + (self.b().powi(2) + 4.0 * self.a() * x.abs()).sqrt() + } - pub fn x(&self, u: f32) -> f32 { extend_linear(u, |u| (self.a() * u.abs() + self.b()) * u, self.internal, self.external) } - pub fn u(&self, x: f32) -> f32 { extend_linear(x, |x| 0.5 * x.signum() * (-self.b() + self.root(x)) / self.a(), self.external, self.internal) } - pub fn dx(&self, u: f32) -> f32 { extend_const(u, |u| 2.0 * self.a() * u.abs() + self.b(), self.internal, 1.0) } - pub fn du(&self, x: f32) -> f32 { extend_const(x, |x| 1.0 / self.root(x), self.external, 1.0) } - pub fn d2u(&self, x: f32) -> f32 { extend_const(x, |x| -2.0 * x.signum() * self.a() * self.root(x).powi(-3), self.external, 0.0) } + pub fn x(&self, u: f32) -> f32 { + extend_linear( + u, + |u| (self.a() * u.abs() + self.b()) * u, + self.internal, + self.external, + ) + } + pub fn u(&self, x: f32) -> f32 { + extend_linear( + x, + |x| 0.5 * x.signum() * (-self.b() + self.root(x)) / self.a(), + self.external, + self.internal, + ) + } + pub fn dx(&self, u: f32) -> f32 { + extend_const( + u, + |u| 2.0 * self.a() * u.abs() + self.b(), + self.internal, + 1.0, + ) + } + pub fn du(&self, x: f32) -> f32 { + extend_const(x, |x| 1.0 / self.root(x), self.external, 1.0) + } + pub fn d2u(&self, x: f32) -> f32 { + extend_const( + x, + |x| -2.0 * x.signum() * self.a() * self.root(x).powi(-3), + self.external, + 0.0, + ) + } } #[cfg(test)] mod test { - use approx::{abs_diff_eq, AbsDiffEq, assert_abs_diff_eq}; + use approx::{abs_diff_eq, assert_abs_diff_eq, AbsDiffEq}; use super::*; @@ -93,23 +137,45 @@ mod test { for x in itertools_num::linspace(-mul * max, mul * max, 100) { let df_num = (testee.value(x + δ) - testee.value(x - δ)) / (2. * δ); let df_expl = testee.derivative(x); - assert!(abs_diff_eq!(df_expl, df_num, epsilon = ε), "At x={x}, df/dx:\nnumerical: {df_num}\nexplicit: {df_expl}\n"); + assert!( + abs_diff_eq!(df_expl, df_num, epsilon = ε), + "At x={x}, df/dx:\nnumerical: {df_num}\nexplicit: {df_expl}\n" + ); } } #[test] fn test_smoothstep_limiter() { - test_limiter(SmoothstepLimiter { min: 20.0, max: 30.0 }, 20.0, 30.0, 1.0 / 32.0); + test_limiter( + SmoothstepLimiter { + min: 20.0, + max: 30.0, + }, + 20.0, + 30.0, + 1.0 / 32.0, + ); } #[test] fn test_smootherstep_limiter() { - test_limiter(SmootherstepLimiter { min: 20.0, max: 30.0 }, 20.0, 30.0, 1.0 / 32.0); + test_limiter( + SmootherstepLimiter { + min: 20.0, + max: 30.0, + }, + 20.0, + 30.0, + 1.0 / 32.0, + ); } #[test] fn test_quadratic_accelerator() { - let testee = super::QuadraticAccelerator { internal: 100.0, external: 150.0 }; + let testee = super::QuadraticAccelerator { + internal: 100.0, + external: 150.0, + }; let ε = 1.0e-4f32; let δ = 1.0 / 8.0; // Mathematically, you want this to be small. Computationally, you don’t. let margin = 1.0 / 16.0; @@ -121,30 +187,51 @@ mod test { for x in itertools_num::linspace(-mul * testee.external, mul * testee.external, 100) { let ux = testee.u(x); let xux = testee.x(ux); - assert!(abs_diff_eq!(x, xux, epsilon = ε), "At x={x}:\nu(x): {ux}\nx(u(x)): {xux}\n"); + assert!( + abs_diff_eq!(x, xux, epsilon = ε), + "At x={x}:\nu(x): {ux}\nx(u(x)): {xux}\n" + ); let du_num = (testee.u(x + δ) - testee.u(x - δ)) / (2. * δ); let du_expl = testee.du(x); - assert!(abs_diff_eq!(du_expl, du_num, epsilon = ε), "At x={x}, du/dx:\nnumerical: {du_num}\nexplicit: {du_expl}\n"); + assert!( + abs_diff_eq!(du_expl, du_num, epsilon = ε), + "At x={x}, du/dx:\nnumerical: {du_num}\nexplicit: {du_expl}\n" + ); let dudx = du_expl * testee.dx(ux); - assert!(abs_diff_eq!(dudx, 1.0, epsilon = ε), "At x={x}:\ndu/dx * dx/du: {dudx}\n"); + assert!( + abs_diff_eq!(dudx, 1.0, epsilon = ε), + "At x={x}:\ndu/dx * dx/du: {dudx}\n" + ); let d2u_num = (testee.du(x + δ) - testee.du(x - δ)) / (2. * δ); let d2u_expl = testee.d2u(x); - assert!(abs_diff_eq!(d2u_expl, d2u_num, epsilon = ε), "At x={x}, d^2u/dx^2:\nnumerical: {d2u_num}\nexplicit: {d2u_expl}\n"); + assert!( + abs_diff_eq!(d2u_expl, d2u_num, epsilon = ε), + "At x={x}, d^2u/dx^2:\nnumerical: {d2u_num}\nexplicit: {d2u_expl}\n" + ); } for u in itertools_num::linspace(-mul * testee.internal, mul * testee.internal, 100) { let xu = testee.x(u); let uxu = testee.u(xu); - assert!(abs_diff_eq!(u, uxu, epsilon = ε), "At u={u}:\nx(u): {xu}\nu(x(u)): {uxu}\n"); + assert!( + abs_diff_eq!(u, uxu, epsilon = ε), + "At u={u}:\nx(u): {xu}\nu(x(u)): {uxu}\n" + ); let dx_num = (testee.x(u + δ) - testee.x(u - δ)) / (2. * δ); let dx_expl = testee.dx(u); - assert!(abs_diff_eq!(dx_expl, dx_num, epsilon = ε), "At u={u}, dx/du:\nnumerical: {dx_num}\nexplicit: {dx_expl}\n"); + assert!( + abs_diff_eq!(dx_expl, dx_num, epsilon = ε), + "At u={u}, dx/du:\nnumerical: {dx_num}\nexplicit: {dx_expl}\n" + ); let dudx = testee.du(xu) * dx_expl; - assert!(abs_diff_eq!(dudx, 1.0, epsilon = ε), "At u={u}:\ndu/dx * dx/du: {dudx}\n"); + assert!( + abs_diff_eq!(dudx, 1.0, epsilon = ε), + "At u={u}:\ndu/dx * dx/du: {dudx}\n" + ); } } } diff --git a/src/bin/flat/main.rs b/src/bin/flat/main.rs index ee853a4..65f8b50 100644 --- a/src/bin/flat/main.rs +++ b/src/bin/flat/main.rs @@ -4,24 +4,26 @@ use flo_canvas::*; use flo_draw::*; use glam::*; -use riemann::{Metric, trace_iter}; +use crate::types::FlatTraceResult; +use riemann::{trace_iter, Metric}; use tube::metric::Tube; use tube::Space; use tube::Subspace::{Boundary, Inner, Outer}; use types::{Location, Object, Ray}; -use crate::types::FlatTraceResult; -mod riemann; -mod fns; mod float_fun; +mod fns; +mod riemann; mod tube; mod types; const DT: f32 = 0.1; -fn draw_loop(gc: &mut Vec, mut pts: impl Iterator) { +fn draw_loop(gc: &mut Vec, mut pts: impl Iterator) { gc.new_path(); - let Some(first) = pts.next() else { return; }; + let Some(first) = pts.next() else { + return; + }; gc.move_to(first.x, first.y); for pt in pts { gc.line_to(pt.x, pt.y); @@ -69,10 +71,21 @@ pub fn main() { gc.stroke_color(Color::Rgba(1.0, 0.5, 0.0, 1.0)); draw_fan_2(gc, &space, vec2(-500.0, 0.0), vec2(1.0, 0.0), 1.0); gc.stroke_color(Color::Rgba(0.5, 1.0, 0.0, 1.0)); - draw_fan_2(gc, &space, vec2(-2.5 * tube.outer_radius, 1.25 * tube.external_halflength), vec2(1.0, -1.0), 1.0); + draw_fan_2( + gc, + &space, + vec2(-2.5 * tube.outer_radius, 1.25 * tube.external_halflength), + vec2(1.0, -1.0), + 1.0, + ); draw_track(gc, &space, vec2(-500.0, 0.0), vec2(1.0, 0.2)); draw_track(gc, &space, vec2(-500.0, 0.0), vec2(1.0, 0.5)); - draw_track(gc, &space, vec2(-0.5 * tube.inner_radius, -1.25 * tube.external_halflength), vec2(0.1, 1.0)); + draw_track( + gc, + &space, + vec2(-0.5 * tube.inner_radius, -1.25 * tube.external_halflength), + vec2(0.1, 1.0), + ); let circle_segments = 47; for obj in &space.objs { @@ -83,25 +96,44 @@ pub fn main() { gc.fill(); gc.stroke_color(Color::Rgba(0.0, 0.0, 0.0, 0.5)); - draw_loop(gc, itertools_num::linspace(0.0, 2.0 * PI, circle_segments).skip(1).map(|φ| { - let dir = Vec2::from_angle(φ) * obj.r; - let dir = obj.loc.rot * dir; - pos + dir - })); + draw_loop( + gc, + itertools_num::linspace(0.0, 2.0 * PI, circle_segments) + .skip(1) + .map(|φ| { + let dir = Vec2::from_angle(φ) * obj.r; + let dir = obj.loc.rot * dir; + pos + dir + }), + ); gc.stroke_color(Color::Rgba(0.0, 0.5, 1.0, 0.5)); - draw_loop(gc, itertools_num::linspace(0.0, 2.0 * PI, circle_segments).skip(1).map(|φ| { - let dir = Vec2::from_angle(φ) * obj.r; - let dir = obj.loc.rot * dir; - space.trace_step(Ray { pos, dir }).pos - })); + draw_loop( + gc, + itertools_num::linspace(0.0, 2.0 * PI, circle_segments) + .skip(1) + .map(|φ| { + let dir = Vec2::from_angle(φ) * obj.r; + let dir = obj.loc.rot * dir; + space.trace_step(Ray { pos, dir }).pos + }), + ); gc.stroke_color(Color::Rgba(0.5, 0.0, 1.0, 1.0)); - draw_loop(gc, itertools_num::linspace(0.0, 2.0 * PI, circle_segments).skip(1).map(|φ| { - let n = obj.r.floor(); - let d = obj.r / n; - let dir = Vec2::from_angle(φ); - let dir = obj.loc.rot * dir * d; - space.trace_iter(Ray { pos, dir }).nth(n as usize).unwrap().pos - })); + draw_loop( + gc, + itertools_num::linspace(0.0, 2.0 * PI, circle_segments) + .skip(1) + .map(|φ| { + let n = obj.r.floor(); + let d = obj.r / n; + let dir = Vec2::from_angle(φ); + let dir = obj.loc.rot * dir * d; + space + .trace_iter(Ray { pos, dir }) + .nth(n as usize) + .unwrap() + .pos + }), + ); } }); }); @@ -109,7 +141,9 @@ pub fn main() { fn rel_to_abs(space: &impl Metric, base: &Location, rel: Vec2, steps: usize) -> Vec2 { let c = 1.0 / (steps as f32); - trace_iter(space, base.pos, base.rot * rel, c * rel.length()).nth(steps - 1).unwrap() + trace_iter(space, base.pos, base.rot * rel, c * rel.length()) + .nth(steps - 1) + .unwrap() } fn draw_cross(gc: &mut Vec, pos: Vec2, r: f32) { @@ -128,7 +162,7 @@ fn draw_ray_2(gc: &mut Vec, space: &Space, base: Vec2, dir: Vec2) { Outer => return (ray, space.trace_outer(ray)), Boundary => continue, }; - }; + } unreachable!("Space::trace_iter terminated!") } @@ -136,7 +170,10 @@ fn draw_ray_2(gc: &mut Vec, space: &Space, base: Vec2, dir: Vec2) { let dir = space.tube.globalize(base, dir); gc.new_path(); gc.move_to(base.x, base.y); - let mut ray = Ray { pos: base, dir: space.tube.normalize_vec_at(base, dir) * DT }; + let mut ray = Ray { + pos: base, + dir: space.tube.normalize_vec_at(base, dir) * DT, + }; for _ in 0..100 { let ret; (ray, ret) = trace_to_flat(gc, space, ray); @@ -150,7 +187,8 @@ fn draw_ray_2(gc: &mut Vec, space: &Space, base: Vec2, dir: Vec2) { let apx_hit_pos = rel_to_abs(&space.tube, &obj.loc, hit.rel.pos, 128); // assert_abs_diff_eq!(apx_hit_pos, hit.pos, epsilon=1.0); let Ray { pos: rel, dir } = hit.rel; - let diff = rel.dot(dir).powi(2) - dir.length_squared() * (rel.length_squared() - obj.r.powi(2)); + let diff = rel.dot(dir).powi(2) + - dir.length_squared() * (rel.length_squared() - obj.r.powi(2)); assert!(diff >= 0.0); let t = (-rel.dot(dir) + diff.sqrt()) / dir.length_squared(); let rel2 = hit.rel.forward(t).pos; @@ -207,7 +245,10 @@ fn draw_track(gc: &mut Vec, space: &Space, start: Vec2, dir: Vec2) { // let mut loc = Location { pos: start, rot: Mat2::IDENTITY }; // let dir = space.tube.globalize(start, dir); // let v = space.tube.normalize(start, dir); - let mut loc = Location { pos: start, rot: mat2(dir, vec2(-dir.y, dir.x)) }; + let mut loc = Location { + pos: start, + rot: mat2(dir, vec2(-dir.y, dir.x)), + }; let v = vec2(1.0, 0.0); let mut draw = |loc: &Location| { let p = loc.pos; @@ -249,8 +290,18 @@ trait Renderable { impl Renderable for Tube { fn render(&self, gc: &mut Vec) { gc.new_path(); - gc.rect(-self.outer_radius, -self.external_halflength, self.outer_radius, self.external_halflength); - gc.rect(-self.inner_radius, -self.external_halflength, self.inner_radius, self.external_halflength); + gc.rect( + -self.outer_radius, + -self.external_halflength, + self.outer_radius, + self.external_halflength, + ); + gc.rect( + -self.inner_radius, + -self.external_halflength, + self.inner_radius, + self.external_halflength, + ); gc.winding_rule(WindingRule::EvenOdd); gc.fill_color(Color::Rgba(0.8, 0.8, 0.8, 1.0)); gc.fill(); diff --git a/src/bin/flat/riemann.rs b/src/bin/flat/riemann.rs index ef07959..2acbe50 100644 --- a/src/bin/flat/riemann.rs +++ b/src/bin/flat/riemann.rs @@ -88,7 +88,11 @@ pub fn krist(space: &impl Metric, pos: Vec2) -> Tens2 { let g = &space.inverse_at(pos); // с верхними индексами let d = space.part_derivs_at(pos); // ret[i][l][k] = sum((m) => .5f * g[m][i] * (d[k][l][m] + d[l][k][m] - d[m][k][l])) - make_tens2(|i, l, k| 0.5 * (0..2).map(|m| g.col(m)[i] * (d[l].col(k)[m] + d[k].col(m)[l] - d[m].col(k)[l])).sum::()) + make_tens2(|i, l, k| { + 0.5 * (0..2) + .map(|m| g.col(m)[i] * (d[l].col(k)[m] + d[k].col(m)[l] - d[m].col(k)[l])) + .sum::() + }) } fn dir_deriv(f: impl Fn(Vec2) -> Mat2, pos: Vec2, delta: Vec2) -> Mat2 { @@ -167,40 +171,109 @@ pub mod samples { #[cfg(test)] mod tests { - use approx::assert_abs_diff_eq; use super::*; + use approx::assert_abs_diff_eq; - use glam::{Mat2, mat2, vec2}; + use glam::{mat2, vec2, Mat2}; use rand::{Rng, SeedableRng}; #[test] fn uniform_scaled_metric() { let mut rng = rand_pcg::Pcg64Mcg::seed_from_u64(17); - let metric = samples::ScaledMetric { scale: vec2(3., 4.) }; - assert_eq!(metric.sqrt_at(rng.gen()), Decomp2 { ortho: Mat2::IDENTITY, diag: vec2(3., 4.) }); - assert_eq!(metric.at(rng.gen()), Mat2::from_cols_array(&[9., 0., 0., 16.])); - assert_eq!(metric.inverse_at(rng.gen()), Mat2::from_cols_array(&[1. / 9., 0., 0., 1. / 16.])); + let metric = samples::ScaledMetric { + scale: vec2(3., 4.), + }; + assert_eq!( + metric.sqrt_at(rng.gen()), + Decomp2 { + ortho: Mat2::IDENTITY, + diag: vec2(3., 4.) + } + ); + assert_eq!( + metric.at(rng.gen()), + Mat2::from_cols_array(&[9., 0., 0., 16.]) + ); + assert_eq!( + metric.inverse_at(rng.gen()), + Mat2::from_cols_array(&[1. / 9., 0., 0., 1. / 16.]) + ); assert_eq!(metric.part_derivs_at(rng.gen()), [Mat2::ZERO, Mat2::ZERO]); assert_eq!(metric.vec_length_at(rng.gen(), vec2(1., 0.)), 3.); assert_eq!(metric.vec_length_at(rng.gen(), vec2(0., 1.)), 4.); assert_eq!(metric.vec_length_at(rng.gen(), vec2(1., 1.)), 5.); - assert_eq!(metric.normalize_vec_at(rng.gen(), vec2(1., 0.)), vec2(1. / 3., 0.)); - assert_eq!(metric.normalize_vec_at(rng.gen(), vec2(0., 1.)), vec2(0., 1. / 4.)); - assert_eq!(metric.normalize_vec_at(rng.gen(), vec2(1., 1.)), vec2(1. / 5., 1. / 5.)); + assert_eq!( + metric.normalize_vec_at(rng.gen(), vec2(1., 0.)), + vec2(1. / 3., 0.) + ); + assert_eq!( + metric.normalize_vec_at(rng.gen(), vec2(0., 1.)), + vec2(0., 1. / 4.) + ); + assert_eq!( + metric.normalize_vec_at(rng.gen(), vec2(1., 1.)), + vec2(1. / 5., 1. / 5.) + ); assert_eq!(metric.globalize(rng.gen(), vec2(1., 0.)), vec2(1. / 3., 0.)); assert_eq!(metric.globalize(rng.gen(), vec2(0., 1.)), vec2(0., 1. / 4.)); - assert_eq!(metric.globalize(rng.gen(), vec2(1., 1.)), vec2(1. / 3., 1. / 4.)); + assert_eq!( + metric.globalize(rng.gen(), vec2(1., 1.)), + vec2(1. / 3., 1. / 4.) + ); } #[test] fn test_trace_iter() { - let metric = samples::ScaledMetric { scale: vec2(2., 4.) }; - assert_eq!(trace_iter(&metric, vec2(3., 5.), vec2(1., 0.), 1.).nth(7).unwrap(), vec2(7., 5.)); - assert_eq!(trace_iter(&metric, vec2(3., 5.), vec2(2., 0.), 1.).nth(7).unwrap(), vec2(7., 5.)); - assert_eq!(trace_iter(&metric, vec2(3., 5.), vec2(1., 0.), 0.5).nth(7).unwrap(), vec2(5., 5.)); - assert_eq!(trace_iter(&metric, vec2(3., 5.), vec2(0., 1.), 1.).nth(9).unwrap(), vec2(3., 7.5)); - assert_eq!(trace_iter(&metric, vec2(3., 5.), vec2(0., 4.), 1.).nth(9).unwrap(), vec2(3., 7.5)); - assert_eq!(trace_iter(&metric, vec2(3., 5.), vec2(0., 1.), 0.5).nth(9).unwrap(), vec2(3., 6.25)); - assert_abs_diff_eq!(trace_iter(&metric, vec2(3., 5.), vec2(0.5, 0.25), std::f32::consts::SQRT_2).nth(7).unwrap(), vec2(7., 7.), epsilon=1e-5); + let metric = samples::ScaledMetric { + scale: vec2(2., 4.), + }; + assert_eq!( + trace_iter(&metric, vec2(3., 5.), vec2(1., 0.), 1.) + .nth(7) + .unwrap(), + vec2(7., 5.) + ); + assert_eq!( + trace_iter(&metric, vec2(3., 5.), vec2(2., 0.), 1.) + .nth(7) + .unwrap(), + vec2(7., 5.) + ); + assert_eq!( + trace_iter(&metric, vec2(3., 5.), vec2(1., 0.), 0.5) + .nth(7) + .unwrap(), + vec2(5., 5.) + ); + assert_eq!( + trace_iter(&metric, vec2(3., 5.), vec2(0., 1.), 1.) + .nth(9) + .unwrap(), + vec2(3., 7.5) + ); + assert_eq!( + trace_iter(&metric, vec2(3., 5.), vec2(0., 4.), 1.) + .nth(9) + .unwrap(), + vec2(3., 7.5) + ); + assert_eq!( + trace_iter(&metric, vec2(3., 5.), vec2(0., 1.), 0.5) + .nth(9) + .unwrap(), + vec2(3., 6.25) + ); + assert_abs_diff_eq!( + trace_iter( + &metric, + vec2(3., 5.), + vec2(0.5, 0.25), + std::f32::consts::SQRT_2 + ) + .nth(7) + .unwrap(), + vec2(7., 7.), + epsilon = 1e-5 + ); } } diff --git a/src/bin/flat/tube/coords.rs b/src/bin/flat/tube/coords.rs index 17fc479..edcc0f1 100644 --- a/src/bin/flat/tube/coords.rs +++ b/src/bin/flat/tube/coords.rs @@ -1,4 +1,4 @@ -use glam::{Mat2, Vec2, vec2}; +use glam::{vec2, Mat2, Vec2}; use crate::riemann::Metric; use crate::types::{Location, Ray}; @@ -10,15 +10,22 @@ pub trait FlatCoordinateSystem { fn global_to_flat(&self, v: T) -> T; } -pub trait FlatRegion: FlatCoordinateSystem + FlatCoordinateSystem + FlatCoordinateSystem { +pub trait FlatRegion: + FlatCoordinateSystem + FlatCoordinateSystem + FlatCoordinateSystem +{ // Измеряет расстояние до выхода за пределы области вдоль луча ray. Луч задаётся в плоской СК. - fn distance_to_boundary(&self, _ray: Ray) -> Option { None } + fn distance_to_boundary(&self, _ray: Ray) -> Option { + None + } } trait MetricCS: FlatCoordinateSystem { fn global_metric(&self) -> &impl Metric; fn flat_to_global_tfm_at(&self, pos: Vec2) -> Mat2 { - self.global_metric().sqrt_at(self.flat_to_global(pos)).inverse().into() + self.global_metric() + .sqrt_at(self.flat_to_global(pos)) + .inverse() + .into() } fn global_to_flat_tfm_at(&self, pos: Vec2) -> Mat2 { self.global_metric().sqrt_at(pos).into() @@ -59,7 +66,11 @@ impl + MetricCS> FlatCoordinateSystem fo pub struct InnerCS(pub Tube); -impl MetricCS for InnerCS { fn global_metric(&self) -> &impl Metric { &self.0 } } +impl MetricCS for InnerCS { + fn global_metric(&self) -> &impl Metric { + &self.0 + } +} impl FlatCoordinateSystem for InnerCS { fn flat_to_global(&self, pos: Vec2) -> Vec2 { @@ -74,20 +85,31 @@ impl FlatCoordinateSystem for InnerCS { impl FlatRegion for InnerCS { fn distance_to_boundary(&self, ray: Ray) -> Option { - Rect { size: vec2(self.0.inner_radius, self.0.internal_halflength) }.trace_out_of(ray) + Rect { + size: vec2(self.0.inner_radius, self.0.internal_halflength), + } + .trace_out_of(ray) } } pub struct OuterCS(pub Tube); -impl MetricCS for OuterCS { fn global_metric(&self) -> &impl Metric { &self.0 } } +impl MetricCS for OuterCS { + fn global_metric(&self) -> &impl Metric { + &self.0 + } +} impl FlatCoordinateSystem for OuterCS { fn flat_to_global(&self, pos: Vec2) -> Vec2 { - let inner = Rect { size: vec2(self.0.inner_radius + 1.0, self.0.external_halflength) }; + let inner = Rect { + size: vec2(self.0.inner_radius + 1.0, self.0.external_halflength), + }; if inner.is_inside(pos) { let Vec2 { x, y: v } = pos; - let y = self.0.y(v - v.signum() * (self.0.external_halflength - self.0.internal_halflength)); + let y = self + .0 + .y(v - v.signum() * (self.0.external_halflength - self.0.internal_halflength)); vec2(x, y) } else { pos @@ -95,10 +117,13 @@ impl FlatCoordinateSystem for OuterCS { } fn global_to_flat(&self, pos: Vec2) -> Vec2 { - let inner = Rect { size: vec2(self.0.inner_radius + 1.0, self.0.external_halflength) }; + let inner = Rect { + size: vec2(self.0.inner_radius + 1.0, self.0.external_halflength), + }; if inner.is_inside(pos) { let Vec2 { x: u, y } = pos; // в основной СК - let v = self.0.v(y) + y.signum() * (self.0.external_halflength - self.0.internal_halflength); + let v = self.0.v(y) + + y.signum() * (self.0.external_halflength - self.0.internal_halflength); vec2(u, v) // в плоском продолжении СК Outer на область Inner } else { pos @@ -108,14 +133,17 @@ impl FlatCoordinateSystem for OuterCS { impl FlatRegion for OuterCS { fn distance_to_boundary(&self, ray: Ray) -> Option { - Rect { size: vec2(self.0.outer_radius, self.0.external_halflength) }.trace_into(ray) + Rect { + size: vec2(self.0.outer_radius, self.0.external_halflength), + } + .trace_into(ray) } } #[cfg(test)] mod tests { - use approx::{AbsDiffEq, assert_abs_diff_eq}; - use glam::{Mat2, mat2, vec2, Vec2}; + use approx::{assert_abs_diff_eq, AbsDiffEq}; + use glam::{mat2, vec2, Mat2, Vec2}; use itertools_num::linspace; use crate::riemann::samples; @@ -126,22 +154,73 @@ mod tests { fn uniform_scaled_metric() { struct Scaled(samples::ScaledMetric, Vec2); impl FlatCoordinateSystem for Scaled { - fn flat_to_global(&self, pos: Vec2) -> Vec2 { (pos - self.1) / self.0.scale } - fn global_to_flat(&self, pos: Vec2) -> Vec2 { pos * self.0.scale + self.1 } + fn flat_to_global(&self, pos: Vec2) -> Vec2 { + (pos - self.1) / self.0.scale + } + fn global_to_flat(&self, pos: Vec2) -> Vec2 { + pos * self.0.scale + self.1 + } } impl MetricCS for Scaled { - fn global_metric(&self) -> &impl Metric { &self.0 } + fn global_metric(&self) -> &impl Metric { + &self.0 + } } - let cs = Scaled(samples::ScaledMetric { scale: vec2(3., 4.) }, vec2(2., 3.)); + let cs = Scaled( + samples::ScaledMetric { + scale: vec2(3., 4.), + }, + vec2(2., 3.), + ); assert_eq!(cs.global_to_flat(vec2(7., 3.)), vec2(23., 15.)); assert_eq!(cs.flat_to_global(vec2(8., 7.)), vec2(2., 1.)); - assert_eq!(cs.global_to_flat(Ray { pos: vec2(7., 3.), dir: vec2(3., 2.) }), Ray { pos: vec2(23., 15.), dir: vec2(9., 8.) }); - assert_eq!(cs.flat_to_global(Ray { pos: vec2(23., 15.), dir: vec2(9., 8.) }), Ray { pos: vec2(7., 3.), dir: vec2(3., 2.) }); - assert_eq!(cs.global_to_flat(Location { pos: vec2(2., 1.), rot: mat2(vec2(0., 1.), vec2(-1., 0.)) }), Location { pos: vec2(8., 7.), rot: mat2(vec2(0., 4.), vec2(-3., 0.)) }); - assert_eq!(cs.flat_to_global(Location { pos: vec2(2., 1.), rot: mat2(vec2(0., 1.), vec2(-1., 0.)) }), Location { pos: vec2(0., -0.5), rot: mat2(vec2(0., 0.25), vec2(-1. / 3., 0.)) }); + assert_eq!( + cs.global_to_flat(Ray { + pos: vec2(7., 3.), + dir: vec2(3., 2.) + }), + Ray { + pos: vec2(23., 15.), + dir: vec2(9., 8.) + } + ); + assert_eq!( + cs.flat_to_global(Ray { + pos: vec2(23., 15.), + dir: vec2(9., 8.) + }), + Ray { + pos: vec2(7., 3.), + dir: vec2(3., 2.) + } + ); + assert_eq!( + cs.global_to_flat(Location { + pos: vec2(2., 1.), + rot: mat2(vec2(0., 1.), vec2(-1., 0.)) + }), + Location { + pos: vec2(8., 7.), + rot: mat2(vec2(0., 4.), vec2(-3., 0.)) + } + ); + assert_eq!( + cs.flat_to_global(Location { + pos: vec2(2., 1.), + rot: mat2(vec2(0., 1.), vec2(-1., 0.)) + }), + Location { + pos: vec2(0., -0.5), + rot: mat2(vec2(0., 0.25), vec2(-1. / 3., 0.)) + } + ); } - fn test_flat_region(region: &impl FlatRegion, range_global: (Vec2, Vec2), range_flat: (Vec2, Vec2)) { + fn test_flat_region( + region: &impl FlatRegion, + range_global: (Vec2, Vec2), + range_flat: (Vec2, Vec2), + ) { #[allow(non_upper_case_globals)] const ε: f32 = 1e-3; macro_rules! assert_eq_at { @@ -149,35 +228,105 @@ mod tests { let at = $at; let left = $left; let right = $right; - assert!(left.abs_diff_eq(right, ε), "Assertion failed at {at}:\n left: {left} = {}\n right: {right} = {}", stringify!($left), stringify!($right)); + assert!( + left.abs_diff_eq(right, ε), + "Assertion failed at {at}:\n left: {left} = {}\n right: {right} = {}", + stringify!($left), + stringify!($right) + ); }; } - fn check_range(name_a: &str, a: Vec2, range_a: (Vec2, Vec2), name_b: &str, b: Vec2, range_b: (Vec2, Vec2)) { + fn check_range( + name_a: &str, + a: Vec2, + range_a: (Vec2, Vec2), + name_b: &str, + b: Vec2, + range_b: (Vec2, Vec2), + ) { assert!(b.cmpge(range_b.0 - ε).all() && b.cmple(range_b.1 + ε).all(), "Assertion failed:\nAt {name_a}: {a}, from range: {range_a:?}\nGot {name_b}: {b}, which is out of range {range_b:?}"); // TODO sort out when to check these conditions: - if a.x.abs_diff_eq(&range_a.0.x, ε) { assert_abs_diff_eq!(b.x, range_b.0.x, epsilon=ε); } - if a.y.abs_diff_eq(&range_a.0.y, ε) { assert_abs_diff_eq!(b.y, range_b.0.y, epsilon=ε); } - if a.x.abs_diff_eq(&range_a.1.x, ε) { assert_abs_diff_eq!(b.x, range_b.1.x, epsilon=ε); } - if a.y.abs_diff_eq(&range_a.1.y, ε) { assert_abs_diff_eq!(b.y, range_b.1.y, epsilon=ε); } + if a.x.abs_diff_eq(&range_a.0.x, ε) { + assert_abs_diff_eq!(b.x, range_b.0.x, epsilon = ε); + } + if a.y.abs_diff_eq(&range_a.0.y, ε) { + assert_abs_diff_eq!(b.y, range_b.0.y, epsilon = ε); + } + if a.x.abs_diff_eq(&range_a.1.x, ε) { + assert_abs_diff_eq!(b.x, range_b.1.x, epsilon = ε); + } + if a.y.abs_diff_eq(&range_a.1.y, ε) { + assert_abs_diff_eq!(b.y, range_b.1.y, epsilon = ε); + } } for x in linspace(range_global.0.x, range_global.1.x, 20) { for y in linspace(range_global.0.y, range_global.1.y, 20) { let pos_global = vec2(x, y); let pos_flat = region.global_to_flat(pos_global); - check_range("global", pos_global, range_global, "flat", pos_flat, range_flat); - assert_eq_at!(pos_global, region.global_to_flat(Location { pos: pos_global, rot: Mat2::IDENTITY }).pos, pos_flat); + check_range( + "global", + pos_global, + range_global, + "flat", + pos_flat, + range_flat, + ); + assert_eq_at!( + pos_global, + region + .global_to_flat(Location { + pos: pos_global, + rot: Mat2::IDENTITY + }) + .pos, + pos_flat + ); assert_eq_at!(pos_global, region.flat_to_global(pos_flat), pos_global); - assert_eq_at!(pos_global, region.flat_to_global(region.global_to_flat(Location { pos: pos_global, rot: Mat2::IDENTITY })).rot, Mat2::IDENTITY); + assert_eq_at!( + pos_global, + region + .flat_to_global(region.global_to_flat(Location { + pos: pos_global, + rot: Mat2::IDENTITY + })) + .rot, + Mat2::IDENTITY + ); } } for x in linspace(range_flat.0.x, range_flat.1.x, 20) { for y in linspace(range_flat.0.y, range_flat.1.y, 20) { let pos_flat = vec2(x, y); let pos_global = region.flat_to_global(pos_flat); - check_range("flat", pos_flat, range_flat, "global", pos_global, range_global); - assert_eq_at!(pos_flat, region.flat_to_global(Location { pos: pos_flat, rot: Mat2::IDENTITY }).pos, pos_global); + check_range( + "flat", + pos_flat, + range_flat, + "global", + pos_global, + range_global, + ); + assert_eq_at!( + pos_flat, + region + .flat_to_global(Location { + pos: pos_flat, + rot: Mat2::IDENTITY + }) + .pos, + pos_global + ); assert_eq_at!(pos_flat, region.global_to_flat(pos_global), pos_flat); - assert_eq_at!(pos_flat, region.global_to_flat(region.flat_to_global(Location { pos: pos_global, rot: Mat2::IDENTITY })).rot, Mat2::IDENTITY); + assert_eq_at!( + pos_flat, + region + .global_to_flat(region.flat_to_global(Location { + pos: pos_global, + rot: Mat2::IDENTITY + })) + .rot, + Mat2::IDENTITY + ); } } } @@ -190,9 +339,21 @@ mod tests { internal_halflength: 100.0, external_halflength: 300.0, }); - test_flat_region(&mapper, (vec2(-30.0, -300.0), vec2(30.0, 300.0)), (vec2(-30.0, -100.0), vec2(30.0, 100.0))); - test_flat_region(&mapper, (vec2(-60.0, -400.0), vec2(60.0, -300.0)), (vec2(-60.0, -200.0), vec2(60.0, -100.0))); - test_flat_region(&mapper, (vec2(-60.0, 300.0), vec2(60.0, 400.0)), (vec2(-60.0, 100.0), vec2(60.0, 200.0))); + test_flat_region( + &mapper, + (vec2(-30.0, -300.0), vec2(30.0, 300.0)), + (vec2(-30.0, -100.0), vec2(30.0, 100.0)), + ); + test_flat_region( + &mapper, + (vec2(-60.0, -400.0), vec2(60.0, -300.0)), + (vec2(-60.0, -200.0), vec2(60.0, -100.0)), + ); + test_flat_region( + &mapper, + (vec2(-60.0, 300.0), vec2(60.0, 400.0)), + (vec2(-60.0, 100.0), vec2(60.0, 200.0)), + ); } #[test] @@ -204,23 +365,68 @@ mod tests { external_halflength: 300.0, }); // TODO replace 200.20016 with something sane - test_flat_region(&mapper, (vec2(-30.0, -300.0), vec2(30.0, -1.0)), (vec2(-30.0, -300.0), vec2(30.0, -200.20016))); - test_flat_region(&mapper, (vec2(-30.0, 1.0), vec2(30.0, 300.0)), (vec2(-30.0, 200.20016), vec2(30.0, 300.0))); - test_flat_region(&mapper, (vec2(-60.0, -400.0), vec2(60.0, -300.0)), (vec2(-60.0, -400.0), vec2(60.0, -300.0))); - test_flat_region(&mapper, (vec2(-60.0, 300.0), vec2(60.0, 400.0)), (vec2(-60.0, 300.0), vec2(60.0, 400.0))); + test_flat_region( + &mapper, + (vec2(-30.0, -300.0), vec2(30.0, -1.0)), + (vec2(-30.0, -300.0), vec2(30.0, -200.20016)), + ); + test_flat_region( + &mapper, + (vec2(-30.0, 1.0), vec2(30.0, 300.0)), + (vec2(-30.0, 200.20016), vec2(30.0, 300.0)), + ); + test_flat_region( + &mapper, + (vec2(-60.0, -400.0), vec2(60.0, -300.0)), + (vec2(-60.0, -400.0), vec2(60.0, -300.0)), + ); + test_flat_region( + &mapper, + (vec2(-60.0, 300.0), vec2(60.0, 400.0)), + (vec2(-60.0, 300.0), vec2(60.0, 400.0)), + ); // straight for x in linspace(-60., 60., 20) { for y in linspace(-320., 320., 20) { - assert_eq!(mapper.global_to_flat(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos.x, x); + assert_eq!( + mapper + .global_to_flat(Location { + pos: vec2(x, y), + rot: Mat2::IDENTITY + }) + .pos + .x, + x + ); } } // symmetrical for x in linspace(0., 60., 20) { for y in linspace(0., 320., 20) { - let pp = mapper.global_to_flat(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos; - let np = mapper.global_to_flat(Location { pos: vec2(-x, y), rot: Mat2::IDENTITY }).pos; - let pn = mapper.global_to_flat(Location { pos: vec2(x, -y), rot: Mat2::IDENTITY }).pos; - let nn = mapper.global_to_flat(Location { pos: vec2(-x, -y), rot: Mat2::IDENTITY }).pos; + let pp = mapper + .global_to_flat(Location { + pos: vec2(x, y), + rot: Mat2::IDENTITY, + }) + .pos; + let np = mapper + .global_to_flat(Location { + pos: vec2(-x, y), + rot: Mat2::IDENTITY, + }) + .pos; + let pn = mapper + .global_to_flat(Location { + pos: vec2(x, -y), + rot: Mat2::IDENTITY, + }) + .pos; + let nn = mapper + .global_to_flat(Location { + pos: vec2(-x, -y), + rot: Mat2::IDENTITY, + }) + .pos; assert_eq!(np, vec2(-pp.x, pp.y)); assert_eq!(pn, vec2(pp.x, -pp.y)); assert_eq!(nn, vec2(-pp.x, -pp.y)); @@ -229,18 +435,42 @@ mod tests { // clean boundary for x in linspace(50., 60., 20) { for y in linspace(0., 320., 20) { - assert_eq!(mapper.global_to_flat(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos.y, y); + assert_eq!( + mapper + .global_to_flat(Location { + pos: vec2(x, y), + rot: Mat2::IDENTITY + }) + .pos + .y, + y + ); } } for x in linspace(0., 60., 20) { for y in linspace(300., 320., 20) { - assert_eq!(mapper.global_to_flat(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos.y, y); + assert_eq!( + mapper + .global_to_flat(Location { + pos: vec2(x, y), + rot: Mat2::IDENTITY + }) + .pos + .y, + y + ); } } // accelerating for x in linspace(-29., 29., 20) { for y in linspace(1., 299., 20) { - let v = mapper.global_to_flat(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos.y; + let v = mapper + .global_to_flat(Location { + pos: vec2(x, y), + rot: Mat2::IDENTITY, + }) + .pos + .y; assert!(v > 200.0); assert!(v > y); } diff --git a/src/bin/flat/tube/metric.rs b/src/bin/flat/tube/metric.rs index 308a8ee..f88b69b 100644 --- a/src/bin/flat/tube/metric.rs +++ b/src/bin/flat/tube/metric.rs @@ -1,4 +1,4 @@ -use glam::{f32, Mat2, Vec2, vec2}; +use glam::{f32, vec2, Mat2, Vec2}; use crate::fns::{self, Limiter}; use crate::riemann::{Decomp2, Metric, Tens2}; @@ -12,13 +12,31 @@ pub struct Tube { } impl Tube { - fn fx(&self) -> impl Limiter { fns::SmootherstepLimiter { min: self.inner_radius, max: self.outer_radius } } - fn fy(&self) -> fns::QuadraticAccelerator { fns::QuadraticAccelerator { internal: self.internal_halflength, external: self.external_halflength } } + fn fx(&self) -> impl Limiter { + fns::SmootherstepLimiter { + min: self.inner_radius, + max: self.outer_radius, + } + } + fn fy(&self) -> fns::QuadraticAccelerator { + fns::QuadraticAccelerator { + internal: self.internal_halflength, + external: self.external_halflength, + } + } - pub fn y(&self, v: f32) -> f32 { self.fy().x(v) } - pub fn v(&self, y: f32) -> f32 { self.fy().u(y) } - pub fn dy(&self, v: f32) -> f32 { self.fy().dx(v) } - pub fn dv(&self, y: f32) -> f32 { self.fy().du(y) } + pub fn y(&self, v: f32) -> f32 { + self.fy().x(v) + } + pub fn v(&self, y: f32) -> f32 { + self.fy().u(y) + } + pub fn dy(&self, v: f32) -> f32 { + self.fy().dx(v) + } + pub fn dv(&self, y: f32) -> f32 { + self.fy().du(y) + } } impl Metric for Tube { @@ -53,7 +71,7 @@ impl Metric for Tube { #[cfg(test)] mod test { use approx::assert_abs_diff_eq; - use glam::{Vec2, vec2}; + use glam::{vec2, Vec2}; use itertools_num::linspace; use crate::riemann::{Decomp2, Metric}; @@ -66,7 +84,9 @@ mod test { fn test_tube_metric_derivs() { struct Approx(Tube); impl Metric for Approx { - fn sqrt_at(&self, pos: Vec2) -> Decomp2 { self.0.sqrt_at(pos) } + fn sqrt_at(&self, pos: Vec2) -> Decomp2 { + self.0.sqrt_at(pos) + } } let testee = Tube { inner_radius: 30.0, @@ -78,8 +98,13 @@ mod test { let epsilon = 1.0e-3; let margin = 1.0 / 16.0; let mul = 1.0 + margin; - for x in itertools_num::linspace(-mul * testee.outer_radius, mul * testee.outer_radius, 100) { - for y in itertools_num::linspace(-mul * testee.external_halflength, mul * testee.external_halflength, 100) { + for x in itertools_num::linspace(-mul * testee.outer_radius, mul * testee.outer_radius, 100) + { + for y in itertools_num::linspace( + -mul * testee.external_halflength, + mul * testee.external_halflength, + 100, + ) { let pos = vec2(x, y); let computed = testee.part_derivs_at(pos); let reference = approx.part_derivs_at(pos); @@ -110,10 +135,10 @@ mod test { let Δ = vec2(bx - ax, 2.0 * (space.tube.internal_halflength + off)); let dir = Δ / (steps as f32); let traced = space.trace_iter(Ray { pos: a, dir }).nth(steps).unwrap(); - assert_abs_diff_eq!(traced.pos.x, b.x, epsilon=1.0e-2); - assert_abs_diff_eq!(traced.pos.y, b.y, epsilon=1.0e1); - assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon=1.0e-3); - assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon=1.0e-2); + assert_abs_diff_eq!(traced.pos.x, b.x, epsilon = 1.0e-2); + assert_abs_diff_eq!(traced.pos.y, b.y, epsilon = 1.0e1); + assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon = 1.0e-3); + assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon = 1.0e-2); } } } @@ -133,17 +158,25 @@ mod test { let ε = 1e-3; let off = 10.0; let steps = 4096; - for ax in linspace(-space.tube.inner_radius + ε, space.tube.inner_radius - ε, 20) { - for bx in linspace(-space.tube.inner_radius + ε, space.tube.inner_radius - ε, 20) { + for ax in linspace( + -space.tube.inner_radius + ε, + space.tube.inner_radius - ε, + 20, + ) { + for bx in linspace( + -space.tube.inner_radius + ε, + space.tube.inner_radius - ε, + 20, + ) { let a = vec2(ax, -(space.tube.external_halflength + off)); let b = vec2(bx, space.tube.external_halflength + off); let Δ = vec2(bx - ax, 2.0 * (space.tube.internal_halflength + off)); let dir = Δ / (steps as f32); let traced = space.trace_iter(Ray { pos: a, dir }).nth(steps).unwrap(); - assert_abs_diff_eq!(traced.pos.x, b.x, epsilon=1.0e-2); - assert_abs_diff_eq!(traced.pos.y, b.y, epsilon=1.0e0); - assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon=1.0e-3); - assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon=1.0e-3); + assert_abs_diff_eq!(traced.pos.x, b.x, epsilon = 1.0e-2); + assert_abs_diff_eq!(traced.pos.y, b.y, epsilon = 1.0e0); + assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon = 1.0e-3); + assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon = 1.0e-3); } } } @@ -168,10 +201,10 @@ mod test { let Δ = vec2(0.0, 2.0 * (space.tube.internal_halflength + off)); let dir = Δ / (steps as f32); let traced = space.trace_iter(Ray { pos: a, dir }).nth(steps).unwrap(); - assert_abs_diff_eq!(traced.pos.x, b.x, epsilon=1.0e-1); - assert_abs_diff_eq!(traced.pos.y, b.y, epsilon=1.0e0); - assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon=1.0e-2); - assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon=1.0e-2); + assert_abs_diff_eq!(traced.pos.x, b.x, epsilon = 1.0e-1); + assert_abs_diff_eq!(traced.pos.y, b.y, epsilon = 1.0e0); + assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon = 1.0e-2); + assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon = 1.0e-2); } } @@ -195,10 +228,10 @@ mod test { let Δ = vec2(0.0, 2.0 * (space.tube.external_halflength + off)); let dir = Δ / (steps as f32); let traced = space.trace_iter(Ray { pos: a, dir }).nth(steps).unwrap(); - assert_abs_diff_eq!(traced.pos.x, b.x, epsilon=2.0e0); - assert_abs_diff_eq!(traced.pos.y, b.y, epsilon=1.0e0); - assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon=1.0e-2); - assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon=1.0e-2); + assert_abs_diff_eq!(traced.pos.x, b.x, epsilon = 2.0e0); + assert_abs_diff_eq!(traced.pos.y, b.y, epsilon = 1.0e0); + assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon = 1.0e-2); + assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon = 1.0e-2); } } } diff --git a/src/bin/flat/tube/mod.rs b/src/bin/flat/tube/mod.rs index 722f98b..29bb9b5 100644 --- a/src/bin/flat/tube/mod.rs +++ b/src/bin/flat/tube/mod.rs @@ -1,4 +1,4 @@ -use glam::{bool, f32, Mat2, Vec2, vec2}; +use glam::{bool, f32, vec2, Mat2, Vec2}; use coords::{FlatCoordinateSystem, InnerCS, OuterCS}; use metric::Tube; @@ -8,8 +8,8 @@ use crate::riemann; use crate::tube::coords::FlatRegion; use crate::types::{FlatTraceResult, Hit, Location, Object, Ray}; -pub mod metric; mod coords; +pub mod metric; pub struct Space { pub tube: Tube, @@ -48,12 +48,16 @@ impl Space { /// Выполняет один шаг перемещения. Работает в любой части пространства. /// off задаётся в локальной СК. Рекомендуется считать небольшими шагами. pub fn move_step(&self, loc: Location, off: Vec2) -> Location { - let corr = Mat2::IDENTITY - riemann::contract(riemann::krist(&self.tube, loc.pos), loc.rot * off); + let corr = + Mat2::IDENTITY - riemann::contract(riemann::krist(&self.tube, loc.pos), loc.rot * off); let p = loc.pos + corr * loc.rot * off; - Location { pos: p, rot: corr * loc.rot } + Location { + pos: p, + rot: corr * loc.rot, + } } - pub fn trace_iter(&self, ray: Ray) -> impl Iterator + '_ { + pub fn trace_iter(&self, ray: Ray) -> impl Iterator + '_ { std::iter::successors(Some(ray), |&ray| Some(self.trace_step(ray))) } @@ -85,15 +89,28 @@ impl Space { } fn list_objects(&self, tfm: impl Fn(Location) -> Location) -> Vec { - self.objs.iter().map(|&Object { id, loc, r }| Object { id, loc: tfm(loc), r }).collect() + self.objs + .iter() + .map(|&Object { id, loc, r }| Object { + id, + loc: tfm(loc), + r, + }) + .collect() } - fn hit_objects(objs: &[Object], ray: Ray, limit: Option, globalize: impl Fn(Vec2) -> Vec2) -> Vec { + fn hit_objects( + objs: &[Object], + ray: Ray, + limit: Option, + globalize: impl Fn(Vec2) -> Vec2, + ) -> Vec { let limit = limit.unwrap_or(f32::INFINITY); objs.iter() .filter_map(|obj| { let rel = ray.pos - obj.loc.pos; - let diff = rel.dot(ray.dir).powi(2) - ray.dir.length_squared() * (rel.length_squared() - obj.r.powi(2)); + let diff = rel.dot(ray.dir).powi(2) + - ray.dir.length_squared() * (rel.length_squared() - obj.r.powi(2)); if diff > 0.0 { let t = (-rel.dot(ray.dir) - diff.sqrt()) / ray.dir.length_squared(); Some((obj, t)) @@ -104,8 +121,17 @@ impl Space { .filter(|&(_, t)| t >= 0.0 && t < limit) .map(|(obj, t)| { let pos = ray.forward(t).pos; - let rel = obj.loc.rot.inverse() * Ray { pos: pos - obj.loc.pos, dir: ray.dir }; - Hit { id: obj.id, distance: t, pos: globalize(pos), rel } + let rel = obj.loc.rot.inverse() + * Ray { + pos: pos - obj.loc.pos, + dir: ray.dir, + }; + Hit { + id: obj.id, + distance: t, + pos: globalize(pos), + rel, + } }) .collect() } @@ -118,7 +144,9 @@ impl Space { let n = ((b - a).length() / step) as usize + 1; let a = cs.global_to_flat(a); let b = cs.global_to_flat(b); - (1..=n).map(|k| cs.flat_to_global(a.lerp(b, k as f32 / n as f32))).collect() + (1..=n) + .map(|k| cs.flat_to_global(a.lerp(b, k as f32 / n as f32))) + .collect() } Boundary => panic!("Can't draw a line here!"), } @@ -132,7 +160,10 @@ struct Rect { impl Rect { /// Отражает луч, чтобы все координаты направления были положительны (допустимо благодаря симметрии Rect). fn flip_ray(ray: Ray) -> Ray { - Ray { pos: ray.pos * ray.dir.signum(), dir: ray.dir.abs() } + Ray { + pos: ray.pos * ray.dir.signum(), + dir: ray.dir.abs(), + } } fn is_inside(&self, pt: Vec2) -> bool { @@ -145,8 +176,12 @@ impl Rect { let ts = (-self.size - ray.pos) / ray.dir; let t = ts.max_element(); let pt = ray.pos + t * ray.dir; - if t < 0.0 { return None; } - if pt.cmpgt(self.size).any() { return None; } + if t < 0.0 { + return None; + } + if pt.cmpgt(self.size).any() { + return None; + } Some(t) } @@ -161,27 +196,149 @@ impl Rect { #[test] fn test_rect() { - assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 5.0) }), Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 5.0) }); - assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(-4.0, 5.0) }), Ray { pos: vec2(-2.0, 3.0), dir: vec2(4.0, 5.0) }); - assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, -5.0) }), Ray { pos: vec2(2.0, -3.0), dir: vec2(4.0, 5.0) }); - assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 0.0) }), Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 0.0) }); + assert_eq!( + Rect::flip_ray(Ray { + pos: vec2(2.0, 3.0), + dir: vec2(4.0, 5.0) + }), + Ray { + pos: vec2(2.0, 3.0), + dir: vec2(4.0, 5.0) + } + ); + assert_eq!( + Rect::flip_ray(Ray { + pos: vec2(2.0, 3.0), + dir: vec2(-4.0, 5.0) + }), + Ray { + pos: vec2(-2.0, 3.0), + dir: vec2(4.0, 5.0) + } + ); + assert_eq!( + Rect::flip_ray(Ray { + pos: vec2(2.0, 3.0), + dir: vec2(4.0, -5.0) + }), + Ray { + pos: vec2(2.0, -3.0), + dir: vec2(4.0, 5.0) + } + ); + assert_eq!( + Rect::flip_ray(Ray { + pos: vec2(2.0, 3.0), + dir: vec2(4.0, 0.0) + }), + Ray { + pos: vec2(2.0, 3.0), + dir: vec2(4.0, 0.0) + } + ); - let r = Rect { size: vec2(2.0, 3.0) }; + let r = Rect { + size: vec2(2.0, 3.0), + }; - assert_eq!(r.trace_into(Ray { pos: vec2(3.0, 3.0), dir: vec2(1.0, 1.0) }), None); - assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 2.0), dir: vec2(1.0, 0.0) }), Some(1.0)); - assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 2.0), dir: vec2(-1.0, 0.0) }), None); - assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 1.0), dir: vec2(2.0, 2.0) }), Some(0.5)); - assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 2.1), dir: vec2(2.0, 2.0) }), None); + assert_eq!( + r.trace_into(Ray { + pos: vec2(3.0, 3.0), + dir: vec2(1.0, 1.0) + }), + None + ); + assert_eq!( + r.trace_into(Ray { + pos: vec2(-3.0, 2.0), + dir: vec2(1.0, 0.0) + }), + Some(1.0) + ); + assert_eq!( + r.trace_into(Ray { + pos: vec2(-3.0, 2.0), + dir: vec2(-1.0, 0.0) + }), + None + ); + assert_eq!( + r.trace_into(Ray { + pos: vec2(-3.0, 1.0), + dir: vec2(2.0, 2.0) + }), + Some(0.5) + ); + assert_eq!( + r.trace_into(Ray { + pos: vec2(-3.0, 2.1), + dir: vec2(2.0, 2.0) + }), + None + ); - assert_eq!(r.trace_into(Ray { pos: vec2(2.0, 3.0), dir: vec2(1.0, 1.0) }), None); - assert_eq!(r.trace_into(Ray { pos: vec2(-2.0, 3.0), dir: vec2(-1.0, 1.0) }), None); - assert_eq!(r.trace_into(Ray { pos: vec2(2.0, 3.0), dir: vec2(-1.0, -1.0) }), Some(0.0)); - assert_eq!(r.trace_into(Ray { pos: vec2(2.0, -3.0), dir: vec2(-1.0, 1.0) }), Some(0.0)); + assert_eq!( + r.trace_into(Ray { + pos: vec2(2.0, 3.0), + dir: vec2(1.0, 1.0) + }), + None + ); + assert_eq!( + r.trace_into(Ray { + pos: vec2(-2.0, 3.0), + dir: vec2(-1.0, 1.0) + }), + None + ); + assert_eq!( + r.trace_into(Ray { + pos: vec2(2.0, 3.0), + dir: vec2(-1.0, -1.0) + }), + Some(0.0) + ); + assert_eq!( + r.trace_into(Ray { + pos: vec2(2.0, -3.0), + dir: vec2(-1.0, 1.0) + }), + Some(0.0) + ); - assert_eq!(r.trace_out_of(Ray { pos: vec2(0.0, 0.0), dir: vec2(1.0, 1.0) }), Some(2.0)); - assert_eq!(r.trace_out_of(Ray { pos: vec2(0.0, 0.0), dir: vec2(0.0, 1.0) }), Some(3.0)); - assert_eq!(r.trace_out_of(Ray { pos: vec2(0.0, 1.0), dir: vec2(0.0, -1.0) }), Some(4.0)); - assert_eq!(r.trace_out_of(Ray { pos: vec2(1.0, 1.0), dir: vec2(0.0, -1.0) }), Some(4.0)); - assert_eq!(r.trace_out_of(Ray { pos: vec2(2.0, 3.0), dir: vec2(1.0, 1.0) }), Some(0.0)); + assert_eq!( + r.trace_out_of(Ray { + pos: vec2(0.0, 0.0), + dir: vec2(1.0, 1.0) + }), + Some(2.0) + ); + assert_eq!( + r.trace_out_of(Ray { + pos: vec2(0.0, 0.0), + dir: vec2(0.0, 1.0) + }), + Some(3.0) + ); + assert_eq!( + r.trace_out_of(Ray { + pos: vec2(0.0, 1.0), + dir: vec2(0.0, -1.0) + }), + Some(4.0) + ); + assert_eq!( + r.trace_out_of(Ray { + pos: vec2(1.0, 1.0), + dir: vec2(0.0, -1.0) + }), + Some(4.0) + ); + assert_eq!( + r.trace_out_of(Ray { + pos: vec2(2.0, 3.0), + dir: vec2(1.0, 1.0) + }), + Some(0.0) + ); } diff --git a/src/bin/flat/types.rs b/src/bin/flat/types.rs index 2edc6d8..72e88dc 100644 --- a/src/bin/flat/types.rs +++ b/src/bin/flat/types.rs @@ -8,7 +8,10 @@ pub struct Ray { impl Ray { pub fn forward(&self, dist: f32) -> Ray { - Ray { pos: self.pos + self.dir * dist, dir: self.dir } + Ray { + pos: self.pos + self.dir * dist, + dir: self.dir, + } } } @@ -16,7 +19,10 @@ impl std::ops::Mul for Mat2 { type Output = Ray; fn mul(self, rhs: Ray) -> Self::Output { - Ray { pos: self * rhs.pos, dir: self * rhs.dir } + Ray { + pos: self * rhs.pos, + dir: self * rhs.dir, + } } } @@ -40,7 +46,7 @@ pub struct Hit { pub distance: f32, pub id: i32, pub pos: Vec2, // положение в основной СК - pub rel: Ray, // в локальной ортонормированной СК объекта + pub rel: Ray, // в локальной ортонормированной СК объекта } #[derive(Clone, Debug)] diff --git a/src/bin/mesh/main.rs b/src/bin/mesh/main.rs index 04a44d6..1a4a1a0 100644 --- a/src/bin/mesh/main.rs +++ b/src/bin/mesh/main.rs @@ -1,13 +1,13 @@ mod mesh_loader; -use std::fs::File; -use std::{env}; -use std::error::Error; -use std::f32::consts::PI; -use std::io::{BufReader}; +use crate::mesh_loader::{load_mesh, Face}; use glam::*; use show_image::{ImageInfo, ImageView, WindowOptions}; -use crate::mesh_loader::{Face, load_mesh}; +use std::env; +use std::error::Error; +use std::f32::consts::PI; +use std::fs::File; +use std::io::BufReader; const W: i32 = 320; const H: i32 = 240; @@ -38,19 +38,26 @@ impl Image { } fn ypr_to_mat(ypr: Vec3) -> Mat3 { - let Vec3 { x: yaw, y: pitch, z: roll } = ypr; + let Vec3 { + x: yaw, + y: pitch, + z: roll, + } = ypr; let m_roll = mat3( vec3(roll.cos(), roll.sin(), 0.0), vec3(-roll.sin(), roll.cos(), 0.0), - vec3(0.0, 0.0, 1.0)); + vec3(0.0, 0.0, 1.0), + ); let m_yaw = mat3( vec3(yaw.cos(), 0.0, yaw.sin()), vec3(0.0, 1.0, 0.0), - vec3(-yaw.sin(), 0.0, yaw.cos())); + vec3(-yaw.sin(), 0.0, yaw.cos()), + ); let m_pitch = mat3( vec3(1.0, 0.0, 0.0), vec3(0.0, pitch.cos(), -pitch.sin()), - vec3(0.0, pitch.sin(), pitch.cos())); + vec3(0.0, pitch.sin(), pitch.cos()), + ); m_roll * m_pitch * m_yaw } @@ -67,7 +74,11 @@ fn trace_to_mesh(mesh: &Mesh, base: Vec3, ray: Vec3) -> Option { for f in mesh { let fs = (0..3).map(|k| edge_dist(f.vertices[k], f.vertices[(k + 1) % 3], base, ray)); if fs.into_iter().all(|f| f >= 0.0) { - let m = mat3(f.vertices[1] - f.vertices[0], f.vertices[2] - f.vertices[0], -ray); + let m = mat3( + f.vertices[1] - f.vertices[0], + f.vertices[2] - f.vertices[0], + -ray, + ); let m = m.inverse(); let rel = m * (base - f.vertices[0]); if rel.z > dist { @@ -107,7 +118,9 @@ fn render(mesh: &Mesh, camera: impl Fn(Vec2) -> (Vec3, Vec3)) -> Image { } else { bkg }; - let color = (color * 255.0).as_ivec3().clamp(IVec3::splat(0), IVec3::splat(255)); + let color = (color * 255.0) + .as_ivec3() + .clamp(IVec3::splat(0), IVec3::splat(255)); img.put_pixel(x, y, Color(color.x as u8, color.y as u8, color.z as u8)); } } @@ -125,7 +138,11 @@ fn main() -> Result<(), Box> { let window = show_image::create_window("Raytracing", WindowOptions::default())?; loop { for phi in 0..360 { - let m_view = ypr_to_mat(vec3((135.0 + phi as f32) * PI / 180.0, -30.0 * PI / 180.0, 0.0f32)); + let m_view = ypr_to_mat(vec3( + (135.0 + phi as f32) * PI / 180.0, + -30.0 * PI / 180.0, + 0.0f32, + )); let m_camera = m_view.transpose(); let img = render(mesh.as_slice(), |off| { // perspective projection diff --git a/src/bin/mesh/mesh_loader.rs b/src/bin/mesh/mesh_loader.rs index 57fc99b..f9d5a8f 100644 --- a/src/bin/mesh/mesh_loader.rs +++ b/src/bin/mesh/mesh_loader.rs @@ -1,5 +1,5 @@ -use std::io; use glam::{vec2, vec3, Vec2, Vec3}; +use std::io; #[derive(Copy, Clone, Debug)] struct ObjVertex { @@ -35,14 +35,23 @@ impl ObjMesh { } fn parse_fv(desc: &&str) -> ObjVertex { - let tokens: Vec<_> = desc.split('/').map(|s| s.parse::().unwrap() - 1).collect(); + let tokens: Vec<_> = desc + .split('/') + .map(|s| s.parse::().unwrap() - 1) + .collect(); assert_eq!(tokens.len(), 3); - ObjVertex { vertex: tokens[0], tex_coord: tokens[1], normal: tokens[2] } + ObjVertex { + vertex: tokens[0], + tex_coord: tokens[1], + normal: tokens[2], + } } fn parse_f(tokens: &[&str]) -> ObjFace { let vertices: Vec<_> = tokens.iter().map(ObjMesh::parse_fv).collect(); - ObjFace { vertices: vertices.as_slice().try_into().unwrap() } + ObjFace { + vertices: vertices.as_slice().try_into().unwrap(), + } } fn read(f: &mut impl io::BufRead) -> io::Result { @@ -57,13 +66,7 @@ impl ObjMesh { if f.read_line(&mut line)? == 0 { break; } - let tokens: Vec<&str> = line - .trim() - .split('#') - .next() - .unwrap() - .split(' ') - .collect(); + let tokens: Vec<&str> = line.trim().split('#').next().unwrap().split(' ').collect(); match tokens[0] { "v" => result.vertices.push(Self::parse_v3(&tokens[1..])), "vn" => result.normals.push(Self::parse_v3(&tokens[1..])), @@ -76,12 +79,13 @@ impl ObjMesh { } fn flatten(&self) -> Vec { - self.faces.iter().map(|face| { - Face { + self.faces + .iter() + .map(|face| Face { vertices: face.vertices.map(|iv| self.vertices[iv.vertex]), normal: self.normals[face.vertices[0].normal], - } - }).collect() + }) + .collect() } }