Compare commits
No commits in common. "cbce6ccc444f0ff8d402247c61b7f3f1894d0ea7" and "afc497002353c9069df30c2e6f0700515305e8a8" have entirely different histories.
cbce6ccc44
...
afc4970023
3827
Cargo.lock
generated
3827
Cargo.lock
generated
File diff suppressed because it is too large
Load Diff
15
Cargo.toml
15
Cargo.toml
|
|
@ -1,25 +1,14 @@
|
|||
[package]
|
||||
name = "refraction"
|
||||
name = "hello"
|
||||
version = "0.1.0"
|
||||
edition = "2021"
|
||||
|
||||
# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
|
||||
|
||||
[profile.dev]
|
||||
panic = 'abort'
|
||||
|
||||
[profile.dev.package."*"]
|
||||
opt-level = 3
|
||||
|
||||
[dependencies]
|
||||
rand = "0.8.5"
|
||||
glam = { version = "0.27.0", features = ["approx", "fast-math", "rand"] }
|
||||
glm = "0.2.3"
|
||||
show-image = "0.14.0"
|
||||
flo_draw = "0.3.1"
|
||||
flo_canvas = "0.3.1"
|
||||
itertools-num = "0.1.3"
|
||||
|
||||
[dev-dependencies]
|
||||
approx = "0.5.1"
|
||||
rand = "0.8.5"
|
||||
rand_pcg = "0.3.1"
|
||||
|
|
|
|||
|
|
@ -1,36 +0,0 @@
|
|||
use glam::FloatExt;
|
||||
|
||||
mod bounds {
|
||||
pub trait Pair<T> {}
|
||||
|
||||
impl<T> Pair<T> for (T, T) {}
|
||||
}
|
||||
|
||||
pub trait FloatExt2<T>: bounds::Pair<T> {
|
||||
fn lerp(self, t: T) -> T;
|
||||
fn inverse_lerp(self, y: T) -> T;
|
||||
}
|
||||
|
||||
impl<F: FloatExt> FloatExt2<F> for (F, F) {
|
||||
fn lerp(self, t: F) -> F { F::lerp(self.0, self.1, t) }
|
||||
fn inverse_lerp(self, y: F) -> F { F::inverse_lerp(self.0, self.1, y) }
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod test {
|
||||
use super::FloatExt2;
|
||||
|
||||
#[test]
|
||||
fn test_lerp() {
|
||||
assert_eq!((3., 7.).lerp(-0.5), 1.);
|
||||
assert_eq!((3., 7.).lerp(0.0), 3.);
|
||||
assert_eq!((3., 7.).lerp(0.5), 5.);
|
||||
assert_eq!((3., 7.).lerp(1.0), 7.);
|
||||
assert_eq!((3., 7.).lerp(1.5), 9.);
|
||||
assert_eq!((3., 7.).inverse_lerp(1.), -0.5);
|
||||
assert_eq!((3., 7.).inverse_lerp(3.), 0.0);
|
||||
assert_eq!((3., 7.).inverse_lerp(5.), 0.5);
|
||||
assert_eq!((3., 7.).inverse_lerp(7.), 1.0);
|
||||
assert_eq!((3., 7.).inverse_lerp(9.), 1.5);
|
||||
}
|
||||
}
|
||||
|
|
@ -1,150 +0,0 @@
|
|||
use crate::float_fun::FloatExt2;
|
||||
|
||||
pub trait Limiter {
|
||||
fn value(&self, x: f32) -> f32;
|
||||
fn derivative(&self, x: f32) -> f32;
|
||||
}
|
||||
|
||||
pub struct SmoothstepLimiter {
|
||||
pub min: f32,
|
||||
pub max: f32,
|
||||
}
|
||||
|
||||
impl Limiter for SmoothstepLimiter {
|
||||
fn value(&self, x: f32) -> f32 {
|
||||
let y = (self.min, self.max).inverse_lerp(x.abs()).clamp(0.0, 1.0);
|
||||
3.0 * y * y - 2.0 * y * y * y
|
||||
}
|
||||
fn derivative(&self, x: f32) -> f32 {
|
||||
if x.abs() > self.min && x.abs() < self.max {
|
||||
let t = (self.min, self.max).inverse_lerp(x.abs());
|
||||
6.0 * x.signum() * t * (1.0 - t) / (self.max - self.min)
|
||||
} else {
|
||||
0.0
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
pub struct SmootherstepLimiter {
|
||||
pub min: f32,
|
||||
pub max: f32,
|
||||
}
|
||||
|
||||
impl Limiter for SmootherstepLimiter {
|
||||
fn value(&self, x: f32) -> f32 {
|
||||
let y = (self.min, self.max).inverse_lerp(x.abs()).clamp(0.0, 1.0);
|
||||
6.0 * y.powi(5) - 15.0 * y.powi(4) + 10.0 * y.powi(3)
|
||||
}
|
||||
fn derivative(&self, x: f32) -> f32 {
|
||||
if x.abs() > self.min && x.abs() < self.max {
|
||||
let t = (self.min, self.max).inverse_lerp(x.abs());
|
||||
30.0 * (t * (1.0 - t)).powi(2) * x.signum() / (self.max - self.min)
|
||||
} else {
|
||||
0.0
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
pub struct QuadraticAccelerator {
|
||||
pub internal: f32,
|
||||
pub external: f32,
|
||||
}
|
||||
|
||||
/// Продолжает функцию f с [-lim, lim] линейно в предположении f(±lim) = ±val, f'(±lim) = 1.
|
||||
fn extend_linear(t: f32, f: impl FnOnce(f32) -> f32, lim: f32, val: f32) -> f32 {
|
||||
if t.abs() <= lim { f(t) } else { t + t.signum() * (val - lim) }
|
||||
}
|
||||
|
||||
/// Продолжает функцию f с [-lim, lim] константой в предположении f(±lim) = val, f'(±lim) = 0.
|
||||
fn extend_const(t: f32, f: impl FnOnce(f32) -> f32, lim: f32, val: f32) -> f32 {
|
||||
if t.abs() <= lim { f(t) } else { val }
|
||||
}
|
||||
|
||||
impl QuadraticAccelerator {
|
||||
fn a(&self) -> f32 { -(self.external - self.internal) / self.internal.powi(2) }
|
||||
fn b(&self) -> f32 { 2.0 * self.external / self.internal - 1.0 }
|
||||
fn root(&self, x: f32) -> f32 { (self.b().powi(2) + 4.0 * self.a() * x.abs()).sqrt() }
|
||||
|
||||
pub fn x(&self, u: f32) -> f32 { extend_linear(u, |u| (self.a() * u.abs() + self.b()) * u, self.internal, self.external) }
|
||||
pub fn u(&self, x: f32) -> f32 { extend_linear(x, |x| 0.5 * x.signum() * (-self.b() + self.root(x)) / self.a(), self.external, self.internal) }
|
||||
pub fn dx(&self, u: f32) -> f32 { extend_const(u, |u| 2.0 * self.a() * u.abs() + self.b(), self.internal, 1.0) }
|
||||
pub fn du(&self, x: f32) -> f32 { extend_const(x, |x| 1.0 / self.root(x), self.external, 1.0) }
|
||||
pub fn d2u(&self, x: f32) -> f32 { extend_const(x, |x| -2.0 * x.signum() * self.a() * self.root(x).powi(-3), self.external, 0.0) }
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod test {
|
||||
use approx::{abs_diff_eq, AbsDiffEq, assert_abs_diff_eq};
|
||||
|
||||
use super::*;
|
||||
|
||||
fn test_limiter(testee: impl Limiter, min: f32, max: f32, δ: f32) {
|
||||
let ε = 1.0e-4f32;
|
||||
let margin = 1.0 / 16.0;
|
||||
let mul = 1.0 + margin;
|
||||
for x in itertools_num::linspace(0., min, 10) {
|
||||
assert_abs_diff_eq!(testee.value(x), 0., epsilon = ε);
|
||||
assert_abs_diff_eq!(testee.value(-x), 0., epsilon = ε);
|
||||
}
|
||||
for x in itertools_num::linspace(max, mul * max, 10) {
|
||||
assert_abs_diff_eq!(testee.value(x), 1., epsilon = ε);
|
||||
assert_abs_diff_eq!(testee.value(-x), 1., epsilon = ε);
|
||||
}
|
||||
for x in itertools_num::linspace(-mul * max, mul * max, 100) {
|
||||
let df_num = (testee.value(x + δ) - testee.value(x - δ)) / (2. * δ);
|
||||
let df_expl = testee.derivative(x);
|
||||
assert!(abs_diff_eq!(df_expl, df_num, epsilon = ε), "At x={x}, df/dx:\nnumerical: {df_num}\nexplicit: {df_expl}\n");
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_smoothstep_limiter() {
|
||||
test_limiter(SmoothstepLimiter { min: 20.0, max: 30.0 }, 20.0, 30.0, 1.0 / 32.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_smootherstep_limiter() {
|
||||
test_limiter(SmootherstepLimiter { min: 20.0, max: 30.0 }, 20.0, 30.0, 1.0 / 32.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_quadratic_accelerator() {
|
||||
let testee = super::QuadraticAccelerator { internal: 100.0, external: 150.0 };
|
||||
let ε = 1.0e-4f32;
|
||||
let δ = 1.0 / 8.0; // Mathematically, you want this to be small. Computationally, you don’t.
|
||||
let margin = 1.0 / 16.0;
|
||||
let mul = 1.0 + margin;
|
||||
assert_abs_diff_eq!(testee.u(testee.external), testee.internal, epsilon = ε);
|
||||
assert_abs_diff_eq!(testee.u(-testee.external), -testee.internal, epsilon = ε);
|
||||
assert_abs_diff_eq!(testee.du(testee.external), 1., epsilon = ε);
|
||||
assert_abs_diff_eq!(testee.du(-testee.external), 1., epsilon = ε);
|
||||
for x in itertools_num::linspace(-mul * testee.external, mul * testee.external, 100) {
|
||||
let ux = testee.u(x);
|
||||
let xux = testee.x(ux);
|
||||
assert!(abs_diff_eq!(x, xux, epsilon = ε), "At x={x}:\nu(x): {ux}\nx(u(x)): {xux}\n");
|
||||
|
||||
let du_num = (testee.u(x + δ) - testee.u(x - δ)) / (2. * δ);
|
||||
let du_expl = testee.du(x);
|
||||
assert!(abs_diff_eq!(du_expl, du_num, epsilon = ε), "At x={x}, du/dx:\nnumerical: {du_num}\nexplicit: {du_expl}\n");
|
||||
|
||||
let dudx = du_expl * testee.dx(ux);
|
||||
assert!(abs_diff_eq!(dudx, 1.0, epsilon = ε), "At x={x}:\ndu/dx * dx/du: {dudx}\n");
|
||||
|
||||
let d2u_num = (testee.du(x + δ) - testee.du(x - δ)) / (2. * δ);
|
||||
let d2u_expl = testee.d2u(x);
|
||||
assert!(abs_diff_eq!(d2u_expl, d2u_num, epsilon = ε), "At x={x}, d^2u/dx^2:\nnumerical: {d2u_num}\nexplicit: {d2u_expl}\n");
|
||||
}
|
||||
for u in itertools_num::linspace(-mul * testee.internal, mul * testee.internal, 100) {
|
||||
let xu = testee.x(u);
|
||||
let uxu = testee.u(xu);
|
||||
assert!(abs_diff_eq!(u, uxu, epsilon = ε), "At u={u}:\nx(u): {xu}\nu(x(u)): {uxu}\n");
|
||||
|
||||
let dx_num = (testee.x(u + δ) - testee.x(u - δ)) / (2. * δ);
|
||||
let dx_expl = testee.dx(u);
|
||||
assert!(abs_diff_eq!(dx_expl, dx_num, epsilon = ε), "At u={u}, dx/du:\nnumerical: {dx_num}\nexplicit: {dx_expl}\n");
|
||||
|
||||
let dudx = testee.du(xu) * dx_expl;
|
||||
assert!(abs_diff_eq!(dudx, 1.0, epsilon = ε), "At u={u}:\ndu/dx * dx/du: {dudx}\n");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -1,252 +0,0 @@
|
|||
use std::f32::consts::{FRAC_PI_2, PI};
|
||||
|
||||
use flo_canvas::*;
|
||||
use flo_draw::*;
|
||||
use glam::*;
|
||||
|
||||
use riemann::{Metric, trace_iter};
|
||||
use tube::metric::Tube;
|
||||
use tube::Space;
|
||||
use tube::Subspace::{Boundary, Inner, Outer};
|
||||
use types::{Location, Object, Ray};
|
||||
|
||||
mod riemann;
|
||||
mod fns;
|
||||
mod float_fun;
|
||||
mod tube;
|
||||
mod types;
|
||||
|
||||
const DT: f32 = 0.1;
|
||||
|
||||
fn draw_loop(gc: &mut Vec<Draw>, mut pts: impl Iterator<Item=Vec2>) {
|
||||
gc.new_path();
|
||||
let Some(first) = pts.next() else { return; };
|
||||
gc.move_to(first.x, first.y);
|
||||
for pt in pts {
|
||||
gc.line_to(pt.x, pt.y);
|
||||
}
|
||||
gc.close_path();
|
||||
gc.stroke();
|
||||
}
|
||||
|
||||
pub fn main() {
|
||||
with_2d_graphics(move || {
|
||||
let canvas = create_drawing_window("Refraction");
|
||||
canvas.draw(|gc| {
|
||||
let tube = Tube {
|
||||
inner_radius: 30.0,
|
||||
outer_radius: 50.0,
|
||||
internal_halflength: 100.0,
|
||||
external_halflength: 300.0,
|
||||
};
|
||||
|
||||
let objs: Vec<_> = [-1.25, -1.00, -0.85, -0.50, 0.00, 0.40, 0.70, 0.95, 1.05]
|
||||
.iter()
|
||||
.enumerate()
|
||||
.map(|(k, &y)| Object {
|
||||
id: k as i32,
|
||||
loc: {
|
||||
let pos = vec2(0.0, y * tube.external_halflength);
|
||||
let adj: Mat2 = tube.sqrt_at(pos).inverse().into();
|
||||
let rot = Mat2::from_angle(y);
|
||||
Location {
|
||||
pos,
|
||||
rot: adj * rot,
|
||||
}
|
||||
},
|
||||
r: 20.0,
|
||||
})
|
||||
.collect();
|
||||
let space = Space { tube, objs };
|
||||
|
||||
gc.canvas_height(500.0);
|
||||
gc.transform(Transform2D::rotate(FRAC_PI_2));
|
||||
tube.render(gc);
|
||||
gc.line_width(0.5);
|
||||
// gc.stroke_color(Color::Rgba(1.0, 0.5, 0.0, 0.5));
|
||||
// draw_fan(gc, &tube, vec2(-500.0, 0.0), vec2(1.0, 0.0), 1.0);
|
||||
gc.stroke_color(Color::Rgba(1.0, 0.5, 0.0, 1.0));
|
||||
draw_fan_2(gc, &space, vec2(-500.0, 0.0), vec2(1.0, 0.0), 1.0);
|
||||
gc.stroke_color(Color::Rgba(0.5, 1.0, 0.0, 1.0));
|
||||
draw_fan_2(gc, &space, vec2(-2.5 * tube.outer_radius, 1.25 * tube.external_halflength), vec2(1.0, -1.0), 1.0);
|
||||
draw_track(gc, &space, vec2(-500.0, 0.0), vec2(1.0, 0.2));
|
||||
draw_track(gc, &space, vec2(-500.0, 0.0), vec2(1.0, 0.5));
|
||||
draw_track(gc, &space, vec2(-0.5 * tube.inner_radius, -1.25 * tube.external_halflength), vec2(0.1, 1.0));
|
||||
|
||||
let circle_segments = 47;
|
||||
for obj in &space.objs {
|
||||
let pos = obj.loc.pos;
|
||||
gc.new_path();
|
||||
gc.circle(pos.x, pos.y, 5.0);
|
||||
gc.fill_color(Color::Rgba(0.0, 0.5, 1.0, 1.0));
|
||||
gc.fill();
|
||||
|
||||
gc.stroke_color(Color::Rgba(0.0, 0.0, 0.0, 0.5));
|
||||
draw_loop(gc, itertools_num::linspace(0.0, 2.0 * PI, circle_segments).skip(1).map(|φ| {
|
||||
let dir = Vec2::from_angle(φ) * obj.r;
|
||||
let dir = obj.loc.rot * dir;
|
||||
pos + dir
|
||||
}));
|
||||
gc.stroke_color(Color::Rgba(0.0, 0.5, 1.0, 0.5));
|
||||
draw_loop(gc, itertools_num::linspace(0.0, 2.0 * PI, circle_segments).skip(1).map(|φ| {
|
||||
let dir = Vec2::from_angle(φ) * obj.r;
|
||||
let dir = obj.loc.rot * dir;
|
||||
space.trace_step(Ray { pos, dir }).pos
|
||||
}));
|
||||
gc.stroke_color(Color::Rgba(0.5, 0.0, 1.0, 1.0));
|
||||
draw_loop(gc, itertools_num::linspace(0.0, 2.0 * PI, circle_segments).skip(1).map(|φ| {
|
||||
let n = obj.r.floor();
|
||||
let d = obj.r / n;
|
||||
let dir = Vec2::from_angle(φ);
|
||||
let dir = obj.loc.rot * dir * d;
|
||||
space.trace_iter(Ray { pos, dir }).nth(n as usize).unwrap().pos
|
||||
}));
|
||||
}
|
||||
});
|
||||
});
|
||||
}
|
||||
|
||||
fn draw_ray_2(gc: &mut Vec<Draw>, space: &Space, base: Vec2, dir: Vec2) {
|
||||
let mut hits = Vec::<Draw>::new();
|
||||
let dir = space.tube.globalize(base, dir);
|
||||
gc.new_path();
|
||||
gc.move_to(base.x, base.y);
|
||||
let mut ray = Ray { pos: base, dir: space.tube.normalize_vec_at(base, dir) * DT };
|
||||
for _ in 0..10000 {
|
||||
ray = space.trace_step(ray);
|
||||
gc.line_to(ray.pos.x, ray.pos.y);
|
||||
if ray.pos.abs().cmpgt(Vec2::splat(1000.0)).any() {
|
||||
break;
|
||||
}
|
||||
let sub = space.which_subspace(ray.pos);
|
||||
if sub == Boundary {
|
||||
continue;
|
||||
}
|
||||
gc.stroke();
|
||||
gc.new_dash_pattern();
|
||||
// gc.dash_length(6.0);
|
||||
gc.new_path();
|
||||
gc.move_to(ray.pos.x, ray.pos.y);
|
||||
let ret = match sub {
|
||||
Inner => space.trace_inner(ray),
|
||||
Outer => space.trace_outer(ray),
|
||||
Boundary => panic!(),
|
||||
};
|
||||
for hit in ret.objects {
|
||||
let obj = space.objs[hit.id as usize];
|
||||
hits.move_to(obj.loc.pos.x, obj.loc.pos.y);
|
||||
for pt in trace_iter(&space.tube, obj.loc.pos, obj.loc.rot * hit.rel.pos, hit.rel.pos.length() / 100.0).take(100) {
|
||||
hits.line_to(pt.x, pt.y);
|
||||
}
|
||||
hits.circle(hit.pos.x, hit.pos.y, 1.5);
|
||||
let Ray { pos: rel, dir } = hit.rel;
|
||||
let diff = rel.dot(dir).powi(2) - dir.length_squared() * (rel.length_squared() - obj.r.powi(2));
|
||||
assert!(diff >= 0.0);
|
||||
let t = (-rel.dot(dir) + diff.sqrt()) / dir.length_squared();
|
||||
let rel2 = hit.rel.forward(t).pos;
|
||||
let pos2 = trace_iter(&space.tube, obj.loc.pos, obj.loc.rot * rel2, rel2.length() / 100.0).nth(100).unwrap();
|
||||
hits.move_to(pos2.x - 1.0, pos2.y - 1.0);
|
||||
hits.line_to(pos2.x + 1.0, pos2.y + 1.0);
|
||||
hits.move_to(pos2.x - 1.0, pos2.y + 1.0);
|
||||
hits.line_to(pos2.x + 1.0, pos2.y - 1.0);
|
||||
}
|
||||
let a = ray.pos;
|
||||
ray = match ret.end {
|
||||
Some(r) => r,
|
||||
None => {
|
||||
ray = ray.forward(1000.0 / DT);
|
||||
gc.line_to(ray.pos.x, ray.pos.y);
|
||||
break;
|
||||
}
|
||||
};
|
||||
for p in space.line(a, ray.pos, 10.0) {
|
||||
gc.line_to(p.x, p.y);
|
||||
}
|
||||
gc.stroke();
|
||||
gc.new_dash_pattern();
|
||||
gc.new_path();
|
||||
gc.move_to(ray.pos.x, ray.pos.y);
|
||||
}
|
||||
gc.stroke();
|
||||
gc.new_path();
|
||||
gc.new_dash_pattern();
|
||||
gc.append(&mut hits);
|
||||
gc.stroke();
|
||||
}
|
||||
|
||||
fn draw_fan_2(gc: &mut Vec<Draw>, space: &Space, base: Vec2, dir: Vec2, spread: f32) {
|
||||
let dir = dir.normalize();
|
||||
let v = vec2(-dir.y, dir.x);
|
||||
for y in itertools_num::linspace(-spread, spread, 101) {
|
||||
draw_ray_2(gc, space, base, dir + y * v);
|
||||
}
|
||||
}
|
||||
|
||||
fn draw_ray(gc: &mut Vec<Draw>, space: &impl Metric, base: Vec2, dir: Vec2) {
|
||||
let dir = space.globalize(base, dir);
|
||||
gc.new_path();
|
||||
gc.move_to(base.x, base.y);
|
||||
for pt in trace_iter(space, base, dir, DT).take(10000) {
|
||||
gc.line_to(pt.x, pt.y);
|
||||
if pt.abs().cmpgt(Vec2::splat(1000.0)).any() {
|
||||
break;
|
||||
}
|
||||
}
|
||||
gc.stroke();
|
||||
}
|
||||
|
||||
fn draw_track(gc: &mut Vec<Draw>, space: &Space, start: Vec2, dir: Vec2) {
|
||||
const SCALE: f32 = 5.0;
|
||||
const STEP: f32 = 2.0 * SCALE;
|
||||
// let mut loc = Location { pos: start, rot: Mat2::IDENTITY };
|
||||
// let dir = space.tube.globalize(start, dir);
|
||||
// let v = space.tube.normalize(start, dir);
|
||||
let mut loc = Location { pos: start, rot: mat2(dir, vec2(-dir.y, dir.x)) };
|
||||
let v = vec2(1.0, 0.0);
|
||||
let mut draw = |loc: &Location| {
|
||||
let p = loc.pos;
|
||||
let ax = p + loc.rot.x_axis * SCALE;
|
||||
let ay = p + loc.rot.y_axis * SCALE;
|
||||
gc.new_path();
|
||||
gc.stroke_color(Color::Rgba(0.7, 0.0, 0.0, 1.0));
|
||||
gc.move_to(p.x, p.y);
|
||||
gc.line_to(ax.x, ax.y);
|
||||
gc.stroke();
|
||||
gc.new_path();
|
||||
gc.stroke_color(Color::Rgba(0.0, 0.7, 0.0, 1.0));
|
||||
gc.move_to(p.x, p.y);
|
||||
gc.line_to(ay.x, ay.y);
|
||||
gc.stroke();
|
||||
};
|
||||
draw(&loc);
|
||||
for _ in 0..1000 {
|
||||
let N = (STEP / DT).floor() as i32;
|
||||
for _ in 0..N {
|
||||
loc = space.move_step(loc, v * DT);
|
||||
}
|
||||
draw(&loc);
|
||||
}
|
||||
}
|
||||
|
||||
fn draw_fan(gc: &mut Vec<Draw>, space: &impl Metric, base: Vec2, dir: Vec2, spread: f32) {
|
||||
let dir = dir.normalize();
|
||||
let v = vec2(-dir.y, dir.x);
|
||||
for y in itertools_num::linspace(-spread, spread, 101) {
|
||||
draw_ray(gc, space, base, dir + y * v);
|
||||
}
|
||||
}
|
||||
|
||||
trait Renderable {
|
||||
fn render(&self, gc: &mut Vec<Draw>);
|
||||
}
|
||||
|
||||
impl Renderable for Tube {
|
||||
fn render(&self, gc: &mut Vec<Draw>) {
|
||||
gc.new_path();
|
||||
gc.rect(-self.outer_radius, -self.external_halflength, self.outer_radius, self.external_halflength);
|
||||
gc.rect(-self.inner_radius, -self.external_halflength, self.inner_radius, self.external_halflength);
|
||||
gc.winding_rule(WindingRule::EvenOdd);
|
||||
gc.fill_color(Color::Rgba(0.8, 0.8, 0.8, 1.0));
|
||||
gc.fill();
|
||||
}
|
||||
}
|
||||
|
|
@ -1,195 +0,0 @@
|
|||
use glam::*;
|
||||
|
||||
#[derive(Copy, Clone, Debug, PartialEq)]
|
||||
pub struct Decomp2 {
|
||||
pub ortho: Mat2,
|
||||
pub diag: Vec2,
|
||||
}
|
||||
|
||||
impl Decomp2 {
|
||||
fn square(&self) -> Self {
|
||||
Self {
|
||||
ortho: self.ortho,
|
||||
diag: self.diag * self.diag,
|
||||
}
|
||||
}
|
||||
|
||||
pub fn inverse(&self) -> Self {
|
||||
Self {
|
||||
ortho: self.ortho,
|
||||
diag: Vec2::splat(1.0) / self.diag,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl From<Decomp2> for Mat2 {
|
||||
fn from(value: Decomp2) -> Self {
|
||||
value.ortho.transpose() * Mat2::from_diagonal(value.diag) * value.ortho
|
||||
}
|
||||
}
|
||||
|
||||
pub type Tens2 = [Mat2; 2];
|
||||
|
||||
pub trait Metric {
|
||||
fn sqrt_at(&self, pos: Vec2) -> Decomp2;
|
||||
|
||||
fn at(&self, pos: Vec2) -> Mat2 {
|
||||
self.sqrt_at(pos).square().into()
|
||||
}
|
||||
|
||||
fn inverse_at(&self, pos: Vec2) -> Mat2 {
|
||||
self.sqrt_at(pos).square().inverse().into()
|
||||
}
|
||||
|
||||
fn part_derivs_at(&self, pos: Vec2) -> Tens2 {
|
||||
part_deriv(|p| self.at(p), pos, 1.0 / 1024.0) // division by such eps is exact which is good for overall precision
|
||||
}
|
||||
|
||||
fn vec_length_at(&self, at: Vec2, v: Vec2) -> f32 {
|
||||
v.dot(self.at(at) * v).sqrt()
|
||||
}
|
||||
|
||||
fn normalize_vec_at(&self, at: Vec2, v: Vec2) -> Vec2 {
|
||||
v / self.vec_length_at(at, v)
|
||||
}
|
||||
|
||||
fn globalize(&self, at: Vec2, v: Vec2) -> Vec2 {
|
||||
Mat2::from(self.sqrt_at(at).inverse()) * v
|
||||
}
|
||||
}
|
||||
|
||||
pub struct TraceIter<'a, M: Metric> {
|
||||
space: &'a M,
|
||||
p: Vec2,
|
||||
v: Vec2,
|
||||
dt: f32,
|
||||
}
|
||||
|
||||
impl<'a, M: Metric> Iterator for TraceIter<'a, M> {
|
||||
type Item = Vec2;
|
||||
|
||||
fn next(&mut self) -> Option<Self::Item> {
|
||||
let a: Vec2 = -contract2(krist(self.space, self.p), self.v);
|
||||
self.v = self.v + a * self.dt;
|
||||
self.p = self.p + self.v * self.dt;
|
||||
Some(self.p)
|
||||
}
|
||||
}
|
||||
|
||||
pub fn trace_iter<M: Metric>(space: &M, base: Vec2, dir: Vec2, dt: f32) -> TraceIter<M> {
|
||||
TraceIter {
|
||||
space,
|
||||
p: base,
|
||||
v: space.normalize_vec_at(base, dir),
|
||||
dt,
|
||||
}
|
||||
}
|
||||
|
||||
pub fn krist(space: &impl Metric, pos: Vec2) -> Tens2 {
|
||||
// Γ^i_k_l = .5 * g^i^m * (g_m_k,l + g_m_l,k - g_k_l,m)
|
||||
let g = &space.inverse_at(pos); // с верхними индексами
|
||||
let d = space.part_derivs_at(pos);
|
||||
// ret[i][l][k] = sum((m) => .5f * g[m][i] * (d[k][l][m] + d[l][k][m] - d[m][k][l]))
|
||||
make_tens2(|i, l, k| 0.5 * (0..2).map(|m| g.col(m)[i] * (d[l].col(k)[m] + d[k].col(m)[l] - d[m].col(k)[l])).sum::<f32>())
|
||||
}
|
||||
|
||||
fn dir_deriv(f: impl Fn(Vec2) -> Mat2, pos: Vec2, delta: Vec2) -> Mat2 {
|
||||
(f(pos + delta) - f(pos - delta)) / (2.0 * delta.length())
|
||||
}
|
||||
|
||||
fn part_deriv(f: impl Fn(Vec2) -> Mat2, pos: Vec2, eps: f32) -> Tens2 {
|
||||
[
|
||||
dir_deriv(&f, pos, vec2(eps, 0.0)),
|
||||
dir_deriv(&f, pos, vec2(0.0, eps)),
|
||||
]
|
||||
}
|
||||
|
||||
/// Сворачивает тензор t с вектором u
|
||||
pub fn contract(t: Tens2, u: Vec2) -> Mat2 {
|
||||
mat2(t[0] * u, t[1] * u).transpose()
|
||||
}
|
||||
|
||||
/// Сворачивает тензор t с вектором v дважды, по второму и третьему индексам.
|
||||
pub fn contract2(t: Tens2, v: Vec2) -> Vec2 {
|
||||
contract(t, v) * v
|
||||
}
|
||||
|
||||
fn make_vec2(f: impl Fn(usize) -> f32) -> Vec2 {
|
||||
Vec2::from_array(std::array::from_fn(|i| f(i)))
|
||||
}
|
||||
|
||||
fn make_mat2(f: impl Fn(usize, usize) -> f32) -> Mat2 {
|
||||
Mat2::from_cols_array_2d(&std::array::from_fn(|i| std::array::from_fn(|j| f(i, j))))
|
||||
}
|
||||
|
||||
fn make_tens2(f: impl Fn(usize, usize, usize) -> f32) -> Tens2 {
|
||||
std::array::from_fn(|i| make_mat2(|j, k| f(i, j, k)))
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn m2() {
|
||||
let m = make_mat2(|i, j| (i + 2 * j) as f32);
|
||||
assert_eq!(m.col(0)[0], 0.0);
|
||||
assert_eq!(m.col(1)[0], 1.0);
|
||||
assert_eq!(m.col(0)[1], 2.0);
|
||||
assert_eq!(m.col(1)[1], 3.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn t2() {
|
||||
let t = make_tens2(|i, j, k| (i + 2 * j + 4 * k) as f32);
|
||||
assert_eq!(t[0].col(0)[0], 0.0);
|
||||
assert_eq!(t[1].col(0)[0], 1.0);
|
||||
assert_eq!(t[0].col(1)[0], 2.0);
|
||||
assert_eq!(t[1].col(1)[0], 3.0);
|
||||
assert_eq!(t[0].col(0)[1], 4.0);
|
||||
assert_eq!(t[1].col(0)[1], 5.0);
|
||||
assert_eq!(t[0].col(1)[1], 6.0);
|
||||
assert_eq!(t[1].col(1)[1], 7.0);
|
||||
}
|
||||
|
||||
pub mod samples {
|
||||
use glam::{Mat2, Vec2};
|
||||
|
||||
use super::{Decomp2, Metric};
|
||||
|
||||
pub struct ScaledMetric {
|
||||
pub scale: Vec2,
|
||||
}
|
||||
|
||||
impl Metric for ScaledMetric {
|
||||
fn sqrt_at(&self, _pos: Vec2) -> Decomp2 {
|
||||
Decomp2 {
|
||||
diag: self.scale,
|
||||
ortho: Mat2::IDENTITY,
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use glam::{Mat2, mat2, vec2};
|
||||
use rand::{Rng, SeedableRng};
|
||||
|
||||
use super::{Decomp2, Metric, samples};
|
||||
|
||||
#[test]
|
||||
fn uniform_scaled_metric() {
|
||||
let mut rng = rand_pcg::Pcg64Mcg::seed_from_u64(17);
|
||||
let metric = samples::ScaledMetric { scale: vec2(3., 4.) };
|
||||
assert_eq!(metric.sqrt_at(rng.gen()), Decomp2 { ortho: Mat2::IDENTITY, diag: vec2(3., 4.) });
|
||||
assert_eq!(metric.at(rng.gen()), Mat2::from_cols_array(&[9., 0., 0., 16.]));
|
||||
assert_eq!(metric.inverse_at(rng.gen()), Mat2::from_cols_array(&[1. / 9., 0., 0., 1. / 16.]));
|
||||
assert_eq!(metric.part_derivs_at(rng.gen()), [Mat2::ZERO, Mat2::ZERO]);
|
||||
assert_eq!(metric.vec_length_at(rng.gen(), vec2(1., 0.)), 3.);
|
||||
assert_eq!(metric.vec_length_at(rng.gen(), vec2(0., 1.)), 4.);
|
||||
assert_eq!(metric.vec_length_at(rng.gen(), vec2(1., 1.)), 5.);
|
||||
assert_eq!(metric.normalize_vec_at(rng.gen(), vec2(1., 0.)), vec2(1. / 3., 0.));
|
||||
assert_eq!(metric.normalize_vec_at(rng.gen(), vec2(0., 1.)), vec2(0., 1. / 4.));
|
||||
assert_eq!(metric.normalize_vec_at(rng.gen(), vec2(1., 1.)), vec2(1. / 5., 1. / 5.));
|
||||
assert_eq!(metric.globalize(rng.gen(), vec2(1., 0.)), vec2(1. / 3., 0.));
|
||||
assert_eq!(metric.globalize(rng.gen(), vec2(0., 1.)), vec2(0., 1. / 4.));
|
||||
assert_eq!(metric.globalize(rng.gen(), vec2(1., 1.)), vec2(1. / 3., 1. / 4.));
|
||||
}
|
||||
}
|
||||
|
|
@ -1,248 +0,0 @@
|
|||
use glam::{Mat2, Vec2, vec2};
|
||||
|
||||
use crate::riemann::Metric;
|
||||
use crate::types::{Location, Ray};
|
||||
|
||||
use super::{Rect, Tube};
|
||||
|
||||
pub trait FlatCoordinateSystem<T> {
|
||||
fn flat_to_global(&self, v: T) -> T;
|
||||
fn global_to_flat(&self, v: T) -> T;
|
||||
}
|
||||
|
||||
pub trait FlatRegion: FlatCoordinateSystem<Vec2> + FlatCoordinateSystem<Ray> + FlatCoordinateSystem<Location> {
|
||||
// Измеряет расстояние до выхода за пределы области вдоль луча ray. Луч задаётся в плоской СК.
|
||||
fn distance_to_boundary(&self, _ray: Ray) -> Option<f32> { None }
|
||||
}
|
||||
|
||||
trait MetricCS: FlatCoordinateSystem<Vec2> {
|
||||
fn global_metric(&self) -> &impl Metric;
|
||||
fn flat_to_global_tfm_at(&self, pos: Vec2) -> Mat2 {
|
||||
self.global_metric().sqrt_at(self.flat_to_global(pos)).inverse().into()
|
||||
}
|
||||
fn global_to_flat_tfm_at(&self, pos: Vec2) -> Mat2 {
|
||||
self.global_metric().sqrt_at(pos).into()
|
||||
}
|
||||
}
|
||||
|
||||
impl<T: FlatCoordinateSystem<Vec2> + MetricCS> FlatCoordinateSystem<Ray> for T {
|
||||
fn flat_to_global(&self, ray: Ray) -> Ray {
|
||||
Ray {
|
||||
pos: self.flat_to_global(ray.pos),
|
||||
dir: self.flat_to_global_tfm_at(ray.pos) * ray.dir,
|
||||
}
|
||||
}
|
||||
|
||||
fn global_to_flat(&self, ray: Ray) -> Ray {
|
||||
Ray {
|
||||
pos: self.global_to_flat(ray.pos),
|
||||
dir: self.global_to_flat_tfm_at(ray.pos) * ray.dir,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl<T: FlatCoordinateSystem<Vec2> + MetricCS> FlatCoordinateSystem<Location> for T {
|
||||
fn flat_to_global(&self, loc: Location) -> Location {
|
||||
Location {
|
||||
pos: self.flat_to_global(loc.pos),
|
||||
rot: self.flat_to_global_tfm_at(loc.pos) * loc.rot,
|
||||
}
|
||||
}
|
||||
|
||||
fn global_to_flat(&self, loc: Location) -> Location {
|
||||
Location {
|
||||
pos: self.global_to_flat(loc.pos), // в плоской СК для Inner или её продолжении на Outer
|
||||
rot: self.global_to_flat_tfm_at(loc.pos) * loc.rot,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
pub struct InnerCS(pub Tube);
|
||||
|
||||
impl MetricCS for InnerCS { fn global_metric(&self) -> &impl Metric { &self.0 } }
|
||||
|
||||
impl FlatCoordinateSystem<Vec2> for InnerCS {
|
||||
fn flat_to_global(&self, pos: Vec2) -> Vec2 {
|
||||
vec2(pos.x, self.0.y(pos.y))
|
||||
}
|
||||
|
||||
// Работает только при |pos.x| ≤ inner_radius или |pos.y| ≥ external_halflength.
|
||||
fn global_to_flat(&self, pos: Vec2) -> Vec2 {
|
||||
vec2(pos.x, self.0.v(pos.y))
|
||||
}
|
||||
}
|
||||
|
||||
impl FlatRegion for InnerCS {
|
||||
fn distance_to_boundary(&self, ray: Ray) -> Option<f32> {
|
||||
Rect { size: vec2(self.0.inner_radius, self.0.internal_halflength) }.trace_out_of(ray)
|
||||
}
|
||||
}
|
||||
|
||||
pub struct OuterCS(pub Tube);
|
||||
|
||||
impl MetricCS for OuterCS { fn global_metric(&self) -> &impl Metric { &self.0 } }
|
||||
|
||||
impl FlatCoordinateSystem<Vec2> for OuterCS {
|
||||
fn flat_to_global(&self, pos: Vec2) -> Vec2 {
|
||||
let inner = Rect { size: vec2(self.0.inner_radius + 1.0, self.0.external_halflength) };
|
||||
if inner.is_inside(pos) {
|
||||
let Vec2 { x, y: v } = pos;
|
||||
let y = self.0.y(v - v.signum() * (self.0.external_halflength - self.0.internal_halflength));
|
||||
vec2(x, y)
|
||||
} else {
|
||||
pos
|
||||
}
|
||||
}
|
||||
|
||||
fn global_to_flat(&self, pos: Vec2) -> Vec2 {
|
||||
let inner = Rect { size: vec2(self.0.inner_radius + 1.0, self.0.external_halflength) };
|
||||
if inner.is_inside(pos) {
|
||||
let Vec2 { x: u, y } = pos; // в основной СК
|
||||
let v = self.0.v(y) + y.signum() * (self.0.external_halflength - self.0.internal_halflength);
|
||||
vec2(u, v) // в плоском продолжении СК Outer на область Inner
|
||||
} else {
|
||||
pos
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl FlatRegion for OuterCS {
|
||||
fn distance_to_boundary(&self, ray: Ray) -> Option<f32> {
|
||||
Rect { size: vec2(self.0.outer_radius, self.0.external_halflength) }.trace_into(ray)
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use approx::{AbsDiffEq, assert_abs_diff_eq};
|
||||
use glam::{Mat2, mat2, vec2, Vec2};
|
||||
use itertools_num::linspace;
|
||||
|
||||
use crate::riemann::samples;
|
||||
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn uniform_scaled_metric() {
|
||||
struct Scaled(samples::ScaledMetric, Vec2);
|
||||
impl FlatCoordinateSystem<Vec2> for Scaled {
|
||||
fn flat_to_global(&self, pos: Vec2) -> Vec2 { (pos - self.1) / self.0.scale }
|
||||
fn global_to_flat(&self, pos: Vec2) -> Vec2 { pos * self.0.scale + self.1 }
|
||||
}
|
||||
impl MetricCS for Scaled {
|
||||
fn global_metric(&self) -> &impl Metric { &self.0 }
|
||||
}
|
||||
let cs = Scaled(samples::ScaledMetric { scale: vec2(3., 4.) }, vec2(2., 3.));
|
||||
assert_eq!(cs.global_to_flat(vec2(7., 3.)), vec2(23., 15.));
|
||||
assert_eq!(cs.flat_to_global(vec2(8., 7.)), vec2(2., 1.));
|
||||
assert_eq!(cs.global_to_flat(Ray { pos: vec2(7., 3.), dir: vec2(3., 2.) }), Ray { pos: vec2(23., 15.), dir: vec2(9., 8.) });
|
||||
assert_eq!(cs.flat_to_global(Ray { pos: vec2(23., 15.), dir: vec2(9., 8.) }), Ray { pos: vec2(7., 3.), dir: vec2(3., 2.) });
|
||||
assert_eq!(cs.global_to_flat(Location { pos: vec2(2., 1.), rot: mat2(vec2(0., 1.), vec2(-1., 0.)) }), Location { pos: vec2(8., 7.), rot: mat2(vec2(0., 4.), vec2(-3., 0.)) });
|
||||
assert_eq!(cs.flat_to_global(Location { pos: vec2(2., 1.), rot: mat2(vec2(0., 1.), vec2(-1., 0.)) }), Location { pos: vec2(0., -0.5), rot: mat2(vec2(0., 0.25), vec2(-1. / 3., 0.)) });
|
||||
}
|
||||
|
||||
fn test_flat_region(region: &impl FlatRegion, range_global: (Vec2, Vec2), range_flat: (Vec2, Vec2)) {
|
||||
const ε: f32 = 1e-3;
|
||||
macro_rules! assert_eq_at {
|
||||
($at: expr, $left: expr, $right: expr) => {
|
||||
let at = $at;
|
||||
let left = $left;
|
||||
let right = $right;
|
||||
assert!(left.abs_diff_eq(right, ε), "Assertion failed at {at}:\n left: {left} = {}\n right: {right} = {}", stringify!($left), stringify!($right));
|
||||
};
|
||||
}
|
||||
fn check_range(name_a: &str, a: Vec2, range_a: (Vec2, Vec2), name_b: &str, b: Vec2, range_b: (Vec2, Vec2)) {
|
||||
assert!(b.cmpge(range_b.0 - ε).all() && b.cmple(range_b.1 + ε).all(), "Assertion failed:\nAt {name_a}: {a}, from range: {range_a:?}\nGot {name_b}: {b}, which is out of range {range_b:?}");
|
||||
// TODO sort out when to check these conditions:
|
||||
if a.x.abs_diff_eq(&range_a.0.x, ε) { assert_abs_diff_eq!(b.x, range_b.0.x, epsilon=ε); }
|
||||
if a.y.abs_diff_eq(&range_a.0.y, ε) { assert_abs_diff_eq!(b.y, range_b.0.y, epsilon=ε); }
|
||||
if a.x.abs_diff_eq(&range_a.1.x, ε) { assert_abs_diff_eq!(b.x, range_b.1.x, epsilon=ε); }
|
||||
if a.y.abs_diff_eq(&range_a.1.y, ε) { assert_abs_diff_eq!(b.y, range_b.1.y, epsilon=ε); }
|
||||
}
|
||||
for x in linspace(range_global.0.x, range_global.1.x, 20) {
|
||||
for y in linspace(range_global.0.y, range_global.1.y, 20) {
|
||||
let pos_global = vec2(x, y);
|
||||
let pos_flat = region.global_to_flat(pos_global);
|
||||
check_range("global", pos_global, range_global, "flat", pos_flat, range_flat);
|
||||
assert_eq_at!(pos_global, region.global_to_flat(Location { pos: pos_global, rot: Mat2::IDENTITY }).pos, pos_flat);
|
||||
assert_eq_at!(pos_global, region.flat_to_global(pos_flat), pos_global);
|
||||
assert_eq_at!(pos_global, region.flat_to_global(region.global_to_flat(Location { pos: pos_global, rot: Mat2::IDENTITY })).rot, Mat2::IDENTITY);
|
||||
}
|
||||
}
|
||||
for x in linspace(range_flat.0.x, range_flat.1.x, 20) {
|
||||
for y in linspace(range_flat.0.y, range_flat.1.y, 20) {
|
||||
let pos_flat = vec2(x, y);
|
||||
let pos_global = region.flat_to_global(pos_flat);
|
||||
check_range("flat", pos_flat, range_flat, "global", pos_global, range_global);
|
||||
assert_eq_at!(pos_flat, region.flat_to_global(Location { pos: pos_flat, rot: Mat2::IDENTITY }).pos, pos_global);
|
||||
assert_eq_at!(pos_flat, region.global_to_flat(pos_global), pos_flat);
|
||||
assert_eq_at!(pos_flat, region.global_to_flat(region.flat_to_global(Location { pos: pos_global, rot: Mat2::IDENTITY })).rot, Mat2::IDENTITY);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_mapper_inner() {
|
||||
let mapper = InnerCS(Tube {
|
||||
inner_radius: 30.0,
|
||||
outer_radius: 50.0,
|
||||
internal_halflength: 100.0,
|
||||
external_halflength: 300.0,
|
||||
});
|
||||
test_flat_region(&mapper, (vec2(-30.0, -300.0), vec2(30.0, 300.0)), (vec2(-30.0, -100.0), vec2(30.0, 100.0)));
|
||||
test_flat_region(&mapper, (vec2(-60.0, -400.0), vec2(60.0, -300.0)), (vec2(-60.0, -200.0), vec2(60.0, -100.0)));
|
||||
test_flat_region(&mapper, (vec2(-60.0, 300.0), vec2(60.0, 400.0)), (vec2(-60.0, 100.0), vec2(60.0, 200.0)));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_mapper_outer() {
|
||||
let mapper = OuterCS(Tube {
|
||||
inner_radius: 30.0,
|
||||
outer_radius: 50.0,
|
||||
internal_halflength: 100.0,
|
||||
external_halflength: 300.0,
|
||||
});
|
||||
// TODO replace 200.20016 with something sane
|
||||
test_flat_region(&mapper, (vec2(-30.0, -300.0), vec2(30.0, -1.0)), (vec2(-30.0, -300.0), vec2(30.0, -200.20016)));
|
||||
test_flat_region(&mapper, (vec2(-30.0, 1.0), vec2(30.0, 300.0)), (vec2(-30.0, 200.20016), vec2(30.0, 300.0)));
|
||||
test_flat_region(&mapper, (vec2(-60.0, -400.0), vec2(60.0, -300.0)), (vec2(-60.0, -400.0), vec2(60.0, -300.0)));
|
||||
test_flat_region(&mapper, (vec2(-60.0, 300.0), vec2(60.0, 400.0)), (vec2(-60.0, 300.0), vec2(60.0, 400.0)));
|
||||
// straight
|
||||
for x in linspace(-60., 60., 20) {
|
||||
for y in linspace(-320., 320., 20) {
|
||||
assert_eq!(mapper.global_to_flat(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos.x, x);
|
||||
}
|
||||
}
|
||||
// symmetrical
|
||||
for x in linspace(0., 60., 20) {
|
||||
for y in linspace(0., 320., 20) {
|
||||
let pp = mapper.global_to_flat(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos;
|
||||
let np = mapper.global_to_flat(Location { pos: vec2(-x, y), rot: Mat2::IDENTITY }).pos;
|
||||
let pn = mapper.global_to_flat(Location { pos: vec2(x, -y), rot: Mat2::IDENTITY }).pos;
|
||||
let nn = mapper.global_to_flat(Location { pos: vec2(-x, -y), rot: Mat2::IDENTITY }).pos;
|
||||
assert_eq!(np, vec2(-pp.x, pp.y));
|
||||
assert_eq!(pn, vec2(pp.x, -pp.y));
|
||||
assert_eq!(nn, vec2(-pp.x, -pp.y));
|
||||
}
|
||||
}
|
||||
// clean boundary
|
||||
for x in linspace(50., 60., 20) {
|
||||
for y in linspace(0., 320., 20) {
|
||||
assert_eq!(mapper.global_to_flat(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos.y, y);
|
||||
}
|
||||
}
|
||||
for x in linspace(0., 60., 20) {
|
||||
for y in linspace(300., 320., 20) {
|
||||
assert_eq!(mapper.global_to_flat(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos.y, y);
|
||||
}
|
||||
}
|
||||
// accelerating
|
||||
for x in linspace(-29., 29., 20) {
|
||||
for y in linspace(1., 299., 20) {
|
||||
let v = mapper.global_to_flat(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos.y;
|
||||
assert!(v > 200.0);
|
||||
assert!(v > y);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -1,204 +0,0 @@
|
|||
use glam::{f32, Mat2, Vec2, vec2};
|
||||
|
||||
use crate::fns::{self, Limiter};
|
||||
use crate::riemann::{Decomp2, Metric, Tens2};
|
||||
|
||||
#[derive(Copy, Clone, Debug)]
|
||||
pub struct Tube {
|
||||
pub outer_radius: f32,
|
||||
pub inner_radius: f32,
|
||||
pub external_halflength: f32,
|
||||
pub internal_halflength: f32,
|
||||
}
|
||||
|
||||
impl Tube {
|
||||
fn fx(&self) -> impl Limiter { fns::SmootherstepLimiter { min: self.inner_radius, max: self.outer_radius } }
|
||||
fn fy(&self) -> fns::QuadraticAccelerator { fns::QuadraticAccelerator { internal: self.internal_halflength, external: self.external_halflength } }
|
||||
|
||||
pub fn y(&self, v: f32) -> f32 { self.fy().x(v) }
|
||||
pub fn v(&self, y: f32) -> f32 { self.fy().u(y) }
|
||||
pub fn dy(&self, v: f32) -> f32 { self.fy().dx(v) }
|
||||
pub fn dv(&self, y: f32) -> f32 { self.fy().du(y) }
|
||||
}
|
||||
|
||||
impl Metric for Tube {
|
||||
fn sqrt_at(&self, pos: Vec2) -> Decomp2 {
|
||||
let sx = self.fx().value(pos.x);
|
||||
let sy = self.fy().du(pos.y);
|
||||
let s = sx + sy - sx * sy;
|
||||
assert!(sx.is_finite());
|
||||
assert!(sy.is_finite());
|
||||
assert!(sy > 0.0);
|
||||
Decomp2 {
|
||||
ortho: Mat2::IDENTITY,
|
||||
diag: vec2(1.0, s),
|
||||
}
|
||||
}
|
||||
|
||||
fn part_derivs_at(&self, pos: Vec2) -> Tens2 {
|
||||
let sx = self.fx().value(pos.x);
|
||||
let sy = self.fy().du(pos.y);
|
||||
let s = sx + sy - sx * sy;
|
||||
let dsx_dx = self.fx().derivative(pos.x);
|
||||
let dsy_dy = self.fy().d2u(pos.y);
|
||||
let ds2_dx = 2.0 * s * (1.0 - sy) * dsx_dx;
|
||||
let ds2_dy = 2.0 * s * (1.0 - sx) * dsy_dy;
|
||||
[
|
||||
Mat2::from_cols_array(&[0.0, 0.0, 0.0, ds2_dx]),
|
||||
Mat2::from_cols_array(&[0.0, 0.0, 0.0, ds2_dy]),
|
||||
]
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod test {
|
||||
use approx::assert_abs_diff_eq;
|
||||
use glam::{Vec2, vec2};
|
||||
use itertools_num::linspace;
|
||||
|
||||
use crate::riemann::{Decomp2, Metric};
|
||||
use crate::tube::Space;
|
||||
use crate::types::Ray;
|
||||
|
||||
use super::Tube;
|
||||
|
||||
#[test]
|
||||
fn test_tube_metric_derivs() {
|
||||
struct Approx(Tube);
|
||||
impl Metric for Approx {
|
||||
fn sqrt_at(&self, pos: Vec2) -> Decomp2 { self.0.sqrt_at(pos) }
|
||||
}
|
||||
let testee = Tube {
|
||||
inner_radius: 30.0,
|
||||
outer_radius: 50.0,
|
||||
internal_halflength: 100.0,
|
||||
external_halflength: 300.0,
|
||||
};
|
||||
let approx = Approx(testee);
|
||||
let epsilon = 1.0e-3;
|
||||
let margin = 1.0 / 16.0;
|
||||
let mul = 1.0 + margin;
|
||||
for x in itertools_num::linspace(-mul * testee.outer_radius, mul * testee.outer_radius, 100) {
|
||||
for y in itertools_num::linspace(-mul * testee.external_halflength, mul * testee.external_halflength, 100) {
|
||||
let pos = vec2(x, y);
|
||||
let computed = testee.part_derivs_at(pos);
|
||||
let reference = approx.part_derivs_at(pos);
|
||||
let eq = (0..2).all(|coord| computed[coord].abs_diff_eq(reference[coord], epsilon));
|
||||
assert!(eq, "Bad derivative computation at {pos}:\n explicit: {computed:?}\n numerical: {reference:?}\n");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_accelerator() {
|
||||
let space = Space {
|
||||
tube: Tube {
|
||||
inner_radius: 30.0,
|
||||
outer_radius: 50.0,
|
||||
internal_halflength: 100.0,
|
||||
external_halflength: 300.0,
|
||||
},
|
||||
objs: vec![],
|
||||
};
|
||||
let ε = 1e-3;
|
||||
let off = 10.0;
|
||||
let steps = 1024;
|
||||
for ax in [-30.0 + ε, -25.0, -3.0, 17.0, 30.0 - ε] {
|
||||
for bx in [0.0, ε, 1.0, 7.0, 30.0 - ε] {
|
||||
let a = vec2(ax, -(space.tube.external_halflength + off));
|
||||
let b = vec2(bx, space.tube.external_halflength + off);
|
||||
let Δ = vec2(bx - ax, 2.0 * (space.tube.internal_halflength + off));
|
||||
let dir = Δ / (steps as f32);
|
||||
let traced = space.trace_iter(Ray { pos: a, dir }).nth(steps).unwrap();
|
||||
assert_abs_diff_eq!(traced.pos.x, b.x, epsilon=1.0e-2);
|
||||
assert_abs_diff_eq!(traced.pos.y, b.y, epsilon=1.0e1);
|
||||
assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon=1.0e-3);
|
||||
assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon=1.0e-2);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
#[ignore]
|
||||
fn test_accelerator_slow() {
|
||||
let space = Space {
|
||||
tube: Tube {
|
||||
inner_radius: 30.0,
|
||||
outer_radius: 50.0,
|
||||
internal_halflength: 100.0,
|
||||
external_halflength: 300.0,
|
||||
},
|
||||
objs: vec![],
|
||||
};
|
||||
let ε = 1e-3;
|
||||
let off = 10.0;
|
||||
let steps = 4096;
|
||||
for ax in linspace(-space.tube.inner_radius + ε, space.tube.inner_radius - ε, 20) {
|
||||
for bx in linspace(-space.tube.inner_radius + ε, space.tube.inner_radius - ε, 20) {
|
||||
let a = vec2(ax, -(space.tube.external_halflength + off));
|
||||
let b = vec2(bx, space.tube.external_halflength + off);
|
||||
let Δ = vec2(bx - ax, 2.0 * (space.tube.internal_halflength + off));
|
||||
let dir = Δ / (steps as f32);
|
||||
let traced = space.trace_iter(Ray { pos: a, dir }).nth(steps).unwrap();
|
||||
assert_abs_diff_eq!(traced.pos.x, b.x, epsilon=1.0e-2);
|
||||
assert_abs_diff_eq!(traced.pos.y, b.y, epsilon=1.0e0);
|
||||
assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon=1.0e-3);
|
||||
assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon=1.0e-3);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_accelerator_inner_edge() {
|
||||
let space = Space {
|
||||
tube: Tube {
|
||||
inner_radius: 30.0,
|
||||
outer_radius: 50.0,
|
||||
internal_halflength: 100.0,
|
||||
external_halflength: 300.0,
|
||||
},
|
||||
objs: vec![],
|
||||
};
|
||||
let ε = 1e-3;
|
||||
let off = 10.0;
|
||||
let steps = 10000;
|
||||
for x in [space.tube.inner_radius - ε, space.tube.inner_radius + ε] {
|
||||
let a = vec2(x, -(space.tube.external_halflength + off));
|
||||
let b = vec2(x, space.tube.external_halflength + off);
|
||||
let Δ = vec2(0.0, 2.0 * (space.tube.internal_halflength + off));
|
||||
let dir = Δ / (steps as f32);
|
||||
let traced = space.trace_iter(Ray { pos: a, dir }).nth(steps).unwrap();
|
||||
assert_abs_diff_eq!(traced.pos.x, b.x, epsilon=1.0e-1);
|
||||
assert_abs_diff_eq!(traced.pos.y, b.y, epsilon=1.0e0);
|
||||
assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon=1.0e-2);
|
||||
assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon=1.0e-2);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_accelerator_outer_edge() {
|
||||
let space = Space {
|
||||
tube: Tube {
|
||||
inner_radius: 30.0,
|
||||
outer_radius: 50.0,
|
||||
internal_halflength: 100.0,
|
||||
external_halflength: 300.0,
|
||||
},
|
||||
objs: vec![],
|
||||
};
|
||||
let ε = 1e-3;
|
||||
let off = 10.0;
|
||||
let steps = 4096;
|
||||
for x in [space.tube.outer_radius + ε, space.tube.outer_radius - ε] {
|
||||
let a = vec2(x, -(space.tube.external_halflength + off));
|
||||
let b = vec2(x, space.tube.external_halflength + off);
|
||||
let Δ = vec2(0.0, 2.0 * (space.tube.external_halflength + off));
|
||||
let dir = Δ / (steps as f32);
|
||||
let traced = space.trace_iter(Ray { pos: a, dir }).nth(steps).unwrap();
|
||||
assert_abs_diff_eq!(traced.pos.x, b.x, epsilon=2.0e0);
|
||||
assert_abs_diff_eq!(traced.pos.y, b.y, epsilon=1.0e0);
|
||||
assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon=1.0e-2);
|
||||
assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon=1.0e-2);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -1,187 +0,0 @@
|
|||
use glam::{bool, f32, Mat2, Vec2, vec2};
|
||||
|
||||
use coords::{FlatCoordinateSystem, InnerCS, OuterCS};
|
||||
use metric::Tube;
|
||||
use Subspace::{Boundary, Inner, Outer};
|
||||
|
||||
use crate::riemann;
|
||||
use crate::tube::coords::FlatRegion;
|
||||
use crate::types::{FlatTraceResult, Hit, Location, Object, Ray};
|
||||
|
||||
pub mod metric;
|
||||
mod coords;
|
||||
|
||||
pub struct Space {
|
||||
pub tube: Tube,
|
||||
pub objs: Vec<Object>,
|
||||
}
|
||||
|
||||
#[derive(PartialEq, Eq, Debug)]
|
||||
pub enum Subspace {
|
||||
Outer,
|
||||
Boundary,
|
||||
Inner,
|
||||
}
|
||||
|
||||
impl Space {
|
||||
pub fn which_subspace(&self, pt: Vec2) -> Subspace {
|
||||
if pt.y.abs() > self.tube.external_halflength {
|
||||
Outer
|
||||
} else if pt.x.abs() > self.tube.outer_radius {
|
||||
Outer
|
||||
} else if pt.x.abs() > self.tube.inner_radius {
|
||||
Boundary
|
||||
} else {
|
||||
Inner
|
||||
}
|
||||
}
|
||||
|
||||
/// Выполняет один шаг трассировки. Работает в любой части пространства, но вне Boundary доступны более эффективные методы.
|
||||
/// ray задаётся в основной СК.
|
||||
pub fn trace_step(&self, ray: Ray) -> Ray {
|
||||
let a: Vec2 = -riemann::contract2(riemann::krist(&self.tube, ray.pos), ray.dir);
|
||||
let v = ray.dir + a;
|
||||
let p = ray.pos + v;
|
||||
Ray { pos: p, dir: v }
|
||||
}
|
||||
|
||||
/// Выполняет один шаг перемещения. Работает в любой части пространства.
|
||||
/// off задаётся в локальной СК. Рекомендуется считать небольшими шагами.
|
||||
pub fn move_step(&self, loc: Location, off: Vec2) -> Location {
|
||||
let corr = Mat2::IDENTITY - riemann::contract(riemann::krist(&self.tube, loc.pos), loc.rot * off);
|
||||
let p = loc.pos + corr * loc.rot * off;
|
||||
Location { pos: p, rot: corr * loc.rot }
|
||||
}
|
||||
|
||||
pub fn trace_iter(&self, ray: Ray) -> impl Iterator<Item=Ray> + '_ {
|
||||
std::iter::successors(Some(ray), |&ray| Some(self.trace_step(ray)))
|
||||
}
|
||||
|
||||
pub fn trace_inner(&self, ray: Ray) -> FlatTraceResult {
|
||||
assert_eq!(self.which_subspace(ray.pos), Inner);
|
||||
self.trace_flat(InnerCS(self.tube), ray)
|
||||
}
|
||||
|
||||
pub fn trace_outer(&self, ray: Ray) -> FlatTraceResult {
|
||||
assert_eq!(self.which_subspace(ray.pos), Outer);
|
||||
self.trace_flat(OuterCS(self.tube), ray)
|
||||
}
|
||||
|
||||
fn trace_flat(&self, cs: impl FlatRegion, ray: Ray) -> FlatTraceResult {
|
||||
let ray = cs.global_to_flat(ray);
|
||||
let dist = cs.distance_to_boundary(ray);
|
||||
let objs = self.list_objects(|loc| cs.global_to_flat(loc));
|
||||
FlatTraceResult {
|
||||
end: dist.map(|dist| cs.flat_to_global(ray.forward(dist))),
|
||||
objects: Self::hit_objects(objs.as_slice(), ray, dist, |pos| cs.flat_to_global(pos)),
|
||||
}
|
||||
}
|
||||
|
||||
fn trace_boundary(&self, ray: Ray) -> Ray {
|
||||
assert_eq!(self.which_subspace(ray.pos), Boundary);
|
||||
self.trace_iter(ray)
|
||||
.find(|&ray| self.which_subspace(ray.pos) != Boundary)
|
||||
.expect("Can't get outta the wall!")
|
||||
}
|
||||
|
||||
fn list_objects(&self, tfm: impl Fn(Location) -> Location) -> Vec<Object> {
|
||||
self.objs.iter().map(|&Object { id, loc, r }| Object { id, loc: tfm(loc), r }).collect()
|
||||
}
|
||||
|
||||
fn hit_objects(objs: &[Object], ray: Ray, limit: Option<f32>, globalize: impl Fn(Vec2) -> Vec2) -> Vec<Hit> {
|
||||
let limit = limit.unwrap_or(f32::INFINITY);
|
||||
objs.iter()
|
||||
.filter_map(|obj| {
|
||||
let rel = ray.pos - obj.loc.pos;
|
||||
let diff = rel.dot(ray.dir).powi(2) - ray.dir.length_squared() * (rel.length_squared() - obj.r.powi(2));
|
||||
if diff > 0.0 {
|
||||
let t = (-rel.dot(ray.dir) - diff.sqrt()) / ray.dir.length_squared();
|
||||
Some((obj, t))
|
||||
} else {
|
||||
None
|
||||
}
|
||||
})
|
||||
.filter(|&(_, t)| t >= 0.0 && t < limit)
|
||||
.map(|(obj, t)| {
|
||||
let pos = ray.forward(t).pos;
|
||||
let rel = obj.loc.rot.inverse() * Ray { pos: pos - obj.loc.pos, dir: ray.dir };
|
||||
Hit { id: obj.id, distance: t, pos: globalize(pos), rel }
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
|
||||
pub fn line(&self, a: Vec2, b: Vec2, step: f32) -> Vec<Vec2> {
|
||||
match self.which_subspace(a) {
|
||||
Outer => vec![b],
|
||||
Inner => {
|
||||
let 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);
|
||||
(1..=n).map(|k| cs.flat_to_global(a.lerp(b, k as f32 / n as f32))).collect()
|
||||
}
|
||||
Boundary => panic!("Can't draw a line here!"),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
struct Rect {
|
||||
pub size: Vec2,
|
||||
}
|
||||
|
||||
impl Rect {
|
||||
/// Отражает луч, чтобы все координаты направления были положительны (допустимо благодаря симметрии Rect).
|
||||
fn flip_ray(ray: Ray) -> Ray {
|
||||
Ray { pos: ray.pos * ray.dir.signum(), dir: ray.dir.abs() }
|
||||
}
|
||||
|
||||
fn is_inside(&self, pt: Vec2) -> bool {
|
||||
pt.abs().cmplt(self.size).all()
|
||||
}
|
||||
|
||||
fn trace_into(&self, ray: Ray) -> Option<f32> {
|
||||
let ray = Self::flip_ray(ray);
|
||||
// ray.pos.x + t * ray.dir.x = −size.x
|
||||
let ts = (-self.size - ray.pos) / ray.dir;
|
||||
let t = ts.max_element();
|
||||
let pt = ray.pos + t * ray.dir;
|
||||
if t < 0.0 { return None; }
|
||||
if pt.cmpgt(self.size).any() { return None; }
|
||||
Some(t)
|
||||
}
|
||||
|
||||
fn trace_out_of(&self, ray: Ray) -> Option<f32> {
|
||||
let ray = Self::flip_ray(ray);
|
||||
// ray.pos.x + t * ray.dir.x = +size.x
|
||||
let ts = (self.size - ray.pos) / ray.dir;
|
||||
let t = ts.min_element();
|
||||
Some(t)
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_rect() {
|
||||
assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 5.0) }), Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 5.0) });
|
||||
assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(-4.0, 5.0) }), Ray { pos: vec2(-2.0, 3.0), dir: vec2(4.0, 5.0) });
|
||||
assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, -5.0) }), Ray { pos: vec2(2.0, -3.0), dir: vec2(4.0, 5.0) });
|
||||
assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 0.0) }), Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 0.0) });
|
||||
|
||||
let r = Rect { size: vec2(2.0, 3.0) };
|
||||
|
||||
assert_eq!(r.trace_into(Ray { pos: vec2(3.0, 3.0), dir: vec2(1.0, 1.0) }), None);
|
||||
assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 2.0), dir: vec2(1.0, 0.0) }), Some(1.0));
|
||||
assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 2.0), dir: vec2(-1.0, 0.0) }), None);
|
||||
assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 1.0), dir: vec2(2.0, 2.0) }), Some(0.5));
|
||||
assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 2.1), dir: vec2(2.0, 2.0) }), None);
|
||||
|
||||
assert_eq!(r.trace_into(Ray { pos: vec2(2.0, 3.0), dir: vec2(1.0, 1.0) }), None);
|
||||
assert_eq!(r.trace_into(Ray { pos: vec2(-2.0, 3.0), dir: vec2(-1.0, 1.0) }), None);
|
||||
assert_eq!(r.trace_into(Ray { pos: vec2(2.0, 3.0), dir: vec2(-1.0, -1.0) }), Some(0.0));
|
||||
assert_eq!(r.trace_into(Ray { pos: vec2(2.0, -3.0), dir: vec2(-1.0, 1.0) }), Some(0.0));
|
||||
|
||||
assert_eq!(r.trace_out_of(Ray { pos: vec2(0.0, 0.0), dir: vec2(1.0, 1.0) }), Some(2.0));
|
||||
assert_eq!(r.trace_out_of(Ray { pos: vec2(0.0, 0.0), dir: vec2(0.0, 1.0) }), Some(3.0));
|
||||
assert_eq!(r.trace_out_of(Ray { pos: vec2(0.0, 1.0), dir: vec2(0.0, -1.0) }), Some(4.0));
|
||||
assert_eq!(r.trace_out_of(Ray { pos: vec2(1.0, 1.0), dir: vec2(0.0, -1.0) }), Some(4.0));
|
||||
assert_eq!(r.trace_out_of(Ray { pos: vec2(2.0, 3.0), dir: vec2(1.0, 1.0) }), Some(0.0));
|
||||
}
|
||||
|
|
@ -1,50 +0,0 @@
|
|||
use glam::{f32, i32, Mat2, Vec2};
|
||||
|
||||
#[derive(Copy, Clone, Debug, PartialEq)]
|
||||
pub struct Ray {
|
||||
pub pos: Vec2,
|
||||
pub dir: Vec2,
|
||||
}
|
||||
|
||||
impl Ray {
|
||||
pub fn forward(&self, dist: f32) -> Ray {
|
||||
Ray { pos: self.pos + self.dir * dist, dir: self.dir }
|
||||
}
|
||||
}
|
||||
|
||||
impl std::ops::Mul<Ray> for Mat2 {
|
||||
type Output = Ray;
|
||||
|
||||
fn mul(self, rhs: Ray) -> Self::Output {
|
||||
Ray { pos: self * rhs.pos, dir: self * rhs.dir }
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Copy, Clone, Debug, PartialEq)]
|
||||
pub struct Location {
|
||||
/// Положение в основной СК
|
||||
pub pos: Vec2,
|
||||
/// Преобразование вектора из локальной ортонормированной в основную СК
|
||||
pub rot: Mat2,
|
||||
}
|
||||
|
||||
#[derive(Copy, Clone, Debug)]
|
||||
pub struct Object {
|
||||
pub id: i32,
|
||||
pub loc: Location,
|
||||
pub r: f32,
|
||||
}
|
||||
|
||||
#[derive(Copy, Clone, Debug)]
|
||||
pub struct Hit {
|
||||
pub distance: f32,
|
||||
pub id: i32,
|
||||
pub pos: Vec2, // положение в основной СК
|
||||
pub rel: Ray, // в локальной ортонормированной СК объекта
|
||||
}
|
||||
|
||||
#[derive(Clone, Debug)]
|
||||
pub struct FlatTraceResult {
|
||||
pub end: Option<Ray>,
|
||||
pub objects: Vec<Hit>,
|
||||
}
|
||||
|
|
@ -1,151 +0,0 @@
|
|||
mod mesh_loader;
|
||||
|
||||
use std::fs::File;
|
||||
use std::{env};
|
||||
use std::error::Error;
|
||||
use std::f32::consts::PI;
|
||||
use std::io::{BufReader};
|
||||
use glam::*;
|
||||
use show_image::{ImageInfo, ImageView, WindowOptions};
|
||||
use crate::mesh_loader::{Face, load_mesh};
|
||||
|
||||
const W: i32 = 320;
|
||||
const H: i32 = 240;
|
||||
|
||||
#[derive(Copy, Clone)]
|
||||
struct Color(u8, u8, u8);
|
||||
|
||||
struct Image {
|
||||
w: i32,
|
||||
h: i32,
|
||||
data: Vec<u8>,
|
||||
}
|
||||
|
||||
impl Image {
|
||||
fn data(&self) -> &[u8] {
|
||||
self.data.as_slice()
|
||||
}
|
||||
|
||||
fn put_pixel(&mut self, x: i32, y: i32, color: Color) {
|
||||
if x < 0 || x >= self.w || y < 0 || y > self.h {
|
||||
return;
|
||||
}
|
||||
let index = 3 * (x + self.w * y) as usize;
|
||||
self.data[index] = color.0;
|
||||
self.data[index + 1] = color.1;
|
||||
self.data[index + 2] = color.2;
|
||||
}
|
||||
}
|
||||
|
||||
fn ypr_to_mat(ypr: Vec3) -> Mat3 {
|
||||
let Vec3 { x: yaw, y: pitch, z: roll } = ypr;
|
||||
let m_roll = mat3(
|
||||
vec3(roll.cos(), roll.sin(), 0.0),
|
||||
vec3(-roll.sin(), roll.cos(), 0.0),
|
||||
vec3(0.0, 0.0, 1.0));
|
||||
let m_yaw = mat3(
|
||||
vec3(yaw.cos(), 0.0, yaw.sin()),
|
||||
vec3(0.0, 1.0, 0.0),
|
||||
vec3(-yaw.sin(), 0.0, yaw.cos()));
|
||||
let m_pitch = mat3(
|
||||
vec3(1.0, 0.0, 0.0),
|
||||
vec3(0.0, pitch.cos(), -pitch.sin()),
|
||||
vec3(0.0, pitch.sin(), pitch.cos()));
|
||||
m_roll * m_pitch * m_yaw
|
||||
}
|
||||
|
||||
type Mesh = [Face];
|
||||
|
||||
struct TraceResult {
|
||||
distance: f32,
|
||||
normal: Vec3,
|
||||
}
|
||||
|
||||
fn trace_to_mesh(mesh: &Mesh, base: Vec3, ray: Vec3) -> Option<TraceResult> {
|
||||
let mut ret: Option<TraceResult> = None;
|
||||
let mut dist = f32::INFINITY;
|
||||
for f in mesh {
|
||||
let fs = (0..3).map(|k| edge_dist(f.vertices[k], f.vertices[(k + 1) % 3], base, ray));
|
||||
if fs.into_iter().all(|f| f >= 0.0) {
|
||||
let m = mat3(f.vertices[1] - f.vertices[0], f.vertices[2] - f.vertices[0], -ray);
|
||||
let m = m.inverse();
|
||||
let rel = m * (base - f.vertices[0]);
|
||||
if rel.z > dist {
|
||||
continue;
|
||||
}
|
||||
dist = rel.z;
|
||||
ret = Some(TraceResult {
|
||||
distance: rel.z,
|
||||
normal: f.normal,
|
||||
});
|
||||
}
|
||||
}
|
||||
ret
|
||||
}
|
||||
|
||||
struct Location {
|
||||
pos: Vec3,
|
||||
rot: Vec4,
|
||||
}
|
||||
|
||||
fn render(mesh: &Mesh, camera: impl Fn(Vec2) -> (Vec3, Vec3)) -> Image {
|
||||
let bkg = vec3(0.0, 0.0, 0.0);
|
||||
let mut img = Image {
|
||||
w: W,
|
||||
h: H,
|
||||
data: vec![0; (3 * W * H) as usize],
|
||||
};
|
||||
let img_size = vec2(W as f32, H as f32);
|
||||
for y in 0..H {
|
||||
for x in 0..W {
|
||||
let img_coords = vec2(x as f32, y as f32);
|
||||
let off = (img_coords - img_size * 0.5) / img_size.y;
|
||||
let (base, ray) = camera(off);
|
||||
let color = if let Some(r) = trace_to_mesh(mesh, base, ray.normalize()) {
|
||||
// to_vec3(0.45) * dot(r.normal, normalize(vec3(-1.0, 1.0, -1.0))) + 0.50
|
||||
r.normal * 0.45 + 0.50
|
||||
} else {
|
||||
bkg
|
||||
};
|
||||
let color = (color * 255.0).as_ivec3().clamp(IVec3::splat(0), IVec3::splat(255));
|
||||
img.put_pixel(x, y, Color(color.x as u8, color.y as u8, color.z as u8));
|
||||
}
|
||||
}
|
||||
img
|
||||
}
|
||||
|
||||
#[show_image::main]
|
||||
fn main() -> Result<(), Box<dyn Error>> {
|
||||
let args: Vec<String> = env::args().collect();
|
||||
let mesh = {
|
||||
let f = File::open(&args[1])?;
|
||||
let mut f = BufReader::new(f);
|
||||
load_mesh(&mut f)?
|
||||
};
|
||||
let window = show_image::create_window("Raytracing", WindowOptions::default())?;
|
||||
loop {
|
||||
for phi in 0..360 {
|
||||
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| {
|
||||
// perspective projection
|
||||
let base = vec3(0.0, 0.0, -40.0);
|
||||
let ray = vec3(off.x, off.y, 2.0);
|
||||
|
||||
// orthographic projection
|
||||
// let base = vec3(off.x, off.y, -10.0);
|
||||
// let ray = vec3(0.0, 0.0, 1.0);
|
||||
|
||||
(m_camera * base, m_camera * ray)
|
||||
});
|
||||
|
||||
let image = ImageView::new(ImageInfo::rgb8(W as u32, H as u32), img.data());
|
||||
window.set_image("image", image)?;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
fn edge_dist(a: Vec3, b: Vec3, base: Vec3, dir: Vec3) -> f32 {
|
||||
// Note: given that the input is not arbitrary but comes from a cartesian product of certain (a, b) pairs and certain (base, dir) pairs, this can be optimized from Cnm to an+bm+cnm with c<C.
|
||||
mat3(b - a, base - a, -dir).determinant()
|
||||
}
|
||||
127
src/main.rs
Normal file
127
src/main.rs
Normal file
|
|
@ -0,0 +1,127 @@
|
|||
mod mesh_loader;
|
||||
|
||||
use std::fs::File;
|
||||
use std::{env};
|
||||
use std::error::Error;
|
||||
use std::f32::consts::PI;
|
||||
use std::io::{BufRead, BufReader};
|
||||
use std::io::Write;
|
||||
use rand::Rng;
|
||||
use glm::*;
|
||||
use show_image::{ImageInfo, ImageView, WindowOptions};
|
||||
use crate::mesh_loader::{Face, load_mesh};
|
||||
|
||||
const W: i32 = 320;
|
||||
const H: i32 = 240;
|
||||
const SCALE: f32 = 30.0;
|
||||
|
||||
#[derive(Copy, Clone)]
|
||||
struct Color(u8, u8, u8);
|
||||
|
||||
struct Image {
|
||||
w: i32,
|
||||
h: i32,
|
||||
data: Vec<u8>,
|
||||
}
|
||||
|
||||
impl Image {
|
||||
fn data(&self) -> &[u8] {
|
||||
self.data.as_slice()
|
||||
}
|
||||
|
||||
fn put_pixel(&mut self, x: i32, y: i32, color: Color) {
|
||||
if x < 0 || x >= self.w || y < 0 || y > self.h {
|
||||
return;
|
||||
}
|
||||
let index = 3 * (x + self.w * y) as usize;
|
||||
self.data[index] = color.0;
|
||||
self.data[index + 1] = color.1;
|
||||
self.data[index + 2] = color.2;
|
||||
}
|
||||
}
|
||||
|
||||
fn ypr_to_mat(ypr: Vec3) -> Mat3 {
|
||||
let Vec3 { x: yaw, y: pitch, z: roll } = ypr;
|
||||
let m_roll = mat3(
|
||||
roll.cos(), roll.sin(), 0.0,
|
||||
-roll.sin(), roll.cos(), 0.0,
|
||||
0.0, 0.0, 1.0);
|
||||
let m_yaw = mat3(
|
||||
yaw.cos(), 0.0, yaw.sin(),
|
||||
0.0, 1.0, 0.0,
|
||||
-yaw.sin(), 0.0, yaw.cos());
|
||||
let m_pitch = mat3(
|
||||
1.0, 0.0, 0.0,
|
||||
0.0, pitch.cos(), -pitch.sin(),
|
||||
0.0, pitch.sin(), pitch.cos());
|
||||
m_roll * m_pitch * m_yaw
|
||||
}
|
||||
|
||||
fn render(mesh: &[Face], camera: impl Fn(Vec2) -> (Vec3, Vec3)) -> Image {
|
||||
let mut img = Image {
|
||||
w: W,
|
||||
h: H,
|
||||
data: vec![0; (3 * W * H) as usize],
|
||||
};
|
||||
let img_size = vec2(W as f32, H as f32);
|
||||
for y in 0..H {
|
||||
for x in 0..W {
|
||||
let img_coords = vec2(x as f32, y as f32);
|
||||
let off = (img_coords - img_size * 0.5) / img_size.y;
|
||||
let (base, ray) = camera(off);
|
||||
let mut dist = f32::INFINITY;
|
||||
for f in mesh {
|
||||
let color = clamp(to_ivec3(f.normal * 120.0 + 128.0), ivec3(0, 0, 0), ivec3(255, 255, 255));
|
||||
let fs = (0..3).map(|k| edge_dist(f.vertices[k], f.vertices[(k + 1) % 3], base, ray));
|
||||
if fs.into_iter().all(|f| f >= 0.0) {
|
||||
let m = Mat3 { c0: f.vertices[1] - f.vertices[0], c1: f.vertices[2] - f.vertices[0], c2: -ray };
|
||||
if let Some(m) = m.inverse() {
|
||||
let rel = m * (base - f.vertices[0]);
|
||||
if rel.z > dist {
|
||||
continue;
|
||||
}
|
||||
dist = rel.z;
|
||||
} else {
|
||||
continue;
|
||||
}
|
||||
img.put_pixel(x, y, Color(color.x as u8, color.y as u8, color.z as u8));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
img
|
||||
}
|
||||
|
||||
#[show_image::main]
|
||||
fn main() -> Result<(), Box<dyn Error>> {
|
||||
let args: Vec<String> = env::args().collect();
|
||||
let mesh = {
|
||||
let f = File::open(&args[1])?;
|
||||
let mut f = BufReader::new(f);
|
||||
load_mesh(&mut f)?
|
||||
};
|
||||
let m_view = ypr_to_mat(vec3(135.0 * PI / 180.0, -30.0 * PI / 180.0, 0.0f32));
|
||||
let m_camera = transpose(&m_view);
|
||||
let img = render(mesh.as_slice(), |off| {
|
||||
// perspective projection
|
||||
let base = vec3(0.0, 0.0, -40.0);
|
||||
let ray = vec3(off.x, off.y, 2.0);
|
||||
|
||||
// orthographic projection
|
||||
// let base = vec3(off.x, off.y, -10.0);
|
||||
// let ray = vec3(0.0, 0.0, 1.0);
|
||||
|
||||
(m_camera * base, m_camera * ray)
|
||||
});
|
||||
|
||||
let window = show_image::create_window("Raytracing", WindowOptions::default())?;
|
||||
let image = ImageView::new(ImageInfo::rgb8(W as u32, H as u32), img.data());
|
||||
window.set_image("image", image)?;
|
||||
window.wait_until_destroyed()?;
|
||||
Ok(())
|
||||
}
|
||||
|
||||
fn edge_dist(a: Vec3, b: Vec3, base: Vec3, dir: Vec3) -> f32 {
|
||||
// Note: given that the input is not arbitrary but comes from a cartesian product of certain (a, b) pairs and certain (base, dir) pairs, this can be optimized from Cnm to an+bm+cnm with c<C.
|
||||
Mat3 { c0: b - a, c1: base - a, c2: -dir }.determinant()
|
||||
}
|
||||
|
|
@ -1,5 +1,5 @@
|
|||
use std::io;
|
||||
use glam::{vec2, vec3, Vec2, Vec3};
|
||||
use glm::{vec2, vec3, Vec2, Vec3};
|
||||
|
||||
#[derive(Copy, Clone, Debug)]
|
||||
struct ObjVertex {
|
||||
Loading…
Reference in New Issue
Block a user