Nice wide lines

Simplicity first!

.kateproject.build

@@ -0,0 +1,18 @@
+{
+    "Auto_generated": "This file is auto-generated. Any extra tags or formatting will be lost",
+    "target_sets": [
+        {
+            "cmake_config": "",
+            "directory": "%B",
+            "loaded_via_cmake": false,
+            "name": "Cargo",
+            "targets": [
+                {
+                    "build_cmd": "cargo build --bin wireframe",
+                    "name": "wireframe",
+                    "run_cmd": "cargo run --bin wireframe"
+                }
+            ]
+        }
+    ]
+}

Cargo.toml

@@ -21,6 +21,11 @@ show-image = "0.14.0"
 flo_draw = "0.3.1"
 flo_canvas = "0.3.1"
 itertools-num = "0.1.3"
+winit = "0.29"
+itertools = "0.13.0"
+wgpu = "22.1.0"
+bytemuck = { version = "1.18.0", features = ["derive"] }
+pollster = "0.3.0"
 
 [dev-dependencies]
 approx = "0.5.1"

src/bin/flat/main.rs

@@ -4,24 +4,14 @@ use flo_canvas::*;
 use flo_draw::*;
 use glam::*;
 
-use crate::ifaces::{DebugTraceable, Traceable};
+use refraction::ifaces::{DebugTraceable, Traceable};
-use crate::types::FlatTraceResult;
+use refraction::tube::metric::Tube;
-use refraction::mathx::MatExt;
+use refraction::tube::Space;
-use riemann::{trace_iter, Metric};
+use refraction::types::{Location, Object, Ray};
-use tube::metric::Tube;
+use refraction::utils::put_object;
-use tube::Space;
+use refraction::DT;
-use tube::Subspace::{Boundary, Inner, Outer};
-use types::{Location, Object, Ray};
 
-mod fns;
+fn draw_loop(gc: &mut Vec<Draw>, mut pts: impl Iterator<Item = Vec3>) {
-mod ifaces;
-mod riemann;
-mod tube;
-mod types;
-
-const DT: f32 = 0.1;
-
-fn draw_loop(gc: &mut Vec<Draw>, mut pts: impl Iterator<Item = Vec2>) {
 	gc.new_path();
 	let Some(first) = pts.next() else {
 		return;
@@ -52,23 +42,31 @@ pub fn main() {
 					id: k as i32,
 					loc: put_object(
 						&tube,
-						vec2(0.0, y * tube.external_halflength),
+						vec3(0.0, y * tube.external_halflength, 0.0),
-						Mat2::from_angle(y),
+						Mat3::from_mat2(Mat2::from_angle(y)),
 					),
 					r: 20.0,
 				})
 				.collect();
 			let space = Space { tube, objs };
-			let cam1 = put_object(&space.tube, vec2(-500., 0.), Mat2::IDENTITY);
+			let cam1 = put_object(&space.tube, vec3(-500., 0., 0.), Mat3::IDENTITY);
 			let cam2 = put_object(
 				&space.tube,
-				vec2(-2.5 * tube.outer_radius, 1.25 * tube.external_halflength),
+				vec3(
-				mat2(vec2(1., -1.), vec2(1., 1.)),
+					-2.5 * tube.outer_radius,
+					1.25 * tube.external_halflength,
+					0.,
+				),
+				mat3(vec3(1., -1., 0.), vec3(1., 1., 0.), vec3(0., 0., 1.)),
 			);
 			let cam3 = put_object(
 				&space.tube,
-				vec2(0.25 * tube.inner_radius, 0.25 * tube.external_halflength),
+				vec3(
-				mat2(vec2(0., -1.), vec2(1., 0.)),
+					0.25 * tube.inner_radius,
+					0.25 * tube.external_halflength,
+					0.,
+				),
+				mat3(vec3(0., -1., 0.), vec3(1., 0., 0.), vec3(0., 0., 1.)),
 			);
 
 			gc.canvas_height(500.0);
@@ -108,6 +106,7 @@ pub fn main() {
 						.skip(1)
 						.map(|φ| {
 							let dir = Vec2::from_angle(φ) * obj.r;
+							let dir = vec3(dir.x, dir.y, 0.);
 							let dir = obj.loc.rot * dir;
 							pos + dir
 						}),
@@ -119,6 +118,7 @@ pub fn main() {
 						.skip(1)
 						.map(|φ| {
 							let dir = Vec2::from_angle(φ) * obj.r;
+							let dir = vec3(dir.x, dir.y, 0.);
 							let dir = obj.loc.rot * dir;
 							space.trace_step(Ray { pos, dir }).pos
 						}),
@@ -132,6 +132,7 @@ pub fn main() {
 							let n = obj.r.floor();
 							let d = obj.r / n;
 							let dir = Vec2::from_angle(φ);
+							let dir = vec3(dir.x, dir.y, 0.);
 							let dir = obj.loc.rot * dir * d;
 							space
 								.trace_iter(Ray { pos, dir })
@@ -145,51 +146,6 @@ 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()
-}
-
-/// Converts a position and a rotation to a [Location]. Only the X direction is preserved from `rot` to ensure the resulting Location describes an orthonormal coordinate system.
-fn put_object(space: &impl Metric, pos: Vec2, rot: Mat2) -> Location {
-	let metric_sqrt = space.sqrt_at(pos);
-	let metric_inv_sqrt = space.sqrt_at(pos).inverse();
-	let rot = metric_inv_sqrt * (metric_sqrt * rot).orthonormalize();
-	Location { pos, rot }
-}
-
-#[test]
-fn test_put_object() {
-	use approx::assert_abs_diff_eq;
-
-	let ε = 1e-5;
-	let m = riemann::samples::ScaledMetric {
-		scale: vec2(3., 4.),
-	};
-
-	let loc = put_object(&m, vec2(1., 2.), mat2(vec2(1., 0.), vec2(0., 1.)));
-	assert_eq!(loc.pos, vec2(1., 2.));
-	assert_abs_diff_eq!(loc.rot * vec2(1., 0.), vec2(1. / 3., 0.), epsilon = ε);
-	assert_abs_diff_eq!(loc.rot * vec2(0., 1.), vec2(0., 1. / 4.), epsilon = ε);
-
-	let loc = put_object(&m, vec2(1., 2.), mat2(vec2(0., 1.), vec2(-1., 0.)));
-	assert_eq!(loc.pos, vec2(1., 2.));
-	assert_abs_diff_eq!(loc.rot * vec2(1., 0.), vec2(0., 1. / 4.), epsilon = ε);
-	assert_abs_diff_eq!(loc.rot * vec2(0., 1.), vec2(-1. / 3., 0.), epsilon = ε);
-
-	let c = 0.5 * std::f32::consts::SQRT_2;
-	let loc = put_object(&m, vec2(1., 2.), mat2(vec2(c, c), vec2(-c, c)));
-	assert_eq!(loc.pos, vec2(1., 2.));
-	assert_abs_diff_eq!(loc.rot * vec2(1., 0.), vec2(1. / 5., 1. / 5.), epsilon = ε);
-	assert_abs_diff_eq!(
-		loc.rot * vec2(0., 1.),
-		vec2(-4. / 15., 3. / 20.),
-		epsilon = ε
-	);
-}
-
 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);
@@ -197,8 +153,8 @@ fn draw_cross(gc: &mut Vec<Draw>, pos: Vec2, r: f32) {
 	gc.line_to(pos.x + r, pos.y - r);
 }
 
-fn draw_ray_2(gc: &mut Vec<Draw>, space: &Space, camera: Location, dir: Vec2) {
+fn draw_ray_2(gc: &mut Vec<Draw>, space: &Space, camera: Location, dir: Vec3) {
-	let pos = vec2(0., 0.);
+	let pos = vec3(0., 0., 0.);
 	let (hits, path) = space.trace_dbg(camera, Ray { pos, dir });
 	let hits2 = space.trace(camera, Ray { pos, dir });
 	for (a, b) in hits.into_iter().zip(hits2.into_iter()) {
@@ -210,20 +166,20 @@ fn draw_ray_2(gc: &mut Vec<Draw>, space: &Space, camera: Location, dir: Vec2) {
 	gc.new_path();
 	gc.move_to(pos.x, pos.y);
 	for pt in &path.points[1..] {
-		gc.line_to(pt.x, pt.y);
+		gc.line_to(pt.pos.x, pt.pos.y);
 	}
 	let end_pos = *path
 		.points
 		.last()
 		.expect("the starting point is always in the path");
-	let dir_pos = end_pos + 1000.0 * path.end_dir;
+	let dir_pos = end_pos.forward(1000. / DT).pos;
 	gc.line_to(dir_pos.x, dir_pos.y);
 	gc.stroke();
 }
 
 fn draw_fan_2(gc: &mut Vec<Draw>, space: &Space, camera: Location, spread: f32) {
 	for y in itertools_num::linspace(-spread, spread, 101) {
-		draw_ray_2(gc, space, camera, vec2(1., y));
+		draw_ray_2(gc, space, camera, vec3(1., y, 0.));
 	}
 }
 
@@ -234,10 +190,14 @@ fn draw_track(gc: &mut Vec<Draw>, space: &Space, start: Vec2, dir: Vec2) {
 	// let dir = space.tube.globalize(start, dir);
 	// let v = space.tube.normalize(start, dir);
 	let mut loc = Location {
-		pos: start,
+		pos: vec3(start.x, start.y, 0.),
-		rot: mat2(dir, vec2(-dir.y, dir.x)),
+		rot: mat3(
+			vec3(dir.x, dir.y, 0.),
+			vec3(-dir.y, dir.x, 0.),
+			vec3(0., 0., 1.),
+		),
 	};
-	let v = vec2(1.0, 0.0);
+	let v = vec3(1., 0., 0.);
 	let mut draw = |loc: &Location| {
 		let p = loc.pos;
 		let ax = p + loc.rot.x_axis * SCALE;

src/bin/flat/riemann.rs

@@ -1,281 +0,0 @@
-use glam::*;
-
-#[derive(Copy, Clone, Debug, PartialEq)]
-pub struct Decomp2 {
-	pub ortho: Mat2,
-	pub diag: Vec2,
-}
-
-impl Decomp2 {
-	fn square(&self) -> Self {
-		Self {
-			ortho: self.ortho,
-			diag: self.diag * self.diag,
-		}
-	}
-
-	pub fn inverse(&self) -> Self {
-		Self {
-			ortho: self.ortho,
-			diag: Vec2::splat(1.0) / self.diag,
-		}
-	}
-}
-
-impl<T> std::ops::Mul<T> for Decomp2
-where
-	Mat2: std::ops::Mul<T>,
-{
-	type Output = <Mat2 as std::ops::Mul<T>>::Output;
-
-	fn mul(self, rhs: T) -> Self::Output {
-		Mat2::from(self) * rhs
-	}
-}
-
-impl From<Decomp2> for Mat2 {
-	fn from(value: Decomp2) -> Self {
-		value.ortho.transpose() * Mat2::from_diagonal(value.diag) * value.ortho
-	}
-}
-
-pub type Tens2 = [Mat2; 2];
-
-pub trait Metric {
-	fn sqrt_at(&self, pos: Vec2) -> Decomp2;
-
-	fn at(&self, pos: Vec2) -> Mat2 {
-		self.sqrt_at(pos).square().into()
-	}
-
-	fn inverse_at(&self, pos: Vec2) -> Mat2 {
-		self.sqrt_at(pos).square().inverse().into()
-	}
-
-	fn part_derivs_at(&self, pos: Vec2) -> Tens2 {
-		part_deriv(|p| self.at(p), pos, 1.0 / 1024.0) // division by such eps is exact which is good for overall precision
-	}
-
-	fn vec_length_at(&self, at: Vec2, v: Vec2) -> f32 {
-		v.dot(self.at(at) * v).sqrt()
-	}
-
-	fn normalize_vec_at(&self, at: Vec2, v: Vec2) -> Vec2 {
-		v / self.vec_length_at(at, v)
-	}
-}
-
-pub struct TraceIter<'a, M: Metric> {
-	space: &'a M,
-	p: Vec2,
-	v: Vec2,
-}
-
-impl<'a, M: Metric> Iterator for TraceIter<'a, M> {
-	type Item = Vec2;
-
-	fn next(&mut self) -> Option<Self::Item> {
-		let a: Vec2 = -contract2(krist(self.space, self.p), self.v);
-		self.v += a;
-		self.p += self.v;
-		Some(self.p)
-	}
-}
-
-pub fn trace_iter<M: Metric>(space: &M, base: Vec2, dir: Vec2, dt: f32) -> TraceIter<M> {
-	TraceIter {
-		space,
-		p: base,
-		v: dt * space.normalize_vec_at(base, dir),
-	}
-}
-
-pub fn krist(space: &impl Metric, pos: Vec2) -> Tens2 {
-	// Γ^i_k_l = .5 * g^i^m * (g_m_k,l + g_m_l,k - g_k_l,m)
-	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::<f32>()
-	})
-}
-
-fn dir_deriv(f: impl Fn(Vec2) -> Mat2, pos: Vec2, delta: Vec2) -> Mat2 {
-	(f(pos + delta) - f(pos - delta)) / (2.0 * delta.length())
-}
-
-fn part_deriv(f: impl Fn(Vec2) -> Mat2, pos: Vec2, eps: f32) -> Tens2 {
-	[
-		dir_deriv(&f, pos, vec2(eps, 0.0)),
-		dir_deriv(&f, pos, vec2(0.0, eps)),
-	]
-}
-
-/// Сворачивает тензор t с вектором u
-pub fn contract(t: Tens2, u: Vec2) -> Mat2 {
-	mat2(t[0] * u, t[1] * u).transpose()
-}
-
-/// Сворачивает тензор t с вектором v дважды, по второму и третьему индексам.
-pub fn contract2(t: Tens2, v: Vec2) -> Vec2 {
-	contract(t, v) * v
-}
-
-fn make_vec2(f: impl Fn(usize) -> f32) -> Vec2 {
-	Vec2::from_array(std::array::from_fn(|i| f(i)))
-}
-
-fn make_mat2(f: impl Fn(usize, usize) -> f32) -> Mat2 {
-	Mat2::from_cols_array_2d(&std::array::from_fn(|i| std::array::from_fn(|j| f(i, j))))
-}
-
-fn make_tens2(f: impl Fn(usize, usize, usize) -> f32) -> Tens2 {
-	std::array::from_fn(|i| make_mat2(|j, k| f(i, j, k)))
-}
-
-#[test]
-fn m2() {
-	let m = make_mat2(|i, j| (i + 2 * j) as f32);
-	assert_eq!(m.col(0)[0], 0.0);
-	assert_eq!(m.col(1)[0], 1.0);
-	assert_eq!(m.col(0)[1], 2.0);
-	assert_eq!(m.col(1)[1], 3.0);
-}
-
-#[test]
-fn t2() {
-	let t = make_tens2(|i, j, k| (i + 2 * j + 4 * k) as f32);
-	assert_eq!(t[0].col(0)[0], 0.0);
-	assert_eq!(t[1].col(0)[0], 1.0);
-	assert_eq!(t[0].col(1)[0], 2.0);
-	assert_eq!(t[1].col(1)[0], 3.0);
-	assert_eq!(t[0].col(0)[1], 4.0);
-	assert_eq!(t[1].col(0)[1], 5.0);
-	assert_eq!(t[0].col(1)[1], 6.0);
-	assert_eq!(t[1].col(1)[1], 7.0);
-}
-
-pub mod samples {
-	use glam::{Mat2, Vec2};
-
-	use super::{Decomp2, Metric};
-
-	pub struct ScaledMetric {
-		/// Specifies unit size in each cardinal direction. E.g. with scale=(2, 3), vector (1, 0) has length 2 while a unit vector with the same direction is (1/2, 0).
-		pub scale: Vec2,
-	}
-
-	impl Metric for ScaledMetric {
-		fn sqrt_at(&self, _pos: Vec2) -> Decomp2 {
-			Decomp2 {
-				diag: self.scale,
-				ortho: Mat2::IDENTITY,
-			}
-		}
-	}
-}
-
-#[cfg(test)]
-mod tests {
-	use super::*;
-	use approx::assert_abs_diff_eq;
-
-	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.])
-		);
-		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.)
-		);
-	}
-
-	#[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,479 +0,0 @@
-use glam::{vec2, Mat2, Vec2};
-
-use crate::riemann::Metric;
-use crate::types::{Location, Ray};
-
-use super::{Rect, Tube};
-
-pub trait FlatCoordinateSystem<T> {
-	fn flat_to_global(&self, v: T) -> T;
-	fn global_to_flat(&self, v: T) -> T;
-}
-
-pub trait FlatRegion:
-	FlatCoordinateSystem<Vec2> + FlatCoordinateSystem<Ray> + FlatCoordinateSystem<Location>
-{
-	// Измеряет расстояние до выхода за пределы области вдоль луча ray. Луч задаётся в плоской СК.
-	fn distance_to_boundary(&self, _ray: Ray) -> Option<f32> {
-		None
-	}
-}
-
-trait MetricCS: FlatCoordinateSystem<Vec2> {
-	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()
-	}
-	fn global_to_flat_tfm_at(&self, pos: Vec2) -> Mat2 {
-		self.global_metric().sqrt_at(pos).into()
-	}
-}
-
-impl<T: FlatCoordinateSystem<Vec2> + MetricCS> FlatCoordinateSystem<Ray> for T {
-	fn flat_to_global(&self, ray: Ray) -> Ray {
-		Ray {
-			pos: self.flat_to_global(ray.pos),
-			dir: self.flat_to_global_tfm_at(ray.pos) * ray.dir,
-		}
-	}
-
-	fn global_to_flat(&self, ray: Ray) -> Ray {
-		Ray {
-			pos: self.global_to_flat(ray.pos),
-			dir: self.global_to_flat_tfm_at(ray.pos) * ray.dir,
-		}
-	}
-}
-
-impl<T: FlatCoordinateSystem<Vec2> + MetricCS> FlatCoordinateSystem<Location> for T {
-	fn flat_to_global(&self, loc: Location) -> Location {
-		Location {
-			pos: self.flat_to_global(loc.pos),
-			rot: self.flat_to_global_tfm_at(loc.pos) * loc.rot,
-		}
-	}
-
-	fn global_to_flat(&self, loc: Location) -> Location {
-		Location {
-			pos: self.global_to_flat(loc.pos), // в плоской СК для Inner или её продолжении на Outer
-			rot: self.global_to_flat_tfm_at(loc.pos) * loc.rot,
-		}
-	}
-}
-
-pub struct InnerCS(pub Tube);
-
-impl MetricCS for InnerCS {
-	fn global_metric(&self) -> &impl Metric {
-		&self.0
-	}
-}
-
-impl FlatCoordinateSystem<Vec2> for InnerCS {
-	fn flat_to_global(&self, pos: Vec2) -> Vec2 {
-		vec2(pos.x, self.0.y(pos.y))
-	}
-
-	// Работает только при |pos.x| ≤ inner_radius или |pos.y| ≥ external_halflength.
-	fn global_to_flat(&self, pos: Vec2) -> Vec2 {
-		vec2(pos.x, self.0.v(pos.y))
-	}
-}
-
-impl FlatRegion for InnerCS {
-	fn distance_to_boundary(&self, ray: Ray) -> Option<f32> {
-		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 FlatCoordinateSystem<Vec2> 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),
-		};
-		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));
-			vec2(x, y)
-		} else {
-			pos
-		}
-	}
-
-	fn global_to_flat(&self, pos: Vec2) -> Vec2 {
-		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);
-			vec2(u, v) // в плоском продолжении СК Outer на область Inner
-		} else {
-			pos
-		}
-	}
-}
-
-impl FlatRegion for OuterCS {
-	fn distance_to_boundary(&self, ray: Ray) -> Option<f32> {
-		Rect {
-			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 glam::{mat2, vec2, Mat2, Vec2};
-	use itertools_num::linspace;
-
-	use crate::riemann::samples;
-
-	use super::*;
-
-	#[test]
-	fn uniform_scaled_metric() {
-		struct Scaled(samples::ScaledMetric, Vec2);
-		impl FlatCoordinateSystem<Vec2> 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
-			}
-		}
-		impl MetricCS for Scaled {
-			fn global_metric(&self) -> &impl Metric {
-				&self.0
-			}
-		}
-		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.))
-			}
-		);
-	}
-
-	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 {
-			($at: expr, $left: expr, $right: expr) => {
-				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)
-				);
-			};
-		}
-		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 = ε);
-			}
-		}
-		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
-				);
-				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
-				);
-			}
-		}
-		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
-				);
-				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
-				);
-			}
-		}
-	}
-
-	#[test]
-	fn test_mapper_inner() {
-		let mapper = InnerCS(Tube {
-			inner_radius: 30.0,
-			outer_radius: 50.0,
-			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]
-	fn test_mapper_outer() {
-		let mapper = OuterCS(Tube {
-			inner_radius: 30.0,
-			outer_radius: 50.0,
-			internal_halflength: 100.0,
-			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)),
-		);
-		// 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
-				);
-			}
-		}
-		// 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;
-				assert_eq!(np, vec2(-pp.x, pp.y));
-				assert_eq!(pn, vec2(pp.x, -pp.y));
-				assert_eq!(nn, vec2(-pp.x, -pp.y));
-			}
-		}
-		// 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
-				);
-			}
-		}
-		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
-				);
-			}
-		}
-		// 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;
-				assert!(v > 200.0);
-				assert!(v > y);
-			}
-		}
-	}
-}

src/bin/wireframe/main.rs

@@ -0,0 +1,570 @@
+use std::{iter, mem, time::Instant};
+
+use glam::{mat4, vec2, vec3, vec4, Mat4, Vec3};
+use wgpu::{util::DeviceExt, ShaderStages};
+use winit::{
+	event::*,
+	event_loop::EventLoop,
+	keyboard::{KeyCode, PhysicalKey},
+	window::{Window, WindowBuilder},
+};
+
+mod scene;
+
+// The coordinate system:
+// * X: forward
+// * Y: left
+// * Z: up
+
+#[repr(C)]
+#[derive(Copy, Clone, Debug, bytemuck::Pod, bytemuck::Zeroable)]
+struct Vertex {
+	position: [f32; 3],
+	tangent: [f32; 3],
+}
+
+struct Wireframe {
+	color: Vec3,
+	data: wgpu::Buffer,
+	size: u32,
+}
+
+fn prepare_scene(device: &wgpu::Device) -> Vec<Wireframe> {
+	scene::build()
+		.into_iter()
+		.map(|line| {
+			let color = line.color;
+			let data: Vec<Vertex> = line
+				.pts
+				.into_iter()
+				.map(|r| Vertex {
+					position: r.pos.to_array(),
+					tangent: r.dir.to_array(),
+				})
+				.collect();
+			let size = data.len() as u32;
+			let data = device.create_buffer_init(&wgpu::util::BufferInitDescriptor {
+				label: Some("Vertex Buffer"),
+				contents: bytemuck::cast_slice(&data),
+				usage: wgpu::BufferUsages::VERTEX,
+			});
+			Wireframe { color, data, size }
+		})
+		.collect()
+}
+
+mod camctl {
+	use glam::{vec3, Mat4, Quat, Vec3};
+
+	pub struct CameraLocation {
+		pos: Vec3,
+		rot: Quat,
+	}
+
+	impl CameraLocation {
+		pub fn new() -> CameraLocation {
+			let rot = Quat::from_euler(glam::EulerRot::ZYX, std::f32::consts::FRAC_PI_4, 0., 0.);
+			let pos = rot * vec3(-200., 0., 50.);
+			CameraLocation { pos, rot }
+		}
+
+		pub fn view_mtx(&self) -> Mat4 {
+			Mat4::from_quat(self.rot.inverse()) * Mat4::from_translation(-self.pos)
+		}
+
+		pub fn move_rel(&mut self, offset: Vec3) {
+			self.pos += self.rot * offset;
+		}
+
+		pub fn rotate_rel_ypr(&mut self, ypr: Vec3) {
+			self.rotate_rel_quat(Quat::from_euler(glam::EulerRot::XYZ, ypr.x, ypr.y, ypr.z));
+		}
+
+		pub fn rotate_rel_quat(&mut self, rot: Quat) {
+			self.rot *= rot;
+		}
+	}
+}
+
+mod keyctl {
+	use std::{collections::HashSet, iter::Sum};
+	use winit::{event::ElementState, keyboard::PhysicalKey};
+
+	pub struct Keyboard {
+		pressed: HashSet<PhysicalKey>,
+	}
+
+	impl Keyboard {
+		pub fn new() -> Self {
+			Keyboard {
+				pressed: Default::default(),
+			}
+		}
+
+		pub fn is_pressed(&self, key: PhysicalKey) -> bool {
+			self.pressed.contains(&key)
+		}
+
+		pub fn set_key_state(&mut self, key: PhysicalKey, state: ElementState) {
+			match state {
+				ElementState::Pressed => self.pressed.insert(key),
+				ElementState::Released => self.pressed.remove(&key),
+			};
+		}
+
+		pub fn control<T: Copy + Sum>(&self, keymap: &[(PhysicalKey, T)]) -> T {
+			keymap
+				.iter()
+				.copied()
+				.filter_map(|(key, ctl)| self.is_pressed(key).then_some(ctl))
+				.sum()
+		}
+	}
+}
+
+static KEYS_MOVE: &'static [(PhysicalKey, Vec3)] = &[
+	(PhysicalKey::Code(KeyCode::KeyW), vec3(1., 0., 0.)),
+	(PhysicalKey::Code(KeyCode::KeyS), vec3(-1., 0., 0.)),
+	(PhysicalKey::Code(KeyCode::KeyA), vec3(0., 1., 0.)),
+	(PhysicalKey::Code(KeyCode::KeyD), vec3(0., -1., 0.)),
+	(PhysicalKey::Code(KeyCode::Space), vec3(0., 0., 1.)),
+	(PhysicalKey::Code(KeyCode::ShiftLeft), vec3(0., 0., -1.)),
+];
+
+static KEYS_ROTATE: &'static [(PhysicalKey, Vec3)] = &[
+	(PhysicalKey::Code(KeyCode::Numpad9), vec3(1., 0., 0.)),
+	(PhysicalKey::Code(KeyCode::Numpad7), vec3(-1., 0., 0.)),
+	(PhysicalKey::Code(KeyCode::Numpad5), vec3(0., 1., 0.)),
+	(PhysicalKey::Code(KeyCode::Numpad8), vec3(0., -1., 0.)),
+	(PhysicalKey::Code(KeyCode::Numpad4), vec3(0., 0., 1.)),
+	(PhysicalKey::Code(KeyCode::Numpad6), vec3(0., 0., -1.)),
+];
+
+#[repr(C)]
+#[derive(Copy, Clone, bytemuck::Pod, bytemuck::Zeroable)]
+struct CameraUniform {
+	mvp: [[f32; 4]; 4],
+	scale: [f32; 2],
+	pad: [u32; 2],
+}
+
+#[repr(C)]
+#[derive(Copy, Clone, bytemuck::Pod, bytemuck::Zeroable)]
+struct LineUniform {
+	color: [f32; 3],
+	_pad: f32,
+}
+
+struct State<'a> {
+	surface: wgpu::Surface<'a>,
+	device: wgpu::Device,
+	queue: wgpu::Queue,
+	config: wgpu::SurfaceConfiguration,
+	size: winit::dpi::PhysicalSize<u32>,
+	render_pipeline: wgpu::RenderPipeline,
+
+	kbd: keyctl::Keyboard,
+	cam: camctl::CameraLocation,
+	t1: Instant,
+
+	camera_buffer: wgpu::Buffer,
+	camera_bind_group: wgpu::BindGroup,
+
+	scene: Vec<Wireframe>,
+
+	window: &'a Window,
+}
+
+impl<'a> State<'a> {
+	async fn new(window: &'a Window) -> State<'a> {
+		let size = window.inner_size();
+
+		let instance = wgpu::Instance::new(wgpu::InstanceDescriptor {
+			backends: wgpu::Backends::PRIMARY,
+			..Default::default()
+		});
+		let surface = instance.create_surface(window).unwrap();
+		let adapter = instance
+			.request_adapter(&wgpu::RequestAdapterOptions {
+				power_preference: wgpu::PowerPreference::default(),
+				compatible_surface: Some(&surface),
+				force_fallback_adapter: false,
+			})
+			.await
+			.unwrap();
+		let (device, queue) = adapter
+			.request_device(
+				&wgpu::DeviceDescriptor {
+					label: None,
+					required_features: wgpu::Features::PUSH_CONSTANTS,
+					required_limits: wgpu::Limits {
+						max_push_constant_size: mem::size_of::<LineUniform>() as u32,
+						..wgpu::Limits::default()
+					},
+					memory_hints: Default::default(),
+				},
+				None, // Trace path
+			)
+			.await
+			.unwrap();
+
+		let surface_caps = surface.get_capabilities(&adapter);
+		let surface_format = surface_caps
+			.formats
+			.iter()
+			.copied()
+			.find(|f| !f.is_srgb())
+			.unwrap_or(surface_caps.formats[0]);
+		let config = wgpu::SurfaceConfiguration {
+			usage: wgpu::TextureUsages::RENDER_ATTACHMENT,
+			format: surface_format,
+			width: size.width,
+			height: size.height,
+			present_mode: surface_caps.present_modes[0],
+			alpha_mode: surface_caps.alpha_modes[0],
+			view_formats: vec![],
+			desired_maximum_frame_latency: 2,
+		};
+
+		let kbd = keyctl::Keyboard::new();
+		let cam = camctl::CameraLocation::new();
+		let t1 = Instant::now();
+
+		let camera_buffer = device.create_buffer(&wgpu::BufferDescriptor {
+			label: Some("Camera Buffer"),
+			usage: wgpu::BufferUsages::UNIFORM | wgpu::BufferUsages::COPY_DST,
+			size: mem::size_of::<CameraUniform>() as u64,
+			mapped_at_creation: false,
+		});
+		let camera_bind_group_layout =
+			device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor {
+				entries: &[wgpu::BindGroupLayoutEntry {
+					binding: 0,
+					visibility: wgpu::ShaderStages::VERTEX,
+					ty: wgpu::BindingType::Buffer {
+						ty: wgpu::BufferBindingType::Uniform,
+						has_dynamic_offset: false,
+						min_binding_size: None,
+					},
+					count: None,
+				}],
+				label: Some("camera_bind_group_layout"),
+			});
+		let camera_bind_group = device.create_bind_group(&wgpu::BindGroupDescriptor {
+			layout: &camera_bind_group_layout,
+			entries: &[wgpu::BindGroupEntry {
+				binding: 0,
+				resource: camera_buffer.as_entire_binding(),
+			}],
+			label: Some("camera_bind_group"),
+		});
+
+		let shader = device.create_shader_module(wgpu::ShaderModuleDescriptor {
+			label: Some("Shader"),
+			source: wgpu::ShaderSource::Wgsl(include_str!("ray.wgsl").into()),
+		});
+
+		let line_push_constant_range = wgpu::PushConstantRange {
+			stages: ShaderStages::VERTEX,
+			range: 0..mem::size_of::<LineUniform>() as u32,
+		};
+
+		let render_pipeline_layout =
+			device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor {
+				label: Some("Render Pipeline Layout"),
+				bind_group_layouts: &[&camera_bind_group_layout],
+				push_constant_ranges: &[line_push_constant_range],
+			});
+
+		let render_pipeline = device.create_render_pipeline(&wgpu::RenderPipelineDescriptor {
+			label: Some("Render Pipeline"),
+			layout: Some(&render_pipeline_layout),
+			vertex: wgpu::VertexState {
+				module: &shader,
+				entry_point: "vs_main",
+				buffers: &[wgpu::VertexBufferLayout {
+					array_stride: mem::size_of::<Vertex>() as wgpu::BufferAddress,
+					step_mode: wgpu::VertexStepMode::Instance,
+					attributes: &[
+						wgpu::VertexAttribute {
+							offset: mem::offset_of!(Vertex, position) as u64,
+							shader_location: 0,
+							format: wgpu::VertexFormat::Float32x3,
+						},
+						wgpu::VertexAttribute {
+							offset: mem::offset_of!(Vertex, tangent) as u64,
+							shader_location: 1,
+							format: wgpu::VertexFormat::Float32x3,
+						},
+						wgpu::VertexAttribute {
+							offset: (mem::size_of::<Vertex>() + mem::offset_of!(Vertex, position))
+								as u64,
+							shader_location: 2,
+							format: wgpu::VertexFormat::Float32x3,
+						},
+						wgpu::VertexAttribute {
+							offset: (mem::size_of::<Vertex>() + mem::offset_of!(Vertex, tangent))
+								as u64,
+							shader_location: 3,
+							format: wgpu::VertexFormat::Float32x3,
+						},
+					],
+				}],
+				compilation_options: Default::default(),
+			},
+			fragment: Some(wgpu::FragmentState {
+				module: &shader,
+				entry_point: "fs_main",
+				targets: &[Some(wgpu::ColorTargetState {
+					format: config.format,
+					blend: Some(wgpu::BlendState {
+						color: wgpu::BlendComponent::OVER,
+						alpha: wgpu::BlendComponent::OVER,
+					}),
+					write_mask: wgpu::ColorWrites::ALL,
+				})],
+				compilation_options: Default::default(),
+			}),
+			primitive: wgpu::PrimitiveState {
+				topology: wgpu::PrimitiveTopology::TriangleStrip,
+				..Default::default()
+			},
+			depth_stencil: None,
+			multisample: wgpu::MultisampleState {
+				count: 1,
+				mask: !0,
+				alpha_to_coverage_enabled: false,
+			},
+			multiview: None,
+			cache: None,
+		});
+
+		let scene = prepare_scene(&device);
+
+		Self {
+			surface,
+			device,
+			queue,
+			config,
+			size,
+			render_pipeline,
+			kbd,
+			cam,
+			t1,
+			scene,
+			camera_buffer,
+			camera_bind_group,
+			window,
+		}
+	}
+
+	pub fn window(&self) -> &Window {
+		&self.window
+	}
+
+	fn resize(&mut self, new_size: winit::dpi::PhysicalSize<u32>) {
+		if new_size.width > 0 && new_size.height > 0 {
+			self.size = new_size;
+			self.config.width = new_size.width;
+			self.config.height = new_size.height;
+			self.surface.configure(&self.device, &self.config);
+		}
+	}
+
+	fn update(&mut self) {
+		let dt = {
+			let t2 = Instant::now();
+			let dt = t2 - self.t1;
+			self.t1 = t2;
+			dt.as_secs_f32()
+		};
+		self.cam.move_rel(100. * dt * self.kbd.control(&KEYS_MOVE));
+		self.cam
+			.rotate_rel_ypr(2. * dt * self.kbd.control(&KEYS_ROTATE));
+
+		let size = vec2(self.config.width as f32, self.config.height as f32);
+		let size = size.normalize() * std::f32::consts::SQRT_2;
+		let proj = make_proj_matrix(vec3(size.x, size.y, 2.), (1., 4096.));
+
+		let my_to_gl = mat4(
+			vec4(0., 0., 1., 0.),
+			vec4(-1., 0., 0., 0.),
+			vec4(0., 1., 0., 0.),
+			vec4(0., 0., 0., 1.),
+		);
+		let view = my_to_gl * self.cam.view_mtx();
+		let mvp = proj * view;
+		let camera_uniform = CameraUniform {
+			mvp: mvp.to_cols_array_2d(),
+			scale: (1. / size).to_array(),
+			pad: [0; 2],
+		};
+		self.queue
+			.write_buffer(&self.camera_buffer, 0, bytemuck::bytes_of(&camera_uniform));
+	}
+
+	fn render(&mut self) -> Result<(), wgpu::SurfaceError> {
+		let output = self.surface.get_current_texture()?;
+		let view = output
+			.texture
+			.create_view(&wgpu::TextureViewDescriptor::default());
+
+		let mut encoder = self
+			.device
+			.create_command_encoder(&wgpu::CommandEncoderDescriptor {
+				label: Some("Render Encoder"),
+			});
+
+		{
+			let mut render_pass = encoder.begin_render_pass(&wgpu::RenderPassDescriptor {
+				label: Some("Render Pass"),
+				color_attachments: &[Some(wgpu::RenderPassColorAttachment {
+					view: &view,
+					resolve_target: None,
+					ops: wgpu::Operations {
+						load: wgpu::LoadOp::Clear(wgpu::Color {
+							r: 0.,
+							g: 0.,
+							b: 0.,
+							a: 1.,
+						}),
+						store: wgpu::StoreOp::Store,
+					},
+				})],
+				depth_stencil_attachment: None,
+				occlusion_query_set: None,
+				timestamp_writes: None,
+			});
+
+			render_pass.set_pipeline(&self.render_pipeline);
+			render_pass.set_bind_group(0, &self.camera_bind_group, &[]);
+			for wireframe in &self.scene {
+				let line = LineUniform {
+					color: wireframe.color.to_array(),
+					_pad: 0.,
+				};
+				render_pass.set_push_constants(ShaderStages::VERTEX, 0, bytemuck::bytes_of(&line));
+				render_pass.set_vertex_buffer(0, wireframe.data.slice(..));
+				render_pass.draw(0..4, 0..wireframe.size - 1);
+			}
+		}
+
+		self.queue.submit(iter::once(encoder.finish()));
+		output.present();
+
+		Ok(())
+	}
+}
+
+pub async fn run() {
+	let event_loop = EventLoop::new().unwrap();
+	let window = WindowBuilder::new()
+		.with_title("Refraction: Wireframe")
+		.build(&event_loop)
+		.unwrap();
+
+	// State::new uses async code, so we're going to wait for it to finish
+	let mut state = State::new(&window).await;
+	let mut surface_configured = false;
+
+	event_loop
+		.run(move |event, control_flow| {
+			match event {
+				Event::WindowEvent {
+					ref event,
+					window_id,
+				} if window_id == state.window().id() => {
+					match event {
+						WindowEvent::KeyboardInput {
+							device_id: _,
+							event,
+							is_synthetic: _,
+						} => {
+							state.kbd.set_key_state(event.physical_key, event.state);
+						}
+						WindowEvent::CloseRequested => control_flow.exit(),
+						WindowEvent::Resized(physical_size) => {
+							surface_configured = true;
+							state.resize(*physical_size);
+						}
+						WindowEvent::RedrawRequested => {
+							// This tells winit that we want another frame after this one
+							state.window().request_redraw();
+
+							if !surface_configured {
+								return;
+							}
+
+							state.update();
+							match state.render() {
+								Ok(_) => {}
+								// Reconfigure the surface if it's lost or outdated
+								Err(wgpu::SurfaceError::Lost | wgpu::SurfaceError::Outdated) => {
+									state.resize(state.size)
+								}
+								// The system is out of memory, we should probably quit
+								Err(wgpu::SurfaceError::OutOfMemory) => {
+									eprintln!("OutOfMemory");
+									control_flow.exit();
+								}
+
+								// This happens when the a frame takes too long to present
+								Err(wgpu::SurfaceError::Timeout) => {
+									eprintln!("Surface timeout")
+								}
+							}
+						}
+						_ => {}
+					}
+				}
+				_ => {}
+			}
+		})
+		.unwrap();
+}
+
+fn main() {
+	pollster::block_on(run());
+}
+
+/// Make a projection matrix, assuming input coordinates are (right, up, forward).
+///
+/// `corner` is a vector that will be mapped to (x=1, y=1) after the perspective division.
+/// `zrange` is the Z range that will be mapped to z∈[-1, 1]. It has no other effect. Both ends have to be positive though.
+fn make_proj_matrix(corner: Vec3, zrange: (f32, f32)) -> Mat4 {
+	let scale = 1.0 / corner;
+	let zspan = zrange.1 - zrange.0;
+	mat4(
+		scale.x * vec4(1., 0., 0., 0.),
+		scale.y * vec4(0., 1., 0., 0.),
+		scale.z * vec4(0., 0., (zrange.0 + zrange.1) / zspan, 1.),
+		scale.z * vec4(0., 0., -2. * zrange.0 * zrange.1 / zspan, 0.),
+	)
+}
+
+#[cfg(test)]
+mod tests {
+	use super::*;
+	use approx::assert_abs_diff_eq;
+	use glam::vec3;
+
+	#[test]
+	fn test_proj_matrix() {
+		let m = make_proj_matrix(vec3(2., 3., 4.), (0.5, 20.0));
+
+		let v = m * vec4(2., 3., 4., 1.);
+		assert_abs_diff_eq!(v.x / v.w, 1.0);
+		assert_abs_diff_eq!(v.y / v.w, 1.0);
+		assert!(-v.w < v.z && v.z < v.w, "z out of range in {v}");
+
+		let v = m * vec4(2., 3., 0.5, 1.);
+		assert_abs_diff_eq!(v.x / v.w, 8.0);
+		assert_abs_diff_eq!(v.y / v.w, 8.0);
+		assert_abs_diff_eq!(v.z / v.w, -1.0);
+
+		let v = m * vec4(2., 3., 20.0, 1.);
+		assert_abs_diff_eq!(v.x / v.w, 0.2);
+		assert_abs_diff_eq!(v.y / v.w, 0.2);
+		assert_abs_diff_eq!(v.z / v.w, 1.0);
+	}
+}

src/bin/wireframe/ray.wgsl

@@ -0,0 +1,60 @@
+struct CameraUniform {
+	mvp: mat4x4<f32>,
+	scale: vec2<f32>,
+}
+@group(0) @binding(0)
+var<uniform> camera: CameraUniform;
+
+struct LineUniform {
+	color: vec3<f32>,
+}
+const width = 0.1;
+var<push_constant> line: LineUniform;
+
+struct SegmentInput {
+	@location(0) a: vec3<f32>,
+	@location(1) ad: vec3<f32>,
+	@location(2) b: vec3<f32>,
+	@location(3) bd: vec3<f32>,
+}
+struct OffsetInput {
+	@builtin(vertex_index) idx: u32,
+}
+
+struct VertexOutput {
+	@builtin(position) clip_position: vec4<f32>,
+	@location(0) vertex_color: vec3<f32>,
+}
+
+@vertex
+fn vs_main(seg: SegmentInput, off: OffsetInput) -> VertexOutput {
+	var out: VertexOutput;
+	out.vertex_color = line.color;
+	var pt: vec3<f32>;
+	var dir: vec3<f32>;
+	switch (off.idx) {
+	case 0u: { pt = seg.a; dir = seg.ad; }
+	case 1u: { pt = seg.a; dir = seg.ad; }
+	case 2u: { pt = seg.b; dir = seg.bd; }
+	case 3u: { pt = seg.b; dir = seg.bd; }
+	default: {}
+	}
+	var sgn: f32;
+	switch (off.idx) {
+	case 0u: { sgn = -1.; }
+	case 1u: { sgn = 1.; }
+	case 2u: { sgn = -1.; }
+	case 3u: { sgn = 1.; }
+	default: {}
+	}
+	let pt_cs = camera.mvp * vec4(pt, 1.);
+	let dir_cs = camera.mvp * vec4(dir, 0.);
+	let normal_cs = camera.scale * normalize(vec2(-dir_cs.y, dir_cs.x));
+	out.clip_position = pt_cs + vec4(sgn * width * normal_cs, 0., 0.);
+	return out;
+}
+
+@fragment
+fn fs_main(in: VertexOutput) -> @location(0) vec4<f32> {
+	return 0.5 * vec4(in.vertex_color, 1.0);
+}

src/bin/wireframe/scene.rs

@@ -0,0 +1,163 @@
+use glam::*;
+use itertools::{chain, iproduct};
+
+use refraction::ifaces::{DebugTraceable, Traceable};
+use refraction::tube::metric::Tube;
+use refraction::tube::Space;
+use refraction::types::{Location, Object, Ray};
+use refraction::utils::put_object;
+
+pub enum Line {
+	Strip(Vec<Ray>),
+	Loop(Vec<Ray>),
+}
+
+pub struct FancyLine {
+	pub color: Vec3,
+	pub pts: Vec<Ray>,
+}
+
+fn paint(onto: &mut Vec<FancyLine>, color: Vec3, lines: Vec<Line>) {
+	onto.extend(lines.into_iter().map(move |line| FancyLine {
+		color,
+		pts: match line {
+			Line::Strip(pts) => pts,
+			Line::Loop(mut pts) => {
+				pts.push(*pts.first().unwrap());
+				pts
+			}
+		},
+	}))
+}
+
+fn draw_line(a: Vec3, b: Vec3) -> Line {
+	let dir = (b - a).normalize();
+	Line::Strip(vec![Ray { pos: a, dir }, Ray { pos: b, dir }])
+}
+
+fn draw_rect(center: Vec3, u: Vec3, v: Vec3) -> Vec<Line> {
+	let a = center - u - v;
+	let b = center + u - v;
+	let c = center + u + v;
+	let d = center - u + v;
+	vec![
+		draw_line(a, b),
+		draw_line(b, c),
+		draw_line(c, d),
+		draw_line(d, a),
+	]
+}
+
+fn draw_ellipse(center: Vec3, u: Vec3, v: Vec3) -> Line {
+	let segments = 47;
+	let step = 2. * std::f32::consts::PI / segments as f32;
+	Line::Loop(
+		(0..segments)
+			.map(|k| k as f32 * step)
+			.map(Vec2::from_angle)
+			.map(|d| Ray {
+				pos: center + d.x * u + d.y * v,
+				dir: -d.y * u + d.x * v,
+			})
+			.collect(),
+	)
+}
+
+pub fn build() -> Vec<FancyLine> {
+	let tube = Tube {
+		inner_radius: 30.0,
+		outer_radius: 50.0,
+		internal_halflength: 100.0,
+		external_halflength: 300.0,
+	};
+
+	let objs: Vec<_> = [-1.25, -1.00, -0.85, -0.50, 0.00, 0.40, 0.70, 0.95, 1.05]
+		.iter()
+		.enumerate()
+		.map(|(k, &y)| Object {
+			id: k as i32,
+			loc: put_object(
+				&tube,
+				vec3(0.0, y * tube.external_halflength, 0.0),
+				Mat3::from_mat2(Mat2::from_angle(y)),
+			),
+			r: 20.0,
+		})
+		.collect();
+	let space = Space { tube, objs };
+	let cam1 = put_object(&space.tube, vec3(-500., 0., 0.), Mat3::IDENTITY);
+	let cam2 = put_object(
+		&space.tube,
+		vec3(
+			-2.5 * tube.outer_radius,
+			1.25 * tube.external_halflength,
+			0.,
+		),
+		mat3(vec3(1., -1., 0.), vec3(1., 1., 0.), vec3(0., 0., 1.)),
+	);
+	let cam3 = put_object(
+		&space.tube,
+		vec3(
+			0.25 * tube.inner_radius,
+			0.25 * tube.external_halflength,
+			0.,
+		),
+		mat3(vec3(0., -1., 0.), vec3(1., 0., 0.), vec3(0., 0., 1.)),
+	);
+
+	let mut gc = vec![];
+	paint(&mut gc, vec3(0.8, 0.8, 0.8), tube.render());
+	paint(&mut gc, vec3(0.0, 0.8, 1.0), draw_fan_2(&space, cam3, 1.0));
+	paint(&mut gc, vec3(0.5, 1.0, 0.0), draw_fan_2(&space, cam2, 1.0));
+	paint(&mut gc, vec3(1.0, 0.5, 0.0), draw_fan_2(&space, cam1, 1.0));
+	gc
+}
+
+fn draw_ray_2(gc: &mut Vec<Line>, space: &Space, camera: Location, dir: Vec3) {
+	let pos = vec3(0., 0., 0.);
+	let (hits, path) = space.trace_dbg(camera, Ray { pos, dir });
+	let hits2 = space.trace(camera, Ray { pos, dir });
+	for (a, b) in hits.into_iter().zip(hits2.into_iter()) {
+		assert_eq!(a.id, b.id);
+		assert_eq!(a.pos, b.pos);
+		assert_eq!(a.rel, b.rel);
+	}
+
+	let mut pts = path.points;
+	let end_pos = *pts
+		.last()
+		.expect("the starting point is always in the path");
+	let dir_pos = end_pos.forward(10000.0);
+	pts.push(dir_pos);
+	gc.push(Line::Strip(pts));
+}
+
+fn draw_fan_2(space: &Space, camera: Location, spread: f32) -> Vec<Line> {
+	let mut gc = vec![];
+	for y in itertools_num::linspace(-spread, spread, 101) {
+		draw_ray_2(&mut gc, space, camera, vec3(1., y, 0.));
+	}
+	gc
+}
+
+trait Renderable {
+	fn render(&self) -> Vec<Line>;
+}
+
+impl Renderable for Tube {
+	fn render(&self) -> Vec<Line> {
+		let lines = 17;
+		let step = 2. * std::f32::consts::PI / lines as f32;
+		let r = 0.5 * (self.outer_radius + self.inner_radius);
+		let w = 0.5 * (self.outer_radius - self.inner_radius);
+		let l = vec3(0., self.external_halflength, 0.);
+		let along = (0..lines)
+			.map(|k| k as f32 * step)
+			.map(Vec2::from_angle)
+			.map(|d| vec3(d.x, 0., d.y))
+			.flat_map(|d| draw_rect(r * d, w * d, l));
+		let caps = iproduct!([self.inner_radius, self.outer_radius], [-l, l])
+			.map(|(r, l)| draw_ellipse(l, vec3(r, 0., 0.), vec3(0., 0., r)));
+		chain!(along, caps).collect()
+	}
+}

src/fns.rs

@@ -1,4 +1,6 @@
-use refraction::mathx::FloatExt2;
+//! Functions used to make metrics.
+
+use crate::mathx::FloatExt2;
 
 pub trait Limiter {
 	fn value(&self, x: f32) -> f32;

src/ifaces.rs

@@ -1,5 +1,4 @@
 use crate::types::{Hit, Location, Ray};
-use glam::Vec2;
 
 pub trait Traceable {
 	/// Traces a ray from a given starting point. `ray` is relative to the camera.
@@ -19,8 +18,7 @@ pub trait OptimizedTraceable: Traceable {
 }
 
 pub struct RayPath {
-	pub points: Vec<Vec2>,
+	pub points: Vec<Ray>,
-	pub end_dir: Vec2,
 }
 
 pub trait DebugTraceable: Traceable {

src/lib.rs

@@ -1,3 +1,11 @@
+mod fns;
+pub mod ifaces;
 pub mod mathx;
 pub mod mesh_loader;
 pub mod mesh_tracer;
+pub mod riemann;
+pub mod tube;
+pub mod types;
+pub mod utils;
+
+pub const DT: f32 = 0.1;

src/mathx.rs

@@ -1,4 +1,4 @@
-use glam::{FloatExt, Mat2, Mat3};
+use glam::{FloatExt, Mat2, Mat3, Vec2, Vec3};
 
 mod bounds {
 	pub trait Pair<T> {}
@@ -44,6 +44,104 @@ impl MatExt for Mat3 {
 	}
 }
 
+/// Represents a 2×2 matrix decomposed as O^T D O, where O is orthogonal and D is diagonal.
+///
+/// Not every matrix can be decomposed like this, only that of a symmetric bilinear function.
+#[derive(Copy, Clone, Debug, PartialEq)]
+pub struct Decomp2 {
+	/// The orthogonal part.
+	///
+	/// Using a non-orthogonal matrix will yield to incorrect results (but no UB).
+	pub ortho: Mat2,
+
+	/// The diagonal part.
+	pub diag: Vec2,
+}
+
+impl Decomp2 {
+	/// Computes the square of this, more efficiently than doing that with a matrix.
+	pub fn square(&self) -> Self {
+		Self {
+			ortho: self.ortho,
+			diag: self.diag * self.diag,
+		}
+	}
+
+	/// Computes the inverse of this, more efficiently than doing that with a matrix.
+	pub fn inverse(&self) -> Self {
+		Self {
+			ortho: self.ortho,
+			diag: Vec2::splat(1.0) / self.diag,
+		}
+	}
+}
+
+impl From<Decomp2> for Mat2 {
+	fn from(value: Decomp2) -> Self {
+		value.ortho.transpose() * Mat2::from_diagonal(value.diag) * value.ortho
+	}
+}
+
+impl<T> std::ops::Mul<T> for Decomp2
+where
+	Mat2: std::ops::Mul<T>,
+{
+	type Output = <Mat2 as std::ops::Mul<T>>::Output;
+
+	fn mul(self, rhs: T) -> Self::Output {
+		Mat2::from(self) * rhs
+	}
+}
+
+/// Represents a 3×3 matrix decomposed as O^T D O, where O is orthogonal and D is diagonal.
+///
+/// Not every matrix can be decomposed like this, only that of a symmetric bilinear function.
+#[derive(Copy, Clone, Debug, PartialEq)]
+pub struct Decomp3 {
+	/// The orthogonal part.
+	///
+	/// Using a non-orthogonal matrix will yield to incorrect results (but no UB).
+	pub ortho: Mat3,
+
+	/// The diagonal part.
+	pub diag: Vec3,
+}
+
+impl Decomp3 {
+	/// Computes the square of this, more efficiently than doing that with a matrix.
+	pub fn square(&self) -> Self {
+		Self {
+			ortho: self.ortho,
+			diag: self.diag * self.diag,
+		}
+	}
+
+	/// Computes the inverse of this, more efficiently than doing that with a matrix.
+	pub fn inverse(&self) -> Self {
+		Self {
+			ortho: self.ortho,
+			diag: Vec3::splat(1.0) / self.diag,
+		}
+	}
+}
+
+impl From<Decomp3> for Mat3 {
+	fn from(value: Decomp3) -> Self {
+		value.ortho.transpose() * Mat3::from_diagonal(value.diag) * value.ortho
+	}
+}
+
+impl<T> std::ops::Mul<T> for Decomp3
+where
+	Mat3: std::ops::Mul<T>,
+{
+	type Output = <Mat3 as std::ops::Mul<T>>::Output;
+
+	fn mul(self, rhs: T) -> Self::Output {
+		Mat3::from(self) * rhs
+	}
+}
+
 #[cfg(test)]
 mod tests {
 	use super::*;

src/riemann.rs

@@ -0,0 +1,248 @@
+use crate::mathx::Decomp3;
+use glam::*;
+
+pub type Tens3 = [Mat3; 3];
+
+pub trait Metric {
+	fn sqrt_at(&self, pos: Vec3) -> Decomp3;
+
+	fn at(&self, pos: Vec3) -> Mat3 {
+		self.sqrt_at(pos).square().into()
+	}
+
+	fn inverse_at(&self, pos: Vec3) -> Mat3 {
+		self.sqrt_at(pos).square().inverse().into()
+	}
+
+	fn part_derivs_at(&self, pos: Vec3) -> Tens3 {
+		part_deriv(|p| self.at(p), pos, 1.0 / 1024.0) // division by such eps is exact which is good for overall precision
+	}
+
+	fn vec_length_at(&self, at: Vec3, v: Vec3) -> f32 {
+		v.dot(self.at(at) * v).sqrt()
+	}
+
+	fn normalize_vec_at(&self, at: Vec3, v: Vec3) -> Vec3 {
+		v / self.vec_length_at(at, v)
+	}
+}
+
+pub struct TraceIter<'a, M: Metric> {
+	space: &'a M,
+	p: Vec3,
+	v: Vec3,
+}
+
+impl<'a, M: Metric> Iterator for TraceIter<'a, M> {
+	type Item = Vec3;
+
+	fn next(&mut self) -> Option<Self::Item> {
+		let a: Vec3 = -contract2(krist(self.space, self.p), self.v);
+		self.v += a;
+		self.p += self.v;
+		Some(self.p)
+	}
+}
+
+pub fn trace_iter<M: Metric>(space: &M, base: Vec3, dir: Vec3, dt: f32) -> TraceIter<M> {
+	TraceIter {
+		space,
+		p: base,
+		v: dt * space.normalize_vec_at(base, dir),
+	}
+}
+
+pub fn krist(space: &impl Metric, pos: Vec3) -> Tens3 {
+	// Γ^i_k_l = .5 * g^i^m * (g_m_k,l + g_m_l,k - g_k_l,m)
+	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_tens3(|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>()
+	})
+}
+
+fn dir_deriv(f: impl Fn(Vec3) -> Mat3, pos: Vec3, delta: Vec3) -> Mat3 {
+	(f(pos + delta) - f(pos - delta)) / (2.0 * delta.length())
+}
+
+fn part_deriv(f: impl Fn(Vec3) -> Mat3, pos: Vec3, eps: f32) -> Tens3 {
+	[
+		dir_deriv(&f, pos, vec3(eps, 0.0, 0.0)),
+		dir_deriv(&f, pos, vec3(0.0, eps, 0.0)),
+		dir_deriv(&f, pos, vec3(0.0, 0.0, eps)),
+	]
+}
+
+/// Сворачивает тензор t с вектором u
+pub fn contract(t: Tens3, u: Vec3) -> Mat3 {
+	mat3(t[0] * u, t[1] * u, t[2] * u).transpose()
+}
+
+/// Сворачивает тензор t с вектором v дважды, по второму и третьему индексам.
+pub fn contract2(t: Tens3, v: Vec3) -> Vec3 {
+	contract(t, v) * v
+}
+
+fn make_vec3(f: impl Fn(usize) -> f32) -> Vec3 {
+	Vec3::from_array(std::array::from_fn(|i| f(i)))
+}
+
+fn make_mat3(f: impl Fn(usize, usize) -> f32) -> Mat3 {
+	Mat3::from_cols_array_2d(&std::array::from_fn(|i| std::array::from_fn(|j| f(i, j))))
+}
+
+fn make_tens3(f: impl Fn(usize, usize, usize) -> f32) -> Tens3 {
+	std::array::from_fn(|i| make_mat3(|j, k| f(i, j, k)))
+}
+
+#[test]
+fn m3() {
+	let m = make_mat3(|i, j| (i + 2 * j) as f32);
+	assert_eq!(m.col(0)[0], 0.0);
+	assert_eq!(m.col(1)[0], 1.0);
+	assert_eq!(m.col(0)[1], 2.0);
+	assert_eq!(m.col(1)[1], 3.0);
+}
+
+#[test]
+fn t3() {
+	let t = make_tens3(|i, j, k| (i + 2 * j + 4 * k) as f32);
+	assert_eq!(t[0].col(0)[0], 0.0);
+	assert_eq!(t[1].col(0)[0], 1.0);
+	assert_eq!(t[0].col(1)[0], 2.0);
+	assert_eq!(t[1].col(1)[0], 3.0);
+	assert_eq!(t[0].col(0)[1], 4.0);
+	assert_eq!(t[1].col(0)[1], 5.0);
+	assert_eq!(t[0].col(1)[1], 6.0);
+	assert_eq!(t[1].col(1)[1], 7.0);
+}
+
+pub mod samples {
+	use glam::{Mat3, Vec3};
+
+	use super::{Decomp3, Metric};
+
+	pub struct ScaledMetric {
+		/// Specifies unit size in each cardinal direction. E.g. with scale=(2, 3), vector (1, 0) has length 2 while a unit vector with the same direction is (1/2, 0).
+		pub scale: Vec3,
+	}
+
+	impl Metric for ScaledMetric {
+		fn sqrt_at(&self, _pos: Vec3) -> Decomp3 {
+			Decomp3 {
+				diag: self.scale,
+				ortho: Mat3::IDENTITY,
+			}
+		}
+	}
+}
+
+#[cfg(test)]
+mod tests {
+	use super::*;
+	use approx::assert_abs_diff_eq;
+
+	use glam::{mat3, vec3, Mat3};
+	use rand::{Rng, SeedableRng};
+
+	#[test]
+	fn uniform_scaled_metric() {
+		let mut rng = rand_pcg::Pcg64Mcg::seed_from_u64(17);
+		let metric = samples::ScaledMetric {
+			scale: vec3(3., 4., 5.),
+		};
+		assert_eq!(
+			metric.sqrt_at(rng.gen()),
+			Decomp3 {
+				ortho: Mat3::IDENTITY,
+				diag: vec3(3., 4., 5.)
+			}
+		);
+		assert_eq!(
+			metric.at(rng.gen()),
+			Mat3::from_cols_array(&[9., 0., 0., 0., 16., 0., 0., 0., 25.])
+		);
+		assert_eq!(
+			metric.inverse_at(rng.gen()),
+			Mat3::from_cols_array(&[1. / 9., 0., 0., 0., 1. / 16., 0., 0., 0., 1. / 25.])
+		);
+		assert_eq!(
+			metric.part_derivs_at(rng.gen()),
+			[Mat3::ZERO, Mat3::ZERO, Mat3::ZERO]
+		);
+		assert_eq!(metric.vec_length_at(rng.gen(), vec3(1., 0., 0.)), 3.);
+		assert_eq!(metric.vec_length_at(rng.gen(), vec3(0., 1., 0.)), 4.);
+		assert_eq!(metric.vec_length_at(rng.gen(), vec3(0., 0., 1.)), 5.);
+		assert_eq!(metric.vec_length_at(rng.gen(), vec3(1., 1., 0.)), 5.);
+		assert_eq!(
+			metric.normalize_vec_at(rng.gen(), vec3(1., 0., 0.)),
+			vec3(1. / 3., 0., 0.)
+		);
+		assert_eq!(
+			metric.normalize_vec_at(rng.gen(), vec3(0., 1., 0.)),
+			vec3(0., 1. / 4., 0.)
+		);
+		assert_eq!(
+			metric.normalize_vec_at(rng.gen(), vec3(1., 1., 0.)),
+			vec3(1. / 5., 1. / 5., 0.)
+		);
+	}
+
+	#[test]
+	fn test_trace_iter() {
+		let metric = samples::ScaledMetric {
+			scale: vec3(2., 4., 3.),
+		};
+		assert_eq!(
+			trace_iter(&metric, vec3(3., 5., 0.), vec3(1., 0., 0.), 1.)
+				.nth(7)
+				.unwrap(),
+			vec3(7., 5., 0.)
+		);
+		assert_eq!(
+			trace_iter(&metric, vec3(3., 5., 0.), vec3(2., 0., 0.), 1.)
+				.nth(7)
+				.unwrap(),
+			vec3(7., 5., 0.)
+		);
+		assert_eq!(
+			trace_iter(&metric, vec3(3., 5., 0.), vec3(1., 0., 0.), 0.5)
+				.nth(7)
+				.unwrap(),
+			vec3(5., 5., 0.)
+		);
+		assert_eq!(
+			trace_iter(&metric, vec3(3., 5., 0.), vec3(0., 1., 0.), 1.)
+				.nth(9)
+				.unwrap(),
+			vec3(3., 7.5, 0.)
+		);
+		assert_eq!(
+			trace_iter(&metric, vec3(3., 5., 0.), vec3(0., 4., 0.), 1.)
+				.nth(9)
+				.unwrap(),
+			vec3(3., 7.5, 0.)
+		);
+		assert_eq!(
+			trace_iter(&metric, vec3(3., 5., 0.), vec3(0., 1., 0.), 0.5)
+				.nth(9)
+				.unwrap(),
+			vec3(3., 6.25, 0.)
+		);
+		assert_abs_diff_eq!(
+			trace_iter(
+				&metric,
+				vec3(3., 5., 0.),
+				vec3(0.5, 0.25, 0.),
+				std::f32::consts::SQRT_2
+			)
+			.nth(7)
+			.unwrap(),
+			vec3(7., 7., 0.),
+			epsilon = 1e-5
+		);
+	}
+}

src/tube/coords.rs

@@ -0,0 +1,513 @@
+use glam::{vec3, Mat3, Vec3};
+
+use crate::riemann::Metric;
+use crate::types::{Location, Ray};
+
+use super::{Rect, Tube};
+
+pub trait FlatCoordinateSystem<T> {
+	fn flat_to_global(&self, v: T) -> T;
+	fn global_to_flat(&self, v: T) -> T;
+}
+
+pub trait FlatRegion:
+	FlatCoordinateSystem<Vec3> + FlatCoordinateSystem<Ray> + FlatCoordinateSystem<Location>
+{
+	// Измеряет расстояние до выхода за пределы области вдоль луча ray. Луч задаётся в плоской СК.
+	fn distance_to_boundary(&self, _ray: Ray) -> Option<f32> {
+		None
+	}
+}
+
+trait MetricCS: FlatCoordinateSystem<Vec3> {
+	fn global_metric(&self) -> &impl Metric;
+	fn flat_to_global_tfm_at(&self, pos: Vec3) -> Mat3 {
+		self.global_metric()
+			.sqrt_at(self.flat_to_global(pos))
+			.inverse()
+			.into()
+	}
+	fn global_to_flat_tfm_at(&self, pos: Vec3) -> Mat3 {
+		self.global_metric().sqrt_at(pos).into()
+	}
+}
+
+impl<T: FlatCoordinateSystem<Vec3> + MetricCS> FlatCoordinateSystem<Ray> for T {
+	fn flat_to_global(&self, ray: Ray) -> Ray {
+		Ray {
+			pos: self.flat_to_global(ray.pos),
+			dir: self.flat_to_global_tfm_at(ray.pos) * ray.dir,
+		}
+	}
+
+	fn global_to_flat(&self, ray: Ray) -> Ray {
+		Ray {
+			pos: self.global_to_flat(ray.pos),
+			dir: self.global_to_flat_tfm_at(ray.pos) * ray.dir,
+		}
+	}
+}
+
+impl<T: FlatCoordinateSystem<Vec3> + MetricCS> FlatCoordinateSystem<Location> for T {
+	fn flat_to_global(&self, loc: Location) -> Location {
+		Location {
+			pos: self.flat_to_global(loc.pos),
+			rot: self.flat_to_global_tfm_at(loc.pos) * loc.rot,
+		}
+	}
+
+	fn global_to_flat(&self, loc: Location) -> Location {
+		Location {
+			pos: self.global_to_flat(loc.pos), // в плоской СК для Inner или её продолжении на Outer
+			rot: self.global_to_flat_tfm_at(loc.pos) * loc.rot,
+		}
+	}
+}
+
+pub struct InnerCS(pub Tube);
+
+impl MetricCS for InnerCS {
+	fn global_metric(&self) -> &impl Metric {
+		&self.0
+	}
+}
+
+impl FlatCoordinateSystem<Vec3> for InnerCS {
+	fn flat_to_global(&self, pos: Vec3) -> Vec3 {
+		vec3(pos.x, self.0.y(pos.y), pos.z)
+	}
+
+	// Работает только при |pos.x| ≤ inner_radius или |pos.y| ≥ external_halflength.
+	fn global_to_flat(&self, pos: Vec3) -> Vec3 {
+		vec3(pos.x, self.0.v(pos.y), pos.z)
+	}
+}
+
+impl FlatRegion for InnerCS {
+	fn distance_to_boundary(&self, ray: Ray) -> Option<f32> {
+		Rect {
+			size: vec3(
+				self.0.inner_radius,
+				self.0.internal_halflength,
+				self.0.inner_radius,
+			),
+		}
+		.trace_out_of(ray)
+	}
+}
+
+pub struct OuterCS(pub Tube);
+
+impl MetricCS for OuterCS {
+	fn global_metric(&self) -> &impl Metric {
+		&self.0
+	}
+}
+
+impl FlatCoordinateSystem<Vec3> for OuterCS {
+	fn flat_to_global(&self, pos: Vec3) -> Vec3 {
+		let inner = Rect {
+			size: vec3(
+				self.0.inner_radius + 1.0,
+				self.0.external_halflength,
+				self.0.inner_radius + 1.0,
+			),
+		};
+		if inner.is_inside(pos) {
+			let Vec3 { x, y: v, z } = pos;
+			let y = self
+				.0
+				.y(v - v.signum() * (self.0.external_halflength - self.0.internal_halflength));
+			vec3(x, y, z)
+		} else {
+			pos
+		}
+	}
+
+	fn global_to_flat(&self, pos: Vec3) -> Vec3 {
+		let inner = Rect {
+			size: vec3(
+				self.0.inner_radius + 1.0,
+				self.0.external_halflength,
+				self.0.inner_radius + 1.0,
+			),
+		};
+		if inner.is_inside(pos) {
+			let Vec3 { x: u, y, z: w } = pos; // в основной СК
+			let v = self.0.v(y)
+				+ y.signum() * (self.0.external_halflength - self.0.internal_halflength);
+			vec3(u, v, w) // в плоском продолжении СК Outer на область Inner
+		} else {
+			pos
+		}
+	}
+}
+
+impl FlatRegion for OuterCS {
+	fn distance_to_boundary(&self, ray: Ray) -> Option<f32> {
+		Rect {
+			size: vec3(
+				self.0.outer_radius,
+				self.0.external_halflength,
+				self.0.outer_radius,
+			),
+		}
+		.trace_into(ray)
+	}
+}
+
+#[cfg(test)]
+mod tests {
+	use approx::{assert_abs_diff_eq, AbsDiffEq};
+	use glam::{mat3, vec3, Mat3, Vec3};
+	use itertools_num::linspace;
+
+	use crate::riemann::samples;
+
+	use super::*;
+
+	#[test]
+	fn uniform_scaled_metric() {
+		struct Scaled(samples::ScaledMetric, Vec3);
+		impl FlatCoordinateSystem<Vec3> for Scaled {
+			fn flat_to_global(&self, pos: Vec3) -> Vec3 {
+				(pos - self.1) / self.0.scale
+			}
+			fn global_to_flat(&self, pos: Vec3) -> Vec3 {
+				pos * self.0.scale + self.1
+			}
+		}
+		impl MetricCS for Scaled {
+			fn global_metric(&self) -> &impl Metric {
+				&self.0
+			}
+		}
+		let cs = Scaled(
+			samples::ScaledMetric {
+				scale: vec3(3., 4., 5.),
+			},
+			vec3(2., 3., 7.),
+		);
+		assert_eq!(cs.global_to_flat(vec3(7., 3., 1.)), vec3(23., 15., 12.));
+		assert_eq!(cs.flat_to_global(vec3(8., 7., 17.)), vec3(2., 1., 2.));
+		assert_eq!(
+			cs.global_to_flat(Ray {
+				pos: vec3(7., 3., 0.),
+				dir: vec3(3., 2., 0.)
+			}),
+			Ray {
+				pos: vec3(23., 15., 7.),
+				dir: vec3(9., 8., 0.)
+			}
+		);
+		assert_eq!(
+			cs.flat_to_global(Ray {
+				pos: vec3(23., 15., 7.),
+				dir: vec3(9., 8., 0.)
+			}),
+			Ray {
+				pos: vec3(7., 3., 0.),
+				dir: vec3(3., 2., 0.)
+			}
+		);
+		assert_eq!(
+			cs.global_to_flat(Location {
+				pos: vec3(2., 1., 0.),
+				rot: mat3(vec3(0., 1., 0.), vec3(-1., 0., 0.), vec3(0., 0., 1.))
+			}),
+			Location {
+				pos: vec3(8., 7., 7.),
+				rot: mat3(vec3(0., 4., 0.), vec3(-3., 0., 0.), vec3(0., 0., 5.))
+			}
+		);
+		assert_eq!(
+			cs.flat_to_global(Location {
+				pos: vec3(2., 1., 7.),
+				rot: mat3(vec3(0., 1., 0.), vec3(-1., 0., 0.), vec3(0., 0., 1.))
+			}),
+			Location {
+				pos: vec3(0., -0.5, 0.),
+				rot: mat3(
+					vec3(0., 0.25, 0.),
+					vec3(-1. / 3., 0., 0.),
+					vec3(0., 0., 0.2)
+				)
+			}
+		);
+	}
+
+	fn test_flat_region(
+		region: &impl FlatRegion,
+		range_global: (Vec3, Vec3),
+		range_flat: (Vec3, Vec3),
+	) {
+		#[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!(
+					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: Vec3,
+			range_a: (Vec3, Vec3),
+			name_b: &str,
+			b: Vec3,
+			range_b: (Vec3, Vec3),
+		) {
+			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 = ε);
+			}
+		}
+		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) {
+				for z in linspace(range_global.0.z, range_global.1.z, 20) {
+					let pos_global = vec3(x, y, z);
+					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: Mat3::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: Mat3::IDENTITY
+							}))
+							.rot,
+						Mat3::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) {
+				for z in linspace(range_flat.0.z, range_flat.1.z, 20) {
+					let pos_flat = vec3(x, y, z);
+					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: Mat3::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: Mat3::IDENTITY
+							}))
+							.rot,
+						Mat3::IDENTITY
+					);
+				}
+			}
+		}
+	}
+
+	#[test]
+	fn test_mapper_inner() {
+		let mapper = InnerCS(Tube {
+			inner_radius: 30.0,
+			outer_radius: 50.0,
+			internal_halflength: 100.0,
+			external_halflength: 300.0,
+		});
+		test_flat_region(
+			&mapper,
+			(vec3(-30.0, -300.0, -30.0), vec3(30.0, 300.0, 30.0)),
+			(vec3(-30.0, -100.0, -30.0), vec3(30.0, 100.0, 30.0)),
+		);
+		test_flat_region(
+			&mapper,
+			(vec3(-60.0, -400.0, -60.0), vec3(60.0, -300.0, 60.0)),
+			(vec3(-60.0, -200.0, -60.0), vec3(60.0, -100.0, 60.0)),
+		);
+		test_flat_region(
+			&mapper,
+			(vec3(-60.0, 300.0, -60.0), vec3(60.0, 400.0, 60.0)),
+			(vec3(-60.0, 100.0, -60.0), vec3(60.0, 200.0, 60.0)),
+		);
+	}
+
+	#[test]
+	fn test_mapper_outer() {
+		let mapper = OuterCS(Tube {
+			inner_radius: 30.0,
+			outer_radius: 50.0,
+			internal_halflength: 100.0,
+			external_halflength: 300.0,
+		});
+		// TODO replace 200.20016 with something sane
+		test_flat_region(
+			&mapper,
+			(vec3(-30.0, -300.0, -30.0), vec3(30.0, -1.0, 30.0)),
+			(vec3(-30.0, -300.0, -30.0), vec3(30.0, -200.20016, 30.0)),
+		);
+		test_flat_region(
+			&mapper,
+			(vec3(-30.0, 1.0, -30.0), vec3(30.0, 300.0, 30.0)),
+			(vec3(-30.0, 200.20016, -30.0), vec3(30.0, 300.0, 30.0)),
+		);
+		test_flat_region(
+			&mapper,
+			(vec3(-60.0, -400.0, -60.0), vec3(60.0, -300.0, 60.0)),
+			(vec3(-60.0, -400.0, -60.0), vec3(60.0, -300.0, 60.0)),
+		);
+		test_flat_region(
+			&mapper,
+			(vec3(-60.0, 300.0, -60.0), vec3(60.0, 400.0, 60.0)),
+			(vec3(-60.0, 300.0, -60.0), vec3(60.0, 400.0, 60.0)),
+		);
+		// straight
+		for x in linspace(-60., 60., 20) {
+			for y in linspace(-320., 320., 20) {
+				for z in linspace(-60., 60., 20) {
+					assert_eq!(
+						mapper
+							.global_to_flat(Location {
+								pos: vec3(x, y, z),
+								rot: Mat3::IDENTITY
+							})
+							.pos
+							.x,
+						x
+					);
+				}
+			}
+		}
+		// symmetrical
+		for x in linspace(0., 60., 20) {
+			for y in linspace(0., 320., 20) {
+				for z in linspace(0., 60., 20) {
+					let pp = mapper
+						.global_to_flat(Location {
+							pos: vec3(x, y, z),
+							rot: Mat3::IDENTITY,
+						})
+						.pos;
+					let np = mapper
+						.global_to_flat(Location {
+							pos: vec3(-x, y, z),
+							rot: Mat3::IDENTITY,
+						})
+						.pos;
+					let pn = mapper
+						.global_to_flat(Location {
+							pos: vec3(x, -y, z),
+							rot: Mat3::IDENTITY,
+						})
+						.pos;
+					let nn = mapper
+						.global_to_flat(Location {
+							pos: vec3(-x, -y, z),
+							rot: Mat3::IDENTITY,
+						})
+						.pos;
+					assert_eq!(np, vec3(-pp.x, pp.y, pp.z));
+					assert_eq!(pn, vec3(pp.x, -pp.y, pp.z));
+					assert_eq!(nn, vec3(-pp.x, -pp.y, pp.z));
+				}
+			}
+		}
+		// clean boundary
+		for x in linspace(50., 60., 20) {
+			for y in linspace(0., 320., 20) {
+				for z in linspace(50., 60., 20) {
+					assert_eq!(
+						mapper
+							.global_to_flat(Location {
+								pos: vec3(x, y, z),
+								rot: Mat3::IDENTITY
+							})
+							.pos
+							.y,
+						y
+					);
+				}
+			}
+		}
+		for x in linspace(0., 60., 20) {
+			for y in linspace(300., 320., 20) {
+				for z in linspace(0., 60., 20) {
+					assert_eq!(
+						mapper
+							.global_to_flat(Location {
+								pos: vec3(x, y, z),
+								rot: Mat3::IDENTITY
+							})
+							.pos
+							.y,
+						y
+					);
+				}
+			}
+		}
+		// accelerating
+		for x in linspace(-29., 29., 20) {
+			for y in linspace(1., 299., 20) {
+				for z in linspace(-29., 29., 20) {
+					let v = mapper
+						.global_to_flat(Location {
+							pos: vec3(x, y, z),
+							rot: Mat3::IDENTITY,
+						})
+						.pos
+						.y;
+					assert!(v > 200.0);
+					assert!(v > y);
+				}
+			}
+		}
+	}
+}

src/tube/metric.rs

@@ -1,7 +1,8 @@
-use glam::{f32, vec2, Mat2, Vec2};
+use glam::{f32, vec3, Mat3, Vec3};
 
 use crate::fns::{self, Limiter};
-use crate::riemann::{Decomp2, Metric, Tens2};
+use crate::mathx::Decomp3;
+use crate::riemann::{Metric, Tens3};
 
 #[derive(Copy, Clone, Debug)]
 pub struct Tube {
@@ -40,20 +41,20 @@ impl Tube {
 }
 
 impl Metric for Tube {
-	fn sqrt_at(&self, pos: Vec2) -> Decomp2 {
+	fn sqrt_at(&self, pos: Vec3) -> Decomp3 {
 		let sx = self.fx().value(pos.x);
 		let sy = self.fy().du(pos.y);
 		let s = sx + sy - sx * sy;
 		assert!(sx.is_finite());
 		assert!(sy.is_finite());
 		assert!(sy > 0.0);
-		Decomp2 {
+		Decomp3 {
-			ortho: Mat2::IDENTITY,
+			ortho: Mat3::IDENTITY,
-			diag: vec2(1.0, s),
+			diag: vec3(1.0, s, 1.0),
 		}
 	}
 
-	fn part_derivs_at(&self, pos: Vec2) -> Tens2 {
+	fn part_derivs_at(&self, pos: Vec3) -> Tens3 {
 		let sx = self.fx().value(pos.x);
 		let sy = self.fy().du(pos.y);
 		let s = sx + sy - sx * sy;
@@ -62,8 +63,9 @@ impl Metric for Tube {
 		let ds2_dx = 2.0 * s * (1.0 - sy) * dsx_dx;
 		let ds2_dy = 2.0 * s * (1.0 - sx) * dsy_dy;
 		[
-			Mat2::from_cols_array(&[0.0, 0.0, 0.0, ds2_dx]),
+			Mat3::from_cols_array(&[0., 0., 0., 0., ds2_dx, 0., 0., 0., 0.]),
-			Mat2::from_cols_array(&[0.0, 0.0, 0.0, ds2_dy]),
+			Mat3::from_cols_array(&[0., 0., 0., 0., ds2_dy, 0., 0., 0., 0.]),
+			Mat3::from_cols_array(&[0., 0., 0., 0., 0., 0., 0., 0., 0.]),
 		]
 	}
 }
@@ -71,10 +73,11 @@ impl Metric for Tube {
 #[cfg(test)]
 mod test {
 	use approx::assert_abs_diff_eq;
-	use glam::{vec2, Vec2};
+	use glam::{vec3, Vec3};
 	use itertools_num::linspace;
 
-	use crate::riemann::{Decomp2, Metric};
+	use crate::mathx::Decomp3;
+	use crate::riemann::Metric;
 	use crate::tube::Space;
 	use crate::types::Ray;
 
@@ -84,7 +87,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: Vec3) -> Decomp3 {
 				self.0.sqrt_at(pos)
 			}
 		}
@@ -98,18 +101,25 @@ 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, 20)
 		{
 			for y in itertools_num::linspace(
 				-mul * testee.external_halflength,
 				mul * testee.external_halflength,
-				100,
+				20,
 			) {
-				let pos = vec2(x, y);
+				for z in itertools_num::linspace(
-				let computed = testee.part_derivs_at(pos);
+					-mul * testee.outer_radius,
-				let reference = approx.part_derivs_at(pos);
+					mul * testee.outer_radius,
-				let eq = (0..2).all(|coord| computed[coord].abs_diff_eq(reference[coord], epsilon));
+					20,
-				assert!(eq, "Bad derivative computation at {pos}:\n  explicit:  {computed:?}\n  numerical: {reference:?}\n");
+				) {
+					let pos = vec3(x, y, z);
+					let computed = testee.part_derivs_at(pos);
+					let reference = approx.part_derivs_at(pos);
+					let eq =
+						(0..2).all(|coord| computed[coord].abs_diff_eq(reference[coord], epsilon));
+					assert!(eq, "Bad derivative computation at {pos}:\n  explicit:  {computed:?}\n  numerical: {reference:?}\n");
+				}
 			}
 		}
 	}
@@ -130,9 +140,9 @@ mod test {
 		let steps = 1024;
 		for ax in [-30.0 + ε, -25.0, -3.0, 17.0, 30.0 - ε] {
 			for bx in [0.0, ε, 1.0, 7.0, 30.0 - ε] {
-				let a = vec2(ax, -(space.tube.external_halflength + off));
+				let a = vec3(ax, -(space.tube.external_halflength + off), 0.);
-				let b = vec2(bx, space.tube.external_halflength + off);
+				let b = vec3(bx, space.tube.external_halflength + off, 0.);
-				let Δ = vec2(bx - ax, 2.0 * (space.tube.internal_halflength + off));
+				let Δ = vec3(bx - ax, 2.0 * (space.tube.internal_halflength + off), 0.);
 				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);
@@ -168,9 +178,9 @@ mod test {
 				space.tube.inner_radius - ε,
 				20,
 			) {
-				let a = vec2(ax, -(space.tube.external_halflength + off));
+				let a = vec3(ax, -(space.tube.external_halflength + off), 0.);
-				let b = vec2(bx, space.tube.external_halflength + off);
+				let b = vec3(bx, space.tube.external_halflength + off, 0.);
-				let Δ = vec2(bx - ax, 2.0 * (space.tube.internal_halflength + off));
+				let Δ = vec3(bx - ax, 2.0 * (space.tube.internal_halflength + off), 0.);
 				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);
@@ -196,9 +206,9 @@ mod test {
 		let off = 10.0;
 		let steps = 10000;
 		for x in [space.tube.inner_radius - ε, space.tube.inner_radius + ε] {
-			let a = vec2(x, -(space.tube.external_halflength + off));
+			let a = vec3(x, -(space.tube.external_halflength + off), 0.);
-			let b = vec2(x, space.tube.external_halflength + off);
+			let b = vec3(x, space.tube.external_halflength + off, 0.);
-			let Δ = vec2(0.0, 2.0 * (space.tube.internal_halflength + off));
+			let Δ = vec3(0.0, 2.0 * (space.tube.internal_halflength + off), 0.);
 			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);
@@ -223,9 +233,9 @@ mod test {
 		let off = 10.0;
 		let steps = 4096;
 		for x in [space.tube.outer_radius + ε, space.tube.outer_radius - ε] {
-			let a = vec2(x, -(space.tube.external_halflength + off));
+			let a = vec3(x, -(space.tube.external_halflength + off), 0.);
-			let b = vec2(x, space.tube.external_halflength + off);
+			let b = vec3(x, space.tube.external_halflength + off, 0.);
-			let Δ = vec2(0.0, 2.0 * (space.tube.external_halflength + off));
+			let Δ = vec3(0.0, 2.0 * (space.tube.external_halflength + off), 0.);
 			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);

src/tube/mod.rs

@@ -1,4 +1,4 @@
-use glam::{bool, f32, vec2, Mat2, Vec2};
+use glam::{bool, f32, vec3, Mat3, Vec3};
 
 use crate::ifaces::{DebugTraceable, RayPath, Traceable};
 use coords::{FlatCoordinateSystem, InnerCS, OuterCS};
@@ -26,7 +26,7 @@ pub enum Subspace {
 }
 
 impl Space {
-	fn which_subspace(&self, pt: Vec2) -> Subspace {
+	fn which_subspace(&self, pt: Vec3) -> Subspace {
 		if pt.y.abs() > self.tube.external_halflength {
 			Outer
 		} else if pt.x.abs() > self.tube.outer_radius {
@@ -41,7 +41,7 @@ impl Space {
 	/// Выполняет один шаг трассировки. Работает в любой части пространства, но вне Boundary доступны более эффективные методы.
 	/// ray задаётся в основной СК.
 	pub fn trace_step(&self, ray: Ray) -> Ray {
-		let a: Vec2 = -riemann::contract2(riemann::krist(&self.tube, ray.pos), ray.dir);
+		let a = -riemann::contract2(riemann::krist(&self.tube, ray.pos), ray.dir);
 		let v = ray.dir + a;
 		let p = ray.pos + v;
 		Ray { pos: p, dir: v }
@@ -49,9 +49,9 @@ impl Space {
 
 	/// Выполняет один шаг перемещения. Работает в любой части пространства.
 	/// off задаётся в локальной СК. Рекомендуется считать небольшими шагами.
-	pub fn move_step(&self, loc: Location, off: Vec2) -> Location {
+	pub fn move_step(&self, loc: Location, off: Vec3) -> Location {
 		let corr =
-			Mat2::IDENTITY - riemann::contract(riemann::krist(&self.tube, loc.pos), loc.rot * off);
+			Mat3::IDENTITY - riemann::contract(riemann::krist(&self.tube, loc.pos), loc.rot * off);
 		let p = loc.pos + corr * loc.rot * off;
 		Location {
 			pos: p,
@@ -73,7 +73,7 @@ impl Space {
 		self.trace_flat(OuterCS(self.tube), ray)
 	}
 
-	fn obj_hitter(&self, pos: Vec2) -> Option<fn(&Self, ray: Ray) -> FlatTraceResult> {
+	fn obj_hitter(&self, pos: Vec3) -> Option<fn(&Self, ray: Ray) -> FlatTraceResult> {
 		match self.which_subspace(pos) {
 			Inner => Some(Self::trace_inner),
 			Outer => Some(Self::trace_outer),
@@ -113,7 +113,7 @@ impl Space {
 		objs: &[Object],
 		ray: Ray,
 		limit: Option<f32>,
-		globalize: impl Fn(Vec2) -> Vec2,
+		globalize: impl Fn(Vec3) -> Vec3,
 	) -> Vec<Hit> {
 		let limit = limit.unwrap_or(f32::INFINITY);
 		objs.iter()
@@ -146,16 +146,25 @@ impl Space {
 			.collect()
 	}
 
-	pub fn line(&self, a: Vec2, b: Vec2, step: f32) -> Vec<Vec2> {
+	pub fn line(&self, a: Vec3, b: Vec3, step: f32) -> Vec<Ray> {
 		match self.which_subspace(a) {
-			Outer => vec![b],
+			Outer => vec![Ray {
+				pos: b,
+				dir: (b - a).normalize(),
+			}],
 			Inner => {
 				let cs = InnerCS(self.tube);
 				let n = ((b - a).length() / step) as usize + 1;
 				let a = cs.global_to_flat(a);
 				let b = cs.global_to_flat(b);
+				let dir = (b - a).normalize();
 				(1..=n)
-					.map(|k| cs.flat_to_global(a.lerp(b, k as f32 / n as f32)))
+					.map(|k| {
+						cs.flat_to_global(Ray {
+							pos: a.lerp(b, k as f32 / n as f32),
+							dir,
+						})
+					})
 					.collect()
 			}
 			Boundary => panic!("Can't draw a line here!"),
@@ -210,9 +219,9 @@ impl DebugTraceable for Space {
 		let mut hits = vec![];
 		let mut ray = self.camera_ray_to_abs(camera, ray);
 
-		let trace_to_flat = |points: &mut Vec<Vec2>, ray| {
+		let trace_to_flat = |points: &mut Vec<Ray>, ray| {
 			for ray in self.trace_iter(ray).skip(1) {
-				points.push(ray.pos);
+				points.push(ray);
 				if let Some(hitter) = self.obj_hitter(ray.pos) {
 					return (ray, hitter(self, ray));
 				}
@@ -220,18 +229,12 @@ impl DebugTraceable for Space {
 			unreachable!("Space::trace_iter terminated!")
 		};
 
-		points.push(ray.pos);
+		points.push(ray);
 		for _ in 0..100 {
 			let (ray_into_flat, ret) = trace_to_flat(&mut points, ray);
 			hits.extend(ret.objects); // TODO fix distance
 			let Some(ray_outta_flat) = ret.end else {
-				return (
+				return (hits, RayPath { points });
-					hits,
-					RayPath {
-						points,
-						end_dir: ray_into_flat.dir.normalize(),
-					},
-				);
 			};
 			points.extend(self.line(ray_into_flat.pos, ray_outta_flat.pos, 10.0));
 			ray = ray_outta_flat;
@@ -241,7 +244,7 @@ impl DebugTraceable for Space {
 }
 
 struct Rect {
-	pub size: Vec2,
+	pub size: Vec3,
 }
 
 impl Rect {
@@ -253,7 +256,7 @@ impl Rect {
 		}
 	}
 
-	fn is_inside(&self, pt: Vec2) -> bool {
+	fn is_inside(&self, pt: Vec3) -> bool {
 		pt.abs().cmplt(self.size).all()
 	}
 
@@ -285,146 +288,146 @@ impl Rect {
 fn test_rect() {
 	assert_eq!(
 		Rect::flip_ray(Ray {
-			pos: vec2(2.0, 3.0),
+			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec2(4.0, 5.0)
+			dir: vec3(4.0, 5.0, 4.0)
 		}),
 		Ray {
-			pos: vec2(2.0, 3.0),
+			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec2(4.0, 5.0)
+			dir: vec3(4.0, 5.0, 4.0)
 		}
 	);
 	assert_eq!(
 		Rect::flip_ray(Ray {
-			pos: vec2(2.0, 3.0),
+			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec2(-4.0, 5.0)
+			dir: vec3(-4.0, 5.0, -4.0)
 		}),
 		Ray {
-			pos: vec2(-2.0, 3.0),
+			pos: vec3(-2.0, 3.0, -2.0),
-			dir: vec2(4.0, 5.0)
+			dir: vec3(4.0, 5.0, 4.0)
 		}
 	);
 	assert_eq!(
 		Rect::flip_ray(Ray {
-			pos: vec2(2.0, 3.0),
+			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec2(4.0, -5.0)
+			dir: vec3(4.0, -5.0, 4.0)
 		}),
 		Ray {
-			pos: vec2(2.0, -3.0),
+			pos: vec3(2.0, -3.0, 2.0),
-			dir: vec2(4.0, 5.0)
+			dir: vec3(4.0, 5.0, 4.0)
 		}
 	);
 	assert_eq!(
 		Rect::flip_ray(Ray {
-			pos: vec2(2.0, 3.0),
+			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec2(4.0, 0.0)
+			dir: vec3(4.0, 0.0, 4.0)
 		}),
 		Ray {
-			pos: vec2(2.0, 3.0),
+			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec2(4.0, 0.0)
+			dir: vec3(4.0, 0.0, 4.0)
 		}
 	);
 
 	let r = Rect {
-		size: vec2(2.0, 3.0),
+		size: vec3(2.0, 3.0, 2.0),
 	};
 
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(3.0, 3.0),
+			pos: vec3(3.0, 3.0, 3.0),
-			dir: vec2(1.0, 1.0)
+			dir: vec3(1.0, 1.0, 1.0)
 		}),
 		None
 	);
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(-3.0, 2.0),
+			pos: vec3(-3.0, 2.0, -3.0),
-			dir: vec2(1.0, 0.0)
+			dir: vec3(1.0, 0.0, 1.0)
 		}),
 		Some(1.0)
 	);
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(-3.0, 2.0),
+			pos: vec3(-3.0, 2.0, -3.0),
-			dir: vec2(-1.0, 0.0)
+			dir: vec3(-1.0, 0.0, -1.0)
 		}),
 		None
 	);
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(-3.0, 1.0),
+			pos: vec3(-3.0, 1.0, -3.0),
-			dir: vec2(2.0, 2.0)
+			dir: vec3(2.0, 2.0, 2.0)
 		}),
 		Some(0.5)
 	);
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(-3.0, 2.1),
+			pos: vec3(-3.0, 2.1, -3.0),
-			dir: vec2(2.0, 2.0)
+			dir: vec3(2.0, 2.0, 2.0)
 		}),
 		None
 	);
 
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(2.0, 3.0),
+			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec2(1.0, 1.0)
+			dir: vec3(1.0, 1.0, 1.0)
 		}),
 		None
 	);
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(-2.0, 3.0),
+			pos: vec3(-2.0, 3.0, -2.0),
-			dir: vec2(-1.0, 1.0)
+			dir: vec3(-1.0, 1.0, -1.0)
 		}),
 		None
 	);
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(2.0, 3.0),
+			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec2(-1.0, -1.0)
+			dir: vec3(-1.0, -1.0, -1.0)
 		}),
 		Some(0.0)
 	);
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(2.0, -3.0),
+			pos: vec3(2.0, -3.0, 2.0),
-			dir: vec2(-1.0, 1.0)
+			dir: vec3(-1.0, 1.0, -1.0)
 		}),
 		Some(0.0)
 	);
 
 	assert_eq!(
 		r.trace_out_of(Ray {
-			pos: vec2(0.0, 0.0),
+			pos: vec3(0.0, 0.0, 0.0),
-			dir: vec2(1.0, 1.0)
+			dir: vec3(1.0, 1.0, 1.0)
 		}),
 		Some(2.0)
 	);
 	assert_eq!(
 		r.trace_out_of(Ray {
-			pos: vec2(0.0, 0.0),
+			pos: vec3(0.0, 0.0, 0.0),
-			dir: vec2(0.0, 1.0)
+			dir: vec3(0.0, 1.0, 0.0)
 		}),
 		Some(3.0)
 	);
 	assert_eq!(
 		r.trace_out_of(Ray {
-			pos: vec2(0.0, 1.0),
+			pos: vec3(0.0, 1.0, 0.0),
-			dir: vec2(0.0, -1.0)
+			dir: vec3(0.0, -1.0, 0.0)
 		}),
 		Some(4.0)
 	);
 	assert_eq!(
 		r.trace_out_of(Ray {
-			pos: vec2(1.0, 1.0),
+			pos: vec3(1.0, 1.0, 1.0),
-			dir: vec2(0.0, -1.0)
+			dir: vec3(0.0, -1.0, 0.0)
 		}),
 		Some(4.0)
 	);
 	assert_eq!(
 		r.trace_out_of(Ray {
-			pos: vec2(2.0, 3.0),
+			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec2(1.0, 1.0)
+			dir: vec3(1.0, 1.0, 1.0)
 		}),
 		Some(0.0)
 	);

src/types.rs

@@ -1,9 +1,9 @@
-use glam::{f32, i32, Mat2, Vec2};
+use glam::{f32, i32, Mat3, Vec3};
 
 #[derive(Copy, Clone, Debug, PartialEq)]
 pub struct Ray {
-	pub pos: Vec2,
+	pub pos: Vec3,
-	pub dir: Vec2,
+	pub dir: Vec3,
 }
 
 impl Ray {
@@ -15,7 +15,7 @@ impl Ray {
 	}
 }
 
-impl std::ops::Mul<Ray> for Mat2 {
+impl std::ops::Mul<Ray> for Mat3 {
 	type Output = Ray;
 
 	fn mul(self, rhs: Ray) -> Self::Output {
@@ -29,9 +29,9 @@ impl std::ops::Mul<Ray> for Mat2 {
 #[derive(Copy, Clone, Debug, PartialEq)]
 pub struct Location {
 	/// Положение в основной СК
-	pub pos: Vec2,
+	pub pos: Vec3,
 	/// Преобразование вектора из локальной ортонормированной в основную СК
-	pub rot: Mat2,
+	pub rot: Mat3,
 }
 
 #[derive(Copy, Clone, Debug)]
@@ -45,7 +45,7 @@ pub struct Object {
 pub struct Hit {
 	pub distance: f32,
 	pub id: i32,
-	pub pos: Vec2, // положение в основной СК
+	pub pos: Vec3, // положение в основной СК
 	pub rel: Ray,  // в локальной ортонормированной СК объекта
 }
 

src/utils.rs

@@ -0,0 +1,92 @@
+//! Utility functions to work with metrics.
+
+use crate::{
+	mathx::MatExt as _,
+	riemann::{trace_iter, Metric},
+	types::Location,
+};
+use glam::{Mat3, Vec3};
+
+pub fn rel_to_abs(space: &impl Metric, base: &Location, rel: Vec3, steps: usize) -> Vec3 {
+	let c = 1.0 / (steps as f32);
+	trace_iter(space, base.pos, base.rot * rel, c * rel.length())
+		.nth(steps - 1)
+		.unwrap()
+}
+
+/// Converts a position and a rotation to a [Location]. Only the X direction is preserved from `rot` to ensure the resulting Location describes an orthonormal coordinate system.
+pub fn put_object(space: &impl Metric, pos: Vec3, rot: Mat3) -> Location {
+	let metric_sqrt = space.sqrt_at(pos);
+	let metric_inv_sqrt = space.sqrt_at(pos).inverse();
+	let rot = metric_inv_sqrt * (metric_sqrt * rot).orthonormalize();
+	Location { pos, rot }
+}
+
+#[cfg(test)]
+mod tests {
+	use super::*;
+	use crate::riemann::samples;
+	use glam::{mat3, vec3};
+
+	#[test]
+	fn test_put_object() {
+		use approx::assert_abs_diff_eq;
+
+		let ε = 1e-5;
+		let m = samples::ScaledMetric {
+			scale: vec3(3., 4., 5.),
+		};
+
+		let loc = put_object(
+			&m,
+			vec3(1., 2., 0.),
+			mat3(vec3(1., 0., 0.), vec3(0., 1., 0.), vec3(0., 0., 1.)),
+		);
+		assert_eq!(loc.pos, vec3(1., 2., 0.));
+		assert_abs_diff_eq!(
+			loc.rot * vec3(1., 0., 0.),
+			vec3(1. / 3., 0., 0.),
+			epsilon = ε
+		);
+		assert_abs_diff_eq!(
+			loc.rot * vec3(0., 1., 0.),
+			vec3(0., 1. / 4., 0.),
+			epsilon = ε
+		);
+
+		let loc = put_object(
+			&m,
+			vec3(1., 2., 0.),
+			mat3(vec3(0., 1., 0.), vec3(-1., 0., 0.), vec3(0., 0., 1.)),
+		);
+		assert_eq!(loc.pos, vec3(1., 2., 0.));
+		assert_abs_diff_eq!(
+			loc.rot * vec3(1., 0., 0.),
+			vec3(0., 1. / 4., 0.),
+			epsilon = ε
+		);
+		assert_abs_diff_eq!(
+			loc.rot * vec3(0., 1., 0.),
+			vec3(-1. / 3., 0., 0.),
+			epsilon = ε
+		);
+
+		let c = 0.5 * std::f32::consts::SQRT_2;
+		let loc = put_object(
+			&m,
+			vec3(1., 2., 0.),
+			mat3(vec3(c, c, 0.), vec3(-c, c, 0.), vec3(0., 0., 1.)),
+		);
+		assert_eq!(loc.pos, vec3(1., 2., 0.));
+		assert_abs_diff_eq!(
+			loc.rot * vec3(1., 0., 0.),
+			vec3(1. / 5., 1. / 5., 0.),
+			epsilon = ε
+		);
+		assert_abs_diff_eq!(
+			loc.rot * vec3(0., 1., 0.),
+			vec3(-4. / 15., 3. / 20., 0.),
+			epsilon = ε
+		);
+	}
+}