Pre-generate normals, for smooth look

Anti-aliasing works better this way

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'
@@ -21,9 +22,11 @@ show-image = "0.14.0"
 flo_draw = "0.3.1"
 flo_canvas = "0.3.1"
 itertools-num = "0.1.3"
-glium = "0.35.0"
+winit = "0.29"
-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 +1,2 @@
 hard_tabs = true
+max_width = 120

src/bin/flat/main.rs

@@ -52,20 +52,12 @@ pub fn main() {
 			let cam1 = put_object(&space.tube, vec3(-500., 0., 0.), Mat3::IDENTITY);
 			let cam2 = put_object(
 				&space.tube,
-				vec3(
+				vec3(-2.5 * tube.outer_radius, 1.25 * tube.external_halflength, 0.),
-					-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(
+				vec3(0.25 * tube.inner_radius, 0.25 * tube.external_halflength, 0.),
-					0.25 * tube.inner_radius,
-					0.25 * tube.external_halflength,
-					0.,
-				),
 				mat3(vec3(0., -1., 0.), vec3(1., 0., 0.), vec3(0., 0., 1.)),
 			);
 
@@ -134,11 +126,7 @@ pub fn main() {
 							let dir = Vec2::from_angle(φ);
 							let dir = vec3(dir.x, dir.y, 0.);
 							let dir = obj.loc.rot * dir * d;
-							space
+							space.trace_iter(Ray { pos, dir }).nth(n as usize).unwrap().pos
-								.trace_iter(Ray { pos, dir })
-								.nth(n as usize)
-								.unwrap()
-								.pos
 						}),
 				);
 			}
@@ -146,6 +134,7 @@ 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);
@@ -166,13 +155,10 @@ 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.x, pt.y);
+		gc.line_to(pt.pos.x, pt.pos.y);
 	}
-	let end_pos = *path
+	let end_pos = *path.points.last().expect("the starting point is always in the path");
-		.points
+	let dir_pos = end_pos.forward(1000. / DT).pos;
-		.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();
 }
@@ -191,11 +177,7 @@ 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(
+		rot: mat3(vec3(dir.x, dir.y, 0.), vec3(-dir.y, dir.x, 0.), vec3(0., 0., 1.)),
-			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| {
@@ -215,8 +197,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;
+		let n = (STEP / DT).floor() as i32;
-		for _ in 0..N {
+		for _ in 0..n {
 			loc = space.move_step(loc, v * DT);
 		}
 		draw(&loc);

src/bin/mesh/main.rs

@@ -1,12 +1,14 @@
 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;
@@ -79,9 +81,7 @@ fn render(mesh: &Mesh, camera: impl Fn(Vec2) -> (Vec3, Vec3)) -> Image {
 			} else {
 				bkg
 			};
-			let color = (color * 255.0)
+			let color = (color * 255.0).as_ivec3().clamp(IVec3::splat(0), IVec3::splat(255));
-				.as_ivec3()
-				.clamp(IVec3::splat(0), IVec3::splat(255));
 			img.put_pixel(x, y, Color(color.x as u8, color.y as u8, color.z as u8));
 		}
 	}
@@ -109,6 +109,10 @@ 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();
@@ -121,15 +125,21 @@ 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 m_view = ypr_to_mat(vec3(
+			let proj = PROJS[PROJ_INDEX.load(Relaxed)];
-				(135.0 + phi as f32) * PI / 180.0,
+			let m_view = ypr_to_mat(vec3((135.0 + phi as f32) * PI / 180.0, -30.0 * PI / 180.0, 0.0f32));
-				-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

@@ -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

@@ -1,73 +1,38 @@
 use std::time::Instant;
 
-use glam::{mat4, vec2, vec3, vec4, Mat4, Vec3};
+use glam::{uvec2, vec3, 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
 
-#[derive(Copy, Clone)]
+fn prepare_scene(device: &wgpu::Device) -> (Vec<meshes::Mesh>, Vec<lines::Line>) {
-struct Vertex {
+	let (meshes, lines) = scene::build();
-	position: [f32; 3],
+	let meshes = meshes
-}
-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| {
+		.map(|mesh| meshes::Mesh::new_list(device, meshes::Attrs { color: mesh.color }, mesh.tris))
-			let color = line.color;
+		.collect();
-			let mode;
+	let lines = lines
-			let data: Vec<Vertex>;
+		.into_iter()
-			match line.line {
+		.map(|line| lines::Line::new_strip(device, lines::Attrs { color: line.color }, line.pts))
-				scene::Line::Lines(_) => todo!(),
+		.collect();
-				scene::Line::Strip(pts) => {
+	(meshes, lines)
-					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};
 
@@ -92,15 +57,48 @@ mod camctl {
 		}
 
 		pub fn rotate_rel_ypr(&mut self, ypr: Vec3) {
-			self.rotate_rel_quat(Quat::from_euler(glam::EulerRot::XYZ, ypr.x, ypr.y, ypr.z));
+			self.rotate_rel_quat(Quat::from_euler(glam::EulerRot::ZYX, ypr.x, ypr.y, ypr.z));
 		}
 
-		pub fn rotate_rel_quat(&mut self, rot: Quat) {
+		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};
@@ -142,162 +140,278 @@ 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::Numpad9), vec3(1., 0., 0.)),
+	(PhysicalKey::Code(KeyCode::Numpad4), vec3(1., 0., 0.)),
-	(PhysicalKey::Code(KeyCode::Numpad7), 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::Numpad4), vec3(0., 0., 1.)),
+	(PhysicalKey::Code(KeyCode::Numpad9), vec3(0., 0., 1.)),
-	(PhysicalKey::Code(KeyCode::Numpad6), vec3(0., 0., -1.)),
+	(PhysicalKey::Code(KeyCode::Numpad7), vec3(0., 0., -1.)),
 ];
 
-fn main() {
+struct State<'a> {
-	let event_loop = EventLoop::builder().build().unwrap();
+	device: wgpu::Device,
-	let cfg = ConfigTemplateBuilder::new().with_multisampling(8);
+	queue: wgpu::Queue,
-	let (window, display) = SimpleWindowBuilder::new()
-		.with_config_template_builder(cfg)
-		.with_title("Refraction: Wireframe")
-		.build(&event_loop);
 
-	let vs_src = include_str!("ray.v.glsl");
+	fps: fps::Counter,
-	let fs_src = include_str!("ray.f.glsl");
+	kbd: keyctl::Keyboard,
-	let program = Program::from_source(&display, vs_src, fs_src, None).unwrap();
+	t1: Instant,
 
-	let scene = prepare_scene(&display);
+	viewport: viewport::Viewport<'a>,
+	cam_loc: camctl::CameraLocation,
+	cam_obj: camera::Camera,
+	line_rend: lines::LineRenderer,
+	mesh_rend: meshes::Renderer,
 
-	let mut kbd = keyctl::Keyboard::new();
+	lines: Vec<lines::Line>,
-	let mut cam = camctl::CameraLocation::new();
+	meshes: Vec<meshes::Mesh>,
-	let mut t1 = Instant::now();
 
-	#[allow(deprecated)]
+	window: &'a Window,
-	event_loop
+}
-		.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));
 
-					let size = window.inner_size();
+impl<'a> State<'a> {
-					let size = vec2(size.width as f32, size.height as f32).normalize()
+	async fn new(window: &'a Window) -> State<'a> {
-						* std::f32::consts::SQRT_2;
+		let size = window.inner_size();
-					let proj = make_proj_matrix(vec3(size.x, size.y, 2.), (1., 4096.));
 
-					let my_to_gl = mat4(
+		let instance = wgpu::Instance::new(wgpu::InstanceDescriptor {
-						vec4(0., 0., 1., 0.),
+			backends: wgpu::Backends::PRIMARY,
-						vec4(-1., 0., 0., 0.),
+			..Default::default()
-						vec4(0., 1., 0., 0.),
+		});
-						vec4(0., 0., 0., 1.),
+		let surface = instance.create_surface(window).unwrap();
-					);
+		let adapter = instance
-					let view = my_to_gl * cam.view_mtx();
+			.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 mut target = display.draw();
+		let viewport = viewport::Viewport::new(&adapter, &device, surface, uvec2(size.width, size.height));
-					target.clear_color(0.0, 0.0, 0.0, 0.0);
 
-					let mvp = proj * view;
+		let kbd = keyctl::Keyboard::new();
-					let params = DrawParameters {
+		let cam_loc = camctl::CameraLocation::new();
-						blend: glium::Blend {
+		let t1 = Instant::now();
-							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 {
+		let cam_obj = camera::Camera::new(&device);
-					device_id: _,
+		let line_rend = lines::LineRenderer::new(
-					event,
+			&device,
-					is_synthetic: _,
+			cam_obj.bind_group_layout(),
-				} => {
+			viewport.format(),
-					kbd.set_key_state(event.physical_key, event.state);
+			Some(wgpu::DepthStencilState {
-				}
+				format: wgpu::TextureFormat::Depth24Plus,
-
+				depth_write_enabled: false,
-				WindowEvent::Resized(window_size) => {
+				depth_compare: wgpu::CompareFunction::LessEqual,
-					display.resize(window_size.into());
+				stencil: wgpu::StencilState::default(),
-				}
+				bias: wgpu::DepthBiasState::default(),
-
+			}),
-				WindowEvent::CloseRequested => {
+			wgpu::MultisampleState {
-					window_target.exit();
+				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,
+			},
+		);
 
-			Event::AboutToWait => {
+		let (meshes, lines) = prepare_scene(&device);
-				window.request_redraw();
+
+		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();
 }
 
-/// Make a projection matrix, assuming input coordinates are (right, up, forward).
+fn main() {
-///
+	pollster::block_on(run());
-/// `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 fps {
-mod tests {
+	use std::time::{Duration, Instant};
-	use super::*;
-	use approx::assert_abs_diff_eq;
-	use glam::vec3;
 
-	#[test]
+	pub struct Counter {
-	fn test_proj_matrix() {
+		fps: f32,
-		let m = make_proj_matrix(vec3(2., 3., 4.), (0.5, 20.0));
+		t1: Instant,
+		frames: u32,
+	}
 
-		let v = m * vec4(2., 3., 4., 1.);
+	impl Counter {
-		assert_abs_diff_eq!(v.x / v.w, 1.0);
+		pub fn new() -> Self {
-		assert_abs_diff_eq!(v.y / v.w, 1.0);
+			Self {
-		assert!(-v.w < v.z && v.z < v.w, "z out of range in {v}");
+				fps: 0.,
+				t1: Instant::now(),
+				frames: 0,
+			}
+		}
 
-		let v = m * vec4(2., 3., 0.5, 1.);
+		pub fn get(&self) -> f32 {
-		assert_abs_diff_eq!(v.x / v.w, 8.0);
+			self.fps
-		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.);
+		pub fn on_frame(&mut self) {
-		assert_abs_diff_eq!(v.x / v.w, 0.2);
+			self.frames += 1;
-		assert_abs_diff_eq!(v.y / v.w, 0.2);
+			let t2 = Instant::now();
-		assert_abs_diff_eq!(v.z / v.w, 1.0);
+			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.f.glsl

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

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

@@ -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

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

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

src/ifaces.rs

@@ -1,5 +1,4 @@
 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.
@@ -19,8 +18,7 @@ pub trait OptimizedTraceable: Traceable {
 }
 
 pub struct RayPath {
-	pub points: Vec<Vec3>,
+	pub points: Vec<Ray>,
-	pub end_dir: Vec3,
 }
 
 pub trait DebugTraceable: Traceable {

src/lib.rs

@@ -4,6 +4,7 @@ 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,10 +36,8 @@ 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
+		let fz = (self.z_axis - self.z_axis.project_onto_normalized(fx) - self.z_axis.project_onto_normalized(fy))
-			- self.z_axis.project_onto_normalized(fx)
+			.normalize();
-			- self.z_axis.project_onto_normalized(fy))
-		.normalize();
 		Self::from_cols(fx, fy, fz)
 	}
 }
@@ -142,6 +140,22 @@ 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,14 +35,11 @@ impl ObjMesh {
 	}
 
 	fn parse_fv(desc: &&str) -> ObjVertex {
-		let tokens: Vec<_> = desc
+		let tokens: Vec<_> = desc.split('/').map(|s| s.parse::<usize>().unwrap() - 1).collect();
-			.split('/')
-			.map(|s| s.parse::<usize>().unwrap() - 1)
-			.collect();
 		assert_eq!(tokens.len(), 3);
 		ObjVertex {
 			vertex: tokens[0],
-			tex_coord: tokens[1],
+			// tex_coord: tokens[1],
 			normal: tokens[2],
 		}
 	}

src/mesh_tracer.rs

@@ -8,21 +8,13 @@ pub struct TraceResult {
 	pub normal: Vec3,
 }
 
-pub fn trace_to_mesh_all(
+pub fn trace_to_mesh_all(mesh: &Mesh, base: Vec3, ray: Vec3) -> impl Iterator<Item = TraceResult> + '_ {
-	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(
+		let m = mat3(f.vertices[1] - f.vertices[0], f.vertices[2] - f.vertices[0], -ray);
-			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::Decomp3;
+use crate::mathx::{make_mat3, 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 / 1024.0) // division by such eps is exact which is good for overall precision
+		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 {
@@ -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..2)
+		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>()
 	})
@@ -86,14 +86,6 @@ 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)))
 }
@@ -145,7 +137,7 @@ mod tests {
 	use super::*;
 	use approx::assert_abs_diff_eq;
 
-	use glam::{mat3, vec3, Mat3};
+	use glam::{vec3, Mat3};
 	use rand::{Rng, SeedableRng};
 
 	#[test]
@@ -169,10 +161,7 @@ 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!(
+		assert_eq!(metric.part_derivs_at(rng.gen()), [Mat3::ZERO, Mat3::ZERO, Mat3::ZERO]);
-			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.);
@@ -233,14 +222,9 @@ mod tests {
 			vec3(3., 6.25, 0.)
 		);
 		assert_abs_diff_eq!(
-			trace_iter(
+			trace_iter(&metric, vec3(3., 5., 0.), vec3(0.5, 0.25, 0.), std::f32::consts::SQRT_2)
-				&metric,
+				.nth(7)
-				vec3(3., 5., 0.),
+				.unwrap(),
-				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

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

src/tube/metric.rs

@@ -42,10 +42,11 @@ impl Tube {
 
 impl Metric for Tube {
 	fn sqrt_at(&self, pos: Vec3) -> Decomp3 {
-		let sx = self.fx().value(pos.x);
+		let r = f32::hypot(pos.x, pos.z);
+		let sr = self.fx().value(r);
 		let sy = self.fy().du(pos.y);
-		let s = sx + sy - sx * sy;
+		let s = sr + sy - sr * sy;
-		assert!(sx.is_finite());
+		assert!(sr.is_finite());
 		assert!(sy.is_finite());
 		assert!(sy > 0.0);
 		Decomp3 {
@@ -55,17 +56,23 @@ impl Metric for Tube {
 	}
 
 	fn part_derivs_at(&self, pos: Vec3) -> Tens3 {
-		let sx = self.fx().value(pos.x);
+		let r = f32::hypot(pos.x, pos.z);
+		let sr = self.fx().value(r);
 		let sy = self.fy().du(pos.y);
-		let s = sx + sy - sx * sy;
+		let s = sr + sy - sr * sy;
-		let dsx_dx = self.fx().derivative(pos.x);
+		let dsr_dr = self.fx().derivative(r);
 		let dsy_dy = self.fy().d2u(pos.y);
-		let ds2_dx = 2.0 * s * (1.0 - sy) * dsx_dx;
+		let ds2_dr = 2.0 * s * (1.0 - sy) * dsr_dr;
-		let ds2_dy = 2.0 * s * (1.0 - sx) * dsy_dy;
+		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., 0., 0., 0., 0., 0.]),
+			Mat3::from_cols_array(&[0., 0., 0., 0., ds2_dz, 0., 0., 0., 0.]),
 		]
 	}
 }
@@ -101,24 +108,17 @@ 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 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 y in itertools_num::linspace(
+				for z in itertools_num::linspace(-mul * testee.outer_radius, mul * testee.outer_radius, 20) {
-				-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 =
+					let eq = (0..3).all(|coord| computed[coord].abs_diff_eq(reference[coord], epsilon));
-						(0..2).all(|coord| computed[coord].abs_diff_eq(reference[coord], epsilon));
+					assert!(
-					assert!(eq, "Bad derivative computation at {pos}:\n  explicit:  {computed:?}\n  numerical: {reference:?}\n");
+						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 δ = vec3(bx - ax, 2.0 * (space.tube.internal_halflength + off), 0.);
-				let dir = Δ / (steps as f32);
+				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,20 +168,12 @@ mod test {
 		let ε = 1e-3;
 		let off = 10.0;
 		let steps = 4096;
-		for ax in linspace(
+		for ax in linspace(-space.tube.inner_radius + ε, space.tube.inner_radius - ε, 20) {
-			-space.tube.inner_radius + ε,
+			for bx in linspace(-space.tube.inner_radius + ε, space.tube.inner_radius - ε, 20) {
-			space.tube.inner_radius - ε,
-			20,
-		) {
-			for bx in linspace(
-				-space.tube.inner_radius + ε,
-				space.tube.inner_radius - ε,
-				20,
-			) {
 				let a = 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 δ = vec3(bx - ax, 2.0 * (space.tube.internal_halflength + off), 0.);
-				let dir = Δ / (steps as f32);
+				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);
@@ -208,8 +200,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 δ = vec3(0.0, 2.0 * (space.tube.internal_halflength + off), 0.);
-			let dir = Δ / (steps as f32);
+			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);
@@ -235,8 +227,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 δ = vec3(0.0, 2.0 * (space.tube.external_halflength + off), 0.);
-			let dir = Δ / (steps as f32);
+			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::{bool, f32, vec3, Mat3, Vec3};
+use glam::{f32, Mat3, Vec3};
 
 use crate::ifaces::{DebugTraceable, RayPath, Traceable};
 use coords::{FlatCoordinateSystem, InnerCS, OuterCS};
@@ -28,10 +28,12 @@ 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 pt.x.abs() > self.tube.outer_radius {
+		} else if r > self.tube.inner_radius {
-			Outer
-		} else if pt.x.abs() > self.tube.inner_radius {
 			Boundary
 		} else {
 			Inner
@@ -50,8 +52,7 @@ impl Space {
 	/// Выполняет один шаг перемещения. Работает в любой части пространства.
 	/// off задаётся в локальной СК. Рекомендуется считать небольшими шагами.
 	pub fn move_step(&self, loc: Location, off: Vec3) -> Location {
-		let corr =
+		let corr = Mat3::IDENTITY - riemann::contract(riemann::krist(&self.tube, loc.pos), loc.rot * off);
-			Mat3::IDENTITY - riemann::contract(riemann::krist(&self.tube, loc.pos), loc.rot * off);
 		let p = loc.pos + corr * loc.rot * off;
 		Location {
 			pos: p,
@@ -91,36 +92,19 @@ 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 {
+			.map(|&Object { id, loc, r }| Object { id, loc: tfm(loc), r })
-				id,
-				loc: tfm(loc),
-				r,
-			})
 			.collect()
 	}
 
-	fn hit_objects(
+	fn hit_objects(objs: &[Object], ray: Ray, limit: Option<f32>, globalize: impl Fn(Vec3) -> Vec3) -> Vec<Hit> {
-		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)
+				let diff = rel.dot(ray.dir).powi(2) - ray.dir.length_squared() * (rel.length_squared() - obj.r.powi(2));
-					- ray.dir.length_squared() * (rel.length_squared() - obj.r.powi(2));
 				if diff > 0.0 {
 					let t = (-rel.dot(ray.dir) - diff.sqrt()) / ray.dir.length_squared();
 					Some((obj, t))
@@ -146,16 +130,25 @@ impl Space {
 			.collect()
 	}
 
-	pub fn line(&self, a: Vec3, b: Vec3, step: f32) -> Vec<Vec3> {
+	pub fn line(&self, a: Vec3, b: Vec3, step: f32) -> Vec<Ray> {
 		match self.which_subspace(a) {
-			Outer => vec![b],
+			Outer => vec![Ray {
+				pos: b,
+				dir: (b - a).normalize(),
+			}],
 			Inner => {
 				let cs = InnerCS(self.tube);
 				let n = ((b - a).length() / step) as usize + 1;
 				let a = cs.global_to_flat(a);
 				let b = cs.global_to_flat(b);
+				let dir = (b - a).normalize();
 				(1..=n)
-					.map(|k| cs.flat_to_global(a.lerp(b, k as f32 / n as f32)))
+					.map(|k| {
+						cs.flat_to_global(Ray {
+							pos: a.lerp(b, k as f32 / n as f32),
+							dir,
+						})
+					})
 					.collect()
 			}
 			Boundary => panic!("Can't draw a line here!"),
@@ -210,9 +203,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<Vec3>, ray| {
+		let trace_to_flat = |points: &mut Vec<Ray>, ray| {
 			for ray in self.trace_iter(ray).skip(1) {
-				points.push(ray.pos);
+				points.push(ray);
 				if let Some(hitter) = self.obj_hitter(ray.pos) {
 					return (ray, hitter(self, ray));
 				}
@@ -220,212 +213,16 @@ impl DebugTraceable for Space {
 			unreachable!("Space::trace_iter terminated!")
 		};
 
-		points.push(ray.pos);
+		points.push(ray);
 		for _ in 0..100 {
 			let (ray_into_flat, ret) = trace_to_flat(&mut points, ray);
 			hits.extend(ret.objects); // TODO fix distance
 			let Some(ray_outta_flat) = ret.end else {
-				return (
+				return (hits, RayPath { points });
-					hits,
-					RayPath {
-						points,
-						end_dir: ray_into_flat.dir.normalize(),
-					},
-				);
 			};
-			points.extend(self.line(ray_into_flat.pos, ray_outta_flat.pos, 10.0));
+			points.extend(self.line(ray_into_flat.pos, ray_outta_flat.pos, 1.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,6 +6,10 @@ 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,16 +43,8 @@ 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!(
+		assert_abs_diff_eq!(loc.rot * vec3(1., 0., 0.), vec3(1. / 3., 0., 0.), epsilon = ε);
-			loc.rot * vec3(1., 0., 0.),
+		assert_abs_diff_eq!(loc.rot * vec3(0., 1., 0.), vec3(0., 1. / 4., 0.), epsilon = ε);
-			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,
@@ -60,16 +52,8 @@ 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!(
+		assert_abs_diff_eq!(loc.rot * vec3(1., 0., 0.), vec3(0., 1. / 4., 0.), epsilon = ε);
-			loc.rot * vec3(1., 0., 0.),
+		assert_abs_diff_eq!(loc.rot * vec3(0., 1., 0.), vec3(-1. / 3., 0., 0.), epsilon = ε);
-			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(
@@ -78,15 +62,7 @@ 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!(
+		assert_abs_diff_eq!(loc.rot * vec3(1., 0., 0.), vec3(1. / 5., 1. / 5., 0.), epsilon = ε);
-			loc.rot * vec3(1., 0., 0.),
+		assert_abs_diff_eq!(loc.rot * vec3(0., 1., 0.), vec3(-4. / 15., 3. / 20., 0.), epsilon = ε);
-			vec3(1. / 5., 1. / 5., 0.),
-			epsilon = ε
-		);
-		assert_abs_diff_eq!(
-			loc.rot * vec3(0., 1., 0.),
-			vec3(-4. / 15., 3. / 20., 0.),
-			epsilon = ε
-		);
 	}
 }