Cargo.toml

@@ -11,9 +11,6 @@ panic = 'abort'
 [profile.dev.package."*"]
 opt-level = 3
 
-[profile.test.package."*"]
-opt-level = 3
-
 [dependencies]
 rand = "0.8.5"
 glam = { version = "0.27.0", features = ["approx", "fast-math", "rand"] }

rustfmt.toml

@@ -1 +0,0 @@
-hard_tabs = true

src/bin/flat/float_fun.rs

@@ -12,12 +12,8 @@ pub trait FloatExt2<T>: bounds::Pair<T> {
 }
 
 impl<F: FloatExt> FloatExt2<F> for (F, F) {
-	fn lerp(self, t: F) -> F {
+	fn lerp(self, t: F) -> F { F::lerp(self.0, self.1, t) }
-		F::lerp(self.0, self.1, t)
+	fn inverse_lerp(self, y: F) -> F { F::inverse_lerp(self.0, self.1, y) }
-	}
-	fn inverse_lerp(self, y: F) -> F {
-		F::inverse_lerp(self.0, self.1, y)
-	}
 }
 
 #[cfg(test)]

src/bin/flat/fns.rs

@@ -52,73 +52,29 @@ 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 {
+	if t.abs() <= lim { f(t) } else { t + t.signum() * (val - 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 {
+	if t.abs() <= lim { f(t) } else { val }
-		f(t)
-	} else {
-		val
-	}
 }
 
 impl QuadraticAccelerator {
-	fn a(&self) -> f32 {
+	fn a(&self) -> f32 { -(self.external - self.internal) / self.internal.powi(2) }
-		-(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 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 {
+	pub fn x(&self, u: f32) -> f32 { extend_linear(u, |u| (self.a() * u.abs() + self.b()) * u, self.internal, self.external) }
-		extend_linear(
+	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) }
-			u,
+	pub fn dx(&self, u: f32) -> f32 { extend_const(u, |u| 2.0 * self.a() * u.abs() + self.b(), self.internal, 1.0) }
-			|u| (self.a() * u.abs() + self.b()) * u,
+	pub fn du(&self, x: f32) -> f32 { extend_const(x, |x| 1.0 / self.root(x), self.external, 1.0) }
-			self.internal,
+	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) }
-			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, assert_abs_diff_eq, AbsDiffEq};
+	use approx::{abs_diff_eq, AbsDiffEq, assert_abs_diff_eq};
 
 	use super::*;
 
@@ -137,45 +93,23 @@ 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!(
+			assert!(abs_diff_eq!(df_expl, df_num, epsilon = ε), "At x={x}, df/dx:\nnumerical: {df_num}\nexplicit:  {df_expl}\n");
-				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(
+		test_limiter(SmoothstepLimiter { min: 20.0, max: 30.0 }, 20.0, 30.0, 1.0 / 32.0);
-			SmoothstepLimiter {
-				min: 20.0,
-				max: 30.0,
-			},
-			20.0,
-			30.0,
-			1.0 / 32.0,
-		);
 	}
 
 	#[test]
 	fn test_smootherstep_limiter() {
-		test_limiter(
+		test_limiter(SmootherstepLimiter { min: 20.0, max: 30.0 }, 20.0, 30.0, 1.0 / 32.0);
-			SmootherstepLimiter {
-				min: 20.0,
-				max: 30.0,
-			},
-			20.0,
-			30.0,
-			1.0 / 32.0,
-		);
 	}
 
 	#[test]
 	fn test_quadratic_accelerator() {
-		let testee = super::QuadraticAccelerator {
+		let testee = super::QuadraticAccelerator { internal: 100.0, external: 150.0 };
-			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;
@@ -187,51 +121,30 @@ 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!(
+			assert!(abs_diff_eq!(x, xux, epsilon = ε), "At x={x}:\nu(x):    {ux}\nx(u(x)): {xux}\n");
-				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!(
+			assert!(abs_diff_eq!(du_expl, du_num, epsilon = ε), "At x={x}, du/dx:\nnumerical: {du_num}\nexplicit:  {du_expl}\n");
-				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!(
+			assert!(abs_diff_eq!(dudx, 1.0, epsilon = ε), "At x={x}:\ndu/dx * dx/du: {dudx}\n");
-				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!(
+			assert!(abs_diff_eq!(d2u_expl, d2u_num, epsilon = ε), "At x={x}, d^2u/dx^2:\nnumerical: {d2u_num}\nexplicit:  {d2u_expl}\n");
-				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!(
+			assert!(abs_diff_eq!(u, uxu, epsilon = ε), "At u={u}:\nx(u):    {xu}\nu(x(u)): {uxu}\n");
-				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!(
+			assert!(abs_diff_eq!(dx_expl, dx_num, epsilon = ε), "At u={u}, dx/du:\nnumerical: {dx_num}\nexplicit:  {dx_expl}\n");
-				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!(
+			assert!(abs_diff_eq!(dudx, 1.0, epsilon = ε), "At u={u}:\ndu/dx * dx/du: {dudx}\n");
-				abs_diff_eq!(dudx, 1.0, epsilon = ε),
-				"At u={u}:\ndu/dx * dx/du: {dudx}\n"
-			);
 		}
 	}
 }

src/bin/flat/main.rs

@@ -4,26 +4,23 @@ use flo_canvas::*;
 use flo_draw::*;
 use glam::*;
 
-use crate::types::FlatTraceResult;
+use riemann::{Metric, trace_iter};
-use riemann::{trace_iter, Metric};
 use tube::metric::Tube;
 use tube::Space;
 use tube::Subspace::{Boundary, Inner, Outer};
 use types::{Location, Object, Ray};
 
-mod float_fun;
-mod fns;
 mod riemann;
+mod fns;
+mod float_fun;
 mod tube;
 mod types;
 
 const DT: f32 = 0.1;
 
-fn draw_loop(gc: &mut Vec<Draw>, mut pts: impl Iterator<Item = Vec2>) {
+fn draw_loop(gc: &mut Vec<Draw>, mut pts: impl Iterator<Item=Vec2>) {
 	gc.new_path();
-	let Some(first) = pts.next() else {
+	let Some(first) = pts.next() else { return; };
-		return;
-	};
 	gc.move_to(first.x, first.y);
 	for pt in pts {
 		gc.line_to(pt.x, pt.y);
@@ -71,21 +68,10 @@ 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(
+			draw_fan_2(gc, &space, vec2(-2.5 * tube.outer_radius, 1.25 * tube.external_halflength), vec2(1.0, -1.0), 1.0);
-				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(
+			draw_track(gc, &space, vec2(-0.5 * tube.inner_radius, -1.25 * tube.external_halflength), vec2(0.1, 1.0));
-				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 {
@@ -96,113 +82,83 @@ pub fn main() {
 				gc.fill();
 
 				gc.stroke_color(Color::Rgba(0.0, 0.0, 0.0, 0.5));
-				draw_loop(
+				draw_loop(gc, itertools_num::linspace(0.0, 2.0 * PI, circle_segments).skip(1).map(|φ| {
-					gc,
+					let dir = Vec2::from_angle(φ) * obj.r;
-					itertools_num::linspace(0.0, 2.0 * PI, circle_segments)
+					let dir = obj.loc.rot * dir;
-						.skip(1)
+					pos + dir
-						.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(
+				draw_loop(gc, itertools_num::linspace(0.0, 2.0 * PI, circle_segments).skip(1).map(|φ| {
-					gc,
+					let dir = Vec2::from_angle(φ) * obj.r;
-					itertools_num::linspace(0.0, 2.0 * PI, circle_segments)
+					let dir = obj.loc.rot * dir;
-						.skip(1)
+					space.trace_step(Ray { pos, dir }).pos
-						.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(
+				draw_loop(gc, itertools_num::linspace(0.0, 2.0 * PI, circle_segments).skip(1).map(|φ| {
-					gc,
+					let n = obj.r.floor();
-					itertools_num::linspace(0.0, 2.0 * PI, circle_segments)
+					let d = obj.r / n;
-						.skip(1)
+					let dir = Vec2::from_angle(φ);
-						.map(|φ| {
+					let dir = obj.loc.rot * dir * d;
-							let n = obj.r.floor();
+					space.trace_iter(Ray { pos, dir }).nth(n as usize).unwrap().pos
-							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
-						}),
-				);
 			}
 		});
 	});
 }
 
-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()
-}
-
-fn draw_cross(gc: &mut Vec<Draw>, pos: Vec2, r: f32) {
-	gc.move_to(pos.x - r, pos.y - r);
-	gc.line_to(pos.x + r, pos.y + r);
-	gc.move_to(pos.x - r, pos.y + r);
-	gc.line_to(pos.x + r, pos.y - r);
-}
-
 fn draw_ray_2(gc: &mut Vec<Draw>, space: &Space, base: Vec2, dir: Vec2) {
-	fn trace_to_flat(gc: &mut Vec<Draw>, space: &Space, ray: Ray) -> (Ray, FlatTraceResult) {
-		for ray in space.trace_iter(ray).skip(1) {
-			gc.line_to(ray.pos.x, ray.pos.y);
-			match space.which_subspace(ray.pos) {
-				Inner => return (ray, space.trace_inner(ray)),
-				Outer => return (ray, space.trace_outer(ray)),
-				Boundary => continue,
-			};
-		}
-		unreachable!("Space::trace_iter terminated!")
-	}
-
 	let mut hits = Vec::<Draw>::new();
 	let dir = space.tube.globalize(base, dir);
 	gc.new_path();
 	gc.move_to(base.x, base.y);
-	let mut ray = Ray {
+	let mut ray = Ray { pos: base, dir: space.tube.normalize_vec_at(base, dir) * DT };
-		pos: base,
+	for _ in 0..10000 {
-		dir: space.tube.normalize_vec_at(base, dir) * DT,
+		ray = space.trace_step(ray);
-	};
+		gc.line_to(ray.pos.x, ray.pos.y);
-	for _ in 0..100 {
+		if ray.pos.abs().cmpgt(Vec2::splat(1000.0)).any() {
-		let ret;
+			break;
-		(ray, ret) = trace_to_flat(gc, space, ray);
+		}
+		let sub = space.which_subspace(ray.pos);
+		if sub == Boundary {
+			continue;
+		}
 		gc.stroke();
 		gc.new_dash_pattern();
 		// gc.dash_length(6.0);
 		gc.new_path();
 		gc.move_to(ray.pos.x, ray.pos.y);
+		let ret = match sub {
+			Inner => space.trace_inner(ray),
+			Outer => space.trace_outer(ray),
+			Boundary => panic!(),
+		};
 		for hit in ret.objects {
 			let obj = space.objs[hit.id as usize];
-			let apx_hit_pos = rel_to_abs(&space.tube, &obj.loc, hit.rel.pos, 128);
+			hits.move_to(obj.loc.pos.x, obj.loc.pos.y);
-			// assert_abs_diff_eq!(apx_hit_pos, hit.pos, epsilon=1.0);
+			for pt in trace_iter(&space.tube, obj.loc.pos, obj.loc.rot * hit.rel.pos, hit.rel.pos.length() / 100.0).take(100) {
+				hits.line_to(pt.x, pt.y);
+			}
+			hits.circle(hit.pos.x, hit.pos.y, 1.5);
 			let Ray { pos: rel, dir } = hit.rel;
-			let diff = rel.dot(dir).powi(2)
+			let diff = rel.dot(dir).powi(2) - dir.length_squared() * (rel.length_squared() - obj.r.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;
-			let pos2 = rel_to_abs(&space.tube, &obj.loc, rel2, 128);
+			let pos2 = trace_iter(&space.tube, obj.loc.pos, obj.loc.rot * rel2, rel2.length() / 100.0).nth(100).unwrap();
-			draw_cross(&mut hits, pos2, 1.0);
+			hits.move_to(pos2.x - 1.0, pos2.y - 1.0);
+			hits.line_to(pos2.x + 1.0, pos2.y + 1.0);
+			hits.move_to(pos2.x - 1.0, pos2.y + 1.0);
+			hits.line_to(pos2.x + 1.0, pos2.y - 1.0);
 		}
 		let a = ray.pos;
-		if let Some(r) = ret.end {
+		ray = match ret.end {
-			ray = r
+			Some(r) => r,
-		} else {
+			None => {
-			ray = ray.forward(1000.0 / DT);
+				ray = ray.forward(1000.0 / DT);
-			gc.line_to(ray.pos.x, ray.pos.y);
+				gc.line_to(ray.pos.x, ray.pos.y);
-			break;
+				break;
-		}
+			}
+		};
 		for p in space.line(a, ray.pos, 10.0) {
 			gc.line_to(p.x, p.y);
 		}
@@ -245,10 +201,7 @@ fn draw_track(gc: &mut Vec<Draw>, 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 {
+	let mut loc = Location { pos: start, rot: mat2(dir, vec2(-dir.y, dir.x)) };
-		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;
@@ -290,18 +243,8 @@ trait Renderable {
 impl Renderable for Tube {
 	fn render(&self, gc: &mut Vec<Draw>) {
 		gc.new_path();
-		gc.rect(
+		gc.rect(-self.outer_radius, -self.external_halflength, self.outer_radius, self.external_halflength);
-			-self.outer_radius,
+		gc.rect(-self.inner_radius, -self.external_halflength, self.inner_radius, self.external_halflength);
-			-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();

src/bin/flat/riemann.rs

@@ -62,6 +62,7 @@ pub struct TraceIter<'a, M: Metric> {
 	space: &'a M,
 	p: Vec2,
 	v: Vec2,
+	dt: f32,
 }
 
 impl<'a, M: Metric> Iterator for TraceIter<'a, M> {
@@ -69,8 +70,8 @@ impl<'a, M: Metric> Iterator for TraceIter<'a, M> {
 
 	fn next(&mut self) -> Option<Self::Item> {
 		let a: Vec2 = -contract2(krist(self.space, self.p), self.v);
-		self.v += a;
+		self.v = self.v + a * self.dt;
-		self.p += self.v;
+		self.p = self.p + self.v * self.dt;
 		Some(self.p)
 	}
 }
@@ -79,7 +80,8 @@ pub fn trace_iter<M: Metric>(space: &M, base: Vec2, dir: Vec2, dt: f32) -> Trace
 	TraceIter {
 		space,
 		p: base,
-		v: dt * space.normalize_vec_at(base, dir),
+		v: space.normalize_vec_at(base, dir),
+		dt,
 	}
 }
 
@@ -88,11 +90,7 @@ 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| {
+	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::<f32>())
-		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::<f32>()
-	})
 }
 
 fn dir_deriv(f: impl Fn(Vec2) -> Mat2, pos: Vec2, delta: Vec2) -> Mat2 {
@@ -171,109 +169,27 @@ pub mod samples {
 
 #[cfg(test)]
 mod tests {
-	use super::*;
+	use glam::{Mat2, mat2, vec2};
-	use approx::assert_abs_diff_eq;
-
-	use glam::{mat2, vec2, Mat2};
 	use rand::{Rng, SeedableRng};
 
+	use super::{Decomp2, Metric, samples};
+
 	#[test]
 	fn uniform_scaled_metric() {
 		let mut rng = rand_pcg::Pcg64Mcg::seed_from_u64(17);
-		let metric = samples::ScaledMetric {
+		let metric = samples::ScaledMetric { scale: vec2(3., 4.) };
-			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!(
+		assert_eq!(metric.inverse_at(rng.gen()), Mat2::from_cols_array(&[1. / 9., 0., 0., 1. / 16.]));
-			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!(
+		assert_eq!(metric.normalize_vec_at(rng.gen(), vec2(1., 0.)), vec2(1. / 3., 0.));
-			metric.normalize_vec_at(rng.gen(), vec2(1., 0.)),
+		assert_eq!(metric.normalize_vec_at(rng.gen(), vec2(0., 1.)), vec2(0., 1. / 4.));
-			vec2(1. / 3., 0.)
+		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(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!(
+		assert_eq!(metric.globalize(rng.gen(), vec2(1., 1.)), vec2(1. / 3., 1. / 4.));
-			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
-		);
 	}
 }

src/bin/flat/tube/coords.rs

@@ -1,4 +1,4 @@
-use glam::{vec2, Mat2, Vec2};
+use glam::{Mat2, Vec2, vec2};
 
 use crate::riemann::Metric;
 use crate::types::{Location, Ray};
@@ -10,22 +10,15 @@ pub trait FlatCoordinateSystem<T> {
 	fn global_to_flat(&self, v: T) -> T;
 }
 
-pub trait FlatRegion:
+pub trait FlatRegion: FlatCoordinateSystem<Vec2> + FlatCoordinateSystem<Ray> + FlatCoordinateSystem<Location> {
-	FlatCoordinateSystem<Vec2> + FlatCoordinateSystem<Ray> + FlatCoordinateSystem<Location>
-{
 	// Измеряет расстояние до выхода за пределы области вдоль луча ray. Луч задаётся в плоской СК.
-	fn distance_to_boundary(&self, _ray: Ray) -> Option<f32> {
+	fn distance_to_boundary(&self, _ray: Ray) -> Option<f32> { None }
-		None
-	}
 }
 
 trait MetricCS: FlatCoordinateSystem<Vec2> {
 	fn global_metric(&self) -> &impl Metric;
 	fn flat_to_global_tfm_at(&self, pos: Vec2) -> Mat2 {
-		self.global_metric()
+		self.global_metric().sqrt_at(self.flat_to_global(pos)).inverse().into()
-			.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()
@@ -66,11 +59,7 @@ impl<T: FlatCoordinateSystem<Vec2> + MetricCS> FlatCoordinateSystem<Location> fo
 
 pub struct InnerCS(pub Tube);
 
-impl MetricCS for InnerCS {
+impl MetricCS for InnerCS { fn global_metric(&self) -> &impl Metric { &self.0 } }
-	fn global_metric(&self) -> &impl Metric {
-		&self.0
-	}
-}
 
 impl FlatCoordinateSystem<Vec2> for InnerCS {
 	fn flat_to_global(&self, pos: Vec2) -> Vec2 {
@@ -85,31 +74,20 @@ impl FlatCoordinateSystem<Vec2> for InnerCS {
 
 impl FlatRegion for InnerCS {
 	fn distance_to_boundary(&self, ray: Ray) -> Option<f32> {
-		Rect {
+		Rect { size: vec2(self.0.inner_radius, self.0.internal_halflength) }.trace_out_of(ray)
-			size: vec2(self.0.inner_radius, self.0.internal_halflength),
-		}
-		.trace_out_of(ray)
 	}
 }
 
 pub struct OuterCS(pub Tube);
 
-impl MetricCS for OuterCS {
+impl MetricCS for OuterCS { fn global_metric(&self) -> &impl Metric { &self.0 } }
-	fn global_metric(&self) -> &impl Metric {
-		&self.0
-	}
-}
 
 impl FlatCoordinateSystem<Vec2> for OuterCS {
 	fn flat_to_global(&self, pos: Vec2) -> Vec2 {
-		let inner = Rect {
+		let inner = Rect { size: vec2(self.0.inner_radius + 1.0, self.0.external_halflength) };
-			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
+			let y = self.0.y(v - v.signum() * (self.0.external_halflength - self.0.internal_halflength));
-				.0
-				.y(v - v.signum() * (self.0.external_halflength - self.0.internal_halflength));
 			vec2(x, y)
 		} else {
 			pos
@@ -117,13 +95,10 @@ impl FlatCoordinateSystem<Vec2> for OuterCS {
 	}
 
 	fn global_to_flat(&self, pos: Vec2) -> Vec2 {
-		let inner = Rect {
+		let inner = Rect { size: vec2(self.0.inner_radius + 1.0, self.0.external_halflength) };
-			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)
+			let v = self.0.v(y) + y.signum() * (self.0.external_halflength - self.0.internal_halflength);
-				+ y.signum() * (self.0.external_halflength - self.0.internal_halflength);
 			vec2(u, v) // в плоском продолжении СК Outer на область Inner
 		} else {
 			pos
@@ -133,17 +108,14 @@ impl FlatCoordinateSystem<Vec2> for OuterCS {
 
 impl FlatRegion for OuterCS {
 	fn distance_to_boundary(&self, ray: Ray) -> Option<f32> {
-		Rect {
+		Rect { size: vec2(self.0.outer_radius, self.0.external_halflength) }.trace_into(ray)
-			size: vec2(self.0.outer_radius, self.0.external_halflength),
-		}
-		.trace_into(ray)
 	}
 }
 
 #[cfg(test)]
 mod tests {
-	use approx::{assert_abs_diff_eq, AbsDiffEq};
+	use approx::{AbsDiffEq, assert_abs_diff_eq};
-	use glam::{mat2, vec2, Mat2, Vec2};
+	use glam::{Mat2, mat2, vec2, Vec2};
 	use itertools_num::linspace;
 
 	use crate::riemann::samples;
@@ -154,179 +126,57 @@ mod tests {
 	fn uniform_scaled_metric() {
 		struct Scaled(samples::ScaledMetric, Vec2);
 		impl FlatCoordinateSystem<Vec2> for Scaled {
-			fn flat_to_global(&self, pos: Vec2) -> Vec2 {
+			fn flat_to_global(&self, pos: Vec2) -> Vec2 { (pos - self.1) / self.0.scale }
-				(pos - self.1) / self.0.scale
+			fn global_to_flat(&self, pos: Vec2) -> Vec2 { pos * self.0.scale + self.1 }
-			}
-			fn global_to_flat(&self, pos: Vec2) -> Vec2 {
-				pos * self.0.scale + self.1
-			}
 		}
 		impl MetricCS for Scaled {
-			fn global_metric(&self) -> &impl Metric {
+			fn global_metric(&self) -> &impl Metric { &self.0 }
-				&self.0
-			}
 		}
-		let cs = Scaled(
+		let cs = Scaled(samples::ScaledMetric { scale: vec2(3., 4.) }, vec2(2., 3.));
-			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!(
+		assert_eq!(cs.global_to_flat(Ray { pos: vec2(7., 3.), dir: vec2(3., 2.) }), Ray { pos: vec2(23., 15.), dir: vec2(9., 8.) });
-			cs.global_to_flat(Ray {
+		assert_eq!(cs.flat_to_global(Ray { pos: vec2(23., 15.), dir: vec2(9., 8.) }), Ray { pos: vec2(7., 3.), dir: vec2(3., 2.) });
-				pos: vec2(7., 3.),
+		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.)) });
-				dir: vec2(3., 2.)
+		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.)) });
-			}),
-			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(
+	fn test_flat_region(region: &impl FlatRegion, range_global: (Vec2, Vec2), range_flat: (Vec2, Vec2)) {
-		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 {
 			($at: expr, $left: expr, $right: expr) => {
 				let at = $at;
 				let left = $left;
 				let right = $right;
-				assert!(
+				assert!(left.abs_diff_eq(right, ε), "Assertion failed at {at}:\n  left: {left} = {}\n right: {right} = {}", stringify!($left), stringify!($right));
-					left.abs_diff_eq(right, ε),
-					"Assertion failed at {at}:\n  left: {left} = {}\n right: {right} = {}",
-					stringify!($left),
-					stringify!($right)
-				);
 			};
 		}
-		fn check_range(
+		fn check_range(name_a: &str, a: Vec2, range_a: (Vec2, Vec2), name_b: &str, b: Vec2, range_b: (Vec2, Vec2)) {
-			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, ε) {
+			if a.x.abs_diff_eq(&range_a.0.x, ε) { assert_abs_diff_eq!(b.x, range_b.0.x, epsilon=ε); }
-				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.0.y, ε) {
+			if a.y.abs_diff_eq(&range_a.1.y, ε) { assert_abs_diff_eq!(b.y, range_b.1.y, epsilon=ε); }
-				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(
+				check_range("global", pos_global, range_global, "flat", pos_flat, range_flat);
-					"global",
+				assert_eq_at!(pos_global, region.global_to_flat(Location { pos: pos_global, rot: Mat2::IDENTITY }).pos, pos_flat);
-					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!(
+				assert_eq_at!(pos_global, region.flat_to_global(region.global_to_flat(Location { pos: pos_global, rot: Mat2::IDENTITY })).rot, Mat2::IDENTITY);
-					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(
+				check_range("flat", pos_flat, range_flat, "global", pos_global, range_global);
-					"flat",
+				assert_eq_at!(pos_flat, region.flat_to_global(Location { pos: pos_flat, rot: Mat2::IDENTITY }).pos, pos_global);
-					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!(
+				assert_eq_at!(pos_flat, region.global_to_flat(region.flat_to_global(Location { pos: pos_global, rot: Mat2::IDENTITY })).rot, Mat2::IDENTITY);
-					pos_flat,
-					region
-						.global_to_flat(region.flat_to_global(Location {
-							pos: pos_global,
-							rot: Mat2::IDENTITY
-						}))
-						.rot,
-					Mat2::IDENTITY
-				);
 			}
 		}
 	}
@@ -339,21 +189,9 @@ mod tests {
 			internal_halflength: 100.0,
 			external_halflength: 300.0,
 		});
-		test_flat_region(
+		test_flat_region(&mapper, (vec2(-30.0, -300.0), vec2(30.0, 300.0)), (vec2(-30.0, -100.0), vec2(30.0, 100.0)));
-			&mapper,
+		test_flat_region(&mapper, (vec2(-60.0, -400.0), vec2(60.0, -300.0)), (vec2(-60.0, -200.0), vec2(60.0, -100.0)));
-			(vec2(-30.0, -300.0), vec2(30.0, 300.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)));
-			(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]
@@ -365,68 +203,23 @@ mod tests {
 			external_halflength: 300.0,
 		});
 		// TODO replace 200.20016 with something sane
-		test_flat_region(
+		test_flat_region(&mapper, (vec2(-30.0, -300.0), vec2(30.0, -1.0)), (vec2(-30.0, -300.0), vec2(30.0, -200.20016)));
-			&mapper,
+		test_flat_region(&mapper, (vec2(-30.0, 1.0), vec2(30.0, 300.0)), (vec2(-30.0, 200.20016), vec2(30.0, 300.0)));
-			(vec2(-30.0, -300.0), vec2(30.0, -1.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)));
-			(vec2(-30.0, -300.0), vec2(30.0, -200.20016)),
+		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, 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!(
+				assert_eq!(mapper.global_to_flat(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos.x, x);
-					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
+				let pp = mapper.global_to_flat(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos;
-					.global_to_flat(Location {
+				let np = mapper.global_to_flat(Location { pos: vec2(-x, y), rot: Mat2::IDENTITY }).pos;
-						pos: vec2(x, y),
+				let pn = mapper.global_to_flat(Location { pos: vec2(x, -y), rot: Mat2::IDENTITY }).pos;
-						rot: Mat2::IDENTITY,
+				let nn = mapper.global_to_flat(Location { pos: vec2(-x, -y), rot: Mat2::IDENTITY }).pos;
-					})
-					.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));
@@ -435,42 +228,18 @@ mod tests {
 		// clean boundary
 		for x in linspace(50., 60., 20) {
 			for y in linspace(0., 320., 20) {
-				assert_eq!(
+				assert_eq!(mapper.global_to_flat(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos.y, y);
-					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!(
+				assert_eq!(mapper.global_to_flat(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos.y, y);
-					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
+				let v = mapper.global_to_flat(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos.y;
-					.global_to_flat(Location {
-						pos: vec2(x, y),
-						rot: Mat2::IDENTITY,
-					})
-					.pos
-					.y;
 				assert!(v > 200.0);
 				assert!(v > y);
 			}

src/bin/flat/tube/metric.rs

@@ -1,4 +1,4 @@
-use glam::{f32, vec2, Mat2, Vec2};
+use glam::{f32, Mat2, Vec2, vec2};
 
 use crate::fns::{self, Limiter};
 use crate::riemann::{Decomp2, Metric, Tens2};
@@ -12,31 +12,13 @@ pub struct Tube {
 }
 
 impl Tube {
-	fn fx(&self) -> impl Limiter {
+	fn fx(&self) -> impl Limiter { fns::SmootherstepLimiter { min: self.inner_radius, max: self.outer_radius } }
-		fns::SmootherstepLimiter {
+	fn fy(&self) -> fns::QuadraticAccelerator { fns::QuadraticAccelerator { internal: self.internal_halflength, external: self.external_halflength } }
-			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 {
+	pub fn y(&self, v: f32) -> f32 { self.fy().x(v) }
-		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 v(&self, y: f32) -> f32 {
+	pub fn dv(&self, y: f32) -> f32 { self.fy().du(y) }
-		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 {
@@ -71,7 +53,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};
@@ -84,9 +66,7 @@ mod test {
 	fn test_tube_metric_derivs() {
 		struct Approx(Tube);
 		impl Metric for Approx {
-			fn sqrt_at(&self, pos: Vec2) -> Decomp2 {
+			fn sqrt_at(&self, pos: Vec2) -> Decomp2 { self.0.sqrt_at(pos) }
-				self.0.sqrt_at(pos)
-			}
 		}
 		let testee = Tube {
 			inner_radius: 30.0,
@@ -98,13 +78,8 @@ 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 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 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);
@@ -135,10 +110,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.x, b.x, epsilon=1.0e-2);
-				assert_abs_diff_eq!(traced.pos.y, b.y, epsilon = 1.0e1);
+				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.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.dir.y, dir.y, epsilon=1.0e-2);
 			}
 		}
 	}
@@ -158,25 +133,17 @@ mod test {
 		let ε = 1e-3;
 		let off = 10.0;
 		let steps = 4096;
-		for ax in linspace(
+		for ax in linspace(-space.tube.inner_radius + ε, space.tube.inner_radius - ε, 20) {
-			-space.tube.inner_radius + ε,
+			for bx in linspace(-space.tube.inner_radius + ε, space.tube.inner_radius - ε, 20) {
-			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.x, b.x, epsilon=1.0e-2);
-				assert_abs_diff_eq!(traced.pos.y, b.y, epsilon = 1.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-3);
+				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.dir.y, dir.y, epsilon=1.0e-3);
 			}
 		}
 	}
@@ -201,10 +168,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.x, b.x, epsilon=1.0e-1);
-			assert_abs_diff_eq!(traced.pos.y, b.y, epsilon = 1.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.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.dir.y, dir.y, epsilon=1.0e-2);
 		}
 	}
 
@@ -228,10 +195,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.x, b.x, epsilon=2.0e0);
-			assert_abs_diff_eq!(traced.pos.y, b.y, epsilon = 1.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.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.dir.y, dir.y, epsilon=1.0e-2);
 		}
 	}
 }

src/bin/flat/tube/mod.rs

@@ -1,4 +1,4 @@
-use glam::{bool, f32, vec2, Mat2, Vec2};
+use glam::{bool, f32, Mat2, Vec2, 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};
 
-mod coords;
 pub mod metric;
+mod coords;
 
 pub struct Space {
 	pub tube: Tube,
@@ -48,16 +48,12 @@ impl Space {
 	/// Выполняет один шаг перемещения. Работает в любой части пространства.
 	/// off задаётся в локальной СК. Рекомендуется считать небольшими шагами.
 	pub fn move_step(&self, loc: Location, off: Vec2) -> Location {
-		let corr =
+		let corr = Mat2::IDENTITY - riemann::contract(riemann::krist(&self.tube, loc.pos), loc.rot * off);
-			Mat2::IDENTITY - riemann::contract(riemann::krist(&self.tube, loc.pos), loc.rot * off);
 		let p = loc.pos + corr * loc.rot * off;
-		Location {
+		Location { pos: p, rot: corr * loc.rot }
-			pos: p,
-			rot: corr * loc.rot,
-		}
 	}
 
-	pub fn trace_iter(&self, ray: Ray) -> impl Iterator<Item = Ray> + '_ {
+	pub fn trace_iter(&self, ray: Ray) -> impl Iterator<Item=Ray> + '_ {
 		std::iter::successors(Some(ray), |&ray| Some(self.trace_step(ray)))
 	}
 
@@ -89,28 +85,15 @@ impl Space {
 	}
 
 	fn list_objects(&self, tfm: impl Fn(Location) -> Location) -> Vec<Object> {
-		self.objs
+		self.objs.iter().map(|&Object { id, loc, r }| Object { id, loc: tfm(loc), r }).collect()
-			.iter()
-			.map(|&Object { id, loc, r }| Object {
-				id,
-				loc: tfm(loc),
-				r,
-			})
-			.collect()
 	}
 
-	fn hit_objects(
+	fn hit_objects(objs: &[Object], ray: Ray, limit: Option<f32>, globalize: impl Fn(Vec2) -> Vec2) -> Vec<Hit> {
-		objs: &[Object],
-		ray: Ray,
-		limit: Option<f32>,
-		globalize: impl Fn(Vec2) -> Vec2,
-	) -> Vec<Hit> {
 		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)
+				let diff = rel.dot(ray.dir).powi(2) - ray.dir.length_squared() * (rel.length_squared() - obj.r.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))
@@ -121,17 +104,8 @@ impl Space {
 			.filter(|&(_, t)| t >= 0.0 && t < limit)
 			.map(|(obj, t)| {
 				let pos = ray.forward(t).pos;
-				let rel = obj.loc.rot.inverse()
+				let rel = obj.loc.rot.inverse() * Ray { pos: pos - obj.loc.pos, dir: ray.dir };
-					* Ray {
+				Hit { id: obj.id, distance: t, pos: globalize(pos), rel }
-						pos: pos - obj.loc.pos,
-						dir: ray.dir,
-					};
-				Hit {
-					id: obj.id,
-					distance: t,
-					pos: globalize(pos),
-					rel,
-				}
 			})
 			.collect()
 	}
@@ -144,9 +118,7 @@ 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)
+				(1..=n).map(|k| cs.flat_to_global(a.lerp(b, k as f32 / n as f32))).collect()
-					.map(|k| cs.flat_to_global(a.lerp(b, k as f32 / n as f32)))
-					.collect()
 			}
 			Boundary => panic!("Can't draw a line here!"),
 		}
@@ -160,10 +132,7 @@ struct Rect {
 impl Rect {
 	/// Отражает луч, чтобы все координаты направления были положительны (допустимо благодаря симметрии Rect).
 	fn flip_ray(ray: Ray) -> Ray {
-		Ray {
+		Ray { pos: ray.pos * ray.dir.signum(), dir: ray.dir.abs() }
-			pos: ray.pos * ray.dir.signum(),
-			dir: ray.dir.abs(),
-		}
 	}
 
 	fn is_inside(&self, pt: Vec2) -> bool {
@@ -176,12 +145,8 @@ 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 {
+		if t < 0.0 { return None; }
-			return None;
+		if pt.cmpgt(self.size).any() { return None; }
-		}
-		if pt.cmpgt(self.size).any() {
-			return None;
-		}
 		Some(t)
 	}
 
@@ -196,149 +161,27 @@ impl Rect {
 
 #[test]
 fn test_rect() {
-	assert_eq!(
+	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) });
-		Rect::flip_ray(Ray {
+	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) });
-			pos: vec2(2.0, 3.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) });
-			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) });
-		}),
-		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 {
+	let r = Rect { size: vec2(2.0, 3.0) };
-		size: vec2(2.0, 3.0),
-	};
 
-	assert_eq!(
+	assert_eq!(r.trace_into(Ray { pos: vec2(3.0, 3.0), dir: vec2(1.0, 1.0) }), None);
-		r.trace_into(Ray {
+	assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 2.0), dir: vec2(1.0, 0.0) }), Some(1.0));
-			pos: vec2(3.0, 3.0),
+	assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 2.0), dir: vec2(-1.0, 0.0) }), None);
-			dir: vec2(1.0, 1.0)
+	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);
-		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!(
+	assert_eq!(r.trace_into(Ray { pos: vec2(2.0, 3.0), dir: vec2(1.0, 1.0) }), None);
-		r.trace_into(Ray {
+	assert_eq!(r.trace_into(Ray { pos: vec2(-2.0, 3.0), dir: vec2(-1.0, 1.0) }), None);
-			pos: vec2(2.0, 3.0),
+	assert_eq!(r.trace_into(Ray { pos: vec2(2.0, 3.0), dir: vec2(-1.0, -1.0) }), Some(0.0));
-			dir: vec2(1.0, 1.0)
+	assert_eq!(r.trace_into(Ray { pos: vec2(2.0, -3.0), dir: vec2(-1.0, 1.0) }), Some(0.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!(
+	assert_eq!(r.trace_out_of(Ray { pos: vec2(0.0, 0.0), dir: vec2(1.0, 1.0) }), Some(2.0));
-		r.trace_out_of(Ray {
+	assert_eq!(r.trace_out_of(Ray { pos: vec2(0.0, 0.0), dir: vec2(0.0, 1.0) }), Some(3.0));
-			pos: vec2(0.0, 0.0),
+	assert_eq!(r.trace_out_of(Ray { pos: vec2(0.0, 1.0), dir: vec2(0.0, -1.0) }), Some(4.0));
-			dir: vec2(1.0, 1.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));
-		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)
-	);
 }

src/bin/flat/types.rs

@@ -8,10 +8,7 @@ pub struct Ray {
 
 impl Ray {
 	pub fn forward(&self, dist: f32) -> Ray {
-		Ray {
+		Ray { pos: self.pos + self.dir * dist, dir: self.dir }
-			pos: self.pos + self.dir * dist,
-			dir: self.dir,
-		}
 	}
 }
 
@@ -19,10 +16,7 @@ impl std::ops::Mul<Ray> for Mat2 {
 	type Output = Ray;
 
 	fn mul(self, rhs: Ray) -> Self::Output {
-		Ray {
+		Ray { pos: self * rhs.pos, dir: self * rhs.dir }
-			pos: self * rhs.pos,
-			dir: self * rhs.dir,
-		}
 	}
 }
 
@@ -46,7 +40,7 @@ pub struct Hit {
 	pub distance: f32,
 	pub id: i32,
 	pub pos: Vec2, // положение в основной СК
-	pub rel: Ray,  // в локальной ортонормированной СК объекта
+	pub rel: Ray, // в локальной ортонормированной СК объекта
 }
 
 #[derive(Clone, Debug)]

src/bin/mesh/main.rs

@@ -1,12 +1,13 @@
-use glam::*;
+mod mesh_loader;
-use refraction::mesh_loader::load_mesh;
+
-use refraction::mesh_tracer::{trace_to_mesh, Mesh};
+use std::fs::File;
-use show_image::{ImageInfo, ImageView, WindowOptions};
+use std::{env};
-use std::env;
 use std::error::Error;
 use std::f32::consts::PI;
-use std::fs::File;
+use std::io::{BufReader};
-use std::io::BufReader;
+use glam::*;
+use show_image::{ImageInfo, ImageView, WindowOptions};
+use crate::mesh_loader::{Face, load_mesh};
 
 const W: i32 = 320;
 const H: i32 = 240;
@@ -37,29 +38,56 @@ impl Image {
 }
 
 fn ypr_to_mat(ypr: Vec3) -> Mat3 {
-	let Vec3 {
+	let Vec3 { x: yaw, y: pitch, z: roll } = ypr;
-		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
 }
 
+type Mesh = [Face];
+
+struct TraceResult {
+	distance: f32,
+	normal: Vec3,
+}
+
+fn trace_to_mesh(mesh: &Mesh, base: Vec3, ray: Vec3) -> Option<TraceResult> {
+	let mut ret: Option<TraceResult> = None;
+	let mut dist = f32::INFINITY;
+	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 = m.inverse();
+			let rel = m * (base - f.vertices[0]);
+			if rel.z > dist {
+				continue;
+			}
+			dist = rel.z;
+			ret = Some(TraceResult {
+				distance: rel.z,
+				normal: f.normal,
+			});
+		}
+	}
+	ret
+}
+
+struct Location {
+	pos: Vec3,
+	rot: Vec4,
+}
+
 fn render(mesh: &Mesh, camera: impl Fn(Vec2) -> (Vec3, Vec3)) -> Image {
 	let bkg = vec3(0.0, 0.0, 0.0);
 	let mut img = Image {
@@ -79,9 +107,7 @@ fn render(mesh: &Mesh, camera: impl Fn(Vec2) -> (Vec3, Vec3)) -> Image {
 			} else {
 				bkg
 			};
-			let color = (color * 255.0)
+			let color = (color * 255.0).as_ivec3().clamp(IVec3::splat(0), IVec3::splat(255));
-				.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));
 		}
 	}
@@ -99,11 +125,7 @@ fn main() -> Result<(), Box<dyn Error>> {
 	let window = show_image::create_window("Raytracing", WindowOptions::default())?;
 	loop {
 		for phi in 0..360 {
-			let m_view = ypr_to_mat(vec3(
+			let m_view = ypr_to_mat(vec3((135.0 + phi as f32) * PI / 180.0, -30.0 * PI / 180.0, 0.0f32));
-				(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
@@ -122,3 +144,8 @@ fn main() -> Result<(), Box<dyn Error>> {
 		}
 	}
 }
+
+fn edge_dist(a: Vec3, b: Vec3, base: Vec3, dir: Vec3) -> f32 {
+	// Note: given that the input is not arbitrary but comes from a cartesian product of certain (a, b) pairs and certain (base, dir) pairs, this can be optimized from Cnm to an+bm+cnm with c<C.
+	mat3(b - a, base - a, -dir).determinant()
+}

src/bin/mesh/mesh_loader.rs

@@ -1,5 +1,5 @@
-use glam::{vec2, vec3, Vec2, Vec3};
 use std::io;
+use glam::{vec2, vec3, Vec2, Vec3};
 
 #[derive(Copy, Clone, Debug)]
 struct ObjVertex {
@@ -35,23 +35,14 @@ impl ObjMesh {
 	}
 
 	fn parse_fv(desc: &&str) -> ObjVertex {
-		let tokens: Vec<_> = desc
+		let tokens: Vec<_> = desc.split('/').map(|s| s.parse::<usize>().unwrap() - 1).collect();
-			.split('/')
-			.map(|s| s.parse::<usize>().unwrap() - 1)
-			.collect();
 		assert_eq!(tokens.len(), 3);
-		ObjVertex {
+		ObjVertex { vertex: tokens[0], tex_coord: tokens[1], normal: tokens[2] }
-			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 {
+		ObjFace { vertices: vertices.as_slice().try_into().unwrap() }
-			vertices: vertices.as_slice().try_into().unwrap(),
-		}
 	}
 
 	fn read(f: &mut impl io::BufRead) -> io::Result<ObjMesh> {
@@ -66,7 +57,13 @@ 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..])),
@@ -79,13 +76,12 @@ impl ObjMesh {
 	}
 
 	fn flatten(&self) -> Vec<Face> {
-		self.faces
+		self.faces.iter().map(|face| {
-			.iter()
+			Face {
-			.map(|face| Face {
 				vertices: face.vertices.map(|iv| self.vertices[iv.vertex]),
 				normal: self.normals[face.vertices[0].normal],
-			})
+			}
-			.collect()
+		}).collect()
 	}
 }
 

src/lib.rs

@@ -1,2 +0,0 @@
-pub mod mesh_loader;
-pub mod mesh_tracer;

src/mesh_tracer.rs

@@ -1,44 +0,0 @@
-use crate::mesh_loader::Face;
-use glam::{mat3, Vec3};
-
-pub type Mesh = [Face];
-
-pub struct TraceResult {
-	pub distance: f32,
-	pub normal: Vec3,
-}
-
-pub fn trace_to_mesh_all(
-	mesh: &Mesh,
-	base: Vec3,
-	ray: Vec3,
-) -> impl Iterator<Item = TraceResult> + '_ {
-	mesh.iter().filter_map(move |f| {
-		let fs = (0..3).map(|k| edge_dist(f.vertices[k], f.vertices[(k + 1) % 3], base, ray));
-		if fs.into_iter().any(|f| f < 0.0) {
-			return None;
-		}
-		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]);
-		Some(TraceResult {
-			distance: rel.z,
-			normal: f.normal,
-		})
-	})
-}
-
-pub fn trace_to_mesh(mesh: &Mesh, base: Vec3, ray: Vec3) -> Option<TraceResult> {
-	trace_to_mesh_all(mesh, base, ray)
-		.filter(|tr| tr.distance >= 0.0)
-		.min_by(|a, b| a.distance.total_cmp(&b.distance))
-}
-
-fn edge_dist(a: Vec3, b: Vec3, base: Vec3, dir: Vec3) -> f32 {
-	// Note: given that the input is not arbitrary but comes from a cartesian product of certain (a, b) pairs and certain (base, dir) pairs, this can be optimized from Cnm to an+bm+cnm with c<C.
-	mat3(b - a, base - a, -dir).determinant()
-}