Pre-generate normals, for smooth look

Anti-aliasing works better this way

.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

@@ -3,7 +3,8 @@ name = "refraction"
 version = "0.1.0"
 edition = "2021"
 
-# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
+[lints.rust]
+mixed_script_confusables = "allow"
 
 [profile.dev]
 panic = 'abort'
@@ -11,11 +12,23 @@ panic = 'abort'
 [profile.dev.package."*"]
 opt-level = 3
 
+[profile.test.package."*"]
+opt-level = 3
+
 [dependencies]
-glam = "0.27.0"
+rand = "0.8.5"
+glam = { version = "0.27.0", features = ["approx", "fast-math", "rand"] }
+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"
+rand = "0.8.5"
+rand_pcg = "0.3.1"

models/spacecraft2.1.obj

@@ -0,0 +1,140 @@
+# Blender 3.6.5
+# www.blender.org
+mtllib spacecraft2.1.mtl
+o Cube
+v -1.000000 0.000000 5.000000
+v -3.000000 0.000000 4.000000
+v -4.000000 0.000000 7.000000
+v -8.000000 0.000000 7.000000
+v -9.000000 0.000000 4.000000
+v -3.000000 0.000000 0.000000
+v -6.000000 0.000000 -7.000000
+v -2.000000 0.000000 -3.000000
+v -7.000000 1.000000 4.000000
+v 0.000000 2.000000 0.000000
+v 0.000000 2.000000 2.000000
+v 0.000000 0.000000 6.000000
+v 0.000000 0.000000 -4.000000
+v -5.000000 0.000000 2.000000
+v 1.000000 0.000000 5.000000
+v 3.000000 0.000000 4.000000
+v 4.000000 0.000000 7.000000
+v 8.000000 0.000000 7.000000
+v 9.000000 0.000000 4.000000
+v 3.000000 0.000000 0.000000
+v 6.000000 0.000000 -7.000000
+v 2.000000 0.000000 -3.000000
+v 7.000000 1.000000 4.000000
+v 5.000000 0.000000 2.000000
+v -7.000000 -1.000000 4.000000
+v 0.000000 -2.000000 0.000000
+v 0.000000 -2.000000 2.000000
+v 7.000000 -1.000000 4.000000
+vn -0.5455 0.8182 -0.1818
+vn -0.4851 0.4851 -0.7276
+vn -0.4439 0.8878 -0.1211
+vn -0.0000 0.9487 0.3162
+vn -0.4423 0.8847 0.1474
+vn 0.2417 0.9670 0.0806
+vn -0.2408 0.8427 0.4815
+vn 0.1455 0.5819 0.8001
+vn -0.1414 0.9899 -0.0000
+vn -0.2182 0.8729 -0.4364
+vn -0.4082 0.8165 0.4082
+vn -0.4851 0.7276 -0.4851
+vn 0.4099 0.9110 -0.0455
+vn -0.0000 -1.0000 -0.0000
+vn 0.5455 0.8182 -0.1818
+vn 0.4851 0.4851 -0.7276
+vn 0.4439 0.8878 -0.1211
+vn 0.4423 0.8847 0.1474
+vn -0.2417 0.9670 0.0806
+vn 0.2408 0.8427 0.4815
+vn -0.1455 0.5819 0.8001
+vn 0.1414 0.9899 -0.0000
+vn 0.2182 0.8729 -0.4364
+vn 0.4082 0.8165 0.4082
+vn 0.4851 0.7276 -0.4851
+vn -0.4099 0.9110 -0.0455
+vn -0.5455 -0.8182 -0.1818
+vn -0.4851 -0.4851 -0.7276
+vn -0.4439 -0.8878 -0.1211
+vn -0.0000 -0.9487 0.3162
+vn -0.4423 -0.8847 0.1474
+vn 0.2417 -0.9670 0.0806
+vn -0.2408 -0.8427 0.4815
+vn 0.1455 -0.5819 0.8001
+vn -0.1414 -0.9899 -0.0000
+vn -0.2182 -0.8729 -0.4364
+vn -0.4082 -0.8165 0.4082
+vn -0.4851 -0.7276 -0.4851
+vn 0.4099 -0.9110 -0.0455
+vn 0.5455 -0.8182 -0.1818
+vn 0.4851 -0.4851 -0.7276
+vn 0.4439 -0.8878 -0.1211
+vn 0.4423 -0.8847 0.1474
+vn -0.2417 -0.9670 0.0806
+vn 0.2408 -0.8427 0.4815
+vn -0.1455 -0.5819 0.8001
+vn 0.1414 -0.9899 -0.0000
+vn 0.2182 -0.8729 -0.4364
+vn 0.4082 -0.8165 0.4082
+vn 0.4851 -0.7276 -0.4851
+vn -0.4099 -0.9110 -0.0455
+vt 0.000000 0.000000
+s 0
+usemtl Material
+f 6/1/1 10/1/1 8/1/1
+f 10/1/2 14/1/2 9/1/2
+f 5/1/3 9/1/3 7/1/3
+f 3/1/4 9/1/4 4/1/4
+f 4/1/5 9/1/5 5/1/5
+f 2/1/6 9/1/6 3/1/6
+f 1/1/7 11/1/7 2/1/7
+f 9/1/8 2/1/8 11/1/8
+f 9/1/9 11/1/9 10/1/9
+f 8/1/10 10/1/10 13/1/10
+f 1/1/11 12/1/11 11/1/11
+f 10/1/12 6/1/12 14/1/12
+f 7/1/13 9/1/13 14/1/13
+f 7/1/14 14/1/14 5/1/14
+f 20/1/15 22/1/15 10/1/15
+f 10/1/16 23/1/16 24/1/16
+f 19/1/17 21/1/17 23/1/17
+f 17/1/4 18/1/4 23/1/4
+f 18/1/18 19/1/18 23/1/18
+f 16/1/19 17/1/19 23/1/19
+f 15/1/20 16/1/20 11/1/20
+f 23/1/21 11/1/21 16/1/21
+f 23/1/22 10/1/22 11/1/22
+f 22/1/23 13/1/23 10/1/23
+f 15/1/24 11/1/24 12/1/24
+f 10/1/25 24/1/25 20/1/25
+f 21/1/26 24/1/26 23/1/26
+f 21/1/14 19/1/14 24/1/14
+f 6/1/27 8/1/27 26/1/27
+f 26/1/28 25/1/28 14/1/28
+f 5/1/29 7/1/29 25/1/29
+f 3/1/30 4/1/30 25/1/30
+f 4/1/31 5/1/31 25/1/31
+f 2/1/32 3/1/32 25/1/32
+f 1/1/33 2/1/33 27/1/33
+f 25/1/34 27/1/34 2/1/34
+f 25/1/35 26/1/35 27/1/35
+f 8/1/36 13/1/36 26/1/36
+f 1/1/37 27/1/37 12/1/37
+f 26/1/38 14/1/38 6/1/38
+f 7/1/39 14/1/39 25/1/39
+f 20/1/40 26/1/40 22/1/40
+f 26/1/41 24/1/41 28/1/41
+f 19/1/42 28/1/42 21/1/42
+f 17/1/30 28/1/30 18/1/30
+f 18/1/43 28/1/43 19/1/43
+f 16/1/44 28/1/44 17/1/44
+f 15/1/45 27/1/45 16/1/45
+f 28/1/46 16/1/46 27/1/46
+f 28/1/47 27/1/47 26/1/47
+f 22/1/48 26/1/48 13/1/48
+f 15/1/49 12/1/49 27/1/49
+f 26/1/50 20/1/50 24/1/50
+f 21/1/51 28/1/51 24/1/51

models/spacecraft2.2.obj

@@ -0,0 +1,77 @@
+# Blender 3.6.5
+# www.blender.org
+mtllib spacecraft2.2.mtl
+o Cube
+v -1.000000 0.000000 5.000000
+v -3.000000 0.000000 4.000000
+v -3.000000 0.000000 0.000000
+v -2.000000 0.000000 -3.000000
+v -7.000000 0.000000 4.000000
+v 0.000000 2.000000 0.000000
+v 0.000000 2.000000 2.000000
+v 0.000000 0.000000 6.000000
+v 0.000000 0.000000 -4.000000
+v 1.000000 0.000000 5.000000
+v 3.000000 0.000000 4.000000
+v 3.000000 0.000000 0.000000
+v 2.000000 0.000000 -3.000000
+v 7.000000 0.000000 4.000000
+v 0.000000 -2.000000 0.000000
+v 0.000000 -2.000000 2.000000
+vn -0.5455 0.8182 -0.1818
+vn -0.4851 0.7276 -0.4851
+vn -0.2408 0.8427 0.4815
+vn -0.0000 0.7071 0.7071
+vn -0.2747 0.9615 -0.0000
+vn -0.2182 0.8729 -0.4364
+vn -0.4082 0.8165 0.4082
+vn 0.5455 0.8182 -0.1818
+vn 0.4851 0.7276 -0.4851
+vn 0.2408 0.8427 0.4815
+vn 0.2747 0.9615 -0.0000
+vn 0.2182 0.8729 -0.4364
+vn 0.4082 0.8165 0.4082
+vn -0.5455 -0.8182 -0.1818
+vn -0.4851 -0.7276 -0.4851
+vn -0.2408 -0.8427 0.4815
+vn -0.0000 -0.7071 0.7071
+vn -0.2747 -0.9615 -0.0000
+vn -0.2182 -0.8729 -0.4364
+vn -0.4082 -0.8165 0.4082
+vn 0.5455 -0.8182 -0.1818
+vn 0.4851 -0.7276 -0.4851
+vn 0.2408 -0.8427 0.4815
+vn 0.2747 -0.9615 -0.0000
+vn 0.2182 -0.8729 -0.4364
+vn 0.4082 -0.8165 0.4082
+vt 0.000000 0.000000
+s 0
+usemtl Material
+f 3/1/1 6/1/1 4/1/1
+f 6/1/2 3/1/2 5/1/2
+f 1/1/3 7/1/3 2/1/3
+f 5/1/4 2/1/4 7/1/4
+f 5/1/5 7/1/5 6/1/5
+f 4/1/6 6/1/6 9/1/6
+f 1/1/7 8/1/7 7/1/7
+f 12/1/8 13/1/8 6/1/8
+f 6/1/9 14/1/9 12/1/9
+f 10/1/10 11/1/10 7/1/10
+f 14/1/4 7/1/4 11/1/4
+f 14/1/11 6/1/11 7/1/11
+f 13/1/12 9/1/12 6/1/12
+f 10/1/13 7/1/13 8/1/13
+f 3/1/14 4/1/14 15/1/14
+f 15/1/15 5/1/15 3/1/15
+f 1/1/16 2/1/16 16/1/16
+f 5/1/17 16/1/17 2/1/17
+f 5/1/18 15/1/18 16/1/18
+f 4/1/19 9/1/19 15/1/19
+f 1/1/20 16/1/20 8/1/20
+f 12/1/21 15/1/21 13/1/21
+f 15/1/22 12/1/22 14/1/22
+f 10/1/23 16/1/23 11/1/23
+f 14/1/17 11/1/17 16/1/17
+f 14/1/24 16/1/24 15/1/24
+f 13/1/25 15/1/25 9/1/25
+f 10/1/26 8/1/26 16/1/26

models/spacecraft2.obj

@@ -0,0 +1,109 @@
+# Blender 3.6.5
+# www.blender.org
+mtllib spacecraft2.mtl
+o Cube
+v -1.000000 0.000000 3.000000
+v -2.000000 0.000000 2.000000
+v -3.000000 0.000000 4.000000
+v -5.000000 0.000000 4.000000
+v -6.000000 0.000000 2.000000
+v -2.000000 0.000000 0.000000
+v -4.000000 0.000000 -6.000000
+v -1.000000 0.000000 -2.000000
+v 1.000000 0.000000 3.000000
+v 2.000000 0.000000 2.000000
+v 3.000000 0.000000 4.000000
+v 5.000000 0.000000 4.000000
+v 6.000000 0.000000 2.000000
+v 4.000000 0.000000 -6.000000
+v 2.000000 0.000000 0.000000
+v 1.000000 0.000000 -2.000000
+v -4.000000 1.000000 2.000000
+v 0.000000 1.000000 0.000000
+v 4.000000 1.000000 2.000000
+v 0.000000 1.000000 1.000000
+v -4.000000 -1.000000 2.000000
+v 0.000000 -1.000000 0.000000
+v 4.000000 -1.000000 2.000000
+v 0.000000 -1.000000 1.000000
+vn -0.0000 0.8944 -0.4472
+vn 0.4364 0.8729 -0.2182
+vn 0.3333 0.6667 -0.6667
+vn -0.3487 0.9300 -0.1162
+vn 0.4444 0.8889 -0.1111
+vn 0.4364 0.8729 0.2182
+vn -0.0000 0.8944 0.4472
+vn -0.4364 0.8729 0.2182
+vn -0.0000 1.0000 -0.0000
+vn 0.3015 0.9045 0.3015
+vn -0.2182 0.4364 0.8729
+vn -0.4364 0.8729 -0.2182
+vn -0.3333 0.6667 -0.6667
+vn 0.3487 0.9300 -0.1162
+vn -0.4444 0.8889 -0.1111
+vn -0.3015 0.9045 0.3015
+vn 0.2182 0.4364 0.8729
+vn -0.0000 -0.8944 -0.4472
+vn 0.4364 -0.8729 -0.2182
+vn 0.3333 -0.6667 -0.6667
+vn -0.3487 -0.9300 -0.1162
+vn 0.4444 -0.8889 -0.1111
+vn 0.4364 -0.8729 0.2182
+vn -0.0000 -0.8944 0.4472
+vn -0.4364 -0.8729 0.2182
+vn -0.0000 -1.0000 -0.0000
+vn 0.3015 -0.9045 0.3015
+vn -0.2182 -0.4364 0.8729
+vn -0.4364 -0.8729 -0.2182
+vn -0.3333 -0.6667 -0.6667
+vn 0.3487 -0.9300 -0.1162
+vn -0.4444 -0.8889 -0.1111
+vn -0.3015 -0.9045 0.3015
+vn 0.2182 -0.4364 0.8729
+vt 0.000000 0.000000
+s 0
+usemtl Material
+f 8/1/1 18/1/1 16/1/1
+f 16/1/2 18/1/2 15/1/2
+f 18/1/3 19/1/3 15/1/3
+f 15/1/4 19/1/4 14/1/4
+f 14/1/5 19/1/5 13/1/5
+f 13/1/6 19/1/6 12/1/6
+f 12/1/7 19/1/7 11/1/7
+f 10/1/8 11/1/8 19/1/8
+f 20/1/9 18/1/9 17/1/9
+f 1/1/7 9/1/7 20/1/7
+f 9/1/10 10/1/10 20/1/10
+f 20/1/11 10/1/11 19/1/11
+f 6/1/12 18/1/12 8/1/12
+f 18/1/13 6/1/13 17/1/13
+f 6/1/14 7/1/14 17/1/14
+f 5/1/15 17/1/15 7/1/15
+f 3/1/7 17/1/7 4/1/7
+f 4/1/8 17/1/8 5/1/8
+f 2/1/6 17/1/6 3/1/6
+f 1/1/16 20/1/16 2/1/16
+f 17/1/17 2/1/17 20/1/17
+f 20/1/9 19/1/9 18/1/9
+f 8/1/18 16/1/18 22/1/18
+f 16/1/19 15/1/19 22/1/19
+f 22/1/20 15/1/20 23/1/20
+f 15/1/21 14/1/21 23/1/21
+f 14/1/22 13/1/22 23/1/22
+f 13/1/23 12/1/23 23/1/23
+f 12/1/24 11/1/24 23/1/24
+f 10/1/25 23/1/25 11/1/25
+f 24/1/26 21/1/26 22/1/26
+f 1/1/24 24/1/24 9/1/24
+f 9/1/27 24/1/27 10/1/27
+f 24/1/28 23/1/28 10/1/28
+f 6/1/29 8/1/29 22/1/29
+f 22/1/30 21/1/30 6/1/30
+f 6/1/31 21/1/31 7/1/31
+f 5/1/32 7/1/32 21/1/32
+f 3/1/24 4/1/24 21/1/24
+f 4/1/25 5/1/25 21/1/25
+f 2/1/23 3/1/23 21/1/23
+f 1/1/33 2/1/33 24/1/33
+f 21/1/34 24/1/34 2/1/34
+f 24/1/26 22/1/26 23/1/26

rustfmt.toml

@@ -0,0 +1,2 @@
+hard_tabs = true
+max_width = 120

src/bin/flat/float_fun.rs

@@ -1,36 +0,0 @@
-use glam::FloatExt;
-
-mod bounds {
-	pub trait Pair<T> {}
-
-	impl<T> Pair<T> for (T, T) {}
-}
-
-pub trait FloatExt2<T>: bounds::Pair<T> {
-	fn lerp(self, t: T) -> T;
-	fn inverse_lerp(self, y: T) -> T;
-}
-
-impl<F: FloatExt> FloatExt2<F> for (F, F) {
-	fn lerp(self, t: F) -> F { F::lerp(self.0, self.1, t) }
-	fn inverse_lerp(self, y: F) -> F { F::inverse_lerp(self.0, self.1, y) }
-}
-
-#[cfg(test)]
-mod test {
-	use super::FloatExt2;
-
-	#[test]
-	fn test_lerp() {
-		assert_eq!((3., 7.).lerp(-0.5), 1.);
-		assert_eq!((3., 7.).lerp(0.0), 3.);
-		assert_eq!((3., 7.).lerp(0.5), 5.);
-		assert_eq!((3., 7.).lerp(1.0), 7.);
-		assert_eq!((3., 7.).lerp(1.5), 9.);
-		assert_eq!((3., 7.).inverse_lerp(1.), -0.5);
-		assert_eq!((3., 7.).inverse_lerp(3.), 0.0);
-		assert_eq!((3., 7.).inverse_lerp(5.), 0.5);
-		assert_eq!((3., 7.).inverse_lerp(7.), 1.0);
-		assert_eq!((3., 7.).inverse_lerp(9.), 1.5);
-	}
-}

src/bin/flat/fns.rs

@@ -1,110 +0,0 @@
-use crate::float_fun::FloatExt2;
-
-pub struct LinearLimiter {
-	pub min: f32,
-	pub max: f32,
-}
-
-impl LinearLimiter {
-	pub fn value(&self, x: f32) -> f32 { (self.min, self.max).inverse_lerp(x.abs()).clamp(0.0, 1.0) }
-	pub fn derivative(&self, x: f32) -> f32 { if x.abs() > self.min && x.abs() < self.max { x.signum() / (self.max - self.min) } else { 0.0 } }
-}
-
-pub struct QuadraticAccelerator {
-	pub internal: f32,
-	pub external: f32,
-}
-
-/// Продолжает функцию f с [-lim, lim] линейно в предположении f(±lim) = ±val, f'(±lim) = 1.
-fn extend_linear(t: f32, f: impl FnOnce(f32) -> f32, lim: f32, val: f32) -> f32 {
-	if t.abs() <= lim { f(t) } else { t + t.signum() * (val - lim) }
-}
-
-/// Продолжает функцию f с [-lim, lim] константой в предположении f(±lim) = val, f'(±lim) = 0.
-fn extend_const(t: f32, f: impl FnOnce(f32) -> f32, lim: f32, val: f32) -> f32 {
-	if t.abs() <= lim { f(t) } else { val }
-}
-
-impl QuadraticAccelerator {
-	fn a(&self) -> f32 { -(self.external - self.internal) / self.internal.powi(2) }
-	fn b(&self) -> f32 { 2.0 * self.external / self.internal - 1.0 }
-	fn root(&self, x: f32) -> f32 { (self.b().powi(2) + 4.0 * self.a() * x.abs()).sqrt() }
-
-	pub fn x(&self, u: f32) -> f32 { extend_linear(u, |u| (self.a() * u.abs() + self.b()) * u, self.internal, self.external) }
-	pub fn u(&self, x: f32) -> f32 { extend_linear(x, |x| 0.5 * x.signum() * (-self.b() + self.root(x)) / self.a(), self.external, self.internal) }
-	pub fn dx(&self, u: f32) -> f32 { extend_const(u, |u| 2.0 * self.a() * u.abs() + self.b(), self.internal, 1.0) }
-	pub fn du(&self, x: f32) -> f32 { extend_const(x, |x| 1.0 / self.root(x), self.external, 1.0) }
-	pub fn d2u(&self, x: f32) -> f32 { extend_const(x, |x| -2.0 * x.signum() * self.a() * self.root(x).powi(-3), self.external, 0.0) }
-}
-
-#[cfg(test)]
-mod test {
-	use approx::{abs_diff_eq, AbsDiffEq, assert_abs_diff_eq};
-
-	#[test]
-	fn test_linear_limiter() {
-		let testee = super::LinearLimiter { min: 20.0, max: 30.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;
-		let mul = 1.0 + margin;
-		for x in itertools_num::linspace(0., testee.min, 10) {
-			assert_abs_diff_eq!(testee.value(x), 0., epsilon = ε);
-			assert_abs_diff_eq!(testee.value(-x), 0., epsilon = ε);
-		}
-		for x in itertools_num::linspace(testee.max, mul * testee.max, 10) {
-			assert_abs_diff_eq!(testee.value(x), 1., epsilon = ε);
-			assert_abs_diff_eq!(testee.value(-x), 1., epsilon = ε);
-		}
-		for x in itertools_num::linspace(-mul * testee.max, mul * testee.max, 100) {
-			// Currently, the derivative is discontinuous at ±min and ±max... let’s just skip these for now.
-			if x.abs().abs_diff_eq(&testee.min, δ) || x.abs().abs_diff_eq(&testee.max, δ) {
-				continue;
-			}
-			let df_num = (testee.value(x + δ) - testee.value(x - δ)) / (2. * δ);
-			let df_expl = testee.derivative(x);
-			assert!(abs_diff_eq!(df_expl, df_num, epsilon = ε), "At x={x}, df/dx:\nnumerical: {df_num}\nexplicit:  {df_expl}\n");
-		}
-	}
-
-	#[test]
-	fn test_quadratic_accelerator() {
-		let testee = super::QuadraticAccelerator { internal: 100.0, external: 150.0 };
-		let ε = 1.0e-4f32;
-		let δ = 1.0 / 8.0; // Mathematically, you want this to be small. Computationally, you don’t.
-		let margin = 1.0 / 16.0;
-		let mul = 1.0 + margin;
-		assert_abs_diff_eq!(testee.u(testee.external), testee.internal, epsilon = ε);
-		assert_abs_diff_eq!(testee.u(-testee.external), -testee.internal, epsilon = ε);
-		assert_abs_diff_eq!(testee.du(testee.external), 1., epsilon = ε);
-		assert_abs_diff_eq!(testee.du(-testee.external), 1., epsilon = ε);
-		for x in itertools_num::linspace(-mul * testee.external, mul * testee.external, 100) {
-			let ux = testee.u(x);
-			let xux = testee.x(ux);
-			assert!(abs_diff_eq!(x, xux, epsilon = ε), "At x={x}:\nu(x):    {ux}\nx(u(x)): {xux}\n");
-
-			let du_num = (testee.u(x + δ) - testee.u(x - δ)) / (2. * δ);
-			let du_expl = testee.du(x);
-			assert!(abs_diff_eq!(du_expl, du_num, epsilon = ε), "At x={x}, du/dx:\nnumerical: {du_num}\nexplicit:  {du_expl}\n");
-
-			let dudx = du_expl * testee.dx(ux);
-			assert!(abs_diff_eq!(dudx, 1.0, epsilon = ε), "At x={x}:\ndu/dx * dx/du: {dudx}\n");
-
-			let d2u_num = (testee.du(x + δ) - testee.du(x - δ)) / (2. * δ);
-			let d2u_expl = testee.d2u(x);
-			assert!(abs_diff_eq!(d2u_expl, d2u_num, epsilon = ε), "At x={x}, d^2u/dx^2:\nnumerical: {d2u_num}\nexplicit:  {d2u_expl}\n");
-		}
-		for u in itertools_num::linspace(-mul * testee.internal, mul * testee.internal, 100) {
-			let xu = testee.x(u);
-			let uxu = testee.u(xu);
-			assert!(abs_diff_eq!(u, uxu, epsilon = ε), "At u={u}:\nx(u):    {xu}\nu(x(u)): {uxu}\n");
-
-			let dx_num = (testee.x(u + δ) - testee.x(u - δ)) / (2. * δ);
-			let dx_expl = testee.dx(u);
-			assert!(abs_diff_eq!(dx_expl, dx_num, epsilon = ε), "At u={u}, dx/du:\nnumerical: {dx_num}\nexplicit:  {dx_expl}\n");
-
-			let dudx = testee.du(xu) * dx_expl;
-			assert!(abs_diff_eq!(dudx, 1.0, epsilon = ε), "At u={u}:\ndu/dx * dx/du: {dudx}\n");
-		}
-	}
-}

src/bin/flat/main.rs

@@ -1,25 +1,21 @@
 use std::f32::consts::{FRAC_PI_2, PI};
-use flo_draw::*;
+
 use flo_canvas::*;
+use flo_draw::*;
 use glam::*;
 
-mod riemann;
+use refraction::ifaces::{DebugTraceable, Traceable};
-mod fns;
+use refraction::tube::metric::Tube;
-mod float_fun;
+use refraction::tube::Space;
-mod tube;
+use refraction::types::{Location, Object, Ray};
-mod types;
+use refraction::utils::put_object;
+use refraction::DT;
 
-use riemann::{Metric, trace_iter};
+fn draw_loop(gc: &mut Vec<Draw>, mut pts: impl Iterator<Item = Vec3>) {
-use tube::Subspace::{Boundary, Inner, Outer};
-use tube::metric::Tube;
-use tube::Space;
-use types::{Location, Object, Ray};
-
-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; };
+	let Some(first) = pts.next() else {
+		return;
+	};
 	gc.move_to(first.x, first.y);
 	for pt in pts {
 		gc.line_to(pt.x, pt.y);
@@ -44,19 +40,26 @@ pub fn main() {
 				.enumerate()
 				.map(|(k, &y)| Object {
 					id: k as i32,
-					loc: {
+					loc: put_object(
-						let pos = vec2(0.0, y * tube.external_halflength);
+						&tube,
-						let adj: Mat2 = tube.sqrt_at(pos).inverse().into();
+						vec3(0.0, y * tube.external_halflength, 0.0),
-						let rot = Mat2::from_angle(y);
+						Mat3::from_mat2(Mat2::from_angle(y)),
-						Location {
+					),
-							pos,
-							rot: adj * rot,
-						}
-					},
 					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.)),
+			);
 
 			gc.canvas_height(500.0);
 			gc.transform(Transform2D::rotate(FRAC_PI_2));
@@ -64,13 +67,21 @@ pub fn main() {
 			gc.line_width(0.5);
 			// gc.stroke_color(Color::Rgba(1.0, 0.5, 0.0, 0.5));
 			// draw_fan(gc, &tube, vec2(-500.0, 0.0), vec2(1.0, 0.0), 1.0);
-			gc.stroke_color(Color::Rgba(1.0, 0.5, 0.0, 1.0));
+			gc.stroke_color(Color::Rgba(0.0, 0.8, 1.0, 1.0));
-			draw_fan_2(gc, &space, vec2(-500.0, 0.0), vec2(1.0, 0.0), 1.0);
+			draw_fan_2(gc, &space, cam3, 1.0);
 			gc.stroke_color(Color::Rgba(0.5, 1.0, 0.0, 1.0));
-			draw_fan_2(gc, &space, vec2(-2.5 * tube.outer_radius, 1.25 * tube.external_halflength), vec2(1.0, -1.0), 1.0);
+			draw_fan_2(gc, &space, cam2, 1.0);
+			gc.stroke_color(Color::Rgba(1.0, 0.5, 0.0, 1.0));
+			draw_fan_2(gc, &space, cam1, 1.0);
+
 			draw_track(gc, &space, vec2(-500.0, 0.0), vec2(1.0, 0.2));
 			draw_track(gc, &space, vec2(-500.0, 0.0), vec2(1.0, 0.5));
-			draw_track(gc, &space, vec2(-0.5 * tube.inner_radius, -1.25 * tube.external_halflength), vec2(0.1, 1.0));
+			draw_track(
+				gc,
+				&space,
+				vec2(-0.5 * tube.inner_radius, -1.25 * tube.external_halflength),
+				vec2(0.1, 1.0),
+			);
 
 			let circle_segments = 47;
 			for obj in &space.objs {
@@ -81,127 +92,94 @@ pub fn main() {
 				gc.fill();
 
 				gc.stroke_color(Color::Rgba(0.0, 0.0, 0.0, 0.5));
-				draw_loop(gc, itertools_num::linspace(0.0, 2.0 * PI, circle_segments).skip(1).map(|φ| {
+				draw_loop(
-					let dir = Vec2::from_angle(φ) * obj.r;
+					gc,
-					let dir = obj.loc.rot * dir;
+					itertools_num::linspace(0.0, 2.0 * PI, circle_segments)
-					pos + dir
+						.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
+						}),
+				);
 				gc.stroke_color(Color::Rgba(0.0, 0.5, 1.0, 0.5));
-				draw_loop(gc, itertools_num::linspace(0.0, 2.0 * PI, circle_segments).skip(1).map(|φ| {
+				draw_loop(
-					let dir = Vec2::from_angle(φ) * obj.r;
+					gc,
-					let dir = obj.loc.rot * dir;
+					itertools_num::linspace(0.0, 2.0 * PI, circle_segments)
-					space.trace_step(Ray { pos, dir }).pos
+						.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
+						}),
+				);
 				gc.stroke_color(Color::Rgba(0.5, 0.0, 1.0, 1.0));
-				draw_loop(gc, itertools_num::linspace(0.0, 2.0 * PI, circle_segments).skip(1).map(|φ| {
+				draw_loop(
-					let n = obj.r.floor();
+					gc,
-					let d = obj.r / n;
+					itertools_num::linspace(0.0, 2.0 * PI, circle_segments)
-					let dir = Vec2::from_angle(φ);
+						.skip(1)
-					let dir = obj.loc.rot * dir * d;
+						.map(|φ| {
-					space.trace_iter(Ray { pos, dir }).nth(n as usize).unwrap().pos
+							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 }).nth(n as usize).unwrap().pos
+						}),
+				);
 			}
 		});
 	});
 }
 
-fn draw_ray_2(gc: &mut Vec<Draw>, space: &Space, base: Vec2, dir: Vec2) {
+#[allow(dead_code)]
-	let mut hits = Vec::<Draw>::new();
+fn draw_cross(gc: &mut Vec<Draw>, pos: Vec2, r: f32) {
-	let dir = space.tube.globalize(base, dir);
+	gc.move_to(pos.x - r, pos.y - r);
-	gc.new_path();
+	gc.line_to(pos.x + r, pos.y + r);
-	gc.move_to(base.x, base.y);
+	gc.move_to(pos.x - r, pos.y + r);
-	let mut ray = Ray { pos: base, dir: space.tube.normalize_vec_at(base, dir) * DT };
+	gc.line_to(pos.x + r, pos.y - r);
-	for _ in 0..10000 {
+}
-		ray = space.trace_step(ray);
+
-		gc.line_to(ray.pos.x, ray.pos.y);
+fn draw_ray_2(gc: &mut Vec<Draw>, space: &Space, camera: Location, dir: Vec3) {
-		if ray.pos.abs().cmpgt(Vec2::splat(1000.0)).any() {
+	let pos = vec3(0., 0., 0.);
-			break;
+	let (hits, path) = space.trace_dbg(camera, Ray { pos, dir });
-		}
+	let hits2 = space.trace(camera, Ray { pos, dir });
-		let sub = space.which_subspace(ray.pos);
+	for (a, b) in hits.into_iter().zip(hits2.into_iter()) {
-		if sub == Boundary {
+		assert_eq!(a.id, b.id);
-			continue;
+		assert_eq!(a.pos, b.pos);
-		}
+		assert_eq!(a.rel, b.rel);
-		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];
-			hits.move_to(obj.loc.pos.x, obj.loc.pos.y);
-			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) - 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 = trace_iter(&space.tube, obj.loc.pos, obj.loc.rot * rel2, rel2.length() / 100.0).nth(100).unwrap();
-			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;
-		ray = match ret.end {
-			Some(r) => r,
-			None => {
-				ray = ray.forward(1000.0 / DT);
-				gc.line_to(ray.pos.x, ray.pos.y);
-				break;
-			}
-		};
-		for p in space.line(a, ray.pos, 10.0) {
-			gc.line_to(p.x, p.y);
-		}
-		gc.stroke();
-		gc.new_dash_pattern();
-		gc.new_path();
-		gc.move_to(ray.pos.x, ray.pos.y);
 	}
-	gc.stroke();
+
 	gc.new_path();
-	gc.new_dash_pattern();
+	gc.move_to(pos.x, pos.y);
-	gc.append(&mut hits);
+	for pt in &path.points[1..] {
+		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.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, base: Vec2, dir: Vec2, spread: f32) {
+fn draw_fan_2(gc: &mut Vec<Draw>, space: &Space, camera: Location, spread: f32) {
-	let dir = dir.normalize();
-	let v = vec2(-dir.y, dir.x);
 	for y in itertools_num::linspace(-spread, spread, 101) {
-		draw_ray_2(gc, space, base, dir + y * v);
+		draw_ray_2(gc, space, camera, vec3(1., y, 0.));
 	}
 }
 
-fn draw_ray(gc: &mut Vec<Draw>, space: &impl Metric, base: Vec2, dir: Vec2) {
-	let dir = space.globalize(base, dir);
-	gc.new_path();
-	gc.move_to(base.x, base.y);
-	for pt in trace_iter(space, base, dir, DT).take(10000) {
-		gc.line_to(pt.x, pt.y);
-		if pt.abs().cmpgt(Vec2::splat(1000.0)).any() {
-			break;
-		}
-	}
-	gc.stroke();
-}
-
 fn draw_track(gc: &mut Vec<Draw>, space: &Space, start: Vec2, dir: Vec2) {
 	const SCALE: f32 = 5.0;
 	const STEP: f32 = 2.0 * SCALE;
 	// let mut loc = Location { pos: start, rot: Mat2::IDENTITY };
 	// let dir = space.tube.globalize(start, dir);
 	// let v = space.tube.normalize(start, dir);
-	let mut loc = Location { pos: start, rot: mat2(dir, vec2(-dir.y, dir.x)) };
+	let mut loc = Location {
-	let v = vec2(1.0, 0.0);
+		pos: vec3(start.x, start.y, 0.),
+		rot: mat3(vec3(dir.x, dir.y, 0.), vec3(-dir.y, dir.x, 0.), vec3(0., 0., 1.)),
+	};
+	let v = vec3(1., 0., 0.);
 	let mut draw = |loc: &Location| {
 		let p = loc.pos;
 		let ax = p + loc.rot.x_axis * SCALE;
@@ -219,22 +197,14 @@ fn draw_track(gc: &mut Vec<Draw>, space: &Space, start: Vec2, dir: Vec2) {
 	};
 	draw(&loc);
 	for _ in 0..1000 {
-		let N = (STEP / DT).floor() as i32;
+		let n = (STEP / DT).floor() as i32;
-		for _ in 0..N {
+		for _ in 0..n {
 			loc = space.move_step(loc, v * DT);
 		}
 		draw(&loc);
 	}
 }
 
-fn draw_fan(gc: &mut Vec<Draw>, space: &impl Metric, base: Vec2, dir: Vec2, spread: f32) {
-	let dir = dir.normalize();
-	let v = vec2(-dir.y, dir.x);
-	for y in itertools_num::linspace(-spread, spread, 101) {
-		draw_ray(gc, space, base, dir + y * v);
-	}
-}
-
 trait Renderable {
 	fn render(&self, gc: &mut Vec<Draw>);
 }
@@ -242,8 +212,18 @@ trait Renderable {
 impl Renderable for Tube {
 	fn render(&self, gc: &mut Vec<Draw>) {
 		gc.new_path();
-		gc.rect(-self.outer_radius, -self.external_halflength, self.outer_radius, self.external_halflength);
+		gc.rect(
-		gc.rect(-self.inner_radius, -self.external_halflength, self.inner_radius, self.external_halflength);
+			-self.outer_radius,
+			-self.external_halflength,
+			self.outer_radius,
+			self.external_halflength,
+		);
+		gc.rect(
+			-self.inner_radius,
+			-self.external_halflength,
+			self.inner_radius,
+			self.external_halflength,
+		);
 		gc.winding_rule(WindingRule::EvenOdd);
 		gc.fill_color(Color::Rgba(0.8, 0.8, 0.8, 1.0));
 		gc.fill();

src/bin/flat/riemann.rs

@@ -1,148 +0,0 @@
-use glam::*;
-
-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(crate) 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
-	}
-}
-
-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)
-	}
-
-	fn globalize(&self, at: Vec2, v: Vec2) -> Vec2 {
-		Mat2::from(self.sqrt_at(at).inverse()) * v
-	}
-}
-
-pub struct TraceIter<'a, M: Metric> {
-	space: &'a M,
-	p: Vec2,
-	v: Vec2,
-	dt: f32,
-}
-
-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 = self.v + a * self.dt;
-		self.p = self.p + self.v * self.dt;
-		Some(self.p)
-	}
-}
-
-pub fn trace_iter<M: Metric>(space: &M, base: Vec2, dir: Vec2, dt: f32) -> TraceIter<M> {
-	TraceIter {
-		space,
-		p: base,
-		v: space.normalize_vec_at(base, dir),
-		dt,
-	}
-}
-
-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);
-}

src/bin/flat/tube/metric.rs

@@ -1,77 +0,0 @@
-use glam::{f32, Mat2, Vec2, vec2};
-use crate::fns;
-use crate::riemann::{Decomp2, Metric, Tens2};
-
-#[derive(Copy, Clone, Debug)]
-pub struct Tube {
-	pub outer_radius: f32,
-	pub inner_radius: f32,
-	pub external_halflength: f32,
-	pub internal_halflength: f32,
-}
-
-impl Tube {
-	fn fx(&self) -> fns::LinearLimiter { fns::LinearLimiter { min: self.inner_radius, max: self.outer_radius } }
-	fn fy(&self) -> fns::QuadraticAccelerator { fns::QuadraticAccelerator { internal: self.internal_halflength, external: self.external_halflength } }
-
-	pub fn y(&self, v: f32) -> f32 { self.fy().x(v) }
-	pub fn v(&self, y: f32) -> f32 { self.fy().u(y) }
-	pub fn dy(&self, v: f32) -> f32 { self.fy().dx(v) }
-	pub fn dv(&self, y: f32) -> f32 { self.fy().du(y) }
-}
-
-impl Metric for Tube {
-	fn sqrt_at(&self, pos: Vec2) -> Decomp2 {
-		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 {
-			ortho: Mat2::IDENTITY,
-			diag: vec2(1.0, s),
-		}
-	}
-
-	fn part_derivs_at(&self, pos: Vec2) -> Tens2 {
-		let sx = self.fx().value(pos.x);
-		let sy = self.fy().du(pos.y);
-		let s = sx + sy - sx * sy;
-		let dsx_dx = self.fx().derivative(pos.x);
-		let dsy_dy = self.fy().d2u(pos.y);
-		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]),
-			Mat2::from_cols_array(&[0.0, 0.0, 0.0, ds2_dy]),
-		]
-	}
-}
-
-#[test]
-fn test_tube_metric_derivs() {
-	struct Approx(Tube);
-	impl Metric for Approx {
-		fn sqrt_at(&self, pos: Vec2) -> Decomp2 { self.0.sqrt_at(pos) }
-	}
-	let testee = Tube {
-		inner_radius: 30.0,
-		outer_radius: 50.0,
-		internal_halflength: 100.0,
-		external_halflength: 300.0,
-	};
-	let approx = Approx(testee);
-	let epsilon = 1.0e-3;
-	let margin = 1.0 / 16.0;
-	let mul = 1.0 + margin;
-	for x in itertools_num::linspace(-mul * testee.outer_radius, mul * testee.outer_radius, 100) {
-		for y in itertools_num::linspace(-mul * testee.external_halflength, mul * testee.external_halflength, 100) {
-			let pos = vec2(x, y);
-			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");
-		}
-	}
-}

src/bin/flat/tube/mod.rs

@@ -1,327 +0,0 @@
-use glam::{bool, f32, Mat2, Vec2, vec2};
-use crate::riemann;
-use crate::riemann::Metric;
-use Subspace::{Boundary, Inner, Outer};
-use metric::Tube;
-use coords::{MapperInner, MapperOuter};
-use crate::types::{FlatTraceResult, Hit, Location, Object, Ray};
-
-pub mod metric;
-
-pub struct Space {
-	pub tube: Tube,
-	pub objs: Vec<Object>,
-}
-
-#[derive(PartialEq, Eq, Debug)]
-pub enum Subspace {
-	Outer,
-	Boundary,
-	Inner,
-}
-
-impl Space {
-	pub fn which_subspace(&self, pt: Vec2) -> Subspace {
-		if pt.y.abs() > self.tube.external_halflength {
-			Outer
-		} else if pt.x.abs() > self.tube.outer_radius {
-			Outer
-		} else if pt.x.abs() > self.tube.inner_radius {
-			Boundary
-		} else {
-			Inner
-		}
-	}
-
-	fn flat_to_global(&self, at: Vec2) -> Mat2 {
-		Mat2::from(self.tube.sqrt_at(at).inverse())
-	}
-
-	fn global_to_flat(&self, at: Vec2) -> Mat2 {
-		Mat2::from(self.tube.sqrt_at(at))
-	}
-
-	/// Выполняет один шаг трассировки. Работает в любой части пространства, но вне Boundary доступны более эффективные методы.
-	/// ray задаётся в основной СК.
-	pub fn trace_step(&self, ray: Ray) -> Ray {
-		let a: Vec2 = -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 }
-	}
-
-	/// Выполняет один шаг перемещения. Работает в любой части пространства.
-	/// off задаётся в локальной СК. Рекомендуется считать небольшими шагами.
-	pub fn move_step(&self, loc: Location, off: Vec2) -> Location {
-		let corr = Mat2::IDENTITY - riemann::contract(riemann::krist(&self.tube, loc.pos), loc.rot * off);
-		let p = loc.pos + corr * loc.rot * off;
-		Location { pos: p, rot: corr * loc.rot }
-	}
-
-	pub fn trace_iter(&self, ray: Ray) -> impl Iterator<Item=Ray> + '_ {
-		std::iter::successors(Some(ray), |&ray| Some(self.trace_step(ray)))
-	}
-
-	pub fn trace_inner(&self, ray: Ray) -> FlatTraceResult {
-		assert_eq!(self.which_subspace(ray.pos), Inner);
-		let inner = Rect { size: vec2(self.tube.inner_radius, self.tube.internal_halflength) };
-		let ray = self.tube.global_to_inner(ray);
-		assert!(inner.is_inside(ray.pos));
-		let dist = inner.trace_out_of(ray).expect("Can't get outta here!");
-		let objs = self.list_objects(|loc| self.tube.global_to_inner(loc));
-		FlatTraceResult {
-			end: Some(self.tube.inner_to_global(ray.forward(dist))),
-			objects: Self::hit_objects(objs.as_slice(), ray, Some(dist), |pos| self.tube.inner_to_global(pos)),
-		}
-	}
-
-	pub fn trace_outer(&self, ray: Ray) -> FlatTraceResult {
-		assert_eq!(self.which_subspace(ray.pos), Outer);
-		let outer = Rect { size: vec2(self.tube.outer_radius, self.tube.external_halflength) };
-		let dist = outer.trace_into(ray);
-		let objs = self.list_objects(|loc| self.tube.global_to_outer(loc));
-		FlatTraceResult {
-			end: dist.map(|dist| ray.forward(dist)),
-			objects: Self::hit_objects(objs.as_slice(), ray, dist, |pos| pos),
-		}
-	}
-
-	fn trace_boundary(&self, ray: Ray) -> Ray {
-		assert_eq!(self.which_subspace(ray.pos), Boundary);
-		self.trace_iter(ray)
-			.find(|&ray| self.which_subspace(ray.pos) != Boundary)
-			.expect("Can't get outta the wall!")
-	}
-
-	fn list_objects(&self, tfm: impl Fn(Location) -> Location) -> Vec<Object> {
-		self.objs.iter().map(|&Object { id, loc, r }| Object { id, loc: tfm(loc), r }).collect()
-	}
-
-	fn hit_objects(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) - 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))
-				} else {
-					None
-				}
-			})
-			.filter(|&(_, t)| t >= 0.0 && t < limit)
-			.map(|(obj, t)| {
-				let pos = ray.forward(t).pos;
-				let rel = obj.loc.rot.inverse() * Ray { pos: pos - obj.loc.pos, dir: ray.dir };
-				Hit { id: obj.id, distance: t, pos: globalize(pos), rel }
-			})
-			.collect()
-	}
-
-	pub fn line(&self, a: Vec2, b: Vec2, step: f32) -> Vec<Vec2> {
-		match self.which_subspace(a) {
-			Outer => vec![b],
-			Inner => {
-				let n = ((b - a).length() / step) as usize + 1;
-				let a = self.tube.global_to_inner(a);
-				let b = self.tube.global_to_inner(b);
-				(1..=n).map(|k| self.tube.inner_to_global(a.lerp(b, k as f32 / n as f32))).collect()
-			}
-			Boundary => panic!("Can't draw a line here!"),
-		}
-	}
-}
-
-struct Rect {
-	pub size: Vec2,
-}
-
-impl Rect {
-	/// Отражает луч, чтобы все координаты направления были положительны (допустимо благодаря симметрии Rect).
-	fn flip_ray(ray: Ray) -> Ray {
-		Ray { pos: ray.pos * ray.dir.signum(), dir: ray.dir.abs() }
-	}
-
-	fn is_inside(&self, pt: Vec2) -> bool {
-		pt.abs().cmplt(self.size).all()
-	}
-
-	fn trace_into(&self, ray: Ray) -> Option<f32> {
-		let ray = Self::flip_ray(ray);
-		// ray.pos.x + t * ray.dir.x = −size.x
-		let ts = (-self.size - ray.pos) / ray.dir;
-		let t = ts.max_element();
-		let pt = ray.pos + t * ray.dir;
-		if t < 0.0 { return None; }
-		if pt.cmpgt(self.size).any() { return None; }
-		Some(t)
-	}
-
-	fn trace_out_of(&self, ray: Ray) -> Option<f32> {
-		let ray = Self::flip_ray(ray);
-		// ray.pos.x + t * ray.dir.x = +size.x
-		let ts = (self.size - ray.pos) / ray.dir;
-		let t = ts.min_element();
-		Some(t)
-	}
-}
-
-#[test]
-fn test_rect() {
-	assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 5.0) }), Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 5.0) });
-	assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(-4.0, 5.0) }), Ray { pos: vec2(-2.0, 3.0), dir: vec2(4.0, 5.0) });
-	assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, -5.0) }), Ray { pos: vec2(2.0, -3.0), dir: vec2(4.0, 5.0) });
-	assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 0.0) }), Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 0.0) });
-
-	let r = Rect { size: vec2(2.0, 3.0) };
-
-	assert_eq!(r.trace_into(Ray { pos: vec2(3.0, 3.0), dir: vec2(1.0, 1.0) }), None);
-	assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 2.0), dir: vec2(1.0, 0.0) }), Some(1.0));
-	assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 2.0), dir: vec2(-1.0, 0.0) }), None);
-	assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 1.0), dir: vec2(2.0, 2.0) }), Some(0.5));
-	assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 2.1), dir: vec2(2.0, 2.0) }), None);
-
-	assert_eq!(r.trace_into(Ray { pos: vec2(2.0, 3.0), dir: vec2(1.0, 1.0) }), None);
-	assert_eq!(r.trace_into(Ray { pos: vec2(-2.0, 3.0), dir: vec2(-1.0, 1.0) }), None);
-	assert_eq!(r.trace_into(Ray { pos: vec2(2.0, 3.0), dir: vec2(-1.0, -1.0) }), Some(0.0));
-	assert_eq!(r.trace_into(Ray { pos: vec2(2.0, -3.0), dir: vec2(-1.0, 1.0) }), Some(0.0));
-
-	assert_eq!(r.trace_out_of(Ray { pos: vec2(0.0, 0.0), dir: vec2(1.0, 1.0) }), Some(2.0));
-	assert_eq!(r.trace_out_of(Ray { pos: vec2(0.0, 0.0), dir: vec2(0.0, 1.0) }), Some(3.0));
-	assert_eq!(r.trace_out_of(Ray { pos: vec2(0.0, 1.0), dir: vec2(0.0, -1.0) }), Some(4.0));
-	assert_eq!(r.trace_out_of(Ray { pos: vec2(1.0, 1.0), dir: vec2(0.0, -1.0) }), Some(4.0));
-	assert_eq!(r.trace_out_of(Ray { pos: vec2(2.0, 3.0), dir: vec2(1.0, 1.0) }), Some(0.0));
-}
-
-mod coords {
-	use glam::{Mat2, Vec2, vec2};
-	use crate::riemann::Metric;
-	use crate::types::{Location, Ray};
-	use super::{Rect, Tube};
-
-	pub trait MapperInner<T> {
-		fn inner_to_global(self, v: T) -> T;
-		fn global_to_inner(self, v: T) -> T;
-	}
-
-	pub trait MapperOuter<T> {
-		fn outer_to_global(self, v: T) -> T;
-		fn global_to_outer(self, v: T) -> T;
-	}
-
-	impl MapperInner<Vec2> for Tube {
-		fn inner_to_global(self, pos: Vec2) -> Vec2 {
-			vec2(pos.x, self.y(pos.y))
-		}
-
-		fn global_to_inner(self, pos: Vec2) -> Vec2 {
-			vec2(pos.x, self.v(pos.y))
-		}
-	}
-
-	impl MapperInner<Ray> for Tube {
-		fn inner_to_global(self, ray: Ray) -> Ray {
-			Ray {
-				pos: self.inner_to_global(ray.pos),
-				dir: vec2(ray.dir.x, self.dy(ray.pos.y) * ray.dir.y),
-			}
-		}
-
-		fn global_to_inner(self, ray: Ray) -> Ray {
-			Ray {
-				pos: self.global_to_inner(ray.pos),
-				dir: vec2(ray.dir.x, self.dv(ray.pos.y) * ray.dir.y),
-			}
-		}
-	}
-
-	impl MapperInner<Location> for Tube {
-		fn inner_to_global(self, loc: Location) -> Location {
-			todo!()
-		}
-
-		// NB: не работает для частей Outer с |y| < external_halflength. Но они и не нужны.
-		fn global_to_inner(self, loc: Location) -> Location {
-			Location {
-				pos: vec2(loc.pos.x, self.v(loc.pos.y)), // в плоской СК для Inner или её продолжении на Outer
-				rot: Mat2::from(self.sqrt_at(loc.pos)) * loc.rot,
-			}
-		}
-	}
-
-	impl MapperOuter<Location> for Tube {
-		fn outer_to_global(self, loc: Location) -> Location {
-			todo!()
-		}
-
-		// NB: имеет разрыв в области Inner на y = 0.
-		fn global_to_outer(self, loc: Location) -> Location {
-			let inner = Rect { size: vec2(self.inner_radius, self.external_halflength) };
-			if inner.is_inside(loc.pos) {
-				let Vec2 { x: u, y } = loc.pos; // в основной СК
-				let v = self.v(y) + y.signum() * (self.external_halflength - self.internal_halflength);
-				Location {
-					pos: vec2(u, v), // в плоском продолжении СК Outer на область Inner
-					rot: Mat2::from(self.sqrt_at(loc.pos)) * loc.rot,
-				}
-			} else {
-				loc
-			}
-		}
-	}
-
-	#[cfg(test)]
-	mod test {
-		use super::{Location, Tube, MapperOuter};
-		use glam::{Mat2, vec2};
-		use itertools_num::linspace;
-
-		#[test]
-		fn test_mapper_outer() {
-			let mapper = Tube {
-				inner_radius: 30.0,
-				outer_radius: 50.0,
-				internal_halflength: 100.0,
-				external_halflength: 300.0,
-			};
-			// straight
-			for x in linspace(-60., 60., 20) {
-				for y in linspace(-320., 320., 20) {
-					assert_eq!(mapper.global_to_outer(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_outer(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos;
-					let np = mapper.global_to_outer(Location { pos: vec2(-x, y), rot: Mat2::IDENTITY }).pos;
-					let pn = mapper.global_to_outer(Location { pos: vec2(x, -y), rot: Mat2::IDENTITY }).pos;
-					let nn = mapper.global_to_outer(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_outer(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_outer(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_outer(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos.y;
-					assert!(v > 200.0);
-					assert!(v > y);
-				}
-			}
-		}
-	}
-}

src/bin/mesh/main.rs

@@ -0,0 +1,153 @@
+use glam::*;
+use refraction::mesh_loader::load_mesh;
+use refraction::mesh_tracer::{trace_to_mesh, Mesh};
+use show_image::event::{ElementState, VirtualKeyCode, WindowEvent};
+use show_image::{exit, ImageInfo, ImageView, WindowOptions};
+use std::env;
+use std::error::Error;
+use std::f32::consts::PI;
+use std::fs::File;
+use std::io::BufReader;
+use std::sync::atomic::{AtomicUsize, Ordering::Relaxed};
+
+const W: i32 = 320;
+const H: i32 = 240;
+
+#[derive(Copy, Clone)]
+struct Color(u8, u8, u8);
+
+struct Image {
+	w: i32,
+	h: i32,
+	data: Vec<u8>,
+}
+
+impl Image {
+	fn data(&self) -> &[u8] {
+		self.data.as_slice()
+	}
+
+	fn put_pixel(&mut self, x: i32, y: i32, color: Color) {
+		if x < 0 || x >= self.w || y < 0 || y > self.h {
+			return;
+		}
+		let index = 3 * (x + self.w * y) as usize;
+		self.data[index] = color.0;
+		self.data[index + 1] = color.1;
+		self.data[index + 2] = color.2;
+	}
+}
+
+fn ypr_to_mat(ypr: Vec3) -> Mat3 {
+	let Vec3 {
+		x: yaw,
+		y: pitch,
+		z: roll,
+	} = ypr;
+	let m_roll = mat3(
+		vec3(roll.cos(), roll.sin(), 0.0),
+		vec3(-roll.sin(), roll.cos(), 0.0),
+		vec3(0.0, 0.0, 1.0),
+	);
+	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()),
+	);
+	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()),
+	);
+	m_roll * m_pitch * m_yaw
+}
+
+fn render(mesh: &Mesh, camera: impl Fn(Vec2) -> (Vec3, Vec3)) -> Image {
+	let bkg = vec3(0.0, 0.0, 0.0);
+	let mut img = Image {
+		w: W,
+		h: H,
+		data: vec![0; (3 * W * H) as usize],
+	};
+	let img_size = vec2(W as f32, H as f32);
+	for y in 0..H {
+		for x in 0..W {
+			let img_coords = vec2(x as f32, y as f32);
+			let off = (img_coords - img_size * 0.5) / img_size.y;
+			let (base, ray) = camera(off);
+			let color = if let Some(r) = trace_to_mesh(mesh, base, ray.normalize()) {
+				// to_vec3(0.45) * dot(r.normal, normalize(vec3(-1.0, 1.0, -1.0))) + 0.50
+				r.normal * 0.45 + 0.50
+			} else {
+				bkg
+			};
+			let color = (color * 255.0).as_ivec3().clamp(IVec3::splat(0), IVec3::splat(255));
+			img.put_pixel(x, y, Color(color.x as u8, color.y as u8, color.z as u8));
+		}
+	}
+	img
+}
+
+fn persp(dist: f32, off: Vec2) -> (Vec3, Vec3) {
+	(vec3(0., 0., -dist), vec3(off.x, off.y, dist))
+}
+
+fn ortho(dist: f32, off: Vec2) -> (Vec3, Vec3) {
+	(vec3(off.x, off.y, -dist), vec3(0., 0., 1.))
+}
+
+#[test]
+fn test_projs() {
+	fn check(f: fn(dist: f32, off: Vec2) -> (Vec3, Vec3), x: f32, y: f32, z: f32) {
+		let (base, ray) = f(z, vec2(x, y));
+		let at_dist = base + ray * (z / ray.z);
+		assert_eq!(at_dist, vec3(x, y, 0.));
+	}
+	check(persp, 1., 2., 3.);
+	check(ortho, 1., 2., 3.);
+	check(persp, 5., 3., 7.);
+	check(ortho, 9., 1., 8.);
+}
+
+// add_event_handler wants 'static + Send. Let it be so.
+static PROJ_INDEX: AtomicUsize = AtomicUsize::new(0);
+static PROJS: [fn(dist: f32, off: Vec2) -> (Vec3, Vec3); 2] = [persp, ortho];
+
+#[show_image::main]
+fn main() -> Result<(), Box<dyn Error>> {
+	let args: Vec<String> = env::args().collect();
+	if args.len() != 2 {
+		println!("Usage: {} path/to/model.obj", args[0]);
+		exit(1);
+	}
+	let mesh = {
+		let f = File::open(&args[1])?;
+		let mut f = BufReader::new(f);
+		load_mesh(&mut f)?
+	};
+	let window = show_image::create_window("Raytracing", WindowOptions::default())?;
+	window.add_event_handler(|_wnd, ev, _ctl| {
+		if let WindowEvent::KeyboardInput(ev) = ev {
+			if ev.input.state != ElementState::Pressed {
+				return;
+			}
+			if let Some(VirtualKeyCode::Tab) = ev.input.key_code {
+				PROJ_INDEX.store((PROJ_INDEX.load(Relaxed) + 1) % PROJS.len(), Relaxed);
+			}
+		}
+	})?;
+	loop {
+		for phi in 0..360 {
+			let proj = PROJS[PROJ_INDEX.load(Relaxed)];
+			let m_view = ypr_to_mat(vec3((135.0 + phi as f32) * PI / 180.0, -30.0 * PI / 180.0, 0.0f32));
+			let m_camera = m_view.transpose();
+			let img = render(mesh.as_slice(), |off| {
+				let (base, ray) = proj(40., 20. * off);
+				(m_camera * base, m_camera * ray)
+			});
+
+			let image = ImageView::new(ImageInfo::rgb8(W as u32, H as u32), img.data());
+			window.set_image("image", image)?;
+		}
+	}
+}

src/bin/wireframe/camera.rs

@@ -0,0 +1,123 @@
+use std::mem;
+
+use glam::{mat4, vec3, vec4, Mat4, Vec2, Vec3};
+
+#[repr(C)]
+#[derive(Copy, Clone, bytemuck::Pod, bytemuck::Zeroable)]
+struct CameraUniform {
+	mvp: [[f32; 4]; 4],
+	scale: [f32; 2],
+	pad: [u32; 2],
+}
+
+pub struct Camera {
+	buf: wgpu::Buffer,
+	bind: wgpu::BindGroup,
+	layout: wgpu::BindGroupLayout,
+}
+
+impl Camera {
+	pub fn new(device: &wgpu::Device) -> Camera {
+		let buf = 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 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 BindGroupLayout"),
+		});
+		let bind = device.create_bind_group(&wgpu::BindGroupDescriptor {
+			layout: &layout,
+			entries: &[wgpu::BindGroupEntry {
+				binding: 0,
+				resource: buf.as_entire_binding(),
+			}],
+			label: Some("Camera BindGroup"),
+		});
+		Camera { buf, bind, layout }
+	}
+
+	pub fn bind_group_layout(&self) -> &wgpu::BindGroupLayout {
+		&self.layout
+	}
+
+	pub fn bind_group(&self) -> &wgpu::BindGroup {
+		&self.bind
+	}
+
+	fn uniform(view_mtx: Mat4, view_size: Vec2) -> CameraUniform {
+		const M: Mat4 = mat4(
+			vec4(0., 0., 1., 0.),
+			vec4(-1., 0., 0., 0.),
+			vec4(0., 1., 0., 0.),
+			vec4(0., 0., 0., 1.),
+		);
+		let size = view_size.normalize() * std::f32::consts::SQRT_2;
+		let proj = make_proj_matrix(vec3(size.x, size.y, 2.), (1., (2f32).powi(16) + 1.)) * M;
+
+		let mvp = proj * view_mtx;
+		CameraUniform {
+			mvp: mvp.to_cols_array_2d(),
+			scale: (1. / size).to_array(),
+			pad: [0; 2],
+		}
+	}
+
+	pub fn set(&self, queue: &wgpu::Queue, view_mtx: Mat4, view_size: Vec2) {
+		let uniform = Self::uniform(view_mtx, view_size);
+		queue.write_buffer(&self.buf, 0, bytemuck::bytes_of(&uniform));
+	}
+}
+
+/// 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/lines.rs

@@ -0,0 +1,168 @@
+use std::mem;
+
+use bytemuck::{bytes_of, cast_slice, Pod, Zeroable};
+use glam::Vec3;
+use refraction::types::Ray;
+use wgpu::util::DeviceExt as _;
+
+#[repr(C)]
+#[derive(Copy, Clone, Debug, Pod, Zeroable)]
+struct Vertex {
+	pub position: [f32; 3],
+	pub tangent: [f32; 3],
+}
+
+#[repr(C)]
+#[derive(Copy, Clone, Pod, Zeroable)]
+struct PushConsts {
+	pub color: [f32; 3],
+	pub _pad: f32,
+}
+
+#[derive(Copy, Clone)]
+pub struct Attrs {
+	pub color: Vec3,
+}
+
+impl Attrs {
+	fn consts(&self) -> PushConsts {
+		PushConsts {
+			color: self.color.to_array(),
+			_pad: 0.,
+		}
+	}
+}
+
+pub struct Line {
+	consts: PushConsts,
+	npoints: u32,
+	buf: wgpu::Buffer,
+}
+
+impl Line {
+	pub fn new_strip(device: &wgpu::Device, attrs: Attrs, points: Vec<Ray>) -> Line {
+		let data: Vec<Vertex> = points
+			.into_iter()
+			.map(|r| Vertex {
+				position: r.pos.to_array(),
+				tangent: r.dir.to_array(),
+			})
+			.collect();
+		let buf = device.create_buffer_init(&wgpu::util::BufferInitDescriptor {
+			label: Some("Vertex Buffer"),
+			contents: cast_slice(&data),
+			usage: wgpu::BufferUsages::VERTEX,
+		});
+		Line {
+			consts: attrs.consts(),
+			npoints: data.len() as u32,
+			buf,
+		}
+	}
+}
+
+pub struct LineRenderer {
+	pipeline: wgpu::RenderPipeline,
+}
+
+static SHADER: &'static str = include_str!("ray.wgsl");
+
+impl LineRenderer {
+	pub fn new(
+		device: &wgpu::Device,
+		cam_layout: &wgpu::BindGroupLayout,
+		target_format: wgpu::TextureFormat,
+		depth_stencil: Option<wgpu::DepthStencilState>,
+		multisample: wgpu::MultisampleState,
+	) -> LineRenderer {
+		let shader = device.create_shader_module(wgpu::ShaderModuleDescriptor {
+			label: Some("Line Shader"),
+			source: wgpu::ShaderSource::Wgsl(SHADER.into()),
+		});
+
+		let consts_range = wgpu::PushConstantRange {
+			stages: wgpu::ShaderStages::VERTEX,
+			range: 0..mem::size_of::<PushConsts>() as u32,
+		};
+
+		let layout = device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor {
+			label: Some("Line RenderPipelineLayout"),
+			bind_group_layouts: &[cam_layout],
+			push_constant_ranges: &[consts_range],
+		});
+
+		let pipeline = device.create_render_pipeline(&wgpu::RenderPipelineDescriptor {
+			label: Some("Line RenderPipeline"),
+			layout: Some(&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: target_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,
+			multisample,
+			multiview: None,
+			cache: None,
+		});
+
+		LineRenderer { pipeline }
+	}
+
+	pub fn render<'a>(
+		&self,
+		pass: &mut wgpu::RenderPass,
+		cam_bind: &wgpu::BindGroup,
+		lines: impl Iterator<Item = &'a Line>,
+	) {
+		pass.set_pipeline(&self.pipeline);
+		pass.set_bind_group(0, cam_bind, &[]);
+		for line in lines {
+			pass.set_push_constants(wgpu::ShaderStages::VERTEX, 0, bytes_of(&line.consts));
+			pass.set_vertex_buffer(0, line.buf.slice(..));
+			pass.draw(0..4, 0..line.npoints - 1);
+		}
+	}
+}

src/bin/wireframe/main.rs

@@ -0,0 +1,417 @@
+use std::time::Instant;
+
+use glam::{uvec2, vec3, Vec3};
+use winit::{
+	event::*,
+	event_loop::EventLoop,
+	keyboard::{KeyCode, PhysicalKey},
+	window::{Window, WindowBuilder},
+};
+
+mod camera;
+mod lines;
+mod meshes;
+mod scene;
+mod viewport;
+
+// The coordinate system:
+// * X: forward
+// * Y: left
+// * Z: up
+
+fn prepare_scene(device: &wgpu::Device) -> (Vec<meshes::Mesh>, Vec<lines::Line>) {
+	let (meshes, lines) = scene::build();
+	let meshes = meshes
+		.into_iter()
+		.map(|mesh| meshes::Mesh::new_list(device, meshes::Attrs { color: mesh.color }, mesh.tris))
+		.collect();
+	let lines = lines
+		.into_iter()
+		.map(|line| lines::Line::new_strip(device, lines::Attrs { color: line.color }, line.pts))
+		.collect();
+	(meshes, lines)
+}
+
+#[cfg(any())]
+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::ZYX, ypr.x, ypr.y, ypr.z));
+		}
+
+		fn rotate_rel_quat(&mut self, rot: Quat) {
+			self.rot *= rot;
+		}
+	}
+}
+
+mod camctl {
+	use glam::{vec3, Mat4, Quat, Vec3};
+
+	pub struct CameraLocation {
+		pos: Vec3,
+		rot: Vec3,
+	}
+
+	fn rot_quat(rot: Vec3) -> Quat {
+		Quat::from_euler(glam::EulerRot::XYZ, rot.z, rot.y, rot.x)
+	}
+
+	impl CameraLocation {
+		pub fn new() -> CameraLocation {
+			let rot = vec3(std::f32::consts::FRAC_PI_4, 0., 0.);
+			let pos = rot_quat(rot) * vec3(-200., 0., 50.);
+			CameraLocation { pos, rot }
+		}
+
+		pub fn view_mtx(&self) -> Mat4 {
+			Mat4::from_quat(rot_quat(-self.rot)) * Mat4::from_translation(-self.pos)
+		}
+
+		pub fn move_rel(&mut self, offset: Vec3) {
+			self.pos += rot_quat(vec3(self.rot.x, 0., 0.)) * offset;
+		}
+
+		pub fn rotate_rel_ypr(&mut self, ypr: Vec3) {
+			self.rot += ypr;
+		}
+	}
+}
+
+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::KeyE), vec3(0., 0., 1.)),
+	(PhysicalKey::Code(KeyCode::KeyQ), vec3(0., 0., -1.)),
+	(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::Numpad4), vec3(1., 0., 0.)),
+	(PhysicalKey::Code(KeyCode::Numpad6), vec3(-1., 0., 0.)),
+	(PhysicalKey::Code(KeyCode::Numpad5), vec3(0., 1., 0.)),
+	(PhysicalKey::Code(KeyCode::Numpad8), vec3(0., -1., 0.)),
+	(PhysicalKey::Code(KeyCode::Numpad9), vec3(0., 0., 1.)),
+	(PhysicalKey::Code(KeyCode::Numpad7), vec3(0., 0., -1.)),
+];
+
+struct State<'a> {
+	device: wgpu::Device,
+	queue: wgpu::Queue,
+
+	fps: fps::Counter,
+	kbd: keyctl::Keyboard,
+	t1: Instant,
+
+	viewport: viewport::Viewport<'a>,
+	cam_loc: camctl::CameraLocation,
+	cam_obj: camera::Camera,
+	line_rend: lines::LineRenderer,
+	mesh_rend: meshes::Renderer,
+
+	lines: Vec<lines::Line>,
+	meshes: Vec<meshes::Mesh>,
+
+	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
+						| wgpu::Features::TEXTURE_ADAPTER_SPECIFIC_FORMAT_FEATURES,
+					required_limits: wgpu::Limits {
+						max_push_constant_size: 16,
+						..wgpu::Limits::default()
+					},
+					memory_hints: Default::default(),
+				},
+				None, // Trace path
+			)
+			.await
+			.unwrap();
+
+		let viewport = viewport::Viewport::new(&adapter, &device, surface, uvec2(size.width, size.height));
+
+		let kbd = keyctl::Keyboard::new();
+		let cam_loc = camctl::CameraLocation::new();
+		let t1 = Instant::now();
+
+		let cam_obj = camera::Camera::new(&device);
+		let line_rend = lines::LineRenderer::new(
+			&device,
+			cam_obj.bind_group_layout(),
+			viewport.format(),
+			Some(wgpu::DepthStencilState {
+				format: wgpu::TextureFormat::Depth24Plus,
+				depth_write_enabled: false,
+				depth_compare: wgpu::CompareFunction::LessEqual,
+				stencil: wgpu::StencilState::default(),
+				bias: wgpu::DepthBiasState::default(),
+			}),
+			wgpu::MultisampleState {
+				count: viewport.sample_count(),
+				mask: !0,
+				alpha_to_coverage_enabled: false,
+			},
+		);
+		let mesh_rend = meshes::Renderer::new(
+			&device,
+			cam_obj.bind_group_layout(),
+			viewport.format(),
+			Some(wgpu::DepthStencilState {
+				format: wgpu::TextureFormat::Depth24Plus,
+				depth_write_enabled: true,
+				depth_compare: wgpu::CompareFunction::LessEqual,
+				stencil: wgpu::StencilState::default(),
+				bias: wgpu::DepthBiasState::default(),
+			}),
+			wgpu::MultisampleState {
+				count: viewport.sample_count(),
+				mask: !0,
+				alpha_to_coverage_enabled: false,
+			},
+		);
+
+		let (meshes, lines) = prepare_scene(&device);
+
+		let fps = fps::Counter::new();
+
+		Self {
+			device,
+			queue,
+			viewport,
+			line_rend,
+			mesh_rend,
+			kbd,
+			fps,
+			cam_loc,
+			cam_obj,
+			t1,
+			lines,
+			meshes,
+			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.viewport
+				.resize(&self.device, uvec2(new_size.width, new_size.height));
+		}
+	}
+
+	fn update(&mut self) {
+		let dt = {
+			let t2 = Instant::now();
+			let dt = t2 - self.t1;
+			self.t1 = t2;
+			dt.as_secs_f32()
+		};
+		let size = self.viewport.size().as_vec2();
+
+		self.cam_loc.move_rel(100. * dt * self.kbd.control(&KEYS_MOVE));
+		self.cam_loc.rotate_rel_ypr(2. * dt * self.kbd.control(&KEYS_ROTATE));
+		self.cam_obj.set(&self.queue, self.cam_loc.view_mtx(), size);
+	}
+
+	fn render(&mut self) -> Result<(), wgpu::SurfaceError> {
+		self.fps.on_frame();
+		self.window
+			.set_title(&format!("Space Refraction ({:.1} FPS)", self.fps.get()));
+		self.viewport
+			.render_single_pass(&self.device, &self.queue, |mut render_pass| {
+				self.mesh_rend
+					.render(&mut render_pass, self.cam_obj.bind_group(), self.meshes.iter());
+				self.line_rend
+					.render(&mut render_pass, self.cam_obj.bind_group(), self.lines.iter());
+			})
+	}
+}
+
+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.viewport.configure(&state.device);
+								}
+								// 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());
+}
+
+mod fps {
+	use std::time::{Duration, Instant};
+
+	pub struct Counter {
+		fps: f32,
+		t1: Instant,
+		frames: u32,
+	}
+
+	impl Counter {
+		pub fn new() -> Self {
+			Self {
+				fps: 0.,
+				t1: Instant::now(),
+				frames: 0,
+			}
+		}
+
+		pub fn get(&self) -> f32 {
+			self.fps
+		}
+
+		pub fn on_frame(&mut self) {
+			self.frames += 1;
+			let t2 = Instant::now();
+			let dt = t2 - self.t1;
+			if dt >= Duration::from_secs(1) {
+				*self = Self {
+					fps: self.frames as f32 / dt.as_secs_f32(),
+					t1: t2,
+					frames: 0,
+				}
+			}
+		}
+	}
+}

src/bin/wireframe/mesh.wgsl

@@ -0,0 +1,63 @@
+struct CameraUniform {
+	mvp: mat4x4<f32>,
+	scale: vec2<f32>,
+}
+@group(0) @binding(0)
+var<uniform> camera: CameraUniform;
+
+struct MeshUniform {
+	color: vec4<f32>,
+}
+var<push_constant> mesh: MeshUniform;
+
+struct VertexInput {
+	@location(0) position: vec3<f32>,
+	@location(1) normal: vec3<f32>,
+}
+
+struct VertexOutput {
+	@builtin(position) clip_position: vec4<f32>,
+	@location(0) vertex_color: vec4<f32>,
+}
+
+struct FragmentOutput {
+	@location(0) color: vec4<f32>,
+	@builtin(sample_mask) mask: u32,
+}
+
+fn hash(key : u32) -> u32 {
+	var v = key;
+	v *= 0xb384af1bu;
+	v ^= v >> 15u;
+	return v;
+}
+
+@vertex
+fn vs_main(ver: VertexInput) -> VertexOutput {
+	var out: VertexOutput;
+	var light = dot(ver.normal, normalize(vec3(1., 1., 1.)));
+	light = .7 + .3 * light;
+	out.vertex_color = vec4(light * mesh.color.xyz, mesh.color.w);
+	out.clip_position = camera.mvp * vec4(ver.position, 1.);
+	return out;
+}
+
+@fragment
+fn fs_main(in: VertexOutput) -> FragmentOutput {
+	var out: FragmentOutput;
+	out.color = vec4(in.vertex_color.xyz, 1.);
+	var x = bitcast<u32>(in.clip_position.x);
+	var y = bitcast<u32>(in.clip_position.y);
+	var z = bitcast<u32>(in.clip_position.z);
+	var alpha = in.vertex_color.w;
+	var seed = hash(hash(hash(x) ^ y) ^ z);
+	var mask = 0u;
+	for (var sample = 0u; sample < 8u; sample++) {
+		var threshold = f32(hash(seed ^ sample)) / 0x1p32;
+		if (alpha > threshold) {
+			mask |= 1u << sample;
+		}
+	}
+	out.mask = mask;
+	return out;
+}

src/bin/wireframe/meshes.rs

@@ -0,0 +1,160 @@
+use std::mem;
+
+use bytemuck::{bytes_of, cast_slice, Pod, Zeroable};
+use glam::{Vec3, Vec4};
+use wgpu::util::DeviceExt as _;
+
+use crate::scene::Face;
+
+#[repr(C)]
+#[derive(Copy, Clone, Debug, Pod, Zeroable)]
+struct Vertex {
+	pub position: [f32; 3],
+	pub normal: [f32; 3],
+}
+
+impl From<crate::scene::Vertex> for Vertex {
+	fn from(value: crate::scene::Vertex) -> Self {
+		Vertex {
+			position: value.position.into(),
+			normal: value.normal.into(),
+		}
+	}
+}
+
+#[repr(C)]
+#[derive(Copy, Clone, Pod, Zeroable)]
+struct PushConsts {
+	pub color: [f32; 4],
+}
+
+#[derive(Copy, Clone)]
+pub struct Attrs {
+	pub color: Vec4,
+}
+
+impl Attrs {
+	fn consts(&self) -> PushConsts {
+		PushConsts {
+			color: self.color.to_array(),
+		}
+	}
+}
+
+pub struct Mesh {
+	consts: PushConsts,
+	npoints: u32,
+	buf: wgpu::Buffer,
+}
+
+impl Mesh {
+	pub fn new_list(device: &wgpu::Device, attrs: Attrs, tris: Vec<Face>) -> Self {
+		let data: Vec<Vertex> = tris.into_iter().flat_map(|face| face.map(Into::into)).collect();
+		let buf = device.create_buffer_init(&wgpu::util::BufferInitDescriptor {
+			label: Some("Mesh Vertex Buffer"),
+			contents: cast_slice(&data),
+			usage: wgpu::BufferUsages::VERTEX,
+		});
+		Mesh {
+			consts: attrs.consts(),
+			npoints: data.len() as u32,
+			buf,
+		}
+	}
+}
+
+pub struct Renderer {
+	pipeline: wgpu::RenderPipeline,
+}
+
+static SHADER: &'static str = include_str!("mesh.wgsl");
+
+impl Renderer {
+	pub fn new(
+		device: &wgpu::Device,
+		cam_layout: &wgpu::BindGroupLayout,
+		target_format: wgpu::TextureFormat,
+		depth_stencil: Option<wgpu::DepthStencilState>,
+		multisample: wgpu::MultisampleState,
+	) -> Renderer {
+		let shader = device.create_shader_module(wgpu::ShaderModuleDescriptor {
+			label: Some("Mesh Shader"),
+			source: wgpu::ShaderSource::Wgsl(SHADER.into()),
+		});
+
+		let consts_range = wgpu::PushConstantRange {
+			stages: wgpu::ShaderStages::VERTEX,
+			range: 0..mem::size_of::<PushConsts>() as u32,
+		};
+
+		let layout = device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor {
+			label: Some("Mesh RenderPipelineLayout"),
+			bind_group_layouts: &[cam_layout],
+			push_constant_ranges: &[consts_range],
+		});
+
+		let pipeline = device.create_render_pipeline(&wgpu::RenderPipelineDescriptor {
+			label: Some("Mesh RenderPipeline"),
+			layout: Some(&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::Vertex,
+					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, normal) as u64,
+							shader_location: 1,
+							format: wgpu::VertexFormat::Float32x3,
+						},
+					],
+				}],
+				compilation_options: Default::default(),
+			},
+			fragment: Some(wgpu::FragmentState {
+				module: &shader,
+				entry_point: "fs_main",
+				targets: &[Some(wgpu::ColorTargetState {
+					format: target_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::TriangleList,
+				..Default::default()
+			},
+			depth_stencil,
+			multisample,
+			multiview: None,
+			cache: None,
+		});
+
+		Renderer { pipeline }
+	}
+
+	pub fn render<'a>(
+		&self,
+		pass: &mut wgpu::RenderPass,
+		cam_bind: &wgpu::BindGroup,
+		meshes: impl Iterator<Item = &'a Mesh>,
+	) {
+		pass.set_pipeline(&self.pipeline);
+		pass.set_bind_group(0, cam_bind, &[]);
+		for mesh in meshes {
+			pass.set_push_constants(wgpu::ShaderStages::VERTEX, 0, bytes_of(&mesh.consts));
+			pass.set_vertex_buffer(0, mesh.buf.slice(..));
+			pass.draw(0..mesh.npoints, 0..1);
+		}
+	}
+}

src/bin/wireframe/ray.wgsl

@@ -0,0 +1,61 @@
+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_cs0 = camera.mvp * vec4(dir, 0.);
+	let dir_cs = dir_cs0.xyz * pt_cs.w - pt_cs.xyz * dir_cs0.w;
+	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,210 @@
+use glam::*;
+use itertools::chain;
+
+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>,
+}
+
+pub struct Vertex {
+	pub position: Vec3,
+	pub normal: Vec3,
+}
+
+pub type Face = [Vertex; 3];
+
+pub enum Mesh {
+	List(Vec<Face>),
+}
+
+pub struct FancyMesh {
+	pub color: Vec4,
+	pub tris: Vec<Face>,
+}
+
+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_mark(pos: Vec3) -> Vec<Line> {
+	[
+		vec3(1., 1., 1.),
+		vec3(1., 1., -1.),
+		vec3(1., -1., 1.),
+		vec3(1., -1., -1.),
+	]
+	.into_iter()
+	.map(|off| draw_line(pos - off, pos + off))
+	.collect()
+}
+
+pub fn build() -> (Vec<FancyMesh>, 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 cam2l = put_object(
+		&space.tube,
+		vec3(-2.5 * tube.outer_radius, 1.25 * tube.external_halflength, 0.),
+		mat3(vec3(1., -0.825, 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![];
+	gc.push(FancyMesh {
+		color: vec4(0.10, 0.12, 0.15, 0.8),
+		tris: tube.render(),
+	});
+	let meshes = gc;
+
+	let mut gc = vec![];
+	paint(&mut gc, vec3(0.0, 0.6, 1.0), draw_fan_2(&space, cam3, vec3(0., 1., 0.)));
+	paint(&mut gc, vec3(0.2, 1.0, 0.0), draw_fan_2(&space, cam2, vec3(0., 1., 0.)));
+	paint(
+		&mut gc,
+		vec3(0.0, 1.0, 0.6),
+		draw_fan_2(&space, cam2l, vec3(0., 0., 1.)),
+	);
+	paint(&mut gc, vec3(1.0, 0.2, 0.0), draw_fan_2(&space, cam1, vec3(0., 1., 0.)));
+	let lines = gc;
+
+	(meshes, lines)
+}
+
+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 });
+	if true {
+		let hits2 = space.trace(camera, Ray { pos, dir });
+		for (a, b) in hits.iter().zip(hits2.into_iter()) {
+			assert_eq!(a.id, b.id);
+			assert_eq!(a.pos, b.pos);
+			assert_eq!(a.rel, b.rel);
+		}
+	}
+
+	for hit in hits {
+		gc.extend(draw_mark(hit.pos));
+	}
+
+	let mut pts = path.points;
+	let end_pos = *pts.last().expect("the starting point is always in the path");
+	pts.extend(itertools::iterate(end_pos, |r| r.forward(100.0)).take(1000));
+	gc.push(Line::Strip(pts));
+}
+
+fn draw_fan_2(space: &Space, camera: Location, spread: Vec3) -> Vec<Line> {
+	let mut gc = vec![];
+	for δ in itertools_num::linspace(-1., 1., 101) {
+		draw_ray_2(&mut gc, space, camera, vec3(1., 0., 0.) + δ * spread);
+	}
+	gc
+}
+
+trait Renderable {
+	fn render(&self) -> Vec<Face>;
+}
+
+impl Renderable for Tube {
+	fn render(&self) -> Vec<Face> {
+		let sides = 42;
+		let step = 2. * std::f32::consts::PI / sides as f32;
+		let dir = |k| {
+			let d = Vec2::from_angle(k as f32 * step);
+			vec3(d.x, 0., d.y)
+		};
+
+		let side = vec3(0., self.external_halflength, 0.);
+		let r1 = self.inner_radius;
+		let r2 = self.outer_radius;
+
+		let vertex = |k, r, h| {
+			let n = dir(k);
+			Vertex {
+				position: r * n + h * side,
+				normal: n,
+			}
+		};
+		let inner = (0..sides).flat_map(|k| {
+			[
+				[vertex(k, -r1, -1.), vertex(k, -r1, 1.), vertex(k + 1, -r1, -1.)],
+				[vertex(k + 1, -r1, -1.), vertex(k, -r1, 1.), vertex(k + 1, -r1, 1.)],
+			]
+		});
+		let outer = (0..sides).flat_map(|k| {
+			[
+				[vertex(k, r2, -1.), vertex(k + 1, r2, -1.), vertex(k, r2, 1.)],
+				[vertex(k + 1, r2, -1.), vertex(k + 1, r2, 1.), vertex(k, r2, 1.)],
+			]
+		});
+
+		let vertex = |k, r, h| {
+			let n = dir(k);
+			Vertex {
+				position: r * n + h * side,
+				normal: h * Vec3::Y,
+			}
+		};
+		let cap = (0..sides).flat_map(|k| {
+			[
+				[vertex(k, r1, -1.), vertex(k + 1, r1, -1.), vertex(k + 1, r2, -1.)],
+				[vertex(k, r1, -1.), vertex(k + 1, r2, -1.), vertex(k, r2, -1.)],
+				[vertex(k, r1, 1.), vertex(k + 1, r1, 1.), vertex(k + 1, r2, 1.)],
+				[vertex(k, r1, 1.), vertex(k + 1, r2, 1.), vertex(k, r2, 1.)],
+			]
+		});
+		chain!(inner, outer, cap).collect()
+	}
+}

src/bin/wireframe/viewport.rs

@@ -0,0 +1,147 @@
+use glam::{uvec2, UVec2};
+
+pub struct Viewport<'a> {
+	surface: wgpu::Surface<'a>,
+	config: wgpu::SurfaceConfiguration,
+	sample_count: u32,
+	multisample: Multisample,
+}
+
+impl<'a> Viewport<'a> {
+	pub fn new(adapter: &wgpu::Adapter, device: &wgpu::Device, surface: wgpu::Surface<'a>, size: UVec2) -> Self {
+		let caps = surface.get_capabilities(adapter);
+		let format = wgpu::TextureFormat::Bgra8UnormSrgb;
+		let sample_count = adapter
+			.get_texture_format_features(format)
+			.flags
+			.supported_sample_counts()
+			.into_iter()
+			.max()
+			.unwrap();
+		eprintln!("Using x{sample_count} mutlisampling");
+		let config = wgpu::SurfaceConfiguration {
+			usage: wgpu::TextureUsages::RENDER_ATTACHMENT,
+			format,
+			width: size.x,
+			height: size.y,
+			present_mode: wgpu::PresentMode::Fifo,
+			alpha_mode: caps.alpha_modes[0],
+			view_formats: vec![],
+			desired_maximum_frame_latency: 2,
+		};
+		let multisample = Multisample::new(device, format, size, sample_count);
+		Self {
+			surface,
+			config,
+			sample_count,
+			multisample,
+		}
+	}
+
+	pub fn configure(&mut self, device: &wgpu::Device) {
+		self.surface.configure(&device, &self.config);
+		self.multisample = Multisample::new(device, self.format(), self.size(), self.sample_count);
+	}
+
+	pub fn resize(&mut self, device: &wgpu::Device, size: UVec2) {
+		self.config.width = size.x;
+		self.config.height = size.y;
+		self.configure(&device);
+	}
+
+	pub fn size(&self) -> UVec2 {
+		uvec2(self.config.width, self.config.height)
+	}
+
+	pub fn format(&self) -> wgpu::TextureFormat {
+		self.config.format
+	}
+
+	pub fn sample_count(&self) -> u32 {
+		self.sample_count
+	}
+
+	pub fn render_single_pass(
+		&mut self,
+		device: &wgpu::Device,
+		queue: &wgpu::Queue,
+		f: impl FnOnce(wgpu::RenderPass),
+	) -> Result<(), wgpu::SurfaceError> {
+		let output = self.surface.get_current_texture()?;
+		let view = output.texture.create_view(&wgpu::TextureViewDescriptor::default());
+		let mut encoder = device.create_command_encoder(&wgpu::CommandEncoderDescriptor {
+			label: Some("Render CommandEncoder"),
+		});
+		let render_pass = encoder.begin_render_pass(&wgpu::RenderPassDescriptor {
+			label: Some("RenderPass"),
+			color_attachments: &[Some(wgpu::RenderPassColorAttachment {
+				view: &self.multisample.color,
+				resolve_target: Some(&view),
+				ops: wgpu::Operations {
+					load: wgpu::LoadOp::Clear(wgpu::Color {
+						r: 0.,
+						g: 0.,
+						b: 0.,
+						a: 1.,
+					}),
+					store: wgpu::StoreOp::Store,
+				},
+			})],
+			depth_stencil_attachment: Some(wgpu::RenderPassDepthStencilAttachment {
+				view: &self.multisample.depth,
+				depth_ops: Some(wgpu::Operations {
+					load: wgpu::LoadOp::Clear(1.),
+					store: wgpu::StoreOp::Discard,
+				}),
+				stencil_ops: None,
+			}),
+			occlusion_query_set: None,
+			timestamp_writes: None,
+		});
+		f(render_pass);
+		queue.submit(std::iter::once(encoder.finish()));
+		output.present();
+		Ok(())
+	}
+}
+
+struct Multisample {
+	color: wgpu::TextureView,
+	depth: wgpu::TextureView,
+}
+
+impl Multisample {
+	fn new(device: &wgpu::Device, format: wgpu::TextureFormat, size: UVec2, sample_count: u32) -> Multisample {
+		let color = device.create_texture(&wgpu::TextureDescriptor {
+			label: Some("Multisample color texture"),
+			size: wgpu::Extent3d {
+				width: size.x,
+				height: size.y,
+				depth_or_array_layers: 1,
+			},
+			mip_level_count: 1,
+			sample_count,
+			dimension: wgpu::TextureDimension::D2,
+			format,
+			usage: wgpu::TextureUsages::RENDER_ATTACHMENT,
+			view_formats: &[],
+		});
+		let color = color.create_view(&wgpu::TextureViewDescriptor::default());
+		let depth = device.create_texture(&wgpu::TextureDescriptor {
+			label: Some("Multisample depth texture"),
+			size: wgpu::Extent3d {
+				width: size.x,
+				height: size.y,
+				depth_or_array_layers: 1,
+			},
+			mip_level_count: 1,
+			sample_count,
+			dimension: wgpu::TextureDimension::D2,
+			format: wgpu::TextureFormat::Depth24Plus,
+			usage: wgpu::TextureUsages::RENDER_ATTACHMENT,
+			view_formats: &[],
+		});
+		let depth = depth.create_view(&wgpu::TextureViewDescriptor::default());
+		Multisample { color, depth }
+	}
+}

src/fns.rs

@@ -0,0 +1,213 @@
+//! Functions used to make metrics.
+
+use crate::mathx::FloatExt2;
+
+pub trait Limiter {
+	fn value(&self, x: f32) -> f32;
+	fn derivative(&self, x: f32) -> f32;
+}
+
+pub struct SmoothstepLimiter {
+	pub min: f32,
+	pub max: f32,
+}
+
+impl Limiter for SmoothstepLimiter {
+	fn value(&self, x: f32) -> f32 {
+		let y = (self.min, self.max).inverse_lerp(x.abs()).clamp(0.0, 1.0);
+		3.0 * y * y - 2.0 * y * y * y
+	}
+	fn derivative(&self, x: f32) -> f32 {
+		if x.abs() > self.min && x.abs() < self.max {
+			let t = (self.min, self.max).inverse_lerp(x.abs());
+			6.0 * x.signum() * t * (1.0 - t) / (self.max - self.min)
+		} else {
+			0.0
+		}
+	}
+}
+
+pub struct SmootherstepLimiter {
+	pub min: f32,
+	pub max: f32,
+}
+
+impl Limiter for SmootherstepLimiter {
+	fn value(&self, x: f32) -> f32 {
+		let y = (self.min, self.max).inverse_lerp(x.abs()).clamp(0.0, 1.0);
+		6.0 * y.powi(5) - 15.0 * y.powi(4) + 10.0 * y.powi(3)
+	}
+	fn derivative(&self, x: f32) -> f32 {
+		if x.abs() > self.min && x.abs() < self.max {
+			let t = (self.min, self.max).inverse_lerp(x.abs());
+			30.0 * (t * (1.0 - t)).powi(2) * x.signum() / (self.max - self.min)
+		} else {
+			0.0
+		}
+	}
+}
+
+pub struct QuadraticAccelerator {
+	pub internal: f32,
+	pub external: f32,
+}
+
+/// Продолжает функцию f с [-lim, lim] линейно в предположении f(±lim) = ±val, f'(±lim) = 1.
+fn extend_linear(t: f32, f: impl FnOnce(f32) -> f32, lim: f32, val: f32) -> f32 {
+	if t.abs() <= lim {
+		f(t)
+	} else {
+		t + t.signum() * (val - lim)
+	}
+}
+
+/// Продолжает функцию f с [-lim, lim] константой в предположении f(±lim) = val, f'(±lim) = 0.
+fn extend_const(t: f32, f: impl FnOnce(f32) -> f32, lim: f32, val: f32) -> f32 {
+	if t.abs() <= lim {
+		f(t)
+	} else {
+		val
+	}
+}
+
+impl QuadraticAccelerator {
+	fn a(&self) -> f32 {
+		-(self.external - self.internal) / self.internal.powi(2)
+	}
+	fn b(&self) -> f32 {
+		2.0 * self.external / self.internal - 1.0
+	}
+	fn root(&self, x: f32) -> f32 {
+		(self.b().powi(2) + 4.0 * self.a() * x.abs()).sqrt()
+	}
+
+	pub fn x(&self, u: f32) -> f32 {
+		extend_linear(u, |u| (self.a() * u.abs() + self.b()) * u, self.internal, self.external)
+	}
+	pub fn u(&self, x: f32) -> f32 {
+		extend_linear(
+			x,
+			|x| 0.5 * x.signum() * (-self.b() + self.root(x)) / self.a(),
+			self.external,
+			self.internal,
+		)
+	}
+	pub fn dx(&self, u: f32) -> f32 {
+		extend_const(u, |u| 2.0 * self.a() * u.abs() + self.b(), self.internal, 1.0)
+	}
+	pub fn du(&self, x: f32) -> f32 {
+		extend_const(x, |x| 1.0 / self.root(x), self.external, 1.0)
+	}
+	pub fn d2u(&self, x: f32) -> f32 {
+		extend_const(
+			x,
+			|x| -2.0 * x.signum() * self.a() * self.root(x).powi(-3),
+			self.external,
+			0.0,
+		)
+	}
+}
+
+#[cfg(test)]
+mod test {
+	use approx::{abs_diff_eq, assert_abs_diff_eq};
+
+	use super::*;
+
+	fn test_limiter(testee: impl Limiter, min: f32, max: f32, δ: f32) {
+		let ε = 1.0e-4f32;
+		let margin = 1.0 / 16.0;
+		let mul = 1.0 + margin;
+		for x in itertools_num::linspace(0., min, 10) {
+			assert_abs_diff_eq!(testee.value(x), 0., epsilon = ε);
+			assert_abs_diff_eq!(testee.value(-x), 0., epsilon = ε);
+		}
+		for x in itertools_num::linspace(max, mul * max, 10) {
+			assert_abs_diff_eq!(testee.value(x), 1., epsilon = ε);
+			assert_abs_diff_eq!(testee.value(-x), 1., epsilon = ε);
+		}
+		for x in itertools_num::linspace(-mul * max, mul * max, 100) {
+			let df_num = (testee.value(x + δ) - testee.value(x - δ)) / (2. * δ);
+			let df_expl = testee.derivative(x);
+			assert!(
+				abs_diff_eq!(df_expl, df_num, epsilon = ε),
+				"At x={x}, df/dx:\nnumerical: {df_num}\nexplicit:  {df_expl}\n"
+			);
+		}
+	}
+
+	#[test]
+	fn test_smoothstep_limiter() {
+		test_limiter(SmoothstepLimiter { min: 20.0, max: 30.0 }, 20.0, 30.0, 1.0 / 32.0);
+	}
+
+	#[test]
+	fn test_smootherstep_limiter() {
+		test_limiter(SmootherstepLimiter { min: 20.0, max: 30.0 }, 20.0, 30.0, 1.0 / 32.0);
+	}
+
+	#[test]
+	fn test_quadratic_accelerator() {
+		let testee = super::QuadraticAccelerator {
+			internal: 100.0,
+			external: 150.0,
+		};
+		let ε = 1.0e-4f32;
+		let δ = 1.0 / 8.0; // Mathematically, you want this to be small. Computationally, you don’t.
+		let margin = 1.0 / 16.0;
+		let mul = 1.0 + margin;
+		assert_abs_diff_eq!(testee.u(testee.external), testee.internal, epsilon = ε);
+		assert_abs_diff_eq!(testee.u(-testee.external), -testee.internal, epsilon = ε);
+		assert_abs_diff_eq!(testee.du(testee.external), 1., epsilon = ε);
+		assert_abs_diff_eq!(testee.du(-testee.external), 1., epsilon = ε);
+		for x in itertools_num::linspace(-mul * testee.external, mul * testee.external, 100) {
+			let ux = testee.u(x);
+			let xux = testee.x(ux);
+			assert!(
+				abs_diff_eq!(x, xux, epsilon = ε),
+				"At x={x}:\nu(x):    {ux}\nx(u(x)): {xux}\n"
+			);
+
+			let du_num = (testee.u(x + δ) - testee.u(x - δ)) / (2. * δ);
+			let du_expl = testee.du(x);
+			assert!(
+				abs_diff_eq!(du_expl, du_num, epsilon = ε),
+				"At x={x}, du/dx:\nnumerical: {du_num}\nexplicit:  {du_expl}\n"
+			);
+
+			let dudx = du_expl * testee.dx(ux);
+			assert!(
+				abs_diff_eq!(dudx, 1.0, epsilon = ε),
+				"At x={x}:\ndu/dx * dx/du: {dudx}\n"
+			);
+
+			let d2u_num = (testee.du(x + δ) - testee.du(x - δ)) / (2. * δ);
+			let d2u_expl = testee.d2u(x);
+			assert!(
+				abs_diff_eq!(d2u_expl, d2u_num, epsilon = ε),
+				"At x={x}, d^2u/dx^2:\nnumerical: {d2u_num}\nexplicit:  {d2u_expl}\n"
+			);
+		}
+		for u in itertools_num::linspace(-mul * testee.internal, mul * testee.internal, 100) {
+			let xu = testee.x(u);
+			let uxu = testee.u(xu);
+			assert!(
+				abs_diff_eq!(u, uxu, epsilon = ε),
+				"At u={u}:\nx(u):    {xu}\nu(x(u)): {uxu}\n"
+			);
+
+			let dx_num = (testee.x(u + δ) - testee.x(u - δ)) / (2. * δ);
+			let dx_expl = testee.dx(u);
+			assert!(
+				abs_diff_eq!(dx_expl, dx_num, epsilon = ε),
+				"At u={u}, dx/du:\nnumerical: {dx_num}\nexplicit:  {dx_expl}\n"
+			);
+
+			let dudx = testee.du(xu) * dx_expl;
+			assert!(
+				abs_diff_eq!(dudx, 1.0, epsilon = ε),
+				"At u={u}:\ndu/dx * dx/du: {dudx}\n"
+			);
+		}
+	}
+}

src/ifaces.rs

@@ -0,0 +1,27 @@
+use crate::types::{Hit, Location, Ray};
+
+pub trait Traceable {
+	/// Traces a ray from a given starting point. `ray` is relative to the camera.
+	///
+	/// Returns all objects the ray touched. `Hit::distance` is along the ray; it may be a rough estimate except, when two objects overlap the difference of corresponding `distance`s shall be correct.
+	fn trace(&self, camera: Location, ray: Ray) -> Vec<Hit>;
+}
+
+pub trait OptimizedTraceable: Traceable {
+	type State;
+
+	/// Prepares tracing from a given starting point. `ray` is relative to the camera.
+	fn init(&self, camera: Location, ray: Ray) -> Self::State;
+
+	/// Similar to [`Traceable::trace`] but allows stopping early.
+	fn trace(&self, state: Self::State) -> (Option<Self::State>, Vec<Hit>);
+}
+
+pub struct RayPath {
+	pub points: Vec<Ray>,
+}
+
+pub trait DebugTraceable: Traceable {
+	/// Identical to [`Traceable::trace`], except also returns the ray path.
+	fn trace_dbg(&self, camera: Location, ray: Ray) -> (Vec<Hit>, RayPath);
+}

src/lib.rs

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

src/mathx.rs

@@ -0,0 +1,219 @@
+use glam::{FloatExt, Mat2, Mat3, Vec2, Vec3};
+
+mod bounds {
+	pub trait Pair<T> {}
+
+	impl<T> Pair<T> for (T, T) {}
+}
+
+pub trait FloatExt2<T>: bounds::Pair<T> {
+	fn lerp(self, t: T) -> T;
+	fn inverse_lerp(self, y: T) -> T;
+}
+
+impl<F: FloatExt> FloatExt2<F> for (F, F) {
+	fn lerp(self, t: F) -> F {
+		F::lerp(self.0, self.1, t)
+	}
+	fn inverse_lerp(self, y: F) -> F {
+		F::inverse_lerp(self.0, self.1, y)
+	}
+}
+
+pub trait MatExt {
+	fn orthonormalize(&self) -> Self;
+}
+
+impl MatExt for Mat2 {
+	fn orthonormalize(&self) -> Self {
+		let fx = self.x_axis.normalize();
+		let fy = (self.y_axis - self.y_axis.project_onto_normalized(fx)).normalize();
+		Self::from_cols(fx, fy)
+	}
+}
+
+impl MatExt for Mat3 {
+	fn orthonormalize(&self) -> Self {
+		let fx = self.x_axis.normalize();
+		let fy = (self.y_axis - self.y_axis.project_onto_normalized(fx)).normalize();
+		let fz = (self.z_axis - self.z_axis.project_onto_normalized(fx) - self.z_axis.project_onto_normalized(fy))
+			.normalize();
+		Self::from_cols(fx, fy, fz)
+	}
+}
+
+/// 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
+	}
+}
+
+pub fn make_vec2(f: impl Fn(usize) -> f32) -> Vec2 {
+	Vec2::from_array(std::array::from_fn(|i| f(i)))
+}
+
+pub 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))))
+}
+
+pub fn make_vec3(f: impl Fn(usize) -> f32) -> Vec3 {
+	Vec3::from_array(std::array::from_fn(|i| f(i)))
+}
+
+pub 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))))
+}
+
+#[cfg(test)]
+mod tests {
+	use super::*;
+
+	use approx::assert_abs_diff_eq;
+	use glam::{mat2, mat3, vec2, vec3};
+
+	#[test]
+	fn test_lerp() {
+		assert_eq!((3., 7.).lerp(-0.5), 1.);
+		assert_eq!((3., 7.).lerp(0.0), 3.);
+		assert_eq!((3., 7.).lerp(0.5), 5.);
+		assert_eq!((3., 7.).lerp(1.0), 7.);
+		assert_eq!((3., 7.).lerp(1.5), 9.);
+		assert_eq!((3., 7.).inverse_lerp(1.), -0.5);
+		assert_eq!((3., 7.).inverse_lerp(3.), 0.0);
+		assert_eq!((3., 7.).inverse_lerp(5.), 0.5);
+		assert_eq!((3., 7.).inverse_lerp(7.), 1.0);
+		assert_eq!((3., 7.).inverse_lerp(9.), 1.5);
+	}
+
+	#[test]
+	fn test_orthonormalize_2d() {
+		assert_abs_diff_eq!(
+			mat2(vec2(1., 0.), vec2(0., 1.)).orthonormalize(),
+			mat2(vec2(1., 0.), vec2(0., 1.)),
+		);
+		assert_abs_diff_eq!(
+			mat2(vec2(2., 0.), vec2(3., 5.)).orthonormalize(),
+			mat2(vec2(1., 0.), vec2(0., 1.)),
+		);
+		assert_abs_diff_eq!(
+			mat2(vec2(0., -3.), vec2(5., 1.)).orthonormalize(),
+			mat2(vec2(0., -1.), vec2(1., 0.)),
+		);
+		assert_abs_diff_eq!(
+			mat2(vec2(3., 4.), vec2(5., 1.)).orthonormalize(),
+			mat2(vec2(0.6, 0.8), vec2(0.8, -0.6)),
+		);
+		assert_abs_diff_eq!(
+			mat2(vec2(3., 4.), vec2(1., 5.)).orthonormalize(),
+			mat2(vec2(0.6, 0.8), vec2(-0.8, 0.6)),
+		);
+	}
+
+	#[test]
+	fn test_orthonormalize_3d() {
+		assert_abs_diff_eq!(
+			mat3(vec3(1., 0., 0.), vec3(0., 1., 0.), vec3(0., 0., 1.)).orthonormalize(),
+			mat3(vec3(1., 0., 0.), vec3(0., 1., 0.), vec3(0., 0., 1.)),
+		);
+		assert_abs_diff_eq!(
+			mat3(vec3(2., 0., 0.), vec3(3., 4., 0.), vec3(5., 6., 7.)).orthonormalize(),
+			mat3(vec3(1., 0., 0.), vec3(0., 1., 0.), vec3(0., 0., 1.)),
+		);
+		assert_abs_diff_eq!(
+			mat3(vec3(0., 5., 0.), vec3(0., 7., 6.), vec3(9., 2., 3.)).orthonormalize(),
+			mat3(vec3(0., 1., 0.), vec3(0., 0., 1.), vec3(1., 0., 0.)),
+		);
+	}
+}

src/mesh_loader.rs

@@ -0,0 +1,97 @@
+use glam::{vec2, vec3, Vec2, Vec3};
+use std::io;
+
+#[derive(Copy, Clone, Debug)]
+struct ObjVertex {
+	vertex: usize,
+	normal: usize,
+	// tex_coord: usize,
+}
+
+#[derive(Copy, Clone, Debug)]
+struct ObjFace {
+	vertices: [ObjVertex; 3],
+}
+
+#[derive(Debug)]
+struct ObjMesh {
+	vertices: Vec<Vec3>,
+	normals: Vec<Vec3>,
+	tex_coords: Vec<Vec2>,
+	faces: Vec<ObjFace>,
+}
+
+impl ObjMesh {
+	fn parse_v2(tokens: &[&str]) -> Vec2 {
+		assert_eq!(tokens.len(), 2);
+		let tokens: Vec<_> = tokens.iter().map(|&t| t.parse().unwrap()).collect();
+		vec2(tokens[0], tokens[1])
+	}
+
+	fn parse_v3(tokens: &[&str]) -> Vec3 {
+		assert_eq!(tokens.len(), 3);
+		let tokens: Vec<_> = tokens.iter().map(|&t| t.parse().unwrap()).collect();
+		vec3(tokens[0], tokens[1], tokens[2])
+	}
+
+	fn parse_fv(desc: &&str) -> ObjVertex {
+		let tokens: Vec<_> = desc.split('/').map(|s| s.parse::<usize>().unwrap() - 1).collect();
+		assert_eq!(tokens.len(), 3);
+		ObjVertex {
+			vertex: tokens[0],
+			// tex_coord: tokens[1],
+			normal: tokens[2],
+		}
+	}
+
+	fn parse_f(tokens: &[&str]) -> ObjFace {
+		let vertices: Vec<_> = tokens.iter().map(ObjMesh::parse_fv).collect();
+		ObjFace {
+			vertices: vertices.as_slice().try_into().unwrap(),
+		}
+	}
+
+	fn read(f: &mut impl io::BufRead) -> io::Result<ObjMesh> {
+		let mut result = ObjMesh {
+			vertices: Vec::new(),
+			normals: Vec::new(),
+			tex_coords: Vec::new(),
+			faces: Vec::new(),
+		};
+		loop {
+			let mut line = String::new();
+			if f.read_line(&mut line)? == 0 {
+				break;
+			}
+			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..])),
+				"vt" => result.tex_coords.push(Self::parse_v2(&tokens[1..])),
+				"f" => result.faces.push(Self::parse_f(&tokens[1..])),
+				_ => (),
+			}
+		}
+		Ok(result)
+	}
+
+	fn flatten(&self) -> Vec<Face> {
+		self.faces
+			.iter()
+			.map(|face| Face {
+				vertices: face.vertices.map(|iv| self.vertices[iv.vertex]),
+				normal: self.normals[face.vertices[0].normal],
+			})
+			.collect()
+	}
+}
+
+#[derive(Copy, Clone, Debug)]
+pub struct Face {
+	pub vertices: [Vec3; 3],
+	pub normal: Vec3,
+}
+
+pub fn load_mesh(f: &mut impl io::BufRead) -> io::Result<Vec<Face>> {
+	Ok(ObjMesh::read(f)?.flatten())
+}

src/mesh_tracer.rs

@@ -0,0 +1,36 @@
+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()
+}

src/riemann.rs

@@ -0,0 +1,232 @@
+use crate::mathx::{make_mat3, 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 / 64.0) // there just isn’t enough precision for a small ε.
+	}
+
+	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..3)
+			.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_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::{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/shape/cylinder.rs

@@ -0,0 +1,240 @@
+use glam::{Vec3, Vec3Swizzles as _};
+
+use crate::types::Ray;
+
+pub struct Cylinder {
+	pub center: Vec3,
+	pub semiaxis: Vec3,
+	pub radius: f32,
+}
+
+impl Cylinder {
+	/// Split a vector into a component along the axis and one orthogonal to it.
+	fn split(&self, dir: Vec3) -> (f32, Vec3) {
+		let along = dir.dot(self.semiaxis) / self.semiaxis.length_squared();
+		(along, dir - along * self.semiaxis)
+	}
+
+	fn cap_inout(p: f32, d: f32) -> (f32, f32) {
+		((-d.signum() - p) / d.abs(), (d.signum() - p) / d.abs())
+	}
+
+	pub fn is_inside(&self, pt: Vec3) -> bool {
+		let (along, ortho) = self.split(pt - self.center);
+		along.abs() < 1. && ortho.length_squared() < self.radius.powi(2)
+	}
+
+	fn trace_inout(&self, ray: Ray) -> Option<(f32, f32)> {
+		let (pos_along, pos_ortho) = self.split(ray.pos - self.center);
+		let (dir_along, dir_ortho) = self.split(ray.dir);
+		let (t_cap_in, t_cap_out) = Self::cap_inout(pos_along, dir_along);
+		if dir_ortho.length_squared() < 1e-3 {
+			if pos_ortho.length_squared() >= self.radius.powi(2) {
+				return None;
+			}
+			return Some((t_cap_in, t_cap_out));
+		}
+		let (t_side_in, t_side_out) = solve_quadratic(
+			dir_ortho.length_squared(),
+			pos_ortho.dot(dir_ortho),
+			pos_ortho.length_squared() - self.radius.powi(2),
+		)?;
+		let t_in = f32::max(t_cap_in, t_side_in);
+		let t_out = f32::min(t_cap_out, t_side_out);
+		if t_out <= t_in {
+			return None;
+		}
+		Some((t_in, t_out))
+	}
+
+	pub fn trace_into(&self, ray: Ray) -> Option<f32> {
+		let (t, _) = self.trace_inout(ray)?;
+		if t < 0. {
+			return None;
+		}
+		Some(t)
+	}
+
+	pub fn trace_out_of(&self, ray: Ray) -> Option<f32> {
+		let (_, t) = self
+			.trace_inout(ray)
+			.expect("the ray starts inside so *has* to cross the boundary");
+		Some(t)
+	}
+}
+
+/// Цилиндр с центром в начале координат и осью вдоль оси Y.
+pub struct YCylinder {
+	pub half_length: f32,
+	pub radius: f32,
+}
+
+impl YCylinder {
+	/// Отражает луч, чтобы все координаты направления были положительны (допустимо благодаря симметрии YCylinder).
+	fn flip_ray(ray: Ray) -> Ray {
+		Ray {
+			pos: ray.pos * ray.dir.signum(),
+			dir: ray.dir.abs(),
+		}
+	}
+
+	pub fn is_inside(&self, pt: Vec3) -> bool {
+		let r = f32::hypot(pt.x, pt.z);
+		pt.y.abs() < self.half_length && r < self.radius
+	}
+
+	pub fn trace_into(&self, ray: Ray) -> Option<f32> {
+		let ray = Self::flip_ray(ray);
+
+		// 1. ray.pos.y + t * ray.dir.y = −half_length
+		let t_cap_in = (-self.half_length - ray.pos.y) / ray.dir.y;
+		let t_cap_out = (self.half_length - ray.pos.y) / ray.dir.y;
+
+		// 2. (ray.pos.x + t * ray.dir.x)² + (ray.pos.z + t * ray.dir.z)² = radius²
+		let pos = ray.pos.xz();
+		let dir = ray.dir.xz();
+		if dir.length_squared() < 1e-6 * ray.dir.length_squared() {
+			if pos.length_squared() >= self.radius.powi(2) {
+				return None;
+			}
+			return Some(t_cap_in).filter(|&t| t > 0.);
+		}
+		let (t_side_in, t_side_out) = solve_quadratic(
+			dir.length_squared(),
+			pos.dot(dir),
+			pos.length_squared() - self.radius.powi(2),
+		)?;
+		let t = f32::max(t_cap_in, t_side_in);
+		if t < 0. {
+			return None;
+		}
+		if t >= t_cap_out || t >= t_side_out {
+			return None;
+		}
+		Some(t)
+	}
+
+	pub fn trace_out_of(&self, ray: Ray) -> Option<f32> {
+		let ray = Self::flip_ray(ray);
+		let t_cap_out = (self.half_length - ray.pos.y) / ray.dir.y;
+		let pos = ray.pos.xz();
+		let dir = ray.dir.xz();
+		if dir.length_squared() < 1e-3 {
+			return Some(t_cap_out);
+		}
+		let (_t_side_in, t_side_out) = solve_quadratic(
+			dir.length_squared(),
+			pos.dot(dir),
+			pos.length_squared() - self.radius.powi(2),
+		)
+		.expect("the ray starts inside and is not along the axis so *has* to cross the side");
+		Some(t_side_out)
+	}
+}
+
+fn solve_quadratic(a: f32, half_b: f32, c: f32) -> Option<(f32, f32)> {
+	let base = -half_b / a;
+	let d = base * base - c / a;
+	if d < 0. {
+		None
+	} else {
+		let δ = d.sqrt();
+		Some((base - δ, base + δ))
+	}
+}
+
+#[cfg(test)]
+mod tests {
+	use super::*;
+	use crate::types::ray;
+	use approx::assert_abs_diff_eq;
+	use glam::vec3;
+	use rand::{Rng, SeedableRng};
+
+	#[test]
+	fn test_split() {
+		let mut rng = rand_pcg::Pcg64Mcg::seed_from_u64(17);
+		let cyl = Cylinder {
+			center: vec3(1., 2., 3.),
+			semiaxis: vec3(4., 5., 6.),
+			radius: 7.,
+		};
+		for _ in 0..100 {
+			let dir = vec3(rng.gen(), rng.gen(), rng.gen());
+			let (along, ortho) = cyl.split(dir);
+			assert_abs_diff_eq!(along * cyl.semiaxis + ortho, dir, epsilon = 1e-5);
+			assert_abs_diff_eq!(cyl.semiaxis.dot(ortho), 0., epsilon = 1e-5);
+		}
+	}
+
+	#[test]
+	fn test_cylinder() {
+		assert_eq!(
+			YCylinder::flip_ray(ray(vec3(2., 3., 2.), vec3(4., 5., 4.))),
+			ray(vec3(2., 3., 2.), vec3(4., 5., 4.)),
+		);
+		assert_eq!(
+			YCylinder::flip_ray(ray(vec3(2., 3., 2.), vec3(-4., 5., -4.))),
+			ray(vec3(-2., 3., -2.), vec3(4., 5., 4.)),
+		);
+		assert_eq!(
+			YCylinder::flip_ray(ray(vec3(2., 3., 2.), vec3(4., -5., 4.))),
+			ray(vec3(2., -3., 2.), vec3(4., 5., 4.)),
+		);
+		assert_eq!(
+			YCylinder::flip_ray(ray(vec3(2., 3., 2.), vec3(4., 0., 4.))),
+			ray(vec3(2., 3., 2.), vec3(4., 0., 4.)),
+		);
+
+		let r = YCylinder {
+			half_length: 3.,
+			radius: 2.,
+		};
+
+		assert_eq!(r.trace_into(ray(vec3(3., 4., 3.), vec3(0., -1., 0.))), None);
+		assert_eq!(r.trace_into(ray(vec3(1., 4., 1.), vec3(0., -1., 0.))), Some(1.));
+		assert_eq!(r.trace_into(ray(vec3(3., 3., 3.), vec3(1., 1., 1.))), None);
+		assert_abs_diff_eq!(
+			r.trace_into(ray(vec3(-3., 2., -3.), vec3(1., 0., 1.))).unwrap(),
+			1.5857864
+		);
+		assert_eq!(r.trace_into(ray(vec3(-3., 2., -3.), vec3(-1., 0., -1.))), None);
+		assert_abs_diff_eq!(
+			r.trace_into(ray(vec3(-3., 1., -3.), vec3(2., 2., 2.))).unwrap(),
+			0.7928932
+		);
+		assert_eq!(r.trace_into(ray(vec3(-3., 2.1, -3.), vec3(2., 2., 2.))), None);
+
+		assert_eq!(r.trace_into(ray(vec3(2., 3., 2.), vec3(1., 1., 1.))), None);
+		assert_eq!(r.trace_into(ray(vec3(-2., 3., -2.), vec3(-1., 1., -1.))), None);
+		assert_eq!(
+			r.trace_into(ray(vec3(1.4142135, 3., 1.4142135), vec3(-1., -1., -1.))),
+			Some(0.)
+		);
+		assert_eq!(
+			r.trace_into(ray(vec3(1.4142135, -3., 1.4142135), vec3(-1., 1., -1.))),
+			Some(0.)
+		);
+		assert_eq!(
+			YCylinder {
+				half_length: 300.,
+				radius: 50.
+			}
+			.trace_into(ray(vec3(-125., 375., 0.), vec3(3., -11., 0.) / 1024.)),
+			Some(25600.)
+		);
+
+		assert_abs_diff_eq!(
+			r.trace_out_of(ray(vec3(0., 0., 0.), vec3(1., 1., 1.))).unwrap(),
+			1.4142135
+		);
+		assert_eq!(r.trace_out_of(ray(vec3(0., 0., 0.), vec3(0., 1., 0.))), Some(3.));
+		assert_eq!(r.trace_out_of(ray(vec3(0., 1., 0.), vec3(0., -1., 0.))), Some(4.));
+		assert_eq!(r.trace_out_of(ray(vec3(1., 1., 1.), vec3(0., -1., 0.))), Some(4.));
+		assert_abs_diff_eq!(
+			r.trace_out_of(ray(vec3(1.4142135, 3., 1.4142135), vec3(1., 1., 1.)))
+				.unwrap(),
+			0.
+		);
+	}
+}

src/shape/mod.rs

@@ -0,0 +1,5 @@
+pub mod cylinder;
+pub mod rect;
+
+pub use cylinder::YCylinder;
+pub use rect::Rect;

src/shape/rect.rs

@@ -0,0 +1,90 @@
+use glam::Vec3;
+
+use crate::types::Ray;
+
+pub struct Rect {
+	pub size: Vec3,
+}
+
+impl Rect {
+	/// Отражает луч, чтобы все координаты направления были положительны (допустимо благодаря симметрии Rect).
+	fn flip_ray(ray: Ray) -> Ray {
+		Ray {
+			pos: ray.pos * ray.dir.signum(),
+			dir: ray.dir.abs(),
+		}
+	}
+
+	pub fn is_inside(&self, pt: Vec3) -> bool {
+		pt.abs().cmplt(self.size).all()
+	}
+
+	pub fn trace_into(&self, ray: Ray) -> Option<f32> {
+		let ray = Self::flip_ray(ray);
+		// ray.pos.x + t * ray.dir.x = −size.x
+		let ts = (-self.size - ray.pos) / ray.dir;
+		let t = ts.max_element();
+		let pt = ray.pos + t * ray.dir;
+		if t < 0.0 {
+			return None;
+		}
+		if pt.cmpgt(self.size).any() {
+			return None;
+		}
+		Some(t)
+	}
+
+	pub fn trace_out_of(&self, ray: Ray) -> Option<f32> {
+		let ray = Self::flip_ray(ray);
+		// ray.pos.x + t * ray.dir.x = +size.x
+		let ts = (self.size - ray.pos) / ray.dir;
+		let t = ts.min_element();
+		Some(t)
+	}
+}
+
+#[cfg(test)]
+mod tests {
+	use super::*;
+	use crate::types::ray;
+	use glam::vec3;
+
+	#[test]
+	fn test_rect() {
+		assert_eq!(
+			Rect::flip_ray(ray(vec3(2., 3., 2.), vec3(4., 5., 4.))),
+			ray(vec3(2., 3., 2.), vec3(4., 5., 4.)),
+		);
+		assert_eq!(
+			Rect::flip_ray(ray(vec3(2., 3., 2.), vec3(-4., 5., -4.))),
+			ray(vec3(-2., 3., -2.), vec3(4., 5., 4.)),
+		);
+		assert_eq!(
+			Rect::flip_ray(ray(vec3(2., 3., 2.), vec3(4., -5., 4.))),
+			ray(vec3(2., -3., 2.), vec3(4., 5., 4.)),
+		);
+		assert_eq!(
+			Rect::flip_ray(ray(vec3(2., 3., 2.), vec3(4., 0., 4.))),
+			ray(vec3(2., 3., 2.), vec3(4., 0., 4.)),
+		);
+
+		let r = Rect { size: vec3(2., 3., 2.) };
+
+		assert_eq!(r.trace_into(ray(vec3(3., 3., 3.), vec3(1., 1., 1.))), None);
+		assert_eq!(r.trace_into(ray(vec3(-3., 2., -3.), vec3(1., 0., 1.))), Some(1.));
+		assert_eq!(r.trace_into(ray(vec3(-3., 2., -3.), vec3(-1., 0., -1.))), None);
+		assert_eq!(r.trace_into(ray(vec3(-3., 1., -3.), vec3(2., 2., 2.))), Some(0.5));
+		assert_eq!(r.trace_into(ray(vec3(-3., 2.1, -3.), vec3(2., 2., 2.))), None);
+
+		assert_eq!(r.trace_into(ray(vec3(2., 3., 2.), vec3(1., 1., 1.))), None);
+		assert_eq!(r.trace_into(ray(vec3(-2., 3., -2.), vec3(-1., 1., -1.))), None);
+		assert_eq!(r.trace_into(ray(vec3(2., 3., 2.), vec3(-1., -1., -1.))), Some(0.));
+		assert_eq!(r.trace_into(ray(vec3(2., -3., 2.), vec3(-1., 1., -1.))), Some(0.));
+
+		assert_eq!(r.trace_out_of(ray(vec3(0., 0., 0.), vec3(1., 1., 1.))), Some(2.));
+		assert_eq!(r.trace_out_of(ray(vec3(0., 0., 0.), vec3(0., 1., 0.))), Some(3.));
+		assert_eq!(r.trace_out_of(ray(vec3(0., 1., 0.), vec3(0., -1., 0.))), Some(4.));
+		assert_eq!(r.trace_out_of(ray(vec3(1., 1., 1.), vec3(0., -1., 0.))), Some(4.));
+		assert_eq!(r.trace_out_of(ray(vec3(2., 3., 2.), vec3(1., 1., 1.))), Some(0.));
+	}
+}

src/tube/coords.rs

@@ -0,0 +1,468 @@
+use glam::{vec3, Mat3, Vec3};
+
+use crate::riemann::Metric;
+use crate::shape::YCylinder;
+use crate::types::{Location, Ray};
+
+use super::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> {
+		YCylinder {
+			radius: self.0.inner_radius,
+			half_length: 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<Vec3> for OuterCS {
+	fn flat_to_global(&self, pos: Vec3) -> Vec3 {
+		let inner = YCylinder {
+			radius: self.0.inner_radius + 1.0,
+			half_length: self.0.external_halflength,
+		};
+		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 = YCylinder {
+			radius: self.0.inner_radius + 1.0,
+			half_length: self.0.external_halflength,
+		};
+		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> {
+		YCylinder {
+			radius: self.0.outer_radius,
+			half_length: self.0.external_halflength,
+		}
+		.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
+		// NOTE it can’t test −30..30 as that area intersects the boundary
+		test_flat_region(
+			&mapper,
+			(vec3(-20.0, -300.0, -20.0), vec3(20.0, -1.0, 20.0)),
+			(vec3(-20.0, -300.0, -20.0), vec3(20.0, -200.20016, 20.0)),
+		);
+		test_flat_region(
+			&mapper,
+			(vec3(-20.0, 1.0, -20.0), vec3(20.0, 300.0, 20.0)),
+			(vec3(-20.0, 200.20016, -20.0), vec3(20.0, 300.0, 20.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(-21., 21., 20) {
+			for y in linspace(1., 299., 20) {
+				for z in linspace(-21., 21., 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

@@ -0,0 +1,239 @@
+use glam::{f32, vec3, Mat3, Vec3};
+
+use crate::fns::{self, Limiter};
+use crate::mathx::Decomp3;
+use crate::riemann::{Metric, Tens3};
+
+#[derive(Copy, Clone, Debug)]
+pub struct Tube {
+	pub outer_radius: f32,
+	pub inner_radius: f32,
+	pub external_halflength: f32,
+	pub internal_halflength: f32,
+}
+
+impl Tube {
+	fn fx(&self) -> impl Limiter {
+		fns::SmootherstepLimiter {
+			min: self.inner_radius,
+			max: self.outer_radius,
+		}
+	}
+	fn fy(&self) -> fns::QuadraticAccelerator {
+		fns::QuadraticAccelerator {
+			internal: self.internal_halflength,
+			external: self.external_halflength,
+		}
+	}
+
+	pub fn y(&self, v: f32) -> f32 {
+		self.fy().x(v)
+	}
+	pub fn v(&self, y: f32) -> f32 {
+		self.fy().u(y)
+	}
+	pub fn dy(&self, v: f32) -> f32 {
+		self.fy().dx(v)
+	}
+	pub fn dv(&self, y: f32) -> f32 {
+		self.fy().du(y)
+	}
+}
+
+impl Metric for Tube {
+	fn sqrt_at(&self, pos: Vec3) -> Decomp3 {
+		let r = f32::hypot(pos.x, pos.z);
+		let sr = self.fx().value(r);
+		let sy = self.fy().du(pos.y);
+		let s = sr + sy - sr * sy;
+		assert!(sr.is_finite());
+		assert!(sy.is_finite());
+		assert!(sy > 0.0);
+		Decomp3 {
+			ortho: Mat3::IDENTITY,
+			diag: vec3(1.0, s, 1.0),
+		}
+	}
+
+	fn part_derivs_at(&self, pos: Vec3) -> Tens3 {
+		let r = f32::hypot(pos.x, pos.z);
+		let sr = self.fx().value(r);
+		let sy = self.fy().du(pos.y);
+		let s = sr + sy - sr * sy;
+		let dsr_dr = self.fx().derivative(r);
+		let dsy_dy = self.fy().d2u(pos.y);
+		let ds2_dr = 2.0 * s * (1.0 - sy) * dsr_dr;
+		let ds2_dy = 2.0 * s * (1.0 - sr) * dsy_dy;
+		let (ds2_dx, ds2_dz) = if ds2_dr.abs() < 1e-5 {
+			(0., 0.)
+		} else {
+			(ds2_dr * pos.x / r, ds2_dr * pos.z / r)
+		};
+		[
+			Mat3::from_cols_array(&[0., 0., 0., 0., ds2_dx, 0., 0., 0., 0.]),
+			Mat3::from_cols_array(&[0., 0., 0., 0., ds2_dy, 0., 0., 0., 0.]),
+			Mat3::from_cols_array(&[0., 0., 0., 0., ds2_dz, 0., 0., 0., 0.]),
+		]
+	}
+}
+
+#[cfg(test)]
+mod test {
+	use approx::assert_abs_diff_eq;
+	use glam::{vec3, Vec3};
+	use itertools_num::linspace;
+
+	use crate::mathx::Decomp3;
+	use crate::riemann::Metric;
+	use crate::tube::Space;
+	use crate::types::Ray;
+
+	use super::Tube;
+
+	#[test]
+	fn test_tube_metric_derivs() {
+		struct Approx(Tube);
+		impl Metric for Approx {
+			fn sqrt_at(&self, pos: Vec3) -> Decomp3 {
+				self.0.sqrt_at(pos)
+			}
+		}
+		let testee = Tube {
+			inner_radius: 30.0,
+			outer_radius: 50.0,
+			internal_halflength: 100.0,
+			external_halflength: 300.0,
+		};
+		let approx = Approx(testee);
+		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, 20) {
+			for y in itertools_num::linspace(-mul * testee.external_halflength, mul * testee.external_halflength, 20) {
+				for z in itertools_num::linspace(-mul * testee.outer_radius, mul * testee.outer_radius, 20) {
+					let pos = vec3(x, y, z);
+					let computed = testee.part_derivs_at(pos);
+					let reference = approx.part_derivs_at(pos);
+					let eq = (0..3).all(|coord| computed[coord].abs_diff_eq(reference[coord], epsilon));
+					assert!(
+						eq,
+						"Bad derivative computation at {pos}:\n  explicit:  {computed:?}\n  numerical: {reference:?}\n"
+					);
+				}
+			}
+		}
+	}
+
+	#[test]
+	fn test_accelerator() {
+		let space = Space {
+			tube: Tube {
+				inner_radius: 30.0,
+				outer_radius: 50.0,
+				internal_halflength: 100.0,
+				external_halflength: 300.0,
+			},
+			objs: vec![],
+		};
+		let ε = 1e-3;
+		let off = 10.0;
+		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 = vec3(ax, -(space.tube.external_halflength + off), 0.);
+				let b = vec3(bx, space.tube.external_halflength + off, 0.);
+				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);
+				assert_abs_diff_eq!(traced.pos.y, b.y, epsilon = 1.0e1);
+				assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon = 1.0e-3);
+				assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon = 1.0e-2);
+			}
+		}
+	}
+
+	#[test]
+	#[ignore]
+	fn test_accelerator_slow() {
+		let space = Space {
+			tube: Tube {
+				inner_radius: 30.0,
+				outer_radius: 50.0,
+				internal_halflength: 100.0,
+				external_halflength: 300.0,
+			},
+			objs: vec![],
+		};
+		let ε = 1e-3;
+		let off = 10.0;
+		let steps = 4096;
+		for ax in linspace(-space.tube.inner_radius + ε, space.tube.inner_radius - ε, 20) {
+			for bx in linspace(-space.tube.inner_radius + ε, space.tube.inner_radius - ε, 20) {
+				let a = vec3(ax, -(space.tube.external_halflength + off), 0.);
+				let b = vec3(bx, space.tube.external_halflength + off, 0.);
+				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);
+				assert_abs_diff_eq!(traced.pos.y, b.y, epsilon = 1.0e0);
+				assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon = 1.0e-3);
+				assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon = 1.0e-3);
+			}
+		}
+	}
+
+	#[test]
+	fn test_accelerator_inner_edge() {
+		let space = Space {
+			tube: Tube {
+				inner_radius: 30.0,
+				outer_radius: 50.0,
+				internal_halflength: 100.0,
+				external_halflength: 300.0,
+			},
+			objs: vec![],
+		};
+		let ε = 1e-3;
+		let off = 10.0;
+		let steps = 10000;
+		for x in [space.tube.inner_radius - ε, space.tube.inner_radius + ε] {
+			let a = vec3(x, -(space.tube.external_halflength + off), 0.);
+			let b = vec3(x, space.tube.external_halflength + off, 0.);
+			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);
+			assert_abs_diff_eq!(traced.pos.y, b.y, epsilon = 1.0e0);
+			assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon = 1.0e-2);
+			assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon = 1.0e-2);
+		}
+	}
+
+	#[test]
+	fn test_accelerator_outer_edge() {
+		let space = Space {
+			tube: Tube {
+				inner_radius: 30.0,
+				outer_radius: 50.0,
+				internal_halflength: 100.0,
+				external_halflength: 300.0,
+			},
+			objs: vec![],
+		};
+		let ε = 1e-3;
+		let off = 10.0;
+		let steps = 4096;
+		for x in [space.tube.outer_radius + ε, space.tube.outer_radius - ε] {
+			let a = vec3(x, -(space.tube.external_halflength + off), 0.);
+			let b = vec3(x, space.tube.external_halflength + off, 0.);
+			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);
+			assert_abs_diff_eq!(traced.pos.y, b.y, epsilon = 1.0e0);
+			assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon = 1.0e-2);
+			assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon = 1.0e-2);
+		}
+	}
+}

src/tube/mod.rs

@@ -0,0 +1,228 @@
+use glam::{f32, Mat3, Vec3};
+
+use crate::ifaces::{DebugTraceable, RayPath, Traceable};
+use coords::{FlatCoordinateSystem, InnerCS, OuterCS};
+use metric::Tube;
+use Subspace::{Boundary, Inner, Outer};
+
+use crate::riemann::Metric;
+use crate::tube::coords::FlatRegion;
+use crate::types::{FlatTraceResult, Hit, Location, Object, Ray};
+use crate::{riemann, DT};
+
+mod coords;
+pub mod metric;
+
+pub struct Space {
+	pub tube: Tube,
+	pub objs: Vec<Object>,
+}
+
+#[derive(PartialEq, Eq, Debug)]
+pub enum Subspace {
+	Outer,
+	Boundary,
+	Inner,
+}
+
+impl Space {
+	fn which_subspace(&self, pt: Vec3) -> Subspace {
+		if pt.y.abs() > self.tube.external_halflength {
+			return Outer;
+		}
+		let r = f32::hypot(pt.x, pt.z);
+		if r > self.tube.outer_radius {
+			Outer
+		} else if r > self.tube.inner_radius {
+			Boundary
+		} else {
+			Inner
+		}
+	}
+
+	/// Выполняет один шаг трассировки. Работает в любой части пространства, но вне Boundary доступны более эффективные методы.
+	/// ray задаётся в основной СК.
+	pub fn trace_step(&self, ray: Ray) -> Ray {
+		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 }
+	}
+
+	/// Выполняет один шаг перемещения. Работает в любой части пространства.
+	/// off задаётся в локальной СК. Рекомендуется считать небольшими шагами.
+	pub fn move_step(&self, loc: Location, off: Vec3) -> Location {
+		let corr = Mat3::IDENTITY - riemann::contract(riemann::krist(&self.tube, loc.pos), loc.rot * off);
+		let p = loc.pos + corr * loc.rot * off;
+		Location {
+			pos: p,
+			rot: corr * loc.rot,
+		}
+	}
+
+	pub fn trace_iter(&self, ray: Ray) -> impl Iterator<Item = Ray> + '_ {
+		std::iter::successors(Some(ray), |&ray| Some(self.trace_step(ray)))
+	}
+
+	fn trace_inner(&self, ray: Ray) -> FlatTraceResult {
+		assert_eq!(self.which_subspace(ray.pos), Inner);
+		self.trace_flat(InnerCS(self.tube), ray)
+	}
+
+	fn trace_outer(&self, ray: Ray) -> FlatTraceResult {
+		assert_eq!(self.which_subspace(ray.pos), Outer);
+		self.trace_flat(OuterCS(self.tube), ray)
+	}
+
+	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),
+			Boundary => None,
+		}
+	}
+
+	fn trace_flat(&self, cs: impl FlatRegion, ray: Ray) -> FlatTraceResult {
+		let ray = cs.global_to_flat(ray);
+		let dist = cs.distance_to_boundary(ray);
+		let objs = self.list_objects(|loc| cs.global_to_flat(loc));
+		FlatTraceResult {
+			end: dist.map(|dist| cs.flat_to_global(ray.forward(dist))),
+			objects: Self::hit_objects(objs.as_slice(), ray, dist, |pos| cs.flat_to_global(pos)),
+		}
+	}
+
+	fn list_objects(&self, tfm: impl Fn(Location) -> Location) -> Vec<Object> {
+		self.objs
+			.iter()
+			.map(|&Object { id, loc, r }| Object { id, loc: tfm(loc), r })
+			.collect()
+	}
+
+	fn hit_objects(objs: &[Object], ray: Ray, limit: Option<f32>, globalize: impl Fn(Vec3) -> Vec3) -> 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) - 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))
+				} else {
+					None
+				}
+			})
+			.filter(|&(_, t)| t >= 0.0 && t < limit)
+			.map(|(obj, t)| {
+				let pos = ray.forward(t).pos;
+				let rel = obj.loc.rot.inverse()
+					* Ray {
+						pos: pos - obj.loc.pos,
+						dir: ray.dir,
+					};
+				Hit {
+					id: obj.id,
+					distance: t,
+					pos: globalize(pos),
+					rel,
+				}
+			})
+			.collect()
+	}
+
+	pub fn line(&self, a: Vec3, b: Vec3, step: f32) -> Vec<Ray> {
+		match self.which_subspace(a) {
+			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(Ray {
+							pos: a.lerp(b, k as f32 / n as f32),
+							dir,
+						})
+					})
+					.collect()
+			}
+			Boundary => panic!("Can't draw a line here!"),
+		}
+	}
+
+	fn camera_ray_to_abs(&self, camera: Location, ray: Ray) -> Ray {
+		let pos = camera.pos;
+		let dir = camera.rot * ray.dir;
+		// TODO account for ray.pos
+		let dir = DT * self.tube.normalize_vec_at(pos, dir);
+		Ray { pos, dir }
+	}
+}
+
+/// Like [`std::iter::successors`] but with an upper limit on iteration count.
+///
+/// # Panics
+///
+/// Panics if the sequence doesn’t terminate in `max_iters` calls of `succ`.
+fn iterate_with_limit<T>(max_iters: usize, init: T, mut succ: impl FnMut(T) -> Option<T>) {
+	let mut state = init;
+	for _ in 0..max_iters {
+		match succ(state) {
+			Some(next) => state = next,
+			None => return,
+		}
+	}
+	panic!("iteration limit exceeded");
+}
+
+impl Traceable for Space {
+	fn trace(&self, camera: Location, ray: Ray) -> Vec<Hit> {
+		let ray = self.camera_ray_to_abs(camera, ray);
+		let mut hits = vec![];
+		iterate_with_limit(100, ray, |ray| {
+			let ret = self
+				.trace_iter(ray)
+				.skip(1)
+				.find_map(|ray| self.obj_hitter(ray.pos).map(|hitter| hitter(self, ray)))
+				.expect("Space::trace_iter does not terminate");
+			hits.extend(ret.objects); // TODO fix distance
+			ret.end
+		});
+		hits
+	}
+}
+
+impl DebugTraceable for Space {
+	fn trace_dbg(&self, camera: Location, ray: Ray) -> (Vec<Hit>, RayPath) {
+		let mut points = vec![];
+		let mut hits = vec![];
+		let mut ray = self.camera_ray_to_abs(camera, ray);
+
+		let trace_to_flat = |points: &mut Vec<Ray>, ray| {
+			for ray in self.trace_iter(ray).skip(1) {
+				points.push(ray);
+				if let Some(hitter) = self.obj_hitter(ray.pos) {
+					return (ray, hitter(self, ray));
+				}
+			}
+			unreachable!("Space::trace_iter terminated!")
+		};
+
+		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 (hits, RayPath { points });
+			};
+			points.extend(self.line(ray_into_flat.pos, ray_outta_flat.pos, 1.0));
+			ray = ray_outta_flat;
+		}
+		panic!("tracing didn't terminate");
+	}
+}

src/types.rs

@@ -1,31 +1,41 @@
-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,
+}
+
+pub fn ray(pos: Vec3, dir: Vec3) -> Ray {
+	Ray { pos, dir }
 }
 
 impl Ray {
 	pub fn forward(&self, dist: f32) -> Ray {
-		Ray { pos: self.pos + self.dir * dist, dir: self.dir }
+		Ray {
+			pos: self.pos + self.dir * dist,
+			dir: self.dir,
+		}
 	}
 }
 
-impl std::ops::Mul<Ray> for Mat2 {
+impl std::ops::Mul<Ray> for Mat3 {
 	type Output = Ray;
 
 	fn mul(self, rhs: Ray) -> Self::Output {
-		Ray { pos: self * rhs.pos, dir: self * rhs.dir }
+		Ray {
+			pos: self * rhs.pos,
+			dir: self * rhs.dir,
+		}
 	}
 }
 
-#[derive(Copy, Clone, Debug)]
+#[derive(Copy, Clone, Debug, PartialEq)]
 pub struct Location {
 	/// Положение в основной СК
-	pub pos: Vec2,
+	pub pos: Vec3,
 	/// Преобразование вектора из локальной ортонормированной в основную СК
-	pub rot: Mat2,
+	pub rot: Mat3,
 }
 
 #[derive(Copy, Clone, Debug)]
@@ -39,8 +49,8 @@ pub struct Object {
 pub struct Hit {
 	pub distance: f32,
 	pub id: i32,
-	pub pos: Vec2, // положение в основной СК
+	pub pos: Vec3, // положение в основной СК
-	pub rel: Ray, // в локальной ортонормированной СК объекта
+	pub rel: Ray,  // в локальной ортонормированной СК объекта
 }
 
 #[derive(Clone, Debug)]

src/utils.rs

@@ -0,0 +1,68 @@
+//! 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 = ε);
+	}
+}