Compare commits
No commits in common. "97286085ab6975dac32e91fcafa16f498ae6e698" and "cbce6ccc444f0ff8d402247c61b7f3f1894d0ea7" have entirely different histories.
97286085ab
...
cbce6ccc44
|
|
@ -11,9 +11,6 @@ panic = 'abort'
|
|||
[profile.dev.package."*"]
|
||||
opt-level = 3
|
||||
|
||||
[profile.test.package."*"]
|
||||
opt-level = 3
|
||||
|
||||
[dependencies]
|
||||
rand = "0.8.5"
|
||||
glam = { version = "0.27.0", features = ["approx", "fast-math", "rand"] }
|
||||
|
|
|
|||
|
|
@ -1 +0,0 @@
|
|||
hard_tabs = true
|
||||
|
|
@ -12,12 +12,8 @@ pub trait FloatExt2<T>: bounds::Pair<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)
|
||||
}
|
||||
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)]
|
||||
|
|
|
|||
|
|
@ -52,73 +52,29 @@ pub struct QuadraticAccelerator {
|
|||
|
||||
/// Продолжает функцию 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)
|
||||
}
|
||||
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
|
||||
}
|
||||
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()
|
||||
}
|
||||
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,
|
||||
)
|
||||
}
|
||||
pub fn x(&self, u: f32) -> f32 { extend_linear(u, |u| (self.a() * u.abs() + self.b()) * u, self.internal, self.external) }
|
||||
pub fn u(&self, x: f32) -> f32 { extend_linear(x, |x| 0.5 * x.signum() * (-self.b() + self.root(x)) / self.a(), self.external, self.internal) }
|
||||
pub fn dx(&self, u: f32) -> f32 { extend_const(u, |u| 2.0 * self.a() * u.abs() + self.b(), self.internal, 1.0) }
|
||||
pub fn du(&self, x: f32) -> f32 { extend_const(x, |x| 1.0 / self.root(x), self.external, 1.0) }
|
||||
pub fn d2u(&self, x: f32) -> f32 { extend_const(x, |x| -2.0 * x.signum() * self.a() * self.root(x).powi(-3), self.external, 0.0) }
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod test {
|
||||
use approx::{abs_diff_eq, assert_abs_diff_eq, AbsDiffEq};
|
||||
use approx::{abs_diff_eq, AbsDiffEq, assert_abs_diff_eq};
|
||||
|
||||
use super::*;
|
||||
|
||||
|
|
@ -137,45 +93,23 @@ mod test {
|
|||
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"
|
||||
);
|
||||
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_limiter(SmoothstepLimiter { min: 20.0, max: 30.0 }, 20.0, 30.0, 1.0 / 32.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_smootherstep_limiter() {
|
||||
test_limiter(
|
||||
SmootherstepLimiter {
|
||||
min: 20.0,
|
||||
max: 30.0,
|
||||
},
|
||||
20.0,
|
||||
30.0,
|
||||
1.0 / 32.0,
|
||||
);
|
||||
test_limiter(SmootherstepLimiter { min: 20.0, max: 30.0 }, 20.0, 30.0, 1.0 / 32.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_quadratic_accelerator() {
|
||||
let testee = super::QuadraticAccelerator {
|
||||
internal: 100.0,
|
||||
external: 150.0,
|
||||
};
|
||||
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;
|
||||
|
|
@ -187,51 +121,30 @@ mod test {
|
|||
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"
|
||||
);
|
||||
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"
|
||||
);
|
||||
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"
|
||||
);
|
||||
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"
|
||||
);
|
||||
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"
|
||||
);
|
||||
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"
|
||||
);
|
||||
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"
|
||||
);
|
||||
assert!(abs_diff_eq!(dudx, 1.0, epsilon = ε), "At u={u}:\ndu/dx * dx/du: {dudx}\n");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -4,16 +4,15 @@ use flo_canvas::*;
|
|||
use flo_draw::*;
|
||||
use glam::*;
|
||||
|
||||
use crate::types::FlatTraceResult;
|
||||
use riemann::{trace_iter, Metric};
|
||||
use riemann::{Metric, trace_iter};
|
||||
use tube::metric::Tube;
|
||||
use tube::Space;
|
||||
use tube::Subspace::{Boundary, Inner, Outer};
|
||||
use types::{Location, Object, Ray};
|
||||
|
||||
mod float_fun;
|
||||
mod fns;
|
||||
mod riemann;
|
||||
mod fns;
|
||||
mod float_fun;
|
||||
mod tube;
|
||||
mod types;
|
||||
|
||||
|
|
@ -21,9 +20,7 @@ const DT: f32 = 0.1;
|
|||
|
||||
fn draw_loop(gc: &mut Vec<Draw>, mut pts: impl Iterator<Item=Vec2>) {
|
||||
gc.new_path();
|
||||
let Some(first) = pts.next() else {
|
||||
return;
|
||||
};
|
||||
let Some(first) = pts.next() else { return; };
|
||||
gc.move_to(first.x, first.y);
|
||||
for pt in pts {
|
||||
gc.line_to(pt.x, pt.y);
|
||||
|
|
@ -71,21 +68,10 @@ pub fn main() {
|
|||
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_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),
|
||||
);
|
||||
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 {
|
||||
|
|
@ -96,113 +82,83 @@ pub fn main() {
|
|||
gc.fill();
|
||||
|
||||
gc.stroke_color(Color::Rgba(0.0, 0.0, 0.0, 0.5));
|
||||
draw_loop(
|
||||
gc,
|
||||
itertools_num::linspace(0.0, 2.0 * PI, circle_segments)
|
||||
.skip(1)
|
||||
.map(|φ| {
|
||||
draw_loop(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(|φ| {
|
||||
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(|φ| {
|
||||
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
|
||||
}),
|
||||
);
|
||||
space.trace_iter(Ray { pos, dir }).nth(n as usize).unwrap().pos
|
||||
}));
|
||||
}
|
||||
});
|
||||
});
|
||||
}
|
||||
|
||||
fn rel_to_abs(space: &impl Metric, base: &Location, rel: Vec2, steps: usize) -> Vec2 {
|
||||
let c = 1.0 / (steps as f32);
|
||||
trace_iter(space, base.pos, base.rot * rel, c * rel.length())
|
||||
.nth(steps - 1)
|
||||
.unwrap()
|
||||
}
|
||||
|
||||
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);
|
||||
gc.move_to(pos.x - r, pos.y + r);
|
||||
gc.line_to(pos.x + r, pos.y - r);
|
||||
}
|
||||
|
||||
fn draw_ray_2(gc: &mut Vec<Draw>, space: &Space, base: Vec2, dir: Vec2) {
|
||||
fn trace_to_flat(gc: &mut Vec<Draw>, space: &Space, ray: Ray) -> (Ray, FlatTraceResult) {
|
||||
for ray in space.trace_iter(ray).skip(1) {
|
||||
gc.line_to(ray.pos.x, ray.pos.y);
|
||||
match space.which_subspace(ray.pos) {
|
||||
Inner => return (ray, space.trace_inner(ray)),
|
||||
Outer => return (ray, space.trace_outer(ray)),
|
||||
Boundary => continue,
|
||||
};
|
||||
}
|
||||
unreachable!("Space::trace_iter terminated!")
|
||||
}
|
||||
|
||||
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..100 {
|
||||
let ret;
|
||||
(ray, ret) = trace_to_flat(gc, space, ray);
|
||||
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];
|
||||
let apx_hit_pos = rel_to_abs(&space.tube, &obj.loc, hit.rel.pos, 128);
|
||||
// assert_abs_diff_eq!(apx_hit_pos, hit.pos, epsilon=1.0);
|
||||
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));
|
||||
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 = rel_to_abs(&space.tube, &obj.loc, rel2, 128);
|
||||
draw_cross(&mut hits, pos2, 1.0);
|
||||
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;
|
||||
if let Some(r) = ret.end {
|
||||
ray = r
|
||||
} else {
|
||||
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);
|
||||
}
|
||||
|
|
@ -245,10 +201,7 @@ fn draw_track(gc: &mut Vec<Draw>, space: &Space, start: Vec2, dir: Vec2) {
|
|||
// let mut loc = Location { pos: start, rot: Mat2::IDENTITY };
|
||||
// let dir = space.tube.globalize(start, dir);
|
||||
// let v = space.tube.normalize(start, dir);
|
||||
let mut loc = Location {
|
||||
pos: start,
|
||||
rot: mat2(dir, vec2(-dir.y, dir.x)),
|
||||
};
|
||||
let mut loc = Location { 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;
|
||||
|
|
@ -290,18 +243,8 @@ trait Renderable {
|
|||
impl Renderable for Tube {
|
||||
fn render(&self, gc: &mut Vec<Draw>) {
|
||||
gc.new_path();
|
||||
gc.rect(
|
||||
-self.outer_radius,
|
||||
-self.external_halflength,
|
||||
self.outer_radius,
|
||||
self.external_halflength,
|
||||
);
|
||||
gc.rect(
|
||||
-self.inner_radius,
|
||||
-self.external_halflength,
|
||||
self.inner_radius,
|
||||
self.external_halflength,
|
||||
);
|
||||
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();
|
||||
|
|
|
|||
|
|
@ -62,6 +62,7 @@ pub struct TraceIter<'a, M: Metric> {
|
|||
space: &'a M,
|
||||
p: Vec2,
|
||||
v: Vec2,
|
||||
dt: f32,
|
||||
}
|
||||
|
||||
impl<'a, M: Metric> Iterator for TraceIter<'a, M> {
|
||||
|
|
@ -69,8 +70,8 @@ impl<'a, M: Metric> Iterator for TraceIter<'a, M> {
|
|||
|
||||
fn next(&mut self) -> Option<Self::Item> {
|
||||
let a: Vec2 = -contract2(krist(self.space, self.p), self.v);
|
||||
self.v += a;
|
||||
self.p += self.v;
|
||||
self.v = self.v + a * self.dt;
|
||||
self.p = self.p + self.v * self.dt;
|
||||
Some(self.p)
|
||||
}
|
||||
}
|
||||
|
|
@ -79,7 +80,8 @@ pub fn trace_iter<M: Metric>(space: &M, base: Vec2, dir: Vec2, dt: f32) -> Trace
|
|||
TraceIter {
|
||||
space,
|
||||
p: base,
|
||||
v: dt * space.normalize_vec_at(base, dir),
|
||||
v: space.normalize_vec_at(base, dir),
|
||||
dt,
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -88,11 +90,7 @@ pub fn krist(space: &impl Metric, pos: Vec2) -> Tens2 {
|
|||
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>()
|
||||
})
|
||||
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 {
|
||||
|
|
@ -171,109 +169,27 @@ pub mod samples {
|
|||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use approx::assert_abs_diff_eq;
|
||||
|
||||
use glam::{mat2, vec2, Mat2};
|
||||
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.])
|
||||
);
|
||||
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.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.)
|
||||
);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_trace_iter() {
|
||||
let metric = samples::ScaledMetric {
|
||||
scale: vec2(2., 4.),
|
||||
};
|
||||
assert_eq!(
|
||||
trace_iter(&metric, vec2(3., 5.), vec2(1., 0.), 1.)
|
||||
.nth(7)
|
||||
.unwrap(),
|
||||
vec2(7., 5.)
|
||||
);
|
||||
assert_eq!(
|
||||
trace_iter(&metric, vec2(3., 5.), vec2(2., 0.), 1.)
|
||||
.nth(7)
|
||||
.unwrap(),
|
||||
vec2(7., 5.)
|
||||
);
|
||||
assert_eq!(
|
||||
trace_iter(&metric, vec2(3., 5.), vec2(1., 0.), 0.5)
|
||||
.nth(7)
|
||||
.unwrap(),
|
||||
vec2(5., 5.)
|
||||
);
|
||||
assert_eq!(
|
||||
trace_iter(&metric, vec2(3., 5.), vec2(0., 1.), 1.)
|
||||
.nth(9)
|
||||
.unwrap(),
|
||||
vec2(3., 7.5)
|
||||
);
|
||||
assert_eq!(
|
||||
trace_iter(&metric, vec2(3., 5.), vec2(0., 4.), 1.)
|
||||
.nth(9)
|
||||
.unwrap(),
|
||||
vec2(3., 7.5)
|
||||
);
|
||||
assert_eq!(
|
||||
trace_iter(&metric, vec2(3., 5.), vec2(0., 1.), 0.5)
|
||||
.nth(9)
|
||||
.unwrap(),
|
||||
vec2(3., 6.25)
|
||||
);
|
||||
assert_abs_diff_eq!(
|
||||
trace_iter(
|
||||
&metric,
|
||||
vec2(3., 5.),
|
||||
vec2(0.5, 0.25),
|
||||
std::f32::consts::SQRT_2
|
||||
)
|
||||
.nth(7)
|
||||
.unwrap(),
|
||||
vec2(7., 7.),
|
||||
epsilon = 1e-5
|
||||
);
|
||||
assert_eq!(metric.globalize(rng.gen(), vec2(1., 1.)), vec2(1. / 3., 1. / 4.));
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
use glam::{vec2, Mat2, Vec2};
|
||||
use glam::{Mat2, Vec2, vec2};
|
||||
|
||||
use crate::riemann::Metric;
|
||||
use crate::types::{Location, Ray};
|
||||
|
|
@ -10,22 +10,15 @@ pub trait FlatCoordinateSystem<T> {
|
|||
fn global_to_flat(&self, v: T) -> T;
|
||||
}
|
||||
|
||||
pub trait FlatRegion:
|
||||
FlatCoordinateSystem<Vec2> + FlatCoordinateSystem<Ray> + FlatCoordinateSystem<Location>
|
||||
{
|
||||
pub trait FlatRegion: FlatCoordinateSystem<Vec2> + FlatCoordinateSystem<Ray> + FlatCoordinateSystem<Location> {
|
||||
// Измеряет расстояние до выхода за пределы области вдоль луча ray. Луч задаётся в плоской СК.
|
||||
fn distance_to_boundary(&self, _ray: Ray) -> Option<f32> {
|
||||
None
|
||||
}
|
||||
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()
|
||||
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()
|
||||
|
|
@ -66,11 +59,7 @@ impl<T: FlatCoordinateSystem<Vec2> + MetricCS> FlatCoordinateSystem<Location> fo
|
|||
|
||||
pub struct InnerCS(pub Tube);
|
||||
|
||||
impl MetricCS for InnerCS {
|
||||
fn global_metric(&self) -> &impl Metric {
|
||||
&self.0
|
||||
}
|
||||
}
|
||||
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 {
|
||||
|
|
@ -85,31 +74,20 @@ impl FlatCoordinateSystem<Vec2> for InnerCS {
|
|||
|
||||
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)
|
||||
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 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),
|
||||
};
|
||||
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));
|
||||
let y = self.0.y(v - v.signum() * (self.0.external_halflength - self.0.internal_halflength));
|
||||
vec2(x, y)
|
||||
} else {
|
||||
pos
|
||||
|
|
@ -117,13 +95,10 @@ impl FlatCoordinateSystem<Vec2> for OuterCS {
|
|||
}
|
||||
|
||||
fn global_to_flat(&self, pos: Vec2) -> Vec2 {
|
||||
let inner = Rect {
|
||||
size: vec2(self.0.inner_radius + 1.0, self.0.external_halflength),
|
||||
};
|
||||
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);
|
||||
let v = self.0.v(y) + y.signum() * (self.0.external_halflength - self.0.internal_halflength);
|
||||
vec2(u, v) // в плоском продолжении СК Outer на область Inner
|
||||
} else {
|
||||
pos
|
||||
|
|
@ -133,17 +108,14 @@ impl FlatCoordinateSystem<Vec2> for OuterCS {
|
|||
|
||||
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)
|
||||
Rect { size: vec2(self.0.outer_radius, self.0.external_halflength) }.trace_into(ray)
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use approx::{assert_abs_diff_eq, AbsDiffEq};
|
||||
use glam::{mat2, vec2, Mat2, Vec2};
|
||||
use approx::{AbsDiffEq, assert_abs_diff_eq};
|
||||
use glam::{Mat2, mat2, vec2, Vec2};
|
||||
use itertools_num::linspace;
|
||||
|
||||
use crate::riemann::samples;
|
||||
|
|
@ -154,179 +126,57 @@ mod tests {
|
|||
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
|
||||
}
|
||||
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
|
||||
fn global_metric(&self) -> &impl Metric { &self.0 }
|
||||
}
|
||||
}
|
||||
let cs = Scaled(
|
||||
samples::ScaledMetric {
|
||||
scale: vec2(3., 4.),
|
||||
},
|
||||
vec2(2., 3.),
|
||||
);
|
||||
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.))
|
||||
}
|
||||
);
|
||||
assert_eq!(cs.global_to_flat(Ray { pos: vec2(7., 3.), dir: vec2(3., 2.) }), Ray { pos: vec2(23., 15.), dir: vec2(9., 8.) });
|
||||
assert_eq!(cs.flat_to_global(Ray { pos: vec2(23., 15.), dir: vec2(9., 8.) }), Ray { pos: vec2(7., 3.), dir: vec2(3., 2.) });
|
||||
assert_eq!(cs.global_to_flat(Location { pos: vec2(2., 1.), rot: mat2(vec2(0., 1.), vec2(-1., 0.)) }), Location { pos: vec2(8., 7.), rot: mat2(vec2(0., 4.), vec2(-3., 0.)) });
|
||||
assert_eq!(cs.flat_to_global(Location { pos: vec2(2., 1.), rot: mat2(vec2(0., 1.), vec2(-1., 0.)) }), Location { pos: vec2(0., -0.5), rot: mat2(vec2(0., 0.25), vec2(-1. / 3., 0.)) });
|
||||
}
|
||||
|
||||
fn test_flat_region(
|
||||
region: &impl FlatRegion,
|
||||
range_global: (Vec2, Vec2),
|
||||
range_flat: (Vec2, Vec2),
|
||||
) {
|
||||
#[allow(non_upper_case_globals)]
|
||||
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)
|
||||
);
|
||||
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),
|
||||
) {
|
||||
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 = ε);
|
||||
}
|
||||
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
|
||||
);
|
||||
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
|
||||
);
|
||||
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
|
||||
);
|
||||
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
|
||||
);
|
||||
assert_eq_at!(pos_flat, region.global_to_flat(region.flat_to_global(Location { pos: pos_global, rot: Mat2::IDENTITY })).rot, Mat2::IDENTITY);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -339,21 +189,9 @@ mod tests {
|
|||
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_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]
|
||||
|
|
@ -365,68 +203,23 @@ mod tests {
|
|||
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)),
|
||||
);
|
||||
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
|
||||
);
|
||||
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;
|
||||
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));
|
||||
|
|
@ -435,42 +228,18 @@ mod tests {
|
|||
// 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
|
||||
);
|
||||
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
|
||||
);
|
||||
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;
|
||||
let v = mapper.global_to_flat(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos.y;
|
||||
assert!(v > 200.0);
|
||||
assert!(v > y);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
use glam::{f32, vec2, Mat2, Vec2};
|
||||
use glam::{f32, Mat2, Vec2, vec2};
|
||||
|
||||
use crate::fns::{self, Limiter};
|
||||
use crate::riemann::{Decomp2, Metric, Tens2};
|
||||
|
|
@ -12,31 +12,13 @@ pub struct Tube {
|
|||
}
|
||||
|
||||
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,
|
||||
}
|
||||
}
|
||||
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)
|
||||
}
|
||||
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 {
|
||||
|
|
@ -71,7 +53,7 @@ impl Metric for Tube {
|
|||
#[cfg(test)]
|
||||
mod test {
|
||||
use approx::assert_abs_diff_eq;
|
||||
use glam::{vec2, Vec2};
|
||||
use glam::{Vec2, vec2};
|
||||
use itertools_num::linspace;
|
||||
|
||||
use crate::riemann::{Decomp2, Metric};
|
||||
|
|
@ -84,9 +66,7 @@ mod 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)
|
||||
}
|
||||
fn sqrt_at(&self, pos: Vec2) -> Decomp2 { self.0.sqrt_at(pos) }
|
||||
}
|
||||
let testee = Tube {
|
||||
inner_radius: 30.0,
|
||||
|
|
@ -98,13 +78,8 @@ mod test {
|
|||
let epsilon = 1.0e-3;
|
||||
let margin = 1.0 / 16.0;
|
||||
let mul = 1.0 + margin;
|
||||
for x in itertools_num::linspace(-mul * testee.outer_radius, mul * testee.outer_radius, 100)
|
||||
{
|
||||
for y in itertools_num::linspace(
|
||||
-mul * testee.external_halflength,
|
||||
mul * testee.external_halflength,
|
||||
100,
|
||||
) {
|
||||
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);
|
||||
|
|
@ -158,16 +133,8 @@ mod test {
|
|||
let ε = 1e-3;
|
||||
let off = 10.0;
|
||||
let steps = 4096;
|
||||
for ax in linspace(
|
||||
-space.tube.inner_radius + ε,
|
||||
space.tube.inner_radius - ε,
|
||||
20,
|
||||
) {
|
||||
for bx in linspace(
|
||||
-space.tube.inner_radius + ε,
|
||||
space.tube.inner_radius - ε,
|
||||
20,
|
||||
) {
|
||||
for ax in linspace(-space.tube.inner_radius + ε, space.tube.inner_radius - ε, 20) {
|
||||
for bx in linspace(-space.tube.inner_radius + ε, space.tube.inner_radius - ε, 20) {
|
||||
let a = 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));
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
use glam::{bool, f32, vec2, Mat2, Vec2};
|
||||
use glam::{bool, f32, Mat2, Vec2, vec2};
|
||||
|
||||
use coords::{FlatCoordinateSystem, InnerCS, OuterCS};
|
||||
use metric::Tube;
|
||||
|
|
@ -8,8 +8,8 @@ use crate::riemann;
|
|||
use crate::tube::coords::FlatRegion;
|
||||
use crate::types::{FlatTraceResult, Hit, Location, Object, Ray};
|
||||
|
||||
mod coords;
|
||||
pub mod metric;
|
||||
mod coords;
|
||||
|
||||
pub struct Space {
|
||||
pub tube: Tube,
|
||||
|
|
@ -48,13 +48,9 @@ impl Space {
|
|||
/// Выполняет один шаг перемещения. Работает в любой части пространства.
|
||||
/// 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 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,
|
||||
}
|
||||
Location { pos: p, rot: corr * loc.rot }
|
||||
}
|
||||
|
||||
pub fn trace_iter(&self, ray: Ray) -> impl Iterator<Item=Ray> + '_ {
|
||||
|
|
@ -89,28 +85,15 @@ impl Space {
|
|||
}
|
||||
|
||||
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()
|
||||
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> {
|
||||
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));
|
||||
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))
|
||||
|
|
@ -121,17 +104,8 @@ impl Space {
|
|||
.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,
|
||||
}
|
||||
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()
|
||||
}
|
||||
|
|
@ -144,9 +118,7 @@ impl Space {
|
|||
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()
|
||||
(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!"),
|
||||
}
|
||||
|
|
@ -160,10 +132,7 @@ struct Rect {
|
|||
impl Rect {
|
||||
/// Отражает луч, чтобы все координаты направления были положительны (допустимо благодаря симметрии Rect).
|
||||
fn flip_ray(ray: Ray) -> Ray {
|
||||
Ray {
|
||||
pos: ray.pos * ray.dir.signum(),
|
||||
dir: ray.dir.abs(),
|
||||
}
|
||||
Ray { pos: ray.pos * ray.dir.signum(), dir: ray.dir.abs() }
|
||||
}
|
||||
|
||||
fn is_inside(&self, pt: Vec2) -> bool {
|
||||
|
|
@ -176,12 +145,8 @@ impl Rect {
|
|||
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;
|
||||
}
|
||||
if t < 0.0 { return None; }
|
||||
if pt.cmpgt(self.size).any() { return None; }
|
||||
Some(t)
|
||||
}
|
||||
|
||||
|
|
@ -196,149 +161,27 @@ impl Rect {
|
|||
|
||||
#[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)
|
||||
}
|
||||
);
|
||||
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),
|
||||
};
|
||||
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(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_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)
|
||||
);
|
||||
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));
|
||||
}
|
||||
|
|
|
|||
|
|
@ -8,10 +8,7 @@ pub struct Ray {
|
|||
|
||||
impl Ray {
|
||||
pub fn forward(&self, dist: f32) -> Ray {
|
||||
Ray {
|
||||
pos: self.pos + self.dir * dist,
|
||||
dir: self.dir,
|
||||
}
|
||||
Ray { pos: self.pos + self.dir * dist, dir: self.dir }
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -19,10 +16,7 @@ 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,
|
||||
}
|
||||
Ray { pos: self * rhs.pos, dir: self * rhs.dir }
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -1,12 +1,13 @@
|
|||
use glam::*;
|
||||
use refraction::mesh_loader::load_mesh;
|
||||
use refraction::mesh_tracer::{trace_to_mesh, Mesh};
|
||||
use show_image::{ImageInfo, ImageView, WindowOptions};
|
||||
use std::env;
|
||||
mod mesh_loader;
|
||||
|
||||
use std::fs::File;
|
||||
use std::{env};
|
||||
use std::error::Error;
|
||||
use std::f32::consts::PI;
|
||||
use std::fs::File;
|
||||
use std::io::BufReader;
|
||||
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;
|
||||
|
|
@ -37,29 +38,56 @@ impl Image {
|
|||
}
|
||||
|
||||
fn ypr_to_mat(ypr: Vec3) -> Mat3 {
|
||||
let Vec3 {
|
||||
x: yaw,
|
||||
y: pitch,
|
||||
z: roll,
|
||||
} = ypr;
|
||||
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),
|
||||
);
|
||||
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()),
|
||||
);
|
||||
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()),
|
||||
);
|
||||
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 {
|
||||
|
|
@ -79,9 +107,7 @@ fn render(mesh: &Mesh, camera: impl Fn(Vec2) -> (Vec3, Vec3)) -> Image {
|
|||
} else {
|
||||
bkg
|
||||
};
|
||||
let color = (color * 255.0)
|
||||
.as_ivec3()
|
||||
.clamp(IVec3::splat(0), IVec3::splat(255));
|
||||
let color = (color * 255.0).as_ivec3().clamp(IVec3::splat(0), IVec3::splat(255));
|
||||
img.put_pixel(x, y, Color(color.x as u8, color.y as u8, color.z as u8));
|
||||
}
|
||||
}
|
||||
|
|
@ -99,11 +125,7 @@ fn main() -> Result<(), Box<dyn Error>> {
|
|||
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_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
|
||||
|
|
@ -122,3 +144,8 @@ fn main() -> Result<(), Box<dyn Error>> {
|
|||
}
|
||||
}
|
||||
}
|
||||
|
||||
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()
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
use glam::{vec2, vec3, Vec2, Vec3};
|
||||
use std::io;
|
||||
use glam::{vec2, vec3, Vec2, Vec3};
|
||||
|
||||
#[derive(Copy, Clone, Debug)]
|
||||
struct ObjVertex {
|
||||
|
|
@ -35,23 +35,14 @@ impl ObjMesh {
|
|||
}
|
||||
|
||||
fn parse_fv(desc: &&str) -> ObjVertex {
|
||||
let tokens: Vec<_> = desc
|
||||
.split('/')
|
||||
.map(|s| s.parse::<usize>().unwrap() - 1)
|
||||
.collect();
|
||||
let tokens: Vec<_> = desc.split('/').map(|s| s.parse::<usize>().unwrap() - 1).collect();
|
||||
assert_eq!(tokens.len(), 3);
|
||||
ObjVertex {
|
||||
vertex: tokens[0],
|
||||
tex_coord: tokens[1],
|
||||
normal: tokens[2],
|
||||
}
|
||||
ObjVertex { vertex: tokens[0], tex_coord: tokens[1], normal: tokens[2] }
|
||||
}
|
||||
|
||||
fn parse_f(tokens: &[&str]) -> ObjFace {
|
||||
let vertices: Vec<_> = tokens.iter().map(ObjMesh::parse_fv).collect();
|
||||
ObjFace {
|
||||
vertices: vertices.as_slice().try_into().unwrap(),
|
||||
}
|
||||
ObjFace { vertices: vertices.as_slice().try_into().unwrap() }
|
||||
}
|
||||
|
||||
fn read(f: &mut impl io::BufRead) -> io::Result<ObjMesh> {
|
||||
|
|
@ -66,7 +57,13 @@ impl ObjMesh {
|
|||
if f.read_line(&mut line)? == 0 {
|
||||
break;
|
||||
}
|
||||
let tokens: Vec<&str> = line.trim().split('#').next().unwrap().split(' ').collect();
|
||||
let tokens: Vec<&str> = line
|
||||
.trim()
|
||||
.split('#')
|
||||
.next()
|
||||
.unwrap()
|
||||
.split(' ')
|
||||
.collect();
|
||||
match tokens[0] {
|
||||
"v" => result.vertices.push(Self::parse_v3(&tokens[1..])),
|
||||
"vn" => result.normals.push(Self::parse_v3(&tokens[1..])),
|
||||
|
|
@ -79,13 +76,12 @@ impl ObjMesh {
|
|||
}
|
||||
|
||||
fn flatten(&self) -> Vec<Face> {
|
||||
self.faces
|
||||
.iter()
|
||||
.map(|face| Face {
|
||||
self.faces.iter().map(|face| {
|
||||
Face {
|
||||
vertices: face.vertices.map(|iv| self.vertices[iv.vertex]),
|
||||
normal: self.normals[face.vertices[0].normal],
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
}).collect()
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -1,2 +0,0 @@
|
|||
pub mod mesh_loader;
|
||||
pub mod mesh_tracer;
|
||||
|
|
@ -1,44 +0,0 @@
|
|||
use crate::mesh_loader::Face;
|
||||
use glam::{mat3, Vec3};
|
||||
|
||||
pub type Mesh = [Face];
|
||||
|
||||
pub struct TraceResult {
|
||||
pub distance: f32,
|
||||
pub normal: Vec3,
|
||||
}
|
||||
|
||||
pub fn trace_to_mesh_all(
|
||||
mesh: &Mesh,
|
||||
base: Vec3,
|
||||
ray: Vec3,
|
||||
) -> impl Iterator<Item = TraceResult> + '_ {
|
||||
mesh.iter().filter_map(move |f| {
|
||||
let fs = (0..3).map(|k| edge_dist(f.vertices[k], f.vertices[(k + 1) % 3], base, ray));
|
||||
if fs.into_iter().any(|f| f < 0.0) {
|
||||
return None;
|
||||
}
|
||||
let m = mat3(
|
||||
f.vertices[1] - f.vertices[0],
|
||||
f.vertices[2] - f.vertices[0],
|
||||
-ray,
|
||||
);
|
||||
let m = m.inverse();
|
||||
let rel = m * (base - f.vertices[0]);
|
||||
Some(TraceResult {
|
||||
distance: rel.z,
|
||||
normal: f.normal,
|
||||
})
|
||||
})
|
||||
}
|
||||
|
||||
pub fn trace_to_mesh(mesh: &Mesh, base: Vec3, ray: Vec3) -> Option<TraceResult> {
|
||||
trace_to_mesh_all(mesh, base, ray)
|
||||
.filter(|tr| tr.distance >= 0.0)
|
||||
.min_by(|a, b| a.distance.total_cmp(&b.distance))
|
||||
}
|
||||
|
||||
fn edge_dist(a: Vec3, b: Vec3, base: Vec3, dir: Vec3) -> f32 {
|
||||
// Note: given that the input is not arbitrary but comes from a cartesian product of certain (a, b) pairs and certain (base, dir) pairs, this can be optimized from Cnm to an+bm+cnm with c<C.
|
||||
mat3(b - a, base - a, -dir).determinant()
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user