Compare commits
No commits in common. "master" and "shaped" have entirely different histories.
|
|
@ -1,18 +0,0 @@
|
|||
{
|
||||
"Auto_generated": "This file is auto-generated. Any extra tags or formatting will be lost",
|
||||
"target_sets": [
|
||||
{
|
||||
"cmake_config": "",
|
||||
"directory": "%B",
|
||||
"loaded_via_cmake": false,
|
||||
"name": "Cargo",
|
||||
"targets": [
|
||||
{
|
||||
"build_cmd": "cargo build --bin wireframe",
|
||||
"name": "wireframe",
|
||||
"run_cmd": "cargo run --bin wireframe"
|
||||
}
|
||||
]
|
||||
}
|
||||
]
|
||||
}
|
||||
2565
Cargo.lock
generated
2565
Cargo.lock
generated
File diff suppressed because it is too large
Load Diff
19
Cargo.toml
19
Cargo.toml
|
|
@ -3,32 +3,19 @@ name = "refraction"
|
|||
version = "0.1.0"
|
||||
edition = "2021"
|
||||
|
||||
[lints.rust]
|
||||
mixed_script_confusables = "allow"
|
||||
# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
|
||||
|
||||
[profile.dev]
|
||||
panic = 'abort'
|
||||
opt-level = 2
|
||||
|
||||
[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"] }
|
||||
show-image = "0.14.0"
|
||||
glam = "0.27.0"
|
||||
flo_draw = "0.3.1"
|
||||
flo_canvas = "0.3.1"
|
||||
itertools-num = "0.1.3"
|
||||
winit = "0.29"
|
||||
itertools = "0.13.0"
|
||||
wgpu = "22.1.0"
|
||||
bytemuck = { version = "1.18.0", features = ["derive"] }
|
||||
pollster = "0.3.0"
|
||||
|
||||
[dev-dependencies]
|
||||
approx = "0.5.1"
|
||||
rand = "0.8.5"
|
||||
rand_pcg = "0.3.1"
|
||||
|
|
|
|||
|
|
@ -1,140 +0,0 @@
|
|||
# Blender 3.6.5
|
||||
# www.blender.org
|
||||
mtllib spacecraft2.1.mtl
|
||||
o Cube
|
||||
v -1.000000 0.000000 5.000000
|
||||
v -3.000000 0.000000 4.000000
|
||||
v -4.000000 0.000000 7.000000
|
||||
v -8.000000 0.000000 7.000000
|
||||
v -9.000000 0.000000 4.000000
|
||||
v -3.000000 0.000000 0.000000
|
||||
v -6.000000 0.000000 -7.000000
|
||||
v -2.000000 0.000000 -3.000000
|
||||
v -7.000000 1.000000 4.000000
|
||||
v 0.000000 2.000000 0.000000
|
||||
v 0.000000 2.000000 2.000000
|
||||
v 0.000000 0.000000 6.000000
|
||||
v 0.000000 0.000000 -4.000000
|
||||
v -5.000000 0.000000 2.000000
|
||||
v 1.000000 0.000000 5.000000
|
||||
v 3.000000 0.000000 4.000000
|
||||
v 4.000000 0.000000 7.000000
|
||||
v 8.000000 0.000000 7.000000
|
||||
v 9.000000 0.000000 4.000000
|
||||
v 3.000000 0.000000 0.000000
|
||||
v 6.000000 0.000000 -7.000000
|
||||
v 2.000000 0.000000 -3.000000
|
||||
v 7.000000 1.000000 4.000000
|
||||
v 5.000000 0.000000 2.000000
|
||||
v -7.000000 -1.000000 4.000000
|
||||
v 0.000000 -2.000000 0.000000
|
||||
v 0.000000 -2.000000 2.000000
|
||||
v 7.000000 -1.000000 4.000000
|
||||
vn -0.5455 0.8182 -0.1818
|
||||
vn -0.4851 0.4851 -0.7276
|
||||
vn -0.4439 0.8878 -0.1211
|
||||
vn -0.0000 0.9487 0.3162
|
||||
vn -0.4423 0.8847 0.1474
|
||||
vn 0.2417 0.9670 0.0806
|
||||
vn -0.2408 0.8427 0.4815
|
||||
vn 0.1455 0.5819 0.8001
|
||||
vn -0.1414 0.9899 -0.0000
|
||||
vn -0.2182 0.8729 -0.4364
|
||||
vn -0.4082 0.8165 0.4082
|
||||
vn -0.4851 0.7276 -0.4851
|
||||
vn 0.4099 0.9110 -0.0455
|
||||
vn -0.0000 -1.0000 -0.0000
|
||||
vn 0.5455 0.8182 -0.1818
|
||||
vn 0.4851 0.4851 -0.7276
|
||||
vn 0.4439 0.8878 -0.1211
|
||||
vn 0.4423 0.8847 0.1474
|
||||
vn -0.2417 0.9670 0.0806
|
||||
vn 0.2408 0.8427 0.4815
|
||||
vn -0.1455 0.5819 0.8001
|
||||
vn 0.1414 0.9899 -0.0000
|
||||
vn 0.2182 0.8729 -0.4364
|
||||
vn 0.4082 0.8165 0.4082
|
||||
vn 0.4851 0.7276 -0.4851
|
||||
vn -0.4099 0.9110 -0.0455
|
||||
vn -0.5455 -0.8182 -0.1818
|
||||
vn -0.4851 -0.4851 -0.7276
|
||||
vn -0.4439 -0.8878 -0.1211
|
||||
vn -0.0000 -0.9487 0.3162
|
||||
vn -0.4423 -0.8847 0.1474
|
||||
vn 0.2417 -0.9670 0.0806
|
||||
vn -0.2408 -0.8427 0.4815
|
||||
vn 0.1455 -0.5819 0.8001
|
||||
vn -0.1414 -0.9899 -0.0000
|
||||
vn -0.2182 -0.8729 -0.4364
|
||||
vn -0.4082 -0.8165 0.4082
|
||||
vn -0.4851 -0.7276 -0.4851
|
||||
vn 0.4099 -0.9110 -0.0455
|
||||
vn 0.5455 -0.8182 -0.1818
|
||||
vn 0.4851 -0.4851 -0.7276
|
||||
vn 0.4439 -0.8878 -0.1211
|
||||
vn 0.4423 -0.8847 0.1474
|
||||
vn -0.2417 -0.9670 0.0806
|
||||
vn 0.2408 -0.8427 0.4815
|
||||
vn -0.1455 -0.5819 0.8001
|
||||
vn 0.1414 -0.9899 -0.0000
|
||||
vn 0.2182 -0.8729 -0.4364
|
||||
vn 0.4082 -0.8165 0.4082
|
||||
vn 0.4851 -0.7276 -0.4851
|
||||
vn -0.4099 -0.9110 -0.0455
|
||||
vt 0.000000 0.000000
|
||||
s 0
|
||||
usemtl Material
|
||||
f 6/1/1 10/1/1 8/1/1
|
||||
f 10/1/2 14/1/2 9/1/2
|
||||
f 5/1/3 9/1/3 7/1/3
|
||||
f 3/1/4 9/1/4 4/1/4
|
||||
f 4/1/5 9/1/5 5/1/5
|
||||
f 2/1/6 9/1/6 3/1/6
|
||||
f 1/1/7 11/1/7 2/1/7
|
||||
f 9/1/8 2/1/8 11/1/8
|
||||
f 9/1/9 11/1/9 10/1/9
|
||||
f 8/1/10 10/1/10 13/1/10
|
||||
f 1/1/11 12/1/11 11/1/11
|
||||
f 10/1/12 6/1/12 14/1/12
|
||||
f 7/1/13 9/1/13 14/1/13
|
||||
f 7/1/14 14/1/14 5/1/14
|
||||
f 20/1/15 22/1/15 10/1/15
|
||||
f 10/1/16 23/1/16 24/1/16
|
||||
f 19/1/17 21/1/17 23/1/17
|
||||
f 17/1/4 18/1/4 23/1/4
|
||||
f 18/1/18 19/1/18 23/1/18
|
||||
f 16/1/19 17/1/19 23/1/19
|
||||
f 15/1/20 16/1/20 11/1/20
|
||||
f 23/1/21 11/1/21 16/1/21
|
||||
f 23/1/22 10/1/22 11/1/22
|
||||
f 22/1/23 13/1/23 10/1/23
|
||||
f 15/1/24 11/1/24 12/1/24
|
||||
f 10/1/25 24/1/25 20/1/25
|
||||
f 21/1/26 24/1/26 23/1/26
|
||||
f 21/1/14 19/1/14 24/1/14
|
||||
f 6/1/27 8/1/27 26/1/27
|
||||
f 26/1/28 25/1/28 14/1/28
|
||||
f 5/1/29 7/1/29 25/1/29
|
||||
f 3/1/30 4/1/30 25/1/30
|
||||
f 4/1/31 5/1/31 25/1/31
|
||||
f 2/1/32 3/1/32 25/1/32
|
||||
f 1/1/33 2/1/33 27/1/33
|
||||
f 25/1/34 27/1/34 2/1/34
|
||||
f 25/1/35 26/1/35 27/1/35
|
||||
f 8/1/36 13/1/36 26/1/36
|
||||
f 1/1/37 27/1/37 12/1/37
|
||||
f 26/1/38 14/1/38 6/1/38
|
||||
f 7/1/39 14/1/39 25/1/39
|
||||
f 20/1/40 26/1/40 22/1/40
|
||||
f 26/1/41 24/1/41 28/1/41
|
||||
f 19/1/42 28/1/42 21/1/42
|
||||
f 17/1/30 28/1/30 18/1/30
|
||||
f 18/1/43 28/1/43 19/1/43
|
||||
f 16/1/44 28/1/44 17/1/44
|
||||
f 15/1/45 27/1/45 16/1/45
|
||||
f 28/1/46 16/1/46 27/1/46
|
||||
f 28/1/47 27/1/47 26/1/47
|
||||
f 22/1/48 26/1/48 13/1/48
|
||||
f 15/1/49 12/1/49 27/1/49
|
||||
f 26/1/50 20/1/50 24/1/50
|
||||
f 21/1/51 28/1/51 24/1/51
|
||||
|
|
@ -1,77 +0,0 @@
|
|||
# Blender 3.6.5
|
||||
# www.blender.org
|
||||
mtllib spacecraft2.2.mtl
|
||||
o Cube
|
||||
v -1.000000 0.000000 5.000000
|
||||
v -3.000000 0.000000 4.000000
|
||||
v -3.000000 0.000000 0.000000
|
||||
v -2.000000 0.000000 -3.000000
|
||||
v -7.000000 0.000000 4.000000
|
||||
v 0.000000 2.000000 0.000000
|
||||
v 0.000000 2.000000 2.000000
|
||||
v 0.000000 0.000000 6.000000
|
||||
v 0.000000 0.000000 -4.000000
|
||||
v 1.000000 0.000000 5.000000
|
||||
v 3.000000 0.000000 4.000000
|
||||
v 3.000000 0.000000 0.000000
|
||||
v 2.000000 0.000000 -3.000000
|
||||
v 7.000000 0.000000 4.000000
|
||||
v 0.000000 -2.000000 0.000000
|
||||
v 0.000000 -2.000000 2.000000
|
||||
vn -0.5455 0.8182 -0.1818
|
||||
vn -0.4851 0.7276 -0.4851
|
||||
vn -0.2408 0.8427 0.4815
|
||||
vn -0.0000 0.7071 0.7071
|
||||
vn -0.2747 0.9615 -0.0000
|
||||
vn -0.2182 0.8729 -0.4364
|
||||
vn -0.4082 0.8165 0.4082
|
||||
vn 0.5455 0.8182 -0.1818
|
||||
vn 0.4851 0.7276 -0.4851
|
||||
vn 0.2408 0.8427 0.4815
|
||||
vn 0.2747 0.9615 -0.0000
|
||||
vn 0.2182 0.8729 -0.4364
|
||||
vn 0.4082 0.8165 0.4082
|
||||
vn -0.5455 -0.8182 -0.1818
|
||||
vn -0.4851 -0.7276 -0.4851
|
||||
vn -0.2408 -0.8427 0.4815
|
||||
vn -0.0000 -0.7071 0.7071
|
||||
vn -0.2747 -0.9615 -0.0000
|
||||
vn -0.2182 -0.8729 -0.4364
|
||||
vn -0.4082 -0.8165 0.4082
|
||||
vn 0.5455 -0.8182 -0.1818
|
||||
vn 0.4851 -0.7276 -0.4851
|
||||
vn 0.2408 -0.8427 0.4815
|
||||
vn 0.2747 -0.9615 -0.0000
|
||||
vn 0.2182 -0.8729 -0.4364
|
||||
vn 0.4082 -0.8165 0.4082
|
||||
vt 0.000000 0.000000
|
||||
s 0
|
||||
usemtl Material
|
||||
f 3/1/1 6/1/1 4/1/1
|
||||
f 6/1/2 3/1/2 5/1/2
|
||||
f 1/1/3 7/1/3 2/1/3
|
||||
f 5/1/4 2/1/4 7/1/4
|
||||
f 5/1/5 7/1/5 6/1/5
|
||||
f 4/1/6 6/1/6 9/1/6
|
||||
f 1/1/7 8/1/7 7/1/7
|
||||
f 12/1/8 13/1/8 6/1/8
|
||||
f 6/1/9 14/1/9 12/1/9
|
||||
f 10/1/10 11/1/10 7/1/10
|
||||
f 14/1/4 7/1/4 11/1/4
|
||||
f 14/1/11 6/1/11 7/1/11
|
||||
f 13/1/12 9/1/12 6/1/12
|
||||
f 10/1/13 7/1/13 8/1/13
|
||||
f 3/1/14 4/1/14 15/1/14
|
||||
f 15/1/15 5/1/15 3/1/15
|
||||
f 1/1/16 2/1/16 16/1/16
|
||||
f 5/1/17 16/1/17 2/1/17
|
||||
f 5/1/18 15/1/18 16/1/18
|
||||
f 4/1/19 9/1/19 15/1/19
|
||||
f 1/1/20 16/1/20 8/1/20
|
||||
f 12/1/21 15/1/21 13/1/21
|
||||
f 15/1/22 12/1/22 14/1/22
|
||||
f 10/1/23 16/1/23 11/1/23
|
||||
f 14/1/17 11/1/17 16/1/17
|
||||
f 14/1/24 16/1/24 15/1/24
|
||||
f 13/1/25 15/1/25 9/1/25
|
||||
f 10/1/26 8/1/26 16/1/26
|
||||
|
|
@ -1,109 +0,0 @@
|
|||
# Blender 3.6.5
|
||||
# www.blender.org
|
||||
mtllib spacecraft2.mtl
|
||||
o Cube
|
||||
v -1.000000 0.000000 3.000000
|
||||
v -2.000000 0.000000 2.000000
|
||||
v -3.000000 0.000000 4.000000
|
||||
v -5.000000 0.000000 4.000000
|
||||
v -6.000000 0.000000 2.000000
|
||||
v -2.000000 0.000000 0.000000
|
||||
v -4.000000 0.000000 -6.000000
|
||||
v -1.000000 0.000000 -2.000000
|
||||
v 1.000000 0.000000 3.000000
|
||||
v 2.000000 0.000000 2.000000
|
||||
v 3.000000 0.000000 4.000000
|
||||
v 5.000000 0.000000 4.000000
|
||||
v 6.000000 0.000000 2.000000
|
||||
v 4.000000 0.000000 -6.000000
|
||||
v 2.000000 0.000000 0.000000
|
||||
v 1.000000 0.000000 -2.000000
|
||||
v -4.000000 1.000000 2.000000
|
||||
v 0.000000 1.000000 0.000000
|
||||
v 4.000000 1.000000 2.000000
|
||||
v 0.000000 1.000000 1.000000
|
||||
v -4.000000 -1.000000 2.000000
|
||||
v 0.000000 -1.000000 0.000000
|
||||
v 4.000000 -1.000000 2.000000
|
||||
v 0.000000 -1.000000 1.000000
|
||||
vn -0.0000 0.8944 -0.4472
|
||||
vn 0.4364 0.8729 -0.2182
|
||||
vn 0.3333 0.6667 -0.6667
|
||||
vn -0.3487 0.9300 -0.1162
|
||||
vn 0.4444 0.8889 -0.1111
|
||||
vn 0.4364 0.8729 0.2182
|
||||
vn -0.0000 0.8944 0.4472
|
||||
vn -0.4364 0.8729 0.2182
|
||||
vn -0.0000 1.0000 -0.0000
|
||||
vn 0.3015 0.9045 0.3015
|
||||
vn -0.2182 0.4364 0.8729
|
||||
vn -0.4364 0.8729 -0.2182
|
||||
vn -0.3333 0.6667 -0.6667
|
||||
vn 0.3487 0.9300 -0.1162
|
||||
vn -0.4444 0.8889 -0.1111
|
||||
vn -0.3015 0.9045 0.3015
|
||||
vn 0.2182 0.4364 0.8729
|
||||
vn -0.0000 -0.8944 -0.4472
|
||||
vn 0.4364 -0.8729 -0.2182
|
||||
vn 0.3333 -0.6667 -0.6667
|
||||
vn -0.3487 -0.9300 -0.1162
|
||||
vn 0.4444 -0.8889 -0.1111
|
||||
vn 0.4364 -0.8729 0.2182
|
||||
vn -0.0000 -0.8944 0.4472
|
||||
vn -0.4364 -0.8729 0.2182
|
||||
vn -0.0000 -1.0000 -0.0000
|
||||
vn 0.3015 -0.9045 0.3015
|
||||
vn -0.2182 -0.4364 0.8729
|
||||
vn -0.4364 -0.8729 -0.2182
|
||||
vn -0.3333 -0.6667 -0.6667
|
||||
vn 0.3487 -0.9300 -0.1162
|
||||
vn -0.4444 -0.8889 -0.1111
|
||||
vn -0.3015 -0.9045 0.3015
|
||||
vn 0.2182 -0.4364 0.8729
|
||||
vt 0.000000 0.000000
|
||||
s 0
|
||||
usemtl Material
|
||||
f 8/1/1 18/1/1 16/1/1
|
||||
f 16/1/2 18/1/2 15/1/2
|
||||
f 18/1/3 19/1/3 15/1/3
|
||||
f 15/1/4 19/1/4 14/1/4
|
||||
f 14/1/5 19/1/5 13/1/5
|
||||
f 13/1/6 19/1/6 12/1/6
|
||||
f 12/1/7 19/1/7 11/1/7
|
||||
f 10/1/8 11/1/8 19/1/8
|
||||
f 20/1/9 18/1/9 17/1/9
|
||||
f 1/1/7 9/1/7 20/1/7
|
||||
f 9/1/10 10/1/10 20/1/10
|
||||
f 20/1/11 10/1/11 19/1/11
|
||||
f 6/1/12 18/1/12 8/1/12
|
||||
f 18/1/13 6/1/13 17/1/13
|
||||
f 6/1/14 7/1/14 17/1/14
|
||||
f 5/1/15 17/1/15 7/1/15
|
||||
f 3/1/7 17/1/7 4/1/7
|
||||
f 4/1/8 17/1/8 5/1/8
|
||||
f 2/1/6 17/1/6 3/1/6
|
||||
f 1/1/16 20/1/16 2/1/16
|
||||
f 17/1/17 2/1/17 20/1/17
|
||||
f 20/1/9 19/1/9 18/1/9
|
||||
f 8/1/18 16/1/18 22/1/18
|
||||
f 16/1/19 15/1/19 22/1/19
|
||||
f 22/1/20 15/1/20 23/1/20
|
||||
f 15/1/21 14/1/21 23/1/21
|
||||
f 14/1/22 13/1/22 23/1/22
|
||||
f 13/1/23 12/1/23 23/1/23
|
||||
f 12/1/24 11/1/24 23/1/24
|
||||
f 10/1/25 23/1/25 11/1/25
|
||||
f 24/1/26 21/1/26 22/1/26
|
||||
f 1/1/24 24/1/24 9/1/24
|
||||
f 9/1/27 24/1/27 10/1/27
|
||||
f 24/1/28 23/1/28 10/1/28
|
||||
f 6/1/29 8/1/29 22/1/29
|
||||
f 22/1/30 21/1/30 6/1/30
|
||||
f 6/1/31 21/1/31 7/1/31
|
||||
f 5/1/32 7/1/32 21/1/32
|
||||
f 3/1/24 4/1/24 21/1/24
|
||||
f 4/1/25 5/1/25 21/1/25
|
||||
f 2/1/23 3/1/23 21/1/23
|
||||
f 1/1/33 2/1/33 24/1/33
|
||||
f 21/1/34 24/1/34 2/1/34
|
||||
f 24/1/26 22/1/26 23/1/26
|
||||
|
|
@ -1,2 +0,0 @@
|
|||
hard_tabs = true
|
||||
max_width = 120
|
||||
707
src/bin/flat.rs
Normal file
707
src/bin/flat.rs
Normal file
|
|
@ -0,0 +1,707 @@
|
|||
use std::collections::HashMap;
|
||||
use flo_draw::*;
|
||||
use flo_canvas::*;
|
||||
use glam::*;
|
||||
use riemann::{Decomp2, Metric, trace_iter};
|
||||
|
||||
#[cfg(test)]
|
||||
use approx::assert_abs_diff_eq;
|
||||
|
||||
const DT: f32 = 0.1;
|
||||
|
||||
pub fn main() {
|
||||
with_2d_graphics(move || {
|
||||
let canvas = create_drawing_window("Refraction");
|
||||
canvas.draw(|gc| {
|
||||
let space = Coil {
|
||||
scale: 3.0,
|
||||
r: 300.0,
|
||||
w: 50.0,
|
||||
m: 10.0,
|
||||
};
|
||||
let tube = Rect {
|
||||
inner_radius: 30.0,
|
||||
outer_radius: 50.0,
|
||||
internal_halflength: 100.0,
|
||||
external_halflength: 300.0,
|
||||
};
|
||||
let mut grid = Grid {
|
||||
hlines: vec![-tube.external_halflength, tube.external_halflength],
|
||||
vlines: vec![-tube.outer_radius, -tube.inner_radius, tube.inner_radius, tube.outer_radius],
|
||||
cells: HashMap::new(),
|
||||
};
|
||||
fn take(k: i32, data: &[f32]) -> f32 {
|
||||
if k < 0 {
|
||||
-f32::INFINITY
|
||||
} else if k as usize >= data.len() {
|
||||
f32::INFINITY
|
||||
} else {
|
||||
data[k as usize]
|
||||
}
|
||||
}
|
||||
for (i, j) in [(0, 0), (1, 0), (2, 0), (3, 0), (4, 0), (0, 1), (4, 1), (0, 2), (1, 2), (2, 2), (3, 2), (4, 2)] {
|
||||
grid.cells.insert((i as usize, j as usize),
|
||||
Box::new(FlatRect {
|
||||
min: vec2(take(i - 1, grid.vlines.as_slice()), take(j - 1, grid.hlines.as_slice())),
|
||||
max: vec2(take(i, grid.vlines.as_slice()), take(j, grid.hlines.as_slice())),
|
||||
}));
|
||||
}
|
||||
grid.cells.insert((2, 1), Box::new(RectInside { rect: tube }));
|
||||
println!("{:?}", grid.cells);
|
||||
gc.canvas_height(1000.0);
|
||||
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, &tube, &grid, vec2(-500.0, 0.0), vec2(1.0, 0.0), 1.0);
|
||||
gc.stroke_color(Color::Rgba(0.0, 0.5, 1.0, 0.5));
|
||||
draw_fan(gc, &tube, vec2(0.0, -0.5 * tube.internal_halflength), vec2(1.0, 1.0), 1.0);
|
||||
gc.stroke_color(Color::Rgba(0.0, 0.5, 1.0, 1.0));
|
||||
draw_fan_2(gc, &tube, &grid, vec2(0.0, -0.5 * tube.internal_halflength), vec2(1.0, 1.0), 1.0);
|
||||
|
||||
gc.new_path();
|
||||
for ((i, j), cell) in &grid.cells {
|
||||
let (a, b) = cell.local_bounds();
|
||||
gc.rect(a.x.clamp(-1000.0, 1000.0), a.y.clamp(-1000.0, 1000.0), b.x.clamp(-1000.0, 1000.0), b.y.clamp(-1000.0, 1000.0));
|
||||
}
|
||||
gc.fill_color(Color::Rgba(0.0, 1.0, 0.0, 0.3));
|
||||
gc.fill();
|
||||
|
||||
gc.new_dash_pattern();
|
||||
gc.dash_length(6.0);
|
||||
gc.new_dash_pattern();
|
||||
gc.new_path();
|
||||
gc.stroke_color(Color::Rgba(0.0, 0.0, 0.0, 0.5));
|
||||
gc.line_width(1.0);
|
||||
for hline in grid.hlines {
|
||||
gc.move_to(-1000.0, hline);
|
||||
gc.line_to(1000.0, hline);
|
||||
}
|
||||
for vline in grid.vlines {
|
||||
gc.move_to(vline, -1000.0);
|
||||
gc.line_to(vline, 1000.0);
|
||||
}
|
||||
gc.stroke();
|
||||
});
|
||||
});
|
||||
}
|
||||
|
||||
fn find_bucket(val: f32, splits: &[f32]) -> usize {
|
||||
for (k, &split) in splits.iter().enumerate() {
|
||||
if val < split {
|
||||
return k;
|
||||
}
|
||||
}
|
||||
splits.len()
|
||||
}
|
||||
|
||||
fn draw_ray_2(gc: &mut Vec<Draw>, space: &impl Metric, grid: &Grid, base: Vec2, dir: Vec2) {
|
||||
let dir = space.globalize(base, dir);
|
||||
gc.new_path();
|
||||
gc.move_to(base.x, base.y);
|
||||
let dt = DT;
|
||||
let mut p = base;
|
||||
let mut v = space.normalize(base, dir);
|
||||
for _ in 0..10000 {
|
||||
let a: Vec2 = -riemann::convolute(riemann::krist(space, p), v);
|
||||
v = v + a * dt;
|
||||
p = p + v * dt;
|
||||
gc.line_to(p.x, p.y);
|
||||
if p.abs().cmpgt(Vec2::splat(1000.0)).any() {
|
||||
break;
|
||||
}
|
||||
let i = find_bucket(p.x, grid.vlines.as_slice());
|
||||
let j = find_bucket(p.y, grid.hlines.as_slice());
|
||||
if let Some(cell) = grid.cells.get(&(i, j)) {
|
||||
gc.stroke();
|
||||
gc.new_dash_pattern();
|
||||
gc.dash_length(6.0);
|
||||
gc.new_path();
|
||||
gc.move_to(p.x, p.y);
|
||||
let ray = cell.ray_to_local(Ray { pos: p, dir: v });
|
||||
let Some(ray) = cell.to_boundary(ray) else {
|
||||
break;
|
||||
};
|
||||
Ray { pos: p, dir: v } = cell.ray_to_global(ray);
|
||||
gc.line_to(p.x, p.y);
|
||||
gc.stroke();
|
||||
gc.new_dash_pattern();
|
||||
gc.new_path();
|
||||
gc.move_to(p.x, p.y);
|
||||
}
|
||||
}
|
||||
gc.stroke();
|
||||
gc.new_dash_pattern();
|
||||
}
|
||||
|
||||
fn draw_fan_2(gc: &mut Vec<Draw>, space: &impl Metric, grid: &Grid, 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, grid, 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_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 Coil {
|
||||
fn render(&self, gc: &mut Vec<Draw>) {
|
||||
gc.new_path();
|
||||
gc.circle(0.0, 0.0, self.r + self.w + self.m);
|
||||
gc.circle(0.0, 0.0, self.r + self.w);
|
||||
gc.circle(0.0, 0.0, self.r - self.w);
|
||||
gc.circle(0.0, 0.0, self.r - self.w - self.m);
|
||||
gc.winding_rule(WindingRule::EvenOdd);
|
||||
gc.fill_color(Color::Rgba(0.8, 0.8, 0.8, 1.0));
|
||||
gc.fill();
|
||||
gc.line_width(0.5);
|
||||
gc.stroke_color(Color::Rgba(1.0, 0.5, 0.0, 1.0));
|
||||
draw_fan(gc, self, vec2(-500.0, 0.0), vec2(1.0, 0.0), 1.0);
|
||||
gc.stroke_color(Color::Rgba(0.0, 0.5, 1.0, 1.0));
|
||||
draw_fan(gc, self, vec2(0.0, self.r), vec2(1.0, 0.0), 1.0);
|
||||
}
|
||||
}
|
||||
|
||||
impl Renderable for Rect {
|
||||
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();
|
||||
}
|
||||
}
|
||||
|
||||
struct Coil {
|
||||
m: f32,
|
||||
scale: f32,
|
||||
r: f32,
|
||||
w: f32,
|
||||
}
|
||||
|
||||
impl Metric for Coil {
|
||||
fn halfmetric(&self, pos: Vec2) -> Decomp2 {
|
||||
let r = pos.length();
|
||||
let dir = pos.normalize();
|
||||
|
||||
let s = smoothbox(r, vec2(self.r - self.w, self.r + self.w), self.m);
|
||||
let t = 1.0.lerp(self.r / r / self.scale, s);
|
||||
|
||||
Decomp2 {
|
||||
ortho: Mat2::from_cols_array(&[
|
||||
dir.x, -dir.y,
|
||||
dir.y, dir.x,
|
||||
]),
|
||||
diag: vec2(1.0, t),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Debug, PartialEq)]
|
||||
struct Ray {
|
||||
pos: Vec2,
|
||||
dir: Vec2,
|
||||
}
|
||||
|
||||
mod basic_shapes {
|
||||
use glam::{Vec2, vec2};
|
||||
use crate::{Ray, shape};
|
||||
|
||||
pub 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() }
|
||||
}
|
||||
}
|
||||
|
||||
impl shape::Shape for Rect {
|
||||
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)
|
||||
}
|
||||
|
||||
fn visualise(&self) -> Vec<Vec2> {
|
||||
vec![vec2(-self.size.x, -self.size.y), vec2(self.size.x, -self.size.y), vec2(self.size.x, self.size.y), vec2(-self.size.x, self.size.y)]
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
use crate::shape::Shape;
|
||||
|
||||
#[test]
|
||||
fn test_rect() {
|
||||
assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 5.0) }), Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 5.0) });
|
||||
assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(-4.0, 5.0) }), Ray { pos: vec2(-2.0, 3.0), dir: vec2(4.0, 5.0) });
|
||||
assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, -5.0) }), Ray { pos: vec2(2.0, -3.0), dir: vec2(4.0, 5.0) });
|
||||
assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 0.0) }), Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 0.0) });
|
||||
|
||||
let r = Rect { size: vec2(2.0, 3.0) };
|
||||
|
||||
assert_eq!(r.trace_into(Ray { pos: vec2(3.0, 3.0), dir: vec2(1.0, 1.0) }), None);
|
||||
assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 2.0), dir: vec2(1.0, 0.0) }), Some(1.0));
|
||||
assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 2.0), dir: vec2(-1.0, 0.0) }), None);
|
||||
assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 1.0), dir: vec2(2.0, 2.0) }), Some(0.5));
|
||||
assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 2.1), dir: vec2(2.0, 2.0) }), None);
|
||||
|
||||
assert_eq!(r.trace_into(Ray { pos: vec2(2.0, 3.0), dir: vec2(1.0, 1.0) }), None);
|
||||
assert_eq!(r.trace_into(Ray { pos: vec2(-2.0, 3.0), dir: vec2(-1.0, 1.0) }), None);
|
||||
assert_eq!(r.trace_into(Ray { pos: vec2(2.0, 3.0), dir: vec2(-1.0, -1.0) }), Some(0.0));
|
||||
assert_eq!(r.trace_into(Ray { pos: vec2(2.0, -3.0), dir: vec2(-1.0, 1.0) }), Some(0.0));
|
||||
|
||||
assert_eq!(r.trace_out_of(Ray { pos: vec2(0.0, 0.0), dir: vec2(1.0, 1.0) }), Some(2.0));
|
||||
assert_eq!(r.trace_out_of(Ray { pos: vec2(0.0, 0.0), dir: vec2(0.0, 1.0) }), Some(3.0));
|
||||
assert_eq!(r.trace_out_of(Ray { pos: vec2(0.0, 1.0), dir: vec2(0.0, -1.0) }), Some(4.0));
|
||||
assert_eq!(r.trace_out_of(Ray { pos: vec2(1.0, 1.0), dir: vec2(0.0, -1.0) }), Some(4.0));
|
||||
assert_eq!(r.trace_out_of(Ray { pos: vec2(2.0, 3.0), dir: vec2(1.0, 1.0) }), Some(0.0));
|
||||
}
|
||||
}
|
||||
|
||||
mod shape {
|
||||
use glam::{Affine2, Vec2};
|
||||
use crate::Ray;
|
||||
|
||||
pub trait Shape {
|
||||
fn is_inside(&self, pt: Vec2) -> bool;
|
||||
|
||||
/// Ищет ближайшее пересечение луча с границей в направлении внутрь контура. Возвращает расстояние (в ray.dir).
|
||||
fn trace_into(&self, ray: Ray) -> Option<f32>;
|
||||
/// Ищет ближайшее пересечение луча с границей в направлении вовне контура. Возвращает расстояние (в ray.dir).
|
||||
fn trace_out_of(&self, ray: Ray) -> Option<f32>;
|
||||
|
||||
/// Возвращает визуальное представление контура, для отладки.
|
||||
fn visualise(&self) -> Vec<Vec2>;
|
||||
}
|
||||
|
||||
pub struct MovedShape<S: Shape> {
|
||||
pub shape: S,
|
||||
pub transform: Affine2, // transform(координаты контура) = 0
|
||||
}
|
||||
|
||||
impl<S: Shape> MovedShape<S> {
|
||||
fn pt_to_inner(&self, pt: Vec2) -> Vec2 {
|
||||
self.transform.transform_point2(pt)
|
||||
}
|
||||
fn pt_to_outer(&self, pt: Vec2) -> Vec2 {
|
||||
self.transform.inverse().transform_point2(pt)
|
||||
}
|
||||
fn vec_to_inner(&self, vec: Vec2) -> Vec2 {
|
||||
self.transform.transform_vector2(vec)
|
||||
}
|
||||
fn vec_to_outer(&self, vec: Vec2) -> Vec2 {
|
||||
self.transform.inverse().transform_vector2(vec)
|
||||
}
|
||||
fn ray_to_inner(&self, ray: Ray) -> Ray {
|
||||
Ray { pos: self.pt_to_inner(ray.pos), dir: self.vec_to_inner(ray.dir) }
|
||||
}
|
||||
fn ray_to_outer(&self, ray: Ray) -> Ray {
|
||||
Ray { pos: self.pt_to_outer(ray.pos), dir: self.vec_to_outer(ray.dir) }
|
||||
}
|
||||
}
|
||||
|
||||
impl<S: Shape> Shape for MovedShape<S> {
|
||||
fn is_inside(&self, pt: Vec2) -> bool {
|
||||
self.shape.is_inside(self.pt_to_inner(pt))
|
||||
}
|
||||
fn trace_into(&self, ray: Ray) -> Option<f32> {
|
||||
self.shape.trace_into(self.ray_to_inner(ray))
|
||||
}
|
||||
fn trace_out_of(&self, ray: Ray) -> Option<f32> {
|
||||
self.shape.trace_out_of(self.ray_to_inner(ray))
|
||||
}
|
||||
fn visualise(&self) -> Vec<Vec2> {
|
||||
self.shape.visualise().iter().map(|pt| self.pt_to_outer(*pt)).collect()
|
||||
}
|
||||
}
|
||||
|
||||
pub struct InvertedShape<S: Shape> {
|
||||
pub shape: S,
|
||||
}
|
||||
|
||||
impl<S: Shape> Shape for InvertedShape<S> {
|
||||
fn is_inside(&self, pt: Vec2) -> bool {
|
||||
!self.shape.is_inside(pt)
|
||||
}
|
||||
fn trace_into(&self, ray: Ray) -> Option<f32> {
|
||||
self.shape.trace_out_of(ray)
|
||||
}
|
||||
fn trace_out_of(&self, ray: Ray) -> Option<f32> {
|
||||
self.shape.trace_into(ray)
|
||||
}
|
||||
fn visualise(&self) -> Vec<Vec2> {
|
||||
self.shape.visualise()
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
struct Grid {
|
||||
hlines: Vec<f32>,
|
||||
vlines: Vec<f32>,
|
||||
cells: HashMap<(usize, usize), Box<dyn FlatCell>>,
|
||||
}
|
||||
|
||||
trait FlatCell: std::fmt::Debug {
|
||||
fn pos_to_global(&self, pos: Vec2) -> Vec2;
|
||||
fn pos_to_local(&self, pos: Vec2) -> Vec2;
|
||||
fn ray_to_global(&self, ray: Ray) -> Ray;
|
||||
fn ray_to_local(&self, ray: Ray) -> Ray;
|
||||
|
||||
fn is_inside(&self, pos: Vec2) -> bool {
|
||||
let bnd = self.local_bounds();
|
||||
pos.cmpge(bnd.0).all() && pos.cmple(bnd.1).all()
|
||||
}
|
||||
fn local_bounds(&self) -> (Vec2, Vec2);
|
||||
|
||||
fn to_boundary(&self, ray: Ray) -> Option<Ray> {
|
||||
assert!(self.is_inside(ray.pos));
|
||||
let sgn = ray.dir.signum();
|
||||
let p = ray.pos * sgn;
|
||||
let v = ray.dir * sgn;
|
||||
let mut bnd = self.local_bounds();
|
||||
if sgn.x < 0.0 {
|
||||
(bnd.0.x, bnd.1.x) = (-bnd.1.x, -bnd.0.x);
|
||||
}
|
||||
if sgn.y < 0.0 {
|
||||
(bnd.0.y, bnd.1.y) = (-bnd.1.y, -bnd.0.y);
|
||||
}
|
||||
let t = if (bnd.1.x - p.x) * v.y <= (bnd.1.y - p.y) * v.x {
|
||||
(bnd.1.x - p.x) / v.x
|
||||
} else {
|
||||
(bnd.1.y - p.y) / v.y
|
||||
};
|
||||
if t <= 1000.0 {
|
||||
Some(Ray { pos: sgn * (p + v * t), dir: sgn * v })
|
||||
} else {
|
||||
None
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Debug)]
|
||||
struct FlatRect {
|
||||
min: Vec2,
|
||||
max: Vec2,
|
||||
}
|
||||
|
||||
impl FlatCell for FlatRect {
|
||||
fn pos_to_global(&self, pos: Vec2) -> Vec2 { pos }
|
||||
fn pos_to_local(&self, pos: Vec2) -> Vec2 { pos }
|
||||
fn ray_to_global(&self, ray: Ray) -> Ray { ray }
|
||||
fn ray_to_local(&self, ray: Ray) -> Ray { ray }
|
||||
fn local_bounds(&self) -> (Vec2, Vec2) { (self.min, self.max) }
|
||||
}
|
||||
|
||||
#[derive(Debug)]
|
||||
struct RectInside {
|
||||
rect: Rect,
|
||||
}
|
||||
|
||||
impl FlatCell for RectInside {
|
||||
fn pos_to_global(&self, pos: Vec2) -> Vec2 {
|
||||
vec2(pos.x, self.rect.x(pos.y))
|
||||
}
|
||||
|
||||
fn pos_to_local(&self, pos: Vec2) -> Vec2 {
|
||||
vec2(pos.x, self.rect.u(pos.y))
|
||||
}
|
||||
|
||||
fn ray_to_global(&self, ray: Ray) -> Ray {
|
||||
Ray {
|
||||
pos: self.pos_to_global(ray.pos),
|
||||
dir: vec2(ray.dir.x, self.rect.dx(ray.pos.y, ray.dir.y)),
|
||||
}
|
||||
}
|
||||
|
||||
fn ray_to_local(&self, ray: Ray) -> Ray {
|
||||
Ray {
|
||||
pos: self.pos_to_local(ray.pos),
|
||||
dir: vec2(ray.dir.x, self.rect.du(ray.pos.y, ray.dir.y)),
|
||||
}
|
||||
}
|
||||
|
||||
fn local_bounds(&self) -> (Vec2, Vec2) {
|
||||
(vec2(-self.rect.inner_radius, -self.rect.internal_halflength), vec2(self.rect.inner_radius, self.rect.internal_halflength))
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Copy, Clone, Debug)]
|
||||
struct Rect {
|
||||
outer_radius: f32,
|
||||
inner_radius: f32,
|
||||
external_halflength: f32,
|
||||
internal_halflength: f32,
|
||||
}
|
||||
|
||||
impl Rect {
|
||||
fn γ(&self) -> f32 { self.external_halflength / self.internal_halflength }
|
||||
fn ri(&self) -> f32 { self.internal_halflength }
|
||||
fn re(&self) -> f32 { self.external_halflength }
|
||||
fn a(&self) -> f32 { (1.0 - self.γ()) / self.ri() }
|
||||
fn b(&self) -> f32 { 2.0 * self.γ() - 1.0 }
|
||||
|
||||
fn root(&self, x: f32) -> f32 { ((2.0 * self.γ() - 1.0).powi(2) + 4.0 * (1.0 - self.γ()) * x / self.ri()).sqrt() }
|
||||
fn d(&self, u: f32) -> f32 { 2.0 * self.a() * u + self.b() }
|
||||
pub fn x(&self, u: f32) -> f32 { (self.a() * u.abs() + self.b()) * u }
|
||||
pub fn u(&self, x: f32) -> f32 { 0.5 * self.ri() * (1.0 - 2.0 * self.γ() + self.root(x.abs())) / (1.0 - self.γ()) * x.signum() }
|
||||
pub fn dx(&self, u: f32, du: f32) -> f32 { du * self.d(u.abs()) }
|
||||
pub fn du(&self, x: f32, dx: f32) -> f32 { dx / self.root(x.abs()) }
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_rect() {
|
||||
let r = Rect {
|
||||
outer_radius: 50.0,
|
||||
inner_radius: 20.0,
|
||||
external_halflength: 100.0,
|
||||
internal_halflength: 10.0,
|
||||
};
|
||||
assert_abs_diff_eq!(r.x(r.internal_halflength), r.external_halflength, epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.x(-r.internal_halflength), -r.external_halflength, epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.u(r.external_halflength), r.internal_halflength, epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.u(-r.external_halflength), -r.internal_halflength, epsilon = 1.0e-5);
|
||||
|
||||
assert_abs_diff_eq!(r.dx(r.internal_halflength, 3.0), 3.0, epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.dx(-r.internal_halflength, 3.0), 3.0, epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.du(r.external_halflength, 3.0), 3.0, epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.du(-r.external_halflength, 3.0), 3.0, epsilon = 1.0e-5);
|
||||
|
||||
assert_abs_diff_eq!(r.u(r.x(1.0)), 1.0, epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.u(r.x(5.0)), 5.0, epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.u(r.x(-5.0)), -5.0, epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.u(r.x(10.0)), 10.0, epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.x(r.u(10.0)), 10.0, epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.x(r.u(50.0)), 50.0, epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.x(r.u(-50.0)), -50.0, epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.x(r.u(100.0)), 100.0, epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.d(r.u(10.0)), r.root(10.0), epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.d(r.u(50.0)), r.root(50.0), epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.d(r.u(100.0)), r.root(100.0), epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.du(10.0, r.dx(r.u(10.0), 3.0)), 3.0, epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.du(50.0, r.dx(r.u(50.0), 3.0)), 3.0, epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.du(-50.0, r.dx(r.u(-50.0), 3.0)), 3.0, epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.du(100.0, r.dx(r.u(100.0), 3.0)), 3.0, epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.dx(1.0, r.du(r.x(1.0), 3.0)), 3.0, epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.dx(5.0, r.du(r.x(5.0), 3.0)), 3.0, epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.dx(-5.0, r.du(r.x(-5.0), 3.0)), 3.0, epsilon = 1.0e-5);
|
||||
assert_abs_diff_eq!(r.dx(10.0, r.du(r.x(10.0), 3.0)), 3.0, epsilon = 1.0e-5);
|
||||
}
|
||||
|
||||
impl Metric for Rect {
|
||||
fn halfmetric(&self, pos: Vec2) -> Decomp2 {
|
||||
let x = pos.y.abs();
|
||||
let sx = ((pos.x.abs() - self.inner_radius) / (self.outer_radius - self.inner_radius)).clamp(0.0, 1.0);
|
||||
let sy = if x <= self.external_halflength { 1.0 / self.root(x) } else { 1.0 };
|
||||
assert!(sx.is_finite());
|
||||
assert!(sy.is_finite());
|
||||
assert!(sy > 0.0);
|
||||
Decomp2 {
|
||||
ortho: Mat2::IDENTITY,
|
||||
diag: vec2(1.0, sy.lerp(1.0, sx)),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
mod riemann {
|
||||
use glam::*;
|
||||
|
||||
pub struct Decomp2 {
|
||||
pub ortho: Mat2,
|
||||
pub diag: Vec2,
|
||||
}
|
||||
|
||||
impl Decomp2 {
|
||||
fn square(&self) -> Self {
|
||||
Self {
|
||||
ortho: self.ortho,
|
||||
diag: self.diag * self.diag,
|
||||
}
|
||||
}
|
||||
|
||||
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
|
||||
}
|
||||
}
|
||||
|
||||
type Tens2 = [Mat2; 2];
|
||||
|
||||
pub trait Metric {
|
||||
fn halfmetric(&self, pos: Vec2) -> Decomp2;
|
||||
|
||||
fn metric(&self, pos: Vec2) -> Mat2 {
|
||||
self.halfmetric(pos).square().into()
|
||||
}
|
||||
|
||||
fn invmetric(&self, pos: Vec2) -> Mat2 {
|
||||
self.halfmetric(pos).square().inverse().into()
|
||||
}
|
||||
|
||||
fn dmetric(&self, pos: Vec2) -> Tens2 {
|
||||
part_deriv(|p| self.metric(p), pos, 1.0e-3)
|
||||
}
|
||||
|
||||
fn length(&self, at: Vec2, v: Vec2) -> f32 {
|
||||
v.dot(self.metric(at) * v).sqrt()
|
||||
}
|
||||
|
||||
fn normalize(&self, at: Vec2, v: Vec2) -> Vec2 {
|
||||
v / self.length(at, v)
|
||||
}
|
||||
|
||||
fn globalize(&self, at: Vec2, v: Vec2) -> Vec2 {
|
||||
Mat2::from(self.halfmetric(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 = -convolute(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(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.invmetric(pos); // с верхними индексами
|
||||
let d = space.dmetric(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)),
|
||||
]
|
||||
}
|
||||
|
||||
pub fn convolute(t: Tens2, v: Vec2) -> Vec2 {
|
||||
vec2(
|
||||
v.dot(t[0] * v),
|
||||
v.dot(t[1] * 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);
|
||||
}
|
||||
}
|
||||
|
||||
fn smoothstep(x: f32) -> f32 {
|
||||
return 3.0 * x * x - 2.0 * x * x * x;
|
||||
}
|
||||
|
||||
fn smoothbox(val: f32, range: Vec2, pad: f32) -> f32 {
|
||||
let slope1 = 1.0 + (val - range.x) / pad;
|
||||
let slope2 = 1.0 - (val - range.y) / pad;
|
||||
let lin = slope1.min(slope2);
|
||||
smoothstep(lin.clamp(0.0, 1.0))
|
||||
}
|
||||
|
|
@ -1,231 +0,0 @@
|
|||
use std::f32::consts::{FRAC_PI_2, PI};
|
||||
|
||||
use flo_canvas::*;
|
||||
use flo_draw::*;
|
||||
use glam::*;
|
||||
|
||||
use refraction::ifaces::{DebugTraceable, Traceable};
|
||||
use refraction::tube::metric::Tube;
|
||||
use refraction::tube::Space;
|
||||
use refraction::types::{Location, Object, Ray};
|
||||
use refraction::utils::put_object;
|
||||
use refraction::DT;
|
||||
|
||||
fn draw_loop(gc: &mut Vec<Draw>, mut pts: impl Iterator<Item = Vec3>) {
|
||||
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: put_object(
|
||||
&tube,
|
||||
vec3(0.0, y * tube.external_halflength, 0.0),
|
||||
Mat3::from_mat2(Mat2::from_angle(y)),
|
||||
),
|
||||
r: 20.0,
|
||||
})
|
||||
.collect();
|
||||
let space = Space { tube, objs };
|
||||
let cam1 = put_object(&space.tube, vec3(-500., 0., 0.), Mat3::IDENTITY);
|
||||
let cam2 = put_object(
|
||||
&space.tube,
|
||||
vec3(-2.5 * tube.outer_radius, 1.25 * tube.external_halflength, 0.),
|
||||
mat3(vec3(1., -1., 0.), vec3(1., 1., 0.), vec3(0., 0., 1.)),
|
||||
);
|
||||
let cam3 = put_object(
|
||||
&space.tube,
|
||||
vec3(0.25 * tube.inner_radius, 0.25 * tube.external_halflength, 0.),
|
||||
mat3(vec3(0., -1., 0.), vec3(1., 0., 0.), vec3(0., 0., 1.)),
|
||||
);
|
||||
|
||||
gc.canvas_height(500.0);
|
||||
gc.transform(Transform2D::rotate(FRAC_PI_2));
|
||||
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(0.0, 0.8, 1.0, 1.0));
|
||||
draw_fan_2(gc, &space, cam3, 1.0);
|
||||
gc.stroke_color(Color::Rgba(0.5, 1.0, 0.0, 1.0));
|
||||
draw_fan_2(gc, &space, cam2, 1.0);
|
||||
gc.stroke_color(Color::Rgba(1.0, 0.5, 0.0, 1.0));
|
||||
draw_fan_2(gc, &space, cam1, 1.0);
|
||||
|
||||
draw_track(gc, &space, vec2(-500.0, 0.0), vec2(1.0, 0.2));
|
||||
draw_track(gc, &space, vec2(-500.0, 0.0), vec2(1.0, 0.5));
|
||||
draw_track(
|
||||
gc,
|
||||
&space,
|
||||
vec2(-0.5 * tube.inner_radius, -1.25 * tube.external_halflength),
|
||||
vec2(0.1, 1.0),
|
||||
);
|
||||
|
||||
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 = vec3(dir.x, dir.y, 0.);
|
||||
let dir = obj.loc.rot * dir;
|
||||
pos + dir
|
||||
}),
|
||||
);
|
||||
gc.stroke_color(Color::Rgba(0.0, 0.5, 1.0, 0.5));
|
||||
draw_loop(
|
||||
gc,
|
||||
itertools_num::linspace(0.0, 2.0 * PI, circle_segments)
|
||||
.skip(1)
|
||||
.map(|φ| {
|
||||
let dir = Vec2::from_angle(φ) * obj.r;
|
||||
let dir = vec3(dir.x, dir.y, 0.);
|
||||
let dir = obj.loc.rot * dir;
|
||||
space.trace_step(Ray { pos, dir }).pos
|
||||
}),
|
||||
);
|
||||
gc.stroke_color(Color::Rgba(0.5, 0.0, 1.0, 1.0));
|
||||
draw_loop(
|
||||
gc,
|
||||
itertools_num::linspace(0.0, 2.0 * PI, circle_segments)
|
||||
.skip(1)
|
||||
.map(|φ| {
|
||||
let n = obj.r.floor();
|
||||
let d = obj.r / n;
|
||||
let dir = Vec2::from_angle(φ);
|
||||
let dir = vec3(dir.x, dir.y, 0.);
|
||||
let dir = obj.loc.rot * dir * d;
|
||||
space.trace_iter(Ray { pos, dir }).nth(n as usize).unwrap().pos
|
||||
}),
|
||||
);
|
||||
}
|
||||
});
|
||||
});
|
||||
}
|
||||
|
||||
#[allow(dead_code)]
|
||||
fn draw_cross(gc: &mut Vec<Draw>, pos: Vec2, r: f32) {
|
||||
gc.move_to(pos.x - r, pos.y - r);
|
||||
gc.line_to(pos.x + r, pos.y + r);
|
||||
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, camera: Location, dir: Vec3) {
|
||||
let pos = vec3(0., 0., 0.);
|
||||
let (hits, path) = space.trace_dbg(camera, Ray { pos, dir });
|
||||
let hits2 = space.trace(camera, Ray { pos, dir });
|
||||
for (a, b) in hits.into_iter().zip(hits2.into_iter()) {
|
||||
assert_eq!(a.id, b.id);
|
||||
assert_eq!(a.pos, b.pos);
|
||||
assert_eq!(a.rel, b.rel);
|
||||
}
|
||||
|
||||
gc.new_path();
|
||||
gc.move_to(pos.x, pos.y);
|
||||
for pt in &path.points[1..] {
|
||||
gc.line_to(pt.pos.x, pt.pos.y);
|
||||
}
|
||||
let end_pos = *path.points.last().expect("the starting point is always in the path");
|
||||
let dir_pos = end_pos.forward(1000. / DT).pos;
|
||||
gc.line_to(dir_pos.x, dir_pos.y);
|
||||
gc.stroke();
|
||||
}
|
||||
|
||||
fn draw_fan_2(gc: &mut Vec<Draw>, space: &Space, camera: Location, spread: f32) {
|
||||
for y in itertools_num::linspace(-spread, spread, 101) {
|
||||
draw_ray_2(gc, space, camera, vec3(1., y, 0.));
|
||||
}
|
||||
}
|
||||
|
||||
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: vec3(start.x, start.y, 0.),
|
||||
rot: mat3(vec3(dir.x, dir.y, 0.), vec3(-dir.y, dir.x, 0.), vec3(0., 0., 1.)),
|
||||
};
|
||||
let v = vec3(1., 0., 0.);
|
||||
let mut draw = |loc: &Location| {
|
||||
let p = loc.pos;
|
||||
let ax = p + loc.rot.x_axis * SCALE;
|
||||
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);
|
||||
}
|
||||
}
|
||||
|
||||
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,153 +0,0 @@
|
|||
use glam::*;
|
||||
use refraction::mesh_loader::load_mesh;
|
||||
use refraction::mesh_tracer::{trace_to_mesh, Mesh};
|
||||
use show_image::event::{ElementState, VirtualKeyCode, WindowEvent};
|
||||
use show_image::{exit, ImageInfo, ImageView, WindowOptions};
|
||||
use std::env;
|
||||
use std::error::Error;
|
||||
use std::f32::consts::PI;
|
||||
use std::fs::File;
|
||||
use std::io::BufReader;
|
||||
use std::sync::atomic::{AtomicUsize, Ordering::Relaxed};
|
||||
|
||||
const W: i32 = 320;
|
||||
const H: i32 = 240;
|
||||
|
||||
#[derive(Copy, Clone)]
|
||||
struct Color(u8, u8, u8);
|
||||
|
||||
struct Image {
|
||||
w: i32,
|
||||
h: i32,
|
||||
data: Vec<u8>,
|
||||
}
|
||||
|
||||
impl Image {
|
||||
fn data(&self) -> &[u8] {
|
||||
self.data.as_slice()
|
||||
}
|
||||
|
||||
fn put_pixel(&mut self, x: i32, y: i32, color: Color) {
|
||||
if x < 0 || x >= self.w || y < 0 || y > self.h {
|
||||
return;
|
||||
}
|
||||
let index = 3 * (x + self.w * y) as usize;
|
||||
self.data[index] = color.0;
|
||||
self.data[index + 1] = color.1;
|
||||
self.data[index + 2] = color.2;
|
||||
}
|
||||
}
|
||||
|
||||
fn ypr_to_mat(ypr: Vec3) -> Mat3 {
|
||||
let Vec3 {
|
||||
x: yaw,
|
||||
y: pitch,
|
||||
z: roll,
|
||||
} = ypr;
|
||||
let m_roll = mat3(
|
||||
vec3(roll.cos(), roll.sin(), 0.0),
|
||||
vec3(-roll.sin(), roll.cos(), 0.0),
|
||||
vec3(0.0, 0.0, 1.0),
|
||||
);
|
||||
let m_yaw = mat3(
|
||||
vec3(yaw.cos(), 0.0, yaw.sin()),
|
||||
vec3(0.0, 1.0, 0.0),
|
||||
vec3(-yaw.sin(), 0.0, yaw.cos()),
|
||||
);
|
||||
let m_pitch = mat3(
|
||||
vec3(1.0, 0.0, 0.0),
|
||||
vec3(0.0, pitch.cos(), -pitch.sin()),
|
||||
vec3(0.0, pitch.sin(), pitch.cos()),
|
||||
);
|
||||
m_roll * m_pitch * m_yaw
|
||||
}
|
||||
|
||||
fn render(mesh: &Mesh, camera: impl Fn(Vec2) -> (Vec3, Vec3)) -> Image {
|
||||
let bkg = vec3(0.0, 0.0, 0.0);
|
||||
let mut img = Image {
|
||||
w: W,
|
||||
h: H,
|
||||
data: vec![0; (3 * W * H) as usize],
|
||||
};
|
||||
let img_size = vec2(W as f32, H as f32);
|
||||
for y in 0..H {
|
||||
for x in 0..W {
|
||||
let img_coords = vec2(x as f32, y as f32);
|
||||
let off = (img_coords - img_size * 0.5) / img_size.y;
|
||||
let (base, ray) = camera(off);
|
||||
let color = if let Some(r) = trace_to_mesh(mesh, base, ray.normalize()) {
|
||||
// to_vec3(0.45) * dot(r.normal, normalize(vec3(-1.0, 1.0, -1.0))) + 0.50
|
||||
r.normal * 0.45 + 0.50
|
||||
} else {
|
||||
bkg
|
||||
};
|
||||
let color = (color * 255.0).as_ivec3().clamp(IVec3::splat(0), IVec3::splat(255));
|
||||
img.put_pixel(x, y, Color(color.x as u8, color.y as u8, color.z as u8));
|
||||
}
|
||||
}
|
||||
img
|
||||
}
|
||||
|
||||
fn persp(dist: f32, off: Vec2) -> (Vec3, Vec3) {
|
||||
(vec3(0., 0., -dist), vec3(off.x, off.y, dist))
|
||||
}
|
||||
|
||||
fn ortho(dist: f32, off: Vec2) -> (Vec3, Vec3) {
|
||||
(vec3(off.x, off.y, -dist), vec3(0., 0., 1.))
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_projs() {
|
||||
fn check(f: fn(dist: f32, off: Vec2) -> (Vec3, Vec3), x: f32, y: f32, z: f32) {
|
||||
let (base, ray) = f(z, vec2(x, y));
|
||||
let at_dist = base + ray * (z / ray.z);
|
||||
assert_eq!(at_dist, vec3(x, y, 0.));
|
||||
}
|
||||
check(persp, 1., 2., 3.);
|
||||
check(ortho, 1., 2., 3.);
|
||||
check(persp, 5., 3., 7.);
|
||||
check(ortho, 9., 1., 8.);
|
||||
}
|
||||
|
||||
// add_event_handler wants 'static + Send. Let it be so.
|
||||
static PROJ_INDEX: AtomicUsize = AtomicUsize::new(0);
|
||||
static PROJS: [fn(dist: f32, off: Vec2) -> (Vec3, Vec3); 2] = [persp, ortho];
|
||||
|
||||
#[show_image::main]
|
||||
fn main() -> Result<(), Box<dyn Error>> {
|
||||
let args: Vec<String> = env::args().collect();
|
||||
if args.len() != 2 {
|
||||
println!("Usage: {} path/to/model.obj", args[0]);
|
||||
exit(1);
|
||||
}
|
||||
let mesh = {
|
||||
let f = File::open(&args[1])?;
|
||||
let mut f = BufReader::new(f);
|
||||
load_mesh(&mut f)?
|
||||
};
|
||||
let window = show_image::create_window("Raytracing", WindowOptions::default())?;
|
||||
window.add_event_handler(|_wnd, ev, _ctl| {
|
||||
if let WindowEvent::KeyboardInput(ev) = ev {
|
||||
if ev.input.state != ElementState::Pressed {
|
||||
return;
|
||||
}
|
||||
if let Some(VirtualKeyCode::Tab) = ev.input.key_code {
|
||||
PROJ_INDEX.store((PROJ_INDEX.load(Relaxed) + 1) % PROJS.len(), Relaxed);
|
||||
}
|
||||
}
|
||||
})?;
|
||||
loop {
|
||||
for phi in 0..360 {
|
||||
let proj = PROJS[PROJ_INDEX.load(Relaxed)];
|
||||
let m_view = ypr_to_mat(vec3((135.0 + phi as f32) * PI / 180.0, -30.0 * PI / 180.0, 0.0f32));
|
||||
let m_camera = m_view.transpose();
|
||||
let img = render(mesh.as_slice(), |off| {
|
||||
let (base, ray) = proj(40., 20. * off);
|
||||
(m_camera * base, m_camera * ray)
|
||||
});
|
||||
|
||||
let image = ImageView::new(ImageInfo::rgb8(W as u32, H as u32), img.data());
|
||||
window.set_image("image", image)?;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -1,123 +0,0 @@
|
|||
use std::mem;
|
||||
|
||||
use glam::{mat4, vec3, vec4, Mat4, Vec2, Vec3};
|
||||
|
||||
#[repr(C)]
|
||||
#[derive(Copy, Clone, bytemuck::Pod, bytemuck::Zeroable)]
|
||||
struct CameraUniform {
|
||||
mvp: [[f32; 4]; 4],
|
||||
scale: [f32; 2],
|
||||
pad: [u32; 2],
|
||||
}
|
||||
|
||||
pub struct Camera {
|
||||
buf: wgpu::Buffer,
|
||||
bind: wgpu::BindGroup,
|
||||
layout: wgpu::BindGroupLayout,
|
||||
}
|
||||
|
||||
impl Camera {
|
||||
pub fn new(device: &wgpu::Device) -> Camera {
|
||||
let buf = device.create_buffer(&wgpu::BufferDescriptor {
|
||||
label: Some("Camera Buffer"),
|
||||
usage: wgpu::BufferUsages::UNIFORM | wgpu::BufferUsages::COPY_DST,
|
||||
size: mem::size_of::<CameraUniform>() as u64,
|
||||
mapped_at_creation: false,
|
||||
});
|
||||
let layout = device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor {
|
||||
entries: &[wgpu::BindGroupLayoutEntry {
|
||||
binding: 0,
|
||||
visibility: wgpu::ShaderStages::VERTEX,
|
||||
ty: wgpu::BindingType::Buffer {
|
||||
ty: wgpu::BufferBindingType::Uniform,
|
||||
has_dynamic_offset: false,
|
||||
min_binding_size: None,
|
||||
},
|
||||
count: None,
|
||||
}],
|
||||
label: Some("Camera BindGroupLayout"),
|
||||
});
|
||||
let bind = device.create_bind_group(&wgpu::BindGroupDescriptor {
|
||||
layout: &layout,
|
||||
entries: &[wgpu::BindGroupEntry {
|
||||
binding: 0,
|
||||
resource: buf.as_entire_binding(),
|
||||
}],
|
||||
label: Some("Camera BindGroup"),
|
||||
});
|
||||
Camera { buf, bind, layout }
|
||||
}
|
||||
|
||||
pub fn bind_group_layout(&self) -> &wgpu::BindGroupLayout {
|
||||
&self.layout
|
||||
}
|
||||
|
||||
pub fn bind_group(&self) -> &wgpu::BindGroup {
|
||||
&self.bind
|
||||
}
|
||||
|
||||
fn uniform(view_mtx: Mat4, view_size: Vec2) -> CameraUniform {
|
||||
const M: Mat4 = mat4(
|
||||
vec4(0., 0., 1., 0.),
|
||||
vec4(-1., 0., 0., 0.),
|
||||
vec4(0., 1., 0., 0.),
|
||||
vec4(0., 0., 0., 1.),
|
||||
);
|
||||
let size = view_size.normalize() * std::f32::consts::SQRT_2;
|
||||
let proj = make_proj_matrix(vec3(size.x, size.y, 2.), (1., (2f32).powi(16) + 1.)) * M;
|
||||
|
||||
let mvp = proj * view_mtx;
|
||||
CameraUniform {
|
||||
mvp: mvp.to_cols_array_2d(),
|
||||
scale: (1. / size).to_array(),
|
||||
pad: [0; 2],
|
||||
}
|
||||
}
|
||||
|
||||
pub fn set(&self, queue: &wgpu::Queue, view_mtx: Mat4, view_size: Vec2) {
|
||||
let uniform = Self::uniform(view_mtx, view_size);
|
||||
queue.write_buffer(&self.buf, 0, bytemuck::bytes_of(&uniform));
|
||||
}
|
||||
}
|
||||
|
||||
/// Make a projection matrix, assuming input coordinates are (right, up, forward).
|
||||
///
|
||||
/// `corner` is a vector that will be mapped to (x=1, y=1) after the perspective division.
|
||||
/// `zrange` is the Z range that will be mapped to z∈[-1, 1]. It has no other effect. Both ends have to be positive though.
|
||||
fn make_proj_matrix(corner: Vec3, zrange: (f32, f32)) -> Mat4 {
|
||||
let scale = 1.0 / corner;
|
||||
let zspan = zrange.1 - zrange.0;
|
||||
mat4(
|
||||
scale.x * vec4(1., 0., 0., 0.),
|
||||
scale.y * vec4(0., 1., 0., 0.),
|
||||
scale.z * vec4(0., 0., (zrange.0 + zrange.1) / zspan, 1.),
|
||||
scale.z * vec4(0., 0., -2. * zrange.0 * zrange.1 / zspan, 0.),
|
||||
)
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use approx::assert_abs_diff_eq;
|
||||
use glam::vec3;
|
||||
|
||||
#[test]
|
||||
fn test_proj_matrix() {
|
||||
let m = make_proj_matrix(vec3(2., 3., 4.), (0.5, 20.0));
|
||||
|
||||
let v = m * vec4(2., 3., 4., 1.);
|
||||
assert_abs_diff_eq!(v.x / v.w, 1.0);
|
||||
assert_abs_diff_eq!(v.y / v.w, 1.0);
|
||||
assert!(-v.w < v.z && v.z < v.w, "z out of range in {v}");
|
||||
|
||||
let v = m * vec4(2., 3., 0.5, 1.);
|
||||
assert_abs_diff_eq!(v.x / v.w, 8.0);
|
||||
assert_abs_diff_eq!(v.y / v.w, 8.0);
|
||||
assert_abs_diff_eq!(v.z / v.w, -1.0);
|
||||
|
||||
let v = m * vec4(2., 3., 20.0, 1.);
|
||||
assert_abs_diff_eq!(v.x / v.w, 0.2);
|
||||
assert_abs_diff_eq!(v.y / v.w, 0.2);
|
||||
assert_abs_diff_eq!(v.z / v.w, 1.0);
|
||||
}
|
||||
}
|
||||
|
|
@ -1,168 +0,0 @@
|
|||
use std::mem;
|
||||
|
||||
use bytemuck::{bytes_of, cast_slice, Pod, Zeroable};
|
||||
use glam::Vec3;
|
||||
use refraction::types::Ray;
|
||||
use wgpu::util::DeviceExt as _;
|
||||
|
||||
#[repr(C)]
|
||||
#[derive(Copy, Clone, Debug, Pod, Zeroable)]
|
||||
struct Vertex {
|
||||
pub position: [f32; 3],
|
||||
pub tangent: [f32; 3],
|
||||
}
|
||||
|
||||
#[repr(C)]
|
||||
#[derive(Copy, Clone, Pod, Zeroable)]
|
||||
struct PushConsts {
|
||||
pub color: [f32; 3],
|
||||
pub _pad: f32,
|
||||
}
|
||||
|
||||
#[derive(Copy, Clone)]
|
||||
pub struct Attrs {
|
||||
pub color: Vec3,
|
||||
}
|
||||
|
||||
impl Attrs {
|
||||
fn consts(&self) -> PushConsts {
|
||||
PushConsts {
|
||||
color: self.color.to_array(),
|
||||
_pad: 0.,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
pub struct Line {
|
||||
consts: PushConsts,
|
||||
npoints: u32,
|
||||
buf: wgpu::Buffer,
|
||||
}
|
||||
|
||||
impl Line {
|
||||
pub fn new_strip(device: &wgpu::Device, attrs: Attrs, points: Vec<Ray>) -> Line {
|
||||
let data: Vec<Vertex> = points
|
||||
.into_iter()
|
||||
.map(|r| Vertex {
|
||||
position: r.pos.to_array(),
|
||||
tangent: r.dir.to_array(),
|
||||
})
|
||||
.collect();
|
||||
let buf = device.create_buffer_init(&wgpu::util::BufferInitDescriptor {
|
||||
label: Some("Vertex Buffer"),
|
||||
contents: cast_slice(&data),
|
||||
usage: wgpu::BufferUsages::VERTEX,
|
||||
});
|
||||
Line {
|
||||
consts: attrs.consts(),
|
||||
npoints: data.len() as u32,
|
||||
buf,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
pub struct LineRenderer {
|
||||
pipeline: wgpu::RenderPipeline,
|
||||
}
|
||||
|
||||
static SHADER: &'static str = include_str!("ray.wgsl");
|
||||
|
||||
impl LineRenderer {
|
||||
pub fn new(
|
||||
device: &wgpu::Device,
|
||||
cam_layout: &wgpu::BindGroupLayout,
|
||||
target_format: wgpu::TextureFormat,
|
||||
depth_stencil: Option<wgpu::DepthStencilState>,
|
||||
multisample: wgpu::MultisampleState,
|
||||
) -> LineRenderer {
|
||||
let shader = device.create_shader_module(wgpu::ShaderModuleDescriptor {
|
||||
label: Some("Line Shader"),
|
||||
source: wgpu::ShaderSource::Wgsl(SHADER.into()),
|
||||
});
|
||||
|
||||
let consts_range = wgpu::PushConstantRange {
|
||||
stages: wgpu::ShaderStages::VERTEX,
|
||||
range: 0..mem::size_of::<PushConsts>() as u32,
|
||||
};
|
||||
|
||||
let layout = device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor {
|
||||
label: Some("Line RenderPipelineLayout"),
|
||||
bind_group_layouts: &[cam_layout],
|
||||
push_constant_ranges: &[consts_range],
|
||||
});
|
||||
|
||||
let pipeline = device.create_render_pipeline(&wgpu::RenderPipelineDescriptor {
|
||||
label: Some("Line RenderPipeline"),
|
||||
layout: Some(&layout),
|
||||
vertex: wgpu::VertexState {
|
||||
module: &shader,
|
||||
entry_point: "vs_main",
|
||||
buffers: &[wgpu::VertexBufferLayout {
|
||||
array_stride: mem::size_of::<Vertex>() as wgpu::BufferAddress,
|
||||
step_mode: wgpu::VertexStepMode::Instance,
|
||||
attributes: &[
|
||||
wgpu::VertexAttribute {
|
||||
offset: mem::offset_of!(Vertex, position) as u64,
|
||||
shader_location: 0,
|
||||
format: wgpu::VertexFormat::Float32x3,
|
||||
},
|
||||
wgpu::VertexAttribute {
|
||||
offset: mem::offset_of!(Vertex, tangent) as u64,
|
||||
shader_location: 1,
|
||||
format: wgpu::VertexFormat::Float32x3,
|
||||
},
|
||||
wgpu::VertexAttribute {
|
||||
offset: (mem::size_of::<Vertex>() + mem::offset_of!(Vertex, position)) as u64,
|
||||
shader_location: 2,
|
||||
format: wgpu::VertexFormat::Float32x3,
|
||||
},
|
||||
wgpu::VertexAttribute {
|
||||
offset: (mem::size_of::<Vertex>() + mem::offset_of!(Vertex, tangent)) as u64,
|
||||
shader_location: 3,
|
||||
format: wgpu::VertexFormat::Float32x3,
|
||||
},
|
||||
],
|
||||
}],
|
||||
compilation_options: Default::default(),
|
||||
},
|
||||
fragment: Some(wgpu::FragmentState {
|
||||
module: &shader,
|
||||
entry_point: "fs_main",
|
||||
targets: &[Some(wgpu::ColorTargetState {
|
||||
format: target_format,
|
||||
blend: Some(wgpu::BlendState {
|
||||
color: wgpu::BlendComponent::OVER,
|
||||
alpha: wgpu::BlendComponent::OVER,
|
||||
}),
|
||||
write_mask: wgpu::ColorWrites::ALL,
|
||||
})],
|
||||
compilation_options: Default::default(),
|
||||
}),
|
||||
primitive: wgpu::PrimitiveState {
|
||||
topology: wgpu::PrimitiveTopology::TriangleStrip,
|
||||
..Default::default()
|
||||
},
|
||||
depth_stencil,
|
||||
multisample,
|
||||
multiview: None,
|
||||
cache: None,
|
||||
});
|
||||
|
||||
LineRenderer { pipeline }
|
||||
}
|
||||
|
||||
pub fn render<'a>(
|
||||
&self,
|
||||
pass: &mut wgpu::RenderPass,
|
||||
cam_bind: &wgpu::BindGroup,
|
||||
lines: impl Iterator<Item = &'a Line>,
|
||||
) {
|
||||
pass.set_pipeline(&self.pipeline);
|
||||
pass.set_bind_group(0, cam_bind, &[]);
|
||||
for line in lines {
|
||||
pass.set_push_constants(wgpu::ShaderStages::VERTEX, 0, bytes_of(&line.consts));
|
||||
pass.set_vertex_buffer(0, line.buf.slice(..));
|
||||
pass.draw(0..4, 0..line.npoints - 1);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -1,417 +0,0 @@
|
|||
use std::time::Instant;
|
||||
|
||||
use glam::{uvec2, vec3, Vec3};
|
||||
use winit::{
|
||||
event::*,
|
||||
event_loop::EventLoop,
|
||||
keyboard::{KeyCode, PhysicalKey},
|
||||
window::{Window, WindowBuilder},
|
||||
};
|
||||
|
||||
mod camera;
|
||||
mod lines;
|
||||
mod meshes;
|
||||
mod scene;
|
||||
mod viewport;
|
||||
|
||||
// The coordinate system:
|
||||
// * X: forward
|
||||
// * Y: left
|
||||
// * Z: up
|
||||
|
||||
fn prepare_scene(device: &wgpu::Device) -> (Vec<meshes::Mesh>, Vec<lines::Line>) {
|
||||
let (meshes, lines) = scene::build();
|
||||
let meshes = meshes
|
||||
.into_iter()
|
||||
.map(|mesh| meshes::Mesh::new_list(device, meshes::Attrs { color: mesh.color }, mesh.tris))
|
||||
.collect();
|
||||
let lines = lines
|
||||
.into_iter()
|
||||
.map(|line| lines::Line::new_strip(device, lines::Attrs { color: line.color }, line.pts))
|
||||
.collect();
|
||||
(meshes, lines)
|
||||
}
|
||||
|
||||
#[cfg(any())]
|
||||
mod camctl {
|
||||
use glam::{vec3, Mat4, Quat, Vec3};
|
||||
|
||||
pub struct CameraLocation {
|
||||
pos: Vec3,
|
||||
rot: Quat,
|
||||
}
|
||||
|
||||
impl CameraLocation {
|
||||
pub fn new() -> CameraLocation {
|
||||
let rot = Quat::from_euler(glam::EulerRot::ZYX, std::f32::consts::FRAC_PI_4, 0., 0.);
|
||||
let pos = rot * vec3(-200., 0., 50.);
|
||||
CameraLocation { pos, rot }
|
||||
}
|
||||
|
||||
pub fn view_mtx(&self) -> Mat4 {
|
||||
Mat4::from_quat(self.rot.inverse()) * Mat4::from_translation(-self.pos)
|
||||
}
|
||||
|
||||
pub fn move_rel(&mut self, offset: Vec3) {
|
||||
self.pos += self.rot * offset;
|
||||
}
|
||||
|
||||
pub fn rotate_rel_ypr(&mut self, ypr: Vec3) {
|
||||
self.rotate_rel_quat(Quat::from_euler(glam::EulerRot::ZYX, ypr.x, ypr.y, ypr.z));
|
||||
}
|
||||
|
||||
fn rotate_rel_quat(&mut self, rot: Quat) {
|
||||
self.rot *= rot;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
mod camctl {
|
||||
use glam::{vec3, Mat4, Quat, Vec3};
|
||||
|
||||
pub struct CameraLocation {
|
||||
pos: Vec3,
|
||||
rot: Vec3,
|
||||
}
|
||||
|
||||
fn rot_quat(rot: Vec3) -> Quat {
|
||||
Quat::from_euler(glam::EulerRot::XYZ, rot.z, rot.y, rot.x)
|
||||
}
|
||||
|
||||
impl CameraLocation {
|
||||
pub fn new() -> CameraLocation {
|
||||
let rot = vec3(std::f32::consts::FRAC_PI_4, 0., 0.);
|
||||
let pos = rot_quat(rot) * vec3(-200., 0., 50.);
|
||||
CameraLocation { pos, rot }
|
||||
}
|
||||
|
||||
pub fn view_mtx(&self) -> Mat4 {
|
||||
Mat4::from_quat(rot_quat(-self.rot)) * Mat4::from_translation(-self.pos)
|
||||
}
|
||||
|
||||
pub fn move_rel(&mut self, offset: Vec3) {
|
||||
self.pos += rot_quat(vec3(self.rot.x, 0., 0.)) * offset;
|
||||
}
|
||||
|
||||
pub fn rotate_rel_ypr(&mut self, ypr: Vec3) {
|
||||
self.rot += ypr;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
mod keyctl {
|
||||
use std::{collections::HashSet, iter::Sum};
|
||||
use winit::{event::ElementState, keyboard::PhysicalKey};
|
||||
|
||||
pub struct Keyboard {
|
||||
pressed: HashSet<PhysicalKey>,
|
||||
}
|
||||
|
||||
impl Keyboard {
|
||||
pub fn new() -> Self {
|
||||
Keyboard {
|
||||
pressed: Default::default(),
|
||||
}
|
||||
}
|
||||
|
||||
pub fn is_pressed(&self, key: PhysicalKey) -> bool {
|
||||
self.pressed.contains(&key)
|
||||
}
|
||||
|
||||
pub fn set_key_state(&mut self, key: PhysicalKey, state: ElementState) {
|
||||
match state {
|
||||
ElementState::Pressed => self.pressed.insert(key),
|
||||
ElementState::Released => self.pressed.remove(&key),
|
||||
};
|
||||
}
|
||||
|
||||
pub fn control<T: Copy + Sum>(&self, keymap: &[(PhysicalKey, T)]) -> T {
|
||||
keymap
|
||||
.iter()
|
||||
.copied()
|
||||
.filter_map(|(key, ctl)| self.is_pressed(key).then_some(ctl))
|
||||
.sum()
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
static KEYS_MOVE: &'static [(PhysicalKey, Vec3)] = &[
|
||||
(PhysicalKey::Code(KeyCode::KeyW), vec3(1., 0., 0.)),
|
||||
(PhysicalKey::Code(KeyCode::KeyS), vec3(-1., 0., 0.)),
|
||||
(PhysicalKey::Code(KeyCode::KeyA), vec3(0., 1., 0.)),
|
||||
(PhysicalKey::Code(KeyCode::KeyD), vec3(0., -1., 0.)),
|
||||
(PhysicalKey::Code(KeyCode::KeyE), vec3(0., 0., 1.)),
|
||||
(PhysicalKey::Code(KeyCode::KeyQ), vec3(0., 0., -1.)),
|
||||
(PhysicalKey::Code(KeyCode::Space), vec3(0., 0., 1.)),
|
||||
(PhysicalKey::Code(KeyCode::ShiftLeft), vec3(0., 0., -1.)),
|
||||
];
|
||||
|
||||
static KEYS_ROTATE: &'static [(PhysicalKey, Vec3)] = &[
|
||||
(PhysicalKey::Code(KeyCode::Numpad4), vec3(1., 0., 0.)),
|
||||
(PhysicalKey::Code(KeyCode::Numpad6), vec3(-1., 0., 0.)),
|
||||
(PhysicalKey::Code(KeyCode::Numpad5), vec3(0., 1., 0.)),
|
||||
(PhysicalKey::Code(KeyCode::Numpad8), vec3(0., -1., 0.)),
|
||||
(PhysicalKey::Code(KeyCode::Numpad9), vec3(0., 0., 1.)),
|
||||
(PhysicalKey::Code(KeyCode::Numpad7), vec3(0., 0., -1.)),
|
||||
];
|
||||
|
||||
struct State<'a> {
|
||||
device: wgpu::Device,
|
||||
queue: wgpu::Queue,
|
||||
|
||||
fps: fps::Counter,
|
||||
kbd: keyctl::Keyboard,
|
||||
t1: Instant,
|
||||
|
||||
viewport: viewport::Viewport<'a>,
|
||||
cam_loc: camctl::CameraLocation,
|
||||
cam_obj: camera::Camera,
|
||||
line_rend: lines::LineRenderer,
|
||||
mesh_rend: meshes::Renderer,
|
||||
|
||||
lines: Vec<lines::Line>,
|
||||
meshes: Vec<meshes::Mesh>,
|
||||
|
||||
window: &'a Window,
|
||||
}
|
||||
|
||||
impl<'a> State<'a> {
|
||||
async fn new(window: &'a Window) -> State<'a> {
|
||||
let size = window.inner_size();
|
||||
|
||||
let instance = wgpu::Instance::new(wgpu::InstanceDescriptor {
|
||||
backends: wgpu::Backends::PRIMARY,
|
||||
..Default::default()
|
||||
});
|
||||
let surface = instance.create_surface(window).unwrap();
|
||||
let adapter = instance
|
||||
.request_adapter(&wgpu::RequestAdapterOptions {
|
||||
power_preference: wgpu::PowerPreference::default(),
|
||||
compatible_surface: Some(&surface),
|
||||
force_fallback_adapter: false,
|
||||
})
|
||||
.await
|
||||
.unwrap();
|
||||
let (device, queue) = adapter
|
||||
.request_device(
|
||||
&wgpu::DeviceDescriptor {
|
||||
label: None,
|
||||
required_features: wgpu::Features::PUSH_CONSTANTS
|
||||
| wgpu::Features::TEXTURE_ADAPTER_SPECIFIC_FORMAT_FEATURES,
|
||||
required_limits: wgpu::Limits {
|
||||
max_push_constant_size: 16,
|
||||
..wgpu::Limits::default()
|
||||
},
|
||||
memory_hints: Default::default(),
|
||||
},
|
||||
None, // Trace path
|
||||
)
|
||||
.await
|
||||
.unwrap();
|
||||
|
||||
let viewport = viewport::Viewport::new(&adapter, &device, surface, uvec2(size.width, size.height));
|
||||
|
||||
let kbd = keyctl::Keyboard::new();
|
||||
let cam_loc = camctl::CameraLocation::new();
|
||||
let t1 = Instant::now();
|
||||
|
||||
let cam_obj = camera::Camera::new(&device);
|
||||
let line_rend = lines::LineRenderer::new(
|
||||
&device,
|
||||
cam_obj.bind_group_layout(),
|
||||
viewport.format(),
|
||||
Some(wgpu::DepthStencilState {
|
||||
format: wgpu::TextureFormat::Depth24Plus,
|
||||
depth_write_enabled: false,
|
||||
depth_compare: wgpu::CompareFunction::LessEqual,
|
||||
stencil: wgpu::StencilState::default(),
|
||||
bias: wgpu::DepthBiasState::default(),
|
||||
}),
|
||||
wgpu::MultisampleState {
|
||||
count: viewport.sample_count(),
|
||||
mask: !0,
|
||||
alpha_to_coverage_enabled: false,
|
||||
},
|
||||
);
|
||||
let mesh_rend = meshes::Renderer::new(
|
||||
&device,
|
||||
cam_obj.bind_group_layout(),
|
||||
viewport.format(),
|
||||
Some(wgpu::DepthStencilState {
|
||||
format: wgpu::TextureFormat::Depth24Plus,
|
||||
depth_write_enabled: true,
|
||||
depth_compare: wgpu::CompareFunction::LessEqual,
|
||||
stencil: wgpu::StencilState::default(),
|
||||
bias: wgpu::DepthBiasState::default(),
|
||||
}),
|
||||
wgpu::MultisampleState {
|
||||
count: viewport.sample_count(),
|
||||
mask: !0,
|
||||
alpha_to_coverage_enabled: false,
|
||||
},
|
||||
);
|
||||
|
||||
let (meshes, lines) = prepare_scene(&device);
|
||||
|
||||
let fps = fps::Counter::new();
|
||||
|
||||
Self {
|
||||
device,
|
||||
queue,
|
||||
viewport,
|
||||
line_rend,
|
||||
mesh_rend,
|
||||
kbd,
|
||||
fps,
|
||||
cam_loc,
|
||||
cam_obj,
|
||||
t1,
|
||||
lines,
|
||||
meshes,
|
||||
window,
|
||||
}
|
||||
}
|
||||
|
||||
pub fn window(&self) -> &Window {
|
||||
&self.window
|
||||
}
|
||||
|
||||
fn resize(&mut self, new_size: winit::dpi::PhysicalSize<u32>) {
|
||||
if new_size.width > 0 && new_size.height > 0 {
|
||||
self.viewport
|
||||
.resize(&self.device, uvec2(new_size.width, new_size.height));
|
||||
}
|
||||
}
|
||||
|
||||
fn update(&mut self) {
|
||||
let dt = {
|
||||
let t2 = Instant::now();
|
||||
let dt = t2 - self.t1;
|
||||
self.t1 = t2;
|
||||
dt.as_secs_f32()
|
||||
};
|
||||
let size = self.viewport.size().as_vec2();
|
||||
|
||||
self.cam_loc.move_rel(100. * dt * self.kbd.control(&KEYS_MOVE));
|
||||
self.cam_loc.rotate_rel_ypr(2. * dt * self.kbd.control(&KEYS_ROTATE));
|
||||
self.cam_obj.set(&self.queue, self.cam_loc.view_mtx(), size);
|
||||
}
|
||||
|
||||
fn render(&mut self) -> Result<(), wgpu::SurfaceError> {
|
||||
self.fps.on_frame();
|
||||
self.window
|
||||
.set_title(&format!("Space Refraction ({:.1} FPS)", self.fps.get()));
|
||||
self.viewport
|
||||
.render_single_pass(&self.device, &self.queue, |mut render_pass| {
|
||||
self.mesh_rend
|
||||
.render(&mut render_pass, self.cam_obj.bind_group(), self.meshes.iter());
|
||||
self.line_rend
|
||||
.render(&mut render_pass, self.cam_obj.bind_group(), self.lines.iter());
|
||||
})
|
||||
}
|
||||
}
|
||||
|
||||
pub async fn run() {
|
||||
let event_loop = EventLoop::new().unwrap();
|
||||
let window = WindowBuilder::new()
|
||||
.with_title("Refraction: Wireframe")
|
||||
.build(&event_loop)
|
||||
.unwrap();
|
||||
|
||||
// State::new uses async code, so we're going to wait for it to finish
|
||||
let mut state = State::new(&window).await;
|
||||
let mut surface_configured = false;
|
||||
|
||||
event_loop
|
||||
.run(move |event, control_flow| {
|
||||
match event {
|
||||
Event::WindowEvent { ref event, window_id } if window_id == state.window().id() => {
|
||||
match event {
|
||||
WindowEvent::KeyboardInput {
|
||||
device_id: _,
|
||||
event,
|
||||
is_synthetic: _,
|
||||
} => {
|
||||
state.kbd.set_key_state(event.physical_key, event.state);
|
||||
}
|
||||
WindowEvent::CloseRequested => control_flow.exit(),
|
||||
WindowEvent::Resized(physical_size) => {
|
||||
surface_configured = true;
|
||||
state.resize(*physical_size);
|
||||
}
|
||||
WindowEvent::RedrawRequested => {
|
||||
// This tells winit that we want another frame after this one
|
||||
state.window().request_redraw();
|
||||
|
||||
if !surface_configured {
|
||||
return;
|
||||
}
|
||||
|
||||
state.update();
|
||||
match state.render() {
|
||||
Ok(_) => {}
|
||||
// Reconfigure the surface if it's lost or outdated
|
||||
Err(wgpu::SurfaceError::Lost | wgpu::SurfaceError::Outdated) => {
|
||||
state.viewport.configure(&state.device);
|
||||
}
|
||||
// The system is out of memory, we should probably quit
|
||||
Err(wgpu::SurfaceError::OutOfMemory) => {
|
||||
eprintln!("OutOfMemory");
|
||||
control_flow.exit();
|
||||
}
|
||||
|
||||
// This happens when the a frame takes too long to present
|
||||
Err(wgpu::SurfaceError::Timeout) => {
|
||||
eprintln!("Surface timeout")
|
||||
}
|
||||
}
|
||||
}
|
||||
_ => {}
|
||||
}
|
||||
}
|
||||
_ => {}
|
||||
}
|
||||
})
|
||||
.unwrap();
|
||||
}
|
||||
|
||||
fn main() {
|
||||
pollster::block_on(run());
|
||||
}
|
||||
|
||||
mod fps {
|
||||
use std::time::{Duration, Instant};
|
||||
|
||||
pub struct Counter {
|
||||
fps: f32,
|
||||
t1: Instant,
|
||||
frames: u32,
|
||||
}
|
||||
|
||||
impl Counter {
|
||||
pub fn new() -> Self {
|
||||
Self {
|
||||
fps: 0.,
|
||||
t1: Instant::now(),
|
||||
frames: 0,
|
||||
}
|
||||
}
|
||||
|
||||
pub fn get(&self) -> f32 {
|
||||
self.fps
|
||||
}
|
||||
|
||||
pub fn on_frame(&mut self) {
|
||||
self.frames += 1;
|
||||
let t2 = Instant::now();
|
||||
let dt = t2 - self.t1;
|
||||
if dt >= Duration::from_secs(1) {
|
||||
*self = Self {
|
||||
fps: self.frames as f32 / dt.as_secs_f32(),
|
||||
t1: t2,
|
||||
frames: 0,
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -1,63 +0,0 @@
|
|||
struct CameraUniform {
|
||||
mvp: mat4x4<f32>,
|
||||
scale: vec2<f32>,
|
||||
}
|
||||
@group(0) @binding(0)
|
||||
var<uniform> camera: CameraUniform;
|
||||
|
||||
struct MeshUniform {
|
||||
color: vec4<f32>,
|
||||
}
|
||||
var<push_constant> mesh: MeshUniform;
|
||||
|
||||
struct VertexInput {
|
||||
@location(0) position: vec3<f32>,
|
||||
@location(1) normal: vec3<f32>,
|
||||
}
|
||||
|
||||
struct VertexOutput {
|
||||
@builtin(position) clip_position: vec4<f32>,
|
||||
@location(0) vertex_color: vec4<f32>,
|
||||
}
|
||||
|
||||
struct FragmentOutput {
|
||||
@location(0) color: vec4<f32>,
|
||||
@builtin(sample_mask) mask: u32,
|
||||
}
|
||||
|
||||
fn hash(key : u32) -> u32 {
|
||||
var v = key;
|
||||
v *= 0xb384af1bu;
|
||||
v ^= v >> 15u;
|
||||
return v;
|
||||
}
|
||||
|
||||
@vertex
|
||||
fn vs_main(ver: VertexInput) -> VertexOutput {
|
||||
var out: VertexOutput;
|
||||
var light = dot(ver.normal, normalize(vec3(1., 1., 1.)));
|
||||
light = .7 + .3 * light;
|
||||
out.vertex_color = vec4(light * mesh.color.xyz, mesh.color.w);
|
||||
out.clip_position = camera.mvp * vec4(ver.position, 1.);
|
||||
return out;
|
||||
}
|
||||
|
||||
@fragment
|
||||
fn fs_main(in: VertexOutput) -> FragmentOutput {
|
||||
var out: FragmentOutput;
|
||||
out.color = vec4(in.vertex_color.xyz, 1.);
|
||||
var x = bitcast<u32>(in.clip_position.x);
|
||||
var y = bitcast<u32>(in.clip_position.y);
|
||||
var z = bitcast<u32>(in.clip_position.z);
|
||||
var alpha = in.vertex_color.w;
|
||||
var seed = hash(hash(hash(x) ^ y) ^ z);
|
||||
var mask = 0u;
|
||||
for (var sample = 0u; sample < 8u; sample++) {
|
||||
var threshold = f32(hash(seed ^ sample)) / 0x1p32;
|
||||
if (alpha > threshold) {
|
||||
mask |= 1u << sample;
|
||||
}
|
||||
}
|
||||
out.mask = mask;
|
||||
return out;
|
||||
}
|
||||
|
|
@ -1,160 +0,0 @@
|
|||
use std::mem;
|
||||
|
||||
use bytemuck::{bytes_of, cast_slice, Pod, Zeroable};
|
||||
use glam::{Vec3, Vec4};
|
||||
use wgpu::util::DeviceExt as _;
|
||||
|
||||
use crate::scene::Face;
|
||||
|
||||
#[repr(C)]
|
||||
#[derive(Copy, Clone, Debug, Pod, Zeroable)]
|
||||
struct Vertex {
|
||||
pub position: [f32; 3],
|
||||
pub normal: [f32; 3],
|
||||
}
|
||||
|
||||
impl From<crate::scene::Vertex> for Vertex {
|
||||
fn from(value: crate::scene::Vertex) -> Self {
|
||||
Vertex {
|
||||
position: value.position.into(),
|
||||
normal: value.normal.into(),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[repr(C)]
|
||||
#[derive(Copy, Clone, Pod, Zeroable)]
|
||||
struct PushConsts {
|
||||
pub color: [f32; 4],
|
||||
}
|
||||
|
||||
#[derive(Copy, Clone)]
|
||||
pub struct Attrs {
|
||||
pub color: Vec4,
|
||||
}
|
||||
|
||||
impl Attrs {
|
||||
fn consts(&self) -> PushConsts {
|
||||
PushConsts {
|
||||
color: self.color.to_array(),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
pub struct Mesh {
|
||||
consts: PushConsts,
|
||||
npoints: u32,
|
||||
buf: wgpu::Buffer,
|
||||
}
|
||||
|
||||
impl Mesh {
|
||||
pub fn new_list(device: &wgpu::Device, attrs: Attrs, tris: Vec<Face>) -> Self {
|
||||
let data: Vec<Vertex> = tris.into_iter().flat_map(|face| face.map(Into::into)).collect();
|
||||
let buf = device.create_buffer_init(&wgpu::util::BufferInitDescriptor {
|
||||
label: Some("Mesh Vertex Buffer"),
|
||||
contents: cast_slice(&data),
|
||||
usage: wgpu::BufferUsages::VERTEX,
|
||||
});
|
||||
Mesh {
|
||||
consts: attrs.consts(),
|
||||
npoints: data.len() as u32,
|
||||
buf,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
pub struct Renderer {
|
||||
pipeline: wgpu::RenderPipeline,
|
||||
}
|
||||
|
||||
static SHADER: &'static str = include_str!("mesh.wgsl");
|
||||
|
||||
impl Renderer {
|
||||
pub fn new(
|
||||
device: &wgpu::Device,
|
||||
cam_layout: &wgpu::BindGroupLayout,
|
||||
target_format: wgpu::TextureFormat,
|
||||
depth_stencil: Option<wgpu::DepthStencilState>,
|
||||
multisample: wgpu::MultisampleState,
|
||||
) -> Renderer {
|
||||
let shader = device.create_shader_module(wgpu::ShaderModuleDescriptor {
|
||||
label: Some("Mesh Shader"),
|
||||
source: wgpu::ShaderSource::Wgsl(SHADER.into()),
|
||||
});
|
||||
|
||||
let consts_range = wgpu::PushConstantRange {
|
||||
stages: wgpu::ShaderStages::VERTEX,
|
||||
range: 0..mem::size_of::<PushConsts>() as u32,
|
||||
};
|
||||
|
||||
let layout = device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor {
|
||||
label: Some("Mesh RenderPipelineLayout"),
|
||||
bind_group_layouts: &[cam_layout],
|
||||
push_constant_ranges: &[consts_range],
|
||||
});
|
||||
|
||||
let pipeline = device.create_render_pipeline(&wgpu::RenderPipelineDescriptor {
|
||||
label: Some("Mesh RenderPipeline"),
|
||||
layout: Some(&layout),
|
||||
vertex: wgpu::VertexState {
|
||||
module: &shader,
|
||||
entry_point: "vs_main",
|
||||
buffers: &[wgpu::VertexBufferLayout {
|
||||
array_stride: mem::size_of::<Vertex>() as wgpu::BufferAddress,
|
||||
step_mode: wgpu::VertexStepMode::Vertex,
|
||||
attributes: &[
|
||||
wgpu::VertexAttribute {
|
||||
offset: mem::offset_of!(Vertex, position) as u64,
|
||||
shader_location: 0,
|
||||
format: wgpu::VertexFormat::Float32x3,
|
||||
},
|
||||
wgpu::VertexAttribute {
|
||||
offset: mem::offset_of!(Vertex, normal) as u64,
|
||||
shader_location: 1,
|
||||
format: wgpu::VertexFormat::Float32x3,
|
||||
},
|
||||
],
|
||||
}],
|
||||
compilation_options: Default::default(),
|
||||
},
|
||||
fragment: Some(wgpu::FragmentState {
|
||||
module: &shader,
|
||||
entry_point: "fs_main",
|
||||
targets: &[Some(wgpu::ColorTargetState {
|
||||
format: target_format,
|
||||
blend: Some(wgpu::BlendState {
|
||||
color: wgpu::BlendComponent::OVER,
|
||||
alpha: wgpu::BlendComponent::OVER,
|
||||
}),
|
||||
write_mask: wgpu::ColorWrites::ALL,
|
||||
})],
|
||||
compilation_options: Default::default(),
|
||||
}),
|
||||
primitive: wgpu::PrimitiveState {
|
||||
topology: wgpu::PrimitiveTopology::TriangleList,
|
||||
..Default::default()
|
||||
},
|
||||
depth_stencil,
|
||||
multisample,
|
||||
multiview: None,
|
||||
cache: None,
|
||||
});
|
||||
|
||||
Renderer { pipeline }
|
||||
}
|
||||
|
||||
pub fn render<'a>(
|
||||
&self,
|
||||
pass: &mut wgpu::RenderPass,
|
||||
cam_bind: &wgpu::BindGroup,
|
||||
meshes: impl Iterator<Item = &'a Mesh>,
|
||||
) {
|
||||
pass.set_pipeline(&self.pipeline);
|
||||
pass.set_bind_group(0, cam_bind, &[]);
|
||||
for mesh in meshes {
|
||||
pass.set_push_constants(wgpu::ShaderStages::VERTEX, 0, bytes_of(&mesh.consts));
|
||||
pass.set_vertex_buffer(0, mesh.buf.slice(..));
|
||||
pass.draw(0..mesh.npoints, 0..1);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -1,61 +0,0 @@
|
|||
struct CameraUniform {
|
||||
mvp: mat4x4<f32>,
|
||||
scale: vec2<f32>,
|
||||
}
|
||||
@group(0) @binding(0)
|
||||
var<uniform> camera: CameraUniform;
|
||||
|
||||
struct LineUniform {
|
||||
color: vec3<f32>,
|
||||
}
|
||||
const width = 0.1;
|
||||
var<push_constant> line: LineUniform;
|
||||
|
||||
struct SegmentInput {
|
||||
@location(0) a: vec3<f32>,
|
||||
@location(1) ad: vec3<f32>,
|
||||
@location(2) b: vec3<f32>,
|
||||
@location(3) bd: vec3<f32>,
|
||||
}
|
||||
struct OffsetInput {
|
||||
@builtin(vertex_index) idx: u32,
|
||||
}
|
||||
|
||||
struct VertexOutput {
|
||||
@builtin(position) clip_position: vec4<f32>,
|
||||
@location(0) vertex_color: vec3<f32>,
|
||||
}
|
||||
|
||||
@vertex
|
||||
fn vs_main(seg: SegmentInput, off: OffsetInput) -> VertexOutput {
|
||||
var out: VertexOutput;
|
||||
out.vertex_color = line.color;
|
||||
var pt: vec3<f32>;
|
||||
var dir: vec3<f32>;
|
||||
switch (off.idx) {
|
||||
case 0u: { pt = seg.a; dir = seg.ad; }
|
||||
case 1u: { pt = seg.a; dir = seg.ad; }
|
||||
case 2u: { pt = seg.b; dir = seg.bd; }
|
||||
case 3u: { pt = seg.b; dir = seg.bd; }
|
||||
default: {}
|
||||
}
|
||||
var sgn: f32;
|
||||
switch (off.idx) {
|
||||
case 0u: { sgn = -1.; }
|
||||
case 1u: { sgn = 1.; }
|
||||
case 2u: { sgn = -1.; }
|
||||
case 3u: { sgn = 1.; }
|
||||
default: {}
|
||||
}
|
||||
let pt_cs = camera.mvp * vec4(pt, 1.);
|
||||
let dir_cs0 = camera.mvp * vec4(dir, 0.);
|
||||
let dir_cs = dir_cs0.xyz * pt_cs.w - pt_cs.xyz * dir_cs0.w;
|
||||
let normal_cs = camera.scale * normalize(vec2(-dir_cs.y, dir_cs.x));
|
||||
out.clip_position = pt_cs + vec4(sgn * width * normal_cs, 0., 0.);
|
||||
return out;
|
||||
}
|
||||
|
||||
@fragment
|
||||
fn fs_main(in: VertexOutput) -> @location(0) vec4<f32> {
|
||||
return 0.5 * vec4(in.vertex_color, 1.0);
|
||||
}
|
||||
|
|
@ -1,210 +0,0 @@
|
|||
use glam::*;
|
||||
use itertools::chain;
|
||||
|
||||
use refraction::ifaces::{DebugTraceable, Traceable};
|
||||
use refraction::tube::metric::Tube;
|
||||
use refraction::tube::Space;
|
||||
use refraction::types::{Location, Object, Ray};
|
||||
use refraction::utils::put_object;
|
||||
|
||||
pub enum Line {
|
||||
Strip(Vec<Ray>),
|
||||
Loop(Vec<Ray>),
|
||||
}
|
||||
|
||||
pub struct FancyLine {
|
||||
pub color: Vec3,
|
||||
pub pts: Vec<Ray>,
|
||||
}
|
||||
|
||||
pub struct Vertex {
|
||||
pub position: Vec3,
|
||||
pub normal: Vec3,
|
||||
}
|
||||
|
||||
pub type Face = [Vertex; 3];
|
||||
|
||||
pub enum Mesh {
|
||||
List(Vec<Face>),
|
||||
}
|
||||
|
||||
pub struct FancyMesh {
|
||||
pub color: Vec4,
|
||||
pub tris: Vec<Face>,
|
||||
}
|
||||
|
||||
fn paint(onto: &mut Vec<FancyLine>, color: Vec3, lines: Vec<Line>) {
|
||||
onto.extend(lines.into_iter().map(move |line| FancyLine {
|
||||
color,
|
||||
pts: match line {
|
||||
Line::Strip(pts) => pts,
|
||||
Line::Loop(mut pts) => {
|
||||
pts.push(*pts.first().unwrap());
|
||||
pts
|
||||
}
|
||||
},
|
||||
}))
|
||||
}
|
||||
|
||||
fn draw_line(a: Vec3, b: Vec3) -> Line {
|
||||
let dir = (b - a).normalize();
|
||||
Line::Strip(vec![Ray { pos: a, dir }, Ray { pos: b, dir }])
|
||||
}
|
||||
|
||||
fn draw_mark(pos: Vec3) -> Vec<Line> {
|
||||
[
|
||||
vec3(1., 1., 1.),
|
||||
vec3(1., 1., -1.),
|
||||
vec3(1., -1., 1.),
|
||||
vec3(1., -1., -1.),
|
||||
]
|
||||
.into_iter()
|
||||
.map(|off| draw_line(pos - off, pos + off))
|
||||
.collect()
|
||||
}
|
||||
|
||||
pub fn build() -> (Vec<FancyMesh>, Vec<FancyLine>) {
|
||||
let tube = Tube {
|
||||
inner_radius: 30.0,
|
||||
outer_radius: 50.0,
|
||||
internal_halflength: 100.0,
|
||||
external_halflength: 300.0,
|
||||
};
|
||||
|
||||
let objs: Vec<_> = [-1.25, -1.00, -0.85, -0.50, 0.00, 0.40, 0.70, 0.95, 1.05]
|
||||
.iter()
|
||||
.enumerate()
|
||||
.map(|(k, &y)| Object {
|
||||
id: k as i32,
|
||||
loc: put_object(
|
||||
&tube,
|
||||
vec3(0.0, y * tube.external_halflength, 0.0),
|
||||
Mat3::from_mat2(Mat2::from_angle(y)),
|
||||
),
|
||||
r: 20.0,
|
||||
})
|
||||
.collect();
|
||||
let space = Space { tube, objs };
|
||||
let cam1 = put_object(&space.tube, vec3(-500., 0., 0.), Mat3::IDENTITY);
|
||||
let cam2 = put_object(
|
||||
&space.tube,
|
||||
vec3(-2.5 * tube.outer_radius, 1.25 * tube.external_halflength, 0.),
|
||||
mat3(vec3(1., -1., 0.), vec3(1., 1., 0.), vec3(0., 0., 1.)),
|
||||
);
|
||||
let cam2l = put_object(
|
||||
&space.tube,
|
||||
vec3(-2.5 * tube.outer_radius, 1.25 * tube.external_halflength, 0.),
|
||||
mat3(vec3(1., -0.825, 0.), vec3(1., 1., 0.), vec3(0., 0., 1.)),
|
||||
);
|
||||
let cam3 = put_object(
|
||||
&space.tube,
|
||||
vec3(0.25 * tube.inner_radius, 0.25 * tube.external_halflength, 0.),
|
||||
mat3(vec3(0., -1., 0.), vec3(1., 0., 0.), vec3(0., 0., 1.)),
|
||||
);
|
||||
|
||||
let mut gc = vec![];
|
||||
gc.push(FancyMesh {
|
||||
color: vec4(0.10, 0.12, 0.15, 0.8),
|
||||
tris: tube.render(),
|
||||
});
|
||||
let meshes = gc;
|
||||
|
||||
let mut gc = vec![];
|
||||
paint(&mut gc, vec3(0.0, 0.6, 1.0), draw_fan_2(&space, cam3, vec3(0., 1., 0.)));
|
||||
paint(&mut gc, vec3(0.2, 1.0, 0.0), draw_fan_2(&space, cam2, vec3(0., 1., 0.)));
|
||||
paint(
|
||||
&mut gc,
|
||||
vec3(0.0, 1.0, 0.6),
|
||||
draw_fan_2(&space, cam2l, vec3(0., 0., 1.)),
|
||||
);
|
||||
paint(&mut gc, vec3(1.0, 0.2, 0.0), draw_fan_2(&space, cam1, vec3(0., 1., 0.)));
|
||||
let lines = gc;
|
||||
|
||||
(meshes, lines)
|
||||
}
|
||||
|
||||
fn draw_ray_2(gc: &mut Vec<Line>, space: &Space, camera: Location, dir: Vec3) {
|
||||
let pos = vec3(0., 0., 0.);
|
||||
let (hits, path) = space.trace_dbg(camera, Ray { pos, dir });
|
||||
if true {
|
||||
let hits2 = space.trace(camera, Ray { pos, dir });
|
||||
for (a, b) in hits.iter().zip(hits2.into_iter()) {
|
||||
assert_eq!(a.id, b.id);
|
||||
assert_eq!(a.pos, b.pos);
|
||||
assert_eq!(a.rel, b.rel);
|
||||
}
|
||||
}
|
||||
|
||||
for hit in hits {
|
||||
gc.extend(draw_mark(hit.pos));
|
||||
}
|
||||
|
||||
let mut pts = path.points;
|
||||
let end_pos = *pts.last().expect("the starting point is always in the path");
|
||||
pts.extend(itertools::iterate(end_pos, |r| r.forward(100.0)).take(1000));
|
||||
gc.push(Line::Strip(pts));
|
||||
}
|
||||
|
||||
fn draw_fan_2(space: &Space, camera: Location, spread: Vec3) -> Vec<Line> {
|
||||
let mut gc = vec![];
|
||||
for δ in itertools_num::linspace(-1., 1., 101) {
|
||||
draw_ray_2(&mut gc, space, camera, vec3(1., 0., 0.) + δ * spread);
|
||||
}
|
||||
gc
|
||||
}
|
||||
|
||||
trait Renderable {
|
||||
fn render(&self) -> Vec<Face>;
|
||||
}
|
||||
|
||||
impl Renderable for Tube {
|
||||
fn render(&self) -> Vec<Face> {
|
||||
let sides = 42;
|
||||
let step = 2. * std::f32::consts::PI / sides as f32;
|
||||
let dir = |k| {
|
||||
let d = Vec2::from_angle(k as f32 * step);
|
||||
vec3(d.x, 0., d.y)
|
||||
};
|
||||
|
||||
let side = vec3(0., self.external_halflength, 0.);
|
||||
let r1 = self.inner_radius;
|
||||
let r2 = self.outer_radius;
|
||||
|
||||
let vertex = |k, r, h| {
|
||||
let n = dir(k);
|
||||
Vertex {
|
||||
position: r * n + h * side,
|
||||
normal: n,
|
||||
}
|
||||
};
|
||||
let inner = (0..sides).flat_map(|k| {
|
||||
[
|
||||
[vertex(k, -r1, -1.), vertex(k, -r1, 1.), vertex(k + 1, -r1, -1.)],
|
||||
[vertex(k + 1, -r1, -1.), vertex(k, -r1, 1.), vertex(k + 1, -r1, 1.)],
|
||||
]
|
||||
});
|
||||
let outer = (0..sides).flat_map(|k| {
|
||||
[
|
||||
[vertex(k, r2, -1.), vertex(k + 1, r2, -1.), vertex(k, r2, 1.)],
|
||||
[vertex(k + 1, r2, -1.), vertex(k + 1, r2, 1.), vertex(k, r2, 1.)],
|
||||
]
|
||||
});
|
||||
|
||||
let vertex = |k, r, h| {
|
||||
let n = dir(k);
|
||||
Vertex {
|
||||
position: r * n + h * side,
|
||||
normal: h * Vec3::Y,
|
||||
}
|
||||
};
|
||||
let cap = (0..sides).flat_map(|k| {
|
||||
[
|
||||
[vertex(k, r1, -1.), vertex(k + 1, r1, -1.), vertex(k + 1, r2, -1.)],
|
||||
[vertex(k, r1, -1.), vertex(k + 1, r2, -1.), vertex(k, r2, -1.)],
|
||||
[vertex(k, r1, 1.), vertex(k + 1, r1, 1.), vertex(k + 1, r2, 1.)],
|
||||
[vertex(k, r1, 1.), vertex(k + 1, r2, 1.), vertex(k, r2, 1.)],
|
||||
]
|
||||
});
|
||||
chain!(inner, outer, cap).collect()
|
||||
}
|
||||
}
|
||||
|
|
@ -1,147 +0,0 @@
|
|||
use glam::{uvec2, UVec2};
|
||||
|
||||
pub struct Viewport<'a> {
|
||||
surface: wgpu::Surface<'a>,
|
||||
config: wgpu::SurfaceConfiguration,
|
||||
sample_count: u32,
|
||||
multisample: Multisample,
|
||||
}
|
||||
|
||||
impl<'a> Viewport<'a> {
|
||||
pub fn new(adapter: &wgpu::Adapter, device: &wgpu::Device, surface: wgpu::Surface<'a>, size: UVec2) -> Self {
|
||||
let caps = surface.get_capabilities(adapter);
|
||||
let format = wgpu::TextureFormat::Bgra8UnormSrgb;
|
||||
let sample_count = adapter
|
||||
.get_texture_format_features(format)
|
||||
.flags
|
||||
.supported_sample_counts()
|
||||
.into_iter()
|
||||
.max()
|
||||
.unwrap();
|
||||
eprintln!("Using x{sample_count} mutlisampling");
|
||||
let config = wgpu::SurfaceConfiguration {
|
||||
usage: wgpu::TextureUsages::RENDER_ATTACHMENT,
|
||||
format,
|
||||
width: size.x,
|
||||
height: size.y,
|
||||
present_mode: wgpu::PresentMode::Fifo,
|
||||
alpha_mode: caps.alpha_modes[0],
|
||||
view_formats: vec![],
|
||||
desired_maximum_frame_latency: 2,
|
||||
};
|
||||
let multisample = Multisample::new(device, format, size, sample_count);
|
||||
Self {
|
||||
surface,
|
||||
config,
|
||||
sample_count,
|
||||
multisample,
|
||||
}
|
||||
}
|
||||
|
||||
pub fn configure(&mut self, device: &wgpu::Device) {
|
||||
self.surface.configure(&device, &self.config);
|
||||
self.multisample = Multisample::new(device, self.format(), self.size(), self.sample_count);
|
||||
}
|
||||
|
||||
pub fn resize(&mut self, device: &wgpu::Device, size: UVec2) {
|
||||
self.config.width = size.x;
|
||||
self.config.height = size.y;
|
||||
self.configure(&device);
|
||||
}
|
||||
|
||||
pub fn size(&self) -> UVec2 {
|
||||
uvec2(self.config.width, self.config.height)
|
||||
}
|
||||
|
||||
pub fn format(&self) -> wgpu::TextureFormat {
|
||||
self.config.format
|
||||
}
|
||||
|
||||
pub fn sample_count(&self) -> u32 {
|
||||
self.sample_count
|
||||
}
|
||||
|
||||
pub fn render_single_pass(
|
||||
&mut self,
|
||||
device: &wgpu::Device,
|
||||
queue: &wgpu::Queue,
|
||||
f: impl FnOnce(wgpu::RenderPass),
|
||||
) -> Result<(), wgpu::SurfaceError> {
|
||||
let output = self.surface.get_current_texture()?;
|
||||
let view = output.texture.create_view(&wgpu::TextureViewDescriptor::default());
|
||||
let mut encoder = device.create_command_encoder(&wgpu::CommandEncoderDescriptor {
|
||||
label: Some("Render CommandEncoder"),
|
||||
});
|
||||
let render_pass = encoder.begin_render_pass(&wgpu::RenderPassDescriptor {
|
||||
label: Some("RenderPass"),
|
||||
color_attachments: &[Some(wgpu::RenderPassColorAttachment {
|
||||
view: &self.multisample.color,
|
||||
resolve_target: Some(&view),
|
||||
ops: wgpu::Operations {
|
||||
load: wgpu::LoadOp::Clear(wgpu::Color {
|
||||
r: 0.,
|
||||
g: 0.,
|
||||
b: 0.,
|
||||
a: 1.,
|
||||
}),
|
||||
store: wgpu::StoreOp::Store,
|
||||
},
|
||||
})],
|
||||
depth_stencil_attachment: Some(wgpu::RenderPassDepthStencilAttachment {
|
||||
view: &self.multisample.depth,
|
||||
depth_ops: Some(wgpu::Operations {
|
||||
load: wgpu::LoadOp::Clear(1.),
|
||||
store: wgpu::StoreOp::Discard,
|
||||
}),
|
||||
stencil_ops: None,
|
||||
}),
|
||||
occlusion_query_set: None,
|
||||
timestamp_writes: None,
|
||||
});
|
||||
f(render_pass);
|
||||
queue.submit(std::iter::once(encoder.finish()));
|
||||
output.present();
|
||||
Ok(())
|
||||
}
|
||||
}
|
||||
|
||||
struct Multisample {
|
||||
color: wgpu::TextureView,
|
||||
depth: wgpu::TextureView,
|
||||
}
|
||||
|
||||
impl Multisample {
|
||||
fn new(device: &wgpu::Device, format: wgpu::TextureFormat, size: UVec2, sample_count: u32) -> Multisample {
|
||||
let color = device.create_texture(&wgpu::TextureDescriptor {
|
||||
label: Some("Multisample color texture"),
|
||||
size: wgpu::Extent3d {
|
||||
width: size.x,
|
||||
height: size.y,
|
||||
depth_or_array_layers: 1,
|
||||
},
|
||||
mip_level_count: 1,
|
||||
sample_count,
|
||||
dimension: wgpu::TextureDimension::D2,
|
||||
format,
|
||||
usage: wgpu::TextureUsages::RENDER_ATTACHMENT,
|
||||
view_formats: &[],
|
||||
});
|
||||
let color = color.create_view(&wgpu::TextureViewDescriptor::default());
|
||||
let depth = device.create_texture(&wgpu::TextureDescriptor {
|
||||
label: Some("Multisample depth texture"),
|
||||
size: wgpu::Extent3d {
|
||||
width: size.x,
|
||||
height: size.y,
|
||||
depth_or_array_layers: 1,
|
||||
},
|
||||
mip_level_count: 1,
|
||||
sample_count,
|
||||
dimension: wgpu::TextureDimension::D2,
|
||||
format: wgpu::TextureFormat::Depth24Plus,
|
||||
usage: wgpu::TextureUsages::RENDER_ATTACHMENT,
|
||||
view_formats: &[],
|
||||
});
|
||||
let depth = depth.create_view(&wgpu::TextureViewDescriptor::default());
|
||||
Multisample { color, depth }
|
||||
}
|
||||
}
|
||||
213
src/fns.rs
213
src/fns.rs
|
|
@ -1,213 +0,0 @@
|
|||
//! Functions used to make metrics.
|
||||
|
||||
use crate::mathx::FloatExt2;
|
||||
|
||||
pub trait Limiter {
|
||||
fn value(&self, x: f32) -> f32;
|
||||
fn derivative(&self, x: f32) -> f32;
|
||||
}
|
||||
|
||||
pub struct SmoothstepLimiter {
|
||||
pub min: f32,
|
||||
pub max: f32,
|
||||
}
|
||||
|
||||
impl Limiter for SmoothstepLimiter {
|
||||
fn value(&self, x: f32) -> f32 {
|
||||
let y = (self.min, self.max).inverse_lerp(x.abs()).clamp(0.0, 1.0);
|
||||
3.0 * y * y - 2.0 * y * y * y
|
||||
}
|
||||
fn derivative(&self, x: f32) -> f32 {
|
||||
if x.abs() > self.min && x.abs() < self.max {
|
||||
let t = (self.min, self.max).inverse_lerp(x.abs());
|
||||
6.0 * x.signum() * t * (1.0 - t) / (self.max - self.min)
|
||||
} else {
|
||||
0.0
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
pub struct SmootherstepLimiter {
|
||||
pub min: f32,
|
||||
pub max: f32,
|
||||
}
|
||||
|
||||
impl Limiter for SmootherstepLimiter {
|
||||
fn value(&self, x: f32) -> f32 {
|
||||
let y = (self.min, self.max).inverse_lerp(x.abs()).clamp(0.0, 1.0);
|
||||
6.0 * y.powi(5) - 15.0 * y.powi(4) + 10.0 * y.powi(3)
|
||||
}
|
||||
fn derivative(&self, x: f32) -> f32 {
|
||||
if x.abs() > self.min && x.abs() < self.max {
|
||||
let t = (self.min, self.max).inverse_lerp(x.abs());
|
||||
30.0 * (t * (1.0 - t)).powi(2) * x.signum() / (self.max - self.min)
|
||||
} else {
|
||||
0.0
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
pub struct QuadraticAccelerator {
|
||||
pub internal: f32,
|
||||
pub external: f32,
|
||||
}
|
||||
|
||||
/// Продолжает функцию f с [-lim, lim] линейно в предположении f(±lim) = ±val, f'(±lim) = 1.
|
||||
fn extend_linear(t: f32, f: impl FnOnce(f32) -> f32, lim: f32, val: f32) -> f32 {
|
||||
if t.abs() <= lim {
|
||||
f(t)
|
||||
} else {
|
||||
t + t.signum() * (val - lim)
|
||||
}
|
||||
}
|
||||
|
||||
/// Продолжает функцию f с [-lim, lim] константой в предположении f(±lim) = val, f'(±lim) = 0.
|
||||
fn extend_const(t: f32, f: impl FnOnce(f32) -> f32, lim: f32, val: f32) -> f32 {
|
||||
if t.abs() <= lim {
|
||||
f(t)
|
||||
} else {
|
||||
val
|
||||
}
|
||||
}
|
||||
|
||||
impl QuadraticAccelerator {
|
||||
fn a(&self) -> f32 {
|
||||
-(self.external - self.internal) / self.internal.powi(2)
|
||||
}
|
||||
fn b(&self) -> f32 {
|
||||
2.0 * self.external / self.internal - 1.0
|
||||
}
|
||||
fn root(&self, x: f32) -> f32 {
|
||||
(self.b().powi(2) + 4.0 * self.a() * x.abs()).sqrt()
|
||||
}
|
||||
|
||||
pub fn x(&self, u: f32) -> f32 {
|
||||
extend_linear(u, |u| (self.a() * u.abs() + self.b()) * u, self.internal, self.external)
|
||||
}
|
||||
pub fn u(&self, x: f32) -> f32 {
|
||||
extend_linear(
|
||||
x,
|
||||
|x| 0.5 * x.signum() * (-self.b() + self.root(x)) / self.a(),
|
||||
self.external,
|
||||
self.internal,
|
||||
)
|
||||
}
|
||||
pub fn dx(&self, u: f32) -> f32 {
|
||||
extend_const(u, |u| 2.0 * self.a() * u.abs() + self.b(), self.internal, 1.0)
|
||||
}
|
||||
pub fn du(&self, x: f32) -> f32 {
|
||||
extend_const(x, |x| 1.0 / self.root(x), self.external, 1.0)
|
||||
}
|
||||
pub fn d2u(&self, x: f32) -> f32 {
|
||||
extend_const(
|
||||
x,
|
||||
|x| -2.0 * x.signum() * self.a() * self.root(x).powi(-3),
|
||||
self.external,
|
||||
0.0,
|
||||
)
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod test {
|
||||
use approx::{abs_diff_eq, assert_abs_diff_eq};
|
||||
|
||||
use super::*;
|
||||
|
||||
fn test_limiter(testee: impl Limiter, min: f32, max: f32, δ: f32) {
|
||||
let ε = 1.0e-4f32;
|
||||
let margin = 1.0 / 16.0;
|
||||
let mul = 1.0 + margin;
|
||||
for x in itertools_num::linspace(0., min, 10) {
|
||||
assert_abs_diff_eq!(testee.value(x), 0., epsilon = ε);
|
||||
assert_abs_diff_eq!(testee.value(-x), 0., epsilon = ε);
|
||||
}
|
||||
for x in itertools_num::linspace(max, mul * max, 10) {
|
||||
assert_abs_diff_eq!(testee.value(x), 1., epsilon = ε);
|
||||
assert_abs_diff_eq!(testee.value(-x), 1., epsilon = ε);
|
||||
}
|
||||
for x in itertools_num::linspace(-mul * max, mul * max, 100) {
|
||||
let df_num = (testee.value(x + δ) - testee.value(x - δ)) / (2. * δ);
|
||||
let df_expl = testee.derivative(x);
|
||||
assert!(
|
||||
abs_diff_eq!(df_expl, df_num, epsilon = ε),
|
||||
"At x={x}, df/dx:\nnumerical: {df_num}\nexplicit: {df_expl}\n"
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_smoothstep_limiter() {
|
||||
test_limiter(SmoothstepLimiter { min: 20.0, max: 30.0 }, 20.0, 30.0, 1.0 / 32.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_smootherstep_limiter() {
|
||||
test_limiter(SmootherstepLimiter { min: 20.0, max: 30.0 }, 20.0, 30.0, 1.0 / 32.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_quadratic_accelerator() {
|
||||
let testee = super::QuadraticAccelerator {
|
||||
internal: 100.0,
|
||||
external: 150.0,
|
||||
};
|
||||
let ε = 1.0e-4f32;
|
||||
let δ = 1.0 / 8.0; // Mathematically, you want this to be small. Computationally, you don’t.
|
||||
let margin = 1.0 / 16.0;
|
||||
let mul = 1.0 + margin;
|
||||
assert_abs_diff_eq!(testee.u(testee.external), testee.internal, epsilon = ε);
|
||||
assert_abs_diff_eq!(testee.u(-testee.external), -testee.internal, epsilon = ε);
|
||||
assert_abs_diff_eq!(testee.du(testee.external), 1., epsilon = ε);
|
||||
assert_abs_diff_eq!(testee.du(-testee.external), 1., epsilon = ε);
|
||||
for x in itertools_num::linspace(-mul * testee.external, mul * testee.external, 100) {
|
||||
let ux = testee.u(x);
|
||||
let xux = testee.x(ux);
|
||||
assert!(
|
||||
abs_diff_eq!(x, xux, epsilon = ε),
|
||||
"At x={x}:\nu(x): {ux}\nx(u(x)): {xux}\n"
|
||||
);
|
||||
|
||||
let du_num = (testee.u(x + δ) - testee.u(x - δ)) / (2. * δ);
|
||||
let du_expl = testee.du(x);
|
||||
assert!(
|
||||
abs_diff_eq!(du_expl, du_num, epsilon = ε),
|
||||
"At x={x}, du/dx:\nnumerical: {du_num}\nexplicit: {du_expl}\n"
|
||||
);
|
||||
|
||||
let dudx = du_expl * testee.dx(ux);
|
||||
assert!(
|
||||
abs_diff_eq!(dudx, 1.0, epsilon = ε),
|
||||
"At x={x}:\ndu/dx * dx/du: {dudx}\n"
|
||||
);
|
||||
|
||||
let d2u_num = (testee.du(x + δ) - testee.du(x - δ)) / (2. * δ);
|
||||
let d2u_expl = testee.d2u(x);
|
||||
assert!(
|
||||
abs_diff_eq!(d2u_expl, d2u_num, epsilon = ε),
|
||||
"At x={x}, d^2u/dx^2:\nnumerical: {d2u_num}\nexplicit: {d2u_expl}\n"
|
||||
);
|
||||
}
|
||||
for u in itertools_num::linspace(-mul * testee.internal, mul * testee.internal, 100) {
|
||||
let xu = testee.x(u);
|
||||
let uxu = testee.u(xu);
|
||||
assert!(
|
||||
abs_diff_eq!(u, uxu, epsilon = ε),
|
||||
"At u={u}:\nx(u): {xu}\nu(x(u)): {uxu}\n"
|
||||
);
|
||||
|
||||
let dx_num = (testee.x(u + δ) - testee.x(u - δ)) / (2. * δ);
|
||||
let dx_expl = testee.dx(u);
|
||||
assert!(
|
||||
abs_diff_eq!(dx_expl, dx_num, epsilon = ε),
|
||||
"At u={u}, dx/du:\nnumerical: {dx_num}\nexplicit: {dx_expl}\n"
|
||||
);
|
||||
|
||||
let dudx = testee.du(xu) * dx_expl;
|
||||
assert!(
|
||||
abs_diff_eq!(dudx, 1.0, epsilon = ε),
|
||||
"At u={u}:\ndu/dx * dx/du: {dudx}\n"
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -1,27 +0,0 @@
|
|||
use crate::types::{Hit, Location, Ray};
|
||||
|
||||
pub trait Traceable {
|
||||
/// Traces a ray from a given starting point. `ray` is relative to the camera.
|
||||
///
|
||||
/// Returns all objects the ray touched. `Hit::distance` is along the ray; it may be a rough estimate except, when two objects overlap the difference of corresponding `distance`s shall be correct.
|
||||
fn trace(&self, camera: Location, ray: Ray) -> Vec<Hit>;
|
||||
}
|
||||
|
||||
pub trait OptimizedTraceable: Traceable {
|
||||
type State;
|
||||
|
||||
/// Prepares tracing from a given starting point. `ray` is relative to the camera.
|
||||
fn init(&self, camera: Location, ray: Ray) -> Self::State;
|
||||
|
||||
/// Similar to [`Traceable::trace`] but allows stopping early.
|
||||
fn trace(&self, state: Self::State) -> (Option<Self::State>, Vec<Hit>);
|
||||
}
|
||||
|
||||
pub struct RayPath {
|
||||
pub points: Vec<Ray>,
|
||||
}
|
||||
|
||||
pub trait DebugTraceable: Traceable {
|
||||
/// Identical to [`Traceable::trace`], except also returns the ray path.
|
||||
fn trace_dbg(&self, camera: Location, ray: Ray) -> (Vec<Hit>, RayPath);
|
||||
}
|
||||
12
src/lib.rs
12
src/lib.rs
|
|
@ -1,12 +0,0 @@
|
|||
mod fns;
|
||||
pub mod ifaces;
|
||||
pub mod mathx;
|
||||
pub mod mesh_loader;
|
||||
pub mod mesh_tracer;
|
||||
pub mod riemann;
|
||||
pub mod shape;
|
||||
pub mod tube;
|
||||
pub mod types;
|
||||
pub mod utils;
|
||||
|
||||
pub const DT: f32 = 0.1;
|
||||
219
src/mathx.rs
219
src/mathx.rs
|
|
@ -1,219 +0,0 @@
|
|||
use glam::{FloatExt, Mat2, Mat3, Vec2, Vec3};
|
||||
|
||||
mod bounds {
|
||||
pub trait Pair<T> {}
|
||||
|
||||
impl<T> Pair<T> for (T, T) {}
|
||||
}
|
||||
|
||||
pub trait FloatExt2<T>: bounds::Pair<T> {
|
||||
fn lerp(self, t: T) -> T;
|
||||
fn inverse_lerp(self, y: T) -> T;
|
||||
}
|
||||
|
||||
impl<F: FloatExt> FloatExt2<F> for (F, F) {
|
||||
fn lerp(self, t: F) -> F {
|
||||
F::lerp(self.0, self.1, t)
|
||||
}
|
||||
fn inverse_lerp(self, y: F) -> F {
|
||||
F::inverse_lerp(self.0, self.1, y)
|
||||
}
|
||||
}
|
||||
|
||||
pub trait MatExt {
|
||||
fn orthonormalize(&self) -> Self;
|
||||
}
|
||||
|
||||
impl MatExt for Mat2 {
|
||||
fn orthonormalize(&self) -> Self {
|
||||
let fx = self.x_axis.normalize();
|
||||
let fy = (self.y_axis - self.y_axis.project_onto_normalized(fx)).normalize();
|
||||
Self::from_cols(fx, fy)
|
||||
}
|
||||
}
|
||||
|
||||
impl MatExt for Mat3 {
|
||||
fn orthonormalize(&self) -> Self {
|
||||
let fx = self.x_axis.normalize();
|
||||
let fy = (self.y_axis - self.y_axis.project_onto_normalized(fx)).normalize();
|
||||
let fz = (self.z_axis - self.z_axis.project_onto_normalized(fx) - self.z_axis.project_onto_normalized(fy))
|
||||
.normalize();
|
||||
Self::from_cols(fx, fy, fz)
|
||||
}
|
||||
}
|
||||
|
||||
/// Represents a 2×2 matrix decomposed as O^T D O, where O is orthogonal and D is diagonal.
|
||||
///
|
||||
/// Not every matrix can be decomposed like this, only that of a symmetric bilinear function.
|
||||
#[derive(Copy, Clone, Debug, PartialEq)]
|
||||
pub struct Decomp2 {
|
||||
/// The orthogonal part.
|
||||
///
|
||||
/// Using a non-orthogonal matrix will yield to incorrect results (but no UB).
|
||||
pub ortho: Mat2,
|
||||
|
||||
/// The diagonal part.
|
||||
pub diag: Vec2,
|
||||
}
|
||||
|
||||
impl Decomp2 {
|
||||
/// Computes the square of this, more efficiently than doing that with a matrix.
|
||||
pub fn square(&self) -> Self {
|
||||
Self {
|
||||
ortho: self.ortho,
|
||||
diag: self.diag * self.diag,
|
||||
}
|
||||
}
|
||||
|
||||
/// Computes the inverse of this, more efficiently than doing that with a matrix.
|
||||
pub fn inverse(&self) -> Self {
|
||||
Self {
|
||||
ortho: self.ortho,
|
||||
diag: Vec2::splat(1.0) / self.diag,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl From<Decomp2> for Mat2 {
|
||||
fn from(value: Decomp2) -> Self {
|
||||
value.ortho.transpose() * Mat2::from_diagonal(value.diag) * value.ortho
|
||||
}
|
||||
}
|
||||
|
||||
impl<T> std::ops::Mul<T> for Decomp2
|
||||
where
|
||||
Mat2: std::ops::Mul<T>,
|
||||
{
|
||||
type Output = <Mat2 as std::ops::Mul<T>>::Output;
|
||||
|
||||
fn mul(self, rhs: T) -> Self::Output {
|
||||
Mat2::from(self) * rhs
|
||||
}
|
||||
}
|
||||
|
||||
/// Represents a 3×3 matrix decomposed as O^T D O, where O is orthogonal and D is diagonal.
|
||||
///
|
||||
/// Not every matrix can be decomposed like this, only that of a symmetric bilinear function.
|
||||
#[derive(Copy, Clone, Debug, PartialEq)]
|
||||
pub struct Decomp3 {
|
||||
/// The orthogonal part.
|
||||
///
|
||||
/// Using a non-orthogonal matrix will yield to incorrect results (but no UB).
|
||||
pub ortho: Mat3,
|
||||
|
||||
/// The diagonal part.
|
||||
pub diag: Vec3,
|
||||
}
|
||||
|
||||
impl Decomp3 {
|
||||
/// Computes the square of this, more efficiently than doing that with a matrix.
|
||||
pub fn square(&self) -> Self {
|
||||
Self {
|
||||
ortho: self.ortho,
|
||||
diag: self.diag * self.diag,
|
||||
}
|
||||
}
|
||||
|
||||
/// Computes the inverse of this, more efficiently than doing that with a matrix.
|
||||
pub fn inverse(&self) -> Self {
|
||||
Self {
|
||||
ortho: self.ortho,
|
||||
diag: Vec3::splat(1.0) / self.diag,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl From<Decomp3> for Mat3 {
|
||||
fn from(value: Decomp3) -> Self {
|
||||
value.ortho.transpose() * Mat3::from_diagonal(value.diag) * value.ortho
|
||||
}
|
||||
}
|
||||
|
||||
impl<T> std::ops::Mul<T> for Decomp3
|
||||
where
|
||||
Mat3: std::ops::Mul<T>,
|
||||
{
|
||||
type Output = <Mat3 as std::ops::Mul<T>>::Output;
|
||||
|
||||
fn mul(self, rhs: T) -> Self::Output {
|
||||
Mat3::from(self) * rhs
|
||||
}
|
||||
}
|
||||
|
||||
pub fn make_vec2(f: impl Fn(usize) -> f32) -> Vec2 {
|
||||
Vec2::from_array(std::array::from_fn(|i| f(i)))
|
||||
}
|
||||
|
||||
pub fn make_mat2(f: impl Fn(usize, usize) -> f32) -> Mat2 {
|
||||
Mat2::from_cols_array_2d(&std::array::from_fn(|i| std::array::from_fn(|j| f(i, j))))
|
||||
}
|
||||
|
||||
pub fn make_vec3(f: impl Fn(usize) -> f32) -> Vec3 {
|
||||
Vec3::from_array(std::array::from_fn(|i| f(i)))
|
||||
}
|
||||
|
||||
pub fn make_mat3(f: impl Fn(usize, usize) -> f32) -> Mat3 {
|
||||
Mat3::from_cols_array_2d(&std::array::from_fn(|i| std::array::from_fn(|j| f(i, j))))
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
use approx::assert_abs_diff_eq;
|
||||
use glam::{mat2, mat3, vec2, vec3};
|
||||
|
||||
#[test]
|
||||
fn test_lerp() {
|
||||
assert_eq!((3., 7.).lerp(-0.5), 1.);
|
||||
assert_eq!((3., 7.).lerp(0.0), 3.);
|
||||
assert_eq!((3., 7.).lerp(0.5), 5.);
|
||||
assert_eq!((3., 7.).lerp(1.0), 7.);
|
||||
assert_eq!((3., 7.).lerp(1.5), 9.);
|
||||
assert_eq!((3., 7.).inverse_lerp(1.), -0.5);
|
||||
assert_eq!((3., 7.).inverse_lerp(3.), 0.0);
|
||||
assert_eq!((3., 7.).inverse_lerp(5.), 0.5);
|
||||
assert_eq!((3., 7.).inverse_lerp(7.), 1.0);
|
||||
assert_eq!((3., 7.).inverse_lerp(9.), 1.5);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_orthonormalize_2d() {
|
||||
assert_abs_diff_eq!(
|
||||
mat2(vec2(1., 0.), vec2(0., 1.)).orthonormalize(),
|
||||
mat2(vec2(1., 0.), vec2(0., 1.)),
|
||||
);
|
||||
assert_abs_diff_eq!(
|
||||
mat2(vec2(2., 0.), vec2(3., 5.)).orthonormalize(),
|
||||
mat2(vec2(1., 0.), vec2(0., 1.)),
|
||||
);
|
||||
assert_abs_diff_eq!(
|
||||
mat2(vec2(0., -3.), vec2(5., 1.)).orthonormalize(),
|
||||
mat2(vec2(0., -1.), vec2(1., 0.)),
|
||||
);
|
||||
assert_abs_diff_eq!(
|
||||
mat2(vec2(3., 4.), vec2(5., 1.)).orthonormalize(),
|
||||
mat2(vec2(0.6, 0.8), vec2(0.8, -0.6)),
|
||||
);
|
||||
assert_abs_diff_eq!(
|
||||
mat2(vec2(3., 4.), vec2(1., 5.)).orthonormalize(),
|
||||
mat2(vec2(0.6, 0.8), vec2(-0.8, 0.6)),
|
||||
);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_orthonormalize_3d() {
|
||||
assert_abs_diff_eq!(
|
||||
mat3(vec3(1., 0., 0.), vec3(0., 1., 0.), vec3(0., 0., 1.)).orthonormalize(),
|
||||
mat3(vec3(1., 0., 0.), vec3(0., 1., 0.), vec3(0., 0., 1.)),
|
||||
);
|
||||
assert_abs_diff_eq!(
|
||||
mat3(vec3(2., 0., 0.), vec3(3., 4., 0.), vec3(5., 6., 7.)).orthonormalize(),
|
||||
mat3(vec3(1., 0., 0.), vec3(0., 1., 0.), vec3(0., 0., 1.)),
|
||||
);
|
||||
assert_abs_diff_eq!(
|
||||
mat3(vec3(0., 5., 0.), vec3(0., 7., 6.), vec3(9., 2., 3.)).orthonormalize(),
|
||||
mat3(vec3(0., 1., 0.), vec3(0., 0., 1.), vec3(1., 0., 0.)),
|
||||
);
|
||||
}
|
||||
}
|
||||
|
|
@ -1,97 +0,0 @@
|
|||
use glam::{vec2, vec3, Vec2, Vec3};
|
||||
use std::io;
|
||||
|
||||
#[derive(Copy, Clone, Debug)]
|
||||
struct ObjVertex {
|
||||
vertex: usize,
|
||||
normal: usize,
|
||||
// tex_coord: usize,
|
||||
}
|
||||
|
||||
#[derive(Copy, Clone, Debug)]
|
||||
struct ObjFace {
|
||||
vertices: [ObjVertex; 3],
|
||||
}
|
||||
|
||||
#[derive(Debug)]
|
||||
struct ObjMesh {
|
||||
vertices: Vec<Vec3>,
|
||||
normals: Vec<Vec3>,
|
||||
tex_coords: Vec<Vec2>,
|
||||
faces: Vec<ObjFace>,
|
||||
}
|
||||
|
||||
impl ObjMesh {
|
||||
fn parse_v2(tokens: &[&str]) -> Vec2 {
|
||||
assert_eq!(tokens.len(), 2);
|
||||
let tokens: Vec<_> = tokens.iter().map(|&t| t.parse().unwrap()).collect();
|
||||
vec2(tokens[0], tokens[1])
|
||||
}
|
||||
|
||||
fn parse_v3(tokens: &[&str]) -> Vec3 {
|
||||
assert_eq!(tokens.len(), 3);
|
||||
let tokens: Vec<_> = tokens.iter().map(|&t| t.parse().unwrap()).collect();
|
||||
vec3(tokens[0], tokens[1], tokens[2])
|
||||
}
|
||||
|
||||
fn parse_fv(desc: &&str) -> ObjVertex {
|
||||
let tokens: Vec<_> = desc.split('/').map(|s| s.parse::<usize>().unwrap() - 1).collect();
|
||||
assert_eq!(tokens.len(), 3);
|
||||
ObjVertex {
|
||||
vertex: tokens[0],
|
||||
// tex_coord: tokens[1],
|
||||
normal: tokens[2],
|
||||
}
|
||||
}
|
||||
|
||||
fn parse_f(tokens: &[&str]) -> ObjFace {
|
||||
let vertices: Vec<_> = tokens.iter().map(ObjMesh::parse_fv).collect();
|
||||
ObjFace {
|
||||
vertices: vertices.as_slice().try_into().unwrap(),
|
||||
}
|
||||
}
|
||||
|
||||
fn read(f: &mut impl io::BufRead) -> io::Result<ObjMesh> {
|
||||
let mut result = ObjMesh {
|
||||
vertices: Vec::new(),
|
||||
normals: Vec::new(),
|
||||
tex_coords: Vec::new(),
|
||||
faces: Vec::new(),
|
||||
};
|
||||
loop {
|
||||
let mut line = String::new();
|
||||
if f.read_line(&mut line)? == 0 {
|
||||
break;
|
||||
}
|
||||
let tokens: Vec<&str> = line.trim().split('#').next().unwrap().split(' ').collect();
|
||||
match tokens[0] {
|
||||
"v" => result.vertices.push(Self::parse_v3(&tokens[1..])),
|
||||
"vn" => result.normals.push(Self::parse_v3(&tokens[1..])),
|
||||
"vt" => result.tex_coords.push(Self::parse_v2(&tokens[1..])),
|
||||
"f" => result.faces.push(Self::parse_f(&tokens[1..])),
|
||||
_ => (),
|
||||
}
|
||||
}
|
||||
Ok(result)
|
||||
}
|
||||
|
||||
fn flatten(&self) -> Vec<Face> {
|
||||
self.faces
|
||||
.iter()
|
||||
.map(|face| Face {
|
||||
vertices: face.vertices.map(|iv| self.vertices[iv.vertex]),
|
||||
normal: self.normals[face.vertices[0].normal],
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Copy, Clone, Debug)]
|
||||
pub struct Face {
|
||||
pub vertices: [Vec3; 3],
|
||||
pub normal: Vec3,
|
||||
}
|
||||
|
||||
pub fn load_mesh(f: &mut impl io::BufRead) -> io::Result<Vec<Face>> {
|
||||
Ok(ObjMesh::read(f)?.flatten())
|
||||
}
|
||||
|
|
@ -1,36 +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()
|
||||
}
|
||||
232
src/riemann.rs
232
src/riemann.rs
|
|
@ -1,232 +0,0 @@
|
|||
use crate::mathx::{make_mat3, Decomp3};
|
||||
use glam::*;
|
||||
|
||||
pub type Tens3 = [Mat3; 3];
|
||||
|
||||
pub trait Metric {
|
||||
fn sqrt_at(&self, pos: Vec3) -> Decomp3;
|
||||
|
||||
fn at(&self, pos: Vec3) -> Mat3 {
|
||||
self.sqrt_at(pos).square().into()
|
||||
}
|
||||
|
||||
fn inverse_at(&self, pos: Vec3) -> Mat3 {
|
||||
self.sqrt_at(pos).square().inverse().into()
|
||||
}
|
||||
|
||||
fn part_derivs_at(&self, pos: Vec3) -> Tens3 {
|
||||
part_deriv(|p| self.at(p), pos, 1.0 / 64.0) // there just isn’t enough precision for a small ε.
|
||||
}
|
||||
|
||||
fn vec_length_at(&self, at: Vec3, v: Vec3) -> f32 {
|
||||
v.dot(self.at(at) * v).sqrt()
|
||||
}
|
||||
|
||||
fn normalize_vec_at(&self, at: Vec3, v: Vec3) -> Vec3 {
|
||||
v / self.vec_length_at(at, v)
|
||||
}
|
||||
}
|
||||
|
||||
pub struct TraceIter<'a, M: Metric> {
|
||||
space: &'a M,
|
||||
p: Vec3,
|
||||
v: Vec3,
|
||||
}
|
||||
|
||||
impl<'a, M: Metric> Iterator for TraceIter<'a, M> {
|
||||
type Item = Vec3;
|
||||
|
||||
fn next(&mut self) -> Option<Self::Item> {
|
||||
let a: Vec3 = -contract2(krist(self.space, self.p), self.v);
|
||||
self.v += a;
|
||||
self.p += self.v;
|
||||
Some(self.p)
|
||||
}
|
||||
}
|
||||
|
||||
pub fn trace_iter<M: Metric>(space: &M, base: Vec3, dir: Vec3, dt: f32) -> TraceIter<M> {
|
||||
TraceIter {
|
||||
space,
|
||||
p: base,
|
||||
v: dt * space.normalize_vec_at(base, dir),
|
||||
}
|
||||
}
|
||||
|
||||
pub fn krist(space: &impl Metric, pos: Vec3) -> Tens3 {
|
||||
// Γ^i_k_l = .5 * g^i^m * (g_m_k,l + g_m_l,k - g_k_l,m)
|
||||
let g = &space.inverse_at(pos); // с верхними индексами
|
||||
let d = space.part_derivs_at(pos);
|
||||
// ret[i][l][k] = sum((m) => .5f * g[m][i] * (d[k][l][m] + d[l][k][m] - d[m][k][l]))
|
||||
make_tens3(|i, l, k| {
|
||||
0.5 * (0..3)
|
||||
.map(|m| g.col(m)[i] * (d[l].col(k)[m] + d[k].col(m)[l] - d[m].col(k)[l]))
|
||||
.sum::<f32>()
|
||||
})
|
||||
}
|
||||
|
||||
fn dir_deriv(f: impl Fn(Vec3) -> Mat3, pos: Vec3, delta: Vec3) -> Mat3 {
|
||||
(f(pos + delta) - f(pos - delta)) / (2.0 * delta.length())
|
||||
}
|
||||
|
||||
fn part_deriv(f: impl Fn(Vec3) -> Mat3, pos: Vec3, eps: f32) -> Tens3 {
|
||||
[
|
||||
dir_deriv(&f, pos, vec3(eps, 0.0, 0.0)),
|
||||
dir_deriv(&f, pos, vec3(0.0, eps, 0.0)),
|
||||
dir_deriv(&f, pos, vec3(0.0, 0.0, eps)),
|
||||
]
|
||||
}
|
||||
|
||||
/// Сворачивает тензор t с вектором u
|
||||
pub fn contract(t: Tens3, u: Vec3) -> Mat3 {
|
||||
mat3(t[0] * u, t[1] * u, t[2] * u).transpose()
|
||||
}
|
||||
|
||||
/// Сворачивает тензор t с вектором v дважды, по второму и третьему индексам.
|
||||
pub fn contract2(t: Tens3, v: Vec3) -> Vec3 {
|
||||
contract(t, v) * v
|
||||
}
|
||||
|
||||
fn make_tens3(f: impl Fn(usize, usize, usize) -> f32) -> Tens3 {
|
||||
std::array::from_fn(|i| make_mat3(|j, k| f(i, j, k)))
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn m3() {
|
||||
let m = make_mat3(|i, j| (i + 2 * j) as f32);
|
||||
assert_eq!(m.col(0)[0], 0.0);
|
||||
assert_eq!(m.col(1)[0], 1.0);
|
||||
assert_eq!(m.col(0)[1], 2.0);
|
||||
assert_eq!(m.col(1)[1], 3.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn t3() {
|
||||
let t = make_tens3(|i, j, k| (i + 2 * j + 4 * k) as f32);
|
||||
assert_eq!(t[0].col(0)[0], 0.0);
|
||||
assert_eq!(t[1].col(0)[0], 1.0);
|
||||
assert_eq!(t[0].col(1)[0], 2.0);
|
||||
assert_eq!(t[1].col(1)[0], 3.0);
|
||||
assert_eq!(t[0].col(0)[1], 4.0);
|
||||
assert_eq!(t[1].col(0)[1], 5.0);
|
||||
assert_eq!(t[0].col(1)[1], 6.0);
|
||||
assert_eq!(t[1].col(1)[1], 7.0);
|
||||
}
|
||||
|
||||
pub mod samples {
|
||||
use glam::{Mat3, Vec3};
|
||||
|
||||
use super::{Decomp3, Metric};
|
||||
|
||||
pub struct ScaledMetric {
|
||||
/// Specifies unit size in each cardinal direction. E.g. with scale=(2, 3), vector (1, 0) has length 2 while a unit vector with the same direction is (1/2, 0).
|
||||
pub scale: Vec3,
|
||||
}
|
||||
|
||||
impl Metric for ScaledMetric {
|
||||
fn sqrt_at(&self, _pos: Vec3) -> Decomp3 {
|
||||
Decomp3 {
|
||||
diag: self.scale,
|
||||
ortho: Mat3::IDENTITY,
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use approx::assert_abs_diff_eq;
|
||||
|
||||
use glam::{vec3, Mat3};
|
||||
use rand::{Rng, SeedableRng};
|
||||
|
||||
#[test]
|
||||
fn uniform_scaled_metric() {
|
||||
let mut rng = rand_pcg::Pcg64Mcg::seed_from_u64(17);
|
||||
let metric = samples::ScaledMetric {
|
||||
scale: vec3(3., 4., 5.),
|
||||
};
|
||||
assert_eq!(
|
||||
metric.sqrt_at(rng.gen()),
|
||||
Decomp3 {
|
||||
ortho: Mat3::IDENTITY,
|
||||
diag: vec3(3., 4., 5.)
|
||||
}
|
||||
);
|
||||
assert_eq!(
|
||||
metric.at(rng.gen()),
|
||||
Mat3::from_cols_array(&[9., 0., 0., 0., 16., 0., 0., 0., 25.])
|
||||
);
|
||||
assert_eq!(
|
||||
metric.inverse_at(rng.gen()),
|
||||
Mat3::from_cols_array(&[1. / 9., 0., 0., 0., 1. / 16., 0., 0., 0., 1. / 25.])
|
||||
);
|
||||
assert_eq!(metric.part_derivs_at(rng.gen()), [Mat3::ZERO, Mat3::ZERO, Mat3::ZERO]);
|
||||
assert_eq!(metric.vec_length_at(rng.gen(), vec3(1., 0., 0.)), 3.);
|
||||
assert_eq!(metric.vec_length_at(rng.gen(), vec3(0., 1., 0.)), 4.);
|
||||
assert_eq!(metric.vec_length_at(rng.gen(), vec3(0., 0., 1.)), 5.);
|
||||
assert_eq!(metric.vec_length_at(rng.gen(), vec3(1., 1., 0.)), 5.);
|
||||
assert_eq!(
|
||||
metric.normalize_vec_at(rng.gen(), vec3(1., 0., 0.)),
|
||||
vec3(1. / 3., 0., 0.)
|
||||
);
|
||||
assert_eq!(
|
||||
metric.normalize_vec_at(rng.gen(), vec3(0., 1., 0.)),
|
||||
vec3(0., 1. / 4., 0.)
|
||||
);
|
||||
assert_eq!(
|
||||
metric.normalize_vec_at(rng.gen(), vec3(1., 1., 0.)),
|
||||
vec3(1. / 5., 1. / 5., 0.)
|
||||
);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_trace_iter() {
|
||||
let metric = samples::ScaledMetric {
|
||||
scale: vec3(2., 4., 3.),
|
||||
};
|
||||
assert_eq!(
|
||||
trace_iter(&metric, vec3(3., 5., 0.), vec3(1., 0., 0.), 1.)
|
||||
.nth(7)
|
||||
.unwrap(),
|
||||
vec3(7., 5., 0.)
|
||||
);
|
||||
assert_eq!(
|
||||
trace_iter(&metric, vec3(3., 5., 0.), vec3(2., 0., 0.), 1.)
|
||||
.nth(7)
|
||||
.unwrap(),
|
||||
vec3(7., 5., 0.)
|
||||
);
|
||||
assert_eq!(
|
||||
trace_iter(&metric, vec3(3., 5., 0.), vec3(1., 0., 0.), 0.5)
|
||||
.nth(7)
|
||||
.unwrap(),
|
||||
vec3(5., 5., 0.)
|
||||
);
|
||||
assert_eq!(
|
||||
trace_iter(&metric, vec3(3., 5., 0.), vec3(0., 1., 0.), 1.)
|
||||
.nth(9)
|
||||
.unwrap(),
|
||||
vec3(3., 7.5, 0.)
|
||||
);
|
||||
assert_eq!(
|
||||
trace_iter(&metric, vec3(3., 5., 0.), vec3(0., 4., 0.), 1.)
|
||||
.nth(9)
|
||||
.unwrap(),
|
||||
vec3(3., 7.5, 0.)
|
||||
);
|
||||
assert_eq!(
|
||||
trace_iter(&metric, vec3(3., 5., 0.), vec3(0., 1., 0.), 0.5)
|
||||
.nth(9)
|
||||
.unwrap(),
|
||||
vec3(3., 6.25, 0.)
|
||||
);
|
||||
assert_abs_diff_eq!(
|
||||
trace_iter(&metric, vec3(3., 5., 0.), vec3(0.5, 0.25, 0.), std::f32::consts::SQRT_2)
|
||||
.nth(7)
|
||||
.unwrap(),
|
||||
vec3(7., 7., 0.),
|
||||
epsilon = 1e-5
|
||||
);
|
||||
}
|
||||
}
|
||||
|
|
@ -1,240 +0,0 @@
|
|||
use glam::{Vec3, Vec3Swizzles as _};
|
||||
|
||||
use crate::types::Ray;
|
||||
|
||||
pub struct Cylinder {
|
||||
pub center: Vec3,
|
||||
pub semiaxis: Vec3,
|
||||
pub radius: f32,
|
||||
}
|
||||
|
||||
impl Cylinder {
|
||||
/// Split a vector into a component along the axis and one orthogonal to it.
|
||||
fn split(&self, dir: Vec3) -> (f32, Vec3) {
|
||||
let along = dir.dot(self.semiaxis) / self.semiaxis.length_squared();
|
||||
(along, dir - along * self.semiaxis)
|
||||
}
|
||||
|
||||
fn cap_inout(p: f32, d: f32) -> (f32, f32) {
|
||||
((-d.signum() - p) / d.abs(), (d.signum() - p) / d.abs())
|
||||
}
|
||||
|
||||
pub fn is_inside(&self, pt: Vec3) -> bool {
|
||||
let (along, ortho) = self.split(pt - self.center);
|
||||
along.abs() < 1. && ortho.length_squared() < self.radius.powi(2)
|
||||
}
|
||||
|
||||
fn trace_inout(&self, ray: Ray) -> Option<(f32, f32)> {
|
||||
let (pos_along, pos_ortho) = self.split(ray.pos - self.center);
|
||||
let (dir_along, dir_ortho) = self.split(ray.dir);
|
||||
let (t_cap_in, t_cap_out) = Self::cap_inout(pos_along, dir_along);
|
||||
if dir_ortho.length_squared() < 1e-3 {
|
||||
if pos_ortho.length_squared() >= self.radius.powi(2) {
|
||||
return None;
|
||||
}
|
||||
return Some((t_cap_in, t_cap_out));
|
||||
}
|
||||
let (t_side_in, t_side_out) = solve_quadratic(
|
||||
dir_ortho.length_squared(),
|
||||
pos_ortho.dot(dir_ortho),
|
||||
pos_ortho.length_squared() - self.radius.powi(2),
|
||||
)?;
|
||||
let t_in = f32::max(t_cap_in, t_side_in);
|
||||
let t_out = f32::min(t_cap_out, t_side_out);
|
||||
if t_out <= t_in {
|
||||
return None;
|
||||
}
|
||||
Some((t_in, t_out))
|
||||
}
|
||||
|
||||
pub fn trace_into(&self, ray: Ray) -> Option<f32> {
|
||||
let (t, _) = self.trace_inout(ray)?;
|
||||
if t < 0. {
|
||||
return None;
|
||||
}
|
||||
Some(t)
|
||||
}
|
||||
|
||||
pub fn trace_out_of(&self, ray: Ray) -> Option<f32> {
|
||||
let (_, t) = self
|
||||
.trace_inout(ray)
|
||||
.expect("the ray starts inside so *has* to cross the boundary");
|
||||
Some(t)
|
||||
}
|
||||
}
|
||||
|
||||
/// Цилиндр с центром в начале координат и осью вдоль оси Y.
|
||||
pub struct YCylinder {
|
||||
pub half_length: f32,
|
||||
pub radius: f32,
|
||||
}
|
||||
|
||||
impl YCylinder {
|
||||
/// Отражает луч, чтобы все координаты направления были положительны (допустимо благодаря симметрии YCylinder).
|
||||
fn flip_ray(ray: Ray) -> Ray {
|
||||
Ray {
|
||||
pos: ray.pos * ray.dir.signum(),
|
||||
dir: ray.dir.abs(),
|
||||
}
|
||||
}
|
||||
|
||||
pub fn is_inside(&self, pt: Vec3) -> bool {
|
||||
let r = f32::hypot(pt.x, pt.z);
|
||||
pt.y.abs() < self.half_length && r < self.radius
|
||||
}
|
||||
|
||||
pub fn trace_into(&self, ray: Ray) -> Option<f32> {
|
||||
let ray = Self::flip_ray(ray);
|
||||
|
||||
// 1. ray.pos.y + t * ray.dir.y = −half_length
|
||||
let t_cap_in = (-self.half_length - ray.pos.y) / ray.dir.y;
|
||||
let t_cap_out = (self.half_length - ray.pos.y) / ray.dir.y;
|
||||
|
||||
// 2. (ray.pos.x + t * ray.dir.x)² + (ray.pos.z + t * ray.dir.z)² = radius²
|
||||
let pos = ray.pos.xz();
|
||||
let dir = ray.dir.xz();
|
||||
if dir.length_squared() < 1e-6 * ray.dir.length_squared() {
|
||||
if pos.length_squared() >= self.radius.powi(2) {
|
||||
return None;
|
||||
}
|
||||
return Some(t_cap_in).filter(|&t| t > 0.);
|
||||
}
|
||||
let (t_side_in, t_side_out) = solve_quadratic(
|
||||
dir.length_squared(),
|
||||
pos.dot(dir),
|
||||
pos.length_squared() - self.radius.powi(2),
|
||||
)?;
|
||||
let t = f32::max(t_cap_in, t_side_in);
|
||||
if t < 0. {
|
||||
return None;
|
||||
}
|
||||
if t >= t_cap_out || t >= t_side_out {
|
||||
return None;
|
||||
}
|
||||
Some(t)
|
||||
}
|
||||
|
||||
pub fn trace_out_of(&self, ray: Ray) -> Option<f32> {
|
||||
let ray = Self::flip_ray(ray);
|
||||
let t_cap_out = (self.half_length - ray.pos.y) / ray.dir.y;
|
||||
let pos = ray.pos.xz();
|
||||
let dir = ray.dir.xz();
|
||||
if dir.length_squared() < 1e-3 {
|
||||
return Some(t_cap_out);
|
||||
}
|
||||
let (_t_side_in, t_side_out) = solve_quadratic(
|
||||
dir.length_squared(),
|
||||
pos.dot(dir),
|
||||
pos.length_squared() - self.radius.powi(2),
|
||||
)
|
||||
.expect("the ray starts inside and is not along the axis so *has* to cross the side");
|
||||
Some(t_side_out)
|
||||
}
|
||||
}
|
||||
|
||||
fn solve_quadratic(a: f32, half_b: f32, c: f32) -> Option<(f32, f32)> {
|
||||
let base = -half_b / a;
|
||||
let d = base * base - c / a;
|
||||
if d < 0. {
|
||||
None
|
||||
} else {
|
||||
let δ = d.sqrt();
|
||||
Some((base - δ, base + δ))
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use crate::types::ray;
|
||||
use approx::assert_abs_diff_eq;
|
||||
use glam::vec3;
|
||||
use rand::{Rng, SeedableRng};
|
||||
|
||||
#[test]
|
||||
fn test_split() {
|
||||
let mut rng = rand_pcg::Pcg64Mcg::seed_from_u64(17);
|
||||
let cyl = Cylinder {
|
||||
center: vec3(1., 2., 3.),
|
||||
semiaxis: vec3(4., 5., 6.),
|
||||
radius: 7.,
|
||||
};
|
||||
for _ in 0..100 {
|
||||
let dir = vec3(rng.gen(), rng.gen(), rng.gen());
|
||||
let (along, ortho) = cyl.split(dir);
|
||||
assert_abs_diff_eq!(along * cyl.semiaxis + ortho, dir, epsilon = 1e-5);
|
||||
assert_abs_diff_eq!(cyl.semiaxis.dot(ortho), 0., epsilon = 1e-5);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_cylinder() {
|
||||
assert_eq!(
|
||||
YCylinder::flip_ray(ray(vec3(2., 3., 2.), vec3(4., 5., 4.))),
|
||||
ray(vec3(2., 3., 2.), vec3(4., 5., 4.)),
|
||||
);
|
||||
assert_eq!(
|
||||
YCylinder::flip_ray(ray(vec3(2., 3., 2.), vec3(-4., 5., -4.))),
|
||||
ray(vec3(-2., 3., -2.), vec3(4., 5., 4.)),
|
||||
);
|
||||
assert_eq!(
|
||||
YCylinder::flip_ray(ray(vec3(2., 3., 2.), vec3(4., -5., 4.))),
|
||||
ray(vec3(2., -3., 2.), vec3(4., 5., 4.)),
|
||||
);
|
||||
assert_eq!(
|
||||
YCylinder::flip_ray(ray(vec3(2., 3., 2.), vec3(4., 0., 4.))),
|
||||
ray(vec3(2., 3., 2.), vec3(4., 0., 4.)),
|
||||
);
|
||||
|
||||
let r = YCylinder {
|
||||
half_length: 3.,
|
||||
radius: 2.,
|
||||
};
|
||||
|
||||
assert_eq!(r.trace_into(ray(vec3(3., 4., 3.), vec3(0., -1., 0.))), None);
|
||||
assert_eq!(r.trace_into(ray(vec3(1., 4., 1.), vec3(0., -1., 0.))), Some(1.));
|
||||
assert_eq!(r.trace_into(ray(vec3(3., 3., 3.), vec3(1., 1., 1.))), None);
|
||||
assert_abs_diff_eq!(
|
||||
r.trace_into(ray(vec3(-3., 2., -3.), vec3(1., 0., 1.))).unwrap(),
|
||||
1.5857864
|
||||
);
|
||||
assert_eq!(r.trace_into(ray(vec3(-3., 2., -3.), vec3(-1., 0., -1.))), None);
|
||||
assert_abs_diff_eq!(
|
||||
r.trace_into(ray(vec3(-3., 1., -3.), vec3(2., 2., 2.))).unwrap(),
|
||||
0.7928932
|
||||
);
|
||||
assert_eq!(r.trace_into(ray(vec3(-3., 2.1, -3.), vec3(2., 2., 2.))), None);
|
||||
|
||||
assert_eq!(r.trace_into(ray(vec3(2., 3., 2.), vec3(1., 1., 1.))), None);
|
||||
assert_eq!(r.trace_into(ray(vec3(-2., 3., -2.), vec3(-1., 1., -1.))), None);
|
||||
assert_eq!(
|
||||
r.trace_into(ray(vec3(1.4142135, 3., 1.4142135), vec3(-1., -1., -1.))),
|
||||
Some(0.)
|
||||
);
|
||||
assert_eq!(
|
||||
r.trace_into(ray(vec3(1.4142135, -3., 1.4142135), vec3(-1., 1., -1.))),
|
||||
Some(0.)
|
||||
);
|
||||
assert_eq!(
|
||||
YCylinder {
|
||||
half_length: 300.,
|
||||
radius: 50.
|
||||
}
|
||||
.trace_into(ray(vec3(-125., 375., 0.), vec3(3., -11., 0.) / 1024.)),
|
||||
Some(25600.)
|
||||
);
|
||||
|
||||
assert_abs_diff_eq!(
|
||||
r.trace_out_of(ray(vec3(0., 0., 0.), vec3(1., 1., 1.))).unwrap(),
|
||||
1.4142135
|
||||
);
|
||||
assert_eq!(r.trace_out_of(ray(vec3(0., 0., 0.), vec3(0., 1., 0.))), Some(3.));
|
||||
assert_eq!(r.trace_out_of(ray(vec3(0., 1., 0.), vec3(0., -1., 0.))), Some(4.));
|
||||
assert_eq!(r.trace_out_of(ray(vec3(1., 1., 1.), vec3(0., -1., 0.))), Some(4.));
|
||||
assert_abs_diff_eq!(
|
||||
r.trace_out_of(ray(vec3(1.4142135, 3., 1.4142135), vec3(1., 1., 1.)))
|
||||
.unwrap(),
|
||||
0.
|
||||
);
|
||||
}
|
||||
}
|
||||
|
|
@ -1,5 +0,0 @@
|
|||
pub mod cylinder;
|
||||
pub mod rect;
|
||||
|
||||
pub use cylinder::YCylinder;
|
||||
pub use rect::Rect;
|
||||
|
|
@ -1,90 +0,0 @@
|
|||
use glam::Vec3;
|
||||
|
||||
use crate::types::Ray;
|
||||
|
||||
pub struct Rect {
|
||||
pub size: Vec3,
|
||||
}
|
||||
|
||||
impl Rect {
|
||||
/// Отражает луч, чтобы все координаты направления были положительны (допустимо благодаря симметрии Rect).
|
||||
fn flip_ray(ray: Ray) -> Ray {
|
||||
Ray {
|
||||
pos: ray.pos * ray.dir.signum(),
|
||||
dir: ray.dir.abs(),
|
||||
}
|
||||
}
|
||||
|
||||
pub fn is_inside(&self, pt: Vec3) -> bool {
|
||||
pt.abs().cmplt(self.size).all()
|
||||
}
|
||||
|
||||
pub fn trace_into(&self, ray: Ray) -> Option<f32> {
|
||||
let ray = Self::flip_ray(ray);
|
||||
// ray.pos.x + t * ray.dir.x = −size.x
|
||||
let ts = (-self.size - ray.pos) / ray.dir;
|
||||
let t = ts.max_element();
|
||||
let pt = ray.pos + t * ray.dir;
|
||||
if t < 0.0 {
|
||||
return None;
|
||||
}
|
||||
if pt.cmpgt(self.size).any() {
|
||||
return None;
|
||||
}
|
||||
Some(t)
|
||||
}
|
||||
|
||||
pub fn trace_out_of(&self, ray: Ray) -> Option<f32> {
|
||||
let ray = Self::flip_ray(ray);
|
||||
// ray.pos.x + t * ray.dir.x = +size.x
|
||||
let ts = (self.size - ray.pos) / ray.dir;
|
||||
let t = ts.min_element();
|
||||
Some(t)
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use crate::types::ray;
|
||||
use glam::vec3;
|
||||
|
||||
#[test]
|
||||
fn test_rect() {
|
||||
assert_eq!(
|
||||
Rect::flip_ray(ray(vec3(2., 3., 2.), vec3(4., 5., 4.))),
|
||||
ray(vec3(2., 3., 2.), vec3(4., 5., 4.)),
|
||||
);
|
||||
assert_eq!(
|
||||
Rect::flip_ray(ray(vec3(2., 3., 2.), vec3(-4., 5., -4.))),
|
||||
ray(vec3(-2., 3., -2.), vec3(4., 5., 4.)),
|
||||
);
|
||||
assert_eq!(
|
||||
Rect::flip_ray(ray(vec3(2., 3., 2.), vec3(4., -5., 4.))),
|
||||
ray(vec3(2., -3., 2.), vec3(4., 5., 4.)),
|
||||
);
|
||||
assert_eq!(
|
||||
Rect::flip_ray(ray(vec3(2., 3., 2.), vec3(4., 0., 4.))),
|
||||
ray(vec3(2., 3., 2.), vec3(4., 0., 4.)),
|
||||
);
|
||||
|
||||
let r = Rect { size: vec3(2., 3., 2.) };
|
||||
|
||||
assert_eq!(r.trace_into(ray(vec3(3., 3., 3.), vec3(1., 1., 1.))), None);
|
||||
assert_eq!(r.trace_into(ray(vec3(-3., 2., -3.), vec3(1., 0., 1.))), Some(1.));
|
||||
assert_eq!(r.trace_into(ray(vec3(-3., 2., -3.), vec3(-1., 0., -1.))), None);
|
||||
assert_eq!(r.trace_into(ray(vec3(-3., 1., -3.), vec3(2., 2., 2.))), Some(0.5));
|
||||
assert_eq!(r.trace_into(ray(vec3(-3., 2.1, -3.), vec3(2., 2., 2.))), None);
|
||||
|
||||
assert_eq!(r.trace_into(ray(vec3(2., 3., 2.), vec3(1., 1., 1.))), None);
|
||||
assert_eq!(r.trace_into(ray(vec3(-2., 3., -2.), vec3(-1., 1., -1.))), None);
|
||||
assert_eq!(r.trace_into(ray(vec3(2., 3., 2.), vec3(-1., -1., -1.))), Some(0.));
|
||||
assert_eq!(r.trace_into(ray(vec3(2., -3., 2.), vec3(-1., 1., -1.))), Some(0.));
|
||||
|
||||
assert_eq!(r.trace_out_of(ray(vec3(0., 0., 0.), vec3(1., 1., 1.))), Some(2.));
|
||||
assert_eq!(r.trace_out_of(ray(vec3(0., 0., 0.), vec3(0., 1., 0.))), Some(3.));
|
||||
assert_eq!(r.trace_out_of(ray(vec3(0., 1., 0.), vec3(0., -1., 0.))), Some(4.));
|
||||
assert_eq!(r.trace_out_of(ray(vec3(1., 1., 1.), vec3(0., -1., 0.))), Some(4.));
|
||||
assert_eq!(r.trace_out_of(ray(vec3(2., 3., 2.), vec3(1., 1., 1.))), Some(0.));
|
||||
}
|
||||
}
|
||||
|
|
@ -1,468 +0,0 @@
|
|||
use glam::{vec3, Mat3, Vec3};
|
||||
|
||||
use crate::riemann::Metric;
|
||||
use crate::shape::YCylinder;
|
||||
use crate::types::{Location, Ray};
|
||||
|
||||
use super::Tube;
|
||||
|
||||
pub trait FlatCoordinateSystem<T> {
|
||||
fn flat_to_global(&self, v: T) -> T;
|
||||
fn global_to_flat(&self, v: T) -> T;
|
||||
}
|
||||
|
||||
pub trait FlatRegion: FlatCoordinateSystem<Vec3> + FlatCoordinateSystem<Ray> + FlatCoordinateSystem<Location> {
|
||||
// Измеряет расстояние до выхода за пределы области вдоль луча ray. Луч задаётся в плоской СК.
|
||||
fn distance_to_boundary(&self, _ray: Ray) -> Option<f32> {
|
||||
None
|
||||
}
|
||||
}
|
||||
|
||||
trait MetricCS: FlatCoordinateSystem<Vec3> {
|
||||
fn global_metric(&self) -> &impl Metric;
|
||||
fn flat_to_global_tfm_at(&self, pos: Vec3) -> Mat3 {
|
||||
self.global_metric().sqrt_at(self.flat_to_global(pos)).inverse().into()
|
||||
}
|
||||
fn global_to_flat_tfm_at(&self, pos: Vec3) -> Mat3 {
|
||||
self.global_metric().sqrt_at(pos).into()
|
||||
}
|
||||
}
|
||||
|
||||
impl<T: FlatCoordinateSystem<Vec3> + MetricCS> FlatCoordinateSystem<Ray> for T {
|
||||
fn flat_to_global(&self, ray: Ray) -> Ray {
|
||||
Ray {
|
||||
pos: self.flat_to_global(ray.pos),
|
||||
dir: self.flat_to_global_tfm_at(ray.pos) * ray.dir,
|
||||
}
|
||||
}
|
||||
|
||||
fn global_to_flat(&self, ray: Ray) -> Ray {
|
||||
Ray {
|
||||
pos: self.global_to_flat(ray.pos),
|
||||
dir: self.global_to_flat_tfm_at(ray.pos) * ray.dir,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl<T: FlatCoordinateSystem<Vec3> + MetricCS> FlatCoordinateSystem<Location> for T {
|
||||
fn flat_to_global(&self, loc: Location) -> Location {
|
||||
Location {
|
||||
pos: self.flat_to_global(loc.pos),
|
||||
rot: self.flat_to_global_tfm_at(loc.pos) * loc.rot,
|
||||
}
|
||||
}
|
||||
|
||||
fn global_to_flat(&self, loc: Location) -> Location {
|
||||
Location {
|
||||
pos: self.global_to_flat(loc.pos), // в плоской СК для Inner или её продолжении на Outer
|
||||
rot: self.global_to_flat_tfm_at(loc.pos) * loc.rot,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
pub struct InnerCS(pub Tube);
|
||||
|
||||
impl MetricCS for InnerCS {
|
||||
fn global_metric(&self) -> &impl Metric {
|
||||
&self.0
|
||||
}
|
||||
}
|
||||
|
||||
impl FlatCoordinateSystem<Vec3> for InnerCS {
|
||||
fn flat_to_global(&self, pos: Vec3) -> Vec3 {
|
||||
vec3(pos.x, self.0.y(pos.y), pos.z)
|
||||
}
|
||||
|
||||
// Работает только при |pos.x| ≤ inner_radius или |pos.y| ≥ external_halflength.
|
||||
fn global_to_flat(&self, pos: Vec3) -> Vec3 {
|
||||
vec3(pos.x, self.0.v(pos.y), pos.z)
|
||||
}
|
||||
}
|
||||
|
||||
impl FlatRegion for InnerCS {
|
||||
fn distance_to_boundary(&self, ray: Ray) -> Option<f32> {
|
||||
YCylinder {
|
||||
radius: self.0.inner_radius,
|
||||
half_length: self.0.internal_halflength,
|
||||
}
|
||||
.trace_out_of(ray)
|
||||
}
|
||||
}
|
||||
|
||||
pub struct OuterCS(pub Tube);
|
||||
|
||||
impl MetricCS for OuterCS {
|
||||
fn global_metric(&self) -> &impl Metric {
|
||||
&self.0
|
||||
}
|
||||
}
|
||||
|
||||
impl FlatCoordinateSystem<Vec3> for OuterCS {
|
||||
fn flat_to_global(&self, pos: Vec3) -> Vec3 {
|
||||
let inner = YCylinder {
|
||||
radius: self.0.inner_radius + 1.0,
|
||||
half_length: self.0.external_halflength,
|
||||
};
|
||||
if inner.is_inside(pos) {
|
||||
let Vec3 { x, y: v, z } = pos;
|
||||
let y = self
|
||||
.0
|
||||
.y(v - v.signum() * (self.0.external_halflength - self.0.internal_halflength));
|
||||
vec3(x, y, z)
|
||||
} else {
|
||||
pos
|
||||
}
|
||||
}
|
||||
|
||||
fn global_to_flat(&self, pos: Vec3) -> Vec3 {
|
||||
let inner = YCylinder {
|
||||
radius: self.0.inner_radius + 1.0,
|
||||
half_length: self.0.external_halflength,
|
||||
};
|
||||
if inner.is_inside(pos) {
|
||||
let Vec3 { x: u, y, z: w } = pos; // в основной СК
|
||||
let v = self.0.v(y) + y.signum() * (self.0.external_halflength - self.0.internal_halflength);
|
||||
vec3(u, v, w) // в плоском продолжении СК Outer на область Inner
|
||||
} else {
|
||||
pos
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl FlatRegion for OuterCS {
|
||||
fn distance_to_boundary(&self, ray: Ray) -> Option<f32> {
|
||||
YCylinder {
|
||||
radius: self.0.outer_radius,
|
||||
half_length: self.0.external_halflength,
|
||||
}
|
||||
.trace_into(ray)
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use approx::{assert_abs_diff_eq, AbsDiffEq};
|
||||
use glam::{mat3, vec3, Mat3, Vec3};
|
||||
use itertools_num::linspace;
|
||||
|
||||
use crate::riemann::samples;
|
||||
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn uniform_scaled_metric() {
|
||||
struct Scaled(samples::ScaledMetric, Vec3);
|
||||
impl FlatCoordinateSystem<Vec3> for Scaled {
|
||||
fn flat_to_global(&self, pos: Vec3) -> Vec3 {
|
||||
(pos - self.1) / self.0.scale
|
||||
}
|
||||
fn global_to_flat(&self, pos: Vec3) -> Vec3 {
|
||||
pos * self.0.scale + self.1
|
||||
}
|
||||
}
|
||||
impl MetricCS for Scaled {
|
||||
fn global_metric(&self) -> &impl Metric {
|
||||
&self.0
|
||||
}
|
||||
}
|
||||
let cs = Scaled(
|
||||
samples::ScaledMetric {
|
||||
scale: vec3(3., 4., 5.),
|
||||
},
|
||||
vec3(2., 3., 7.),
|
||||
);
|
||||
assert_eq!(cs.global_to_flat(vec3(7., 3., 1.)), vec3(23., 15., 12.));
|
||||
assert_eq!(cs.flat_to_global(vec3(8., 7., 17.)), vec3(2., 1., 2.));
|
||||
assert_eq!(
|
||||
cs.global_to_flat(Ray {
|
||||
pos: vec3(7., 3., 0.),
|
||||
dir: vec3(3., 2., 0.)
|
||||
}),
|
||||
Ray {
|
||||
pos: vec3(23., 15., 7.),
|
||||
dir: vec3(9., 8., 0.)
|
||||
}
|
||||
);
|
||||
assert_eq!(
|
||||
cs.flat_to_global(Ray {
|
||||
pos: vec3(23., 15., 7.),
|
||||
dir: vec3(9., 8., 0.)
|
||||
}),
|
||||
Ray {
|
||||
pos: vec3(7., 3., 0.),
|
||||
dir: vec3(3., 2., 0.)
|
||||
}
|
||||
);
|
||||
assert_eq!(
|
||||
cs.global_to_flat(Location {
|
||||
pos: vec3(2., 1., 0.),
|
||||
rot: mat3(vec3(0., 1., 0.), vec3(-1., 0., 0.), vec3(0., 0., 1.))
|
||||
}),
|
||||
Location {
|
||||
pos: vec3(8., 7., 7.),
|
||||
rot: mat3(vec3(0., 4., 0.), vec3(-3., 0., 0.), vec3(0., 0., 5.))
|
||||
}
|
||||
);
|
||||
assert_eq!(
|
||||
cs.flat_to_global(Location {
|
||||
pos: vec3(2., 1., 7.),
|
||||
rot: mat3(vec3(0., 1., 0.), vec3(-1., 0., 0.), vec3(0., 0., 1.))
|
||||
}),
|
||||
Location {
|
||||
pos: vec3(0., -0.5, 0.),
|
||||
rot: mat3(vec3(0., 0.25, 0.), vec3(-1. / 3., 0., 0.), vec3(0., 0., 0.2))
|
||||
}
|
||||
);
|
||||
}
|
||||
|
||||
fn test_flat_region(region: &impl FlatRegion, range_global: (Vec3, Vec3), range_flat: (Vec3, Vec3)) {
|
||||
#[allow(non_upper_case_globals)]
|
||||
const ε: f32 = 1e-3;
|
||||
macro_rules! assert_eq_at {
|
||||
($at: expr, $left: expr, $right: expr) => {
|
||||
let at = $at;
|
||||
let left = $left;
|
||||
let right = $right;
|
||||
assert!(
|
||||
left.abs_diff_eq(right, ε),
|
||||
"Assertion failed at {at}:\n left: {left} = {}\n right: {right} = {}",
|
||||
stringify!($left),
|
||||
stringify!($right)
|
||||
);
|
||||
};
|
||||
}
|
||||
fn check_range(name_a: &str, a: Vec3, range_a: (Vec3, Vec3), name_b: &str, b: Vec3, range_b: (Vec3, Vec3)) {
|
||||
assert!(b.cmpge(range_b.0 - ε).all() && b.cmple(range_b.1 + ε).all(), "Assertion failed:\nAt {name_a}: {a}, from range: {range_a:?}\nGot {name_b}: {b}, which is out of range {range_b:?}");
|
||||
// TODO sort out when to check these conditions:
|
||||
if a.x.abs_diff_eq(&range_a.0.x, ε) {
|
||||
assert_abs_diff_eq!(b.x, range_b.0.x, epsilon = ε);
|
||||
}
|
||||
if a.y.abs_diff_eq(&range_a.0.y, ε) {
|
||||
assert_abs_diff_eq!(b.y, range_b.0.y, epsilon = ε);
|
||||
}
|
||||
if a.x.abs_diff_eq(&range_a.1.x, ε) {
|
||||
assert_abs_diff_eq!(b.x, range_b.1.x, epsilon = ε);
|
||||
}
|
||||
if a.y.abs_diff_eq(&range_a.1.y, ε) {
|
||||
assert_abs_diff_eq!(b.y, range_b.1.y, epsilon = ε);
|
||||
}
|
||||
}
|
||||
for x in linspace(range_global.0.x, range_global.1.x, 20) {
|
||||
for y in linspace(range_global.0.y, range_global.1.y, 20) {
|
||||
for z in linspace(range_global.0.z, range_global.1.z, 20) {
|
||||
let pos_global = vec3(x, y, z);
|
||||
let pos_flat = region.global_to_flat(pos_global);
|
||||
check_range("global", pos_global, range_global, "flat", pos_flat, range_flat);
|
||||
assert_eq_at!(
|
||||
pos_global,
|
||||
region
|
||||
.global_to_flat(Location {
|
||||
pos: pos_global,
|
||||
rot: Mat3::IDENTITY
|
||||
})
|
||||
.pos,
|
||||
pos_flat
|
||||
);
|
||||
assert_eq_at!(pos_global, region.flat_to_global(pos_flat), pos_global);
|
||||
assert_eq_at!(
|
||||
pos_global,
|
||||
region
|
||||
.flat_to_global(region.global_to_flat(Location {
|
||||
pos: pos_global,
|
||||
rot: Mat3::IDENTITY
|
||||
}))
|
||||
.rot,
|
||||
Mat3::IDENTITY
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
for x in linspace(range_flat.0.x, range_flat.1.x, 20) {
|
||||
for y in linspace(range_flat.0.y, range_flat.1.y, 20) {
|
||||
for z in linspace(range_flat.0.z, range_flat.1.z, 20) {
|
||||
let pos_flat = vec3(x, y, z);
|
||||
let pos_global = region.flat_to_global(pos_flat);
|
||||
check_range("flat", pos_flat, range_flat, "global", pos_global, range_global);
|
||||
assert_eq_at!(
|
||||
pos_flat,
|
||||
region
|
||||
.flat_to_global(Location {
|
||||
pos: pos_flat,
|
||||
rot: Mat3::IDENTITY
|
||||
})
|
||||
.pos,
|
||||
pos_global
|
||||
);
|
||||
assert_eq_at!(pos_flat, region.global_to_flat(pos_global), pos_flat);
|
||||
assert_eq_at!(
|
||||
pos_flat,
|
||||
region
|
||||
.global_to_flat(region.flat_to_global(Location {
|
||||
pos: pos_global,
|
||||
rot: Mat3::IDENTITY
|
||||
}))
|
||||
.rot,
|
||||
Mat3::IDENTITY
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_mapper_inner() {
|
||||
let mapper = InnerCS(Tube {
|
||||
inner_radius: 30.0,
|
||||
outer_radius: 50.0,
|
||||
internal_halflength: 100.0,
|
||||
external_halflength: 300.0,
|
||||
});
|
||||
test_flat_region(
|
||||
&mapper,
|
||||
(vec3(-30.0, -300.0, -30.0), vec3(30.0, 300.0, 30.0)),
|
||||
(vec3(-30.0, -100.0, -30.0), vec3(30.0, 100.0, 30.0)),
|
||||
);
|
||||
test_flat_region(
|
||||
&mapper,
|
||||
(vec3(-60.0, -400.0, -60.0), vec3(60.0, -300.0, 60.0)),
|
||||
(vec3(-60.0, -200.0, -60.0), vec3(60.0, -100.0, 60.0)),
|
||||
);
|
||||
test_flat_region(
|
||||
&mapper,
|
||||
(vec3(-60.0, 300.0, -60.0), vec3(60.0, 400.0, 60.0)),
|
||||
(vec3(-60.0, 100.0, -60.0), vec3(60.0, 200.0, 60.0)),
|
||||
);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_mapper_outer() {
|
||||
let mapper = OuterCS(Tube {
|
||||
inner_radius: 30.0,
|
||||
outer_radius: 50.0,
|
||||
internal_halflength: 100.0,
|
||||
external_halflength: 300.0,
|
||||
});
|
||||
// TODO replace 200.20016 with something sane
|
||||
// NOTE it can’t test −30..30 as that area intersects the boundary
|
||||
test_flat_region(
|
||||
&mapper,
|
||||
(vec3(-20.0, -300.0, -20.0), vec3(20.0, -1.0, 20.0)),
|
||||
(vec3(-20.0, -300.0, -20.0), vec3(20.0, -200.20016, 20.0)),
|
||||
);
|
||||
test_flat_region(
|
||||
&mapper,
|
||||
(vec3(-20.0, 1.0, -20.0), vec3(20.0, 300.0, 20.0)),
|
||||
(vec3(-20.0, 200.20016, -20.0), vec3(20.0, 300.0, 20.0)),
|
||||
);
|
||||
test_flat_region(
|
||||
&mapper,
|
||||
(vec3(-60.0, -400.0, -60.0), vec3(60.0, -300.0, 60.0)),
|
||||
(vec3(-60.0, -400.0, -60.0), vec3(60.0, -300.0, 60.0)),
|
||||
);
|
||||
test_flat_region(
|
||||
&mapper,
|
||||
(vec3(-60.0, 300.0, -60.0), vec3(60.0, 400.0, 60.0)),
|
||||
(vec3(-60.0, 300.0, -60.0), vec3(60.0, 400.0, 60.0)),
|
||||
);
|
||||
// straight
|
||||
for x in linspace(-60., 60., 20) {
|
||||
for y in linspace(-320., 320., 20) {
|
||||
for z in linspace(-60., 60., 20) {
|
||||
assert_eq!(
|
||||
mapper
|
||||
.global_to_flat(Location {
|
||||
pos: vec3(x, y, z),
|
||||
rot: Mat3::IDENTITY
|
||||
})
|
||||
.pos
|
||||
.x,
|
||||
x
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
// symmetrical
|
||||
for x in linspace(0., 60., 20) {
|
||||
for y in linspace(0., 320., 20) {
|
||||
for z in linspace(0., 60., 20) {
|
||||
let pp = mapper
|
||||
.global_to_flat(Location {
|
||||
pos: vec3(x, y, z),
|
||||
rot: Mat3::IDENTITY,
|
||||
})
|
||||
.pos;
|
||||
let np = mapper
|
||||
.global_to_flat(Location {
|
||||
pos: vec3(-x, y, z),
|
||||
rot: Mat3::IDENTITY,
|
||||
})
|
||||
.pos;
|
||||
let pn = mapper
|
||||
.global_to_flat(Location {
|
||||
pos: vec3(x, -y, z),
|
||||
rot: Mat3::IDENTITY,
|
||||
})
|
||||
.pos;
|
||||
let nn = mapper
|
||||
.global_to_flat(Location {
|
||||
pos: vec3(-x, -y, z),
|
||||
rot: Mat3::IDENTITY,
|
||||
})
|
||||
.pos;
|
||||
assert_eq!(np, vec3(-pp.x, pp.y, pp.z));
|
||||
assert_eq!(pn, vec3(pp.x, -pp.y, pp.z));
|
||||
assert_eq!(nn, vec3(-pp.x, -pp.y, pp.z));
|
||||
}
|
||||
}
|
||||
}
|
||||
// clean boundary
|
||||
for x in linspace(50., 60., 20) {
|
||||
for y in linspace(0., 320., 20) {
|
||||
for z in linspace(50., 60., 20) {
|
||||
assert_eq!(
|
||||
mapper
|
||||
.global_to_flat(Location {
|
||||
pos: vec3(x, y, z),
|
||||
rot: Mat3::IDENTITY
|
||||
})
|
||||
.pos
|
||||
.y,
|
||||
y
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
for x in linspace(0., 60., 20) {
|
||||
for y in linspace(300., 320., 20) {
|
||||
for z in linspace(0., 60., 20) {
|
||||
assert_eq!(
|
||||
mapper
|
||||
.global_to_flat(Location {
|
||||
pos: vec3(x, y, z),
|
||||
rot: Mat3::IDENTITY
|
||||
})
|
||||
.pos
|
||||
.y,
|
||||
y
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
// accelerating
|
||||
for x in linspace(-21., 21., 20) {
|
||||
for y in linspace(1., 299., 20) {
|
||||
for z in linspace(-21., 21., 20) {
|
||||
let v = mapper
|
||||
.global_to_flat(Location {
|
||||
pos: vec3(x, y, z),
|
||||
rot: Mat3::IDENTITY,
|
||||
})
|
||||
.pos
|
||||
.y;
|
||||
assert!(v > 200.0);
|
||||
assert!(v > y);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -1,239 +0,0 @@
|
|||
use glam::{f32, vec3, Mat3, Vec3};
|
||||
|
||||
use crate::fns::{self, Limiter};
|
||||
use crate::mathx::Decomp3;
|
||||
use crate::riemann::{Metric, Tens3};
|
||||
|
||||
#[derive(Copy, Clone, Debug)]
|
||||
pub struct Tube {
|
||||
pub outer_radius: f32,
|
||||
pub inner_radius: f32,
|
||||
pub external_halflength: f32,
|
||||
pub internal_halflength: f32,
|
||||
}
|
||||
|
||||
impl Tube {
|
||||
fn fx(&self) -> impl Limiter {
|
||||
fns::SmootherstepLimiter {
|
||||
min: self.inner_radius,
|
||||
max: self.outer_radius,
|
||||
}
|
||||
}
|
||||
fn fy(&self) -> fns::QuadraticAccelerator {
|
||||
fns::QuadraticAccelerator {
|
||||
internal: self.internal_halflength,
|
||||
external: self.external_halflength,
|
||||
}
|
||||
}
|
||||
|
||||
pub fn y(&self, v: f32) -> f32 {
|
||||
self.fy().x(v)
|
||||
}
|
||||
pub fn v(&self, y: f32) -> f32 {
|
||||
self.fy().u(y)
|
||||
}
|
||||
pub fn dy(&self, v: f32) -> f32 {
|
||||
self.fy().dx(v)
|
||||
}
|
||||
pub fn dv(&self, y: f32) -> f32 {
|
||||
self.fy().du(y)
|
||||
}
|
||||
}
|
||||
|
||||
impl Metric for Tube {
|
||||
fn sqrt_at(&self, pos: Vec3) -> Decomp3 {
|
||||
let r = f32::hypot(pos.x, pos.z);
|
||||
let sr = self.fx().value(r);
|
||||
let sy = self.fy().du(pos.y);
|
||||
let s = sr + sy - sr * sy;
|
||||
assert!(sr.is_finite());
|
||||
assert!(sy.is_finite());
|
||||
assert!(sy > 0.0);
|
||||
Decomp3 {
|
||||
ortho: Mat3::IDENTITY,
|
||||
diag: vec3(1.0, s, 1.0),
|
||||
}
|
||||
}
|
||||
|
||||
fn part_derivs_at(&self, pos: Vec3) -> Tens3 {
|
||||
let r = f32::hypot(pos.x, pos.z);
|
||||
let sr = self.fx().value(r);
|
||||
let sy = self.fy().du(pos.y);
|
||||
let s = sr + sy - sr * sy;
|
||||
let dsr_dr = self.fx().derivative(r);
|
||||
let dsy_dy = self.fy().d2u(pos.y);
|
||||
let ds2_dr = 2.0 * s * (1.0 - sy) * dsr_dr;
|
||||
let ds2_dy = 2.0 * s * (1.0 - sr) * dsy_dy;
|
||||
let (ds2_dx, ds2_dz) = if ds2_dr.abs() < 1e-5 {
|
||||
(0., 0.)
|
||||
} else {
|
||||
(ds2_dr * pos.x / r, ds2_dr * pos.z / r)
|
||||
};
|
||||
[
|
||||
Mat3::from_cols_array(&[0., 0., 0., 0., ds2_dx, 0., 0., 0., 0.]),
|
||||
Mat3::from_cols_array(&[0., 0., 0., 0., ds2_dy, 0., 0., 0., 0.]),
|
||||
Mat3::from_cols_array(&[0., 0., 0., 0., ds2_dz, 0., 0., 0., 0.]),
|
||||
]
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod test {
|
||||
use approx::assert_abs_diff_eq;
|
||||
use glam::{vec3, Vec3};
|
||||
use itertools_num::linspace;
|
||||
|
||||
use crate::mathx::Decomp3;
|
||||
use crate::riemann::Metric;
|
||||
use crate::tube::Space;
|
||||
use crate::types::Ray;
|
||||
|
||||
use super::Tube;
|
||||
|
||||
#[test]
|
||||
fn test_tube_metric_derivs() {
|
||||
struct Approx(Tube);
|
||||
impl Metric for Approx {
|
||||
fn sqrt_at(&self, pos: Vec3) -> Decomp3 {
|
||||
self.0.sqrt_at(pos)
|
||||
}
|
||||
}
|
||||
let testee = Tube {
|
||||
inner_radius: 30.0,
|
||||
outer_radius: 50.0,
|
||||
internal_halflength: 100.0,
|
||||
external_halflength: 300.0,
|
||||
};
|
||||
let approx = Approx(testee);
|
||||
let epsilon = 1.0e-3;
|
||||
let margin = 1.0 / 16.0;
|
||||
let mul = 1.0 + margin;
|
||||
for x in itertools_num::linspace(-mul * testee.outer_radius, mul * testee.outer_radius, 20) {
|
||||
for y in itertools_num::linspace(-mul * testee.external_halflength, mul * testee.external_halflength, 20) {
|
||||
for z in itertools_num::linspace(-mul * testee.outer_radius, mul * testee.outer_radius, 20) {
|
||||
let pos = vec3(x, y, z);
|
||||
let computed = testee.part_derivs_at(pos);
|
||||
let reference = approx.part_derivs_at(pos);
|
||||
let eq = (0..3).all(|coord| computed[coord].abs_diff_eq(reference[coord], epsilon));
|
||||
assert!(
|
||||
eq,
|
||||
"Bad derivative computation at {pos}:\n explicit: {computed:?}\n numerical: {reference:?}\n"
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_accelerator() {
|
||||
let space = Space {
|
||||
tube: Tube {
|
||||
inner_radius: 30.0,
|
||||
outer_radius: 50.0,
|
||||
internal_halflength: 100.0,
|
||||
external_halflength: 300.0,
|
||||
},
|
||||
objs: vec![],
|
||||
};
|
||||
let ε = 1e-3;
|
||||
let off = 10.0;
|
||||
let steps = 1024;
|
||||
for ax in [-30.0 + ε, -25.0, -3.0, 17.0, 30.0 - ε] {
|
||||
for bx in [0.0, ε, 1.0, 7.0, 30.0 - ε] {
|
||||
let a = vec3(ax, -(space.tube.external_halflength + off), 0.);
|
||||
let b = vec3(bx, space.tube.external_halflength + off, 0.);
|
||||
let δ = vec3(bx - ax, 2.0 * (space.tube.internal_halflength + off), 0.);
|
||||
let dir = δ / (steps as f32);
|
||||
let traced = space.trace_iter(Ray { pos: a, dir }).nth(steps).unwrap();
|
||||
assert_abs_diff_eq!(traced.pos.x, b.x, epsilon = 1.0e-2);
|
||||
assert_abs_diff_eq!(traced.pos.y, b.y, epsilon = 1.0e1);
|
||||
assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon = 1.0e-3);
|
||||
assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon = 1.0e-2);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
#[ignore]
|
||||
fn test_accelerator_slow() {
|
||||
let space = Space {
|
||||
tube: Tube {
|
||||
inner_radius: 30.0,
|
||||
outer_radius: 50.0,
|
||||
internal_halflength: 100.0,
|
||||
external_halflength: 300.0,
|
||||
},
|
||||
objs: vec![],
|
||||
};
|
||||
let ε = 1e-3;
|
||||
let off = 10.0;
|
||||
let steps = 4096;
|
||||
for ax in linspace(-space.tube.inner_radius + ε, space.tube.inner_radius - ε, 20) {
|
||||
for bx in linspace(-space.tube.inner_radius + ε, space.tube.inner_radius - ε, 20) {
|
||||
let a = vec3(ax, -(space.tube.external_halflength + off), 0.);
|
||||
let b = vec3(bx, space.tube.external_halflength + off, 0.);
|
||||
let δ = vec3(bx - ax, 2.0 * (space.tube.internal_halflength + off), 0.);
|
||||
let dir = δ / (steps as f32);
|
||||
let traced = space.trace_iter(Ray { pos: a, dir }).nth(steps).unwrap();
|
||||
assert_abs_diff_eq!(traced.pos.x, b.x, epsilon = 1.0e-2);
|
||||
assert_abs_diff_eq!(traced.pos.y, b.y, epsilon = 1.0e0);
|
||||
assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon = 1.0e-3);
|
||||
assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon = 1.0e-3);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_accelerator_inner_edge() {
|
||||
let space = Space {
|
||||
tube: Tube {
|
||||
inner_radius: 30.0,
|
||||
outer_radius: 50.0,
|
||||
internal_halflength: 100.0,
|
||||
external_halflength: 300.0,
|
||||
},
|
||||
objs: vec![],
|
||||
};
|
||||
let ε = 1e-3;
|
||||
let off = 10.0;
|
||||
let steps = 10000;
|
||||
for x in [space.tube.inner_radius - ε, space.tube.inner_radius + ε] {
|
||||
let a = vec3(x, -(space.tube.external_halflength + off), 0.);
|
||||
let b = vec3(x, space.tube.external_halflength + off, 0.);
|
||||
let δ = vec3(0.0, 2.0 * (space.tube.internal_halflength + off), 0.);
|
||||
let dir = δ / (steps as f32);
|
||||
let traced = space.trace_iter(Ray { pos: a, dir }).nth(steps).unwrap();
|
||||
assert_abs_diff_eq!(traced.pos.x, b.x, epsilon = 1.0e-1);
|
||||
assert_abs_diff_eq!(traced.pos.y, b.y, epsilon = 1.0e0);
|
||||
assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon = 1.0e-2);
|
||||
assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon = 1.0e-2);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_accelerator_outer_edge() {
|
||||
let space = Space {
|
||||
tube: Tube {
|
||||
inner_radius: 30.0,
|
||||
outer_radius: 50.0,
|
||||
internal_halflength: 100.0,
|
||||
external_halflength: 300.0,
|
||||
},
|
||||
objs: vec![],
|
||||
};
|
||||
let ε = 1e-3;
|
||||
let off = 10.0;
|
||||
let steps = 4096;
|
||||
for x in [space.tube.outer_radius + ε, space.tube.outer_radius - ε] {
|
||||
let a = vec3(x, -(space.tube.external_halflength + off), 0.);
|
||||
let b = vec3(x, space.tube.external_halflength + off, 0.);
|
||||
let δ = vec3(0.0, 2.0 * (space.tube.external_halflength + off), 0.);
|
||||
let dir = δ / (steps as f32);
|
||||
let traced = space.trace_iter(Ray { pos: a, dir }).nth(steps).unwrap();
|
||||
assert_abs_diff_eq!(traced.pos.x, b.x, epsilon = 2.0e0);
|
||||
assert_abs_diff_eq!(traced.pos.y, b.y, epsilon = 1.0e0);
|
||||
assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon = 1.0e-2);
|
||||
assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon = 1.0e-2);
|
||||
}
|
||||
}
|
||||
}
|
||||
228
src/tube/mod.rs
228
src/tube/mod.rs
|
|
@ -1,228 +0,0 @@
|
|||
use glam::{f32, Mat3, Vec3};
|
||||
|
||||
use crate::ifaces::{DebugTraceable, RayPath, Traceable};
|
||||
use coords::{FlatCoordinateSystem, InnerCS, OuterCS};
|
||||
use metric::Tube;
|
||||
use Subspace::{Boundary, Inner, Outer};
|
||||
|
||||
use crate::riemann::Metric;
|
||||
use crate::tube::coords::FlatRegion;
|
||||
use crate::types::{FlatTraceResult, Hit, Location, Object, Ray};
|
||||
use crate::{riemann, DT};
|
||||
|
||||
mod coords;
|
||||
pub mod metric;
|
||||
|
||||
pub struct Space {
|
||||
pub tube: Tube,
|
||||
pub objs: Vec<Object>,
|
||||
}
|
||||
|
||||
#[derive(PartialEq, Eq, Debug)]
|
||||
pub enum Subspace {
|
||||
Outer,
|
||||
Boundary,
|
||||
Inner,
|
||||
}
|
||||
|
||||
impl Space {
|
||||
fn which_subspace(&self, pt: Vec3) -> Subspace {
|
||||
if pt.y.abs() > self.tube.external_halflength {
|
||||
return Outer;
|
||||
}
|
||||
let r = f32::hypot(pt.x, pt.z);
|
||||
if r > self.tube.outer_radius {
|
||||
Outer
|
||||
} else if r > self.tube.inner_radius {
|
||||
Boundary
|
||||
} else {
|
||||
Inner
|
||||
}
|
||||
}
|
||||
|
||||
/// Выполняет один шаг трассировки. Работает в любой части пространства, но вне Boundary доступны более эффективные методы.
|
||||
/// ray задаётся в основной СК.
|
||||
pub fn trace_step(&self, ray: Ray) -> Ray {
|
||||
let a = -riemann::contract2(riemann::krist(&self.tube, ray.pos), ray.dir);
|
||||
let v = ray.dir + a;
|
||||
let p = ray.pos + v;
|
||||
Ray { pos: p, dir: v }
|
||||
}
|
||||
|
||||
/// Выполняет один шаг перемещения. Работает в любой части пространства.
|
||||
/// off задаётся в локальной СК. Рекомендуется считать небольшими шагами.
|
||||
pub fn move_step(&self, loc: Location, off: Vec3) -> Location {
|
||||
let corr = Mat3::IDENTITY - riemann::contract(riemann::krist(&self.tube, loc.pos), loc.rot * off);
|
||||
let p = loc.pos + corr * loc.rot * off;
|
||||
Location {
|
||||
pos: p,
|
||||
rot: corr * loc.rot,
|
||||
}
|
||||
}
|
||||
|
||||
pub fn trace_iter(&self, ray: Ray) -> impl Iterator<Item = Ray> + '_ {
|
||||
std::iter::successors(Some(ray), |&ray| Some(self.trace_step(ray)))
|
||||
}
|
||||
|
||||
fn trace_inner(&self, ray: Ray) -> FlatTraceResult {
|
||||
assert_eq!(self.which_subspace(ray.pos), Inner);
|
||||
self.trace_flat(InnerCS(self.tube), ray)
|
||||
}
|
||||
|
||||
fn trace_outer(&self, ray: Ray) -> FlatTraceResult {
|
||||
assert_eq!(self.which_subspace(ray.pos), Outer);
|
||||
self.trace_flat(OuterCS(self.tube), ray)
|
||||
}
|
||||
|
||||
fn obj_hitter(&self, pos: Vec3) -> Option<fn(&Self, ray: Ray) -> FlatTraceResult> {
|
||||
match self.which_subspace(pos) {
|
||||
Inner => Some(Self::trace_inner),
|
||||
Outer => Some(Self::trace_outer),
|
||||
Boundary => None,
|
||||
}
|
||||
}
|
||||
|
||||
fn trace_flat(&self, cs: impl FlatRegion, ray: Ray) -> FlatTraceResult {
|
||||
let ray = cs.global_to_flat(ray);
|
||||
let dist = cs.distance_to_boundary(ray);
|
||||
let objs = self.list_objects(|loc| cs.global_to_flat(loc));
|
||||
FlatTraceResult {
|
||||
end: dist.map(|dist| cs.flat_to_global(ray.forward(dist))),
|
||||
objects: Self::hit_objects(objs.as_slice(), ray, dist, |pos| cs.flat_to_global(pos)),
|
||||
}
|
||||
}
|
||||
|
||||
fn list_objects(&self, tfm: impl Fn(Location) -> Location) -> Vec<Object> {
|
||||
self.objs
|
||||
.iter()
|
||||
.map(|&Object { id, loc, r }| Object { id, loc: tfm(loc), r })
|
||||
.collect()
|
||||
}
|
||||
|
||||
fn hit_objects(objs: &[Object], ray: Ray, limit: Option<f32>, globalize: impl Fn(Vec3) -> Vec3) -> Vec<Hit> {
|
||||
let limit = limit.unwrap_or(f32::INFINITY);
|
||||
objs.iter()
|
||||
.filter_map(|obj| {
|
||||
let rel = ray.pos - obj.loc.pos;
|
||||
let diff = rel.dot(ray.dir).powi(2) - ray.dir.length_squared() * (rel.length_squared() - obj.r.powi(2));
|
||||
if diff > 0.0 {
|
||||
let t = (-rel.dot(ray.dir) - diff.sqrt()) / ray.dir.length_squared();
|
||||
Some((obj, t))
|
||||
} else {
|
||||
None
|
||||
}
|
||||
})
|
||||
.filter(|&(_, t)| t >= 0.0 && t < limit)
|
||||
.map(|(obj, t)| {
|
||||
let pos = ray.forward(t).pos;
|
||||
let rel = obj.loc.rot.inverse()
|
||||
* Ray {
|
||||
pos: pos - obj.loc.pos,
|
||||
dir: ray.dir,
|
||||
};
|
||||
Hit {
|
||||
id: obj.id,
|
||||
distance: t,
|
||||
pos: globalize(pos),
|
||||
rel,
|
||||
}
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
|
||||
pub fn line(&self, a: Vec3, b: Vec3, step: f32) -> Vec<Ray> {
|
||||
match self.which_subspace(a) {
|
||||
Outer => vec![Ray {
|
||||
pos: b,
|
||||
dir: (b - a).normalize(),
|
||||
}],
|
||||
Inner => {
|
||||
let cs = InnerCS(self.tube);
|
||||
let n = ((b - a).length() / step) as usize + 1;
|
||||
let a = cs.global_to_flat(a);
|
||||
let b = cs.global_to_flat(b);
|
||||
let dir = (b - a).normalize();
|
||||
(1..=n)
|
||||
.map(|k| {
|
||||
cs.flat_to_global(Ray {
|
||||
pos: a.lerp(b, k as f32 / n as f32),
|
||||
dir,
|
||||
})
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
Boundary => panic!("Can't draw a line here!"),
|
||||
}
|
||||
}
|
||||
|
||||
fn camera_ray_to_abs(&self, camera: Location, ray: Ray) -> Ray {
|
||||
let pos = camera.pos;
|
||||
let dir = camera.rot * ray.dir;
|
||||
// TODO account for ray.pos
|
||||
let dir = DT * self.tube.normalize_vec_at(pos, dir);
|
||||
Ray { pos, dir }
|
||||
}
|
||||
}
|
||||
|
||||
/// Like [`std::iter::successors`] but with an upper limit on iteration count.
|
||||
///
|
||||
/// # Panics
|
||||
///
|
||||
/// Panics if the sequence doesn’t terminate in `max_iters` calls of `succ`.
|
||||
fn iterate_with_limit<T>(max_iters: usize, init: T, mut succ: impl FnMut(T) -> Option<T>) {
|
||||
let mut state = init;
|
||||
for _ in 0..max_iters {
|
||||
match succ(state) {
|
||||
Some(next) => state = next,
|
||||
None => return,
|
||||
}
|
||||
}
|
||||
panic!("iteration limit exceeded");
|
||||
}
|
||||
|
||||
impl Traceable for Space {
|
||||
fn trace(&self, camera: Location, ray: Ray) -> Vec<Hit> {
|
||||
let ray = self.camera_ray_to_abs(camera, ray);
|
||||
let mut hits = vec![];
|
||||
iterate_with_limit(100, ray, |ray| {
|
||||
let ret = self
|
||||
.trace_iter(ray)
|
||||
.skip(1)
|
||||
.find_map(|ray| self.obj_hitter(ray.pos).map(|hitter| hitter(self, ray)))
|
||||
.expect("Space::trace_iter does not terminate");
|
||||
hits.extend(ret.objects); // TODO fix distance
|
||||
ret.end
|
||||
});
|
||||
hits
|
||||
}
|
||||
}
|
||||
|
||||
impl DebugTraceable for Space {
|
||||
fn trace_dbg(&self, camera: Location, ray: Ray) -> (Vec<Hit>, RayPath) {
|
||||
let mut points = vec![];
|
||||
let mut hits = vec![];
|
||||
let mut ray = self.camera_ray_to_abs(camera, ray);
|
||||
|
||||
let trace_to_flat = |points: &mut Vec<Ray>, ray| {
|
||||
for ray in self.trace_iter(ray).skip(1) {
|
||||
points.push(ray);
|
||||
if let Some(hitter) = self.obj_hitter(ray.pos) {
|
||||
return (ray, hitter(self, ray));
|
||||
}
|
||||
}
|
||||
unreachable!("Space::trace_iter terminated!")
|
||||
};
|
||||
|
||||
points.push(ray);
|
||||
for _ in 0..100 {
|
||||
let (ray_into_flat, ret) = trace_to_flat(&mut points, ray);
|
||||
hits.extend(ret.objects); // TODO fix distance
|
||||
let Some(ray_outta_flat) = ret.end else {
|
||||
return (hits, RayPath { points });
|
||||
};
|
||||
points.extend(self.line(ray_into_flat.pos, ray_outta_flat.pos, 1.0));
|
||||
ray = ray_outta_flat;
|
||||
}
|
||||
panic!("tracing didn't terminate");
|
||||
}
|
||||
}
|
||||
60
src/types.rs
60
src/types.rs
|
|
@ -1,60 +0,0 @@
|
|||
use glam::{f32, i32, Mat3, Vec3};
|
||||
|
||||
#[derive(Copy, Clone, Debug, PartialEq)]
|
||||
pub struct Ray {
|
||||
pub pos: Vec3,
|
||||
pub dir: Vec3,
|
||||
}
|
||||
|
||||
pub fn ray(pos: Vec3, dir: Vec3) -> Ray {
|
||||
Ray { pos, dir }
|
||||
}
|
||||
|
||||
impl Ray {
|
||||
pub fn forward(&self, dist: f32) -> Ray {
|
||||
Ray {
|
||||
pos: self.pos + self.dir * dist,
|
||||
dir: self.dir,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl std::ops::Mul<Ray> for Mat3 {
|
||||
type Output = Ray;
|
||||
|
||||
fn mul(self, rhs: Ray) -> Self::Output {
|
||||
Ray {
|
||||
pos: self * rhs.pos,
|
||||
dir: self * rhs.dir,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Copy, Clone, Debug, PartialEq)]
|
||||
pub struct Location {
|
||||
/// Положение в основной СК
|
||||
pub pos: Vec3,
|
||||
/// Преобразование вектора из локальной ортонормированной в основную СК
|
||||
pub rot: Mat3,
|
||||
}
|
||||
|
||||
#[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: Vec3, // положение в основной СК
|
||||
pub rel: Ray, // в локальной ортонормированной СК объекта
|
||||
}
|
||||
|
||||
#[derive(Clone, Debug)]
|
||||
pub struct FlatTraceResult {
|
||||
pub end: Option<Ray>,
|
||||
pub objects: Vec<Hit>,
|
||||
}
|
||||
68
src/utils.rs
68
src/utils.rs
|
|
@ -1,68 +0,0 @@
|
|||
//! Utility functions to work with metrics.
|
||||
|
||||
use crate::{
|
||||
mathx::MatExt as _,
|
||||
riemann::{trace_iter, Metric},
|
||||
types::Location,
|
||||
};
|
||||
use glam::{Mat3, Vec3};
|
||||
|
||||
pub fn rel_to_abs(space: &impl Metric, base: &Location, rel: Vec3, steps: usize) -> Vec3 {
|
||||
let c = 1.0 / (steps as f32);
|
||||
trace_iter(space, base.pos, base.rot * rel, c * rel.length())
|
||||
.nth(steps - 1)
|
||||
.unwrap()
|
||||
}
|
||||
|
||||
/// Converts a position and a rotation to a [Location]. Only the X direction is preserved from `rot` to ensure the resulting Location describes an orthonormal coordinate system.
|
||||
pub fn put_object(space: &impl Metric, pos: Vec3, rot: Mat3) -> Location {
|
||||
let metric_sqrt = space.sqrt_at(pos);
|
||||
let metric_inv_sqrt = space.sqrt_at(pos).inverse();
|
||||
let rot = metric_inv_sqrt * (metric_sqrt * rot).orthonormalize();
|
||||
Location { pos, rot }
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use crate::riemann::samples;
|
||||
use glam::{mat3, vec3};
|
||||
|
||||
#[test]
|
||||
fn test_put_object() {
|
||||
use approx::assert_abs_diff_eq;
|
||||
|
||||
let ε = 1e-5;
|
||||
let m = samples::ScaledMetric {
|
||||
scale: vec3(3., 4., 5.),
|
||||
};
|
||||
|
||||
let loc = put_object(
|
||||
&m,
|
||||
vec3(1., 2., 0.),
|
||||
mat3(vec3(1., 0., 0.), vec3(0., 1., 0.), vec3(0., 0., 1.)),
|
||||
);
|
||||
assert_eq!(loc.pos, vec3(1., 2., 0.));
|
||||
assert_abs_diff_eq!(loc.rot * vec3(1., 0., 0.), vec3(1. / 3., 0., 0.), epsilon = ε);
|
||||
assert_abs_diff_eq!(loc.rot * vec3(0., 1., 0.), vec3(0., 1. / 4., 0.), epsilon = ε);
|
||||
|
||||
let loc = put_object(
|
||||
&m,
|
||||
vec3(1., 2., 0.),
|
||||
mat3(vec3(0., 1., 0.), vec3(-1., 0., 0.), vec3(0., 0., 1.)),
|
||||
);
|
||||
assert_eq!(loc.pos, vec3(1., 2., 0.));
|
||||
assert_abs_diff_eq!(loc.rot * vec3(1., 0., 0.), vec3(0., 1. / 4., 0.), epsilon = ε);
|
||||
assert_abs_diff_eq!(loc.rot * vec3(0., 1., 0.), vec3(-1. / 3., 0., 0.), epsilon = ε);
|
||||
|
||||
let c = 0.5 * std::f32::consts::SQRT_2;
|
||||
let loc = put_object(
|
||||
&m,
|
||||
vec3(1., 2., 0.),
|
||||
mat3(vec3(c, c, 0.), vec3(-c, c, 0.), vec3(0., 0., 1.)),
|
||||
);
|
||||
assert_eq!(loc.pos, vec3(1., 2., 0.));
|
||||
assert_abs_diff_eq!(loc.rot * vec3(1., 0., 0.), vec3(1. / 5., 1. / 5., 0.), epsilon = ε);
|
||||
assert_abs_diff_eq!(loc.rot * vec3(0., 1., 0.), vec3(-4. / 15., 3. / 20., 0.), epsilon = ε);
|
||||
}
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user