Cargo.toml

@@ -3,8 +3,7 @@ name = "refraction"
 version = "0.1.0"
 edition = "2021"
 
-[lints.rust]
-mixed_script_confusables = "allow"
+# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
 
 [profile.dev]
 panic = 'abort'
@@ -22,11 +21,9 @@ show-image = "0.14.0"
 flo_draw = "0.3.1"
 flo_canvas = "0.3.1"
 itertools-num = "0.1.3"
-winit = "0.29"
+glium = "0.35.0"
+winit = "0.30.5"
 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"

rustfmt.toml

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

src/bin/flat/main.rs

@@ -52,12 +52,20 @@ pub fn main() {
 			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.),
+				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.),
+				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.)),
 			);
 
@@ -126,7 +134,11 @@ pub fn main() {
 							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
+							space
+								.trace_iter(Ray { pos, dir })
+								.nth(n as usize)
+								.unwrap()
+								.pos
 						}),
 				);
 			}
@@ -134,7 +146,6 @@ pub fn main() {
 	});
 }
 
-#[allow(dead_code)]
 fn draw_cross(gc: &mut Vec<Draw>, pos: Vec2, r: f32) {
 	gc.move_to(pos.x - r, pos.y - r);
 	gc.line_to(pos.x + r, pos.y + r);
@@ -155,10 +166,13 @@ fn draw_ray_2(gc: &mut Vec<Draw>, space: &Space, camera: Location, dir: Vec3) {
 	gc.new_path();
 	gc.move_to(pos.x, pos.y);
 	for pt in &path.points[1..] {
-		gc.line_to(pt.pos.x, pt.pos.y);
+		gc.line_to(pt.x, pt.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;
+	let end_pos = *path
+		.points
+		.last()
+		.expect("the starting point is always in the path");
+	let dir_pos = end_pos + 1000.0 * path.end_dir;
 	gc.line_to(dir_pos.x, dir_pos.y);
 	gc.stroke();
 }
@@ -177,7 +191,11 @@ fn draw_track(gc: &mut Vec<Draw>, space: &Space, start: Vec2, dir: Vec2) {
 	// let v = space.tube.normalize(start, dir);
 	let mut loc = Location {
 		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.)),
+		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| {
@@ -197,8 +215,8 @@ 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;
-		for _ in 0..n {
+		let N = (STEP / DT).floor() as i32;
+		for _ in 0..N {
 			loc = space.move_step(loc, v * DT);
 		}
 		draw(&loc);

src/bin/mesh/main.rs

@@ -1,14 +1,12 @@
 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;
@@ -81,7 +79,9 @@ fn render(mesh: &Mesh, camera: impl Fn(Vec2) -> (Vec3, Vec3)) -> Image {
 			} else {
 				bkg
 			};
-			let color = (color * 255.0).as_ivec3().clamp(IVec3::splat(0), IVec3::splat(255));
+			let color = (color * 255.0)
+				.as_ivec3()
+				.clamp(IVec3::splat(0), IVec3::splat(255));
 			img.put_pixel(x, y, Color(color.x as u8, color.y as u8, color.z as u8));
 		}
 	}
@@ -109,10 +109,6 @@ fn test_projs() {
 	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();
@@ -125,21 +121,15 @@ fn main() -> Result<(), Box<dyn Error>> {
 		let mut f = BufReader::new(f);
 		load_mesh(&mut f)?
 	};
+	let proj = persp;
 	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_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);

src/bin/wireframe/camera.rs

@@ -1,123 +0,0 @@
-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

@@ -1,168 +0,0 @@
-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

@@ -1,38 +1,73 @@
 use std::time::Instant;
 
-use glam::{uvec2, vec3, Vec3};
+use glam::{mat4, vec2, vec3, vec4, Mat4, Vec3};
+use glium::{
+	backend::{glutin::SimpleWindowBuilder, Facade},
+	glutin::config::ConfigTemplateBuilder,
+	implement_vertex,
+	index::PrimitiveType,
+	uniform,
+	winit::event::{Event, WindowEvent},
+	DrawParameters, Program, Surface, VertexBuffer,
+};
 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)
+#[derive(Copy, Clone)]
+struct Vertex {
+	position: [f32; 3],
+}
+implement_vertex!(Vertex, position);
+
+struct Wireframe {
+	color: Vec3,
+	mode: PrimitiveType,
+	data: VertexBuffer<Vertex>,
+}
+
+fn prepare_scene(display: &impl Facade) -> Vec<Wireframe> {
+	scene::build()
+		.into_iter()
+		.map(|line| {
+			let color = line.color;
+			let mode;
+			let data: Vec<Vertex>;
+			match line.line {
+				scene::Line::Lines(_) => todo!(),
+				scene::Line::Strip(pts) => {
+					mode = PrimitiveType::LineStrip;
+					data = pts
+						.into_iter()
+						.map(|p| Vertex {
+							position: p.to_array(),
+						})
+						.collect();
+				}
+				scene::Line::Loop(pts) => {
+					mode = PrimitiveType::LineLoop;
+					data = pts
+						.into_iter()
+						.map(|p| Vertex {
+							position: p.to_array(),
+						})
+						.collect();
+				}
+			};
+			let data = VertexBuffer::new(display, &data).unwrap();
+			Wireframe { color, mode, data }
+		})
+		.collect()
 }
 
-#[cfg(any())]
 mod camctl {
 	use glam::{vec3, Mat4, Quat, Vec3};
 
@@ -57,48 +92,15 @@ mod camctl {
 		}
 
 		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));
+			self.rotate_rel_quat(Quat::from_euler(glam::EulerRot::XYZ, ypr.x, ypr.y, ypr.z));
 		}
 
-		fn rotate_rel_quat(&mut self, rot: Quat) {
+		pub 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};
@@ -140,278 +142,162 @@ static KEYS_MOVE: &'static [(PhysicalKey, Vec3)] = &[
 	(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::Numpad9), vec3(1., 0., 0.)),
+	(PhysicalKey::Code(KeyCode::Numpad7), vec3(-1., 0., 0.)),
 	(PhysicalKey::Code(KeyCode::Numpad5), vec3(0., 1., 0.)),
 	(PhysicalKey::Code(KeyCode::Numpad8), vec3(0., -1., 0.)),
-	(PhysicalKey::Code(KeyCode::Numpad9), vec3(0., 0., 1.)),
-	(PhysicalKey::Code(KeyCode::Numpad7), vec3(0., 0., -1.)),
+	(PhysicalKey::Code(KeyCode::Numpad4), vec3(0., 0., 1.)),
+	(PhysicalKey::Code(KeyCode::Numpad6), 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()
+fn main() {
+	let event_loop = EventLoop::builder().build().unwrap();
+	let cfg = ConfigTemplateBuilder::new().with_multisampling(8);
+	let (window, display) = SimpleWindowBuilder::new()
+		.with_config_template_builder(cfg)
 		.with_title("Refraction: Wireframe")
-		.build(&event_loop)
-		.unwrap();
+		.build(&event_loop);
 
-	// 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;
+	let vs_src = include_str!("ray.v.glsl");
+	let fs_src = include_str!("ray.f.glsl");
+	let program = Program::from_source(&display, vs_src, fs_src, None).unwrap();
 
+	let scene = prepare_scene(&display);
+
+	let mut kbd = keyctl::Keyboard::new();
+	let mut cam = camctl::CameraLocation::new();
+	let mut t1 = Instant::now();
+
+	#[allow(deprecated)]
 	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();
+		.run(move |ev, window_target| match ev {
+			Event::WindowEvent { event, .. } => match event {
+				WindowEvent::RedrawRequested => {
+					let dt = {
+						let t2 = Instant::now();
+						let dt = t2 - t1;
+						t1 = t2;
+						dt.as_secs_f32()
+					};
+					cam.move_rel(100. * dt * kbd.control(&KEYS_MOVE));
+					cam.rotate_rel_ypr(2. * dt * kbd.control(&KEYS_ROTATE));
 
-							if !surface_configured {
-								return;
-							}
+					let size = window.inner_size();
+					let size = vec2(size.width as f32, size.height as f32).normalize()
+						* std::f32::consts::SQRT_2;
+					let proj = make_proj_matrix(vec3(size.x, size.y, 2.), (1., 4096.));
 
-							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();
-								}
+					let my_to_gl = mat4(
+						vec4(0., 0., 1., 0.),
+						vec4(-1., 0., 0., 0.),
+						vec4(0., 1., 0., 0.),
+						vec4(0., 0., 0., 1.),
+					);
+					let view = my_to_gl * cam.view_mtx();
 
-								// This happens when the a frame takes too long to present
-								Err(wgpu::SurfaceError::Timeout) => {
-									eprintln!("Surface timeout")
-								}
-							}
-						}
-						_ => {}
+					let mut target = display.draw();
+					target.clear_color(0.0, 0.0, 0.0, 0.0);
+
+					let mvp = proj * view;
+					let params = DrawParameters {
+						blend: glium::Blend {
+							color: glium::BlendingFunction::Addition {
+								source: glium::LinearBlendingFactor::One,
+								destination: glium::LinearBlendingFactor::OneMinusSourceAlpha,
+							},
+							alpha: glium::BlendingFunction::Addition {
+								source: glium::LinearBlendingFactor::One,
+								destination: glium::LinearBlendingFactor::OneMinusSourceAlpha,
+							},
+							constant_value: (0., 0., 0., 0.),
+						},
+						line_width: Some(3.),
+						smooth: Some(glium::Smooth::Nicest),
+						..Default::default()
+					};
+					for mesh in &scene {
+						let uniforms = uniform! {
+							mvp: mvp.to_cols_array_2d(),
+							color: mesh.color.to_array(),
+						};
+						let indices = glium::index::NoIndices(mesh.mode);
+						target
+							.draw(&mesh.data, &indices, &program, &uniforms, &params)
+							.unwrap();
 					}
+					target.finish().unwrap();
 				}
-				_ => {}
+
+				WindowEvent::KeyboardInput {
+					device_id: _,
+					event,
+					is_synthetic: _,
+				} => {
+					kbd.set_key_state(event.physical_key, event.state);
+				}
+
+				WindowEvent::Resized(window_size) => {
+					display.resize(window_size.into());
+				}
+
+				WindowEvent::CloseRequested => {
+					window_target.exit();
+				}
+
+				_ => (),
+			},
+
+			Event::AboutToWait => {
+				window.request_redraw();
 			}
+			_ => (),
 		})
 		.unwrap();
 }
 
-fn main() {
-	pollster::block_on(run());
+/// Make a projection matrix, assuming input coordinates are (right, up, forward).
+///
+/// `corner` is a vector that will be mapped to (x=1, y=1) after the perspective division.
+/// `zrange` is the Z range that will be mapped to z∈[-1, 1]. It has no other effect. Both ends have to be positive though.
+fn make_proj_matrix(corner: Vec3, zrange: (f32, f32)) -> Mat4 {
+	let scale = 1.0 / corner;
+	let zspan = zrange.1 - zrange.0;
+	mat4(
+		scale.x * vec4(1., 0., 0., 0.),
+		scale.y * vec4(0., 1., 0., 0.),
+		scale.z * vec4(0., 0., (zrange.0 + zrange.1) / zspan, 1.),
+		scale.z * vec4(0., 0., -2. * zrange.0 * zrange.1 / zspan, 0.),
+	)
 }
 
-mod fps {
-	use std::time::{Duration, Instant};
+#[cfg(test)]
+mod tests {
+	use super::*;
+	use approx::assert_abs_diff_eq;
+	use glam::vec3;
 
-	pub struct Counter {
-		fps: f32,
-		t1: Instant,
-		frames: u32,
-	}
+	#[test]
+	fn test_proj_matrix() {
+		let m = make_proj_matrix(vec3(2., 3., 4.), (0.5, 20.0));
 
-	impl Counter {
-		pub fn new() -> Self {
-			Self {
-				fps: 0.,
-				t1: Instant::now(),
-				frames: 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}");
 
-		pub fn get(&self) -> f32 {
-			self.fps
-		}
+		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);
 
-		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,
-				}
-			}
-		}
+		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/mesh.wgsl

@@ -1,63 +0,0 @@
-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

@@ -1,160 +0,0 @@
-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.f.glsl

@@ -0,0 +1,10 @@
+#version 330
+
+in vec3 vertex_color;
+
+out vec4 color;
+
+void main() {
+	float opacity = pow(0.5 - 0.5 * gl_FragCoord.z, 0.25);
+	color = opacity * vec4(vertex_color, 1.0);
+}

src/bin/wireframe/ray.v.glsl

@@ -0,0 +1,13 @@
+#version 330
+
+uniform mat4 mvp;
+uniform vec3 color;
+
+in vec3 position;
+
+out vec3 vertex_color;
+
+void main() {
+	vertex_color = color;
+	gl_Position = mvp * vec4(position, 1.0);
+}

src/bin/wireframe/ray.wgsl

@@ -1,61 +0,0 @@
-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

@@ -1,5 +1,5 @@
 use glam::*;
-use itertools::chain;
+use itertools::{chain, iproduct};
 
 use refraction::ifaces::{DebugTraceable, Traceable};
 use refraction::tube::metric::Tube;
@@ -8,62 +8,42 @@ use refraction::types::{Location, Object, Ray};
 use refraction::utils::put_object;
 
 pub enum Line {
-	Strip(Vec<Ray>),
-	Loop(Vec<Ray>),
+	Lines(Vec<(Vec3, Vec3)>),
+	Strip(Vec<Vec3>),
+	Loop(Vec<Vec3>),
 }
 
 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>,
+	pub line: Line,
 }
 
 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
-			}
-		},
-	}))
+	onto.extend(lines.into_iter().map(move |line| FancyLine { color, line }))
 }
 
-fn draw_line(a: Vec3, b: Vec3) -> Line {
-	let dir = (b - a).normalize();
-	Line::Strip(vec![Ray { pos: a, dir }, Ray { pos: b, dir }])
+fn draw_rect(center: Vec3, u: Vec3, v: Vec3) -> Line {
+	Line::Loop(vec![
+		center - u - v,
+		center + u - v,
+		center + u + v,
+		center - u + v,
+	])
 }
 
-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()
+fn draw_ellipse(center: Vec3, u: Vec3, v: Vec3) -> Line {
+	let segments = 47;
+	let step = 2. * std::f32::consts::PI / segments as f32;
+	Line::Loop(
+		(0..segments)
+			.map(|k| k as f32 * step)
+			.map(Vec2::from_angle)
+			.map(|d| center + d.x * u + d.y * v)
+			.collect(),
+	)
 }
 
-pub fn build() -> (Vec<FancyMesh>, Vec<FancyLine>) {
+pub fn build() -> Vec<FancyLine> {
 	let tube = Tube {
 		inner_radius: 30.0,
 		outer_radius: 50.0,
@@ -88,123 +68,76 @@ pub fn build() -> (Vec<FancyMesh>, Vec<FancyLine>) {
 	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.),
+		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.),
+		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)
+	paint(&mut gc, vec3(0.8, 0.8, 0.8), tube.render());
+	paint(&mut gc, vec3(0.0, 0.8, 1.0), draw_fan_2(&space, cam3, 1.0));
+	paint(&mut gc, vec3(0.5, 1.0, 0.0), draw_fan_2(&space, cam2, 1.0));
+	paint(&mut gc, vec3(1.0, 0.5, 0.0), draw_fan_2(&space, cam1, 1.0));
+	gc
 }
 
 fn draw_ray_2(gc: &mut Vec<Line>, space: &Space, camera: Location, dir: Vec3) {
 	let pos = vec3(0., 0., 0.);
 	let (hits, path) = space.trace_dbg(camera, Ray { pos, dir });
-	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 hits2 = space.trace(camera, Ray { pos, dir });
+	for (a, b) in hits.into_iter().zip(hits2.into_iter()) {
+		assert_eq!(a.id, b.id);
+		assert_eq!(a.pos, b.pos);
+		assert_eq!(a.rel, b.rel);
 	}
 
 	let mut pts = path.points;
-	let end_pos = *pts.last().expect("the starting point is always in the path");
-	pts.extend(itertools::iterate(end_pos, |r| r.forward(100.0)).take(1000));
+	let end_pos = *pts
+		.last()
+		.expect("the starting point is always in the path");
+	let dir_pos = end_pos + 10000.0 * path.end_dir;
+	pts.push(dir_pos);
 	gc.push(Line::Strip(pts));
 }
 
-fn draw_fan_2(space: &Space, camera: Location, spread: Vec3) -> Vec<Line> {
+fn draw_fan_2(space: &Space, camera: Location, spread: f32) -> 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);
+	for y in itertools_num::linspace(-spread, spread, 101) {
+		draw_ray_2(&mut gc, space, camera, vec3(1., y, 0.));
 	}
 	gc
 }
 
 trait Renderable {
-	fn render(&self) -> Vec<Face>;
+	fn render(&self) -> Vec<Line>;
 }
 
 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()
+	fn render(&self) -> Vec<Line> {
+		let lines = 17;
+		let step = 2. * std::f32::consts::PI / lines as f32;
+		let r = 0.5 * (self.outer_radius + self.inner_radius);
+		let w = 0.5 * (self.outer_radius - self.inner_radius);
+		let l = vec3(0., self.external_halflength, 0.);
+		let along = (0..lines)
+			.map(|k| k as f32 * step)
+			.map(Vec2::from_angle)
+			.map(|d| vec3(d.x, 0., d.y))
+			.map(|d| draw_rect(r * d, w * d, l));
+		let caps = iproduct!([self.inner_radius, self.outer_radius], [-l, l])
+			.map(|(r, l)| draw_ellipse(l, vec3(r, 0., 0.), vec3(0., 0., r)));
+		chain!(along, caps).collect()
 	}
 }

src/bin/wireframe/viewport.rs

@@ -1,147 +0,0 @@
-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

@@ -82,7 +82,12 @@ impl QuadraticAccelerator {
 	}
 
 	pub fn x(&self, u: f32) -> f32 {
-		extend_linear(u, |u| (self.a() * u.abs() + self.b()) * u, self.internal, self.external)
+		extend_linear(
+			u,
+			|u| (self.a() * u.abs() + self.b()) * u,
+			self.internal,
+			self.external,
+		)
 	}
 	pub fn u(&self, x: f32) -> f32 {
 		extend_linear(
@@ -93,7 +98,12 @@ impl QuadraticAccelerator {
 		)
 	}
 	pub fn dx(&self, u: f32) -> f32 {
-		extend_const(u, |u| 2.0 * self.a() * u.abs() + self.b(), self.internal, 1.0)
+		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)
@@ -110,7 +120,7 @@ impl QuadraticAccelerator {
 
 #[cfg(test)]
 mod test {
-	use approx::{abs_diff_eq, assert_abs_diff_eq};
+	use approx::{abs_diff_eq, assert_abs_diff_eq, AbsDiffEq};
 
 	use super::*;
 
@@ -138,12 +148,28 @@ mod test {
 
 	#[test]
 	fn test_smoothstep_limiter() {
-		test_limiter(SmoothstepLimiter { min: 20.0, max: 30.0 }, 20.0, 30.0, 1.0 / 32.0);
+		test_limiter(
+			SmoothstepLimiter {
+				min: 20.0,
+				max: 30.0,
+			},
+			20.0,
+			30.0,
+			1.0 / 32.0,
+		);
 	}
 
 	#[test]
 	fn test_smootherstep_limiter() {
-		test_limiter(SmootherstepLimiter { min: 20.0, max: 30.0 }, 20.0, 30.0, 1.0 / 32.0);
+		test_limiter(
+			SmootherstepLimiter {
+				min: 20.0,
+				max: 30.0,
+			},
+			20.0,
+			30.0,
+			1.0 / 32.0,
+		);
 	}
 
 	#[test]

src/ifaces.rs

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

src/lib.rs

@@ -4,7 +4,6 @@ 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;

src/mathx.rs

@@ -36,8 +36,10 @@ 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();
+		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)
 	}
 }
@@ -140,22 +142,6 @@ where
 	}
 }
 
-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::*;

src/mesh_loader.rs

@@ -5,7 +5,7 @@ use std::io;
 struct ObjVertex {
 	vertex: usize,
 	normal: usize,
-	// tex_coord: usize,
+	tex_coord: usize,
 }
 
 #[derive(Copy, Clone, Debug)]
@@ -35,11 +35,14 @@ impl ObjMesh {
 	}
 
 	fn parse_fv(desc: &&str) -> ObjVertex {
-		let tokens: Vec<_> = desc.split('/').map(|s| s.parse::<usize>().unwrap() - 1).collect();
+		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],
+			tex_coord: tokens[1],
 			normal: tokens[2],
 		}
 	}

src/mesh_tracer.rs

@@ -8,13 +8,21 @@ pub struct TraceResult {
 	pub normal: Vec3,
 }
 
-pub fn trace_to_mesh_all(mesh: &Mesh, base: Vec3, ray: Vec3) -> impl Iterator<Item = TraceResult> + '_ {
+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 = 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 {

src/riemann.rs

@@ -1,4 +1,4 @@
-use crate::mathx::{make_mat3, Decomp3};
+use crate::mathx::Decomp3;
 use glam::*;
 
 pub type Tens3 = [Mat3; 3];
@@ -15,7 +15,7 @@ pub trait Metric {
 	}
 
 	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 ε.
+		part_deriv(|p| self.at(p), pos, 1.0 / 1024.0) // division by such eps is exact which is good for overall precision
 	}
 
 	fn vec_length_at(&self, at: Vec3, v: Vec3) -> f32 {
@@ -58,7 +58,7 @@ pub fn krist(space: &impl Metric, pos: Vec3) -> Tens3 {
 	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)
+		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>()
 	})
@@ -86,6 +86,14 @@ pub fn contract2(t: Tens3, v: Vec3) -> Vec3 {
 	contract(t, v) * v
 }
 
+fn make_vec3(f: impl Fn(usize) -> f32) -> Vec3 {
+	Vec3::from_array(std::array::from_fn(|i| f(i)))
+}
+
+fn make_mat3(f: impl Fn(usize, usize) -> f32) -> Mat3 {
+	Mat3::from_cols_array_2d(&std::array::from_fn(|i| std::array::from_fn(|j| f(i, j))))
+}
+
 fn make_tens3(f: impl Fn(usize, usize, usize) -> f32) -> Tens3 {
 	std::array::from_fn(|i| make_mat3(|j, k| f(i, j, k)))
 }
@@ -137,7 +145,7 @@ mod tests {
 	use super::*;
 	use approx::assert_abs_diff_eq;
 
-	use glam::{vec3, Mat3};
+	use glam::{mat3, vec3, Mat3};
 	use rand::{Rng, SeedableRng};
 
 	#[test]
@@ -161,7 +169,10 @@ mod tests {
 			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.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.);
@@ -222,9 +233,14 @@ mod tests {
 			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(),
+			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

@@ -1,240 +0,0 @@
-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

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

src/shape/rect.rs

@@ -1,90 +0,0 @@
-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

@@ -1,17 +1,18 @@
 use glam::{vec3, Mat3, Vec3};
 
 use crate::riemann::Metric;
-use crate::shape::YCylinder;
 use crate::types::{Location, Ray};
 
-use super::Tube;
+use super::{Rect, Tube};
 
 pub trait FlatCoordinateSystem<T> {
 	fn flat_to_global(&self, v: T) -> T;
 	fn global_to_flat(&self, v: T) -> T;
 }
 
-pub trait FlatRegion: FlatCoordinateSystem<Vec3> + FlatCoordinateSystem<Ray> + FlatCoordinateSystem<Location> {
+pub trait FlatRegion:
+	FlatCoordinateSystem<Vec3> + FlatCoordinateSystem<Ray> + FlatCoordinateSystem<Location>
+{
 	// Измеряет расстояние до выхода за пределы области вдоль луча ray. Луч задаётся в плоской СК.
 	fn distance_to_boundary(&self, _ray: Ray) -> Option<f32> {
 		None
@@ -21,7 +22,10 @@ pub trait FlatRegion: FlatCoordinateSystem<Vec3> + FlatCoordinateSystem<Ray> + F
 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()
+		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()
@@ -81,9 +85,12 @@ impl FlatCoordinateSystem<Vec3> for InnerCS {
 
 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,
+		Rect {
+			size: vec3(
+				self.0.inner_radius,
+				self.0.internal_halflength,
+				self.0.inner_radius,
+			),
 		}
 		.trace_out_of(ray)
 	}
@@ -99,9 +106,12 @@ impl MetricCS for OuterCS {
 
 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,
+		let inner = Rect {
+			size: vec3(
+				self.0.inner_radius + 1.0,
+				self.0.external_halflength,
+				self.0.inner_radius + 1.0,
+			),
 		};
 		if inner.is_inside(pos) {
 			let Vec3 { x, y: v, z } = pos;
@@ -115,13 +125,17 @@ impl FlatCoordinateSystem<Vec3> for OuterCS {
 	}
 
 	fn global_to_flat(&self, pos: Vec3) -> Vec3 {
-		let inner = YCylinder {
-			radius: self.0.inner_radius + 1.0,
-			half_length: self.0.external_halflength,
+		let inner = Rect {
+			size: vec3(
+				self.0.inner_radius + 1.0,
+				self.0.external_halflength,
+				self.0.inner_radius + 1.0,
+			),
 		};
 		if inner.is_inside(pos) {
 			let Vec3 { x: u, y, z: w } = pos; // в основной СК
-			let v = self.0.v(y) + y.signum() * (self.0.external_halflength - self.0.internal_halflength);
+			let v = self.0.v(y)
+				+ y.signum() * (self.0.external_halflength - self.0.internal_halflength);
 			vec3(u, v, w) // в плоском продолжении СК Outer на область Inner
 		} else {
 			pos
@@ -131,9 +145,12 @@ impl FlatCoordinateSystem<Vec3> for OuterCS {
 
 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,
+		Rect {
+			size: vec3(
+				self.0.outer_radius,
+				self.0.external_halflength,
+				self.0.outer_radius,
+			),
 		}
 		.trace_into(ray)
 	}
@@ -210,12 +227,20 @@ mod tests {
 			}),
 			Location {
 				pos: vec3(0., -0.5, 0.),
-				rot: mat3(vec3(0., 0.25, 0.), vec3(-1. / 3., 0., 0.), vec3(0., 0., 0.2))
+				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)) {
+	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 {
@@ -231,7 +256,14 @@ mod tests {
 				);
 			};
 		}
-		fn check_range(name_a: &str, a: Vec3, range_a: (Vec3, Vec3), name_b: &str, b: Vec3, range_b: (Vec3, Vec3)) {
+		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, ε) {
@@ -252,7 +284,14 @@ mod tests {
 				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);
+					check_range(
+						"global",
+						pos_global,
+						range_global,
+						"flat",
+						pos_flat,
+						range_flat,
+					);
 					assert_eq_at!(
 						pos_global,
 						region
@@ -282,7 +321,14 @@ mod tests {
 				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);
+					check_range(
+						"flat",
+						pos_flat,
+						range_flat,
+						"global",
+						pos_global,
+						range_global,
+					);
 					assert_eq_at!(
 						pos_flat,
 						region
@@ -343,16 +389,15 @@ mod tests {
 			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)),
+			(vec3(-30.0, -300.0, -30.0), vec3(30.0, -1.0, 30.0)),
+			(vec3(-30.0, -300.0, -30.0), vec3(30.0, -200.20016, 30.0)),
 		);
 		test_flat_region(
 			&mapper,
-			(vec3(-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)),
+			(vec3(-30.0, 1.0, -30.0), vec3(30.0, 300.0, 30.0)),
+			(vec3(-30.0, 200.20016, -30.0), vec3(30.0, 300.0, 30.0)),
 		);
 		test_flat_region(
 			&mapper,
@@ -449,9 +494,9 @@ mod tests {
 			}
 		}
 		// accelerating
-		for x in linspace(-21., 21., 20) {
+		for x in linspace(-29., 29., 20) {
 			for y in linspace(1., 299., 20) {
-				for z in linspace(-21., 21., 20) {
+				for z in linspace(-29., 29., 20) {
 					let v = mapper
 						.global_to_flat(Location {
 							pos: vec3(x, y, z),

src/tube/metric.rs

@@ -42,11 +42,10 @@ impl Tube {
 
 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 sx = self.fx().value(pos.x);
 		let sy = self.fy().du(pos.y);
-		let s = sr + sy - sr * sy;
-		assert!(sr.is_finite());
+		let s = sx + sy - sx * sy;
+		assert!(sx.is_finite());
 		assert!(sy.is_finite());
 		assert!(sy > 0.0);
 		Decomp3 {
@@ -56,23 +55,17 @@ impl Metric for Tube {
 	}
 
 	fn part_derivs_at(&self, pos: Vec3) -> Tens3 {
-		let r = f32::hypot(pos.x, pos.z);
-		let sr = self.fx().value(r);
+		let sx = self.fx().value(pos.x);
 		let sy = self.fy().du(pos.y);
-		let s = sr + sy - sr * sy;
-		let dsr_dr = self.fx().derivative(r);
+		let s = sx + sy - sx * sy;
+		let dsx_dx = self.fx().derivative(pos.x);
 		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)
-		};
+		let ds2_dx = 2.0 * s * (1.0 - sy) * dsx_dx;
+		let ds2_dy = 2.0 * s * (1.0 - sx) * dsy_dy;
 		[
 			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.]),
+			Mat3::from_cols_array(&[0., 0., 0., 0., 0., 0., 0., 0., 0.]),
 		]
 	}
 }
@@ -108,17 +101,24 @@ mod test {
 		let epsilon = 1.0e-3;
 		let margin = 1.0 / 16.0;
 		let mul = 1.0 + margin;
-		for x in itertools_num::linspace(-mul * testee.outer_radius, mul * testee.outer_radius, 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) {
+		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"
-					);
+					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");
 				}
 			}
 		}
@@ -142,8 +142,8 @@ mod test {
 			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 Δ = 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);
@@ -168,12 +168,20 @@ mod test {
 		let ε = 1e-3;
 		let off = 10.0;
 		let steps = 4096;
-		for ax in linspace(-space.tube.inner_radius + ε, space.tube.inner_radius - ε, 20) {
-			for bx in linspace(-space.tube.inner_radius + ε, space.tube.inner_radius - ε, 20) {
+		for ax in linspace(
+			-space.tube.inner_radius + ε,
+			space.tube.inner_radius - ε,
+			20,
+		) {
+			for bx in linspace(
+				-space.tube.inner_radius + ε,
+				space.tube.inner_radius - ε,
+				20,
+			) {
 				let a = 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 Δ = 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);
@@ -200,8 +208,8 @@ mod test {
 		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 Δ = 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);
@@ -227,8 +235,8 @@ mod test {
 		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 Δ = 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);

src/tube/mod.rs

@@ -1,4 +1,4 @@
-use glam::{f32, Mat3, Vec3};
+use glam::{bool, f32, vec3, Mat3, Vec3};
 
 use crate::ifaces::{DebugTraceable, RayPath, Traceable};
 use coords::{FlatCoordinateSystem, InnerCS, OuterCS};
@@ -28,12 +28,10 @@ pub enum Subspace {
 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 {
+		} else if pt.x.abs() > self.tube.outer_radius {
+			Outer
+		} else if pt.x.abs() > self.tube.inner_radius {
 			Boundary
 		} else {
 			Inner
@@ -52,7 +50,8 @@ impl Space {
 	/// Выполняет один шаг перемещения. Работает в любой части пространства.
 	/// 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 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,
@@ -92,19 +91,36 @@ impl Space {
 		}
 	}
 
+	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 })
+			.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> {
+	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));
+				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))
@@ -130,25 +146,16 @@ impl Space {
 			.collect()
 	}
 
-	pub fn line(&self, a: Vec3, b: Vec3, step: f32) -> Vec<Ray> {
+	pub fn line(&self, a: Vec3, b: Vec3, step: f32) -> Vec<Vec3> {
 		match self.which_subspace(a) {
-			Outer => vec![Ray {
-				pos: b,
-				dir: (b - a).normalize(),
-			}],
+			Outer => vec![b],
 			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,
-						})
-					})
+					.map(|k| cs.flat_to_global(a.lerp(b, k as f32 / n as f32)))
 					.collect()
 			}
 			Boundary => panic!("Can't draw a line here!"),
@@ -203,9 +210,9 @@ impl DebugTraceable for Space {
 		let mut hits = vec![];
 		let mut ray = self.camera_ray_to_abs(camera, ray);
 
-		let trace_to_flat = |points: &mut Vec<Ray>, ray| {
+		let trace_to_flat = |points: &mut Vec<Vec3>, ray| {
 			for ray in self.trace_iter(ray).skip(1) {
-				points.push(ray);
+				points.push(ray.pos);
 				if let Some(hitter) = self.obj_hitter(ray.pos) {
 					return (ray, hitter(self, ray));
 				}
@@ -213,16 +220,212 @@ impl DebugTraceable for Space {
 			unreachable!("Space::trace_iter terminated!")
 		};
 
-		points.push(ray);
+		points.push(ray.pos);
 		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 });
+				return (
+					hits,
+					RayPath {
+						points,
+						end_dir: ray_into_flat.dir.normalize(),
+					},
+				);
 			};
-			points.extend(self.line(ray_into_flat.pos, ray_outta_flat.pos, 1.0));
+			points.extend(self.line(ray_into_flat.pos, ray_outta_flat.pos, 10.0));
 			ray = ray_outta_flat;
 		}
 		panic!("tracing didn't terminate");
 	}
 }
+
+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(),
+		}
+	}
+
+	fn is_inside(&self, pt: Vec3) -> 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: vec3(2.0, 3.0, 2.0),
+			dir: vec3(4.0, 5.0, 4.0)
+		}),
+		Ray {
+			pos: vec3(2.0, 3.0, 2.0),
+			dir: vec3(4.0, 5.0, 4.0)
+		}
+	);
+	assert_eq!(
+		Rect::flip_ray(Ray {
+			pos: vec3(2.0, 3.0, 2.0),
+			dir: vec3(-4.0, 5.0, -4.0)
+		}),
+		Ray {
+			pos: vec3(-2.0, 3.0, -2.0),
+			dir: vec3(4.0, 5.0, 4.0)
+		}
+	);
+	assert_eq!(
+		Rect::flip_ray(Ray {
+			pos: vec3(2.0, 3.0, 2.0),
+			dir: vec3(4.0, -5.0, 4.0)
+		}),
+		Ray {
+			pos: vec3(2.0, -3.0, 2.0),
+			dir: vec3(4.0, 5.0, 4.0)
+		}
+	);
+	assert_eq!(
+		Rect::flip_ray(Ray {
+			pos: vec3(2.0, 3.0, 2.0),
+			dir: vec3(4.0, 0.0, 4.0)
+		}),
+		Ray {
+			pos: vec3(2.0, 3.0, 2.0),
+			dir: vec3(4.0, 0.0, 4.0)
+		}
+	);
+
+	let r = Rect {
+		size: vec3(2.0, 3.0, 2.0),
+	};
+
+	assert_eq!(
+		r.trace_into(Ray {
+			pos: vec3(3.0, 3.0, 3.0),
+			dir: vec3(1.0, 1.0, 1.0)
+		}),
+		None
+	);
+	assert_eq!(
+		r.trace_into(Ray {
+			pos: vec3(-3.0, 2.0, -3.0),
+			dir: vec3(1.0, 0.0, 1.0)
+		}),
+		Some(1.0)
+	);
+	assert_eq!(
+		r.trace_into(Ray {
+			pos: vec3(-3.0, 2.0, -3.0),
+			dir: vec3(-1.0, 0.0, -1.0)
+		}),
+		None
+	);
+	assert_eq!(
+		r.trace_into(Ray {
+			pos: vec3(-3.0, 1.0, -3.0),
+			dir: vec3(2.0, 2.0, 2.0)
+		}),
+		Some(0.5)
+	);
+	assert_eq!(
+		r.trace_into(Ray {
+			pos: vec3(-3.0, 2.1, -3.0),
+			dir: vec3(2.0, 2.0, 2.0)
+		}),
+		None
+	);
+
+	assert_eq!(
+		r.trace_into(Ray {
+			pos: vec3(2.0, 3.0, 2.0),
+			dir: vec3(1.0, 1.0, 1.0)
+		}),
+		None
+	);
+	assert_eq!(
+		r.trace_into(Ray {
+			pos: vec3(-2.0, 3.0, -2.0),
+			dir: vec3(-1.0, 1.0, -1.0)
+		}),
+		None
+	);
+	assert_eq!(
+		r.trace_into(Ray {
+			pos: vec3(2.0, 3.0, 2.0),
+			dir: vec3(-1.0, -1.0, -1.0)
+		}),
+		Some(0.0)
+	);
+	assert_eq!(
+		r.trace_into(Ray {
+			pos: vec3(2.0, -3.0, 2.0),
+			dir: vec3(-1.0, 1.0, -1.0)
+		}),
+		Some(0.0)
+	);
+
+	assert_eq!(
+		r.trace_out_of(Ray {
+			pos: vec3(0.0, 0.0, 0.0),
+			dir: vec3(1.0, 1.0, 1.0)
+		}),
+		Some(2.0)
+	);
+	assert_eq!(
+		r.trace_out_of(Ray {
+			pos: vec3(0.0, 0.0, 0.0),
+			dir: vec3(0.0, 1.0, 0.0)
+		}),
+		Some(3.0)
+	);
+	assert_eq!(
+		r.trace_out_of(Ray {
+			pos: vec3(0.0, 1.0, 0.0),
+			dir: vec3(0.0, -1.0, 0.0)
+		}),
+		Some(4.0)
+	);
+	assert_eq!(
+		r.trace_out_of(Ray {
+			pos: vec3(1.0, 1.0, 1.0),
+			dir: vec3(0.0, -1.0, 0.0)
+		}),
+		Some(4.0)
+	);
+	assert_eq!(
+		r.trace_out_of(Ray {
+			pos: vec3(2.0, 3.0, 2.0),
+			dir: vec3(1.0, 1.0, 1.0)
+		}),
+		Some(0.0)
+	);
+}

src/types.rs

@@ -6,10 +6,6 @@ pub struct Ray {
 	pub dir: Vec3,
 }
 
-pub fn ray(pos: Vec3, dir: Vec3) -> Ray {
-	Ray { pos, dir }
-}
-
 impl Ray {
 	pub fn forward(&self, dist: f32) -> Ray {
 		Ray {

src/utils.rs

@@ -43,8 +43,16 @@ mod tests {
 			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 = ε);
+		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,
@@ -52,8 +60,16 @@ mod tests {
 			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 = ε);
+		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(
@@ -62,7 +78,15 @@ mod tests {
 			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 = ε);
+		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 = ε
+		);
 	}
 }