matthew
/
vanrijn
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Compare commits

base: matthew:master
Branches Tags
matthew:master
matthew:glass-bunny
matthew:github-build-action
..
compare: matthew:github-build-action
Branches Tags
matthew:master
matthew:glass-bunny
matthew:github-build-action

5 Commits
master ... github-bui

Author SHA1 Message Date
Matthew Gordon 8c4a0002fe
WIP 2019-11-28 10:49:21 -05:00
Matthew Gordon a700a28d9c
Github build action WIP 2019-11-28 10:39:04 -05:00
Matthew Gordon f52d201ffd
Github build action: WIP 2019-11-28 10:34:35 -05:00
Matthew Gordon ea5e9c3870
Github build action: Install SDL2 2019-11-28 10:33:09 -05:00
Matthew Gordon 32f5939458
Add github action to build and test 2019-11-28 10:21:53 -05:00
73 changed files with 914 additions and 7786 deletions
Show Stats Expand all files Collapse all files

2
.cargo/config
View File

@ -1,2 +0,0 @@
[build]
rustflags = "-C target-cpu=native"

2
.gitattributes vendored
View File

@ -1,2 +0,0 @@
*.png filter=lfs diff=lfs merge=lfs -text
test_data/stanford_bunny.obj filter=lfs diff=lfs merge=lfs -text

BIN
.github/output.png (Stored with Git LFS) vendored
View File

Binary file not shown.

BIN
.github/output2.png (Stored with Git LFS) vendored
View File

Binary file not shown.

BIN
.github/output3.png (Stored with Git LFS) vendored
View File

Binary file not shown.

6
.github/workflows/rust.yml vendored
View File

@ -8,11 +8,7 @@ jobs:
runs-on: ubuntu-latest runs-on: ubuntu-latest
steps: steps:
- uses: actions/checkout@v2 - uses: actions/checkout@v1
- uses: actions-rs/toolchain@v1
with:
profile: minimal
toolchain: nightly
- name: Install dependencies - name: Install dependencies
run: sudo apt-get update && sudo apt-get install libsdl2-dev run: sudo apt-get update && sudo apt-get install libsdl2-dev
- name: Build - name: Build

6
.gitignore vendored
View File

@ -1,9 +1,5 @@
*~ *~
# Linux perf results
perf.data
perf.data.*
# Generated by Cargo # Generated by Cargo
# will have compiled files and executables # will have compiled files and executables
/target/ /target/
@ -21,4 +17,4 @@ Cargo.lock
#already existing elements are commented out #already existing elements are commented out
/target /target
#**/*.rs.bk #**/*.rs.bk

36
Cargo.toml
View File

@ -4,35 +4,13 @@ version = "0.1.0"
authors = ["Matthew Gordon <matthew.scott.gordon@gmail.com>"] authors = ["Matthew Gordon <matthew.scott.gordon@gmail.com>"]
edition = "2018" edition = "2018"
# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
[dependencies] [dependencies]
itertools = "0.14"
obj = "0.9"
quickcheck = "0.9"
quickcheck_macros = "0.9"
rand = "0.7"
rayon = "1.3"
sdl2 = "0.32" sdl2 = "0.32"
csv = "1.1.3" quickcheck = "0.9"
clap = "2.33" quickcheck_macros = "0.8"
png = "0.16"
[dev-dependencies] [dependencies.nalgebra]
criterion = "0.3" version = "0.19"
features = ["arbitrary"]
[[bench]]
name = "simple_scene"
harness = false
[profile.dev]
opt-level = 3
[profile.bench]
opt-level = 3
lto = true
codegen-units = 1
debug = true
[profile.release]
opt-level = 3
lto = true
codegen-units = 1

28
README.md
View File

@ -1,26 +1,2 @@
![](https://github.com/matthewscottgordon/vanrijn/workflows/Rust/badge.svg) # vanrijn
Ray tracer
# Vanrijn
This project is very much a work-in-progress and at this point.
Vanrijn is (or at least will be) a [physically based](https://en.wikipedia.org/wiki/Physically_based_rendering)
[ray tracer](https://en.wikipedia.org/wiki/Ray_tracing_(graphics)). Many thanks to the
authors of the book
["Physically Based Rendering: From Theory to Implementation](https://www.pbrt.org/) from
which many of the algorithms used here are taken. This is, however _not_ a Rust port of
the C++ PBRT rederer described in that book.
This crate is structured as a library; main.rs is just a glorified test harness which
shows an example of using the library to render a scene. It uses SDL2 to display the
rendered image.
On Ubuntu 19.04, if you have the libsdl2-dev package installed you
should be able to run `cargo run --release -- --size 800 600` and see
a window with a test scene rendered into it. In theory it should work
on any platform with SDL2 installed but I've only tested it on Ubuntu
Linux.
![](.github/output3.png?raw=true "Test Image 3")
![](.github/output.png?raw=true "Test Image 1")
![](.github/output2.png?raw=true "Test Image")

53
benches/simple_scene.rs
View File

@ -1,53 +0,0 @@
use criterion::{criterion_group, criterion_main, Criterion};
use vanrijn::colour::{ColourRgbF, NamedColour, Spectrum};
use vanrijn::materials::ReflectiveMaterial;
use vanrijn::math::Vec3;
use vanrijn::mesh::load_obj;
use vanrijn::partial_render_scene;
use vanrijn::raycasting::BoundingVolumeHierarchy;
use vanrijn::scene::Scene;
use vanrijn::util::Tile;
use std::path::Path;
use std::sync::Arc;
fn simple_scene(bencher: &mut Criterion) {
let image_width = 6;
let image_height = 6;
let model_file_path =
Path::new(env!("CARGO_MANIFEST_DIR")).join("test_data/stanford_bunny.obj");
bencher.bench_function("simple_scene", |b| {
let scene = Scene {
camera_location: Vec3::new(-2.0, 1.0, -5.0),
objects: vec![Box::new(BoundingVolumeHierarchy::build(
load_obj(
&model_file_path,
Arc::new(ReflectiveMaterial {
colour: Spectrum::reflection_from_linear_rgb(&ColourRgbF::from_named(
NamedColour::Yellow,
)),
diffuse_strength: 0.05,
reflection_strength: 0.9,
}),
)
.unwrap()
.as_mut_slice(),
))],
};
b.iter(|| {
let tile = Tile {
start_column: 0,
end_column: image_width,
start_row: 0,
end_row: image_height,
};
partial_render_scene(&scene, tile, image_height, image_width);
})
});
}
criterion_group!(benches, simple_scene);
criterion_main!(benches);

1
rust-toolchain
View File

@ -1 +0,0 @@
nightly

328
src/accumulation_buffer.rs
View File

@ -1,328 +0,0 @@
use crate::colour::{ColourXyz, Photon};
use crate::image::{ImageRgbU8, ToneMapper};
use crate::util::{Array2D, Tile};
#[derive(Clone, Debug)]
pub struct AccumulationBuffer {
colour_buffer: Array2D<ColourXyz>,
colour_sum_buffer: Array2D<ColourXyz>,
colour_bias_buffer: Array2D<ColourXyz>,
weight_buffer: Array2D<f64>,
weight_bias_buffer: Array2D<f64>,
}
impl AccumulationBuffer {
pub fn new(height: usize, width: usize) -> AccumulationBuffer {
let colour_buffer = Array2D::new(width, height);
let colour_sum_buffer = Array2D::new(width, height);
let colour_bias_buffer = Array2D::new(width, height);
let weight_buffer = Array2D::new(width, height);
let weight_bias_buffer = Array2D::new(width, height);
AccumulationBuffer {
colour_buffer,
colour_sum_buffer,
colour_bias_buffer,
weight_buffer,
weight_bias_buffer,
}
}
pub fn width(&self) -> usize {
self.colour_buffer.get_width()
}
pub fn height(&self) -> usize {
self.colour_buffer.get_height()
}
pub fn to_image_rgb_u8<Op: ToneMapper<ColourXyz>>(&self, tone_mapper: &Op) -> ImageRgbU8 {
let mut result = ImageRgbU8::new(self.width(), self.height());
tone_mapper.apply_tone_mapping(&self.colour_buffer, &mut result);
result
}
pub fn update_pixel(&mut self, row: usize, column: usize, photon: &Photon, weight: f64) {
let buffer_colour = &mut self.colour_buffer[row][column];
let buffer_colour_sum = &mut self.colour_sum_buffer[row][column];
let buffer_colour_bias = &mut self.colour_bias_buffer[row][column];
let buffer_weight = &mut self.weight_buffer[row][column];
let buffer_weight_bias = &mut self.weight_bias_buffer[row][column];
let photon_colour = ColourXyz::from_photon(photon);
let weight_sum_y = weight - *buffer_weight_bias;
let weight_sum_t = *buffer_weight + weight_sum_y;
*buffer_weight_bias = (weight_sum_t - *buffer_weight) - weight_sum_y;
*buffer_weight = weight_sum_t;
let colour_sum_y = photon_colour.values * weight - buffer_colour_bias.values;
let colour_sum_t = buffer_colour_sum.values + colour_sum_y;
buffer_colour_bias.values = (colour_sum_t - buffer_colour_sum.values) - colour_sum_y;
buffer_colour_sum.values = colour_sum_t;
buffer_colour.values = buffer_colour_sum.values / *buffer_weight;
}
pub fn merge_tile(&mut self, tile: &Tile, src: &AccumulationBuffer) {
assert!(tile.width() == src.width());
assert!(tile.height() == src.height());
for i in 0..tile.height() {
for j in 0..tile.width() {
let dst_colour = &mut self.colour_buffer[tile.start_row + i][tile.start_column + j];
let dst_weight = &mut self.weight_buffer[tile.start_row + i][tile.start_column + j];
*dst_colour = blend(
dst_colour,
*dst_weight,
&src.colour_buffer[i][j],
src.weight_buffer[i][j],
);
*dst_weight += src.weight_buffer[i][j];
}
}
}
}
fn blend(colour1: &ColourXyz, weight1: f64, colour2: &ColourXyz, weight2: f64) -> ColourXyz {
ColourXyz {
values: (colour1.values * weight1 + colour2.values * weight2) * (1.0 / (weight1 + weight2)),
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn has_expected_width() {
let target = AccumulationBuffer::new(16, 12);
assert!(target.width() == 16);
}
#[test]
fn has_expected_height() {
let target = AccumulationBuffer::new(16, 12);
assert!(target.height() == 12);
}
#[test]
fn update_pixel_does_not_panic_inside_bounds() {
let mut target = AccumulationBuffer::new(16, 12);
for i in 0..12 {
for j in 0..16 {
target.update_pixel(i, j, &Default::default(), 1.0);
}
}
}
#[test]
#[should_panic]
fn update_pixel_panics_when_row_to_large() {
let mut target = AccumulationBuffer::new(16, 12);
target.update_pixel(12, 0, &Default::default(), 1.0);
}
#[test]
#[should_panic]
fn update_pixel_panics_when_column_to_large() {
let mut target = AccumulationBuffer::new(16, 12);
target.update_pixel(0, 16, &Default::default(), 1.0);
}
#[test]
fn first_update_sets_expected_value() {
let mut target = AccumulationBuffer::new(16, 12);
let photon = Photon {
wavelength: 589.0,
intensity: 1.5,
};
let row = 4;
let column = 5;
let weight = 0.8;
target.update_pixel(row, column, &photon, weight);
assert!(target.colour_buffer[row][column] == ColourXyz::from_photon(&photon));
assert!(target.weight_buffer[row][column] == weight);
}
#[test]
fn first_update_only_sets_expected_value() {
let mut target = AccumulationBuffer::new(16, 12);
let original = target.clone();
let photon = Photon {
wavelength: 589.0,
intensity: 1.5,
};
let set_row = 4;
let set_column = 5;
target.update_pixel(set_row, set_column, &photon, 0.8);
for i in 0..12 {
for j in 0..16 {
if i != set_row && j != set_column {
assert!(target.colour_buffer[i][j] == original.colour_buffer[i][j]);
assert!(target.weight_buffer[i][j] == original.weight_buffer[i][j]);
}
}
}
}
#[test]
fn second_update_blends_colours() {
let mut target = AccumulationBuffer::new(16, 12);
let photon1 = Photon {
wavelength: 589.0,
intensity: 0.5,
};
let photon2 = Photon {
wavelength: 656.0,
intensity: 1.5,
};
let colour1 = ColourXyz::from_photon(&photon1);
let colour2 = ColourXyz::from_photon(&photon2);
let expected_x = (colour1.x() + colour2.x()) / 2.0;
let expected_y = (colour1.y() + colour2.y()) / 2.0;
let expected_z = (colour1.z() + colour2.z()) / 2.0;
let row = 4;
let column = 5;
target.update_pixel(row, column, &photon1, 1.0);
target.update_pixel(row, column, &photon2, 1.0);
assert!(target.colour_buffer[row][column].x() == expected_x);
assert!(target.colour_buffer[row][column].y() == expected_y);
assert!(target.colour_buffer[row][column].z() == expected_z);
}
#[test]
fn second_update_blends_colours_proportionally() {
let mut target = AccumulationBuffer::new(16, 12);
let photon1 = Photon {
wavelength: 589.0,
intensity: 0.5,
};
let photon2 = Photon {
wavelength: 656.0,
intensity: 1.5,
};
let colour1 = ColourXyz::from_photon(&photon1);
let colour2 = ColourXyz::from_photon(&photon2);
let weight1 = 0.75;
let weight2 = 1.25;
let expected_x = (colour1.x() * weight1 + colour2.x() * weight2) / (weight1 + weight2);
let expected_y = (colour1.y() * weight1 + colour2.y() * weight2) / (weight1 + weight2);
let expected_z = (colour1.z() * weight1 + colour2.z() * weight2) / (weight1 + weight2);
let row = 4;
let column = 5;
target.update_pixel(row, column, &photon1, weight1);
target.update_pixel(row, column, &photon2, weight2);
assert!(target.colour_buffer[row][column].x() == expected_x);
assert!(target.colour_buffer[row][column].y() == expected_y);
assert!(target.colour_buffer[row][column].z() == expected_z);
}
#[test]
fn third_update_blends_colours_proportionally() {
let mut target = AccumulationBuffer::new(16, 12);
let photon1 = Photon {
wavelength: 589.0,
intensity: 0.5,
};
let photon2 = Photon {
wavelength: 656.0,
intensity: 1.5,
};
let photon3 = Photon {
wavelength: 393.0,
intensity: 1.2,
};
let colour1 = ColourXyz::from_photon(&photon1);
let colour2 = ColourXyz::from_photon(&photon2);
let colour3 = ColourXyz::from_photon(&photon3);
let weight1 = 0.75;
let weight2 = 1.25;
let weight3 = 0.5;
let expected_x = (colour1.x() * weight1 + colour2.x() * weight2 + colour3.x() * weight3)
/ (weight1 + weight2 + weight3);
let expected_y = (colour1.y() * weight1 + colour2.y() * weight2 + colour3.y() * weight3)
/ (weight1 + weight2 + weight3);
let expected_z = (colour1.z() * weight1 + colour2.z() * weight2 + colour3.z() * weight3)
/ (weight1 + weight2 + weight3);
let row = 4;
let column = 5;
target.update_pixel(row, column, &photon1, weight1);
target.update_pixel(row, column, &photon2, weight2);
target.update_pixel(row, column, &photon3, weight3);
assert!(target.colour_buffer[row][column].x() == expected_x);
assert!(target.colour_buffer[row][column].y() == expected_y);
assert!(target.colour_buffer[row][column].z() == expected_z);
}
#[test]
fn merge_tile_produces_same_results_as_applying_photons_directly() {
let mut single_buffer = AccumulationBuffer::new(16, 12);
let mut large_buffer = AccumulationBuffer::new(16, 12);
let mut small_buffer = AccumulationBuffer::new(4, 5);
let tile = Tile {
start_column: 3,
end_column: 7,
start_row: 4,
end_row: 9,
};
for i in 0..12 {
for j in 0..16 {
let wavelength = 350.0 + (i * j) as f64;
let intensity = 1.0;
let weight = 0.2 + i as f64 * 0.02 + j as f64 * 0.3;
single_buffer.update_pixel(
i,
j,
&Photon {
wavelength,
intensity,
},
weight,
);
large_buffer.update_pixel(
i,
j,
&Photon {
wavelength,
intensity,
},
weight,
);
}
}
for i in 0..5 {
for j in 0..4 {
let wavelength = 700.0 - (i * j) as f64;
let intensity = 1.0;
let weight = 0.2 + i as f64 * 0.02 + j as f64 * 0.3;
small_buffer.update_pixel(
i,
j,
&Photon {
wavelength,
intensity,
},
weight,
);
single_buffer.update_pixel(
tile.start_row + i,
tile.start_column + j,
&Photon {
wavelength,
intensity,
},
weight,
);
}
}
large_buffer.merge_tile(&tile, &small_buffer);
for i in 0..12 {
for j in 0..16 {
assert!(
(large_buffer.colour_buffer[i][j].values
- single_buffer.colour_buffer[i][j].values)
.norm()
< 0.0000000001
);
assert!(large_buffer.weight_buffer[i][j] == single_buffer.weight_buffer[i][j]);
}
}
}
}

67
src/algebra_utils.rs Normal file
View File

@ -0,0 +1,67 @@
use nalgebra::{Matrix3, RealField, Vector3};
pub fn try_change_of_basis_matrix<T: RealField>(
x: &Vector3<T>,
y: &Vector3<T>,
z: &Vector3<T>,
) -> Option<Matrix3<T>> {
Some(Matrix3::from_rows(&[x.transpose(), y.transpose(), z.transpose()]))
}
#[cfg(test)]
mod tests {
use super::*;
#[cfg(test)]
mod change_of_basis_matrix {
use super::*;
use quickcheck_macros::quickcheck;
#[test]
fn produces_isentity_when_passed_axes() {
let target: Matrix3<f32> = try_change_of_basis_matrix(
&Vector3::x_axis(),
&Vector3::y_axis(),
&Vector3::z_axis(),
)
.unwrap();
assert!(target == Matrix3::identity())
}
#[quickcheck]
fn swap_xy_does_not_change_z(v: Vector3<f32>) {
let target: Matrix3<f32> = try_change_of_basis_matrix(
&Vector3::y_axis(),
&Vector3::x_axis(),
&Vector3::z_axis(),
)
.unwrap();
let v2 = target * v;
assert!(v2.z == v.z)
}
#[quickcheck]
fn swap_xy_copies_y_to_x(v: Vector3<f32>) {
let target: Matrix3<f32> = try_change_of_basis_matrix(
&Vector3::y_axis(),
&Vector3::x_axis(),
&Vector3::z_axis(),
)
.unwrap();
let v2 = target * v;
assert!(v2.x == v.y)
}
#[quickcheck]
fn swap_xy_copies_x_to_y(v: Vector3<f32>) {
let target: Matrix3<f32> = try_change_of_basis_matrix(
&Vector3::y_axis(),
&Vector3::x_axis(),
&Vector3::z_axis(),
)
.unwrap();
let v2 = target * v;
assert!(v2.y == v.x)
}
}
}

186
src/camera.rs
View File

@ -1,31 +1,29 @@
use crate::math::Vec3; use nalgebra::{convert, RealField, Vector3};
use super::accumulation_buffer::AccumulationBuffer; use super::colour::{ColourRgbF, NamedColour};
use super::colour::Photon; use super::image::ImageRgbF;
use super::integrators::{Integrator, SimpleRandomIntegrator}; use super::integrators::{DirectionalLight, Integrator, WhittedIntegrator};
use super::raycasting::Ray; use super::raycasting::Ray;
use super::sampler::Sampler; use super::sampler::Sampler;
use super::scene::Scene; use super::scene::Scene;
use super::util::Tile;
use rand::random; struct ImageSampler<T: RealField> {
image_height_pixels: u32,
image_width_pixels: u32,
struct ImageSampler { film_width: T,
image_height_pixels: usize, film_height: T,
image_width_pixels: usize,
film_width: f64, camera_location: Vector3<T>,
film_height: f64, film_distance: T,
camera_location: Vec3<f64>,
film_distance: f64,
} }
impl ImageSampler { impl<T: RealField> ImageSampler<T> {
pub fn new(width: usize, height: usize, camera_location: Vec3<f64>) -> ImageSampler { pub fn new(width: u32, height: u32, camera_location: Vector3<T>) -> ImageSampler<T> {
let (film_width, film_height) = { let (film_width, film_height) = {
let width = width as f64; let width: T = convert(width as f64);
let height = height as f64; let height: T = convert(height as f64);
let film_size = 1.0; let film_size: T = convert(1.0);
if width > height { if width > height {
(width / height, film_size) (width / height, film_size)
} else { } else {
@ -35,98 +33,87 @@ impl ImageSampler {
ImageSampler { ImageSampler {
image_height_pixels: height, image_height_pixels: height,
image_width_pixels: width, image_width_pixels: width,
film_distance: 1.0, film_distance: convert(1.0),
film_width, film_width,
film_height, film_height,
camera_location, camera_location,
} }
} }
fn scale(i: usize, n: usize, l: f64) -> f64 { fn scale(i: u32, n: u32, l: T) -> T {
let n = n as f64; let one: T = convert(1.0);
let i = i as f64; let n: T = convert(n as f64);
let pixel_size = l * (1.0 / n); let i: T = convert(i as f64);
(i + random::<f64>()) * pixel_size let pixel_size: T = l * (one / n);
(i + convert(0.5)) * pixel_size
} }
fn ray_for_pixel(&self, row: usize, column: usize) -> Ray { fn ray_for_pixel(&self, row: u32, column: u32) -> Ray<T> {
Ray::new( Ray::new(
self.camera_location, self.camera_location,
Vec3::new( Vector3::new(
Self::scale(column, self.image_width_pixels, self.film_width) Self::scale(column, self.image_width_pixels, self.film_width)
- self.film_width * 0.5, - self.film_width * convert(0.5),
Self::scale( Self::scale(row, self.image_height_pixels, self.film_height)
self.image_height_pixels - (row + 1), - self.film_height * convert(0.5),
self.image_height_pixels,
self.film_height,
) - self.film_height * 0.5,
self.film_distance, self.film_distance,
), ),
) )
} }
} }
const RECURSION_LIMIT: u16 = 128; pub fn render_scene<T: RealField>(output_image: &mut ImageRgbF<T>, scene: &Scene<T>) {
partial_render_scene(
output_image,
scene,
0,
output_image.get_height(),
0,
output_image.get_width(),
)
}
/// Render a rectangular section of the image. pub fn partial_render_scene<T: RealField>(
/// output_image: &mut ImageRgbF<T>,
/// The contents and the image, along with the camera, are defined by `scene`. scene: &Scene<T>,
/// row_start: u32,
/// Assuming an overall image size given by `width` and `height`, the part of the image row_end: u32,
/// defined by `tile` is rendered and returned. Rendering a tile at a time allows a partially- column_start: u32,
/// rendered image to be displayed to the user. column_end: u32,
/// ) {
/// # Examples let image_sampler = ImageSampler::new(
// output_image.get_width(),
/// ``` output_image.get_height(),
/// # use vanrijn::math::Vec3; scene.camera_location,
/// # use vanrijn::scene::Scene; );
/// # use vanrijn::util::TileIterator; let ambient_intensity: T = convert(0.0);
/// # use vanrijn::partial_render_scene; let directional_intensity1: T = convert(7.0);
/// # let scene = Scene { camera_location: Vec3::new(0.0, 0.0, 0.0), objects: vec![] }; let directional_intensity2: T = convert(3.0);
/// let image_width = 640; let integrator = WhittedIntegrator::<T> {
/// let image_height = 480; ambient_light: ColourRgbF::from_named(NamedColour::White) * ambient_intensity,
/// let time_size = 32; lights: vec![
/// for tile in TileIterator::new(640, 480, 32) { DirectionalLight {
/// let tile_image = partial_render_scene( &scene, tile, image_height, image_width ); direction: Vector3::new(convert(1.0), convert(1.0), convert(-1.0)).normalize(),
/// // display and/or save tile_image colour: ColourRgbF::from_named(NamedColour::White) * directional_intensity1,
/// } },
/// ``` DirectionalLight {
pub fn partial_render_scene( direction: Vector3::new(convert(-0.5), convert(2.0), convert(-0.5)).normalize(),
scene: &Scene, colour: ColourRgbF::from_named(NamedColour::White) * directional_intensity2,
tile: Tile, },
height: usize, ],
width: usize, };
) -> AccumulationBuffer {
let mut output_image_tile = AccumulationBuffer::new(tile.width(), tile.height());
let image_sampler = ImageSampler::new(width, height, scene.camera_location);
let integrator = SimpleRandomIntegrator {};
let sampler = Sampler { scene }; let sampler = Sampler { scene };
for column in 0..tile.width() { for column in column_start..column_end {
for row in 0..tile.height() { for row in row_start..row_end {
let ray = image_sampler.ray_for_pixel(tile.start_row + row, tile.start_column + column); let ray = image_sampler.ray_for_pixel(row, column);
let hit = sampler.sample(&ray); let hit = sampler.sample(&ray);
let photon = match hit { let colour = match hit {
None => Photon { None => ColourRgbF::from_named(NamedColour::Black),
wavelength: 0.0, Some(intersection_info) => integrator.integrate(&sampler, &intersection_info),
intensity: 0.0,
},
Some(intersection_info) => integrator.integrate(
&sampler,
&intersection_info,
&Photon::random_wavelength(),
RECURSION_LIMIT,
),
}; };
output_image_tile.update_pixel( output_image.set_colour(row, column, colour);
row,
column,
&photon.scale_intensity(Photon::random_wavelength_pdf(photon.wavelength)),
1.0,
);
} }
} }
output_image_tile
} }
#[cfg(test)] #[cfg(test)]
@ -134,7 +121,7 @@ mod tests {
use super::*; use super::*;
use crate::materials::LambertianMaterial; use crate::materials::LambertianMaterial;
use crate::raycasting::{Intersect, IntersectionInfo, Plane}; use crate::raycasting::{Intersect, IntersectionInfo, Plane};
use std::sync::Arc; use std::rc::Rc;
#[cfg(test)] #[cfg(test)]
mod imagesampler { mod imagesampler {
@ -143,23 +130,23 @@ mod tests {
#[test] #[test]
fn scale_returns_correct_value_for_zero() { fn scale_returns_correct_value_for_zero() {
let correct_value = (3.0 / 10.0) / 2.0; let correct_value = (3.0 / 10.0) / 2.0;
assert!((ImageSampler::scale(0, 10, 3.0f64) - correct_value).abs() < 0.5) assert!((ImageSampler::scale(0, 10, 3.0f64) - correct_value).abs() < 0.0000000001)
} }
#[test] #[test]
fn scale_returns_correct_value_for_last_pixel() { fn scale_returns_correct_value_for_last_pixel() {
let correct_value = 3.0 - (3.0 / 10.0) / 2.0; let correct_value = 3.0 - (3.0 / 10.0) / 2.0;
assert!((ImageSampler::scale(9, 10, 3.0f64) - correct_value).abs() < 0.5) assert!((ImageSampler::scale(9, 10, 3.0f64) - correct_value).abs() < 0.0000000001)
} }
#[test] #[test]
fn ray_for_pixel_returns_value_that_intersects_film_plane_at_expected_location() { fn ray_for_pixel_returns_value_that_intersects_film_plane_at_expected_location() {
let target = ImageSampler::new(800, 600, Vec3::new(0.0, 0.0, 0.0)); let target = ImageSampler::new(800, 600, Vector3::new(0.0, 0.0, 0.0));
let ray = target.ray_for_pixel(100, 200); let ray = target.ray_for_pixel(100, 200);
let film_plane = Plane::new( let film_plane = Plane::new(
Vec3::new(0.0, 0.0, 1.0), Vector3::new(0.0, 0.0, 1.0),
target.film_distance, target.film_distance,
Arc::new(LambertianMaterial::new_dummy()), Rc::new(LambertianMaterial::<f64>::new_dummy()),
); );
let point_on_film_plane = match film_plane.intersect(&ray) { let point_on_film_plane = match film_plane.intersect(&ray) {
Some(IntersectionInfo { Some(IntersectionInfo {
@ -173,12 +160,13 @@ mod tests {
}) => location, }) => location,
None => panic!(), None => panic!(),
}; };
let expected_x = target.film_width * (200.0/800.0 - 0.5); let expected_x: f64 =
assert!(point_on_film_plane.x() - expected_x < target.film_width / 800.0); ImageSampler::scale(200, 800, target.film_width) - target.film_width * 0.5;
assert!(point_on_film_plane.x() - expected_x >= 0.0); print!("{}, {}", expected_x, point_on_film_plane);
let expected_y = -target.film_height * (100.0/600.0 - 0.5); assert!((point_on_film_plane.x - expected_x).abs() < 0.0000000001);
assert!(expected_y - point_on_film_plane.y() < target.film_height / 600.0); let expected_y =
assert!(expected_y - point_on_film_plane.y() >= 0.0); ImageSampler::scale(100, 600, target.film_height) - target.film_height * 0.5;
assert!((point_on_film_plane.y - expected_y).abs() < 0.0000000001);
} }
} }
} }

113
src/colour/colour_rgb.rs → src/colour.rs
View File

@ -1,56 +1,59 @@
use crate::math::Vec3; use nalgebra::{convert, RealField, Vector3};
use std::ops::{Add, Mul}; use std::ops::{Add, Mul};
#[derive(Copy, Clone, Debug, Default)] #[derive(Copy, Clone, Debug)]
pub struct ColourRgbF { pub struct ColourRgbF<T: RealField> {
pub values: Vec3<f64>, values: Vector3<T>,
} }
impl ColourRgbF { impl<T: RealField> ColourRgbF<T> {
pub fn new(red: f64, green: f64, blue: f64) -> ColourRgbF { pub fn new(red: T, green: T, blue: T) -> ColourRgbF<T> {
ColourRgbF { ColourRgbF {
values: Vec3::new(red, green, blue), values: Vector3::new(red, green, blue),
} }
} }
pub fn from_named(name: NamedColour) -> ColourRgbF { pub fn from_named(name: NamedColour) -> ColourRgbF<T> {
let zero: T = convert(0.0);
let half: T = convert(0.5);
let one: T = convert(1.0);
match name { match name {
NamedColour::Black => ColourRgbF::new(0.0, 0.0, 0.0), NamedColour::Black => ColourRgbF::new(zero, zero, zero),
NamedColour::White => ColourRgbF::new(1.0, 1.0, 1.0), NamedColour::White => ColourRgbF::new(one, one, one),
NamedColour::Red => ColourRgbF::new(1.0, 0.0, 0.0), NamedColour::Red => ColourRgbF::new(one, zero, zero),
NamedColour::Lime => ColourRgbF::new(0.0, 1.0, 0.0), NamedColour::Lime => ColourRgbF::new(zero, one, zero),
NamedColour::Blue => ColourRgbF::new(0.0, 0.0, 1.0), NamedColour::Blue => ColourRgbF::new(zero, zero, one),
NamedColour::Yellow => ColourRgbF::new(1.0, 1.0, 0.0), NamedColour::Yellow => ColourRgbF::new(one, one, zero),
NamedColour::Cyan => ColourRgbF::new(0.0, 1.0, 1.0), NamedColour::Cyan => ColourRgbF::new(zero, one, one),
NamedColour::Magenta => ColourRgbF::new(1.0, 0.0, 1.0), NamedColour::Magenta => ColourRgbF::new(one, zero, one),
NamedColour::Gray => ColourRgbF::new(0.5, 0.5, 0.5), NamedColour::Gray => ColourRgbF::new(half, half, half),
NamedColour::Maroon => ColourRgbF::new(0.5, 0.0, 0.0), NamedColour::Maroon => ColourRgbF::new(half, zero, zero),
NamedColour::Olive => ColourRgbF::new(0.5, 0.5, 0.0), NamedColour::Olive => ColourRgbF::new(half, half, zero),
NamedColour::Green => ColourRgbF::new(0.0, 0.5, 0.0), NamedColour::Green => ColourRgbF::new(zero, half, zero),
NamedColour::Purple => ColourRgbF::new(0.5, 0.0, 0.5), NamedColour::Purple => ColourRgbF::new(half, zero, half),
NamedColour::Teal => ColourRgbF::new(0.0, 0.5, 0.5), NamedColour::Teal => ColourRgbF::new(zero, half, half),
NamedColour::Navy => ColourRgbF::new(0.0, 0.0, 0.5), NamedColour::Navy => ColourRgbF::new(zero, zero, half),
} }
} }
pub fn from_vec3(v: &Vec3<f64>) -> ColourRgbF { pub fn from_vector3(v: &Vector3<T>) -> ColourRgbF<T> {
ColourRgbF { values: *v } ColourRgbF { values: *v }
} }
pub fn red(&self) -> f64 { pub fn red(&self) -> T {
self.values.x() self.values[0]
} }
pub fn green(&self) -> f64 { pub fn green(&self) -> T {
self.values.y() self.values[1]
} }
pub fn blue(&self) -> f64 { pub fn blue(&self) -> T {
self.values.z() self.values[2]
} }
pub fn as_vec3(&self) -> &Vec3<f64> { pub fn as_vector3(&self) -> &Vector3<T> {
&self.values &self.values
} }
} }
@ -77,27 +80,27 @@ pub enum NamedColour {
Navy, Navy,
} }
impl Add<ColourRgbF> for ColourRgbF { impl<T: RealField> Add<ColourRgbF<T>> for ColourRgbF<T> {
type Output = ColourRgbF; type Output = ColourRgbF<T>;
fn add(self, rhs: ColourRgbF) -> ColourRgbF { fn add(self, rhs: ColourRgbF<T>) -> ColourRgbF<T> {
ColourRgbF { ColourRgbF {
values: self.values + rhs.values, values: self.values + rhs.values,
} }
} }
} }
impl Mul<f64> for ColourRgbF { impl<T: RealField> Mul<T> for ColourRgbF<T> {
type Output = ColourRgbF; type Output = ColourRgbF<T>;
fn mul(self, rhs: f64) -> ColourRgbF { fn mul(self, rhs: T) -> ColourRgbF<T> {
ColourRgbF { ColourRgbF {
values: self.values * rhs, values: self.values * rhs,
} }
} }
} }
impl Mul<ColourRgbF> for ColourRgbF { impl<T: RealField> Mul<ColourRgbF<T>> for ColourRgbF<T> {
type Output = ColourRgbF; type Output = ColourRgbF<T>;
fn mul(self, rhs: ColourRgbF) -> ColourRgbF { fn mul(self, rhs: ColourRgbF<T>) -> ColourRgbF<T> {
ColourRgbF { ColourRgbF {
values: self.values.component_mul(&rhs.values), values: self.values.component_mul(&rhs.values),
} }
@ -112,9 +115,9 @@ mod tests {
use super::*; use super::*;
use quickcheck::{Arbitrary, Gen}; use quickcheck::{Arbitrary, Gen};
use quickcheck_macros::quickcheck; use quickcheck_macros::quickcheck;
impl Arbitrary for ColourRgbF { impl<T: Arbitrary + RealField> Arbitrary for ColourRgbF<T> {
fn arbitrary<G: Gen>(g: &mut G) -> ColourRgbF { fn arbitrary<G: Gen>(g: &mut G) -> ColourRgbF<T> {
let values = <Vec3<f64> as Arbitrary>::arbitrary(g); let values = <Vector3<T> as Arbitrary>::arbitrary(g);
ColourRgbF { values } ColourRgbF { values }
} }
} }
@ -130,12 +133,12 @@ mod tests {
#[test] #[test]
fn as_vector3_returns_expected_vector() { fn as_vector3_returns_expected_vector() {
let target = ColourRgbF::new(1.0, 2.0, 3.0); let target = ColourRgbF::new(1.0, 2.0, 3.0);
let result = target.as_vec3(); let result = target.as_vector3();
assert!(result.x() == 1.0); assert!(result.x == 1.0);
} }
#[quickcheck] #[quickcheck]
fn any_colour_multiplied_by_zero_is_black(colour: ColourRgbF) { fn any_colour_multiplied_by_zero_is_black(colour: ColourRgbF<f64>) {
let target = colour * 0.0; let target = colour * 0.0;
assert!(target.red() == 0.0); assert!(target.red() == 0.0);
assert!(target.green() == 0.0); assert!(target.green() == 0.0);
@ -143,14 +146,17 @@ mod tests {
} }
#[quickcheck] #[quickcheck]
fn red_channel_multiplied_by_scalar_yields_correct_result(colour: ColourRgbF, scalar: f64) { fn red_channel_multiplied_by_scalar_yields_correct_result(
colour: ColourRgbF<f64>,
scalar: f64,
) {
let target = colour * scalar; let target = colour * scalar;
assert!(target.red() == colour.red() * scalar); assert!(target.red() == colour.red() * scalar);
} }
#[quickcheck] #[quickcheck]
fn green_channel_multiplied_by_scalar_yields_correct_result( fn green_channel_multiplied_by_scalar_yields_correct_result(
colour: ColourRgbF, colour: ColourRgbF<f64>,
scalar: f64, scalar: f64,
) { ) {
let target = colour * scalar; let target = colour * scalar;
@ -159,7 +165,7 @@ mod tests {
#[quickcheck] #[quickcheck]
fn blue_channel_multiplied_by_scalar_yields_correct_result( fn blue_channel_multiplied_by_scalar_yields_correct_result(
colour: ColourRgbF, colour: ColourRgbF<f64>,
scalar: f64, scalar: f64,
) { ) {
let target = colour * scalar; let target = colour * scalar;
@ -167,7 +173,10 @@ mod tests {
} }
#[quickcheck] #[quickcheck]
fn adding_colourrgbf_adds_individual_channels(colour1: ColourRgbF, colour2: ColourRgbF) { fn adding_colourrgbf_adds_individual_channels(
colour1: ColourRgbF<f64>,
colour2: ColourRgbF<f64>,
) {
let target = colour1 + colour2; let target = colour1 + colour2;
assert!(target.red() == colour1.red() + colour2.red()); assert!(target.red() == colour1.red() + colour2.red());
assert!(target.green() == colour1.green() + colour2.green()); assert!(target.green() == colour1.green() + colour2.green());
@ -176,8 +185,8 @@ mod tests {
#[quickcheck] #[quickcheck]
fn multiplying_colourrgbf_adds_individual_channels( fn multiplying_colourrgbf_adds_individual_channels(
colour1: ColourRgbF, colour1: ColourRgbF<f64>,
colour2: ColourRgbF, colour2: ColourRgbF<f64>,
) { ) {
let target = colour1 * colour2; let target = colour1 * colour2;
assert!(target.red() == colour1.red() * colour2.red()); assert!(target.red() == colour1.red() * colour2.red());

134
src/colour/colour_xyz.rs
View File

@ -1,134 +0,0 @@
use crate::math::{Mat3, Vec3};
use super::{ColourRgbF, Photon};
/// A CIE XYZ Colour Value
#[derive(Clone, Debug, Default, PartialEq)]
pub struct ColourXyz {
pub values: Vec3<f64>,
}
impl ColourXyz {
/// Construct a ColourXyz with the specified XYZ values
pub fn new(x: f64, y: f64, z: f64) -> ColourXyz {
ColourXyz {
values: Vec3::new(x, y, z),
}
}
/// Calculate the XYZ colour of a laser light with the given wavelength
///
/// The wavelength is in nanometres.
pub fn for_wavelength(wavelength: f64) -> ColourXyz {
let values = Vec3::new(
colour_matching_function_x(wavelength),
colour_matching_function_y(wavelength),
colour_matching_function_z(wavelength),
);
ColourXyz { values }
}
pub fn from_photon(photon: &Photon) -> ColourXyz {
let mut result = Self::for_wavelength(photon.wavelength);
result.values *= photon.intensity;
result
}
pub fn x(&self) -> f64 {
self.values.x()
}
pub fn y(&self) -> f64 {
self.values.y()
}
pub fn z(&self) -> f64 {
self.values.z()
}
pub fn to_linear_rgb(&self) -> ColourRgbF {
let transform = Mat3::from_rows(
&Vec3::new(3.24096994, -1.53738318, -0.49861076),
&Vec3::new(-0.96924364, 1.87596750, 0.04155506),
&Vec3::new(0.05563008, -0.20397696, 1.05697151),
);
ColourRgbF::from_vec3(&(transform * self.values))
}
pub fn from_linear_rgb(rgb: &ColourRgbF) -> ColourXyz {
let transform = Mat3::from_rows(
&Vec3::new(0.41239080, 0.35758434, 0.18048079),
&Vec3::new(0.21263901, 0.71516868, 0.07219232),
&Vec3::new(0.01933082, 0.11919478, 0.95053215),
);
ColourXyz {
values: transform * rgb.values,
}
}
pub fn to_srgb(&self) -> ColourRgbF {
let mut srgb = self.to_linear_rgb();
for element in srgb.values.coords.iter_mut() {
*element = srgb_gamma(*element);
}
srgb
}
}
fn srgb_gamma(u: f64) -> f64 {
if u <= 0.0031308 {
12.98 * u
} else {
1.005 * u.powf(1.0 / 2.4) - 0.055
}
}
fn gaussian(wavelength: f64, alpha: f64, mu: f64, sigma1: f64, sigma2: f64) -> f64 {
let denominator = 2.0 * (if wavelength < mu { sigma1 } else { sigma2 }).powi(2);
alpha * (-(wavelength - mu).powi(2) / denominator).exp()
}
fn colour_matching_function_x(wavelength: f64) -> f64 {
gaussian(wavelength, 1.056, 599.8, 37.9, 31.0)
+ gaussian(wavelength, 0.362, 442.0, 16.0, 26.7)
+ gaussian(wavelength, -0.065, 501.1, 20.4, 26.2)
}
fn colour_matching_function_y(wavelength: f64) -> f64 {
gaussian(wavelength, 0.821, 568.8, 46.9, 40.5) + gaussian(wavelength, 0.286, 530.9, 16.3, 31.1)
}
fn colour_matching_function_z(wavelength: f64) -> f64 {
gaussian(wavelength, 1.217, 437.0, 11.8, 36.0) + gaussian(wavelength, 0.681, 459.0, 26.0, 13.8)
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn x_returns_zero_for_default() {
let target: ColourXyz = Default::default();
assert!(target.x() == 0.0);
}
#[test]
fn x_returns_specified_value_after_constructioni_with_new() {
let target = ColourXyz::new(0.1, 0.2, 0.3);
assert!(target.x() == 0.1);
}
#[test]
fn z_returns_specified_value_after_constructioni_with_new() {
let target = ColourXyz::new(0.1, 0.2, 0.3);
assert!(target.z() == 0.3);
}
#[test]
fn roundtrip_to_linear_rgb_yields_original_values() {
let target = ColourXyz::new(0.1, 0.2, 0.3);
let rgb = target.to_linear_rgb();
let xyz = ColourXyz::from_linear_rgb(&rgb);
assert!((target.values - xyz.values).norm() < 0.00000001);
}
}

14
src/colour/mod.rs
View File

@ -1,14 +0,0 @@
pub mod colour_rgb;
pub use colour_rgb::{ColourRgbF, ColourRgbU8, NamedColour};
pub mod photon;
pub use photon::Photon;
pub mod colour_xyz;
pub use colour_xyz::ColourXyz;
pub mod spectrum;
pub use spectrum::Spectrum;
pub const SHORTEST_VISIBLE_WAVELENGTH: f64 = 380.0;
pub const LONGEST_VISIBLE_WAVELENGTH: f64 = 740.0;

43
src/colour/photon.rs
View File

@ -1,43 +0,0 @@
use crate::colour::{LONGEST_VISIBLE_WAVELENGTH, SHORTEST_VISIBLE_WAVELENGTH};
use rand::random;
/// A quantum of light with a given wavelength and intensity
#[derive(Clone, Default, Debug)]
pub struct Photon {
/// The wavelength in nanometres
pub wavelength: f64,
/// The intensity of the light
///
/// Depending on context, this might represent actual intensity in W/sr,
/// radiant flux in W, irradiance in W/m^2, or radiance in W/(m^2sr).
pub intensity: f64,
}
impl Photon {
pub fn random_wavelength() -> Photon {
Photon {
wavelength: SHORTEST_VISIBLE_WAVELENGTH
+ (LONGEST_VISIBLE_WAVELENGTH - SHORTEST_VISIBLE_WAVELENGTH) * random::<f64>(),
intensity: 0.0,
}
}
pub fn random_wavelength_pdf(_wavelength: f64) -> f64 {
LONGEST_VISIBLE_WAVELENGTH - SHORTEST_VISIBLE_WAVELENGTH
}
pub fn scale_intensity(&self, scale_factor: f64) -> Photon {
Photon {
wavelength: self.wavelength,
intensity: self.intensity * scale_factor,
}
}
pub fn set_intensity(&self, intensity: f64) -> Photon {
Photon {
wavelength: self.wavelength,
intensity,
}
}
}

491
src/colour/spectrum.rs
View File

@ -1,491 +0,0 @@
use crate::colour::{ColourRgbF, Photon, LONGEST_VISIBLE_WAVELENGTH, SHORTEST_VISIBLE_WAVELENGTH};
use itertools::izip;
#[derive(Debug)]
pub struct Spectrum {
shortest_wavelength: f64,
longest_wavelength: f64,
samples: Vec<f64>,
}
impl Spectrum {
pub fn black() -> Spectrum {
Spectrum {
shortest_wavelength: SHORTEST_VISIBLE_WAVELENGTH,
longest_wavelength: LONGEST_VISIBLE_WAVELENGTH,
samples: vec![0.0, 0.0],
}
}
pub fn grey(brightness: f64) -> Spectrum {
Spectrum {
shortest_wavelength: SHORTEST_VISIBLE_WAVELENGTH,
longest_wavelength: LONGEST_VISIBLE_WAVELENGTH,
samples: vec![brightness; 2],
}
}
pub fn diamond_index_of_refraction() -> Spectrum {
Spectrum {
shortest_wavelength: 326.27,
longest_wavelength: 774.9,
samples: vec![
2.505813241,
2.487866556,
2.473323675,
2.464986815,
2.455051934,
2.441251728,
2.431478974,
2.427076431,
2.420857286,
2.411429037,
2.406543164,
2.406202402,
],
}
}
fn wavelength_range(&self) -> f64 {
self.longest_wavelength - self.shortest_wavelength
}
fn index_at_or_before_wavelength(&self, wavelength: f64) -> usize {
((self.samples.len() - 1) as f64
* ((wavelength - self.shortest_wavelength) / self.wavelength_range())) as usize
}
fn wavelength_at_index(&self, index: usize) -> f64 {
(index as f64) / ((self.samples.len() - 1) as f64) * self.wavelength_range()
+ self.shortest_wavelength
}
pub fn intensity_at_wavelength(&self, wavelength: f64) -> f64 {
if wavelength < self.shortest_wavelength || wavelength > self.longest_wavelength {
0.0
} else {
let index_before = self.index_at_or_before_wavelength(wavelength);
let wavelength_before = self.wavelength_at_index(index_before);
if index_before == self.samples.len() - 1 {
self.samples[index_before]
} else {
let wavelength_after = self.wavelength_at_index(index_before + 1);
let delta = wavelength_after - wavelength_before;
let ratio = (wavelength - wavelength_before) / delta;
self.samples[index_before] * (1.0 - ratio) + self.samples[index_before + 1] * ratio
}
}
}
pub fn reflection_from_linear_rgb(colour: &ColourRgbF) -> Spectrum {
Spectrum {
shortest_wavelength: rgb_reference_spectrum::SHORTEST_WAVELENGTH,
longest_wavelength: rgb_reference_spectrum::LONGEST_WAVELENGTH,
#[allow(clippy::collapsible_else_if)]
samples: if colour.red() <= colour.green() && colour.red() <= colour.blue() {
if colour.green() <= colour.blue() {
izip![
rgb_reference_spectrum::reflection::WHITE.iter(),
rgb_reference_spectrum::reflection::CYAN.iter(),
rgb_reference_spectrum::reflection::BLUE.iter()
]
.map(|(white, cyan, blue)| {
colour.red() * white
+ (colour.green() - colour.red()) * cyan
+ (colour.blue() - colour.green()) * blue
})
.collect()
} else {
izip![
rgb_reference_spectrum::reflection::WHITE.iter(),
rgb_reference_spectrum::reflection::CYAN.iter(),
rgb_reference_spectrum::reflection::GREEN.iter()
]
.map(|(white, cyan, green)| {
colour.red() * white
+ (colour.blue() - colour.red()) * cyan
+ (colour.green() - colour.blue()) * green
})
.collect()
}
} else if colour.green() <= colour.red() && colour.green() < colour.blue() {
if colour.red() <= colour.blue() {
izip![
rgb_reference_spectrum::reflection::WHITE.iter(),
rgb_reference_spectrum::reflection::MAGENTA.iter(),
rgb_reference_spectrum::reflection::BLUE.iter()
]
.map(|(white, magenta, blue)| {
colour.green() * white
+ (colour.red() - colour.green()) * magenta
+ (colour.blue() - colour.red()) * blue
})
.collect()
} else {
izip![
rgb_reference_spectrum::reflection::WHITE.iter(),
rgb_reference_spectrum::reflection::MAGENTA.iter(),
rgb_reference_spectrum::reflection::RED.iter()
]
.map(|(white, magenta, red)| {
colour.green() * white
+ (colour.blue() - colour.green()) * magenta
+ (colour.red() - colour.blue()) * red
})
.collect()
}
} else {
if colour.red() <= colour.green() {
izip![
rgb_reference_spectrum::reflection::WHITE.iter(),
rgb_reference_spectrum::reflection::YELLOW.iter(),
rgb_reference_spectrum::reflection::GREEN.iter()
]
.map(|(white, yellow, green)| {
colour.blue() * white
+ (colour.red() - colour.blue()) * yellow
+ (colour.green() - colour.red()) * green
})
.collect()
} else {
izip![
rgb_reference_spectrum::reflection::WHITE.iter(),
rgb_reference_spectrum::reflection::YELLOW.iter(),
rgb_reference_spectrum::reflection::RED.iter()
]
.map(|(white, yellow, red)| {
colour.blue() * white
+ (colour.green() - colour.blue()) * yellow
+ (colour.red() - colour.green()) * red
})
.collect()
}
},
}
}
pub fn scale_photon(&self, photon: &Photon) -> Photon {
let wavelength = photon.wavelength;
photon.scale_intensity(self.intensity_at_wavelength(wavelength))
}
pub fn emit_photon(&self, photon: &Photon) -> Photon {
let wavelength = photon.wavelength;
photon.set_intensity(self.intensity_at_wavelength(wavelength))
}
}
mod rgb_reference_spectrum {
pub const SHORTEST_WAVELENGTH: f64 = 380.0;
pub const LONGEST_WAVELENGTH: f64 = 720.0;
#[allow(clippy::excessive_precision)]
pub mod reflection {
pub const WHITE: [f64; 32] = [
1.0618958571272863e+00,
1.0615019980348779e+00,
1.0614335379927147e+00,
1.0622711654692485e+00,
1.0622036218416742e+00,
1.0625059965187085e+00,
1.0623938486985884e+00,
1.0624706448043137e+00,
1.0625048144827762e+00,
1.0624366131308856e+00,
1.0620694238892607e+00,
1.0613167586932164e+00,
1.0610334029377020e+00,
1.0613868564828413e+00,
1.0614215366116762e+00,
1.0620336151299086e+00,
1.0625497454805051e+00,
1.0624317487992085e+00,
1.0625249140554480e+00,
1.0624277664486914e+00,
1.0624749854090769e+00,
1.0625538581025402e+00,
1.0625326910104864e+00,
1.0623922312225325e+00,
1.0623650980354129e+00,
1.0625256476715284e+00,
1.0612277619533155e+00,
1.0594262608698046e+00,
1.0599810758292072e+00,
1.0602547314449409e+00,
1.0601263046243634e+00,
1.0606565756823634e+00,
];
pub const CYAN: [f64; 32] = [
1.0414628021426751e+00,
1.0328661533771188e+00,
1.0126146228964314e+00,
1.0350460524836209e+00,
1.0078661447098567e+00,
1.0422280385081280e+00,
1.0442596738499825e+00,
1.0535238290294409e+00,
1.0180776226938120e+00,
1.0442729908727713e+00,
1.0529362541920750e+00,
1.0537034271160244e+00,
1.0533901869215969e+00,
1.0537782700979574e+00,
1.0527093770467102e+00,
1.0530449040446797e+00,
1.0550554640191208e+00,
1.0553673610724821e+00,
1.0454306634683976e+00,
6.2348950639230805e-01,
1.8038071613188977e-01,
-7.6303759201984539e-03,
-1.5217847035781367e-04,
-7.5102257347258311e-03,
-2.1708639328491472e-03,
6.5919466602369636e-04,
1.2278815318539780e-02,
-4.4669775637208031e-03,
1.7119799082865147e-02,
4.9211089759759801e-03,
5.8762925143334985e-03,
2.5259399415550079e-02,
];
pub const MAGENTA: [f64; 32] = [
9.9422138151236850e-01,
9.8986937122975682e-01,
9.8293658286116958e-01,
9.9627868399859310e-01,
1.0198955019000133e+00,
1.0166395501210359e+00,
1.0220913178757398e+00,
9.9651666040682441e-01,
1.0097766178917882e+00,
1.0215422470827016e+00,
6.4031953387790963e-01,
2.5012379477078184e-03,
6.5339939555769944e-03,
2.8334080462675826e-03,
-5.1209675389074505e-11,
-9.0592291646646381e-03,
3.3936718323331200e-03,
-3.0638741121828406e-03,
2.2203936168286292e-01,
6.3141140024811970e-01,
9.7480985576500956e-01,
9.7209562333590571e-01,
1.0173770302868150e+00,
9.9875194322734129e-01,
9.4701725739602238e-01,
8.5258623154354796e-01,
9.4897798581660842e-01,
9.4751876096521492e-01,
9.9598944191059791e-01,
8.6301351503809076e-01,
8.9150987853523145e-01,
8.4866492652845082e-01,
];
pub const YELLOW: [f64; 32] = [
5.5740622924920873e-03,
-4.7982831631446787e-03,
-5.2536564298613798e-03,
-6.4571480044499710e-03,
-5.9693514658007013e-03,
-2.1836716037686721e-03,
1.6781120601055327e-02,
9.6096355429062641e-02,
2.1217357081986446e-01,
3.6169133290685068e-01,
5.3961011543232529e-01,
7.4408810492171507e-01,
9.2209571148394054e-01,
1.0460304298411225e+00,
1.0513824989063714e+00,
1.0511991822135085e+00,
1.0510530911991052e+00,
1.0517397230360510e+00,
1.0516043086790485e+00,
1.0511944032061460e+00,
1.0511590325868068e+00,
1.0516612465483031e+00,
1.0514038526836869e+00,
1.0515941029228475e+00,
1.0511460436960840e+00,
1.0515123758830476e+00,
1.0508871369510702e+00,
1.0508923708102380e+00,
1.0477492815668303e+00,
1.0493272144017338e+00,
1.0435963333422726e+00,
1.0392280772051465e+00,
];
pub const RED: [f64; 32] = [
1.6575604867086180e-01,
1.1846442802747797e-01,
1.2408293329637447e-01,
1.1371272058349924e-01,
7.8992434518899132e-02,
3.2205603593106549e-02,
-1.0798365407877875e-02,
1.8051975516730392e-02,
5.3407196598730527e-03,
1.3654918729501336e-02,
-5.9564213545642841e-03,
-1.8444365067353252e-03,
-1.0571884361529504e-02,
-2.9375521078000011e-03,
-1.0790476271835936e-02,
-8.0224306697503633e-03,
-2.2669167702495940e-03,
7.0200240494706634e-03,
-8.1528469000299308e-03,
6.0772866969252792e-01,
9.8831560865432400e-01,
9.9391691044078823e-01,
1.0039338994753197e+00,
9.9234499861167125e-01,
9.9926530858855522e-01,
1.0084621557617270e+00,
9.8358296827441216e-01,
1.0085023660099048e+00,
9.7451138326568698e-01,
9.8543269570059944e-01,
9.3495763980962043e-01,
9.8713907792319400e-01,
];
pub const GREEN: [f64; 32] = [
2.6494153587602255e-03,
-5.0175013429732242e-03,
-1.2547236272489583e-02,
-9.4554964308388671e-03,
-1.2526086181600525e-02,
-7.9170697760437767e-03,
-7.9955735204175690e-03,
-9.3559433444469070e-03,
6.5468611982999303e-02,
3.9572875517634137e-01,
7.5244022299886659e-01,
9.6376478690218559e-01,
9.9854433855162328e-01,
9.9992977025287921e-01,
9.9939086751140449e-01,
9.9994372267071396e-01,
9.9939121813418674e-01,
9.9911237310424483e-01,
9.6019584878271580e-01,
6.3186279338432438e-01,
2.5797401028763473e-01,
9.4014888527335638e-03,
-3.0798345608649747e-03,
-4.5230367033685034e-03,
-6.8933410388274038e-03,
-9.0352195539015398e-03,
-8.5913667165340209e-03,
-8.3690869120289398e-03,
-7.8685832338754313e-03,
-8.3657578711085132e-06,
5.4301225442817177e-03,
-2.7745589759259194e-03,
];
pub const BLUE: [f64; 32] = [
9.9209771469720676e-01,
9.8876426059369127e-01,
9.9539040744505636e-01,
9.9529317353008218e-01,
9.9181447411633950e-01,
1.0002584039673432e+00,
9.9968478437342512e-01,
9.9988120766657174e-01,
9.8504012146370434e-01,
7.9029849053031276e-01,
5.6082198617463974e-01,
3.3133458513996528e-01,
1.3692410840839175e-01,
1.8914906559664151e-02,
-5.1129770932550889e-06,
-4.2395493167891873e-04,
-4.1934593101534273e-04,
1.7473028136486615e-03,
3.7999160177631316e-03,
-5.5101474906588642e-04,
-4.3716662898480967e-05,
7.5874501748732798e-03,
2.5795650780554021e-02,
3.8168376532500548e-02,
4.9489586408030833e-02,
4.9595992290102905e-02,
4.9814819505812249e-02,
3.9840911064978023e-02,
3.0501024937233868e-02,
2.1243054765241080e-02,
6.9596532104356399e-03,
4.1733649330980525e-03,
];
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn intensity_at_wavelength_returns_expected_value_at_minimum_wavelength() {
let target = Spectrum {
shortest_wavelength: 400.5,
longest_wavelength: 700.25,
samples: vec![0.5, 1.0, 0.75, 1.5],
};
assert!(target.intensity_at_wavelength(400.5) == 0.5)
}
#[test]
fn intensity_at_wavelength_returns_expected_value_at_max_wavelength() {
let target = Spectrum {
shortest_wavelength: 400.5,
longest_wavelength: 700.25,
samples: vec![0.5, 1.0, 0.75, 1.5],
};
assert!(target.intensity_at_wavelength(700.25) == 1.5)
}
#[test]
fn intensity_at_wavelength_returns_expected_value_at_interior_sample_wavelength() {
let target = Spectrum {
shortest_wavelength: 400.0,
longest_wavelength: 700.0,
samples: vec![0.5, 1.0, 0.75, 1.5],
};
assert!(target.intensity_at_wavelength(500.0) == 1.0);
assert!(target.intensity_at_wavelength(600.0) == 0.75);
}
#[test]
fn intensity_at_wavelength_returns_expected_value_at_halfway_between_sample_wavelength() {
let target = Spectrum {
shortest_wavelength: 400.0,
longest_wavelength: 700.0,
samples: vec![0.5, 1.0, 0.75, 1.5],
};
assert!(target.intensity_at_wavelength(450.0) == 0.75);
assert!(target.intensity_at_wavelength(550.0) == 0.875);
assert!(target.intensity_at_wavelength(650.0) == 1.125);
}
#[test]
fn intensity_below_minimum_wavelength_is_zero() {
let target = Spectrum {
shortest_wavelength: 400.0,
longest_wavelength: 700.0,
samples: vec![0.5, 1.0, 0.75, 1.5],
};
assert!(target.intensity_at_wavelength(399.9999) == 0.0);
}
#[test]
fn intensity_above_maximum_wavelength_is_zero() {
let target = Spectrum {
shortest_wavelength: 400.0,
longest_wavelength: 700.0,
samples: vec![0.5, 1.0, 0.75, 1.5],
};
assert!(target.intensity_at_wavelength(700.0001) == 0.0);
}
}

222
src/image.rs
View File

@ -1,105 +1,123 @@
use std::fs::File; use std::convert::TryInto;
use std::io::BufWriter;
use std::path::Path;
use crate::colour::{ColourRgbF, ColourRgbU8, ColourXyz}; use nalgebra::{clamp, convert, RealField, Vector3};
use crate::util::Array2D;
use super::colour::{ColourRgbF, ColourRgbU8};
#[derive(Debug)]
pub struct ImageRgbU8 { pub struct ImageRgbU8 {
data: Array2D<[u8; 3]>, pixel_data: Vec<u8>,
width: u32,
height: u32,
} }
impl ImageRgbU8 { impl ImageRgbU8 {
pub fn new(width: usize, height: usize) -> ImageRgbU8 { pub fn new(width: u32, height: u32) -> ImageRgbU8 {
ImageRgbU8 { ImageRgbU8 {
data: Array2D::new(height, width), width: width,
height: height,
pixel_data: vec![0; (width * height * 3) as usize],
} }
} }
pub fn get_colour(&self, row: usize, column: usize) -> ColourRgbU8 { pub fn clear(&mut self) -> &mut ImageRgbU8 {
for byte in self.pixel_data.iter_mut() {
*byte = 0u8;
}
self
}
pub fn get_colour(&self, row: u32, column: u32) -> ColourRgbU8 {
assert!(row < self.height && column < self.width);
let index = self.calculate_index(row, column);
ColourRgbU8 { ColourRgbU8 {
values: self.data[row][column], values: self.pixel_data[index..index + 3]
.try_into()
.expect("Wrong length."),
} }
} }
pub fn set_colour(&mut self, row: usize, column: usize, colour: ColourRgbU8) { pub fn set_colour(&mut self, row: u32, column: u32, colour: ColourRgbU8) {
let slice = &mut self.data[row][column]; assert!(row < self.height && column < self.width);
slice.copy_from_slice(&colour.values[..]); let index = self.calculate_index(row, column);
self.pixel_data[index..index + 3].copy_from_slice(&colour.values[..]);
} }
pub fn get_pixel_data(&self) -> &[u8] { pub fn get_pixel_data(&self) -> &Vec<u8> {
let data = self.data.as_slice(); &self.pixel_data
unsafe { std::slice::from_raw_parts(data[0].as_ptr(), data.len() * 3) }
} }
pub fn get_width(&self) -> usize { pub fn get_width(&self) -> u32 {
self.data.get_width() self.width
} }
pub fn get_height(&self) -> usize { pub fn get_height(&self) -> u32 {
self.data.get_height() self.height
} }
pub fn num_channels() -> usize { pub fn num_channels() -> u32 {
3 3
} }
pub fn update(&mut self, start_row: usize, start_column: usize, image: &ImageRgbU8) { fn calculate_index(&self, row: u32, column: u32) -> usize {
self.data.update_block(start_row, start_column, &image.data); assert!(row < self.height && column < self.width);
} (((self.height - (row + 1)) * self.width + column) * Self::num_channels()) as usize
pub fn write_png(&self, filename: &Path) -> Result<(), std::io::Error> {
let file = File::create(filename)?;
let file_buffer = &mut BufWriter::new(file);
let mut encoder = png::Encoder::new(
file_buffer,
self.get_width() as u32,
self.get_height() as u32,
);
encoder.set_color(png::ColorType::RGB);
encoder.set_depth(png::BitDepth::Eight);
let mut writer = encoder.write_header()?;
writer.write_image_data(self.get_pixel_data())?;
Ok(())
} }
} }
pub struct ImageRgbF { pub struct ImageRgbF<T: RealField> {
pub data: Array2D<ColourRgbF>, pixel_data: Vec<T>,
width: u32,
height: u32,
} }
impl ImageRgbF { impl<T: RealField> ImageRgbF<T> {
pub fn new(width: usize, height: usize) -> ImageRgbF { pub fn new(width: u32, height: u32) -> ImageRgbF<T> {
ImageRgbF { ImageRgbF {
data: Array2D::new(height, width), width: width,
height: height,
pixel_data: vec![convert(0.0); (width * height * 3) as usize],
} }
} }
pub fn clear(&mut self) { pub fn clear(&mut self) -> &mut ImageRgbF<T> {
self.data.clear(); for elem in self.pixel_data.iter_mut() {
*elem = T::zero();
}
self
} }
pub fn get_colour(&self, row: usize, column: usize) -> ColourRgbF { pub fn get_colour(&self, row: u32, column: u32) -> ColourRgbF<T> {
self.data[row][column] assert!(row < self.height && column < self.width);
let index = self.calculate_index(row, column);
ColourRgbF::from_vector3(&Vector3::from_row_slice(&self.pixel_data[index..index + 3]))
} }
pub fn set_colour(&mut self, row: usize, column: usize, colour: ColourRgbF) { pub fn set_colour(&mut self, row: u32, column: u32, colour: ColourRgbF<T>) {
self.data[row][column] = colour; assert!(row < self.height && column < self.width);
let index = self.calculate_index(row, column);
self.pixel_data[index..index + 3].copy_from_slice(&colour.as_vector3().as_slice());
} }
pub fn get_width(&self) -> usize { pub fn get_pixel_data(&self) -> &Vec<T> {
self.data.get_width() &self.pixel_data
} }
pub fn get_height(&self) -> usize { pub fn get_width(&self) -> u32 {
self.data.get_height() self.width
} }
pub fn num_channels() -> usize { pub fn get_height(&self) -> u32 {
self.height
}
pub fn num_channels() -> u32 {
3 3
} }
fn calculate_index(&self, row: u32, column: u32) -> usize {
assert!(row < self.height && column < self.width);
(((self.height - (row + 1)) * self.width + column) * Self::num_channels()) as usize
}
} }
pub trait NormalizedAsByte { pub trait NormalizedAsByte {
@ -109,67 +127,47 @@ pub trait NormalizedAsByte {
impl NormalizedAsByte for f32 { impl NormalizedAsByte for f32 {
fn normalized_to_byte(self) -> u8 { fn normalized_to_byte(self) -> u8 {
(self * (u8::MAX as f32)) as u8 (self * (std::u8::MAX as f32)) as u8
} }
fn byte_to_normalized(byte: u8) -> f32 { fn byte_to_normalized(byte: u8) -> f32 {
(byte as f32) / (u8::MAX as f32) (byte as f32) / (std::u8::MAX as f32)
} }
} }
impl NormalizedAsByte for f64 { impl NormalizedAsByte for f64 {
fn normalized_to_byte(self) -> u8 { fn normalized_to_byte(self) -> u8 {
(self * (u8::MAX as f64)) as u8 (self * (std::u8::MAX as f64)) as u8
} }
fn byte_to_normalized(byte: u8) -> f64 { fn byte_to_normalized(byte: u8) -> f64 {
(byte as f64) / (u8::MAX as f64) (byte as f64) / (std::u8::MAX as f64)
} }
} }
pub trait ToneMapper<SourceType> { pub trait ToneMapper<T: RealField> {
fn apply_tone_mapping(&self, image_in: &Array2D<SourceType>, image_out: &mut ImageRgbU8); fn apply_tone_mapping(&self, image_in: &ImageRgbF<T>, image_out: &mut ImageRgbU8);
} }
#[derive(Default)]
pub struct ClampingToneMapper {} pub struct ClampingToneMapper {}
impl ClampingToneMapper { impl ClampingToneMapper {
fn clamp(v: &f64) -> u8 { pub fn new() -> ClampingToneMapper {
v.clamp(0.0, 1.0).normalized_to_byte() ClampingToneMapper {}
}
fn clamp<T: RealField + NormalizedAsByte>(v: &T) -> u8 {
clamp(v, &T::zero(), &T::one()).normalized_to_byte()
} }
} }
impl ToneMapper<ColourRgbF> for ClampingToneMapper { impl<T: RealField + NormalizedAsByte> ToneMapper<T> for ClampingToneMapper {
fn apply_tone_mapping(&self, image_in: &Array2D<ColourRgbF>, image_out: &mut ImageRgbU8) { fn apply_tone_mapping(&self, image_in: &ImageRgbF<T>, image_out: &mut ImageRgbU8) {
assert!(image_in.get_width() == image_out.get_width()); assert!(image_in.get_width() == image_out.get_width());
assert!(image_in.get_height() == image_out.get_height()); assert!(image_in.get_height() == image_out.get_height());
for column in 0..image_in.get_width() { for column in 0..image_in.get_width() {
for row in 0..image_in.get_height() { for row in 0..image_in.get_height() {
let colour = image_in[row][column]; let colour = image_in.get_colour(row, column);
image_out.set_colour(
row,
column,
ColourRgbU8 {
values: [
Self::clamp(&colour.red()),
Self::clamp(&colour.green()),
Self::clamp(&colour.blue()),
],
},
);
}
}
}
}
impl ToneMapper<ColourXyz> for ClampingToneMapper {
fn apply_tone_mapping(&self, image_in: &Array2D<ColourXyz>, image_out: &mut ImageRgbU8) {
assert!(image_in.get_width() == image_out.get_width());
assert!(image_in.get_height() == image_out.get_height());
for column in 0..image_in.get_width() {
for row in 0..image_in.get_height() {
let colour = image_in[row][column].to_srgb();
image_out.set_colour( image_out.set_colour(
row, row,
column, column,
@ -190,30 +188,6 @@ impl ToneMapper<ColourXyz> for ClampingToneMapper {
mod tests { mod tests {
use super::*; use super::*;
#[test]
fn get_pixel_data_returns_correct_values() {
let mut target = ImageRgbU8::new(4, 3);
for i in 0..3 {
for j in 0..4 {
target.set_colour(
i,
j,
ColourRgbU8 {
values: [i as u8, j as u8, i as u8],
},
)
}
}
for i in 0..3 {
for j in 0..4 {
let index = (i * 4 + j) * 3;
assert!(target.get_pixel_data()[index] == i as u8);
assert!(target.get_pixel_data()[index + 1] == j as u8);
assert!(target.get_pixel_data()[index + 2] == i as u8);
}
}
}
mod normalized_as_byte { mod normalized_as_byte {
use super::*; use super::*;
@ -287,7 +261,7 @@ mod tests {
let mut image_in = ImageRgbF::new(1, 1); let mut image_in = ImageRgbF::new(1, 1);
let mut image_out = ImageRgbU8::new(1, 1); let mut image_out = ImageRgbU8::new(1, 1);
image_in.set_colour(0, 0, ColourRgbF::new(0.0, 0.0, 0.0)); image_in.set_colour(0, 0, ColourRgbF::new(0.0, 0.0, 0.0));
target.apply_tone_mapping(&image_in.data, &mut image_out); target.apply_tone_mapping(&image_in, &mut image_out);
assert!(image_out.get_colour(0, 0).values == [0, 0, 0]); assert!(image_out.get_colour(0, 0).values == [0, 0, 0]);
} }
@ -297,8 +271,8 @@ mod tests {
let mut image_in = ImageRgbF::new(1, 1); let mut image_in = ImageRgbF::new(1, 1);
let mut image_out = ImageRgbU8::new(1, 1); let mut image_out = ImageRgbU8::new(1, 1);
image_in.set_colour(0, 0, ColourRgbF::new(1.0, 1.0, 1.0)); image_in.set_colour(0, 0, ColourRgbF::new(1.0, 1.0, 1.0));
target.apply_tone_mapping(&image_in.data, &mut image_out); target.apply_tone_mapping(&image_in, &mut image_out);
assert!(image_out.get_colour(0, 0).values == [0xff, 0xff, 0xff]); assert!(image_out.get_colour(0, 0).values == [0xff, 0xff, 0xff]);
} }
#[test] #[test]
@ -307,8 +281,8 @@ mod tests {
let mut image_in = ImageRgbF::new(1, 1); let mut image_in = ImageRgbF::new(1, 1);
let mut image_out = ImageRgbU8::new(1, 1); let mut image_out = ImageRgbU8::new(1, 1);
image_in.set_colour(0, 0, ColourRgbF::new(2.0, 2.0, 2.0)); image_in.set_colour(0, 0, ColourRgbF::new(2.0, 2.0, 2.0));
target.apply_tone_mapping(&image_in.data, &mut image_out); target.apply_tone_mapping(&image_in, &mut image_out);
assert!(image_out.get_colour(0, 0).values == [0xff, 0xff, 0xff]); assert!(image_out.get_colour(0, 0).values == [0xff, 0xff, 0xff]);
} }
#[test] #[test]
@ -317,8 +291,8 @@ mod tests {
let mut image_in = ImageRgbF::new(1, 1); let mut image_in = ImageRgbF::new(1, 1);
let mut image_out = ImageRgbU8::new(1, 1); let mut image_out = ImageRgbU8::new(1, 1);
image_in.set_colour(0, 0, ColourRgbF::new(0.0, 2.0, 0.0)); image_in.set_colour(0, 0, ColourRgbF::new(0.0, 2.0, 0.0));
target.apply_tone_mapping(&image_in.data, &mut image_out); target.apply_tone_mapping(&image_in, &mut image_out);
assert!(image_out.get_colour(0, 0).values == [0x0, 0xff, 0x0]); assert!(image_out.get_colour(0, 0).values == [0x0, 0xff, 0x0]);
} }
#[test] #[test]
@ -327,8 +301,8 @@ mod tests {
let mut image_in = ImageRgbF::new(1, 1); let mut image_in = ImageRgbF::new(1, 1);
let mut image_out = ImageRgbU8::new(1, 1); let mut image_out = ImageRgbU8::new(1, 1);
image_in.set_colour(0, 0, ColourRgbF::new(0.5, 0.0, 0.0)); image_in.set_colour(0, 0, ColourRgbF::new(0.5, 0.0, 0.0));
target.apply_tone_mapping(&image_in.data, &mut image_out); target.apply_tone_mapping(&image_in, &mut image_out);
assert!(image_out.get_colour(0, 0).values == [0x7f, 0x0, 0x0]); assert!(image_out.get_colour(0, 0).values == [0x7f, 0x0, 0x0]);
} }
} }
} }

47
src/integrators.rs Normal file
View File

@ -0,0 +1,47 @@
use nalgebra::{convert, RealField, Vector3};
use super::algebra_utils::try_change_of_basis_matrix;
use super::colour::ColourRgbF;
use super::raycasting::{IntersectionInfo, Ray};
use super::sampler::Sampler;
pub trait Integrator<T: RealField> {
fn integrate(&self, sampler: &Sampler<T>, info: &IntersectionInfo<T>) -> ColourRgbF<T>;
}
pub struct DirectionalLight<T: RealField> {
pub direction: Vector3<T>,
pub colour: ColourRgbF<T>,
}
pub struct WhittedIntegrator<T: RealField> {
pub ambient_light: ColourRgbF<T>,
pub lights: Vec<DirectionalLight<T>>,
}
// TODO: Get rid of the magic bias number, which should be calculated base on expected error
// bounds and tangent direction
impl<T: RealField> Integrator<T> for WhittedIntegrator<T> {
fn integrate(&self, sampler: &Sampler<T>, info: &IntersectionInfo<T>) -> ColourRgbF<T> {
self.lights
.iter()
.map(|light| {
let basis_change =
try_change_of_basis_matrix(&info.tangent, &info.cotangent, &info.normal)
.expect("Normal, tangent and cotangent don't for a valid basis.");
match sampler
.sample(&Ray::new(info.location, light.direction).bias(convert(0.000000001)))
{
Some(_) => self.ambient_light,
None => {
info.material.bsdf()(
basis_change * info.retro,
basis_change * light.direction,
light.colour,
) * light.direction.dot(&info.normal).abs()
}
}
})
.fold(self.ambient_light, |a, b| a + b)
}
}

19
src/integrators/mod.rs
View File

@ -1,19 +0,0 @@
use super::colour::Photon;
use super::raycasting::IntersectionInfo;
use super::sampler::Sampler;
mod whitted_integrator;
pub use whitted_integrator::*;
mod simple_random_integrator;
pub use simple_random_integrator::*;
pub trait Integrator {
fn integrate(
&self,
sampler: &Sampler,
info: &IntersectionInfo,
photon: &Photon,
recursion_limit: u16,
) -> Photon;
}

65
src/integrators/simple_random_integrator.rs
View File

@ -1,65 +0,0 @@
use crate::colour::{ColourRgbF, Photon, Spectrum};
use crate::materials::MaterialSampleResult;
use crate::math::Vec3;
use crate::raycasting::{IntersectionInfo, Ray};
use crate::sampler::Sampler;
use crate::util::algebra_utils::try_change_of_basis_matrix;
use super::Integrator;
pub struct SimpleRandomIntegrator {}
impl Integrator for SimpleRandomIntegrator {
fn integrate(
&self,
sampler: &Sampler,
info: &IntersectionInfo,
photon: &Photon,
recursion_limit: u16,
) -> Photon {
if recursion_limit == 0 {
return Photon {
wavelength: 0.0,
intensity: 0.0,
};
}
let world_to_bsdf_space =
try_change_of_basis_matrix(&info.tangent, &info.cotangent, &info.normal)
.expect("Normal, tangent and cotangent don't form a valid basis.");
let bsdf_to_world_space = world_to_bsdf_space
.try_inverse()
.expect("Expected matrix to be invertable.");
let world_space_w_i = info.retro;
let w_i = world_to_bsdf_space * world_space_w_i;
let MaterialSampleResult {
direction: w_o,
pdf: w_o_pdf,
} = info.material.sample(&w_i, photon);
let world_space_w_o = bsdf_to_world_space * w_o;
info.material.bsdf()(
&w_o,
&w_i,
&match sampler.sample(&Ray::new(info.location, world_space_w_o).bias(0.000_000_1)) {
None => photon.set_intensity(test_lighting_environment(
&world_space_w_o,
photon.wavelength,
)),
Some(recursive_hit) => {
self.integrate(sampler, &recursive_hit, photon, recursion_limit - 1)
}
}
.scale_intensity(w_o_pdf)
.scale_intensity(world_space_w_o.dot(&info.normal).abs()),
)
}
}
pub fn test_lighting_environment(w_o: &Vec3<f64>, wavelength: f64) -> f64 {
//let sun_direction = Vec3::new(1.0, 1.0, -1.0).normalize();
//if w_o.dot(&sun_direction) >= 0.99 {
// 300.0
//} else {
let sky_colour = ColourRgbF::new(w_o.y(), w_o.y(), 1.0);
Spectrum::reflection_from_linear_rgb(&sky_colour).intensity_at_wavelength(wavelength)
//}
}

87
src/integrators/whitted_integrator.rs
View File

@ -1,87 +0,0 @@
use crate::colour::{Photon, Spectrum};
use crate::materials::MaterialSampleResult;
use crate::math::Vec3;
use crate::raycasting::{IntersectionInfo, Ray};
use crate::sampler::Sampler;
use crate::util::algebra_utils::try_change_of_basis_matrix;
use super::Integrator;
pub struct DirectionalLight {
pub direction: Vec3<f64>,
pub spectrum: Spectrum,
}
pub struct WhittedIntegrator {
pub ambient_light: Spectrum,
pub lights: Vec<DirectionalLight>,
}
impl Integrator for WhittedIntegrator {
fn integrate(
&self,
sampler: &Sampler,
info: &IntersectionInfo,
photon: &Photon,
recursion_limit: u16,
) -> Photon {
let world_to_bsdf_space =
try_change_of_basis_matrix(&info.tangent, &info.cotangent, &info.normal)
.expect("Normal, tangent and cotangent don't for a valid basis.");
let bsdf_to_world_space = world_to_bsdf_space
.try_inverse()
.expect("Expected matrix to be invertable.");
self.lights
.iter()
.map(|light| {
match sampler.sample(&Ray::new(info.location, light.direction).bias(0.000_000_1)) {
Some(_) => self.ambient_light.emit_photon(photon),
None => info.material.bsdf()(
&(world_to_bsdf_space * info.retro),
&(world_to_bsdf_space * light.direction),
&light
.spectrum
.emit_photon(photon)
.scale_intensity(light.direction.dot(&info.normal).abs()),
),
}
})
.chain(
[info
.material
.sample(&(world_to_bsdf_space * info.retro), photon)]
.iter()
.map(|MaterialSampleResult { direction, pdf: _ }| {
let world_space_direction = bsdf_to_world_space * direction;
match sampler
.sample(&Ray::new(info.location, world_space_direction).bias(0.000_000_1))
{
Some(recursive_hit) => {
if recursion_limit > 0 {
let photon = info.material.bsdf()(
&(world_to_bsdf_space * info.retro),
direction,
&self.integrate(
sampler,
&recursive_hit,
photon,
recursion_limit - 1,
),
);
photon
.scale_intensity(world_space_direction.dot(&info.normal).abs())
} else {
photon.scale_intensity(0.0)
}
}
None => photon.scale_intensity(0.0),
}
}),
)
.fold(photon.clone(), |a, b| {
let mut result = a;
result.intensity += b.intensity;
result
})
}
}

13
src/lib.rs
View File

@ -1,18 +1,9 @@
#![doc = include_str!("../README.md")] pub mod algebra_utils;
pub mod camera;
pub mod accumulation_buffer;
mod camera;
pub mod colour; pub mod colour;
pub mod image; pub mod image;
pub mod integrators; pub mod integrators;
pub mod materials; pub mod materials;
pub mod math;
pub mod mesh;
pub mod random_distributions;
pub mod raycasting; pub mod raycasting;
pub mod realtype;
pub mod sampler; pub mod sampler;
pub mod scene; pub mod scene;
pub mod util;
pub use camera::partial_render_scene;

268
src/main.rs
View File

@ -1,92 +1,35 @@
use rayon::prelude::*;
use sdl2::event::Event; use sdl2::event::Event;
use sdl2::keyboard::Keycode; use sdl2::keyboard::Keycode;
use sdl2::pixels::PixelFormatEnum; use sdl2::pixels::PixelFormatEnum;
use sdl2::rect::Rect;
use sdl2::render::{Canvas, Texture}; use sdl2::render::{Canvas, Texture};
use sdl2::Sdl; use sdl2::Sdl;
use clap::Arg;
use std::path::{Path, PathBuf};
use std::sync::{mpsc, Arc};
use std::time::Duration; use std::time::Duration;
use vanrijn::accumulation_buffer::AccumulationBuffer; use nalgebra::Vector3;
use vanrijn::colour::{ColourRgbF, NamedColour, Spectrum};
use vanrijn::image::{ClampingToneMapper, ImageRgbU8}; use std::cmp::min;
use vanrijn::materials::LambertianMaterial; use std::rc::Rc;
use vanrijn::math::Vec3;
use vanrijn::mesh::load_obj; use vanrijn::camera::partial_render_scene;
use vanrijn::partial_render_scene; use vanrijn::colour::{ColourRgbF, NamedColour};
use vanrijn::raycasting::{BoundingVolumeHierarchy, Plane, Primitive, Sphere}; use vanrijn::image::{ClampingToneMapper, ImageRgbF, ImageRgbU8, ToneMapper};
use vanrijn::materials::{LambertianMaterial, PhongMaterial};
use vanrijn::raycasting::{Plane, Sphere};
use vanrijn::scene::Scene; use vanrijn::scene::Scene;
use vanrijn::util::TileIterator;
#[derive(Debug)]
struct CommandLineParameters {
width: usize,
height: usize,
output_file: Option<PathBuf>,
time: f64,
}
fn parse_args() -> CommandLineParameters {
let matches = clap::App::new("vanrijn")
.version("alpha")
.author("Matthew Gordon <matthew@gordon.earth")
.arg(
Arg::with_name("size")
.long("size")
.value_name("SIZE")
.help("The width and height of the output image, in pixels.")
.takes_value(true)
.number_of_values(2)
.required(true),
)
.arg(
Arg::with_name("output_png")
.long("out")
.value_name("FILENAME")
.help("Filename for output PNG.")
.takes_value(true)
.required(false),
)
.arg(
Arg::with_name("time")
.long("time")
.value_name("SECONDS")
.takes_value(true)
.default_value("0"),
)
.get_matches();
let mut size_iter = matches.values_of("size").unwrap();
let width = size_iter.next().unwrap().parse().unwrap();
let height = size_iter.next().unwrap().parse().unwrap();
let output_file = matches.value_of_os("output_png").map(PathBuf::from);
let time = matches.value_of("time").unwrap().parse().unwrap();
CommandLineParameters {
width,
height,
output_file,
time,
}
}
fn update_texture(image: &ImageRgbU8, texture: &mut Texture) { fn update_texture(image: &ImageRgbU8, texture: &mut Texture) {
texture texture
.update( .update(
Rect::new(0, 0, image.get_width() as u32, image.get_height() as u32), None,
image.get_pixel_data(), image.get_pixel_data().as_slice(),
image.get_width() * ImageRgbU8::num_channels(), (image.get_width() * ImageRgbU8::num_channels()) as usize,
) )
.expect("Couldn't update texture."); .expect("Couldn't update texture.");
} }
fn init_canvas( fn init_canvas(
image_width: usize, image_width: u32,
image_height: usize, image_height: u32,
) -> Result<(Sdl, Canvas<sdl2::video::Window>), Box<dyn std::error::Error>> { ) -> Result<(Sdl, Canvas<sdl2::video::Window>), Box<dyn std::error::Error>> {
let sdl_context = sdl2::init()?; let sdl_context = sdl2::init()?;
let video_subsystem = sdl_context.video()?; let video_subsystem = sdl_context.video()?;
@ -102,11 +45,8 @@ fn init_canvas(
} }
pub fn main() -> Result<(), Box<dyn std::error::Error>> { pub fn main() -> Result<(), Box<dyn std::error::Error>> {
let parameters = parse_args(); let image_width = 1200;
let image_width = parameters.width; let image_height = 900;
let image_height = parameters.height;
let mut rendered_image = AccumulationBuffer::new(image_width, image_height);
let (sdl_context, mut canvas) = init_canvas(image_width, image_height)?; let (sdl_context, mut canvas) = init_canvas(image_width, image_height)?;
@ -116,118 +56,88 @@ pub fn main() -> Result<(), Box<dyn std::error::Error>> {
image_width as u32, image_width as u32,
image_height as u32, image_height as u32,
)?; )?;
let mut output_image = ImageRgbF::<f64>::new(image_width, image_height);
let model_file_path =
Path::new(env!("CARGO_MANIFEST_DIR")).join("test_data/stanford_bunny.obj");
println!("Loading object...");
let mut model_object = load_obj(
&model_file_path,
Arc::new(LambertianMaterial {
colour: Spectrum::reflection_from_linear_rgb(&ColourRgbF::from_named(
NamedColour::Yellow,
)),
diffuse_strength: 0.05,
//reflection_strength: 0.9,
}),
)?;
println!("Building BVH...");
let model_bvh: Box<dyn Primitive> =
Box::new(BoundingVolumeHierarchy::build(model_object.as_mut_slice()));
println!("Constructing Scene...");
let scene = Scene { let scene = Scene {
camera_location: Vec3::new(-2.0, 1.0, -5.0), camera_location: Vector3::new(0.0, 0.0, 0.0),
objects: vec![ objects: vec![
Box::new(vec![ Box::new(Plane::new(
Box::new(Plane::new( Vector3::new(0.0, 1.0, 0.0),
Vec3::new(0.0, 1.0, 0.0), -2.0,
-2.0, Rc::new(LambertianMaterial {
Arc::new(LambertianMaterial { colour: ColourRgbF::new(0.55, 0.27, 0.04),
colour: Spectrum::reflection_from_linear_rgb(&ColourRgbF::new( diffuse_strength: 0.1,
0.55, 0.27, 0.04, }),
)), )),
diffuse_strength: 0.1, Box::new(Sphere::new(
}), Vector3::new(1.25, -0.5, 6.0),
)) as Box<dyn Primitive>, 1.0,
Box::new(Sphere::new( Rc::new(LambertianMaterial {
Vec3::new(-6.25, -0.5, 1.0), colour: ColourRgbF::from_named(NamedColour::Green),
1.0, diffuse_strength: 0.1,
Arc::new(LambertianMaterial { }),
colour: Spectrum::reflection_from_linear_rgb(&ColourRgbF::from_named( )),
NamedColour::Green, Box::new(Sphere::new(
)), Vector3::new(-1.25, -0.5, 6.0),
diffuse_strength: 0.1, 1.0,
}), Rc::new(LambertianMaterial {
)), colour: ColourRgbF::from_named(NamedColour::Blue),
Box::new(Sphere::new( diffuse_strength: 0.1,
Vec3::new(-4.25, -0.5, 2.0), }),
1.0, )),
Arc::new(LambertianMaterial { Box::new(Sphere::new(
colour: Spectrum::reflection_from_linear_rgb(&ColourRgbF::from_named( Vector3::new(0.0, 1.5, 6.0),
NamedColour::Blue, 1.0,
)), Rc::new(PhongMaterial {
diffuse_strength: 0.1, colour: ColourRgbF::from_named(NamedColour::Red),
// diffuse_strength: 0.01, diffuse_strength: 0.05,
// reflection_strength: 0.99, smoothness: 20.0,
}), specular_strength: 250.0,
)), }),
Box::new(Sphere::new( )),
Vec3::new(-5.0, 1.5, 1.0),
1.0,
Arc::new(LambertianMaterial {
colour: Spectrum::reflection_from_linear_rgb(&ColourRgbF::from_named(
NamedColour::Red,
)),
diffuse_strength: 0.05,
//smoothness: 100.0,
//specular_strength: 1.0,
}),
)),
]) as Box<dyn Primitive>,
model_bvh,
], ],
}; };
println!("Done.");
let mut event_pump = sdl_context.event_pump()?; let mut event_pump = sdl_context.event_pump()?;
let mut i = 0;
let (tile_tx, tile_rx) = mpsc::channel();
let mut tile_rx = Some(tile_rx);
let worker_boss = std::thread::spawn(move || {
let end_tx = tile_tx.clone();
TileIterator::new(image_width as usize, image_height as usize, 2048)
.cycle()
.map(move |tile| (tile, tile_tx.clone()))
.par_bridge()
.try_for_each(|(tile, tx)| {
let rendered_tile = partial_render_scene(&scene, tile, image_height, image_width);
// There's nothing we can do if this fails, and we're already
// at the end of the function anyway, so just ignore result.
tx.send(Some((tile, rendered_tile))).ok()
});
end_tx.send(None).ok();
});
'running: loop { 'running: loop {
if let Some(ref tile_rx) = tile_rx { let tile_size = 128;
for message in tile_rx.try_iter() { for tile_row in 0..1 + (output_image.get_height() + 1) / tile_size {
if let Some((tile, tile_accumulation_buffer)) = message { for tile_column in 0..1 + (output_image.get_width() + 1) / tile_size {
rendered_image.merge_tile(&tile, &tile_accumulation_buffer); let row_start = tile_row * tile_size;
let rgb_image = rendered_image.to_image_rgb_u8(&ClampingToneMapper {}); let row_end = min(tile_row * tile_size + tile_size, output_image.get_height());
update_texture(&rgb_image, &mut rendered_image_texture); let column_start = tile_column * tile_size;
canvas.copy(&rendered_image_texture, None, None).unwrap(); let column_end = min(
canvas.present(); tile_column * tile_size + tile_size,
} else if let Some(image_filename) = parameters.output_file { output_image.get_width(),
rendered_image );
.to_image_rgb_u8(&ClampingToneMapper {}) partial_render_scene(
.write_png(&image_filename)?; &mut output_image,
break 'running; &scene,
row_start,
row_end,
column_start,
column_end,
);
let mut output_image_rgbu8 = ImageRgbU8::new(image_width, image_height);
ClampingToneMapper {}.apply_tone_mapping(&output_image, &mut output_image_rgbu8);
update_texture(&output_image_rgbu8, &mut rendered_image_texture);
canvas.copy(&rendered_image_texture, None, None)?;
canvas.present();
for event in event_pump.poll_iter() {
match event {
Event::Quit { .. }
| Event::KeyDown {
keycode: Some(Keycode::Escape),
..
} => break 'running,
_ => {}
}
} }
} }
} }
i = (i + 1) % 255;
for event in event_pump.poll_iter() { for event in event_pump.poll_iter() {
match event { match event {
Event::Quit { .. } Event::Quit { .. }
@ -241,7 +151,5 @@ pub fn main() -> Result<(), Box<dyn std::error::Error>> {
::std::thread::sleep(Duration::new(0, 1_000_000_000u32 / 60)); ::std::thread::sleep(Duration::new(0, 1_000_000_000u32 / 60));
} }
drop(tile_rx.take());
worker_boss.join().expect("Couldn't join worker threads.");
Ok(()) Ok(())
} }

66
src/materials.rs Normal file
View File

@ -0,0 +1,66 @@
use nalgebra::{RealField, Vector3};
use super::colour::{ColourRgbF, NamedColour};
use std::fmt::Debug;
pub trait Material<T: RealField>: Debug {
fn bsdf<'a>(
&'a self,
) -> Box<dyn Fn(Vector3<T>, Vector3<T>, ColourRgbF<T>) -> ColourRgbF<T> + 'a>;
}
#[derive(Debug)]
pub struct LambertianMaterial<T: RealField> {
pub colour: ColourRgbF<T>,
pub diffuse_strength: T,
}
impl<T: RealField> LambertianMaterial<T> {
pub fn new_dummy() -> LambertianMaterial<T> {
LambertianMaterial {
colour: ColourRgbF::new(T::one(), T::one(), T::one()),
diffuse_strength: T::one(),
}
}
}
impl<T: RealField> Material<T> for LambertianMaterial<T> {
fn bsdf<'a>(
&'a self,
) -> Box<dyn Fn(Vector3<T>, Vector3<T>, ColourRgbF<T>) -> ColourRgbF<T> + 'a> {
Box::new(
move |_w_o: Vector3<T>, _w_i: Vector3<T>, colour_in: ColourRgbF<T>| {
self.colour * colour_in * self.diffuse_strength
},
)
}
}
#[derive(Debug)]
pub struct PhongMaterial<T: RealField> {
pub colour: ColourRgbF<T>,
pub diffuse_strength: T,
pub specular_strength: T,
pub smoothness: T,
}
impl<T: RealField> Material<T> for PhongMaterial<T> {
fn bsdf<'a>(
&'a self,
) -> Box<dyn Fn(Vector3<T>, Vector3<T>, ColourRgbF<T>) -> ColourRgbF<T> + 'a> {
Box::new(
move |w_o: Vector3<T>, w_i: Vector3<T>, colour_in: ColourRgbF<T>| {
if w_i.z < T::zero() || w_o.z < T::zero() {
ColourRgbF::from_vector3(&Vector3::zeros())
} else {
let reflection_vector = Vector3::new(-w_i.x, -w_i.y, w_i.z);
self.colour * colour_in * self.diffuse_strength
+ ColourRgbF::from_named(NamedColour::White)
* w_o.dot(&reflection_vector).abs().powf(self.smoothness)
* (self.specular_strength / w_i.dot(&Vector3::z_axis()))
}
},
)
}
}

60
src/materials/lambertian_material.rs
View File

@ -1,60 +0,0 @@
use crate::colour::{Photon, Spectrum};
use crate::math::Vec3;
use super::{Material, MaterialSampleResult};
use rand::distributions::Open01;
use rand::{thread_rng, Rng};
use std::f64::consts::PI;
use std::fmt::Debug;
#[derive(Debug)]
pub struct LambertianMaterial {
pub colour: Spectrum,
pub diffuse_strength: f64,
}
impl LambertianMaterial {
pub fn new_dummy() -> LambertianMaterial {
LambertianMaterial {
colour: Spectrum::black(),
diffuse_strength: 1.0,
}
}
}
impl Material for LambertianMaterial {
fn bsdf<'a>(&'a self) -> Box<dyn Fn(&Vec3<f64>, &Vec3<f64>, &Photon) -> Photon + 'a> {
Box::new(move |_w_o: &Vec3<f64>, _w_i: &Vec3<f64>, photon_in: &Photon| {
let mut result = self.colour.scale_photon(photon_in);
result.intensity *= self.diffuse_strength;
result
})
}
fn sample(&self, _w_i: &Vec3<f64>, _photon: &Photon) -> MaterialSampleResult {
let mut rng = thread_rng();
let mut w_o = Vec3::new(
2.0 * rng.sample::<f64, _>(Open01) - 1.0,
2.0 * rng.sample::<f64, _>(Open01) - 1.0,
0.0,
);
while w_o.norm_squared() > 1.0 {
w_o = Vec3::new(
2.0 * rng.sample::<f64, _>(Open01) - 1.0,
2.0 * rng.sample::<f64, _>(Open01) - 1.0,
0.0,
);
}
w_o.coords[2] = (1.0 - w_o.x() * w_o.x() - w_o.y() * w_o.y())
.sqrt()
.max(0.0);
let cos_theta = w_o.dot(&Vec3::unit_z());
let sin_theta = (1.0 - cos_theta * cos_theta).sqrt();
MaterialSampleResult {
direction: w_o.normalize(),
pdf: (cos_theta * sin_theta) / PI,
}
}
}

36
src/materials/mod.rs
View File

@ -1,36 +0,0 @@
use crate::math::Vec3;
use super::colour::Photon;
use super::random_distributions::{CosineWeightedHemisphere, RandomDistribution};
use std::fmt::Debug;
pub mod lambertian_material;
pub use lambertian_material::LambertianMaterial;
pub mod phong_material;
pub use phong_material::PhongMaterial;
pub mod reflective_material;
pub use reflective_material::ReflectiveMaterial;
pub mod smooth_transparent_dialectric;
pub use smooth_transparent_dialectric::SmoothTransparentDialectric;
pub struct MaterialSampleResult {
pub direction: Vec3<f64>,
pub pdf: f64,
}
type BsdfFunc<'a> = Box<dyn Fn(&Vec3<f64>, &Vec3<f64>, &Photon) -> Photon + 'a>;
pub trait Material: Debug + Sync + Send {
fn bsdf<'a>(&'a self) -> BsdfFunc<'a>;
fn sample(&self, _w_i: &Vec3<f64>, _photon: &Photon) -> MaterialSampleResult {
let distribution = CosineWeightedHemisphere::new();
let direction = distribution.value();
let pdf = distribution.pdf(direction);
MaterialSampleResult { direction, pdf }
}
}

37
src/materials/phong_material.rs
View File

@ -1,37 +0,0 @@
use crate::colour::{Photon, Spectrum};
use crate::math::Vec3;
use std::fmt::Debug;
use super::Material;
#[derive(Debug)]
pub struct PhongMaterial {
pub colour: Spectrum,
pub diffuse_strength: f64,
pub specular_strength: f64,
pub smoothness: f64,
}
impl Material for PhongMaterial {
fn bsdf<'a>(&'a self) -> Box<dyn Fn(&Vec3<f64>, &Vec3<f64>, &Photon) -> Photon + 'a> {
Box::new(move |w_o: &Vec3<f64>, w_i: &Vec3<f64>, photon_in: &Photon| {
if w_i.z() < 0.0 || w_o.z() < 0.0 {
Photon {
wavelength: photon_in.wavelength,
intensity: 0.0,
}
} else {
let reflection_vector = Vec3::new(-w_i.x(), -w_i.y(), w_i.z());
let intensity = self.colour.scale_photon(photon_in).intensity
* self.diffuse_strength
+ w_o.dot(&reflection_vector).abs().powf(self.smoothness)
* (self.specular_strength / w_i.dot(&Vec3::unit_z()));
Photon {
wavelength: photon_in.wavelength,
intensity,
}
}
})
}
}

48
src/materials/reflective_material.rs
View File

@ -1,48 +0,0 @@
use crate::colour::{Photon, Spectrum};
use crate::math::Vec3;
use std::fmt::Debug;
use super::{Material, MaterialSampleResult};
#[derive(Debug)]
pub struct ReflectiveMaterial {
pub colour: Spectrum,
pub diffuse_strength: f64,
pub reflection_strength: f64,
}
impl Material for ReflectiveMaterial {
fn bsdf<'a>(&'a self) -> Box<dyn Fn(&Vec3<f64>, &Vec3<f64>, &Photon) -> Photon + 'a> {
Box::new(move |w_o: &Vec3<f64>, w_i: &Vec3<f64>, photon_in: &Photon| {
if w_i.z() <= 0.0 || w_o.z() <= 0.0 {
Photon {
wavelength: photon_in.wavelength,
intensity: 0.0,
}
} else {
let reflection_vector = Vec3::new(-w_o.x(), -w_o.y(), w_o.z());
let mut photon_out = self.colour.scale_photon(photon_in);
photon_out.intensity *= self.diffuse_strength;
let sigma = 0.05;
let two = 2.0;
// These are normalized vectors, but sometimes rounding errors cause the
// dot product to be slightly above 1 or below 0. The call to clamp
// ensures the values stay within the domain of acos,
let theta = w_i.dot(&reflection_vector).clamp(0.0, 1.0).abs().acos();
let reflection_factor =
self.reflection_strength * (-(theta.powf(two)) / (two * sigma * sigma)).exp();
photon_out.intensity =
photon_out.intensity * (1.0 - reflection_factor) + reflection_factor;
photon_out
}
})
}
fn sample(&self, w_o: &Vec3<f64>, _photon: &Photon) -> MaterialSampleResult {
MaterialSampleResult {
direction: Vec3::new(-w_o.x(), -w_o.y(), w_o.z()),
pdf: 1.0,
}
}
}

84
src/materials/rgb_sampled_bsdf_material.rs
View File

@ -1,84 +0,0 @@
use super::{Bsdf, Material};
use crate::colour::ColourRgbF;
use crate::math::Vec3;
use crate::realtype::NormalizedToU32;
use std::error::Error;
use std::fs::File;
use std::io::BufReader;
use std::sync::Arc;
use std::f64::consts::{FRAC_PI_2, PI};
#[derive(Debug)]
pub struct RgbSampledBsdfMaterial {
lut: Arc<Vec<Vec<Vec<Vec<Vec3>>>>>,
}
fn expand_and_index<T: Clone>(v: &mut Vec<T>, i: usize, default: T) -> &mut T {
if v.len() < i + 1 {
v.resize(i + 1, default);
}
&mut v[i]
}
impl RgbSampledBsdfMaterial {
pub fn from_csv_file(filename: &str) -> Result<RgbSampledBsdfMaterial, Box<dyn Error>> {
let csv_file = File::open(filename)?;
let mut reader = csv::Reader::from_reader(BufReader::new(&csv_file));
let mut lut = Vec::new();
for row_result in reader.records() {
let row = row_result?;
let theta_in_index = row[0].trim().parse::<usize>()?;
let phi_in_index = row[1].trim().parse::<usize>()?;
let theta_out_index = row[2].trim().parse::<usize>()?;
let phi_out_index = row[3].trim().parse::<usize>()?;
let red = row[4].trim().parse::<f64>()?;
let green = row[5].trim().parse::<f64>()?;
let blue = row[6].trim().parse::<f64>()?;
*expand_and_index(
expand_and_index(
expand_and_index(
expand_and_index(&mut lut, theta_in_index, Vec::new()),
phi_in_index,
Vec::new(),
),
theta_out_index,
Vec::new(),
),
phi_out_index,
Vec3::zeros(),
) = Vec3::new(red, green, blue);
}
let lut = Arc::new(lut);
Ok(RgbSampledBsdfMaterial { lut })
}
}
impl<'a> Material for RgbSampledBsdfMaterial {
fn bsdf(&self) -> Bsdf {
let lut = Arc::clone(&self.lut);
Box::new(move |w_in, w_out, colour_in| {
if w_in.z() < 0.0 || w_out.z() < 0.0 {
return ColourRgbF::new(0.0, 0.0, 0.0);
}
let theta_in = w_in.z().acos();
let theta_in_index = (theta_in / FRAC_PI_2).normalized_to_u32(4) as usize;
let phi_in = w_in.y().atan2(w_in.x()) + PI;
let phi_in_index = (phi_in / (2.0 * PI)).normalized_to_u32(6) as usize;
let theta_out = w_out.z().acos();
let theta_out_index = (theta_out / FRAC_PI_2).normalized_to_u32(4) as usize;
let phi_out = w_out.y().atan2(w_out.x()) + PI;
let phi_out_index = (phi_out / (2.0 * PI)).normalized_to_u32(6) as usize;
ColourRgbF::from_vec3(
&colour_in.as_vec3().component_mul(
&lut[theta_in_index][phi_in_index][theta_out_index][phi_out_index],
),
)
})
}
fn sample(&self, w_o: &Vec3) -> Vec<Vec3> {
vec![Vec3::new(-w_o.x(), -w_o.y(), w_o.z())]
}
}

115
src/materials/smooth_transparent_dialectric.rs
View File

@ -1,115 +0,0 @@
use crate::colour::{Photon, Spectrum};
use crate::materials::{Material, MaterialSampleResult};
use crate::math::Vec3;
use rand::random;
#[derive(Debug)]
struct FresnelResult {
reflection_direction: Vec3<f64>,
reflection_strength: f64,
transmission_direction: Vec3<f64>,
transmission_strength: f64,
}
fn fresnel(w_i: &Vec3<f64>, eta1: f64, eta2: f64) -> FresnelResult {
let normal = if w_i.z() > 0.0 {
Vec3::unit_z()
} else {
-Vec3::unit_z()
};
let reflection_direction = Vec3::new(-w_i.x(), -w_i.y(), w_i.z());
let r = eta1 / eta2;
let cos_theta1 = normal.dot(w_i);
let cos_theta2_squared = 1.0 - r * r * (1.0 - cos_theta1 * cos_theta1);
let mut result = if cos_theta2_squared >= 0.0 {
let cos_theta2 = cos_theta2_squared.sqrt();
let reflection_strength_parallel_sqrt =
(eta1 * cos_theta2 - eta2 * cos_theta1) / (eta1 * cos_theta2 + eta2 * cos_theta1);
let reflection_strength_perpendicular_sqrt =
(eta1 * cos_theta1 - eta2 * cos_theta2) / (eta1 * cos_theta1 + eta2 * cos_theta2);
let reflection_strength = 0.5
* (reflection_strength_parallel_sqrt * reflection_strength_parallel_sqrt
+ reflection_strength_perpendicular_sqrt * reflection_strength_perpendicular_sqrt);
let transmission_direction =
(-r * w_i + (r * cos_theta1 - cos_theta2) * normal).normalize();
let transmission_strength = 1.0 - reflection_strength;
FresnelResult {
reflection_direction,
reflection_strength,
transmission_direction,
transmission_strength,
}
} else {
let reflection_strength = 1.0;
let transmission_strength = 0.0;
let transmission_direction = Default::default();
FresnelResult {
reflection_direction,
reflection_strength,
transmission_direction,
transmission_strength,
}
};
if w_i.z() < 0.0 {
result.reflection_direction.coords[2] *= -1.0;
result.transmission_direction.coords[2] *= -1.0;
}
result
}
#[derive(Debug)]
pub struct SmoothTransparentDialectric {
eta: Spectrum,
}
impl SmoothTransparentDialectric {
pub fn new(eta: Spectrum) -> SmoothTransparentDialectric {
SmoothTransparentDialectric { eta }
}
}
impl Material for SmoothTransparentDialectric {
fn bsdf<'a>(&'a self) -> Box<dyn Fn(&Vec3<f64>, &Vec3<f64>, &Photon) -> Photon + 'a> {
Box::new(move |w_o: &Vec3<f64>, w_i: &Vec3<f64>, photon_in: &Photon| {
let (eta1, eta2) = if w_i.z() >= 0.0 {
(1.0, self.eta.intensity_at_wavelength(photon_in.wavelength))
} else {
(self.eta.intensity_at_wavelength(photon_in.wavelength), 1.0)
};
let fresnel = fresnel(w_i, eta1, eta2);
if (*w_o - fresnel.reflection_direction).norm_squared() < 0.0000000001 {
photon_in.scale_intensity(fresnel.reflection_strength)
} else if (*w_o - fresnel.transmission_direction).norm_squared() < 0.0000000001 {
photon_in.scale_intensity(fresnel.transmission_strength)
} else {
photon_in.set_intensity(0.0)
}
})
}
fn sample(&self, w_i: &Vec3<f64>, photon: &Photon) -> MaterialSampleResult {
let (eta1, eta2) = if w_i.z() >= 0.0 {
(1.0, self.eta.intensity_at_wavelength(photon.wavelength))
} else {
(self.eta.intensity_at_wavelength(photon.wavelength), 1.0)
};
let fresnel = fresnel(w_i, eta1, eta2);
if fresnel.transmission_strength <= 0.0000000001 {
MaterialSampleResult {
direction: fresnel.reflection_direction,
pdf: 0.5,
}
} else if fresnel.reflection_strength <= 0.0000000001 || random() {
MaterialSampleResult {
direction: fresnel.transmission_direction,
pdf: 0.5,
}
} else {
MaterialSampleResult {
direction: fresnel.reflection_direction,
pdf: 0.5,
}
}
}
}

160
src/math/affine3.rs
View File

@ -1,160 +0,0 @@
use super::{Float, Mat3, Mat4, Vec3, Vec4};
use std::ops::{Mul, MulAssign};
#[derive(PartialEq, Debug)]
pub struct Affine3<T: Float> {
matrix: Mat4<T>,
}
impl<T: Float> Affine3<T> {
pub fn translation(delta: Vec3<T>) -> Self {
#[rustfmt::skip]
let matrix = Mat4::new(T::one(), T::zero(), T::zero(), delta.x(),
T::zero(), T::one() , T::zero(), delta.y(),
T::zero(), T::zero(), T::one(), delta.z(),
T::zero(), T::zero(), T::zero(), T::one());
Self { matrix }
}
pub fn rotation(axis: Vec3<T>, angle: T) -> Self {
let x = axis.x();
let y = axis.y();
let z = axis.z();
let cos = angle.cos();
let ncos = T::one() - cos;
let sin = angle.sin();
#[rustfmt::skip]
let matrix = Mat4::new(x*x*ncos+cos, y*x*ncos-z*sin, z*x*ncos+y*sin, T::zero(),
x*y*ncos+z*sin, y*y*ncos+cos, z*y*ncos-x*sin, T::zero(),
x*z*ncos-y*sin, y*z*ncos+x*sin, z*z*ncos+cos, T::zero(),
T::zero(), T::zero(), T::zero(), T::one());
Self { matrix }
}
pub fn scale(s: T) -> Self {
#[rustfmt::skip]
let matrix = Mat4::new(s, T::zero(), T::zero(), T::zero(),
T::zero(), s , T::zero(), T::zero(),
T::zero(), T::zero(), s, T::zero(),
T::zero(), T::zero(), T::zero(), T::one());
Self { matrix }
}
pub fn get_element(&self, row: usize, column: usize) -> T {
self.matrix.get_element(row, column)
}
pub fn get_row(&self, row: usize) -> Vec4<T> {
self.matrix.get_row(row)
}
pub fn get_column(&self, column: usize) -> Vec4<T> {
self.matrix.get_column(column)
}
pub fn linear_map(&self) -> Mat3<T> {
Mat3::new(
self.matrix.get_element(0, 0),
self.matrix.get_element(0, 1),
self.matrix.get_element(0, 2),
self.matrix.get_element(1, 0),
self.matrix.get_element(1, 1),
self.matrix.get_element(1, 2),
self.matrix.get_element(2, 0),
self.matrix.get_element(2, 1),
self.matrix.get_element(2, 2),
)
}
pub fn inverse(&self) -> Affine3<T> {
// linear map should always be invertable.
let inner = self.linear_map().try_inverse().unwrap();
let translation = inner * self.matrix.get_column(3).xyz();
#[rustfmt::skip]
let matrix = Mat4::new(
inner.get_element(0,0), inner.get_element(0,1), inner.get_element(0,2), translation.x(),
inner.get_element(1,0), inner.get_element(1,1), inner.get_element(1,2), translation.y(),
inner.get_element(2,0), inner.get_element(2,1), inner.get_element(2,2), translation.z(),
T::zero(), T::zero(), T::zero(), T::one());
Self { matrix }
}
}
impl<T: Float> Mul<Affine3<T>> for Affine3<T> {
type Output = Self;
fn mul(self, rhs: Affine3<T>) -> Affine3<T> {
let matrix = self.matrix * rhs.matrix;
Affine3 { matrix }
}
}
impl<T: Float> Mul<Mat4<T>> for Affine3<T> {
type Output = Mat4<T>;
fn mul(self, rhs: Mat4<T>) -> Mat4<T> {
self.matrix * rhs
}
}
impl<T: Float> Mul<Affine3<T>> for Mat4<T> {
type Output = Mat4<T>;
fn mul(self, rhs: Affine3<T>) -> Mat4<T> {
self * rhs.matrix
}
}
impl<T: Float> MulAssign<Affine3<T>> for Affine3<T> {
fn mul_assign(&mut self, rhs: Affine3<T>) {
self.matrix *= rhs.matrix
}
}
impl<T: Float> Mul<Vec4<T>> for Affine3<T> {
type Output = Vec4<T>;
fn mul(self, rhs: Vec4<T>) -> Vec4<T> {
self.matrix * rhs
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn translate_translates_vector() {
let p = Vec4::new(1.0, 2.0, 3.0, 1.0);
let v = Vec3::new(4.0, 5.0, 6.0);
let target = Affine3::translation(v);
let diff = (target * p).xyz() - (p.xyz() + v);
assert!(diff.norm() < 0.0000000001);
}
#[test]
fn rotate_rotates_vector() {
let x = Vec4::new(1.0, 0.0, 0.0, 1.0);
let y = Vec4::new(0.0, 1.0, 0.0, 1.0);
let z = Vec4::new(0.0, 0.0, 1.0, 1.0);
let target = Affine3::rotation(z.xyz(), std::f64::consts::PI/2.0) * y;
let diff = -x.xyz() - target.xyz();
assert!(diff.norm() < 0.0000000001);
}
#[test]
fn linear_map_is_inner_matrix() {
#[rustfmt::skip]
let target = Affine3{
matrix: Mat4::new(1.0, 2.0, 3.0, 4.0,
5.0, 6.0, 7.0, 8.0,
9.0, 10.0, 11.0, 12.0,
0.0, 0.0, 0.0, 1.0)};
let linear_map = target.linear_map();
for i in 0..2 {
for j in 0..2 {
assert!(linear_map.get_element(i, j) == target.get_element(i, j));
}
}
}
}

35
src/math/mat2.rs
View File

@ -1,35 +0,0 @@
use super::Float;
#[derive(PartialEq, Debug, Copy, Clone)]
pub struct Mat2<T: Float> {
pub elements: [[T; 2]; 2],
}
impl<T: Float> Mat2<T> {
pub fn new(m00: T, m01: T, m10: T, m11: T) -> Mat2<T> {
Mat2 {
elements: [[m00, m01], [m10, m11]],
}
}
pub fn determinant(&self) -> T {
self.elements[0][0] * self.elements[1][1] - self.elements[0][1] * self.elements[1][0]
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn determinant_returns_expected_value() {
let target1 = Mat2::new(1.0, 2.0, 3.0, 4.0);
let target2 = Mat2::new(1.0, -2.0, 3.0, 4.0);
let target3 = Mat2::new(1.0, 1.0, 1.0, 1.0);
let target4 = Mat2::new(21.0, 45.0, -16.0, 0.0);
assert!(target1.determinant() == -2.0);
assert!(target2.determinant() == 10.0);
assert!(target3.determinant() == 0.0);
assert!(target4.determinant() == 720.0);
}
}

351
src/math/mat3.rs
View File

@ -1,351 +0,0 @@
use super::{Float, Mat2, Vec3};
use std::ops::{Mul, MulAssign};
#[derive(PartialEq, Debug, Copy, Clone)]
pub struct Mat3<T: Float> {
elements: [[T; 3]; 3],
}
impl<T: Float> Mat3<T> {
#[allow(clippy::too_many_arguments)]
pub fn new(m00: T, m01: T, m02: T, m10: T, m11: T, m12: T, m20: T, m21: T, m22: T) -> Mat3<T> {
Mat3 {
elements: [[m00, m01, m02], [m10, m11, m12], [m20, m21, m22]],
}
}
pub fn identity() -> Mat3<T> {
Mat3 {
elements: [
[T::one(), T::zero(), T::zero()],
[T::zero(), T::one(), T::zero()],
[T::zero(), T::zero(), T::one()],
],
}
}
pub fn from_rows(r0: &Vec3<T>, r1: &Vec3<T>, r2: &Vec3<T>) -> Mat3<T> {
let mut elements = [[T::zero(); 3]; 3];
for (row, v) in elements.iter_mut().zip([r0, r1, r2].iter()) {
for (it, val) in row.iter_mut().zip(v.coords.iter()) {
*it = *val;
}
}
Mat3 { elements }
}
pub fn get_element(&self, row: usize, column: usize) -> T {
self.elements[row][column]
}
pub fn get_row(&self, row: usize) -> Vec3<T> {
Vec3 {
coords: self.elements[row],
}
}
pub fn get_column(&self, column: usize) -> Vec3<T> {
let mut coords = [T::zero(); 3];
for (coord, row) in coords.iter_mut().zip(self.elements.iter()) {
*coord = row[column];
}
Vec3 { coords }
}
pub fn transpose(&self) -> Mat3<T> {
let mut elements = [[T::zero(); 3]; 3];
for i in 0..3 {
for j in 0..3 {
elements[i][j] = self.elements[j][i];
}
}
Mat3 { elements }
}
pub fn first_minor(&self, row: usize, column: usize) -> T {
let mut elements = [[T::zero(); 2]; 2];
let mut i_dst = 0;
let mut j_dst = 0;
for i_src in 0..3 {
if i_src != row {
for j_src in 0..3 {
if j_src != column {
elements[i_dst][j_dst] = self.get_element(i_src, j_src);
j_dst += 1;
}
}
i_dst += 1;
j_dst = 0;
}
}
let minor_matrix = Mat2 { elements };
minor_matrix.determinant()
}
pub fn cofactor(&self, row: usize, column: usize) -> T {
T::from((-1i32).pow((row + column) as u32)) * self.first_minor(row, column)
}
pub fn cofactor_matrix(&self) -> Mat3<T> {
let mut elements = [[T::zero(); 3]; 3];
for i in 0..3 {
for j in 0..3 {
elements[i][j] = self.cofactor(i, j);
}
}
Mat3 { elements }
}
pub fn determinant(&self) -> T {
self.elements[0][0] * self.first_minor(0, 0) - self.elements[0][1] * self.first_minor(0, 1)
+ self.elements[0][2] * self.first_minor(0, 2)
}
pub fn try_inverse(&self) -> Option<Mat3<T>> {
let determinant = self.determinant();
if determinant == T::zero() {
None
} else {
Some(self.cofactor_matrix().transpose() * determinant)
}
}
}
impl<T: Float> Mul<Mat3<T>> for Mat3<T> {
type Output = Self;
fn mul(self, rhs: Self) -> Self {
let mut elements = [[T::zero(); 3]; 3];
for row in 0..3 {
for column in 0..3 {
elements[row][column] = self.get_row(row).dot(&rhs.get_column(column));
}
}
Mat3 { elements }
}
}
impl<T: Float> MulAssign<Mat3<T>> for Mat3<T> {
fn mul_assign(&mut self, rhs: Self) {
for row in 0..3 {
let mut new_row = [T::zero(); 3];
for column in 0..3 {
new_row[column] = self.get_row(row).dot(&rhs.get_column(column));
}
self.elements[row] = new_row;
}
}
}
impl<T: Float> Mul<Vec3<T>> for Mat3<T> {
type Output = Vec3<T>;
fn mul(self, rhs: Vec3<T>) -> Vec3<T> {
let mut coords = [T::zero(); 3];
for (coord, row) in coords.iter_mut().zip(self.elements.iter()) {
*coord = Vec3 { coords: *row }.dot(&rhs);
}
Vec3 { coords }
}
}
impl<T: Float> Mul<&Vec3<T>> for Mat3<T> {
type Output = Vec3<T>;
fn mul(self, rhs: &Vec3<T>) -> Vec3<T> {
let mut coords = [T::zero(); 3];
for (coord, row) in coords.iter_mut().zip(self.elements.iter()) {
*coord = Vec3 { coords: *row }.dot(rhs);
}
Vec3 { coords }
}
}
impl<T: Float> Mul<T> for Mat3<T> {
type Output = Mat3<T>;
fn mul(self, rhs: T) -> Mat3<T> {
let mut elements = [[T::zero(); 3]; 3];
for i in 0..3 {
for j in 0..3 {
elements[i][j] = self.elements[i][j] * rhs;
}
}
Mat3 { elements }
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn elements_are_in_expected_locations() {
let target = Mat3::new(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0);
assert!(target.get_element(0, 0) == 1.0);
assert!(target.get_element(0, 1) == 2.0);
assert!(target.get_element(0, 2) == 3.0);
assert!(target.get_element(1, 0) == 4.0);
assert!(target.get_element(1, 1) == 5.0);
assert!(target.get_element(1, 2) == 6.0);
assert!(target.get_element(2, 0) == 7.0);
assert!(target.get_element(2, 1) == 8.0);
assert!(target.get_element(2, 2) == 9.0);
}
#[test]
fn from_rows_places_values_in_rows() {
let target = Mat3::from_rows(
&Vec3::new(1.0, 2.0, 3.0),
&Vec3::new(4.0, 5.0, 6.0),
&Vec3::new(7.0, 8.0, 9.0),
);
assert!(target.get_element(0, 0) == 1.0);
assert!(target.get_element(0, 1) == 2.0);
assert!(target.get_element(0, 2) == 3.0);
assert!(target.get_element(1, 0) == 4.0);
assert!(target.get_element(1, 1) == 5.0);
assert!(target.get_element(1, 2) == 6.0);
assert!(target.get_element(2, 0) == 7.0);
assert!(target.get_element(2, 1) == 8.0);
assert!(target.get_element(2, 2) == 9.0);
}
#[test]
fn get_column_returns_expected_value() {
let target = Mat3::from_rows(
&Vec3::new(1.0, 2.0, 3.0),
&Vec3::new(4.0, 5.0, 6.0),
&Vec3::new(7.0, 8.0, 9.0),
);
assert!(target.get_column(0) == Vec3::new(1.0, 4.0, 7.0));
assert!(target.get_column(1) == Vec3::new(2.0, 5.0, 8.0));
assert!(target.get_column(2) == Vec3::new(3.0, 6.0, 9.0));
}
#[test]
fn transpose_returns_expected_result() {
let target = Mat3::from_rows(
&Vec3::new(1.0, 2.0, 3.0),
&Vec3::new(4.0, 5.0, 6.0),
&Vec3::new(7.0, 8.0, 9.0),
);
let expected = Mat3::from_rows(
&Vec3::new(1.0, 4.0, 7.0),
&Vec3::new(2.0, 5.0, 8.0),
&Vec3::new(3.0, 6.0, 9.0),
);
assert!(target.transpose() == expected);
}
#[test]
fn cofactor_matrix_returns_expected_result() {
let target = Mat3::from_rows(
&Vec3::new(1.0, 2.0, 3.0),
&Vec3::new(4.0, 5.0, 6.0),
&Vec3::new(7.0, 8.0, 9.0),
);
let expected = Mat3::from_rows(
&Vec3::new(-3.0, 6.0, -3.0),
&Vec3::new(6.0, -12.0, 6.0),
&Vec3::new(-3.0, 6.0, -3.0),
);
assert!(target.cofactor_matrix() == expected);
}
#[test]
fn determinant_returns_expected_result() {
let target = Mat3::from_rows(
&Vec3::new(1.0, 3.0, 2.0),
&Vec3::new(4.0, 5.0, 6.0),
&Vec3::new(7.0, 8.0, 9.0),
);
assert!(target.determinant() == 9.0);
}
#[test]
fn inverse_of_singular_matrix_is_none_result() {
let target = Mat3::from_rows(
&Vec3::new(1.0, 2.0, 3.0),
&Vec3::new(4.0, 5.0, 6.0),
&Vec3::new(7.0, 8.0, 9.0),
);
let expected = None;
assert!(target.try_inverse() == expected);
}
#[test]
fn inverse_of_identity_is_identity() {
assert!(Mat3::<f64>::identity().try_inverse() == Some(Mat3::identity()));
}
#[test]
fn inverse_returns_expected_result() {
let target = Mat3::from_rows(
&Vec3::new(4.0, -5.0, -2.0),
&Vec3::new(5.0, -6.0, -2.0),
&Vec3::new(-8.0, 9.0, 3.0),
);
let expected = Some(Mat3::from_rows(
&Vec3::new(0.0, -3.0, -2.0),
&Vec3::new(1.0, -4.0, -2.0),
&Vec3::new(-3.0, 4.0, 1.0),
));
assert!(target.try_inverse() == expected);
}
#[test]
fn mul_with_mat3_returns_expected_result() {
let a = Mat3::from_rows(
&Vec3::new(1.0, 2.0, 3.0),
&Vec3::new(4.0, 5.0, 6.0),
&Vec3::new(7.0, 8.0, 9.0),
);
let b = Mat3::from_rows(
&Vec3::new(10.0, 11.0, 12.0),
&Vec3::new(13.0, 14.0, 15.0),
&Vec3::new(16.0, 17.0, 18.0),
);
let c = Mat3::from_rows(
&Vec3::new(84.0, 90.0, 96.0),
&Vec3::new(201.0, 216.0, 231.0),
&Vec3::new(318.0, 342.0, 366.0),
);
assert!(a * b == c);
}
#[test]
fn mul_assign_returns_expected_result() {
let mut a = Mat3::from_rows(
&Vec3::new(1.0, 2.0, 3.0),
&Vec3::new(4.0, 5.0, 6.0),
&Vec3::new(7.0, 8.0, 9.0),
);
let b = Mat3::from_rows(
&Vec3::new(10.0, 11.0, 12.0),
&Vec3::new(13.0, 14.0, 15.0),
&Vec3::new(16.0, 17.0, 18.0),
);
let c = Mat3::from_rows(
&Vec3::new(84.0, 90.0, 96.0),
&Vec3::new(201.0, 216.0, 231.0),
&Vec3::new(318.0, 342.0, 366.0),
);
a *= b;
assert!(a == c);
}
#[test]
fn mul_with_vec3_returns_expected_result() {
let a = Mat3::from_rows(
&Vec3::new(1.0, 2.0, 3.0),
&Vec3::new(4.0, 5.0, 6.0),
&Vec3::new(7.0, 8.0, 9.0),
);
let b = Vec3::new(10.0, 11.0, 12.0);
let c = Vec3::new(68.0, 167.0, 266.0);
assert!(a * b == c);
}
}

234
src/math/mat4.rs
View File

@ -1,234 +0,0 @@
use super::{Float,Vec4};
use std::ops::{Mul, MulAssign};
#[derive(PartialEq, Debug)]
pub struct Mat4<T:Float> {
elements: [[T; 4]; 4],
}
impl<T:Float> Mat4<T> {
#[allow(clippy::too_many_arguments)]
pub fn new(
m00: T,
m01: T,
m02: T,
m03: T,
m10: T,
m11: T,
m12: T,
m13: T,
m20: T,
m21: T,
m22: T,
m23: T,
m30: T,
m31: T,
m32: T,
m33: T,
) -> Mat4<T> {
Mat4 {
elements: [
[m00, m01, m02, m03],
[m10, m11, m12, m13],
[m20, m21, m22, m23],
[m30, m31, m32, m33],
],
}
}
pub fn from_rows(r0: &Vec4<T>, r1: &Vec4<T>, r2: &Vec4<T>, r3: &Vec4<T>) -> Mat4<T> {
let mut elements = [[T::zero(); 4]; 4];
for (row, v) in elements.iter_mut().zip([r0, r1, r2, r3].iter()) {
for (it, val) in row.iter_mut().zip(v.coords.iter()) {
*it = *val;
}
}
Mat4 { elements }
}
pub fn get_element(&self, row: usize, column: usize) -> T {
self.elements[row][column]
}
pub fn get_row(&self, row: usize) -> Vec4<T> {
Vec4 {
coords: self.elements[row],
}
}
pub fn get_column(&self, column: usize) -> Vec4<T> {
let mut coords = [T::zero(); 4];
for (coord, row) in coords.iter_mut().zip(self.elements.iter()) {
*coord = row[column];
}
Vec4 { coords }
}
}
impl<T:Float> Mul<Mat4<T>> for Mat4<T> {
type Output = Self;
fn mul(self, rhs: Self) -> Self {
let mut elements = [[T::zero(); 4]; 4];
for row in 0..4 {
for column in 0..4 {
elements[row][column] = self.get_row(row).dot(&rhs.get_column(column));
}
}
Mat4 { elements }
}
}
impl<T:Float> MulAssign<Mat4<T>> for Mat4<T> {
fn mul_assign(&mut self, rhs: Self) {
for row in 0..4 {
let mut new_row = [T::zero(); 4];
for column in 0..4 {
new_row[column] = self.get_row(row).dot(&rhs.get_column(column));
}
self.elements[row] = new_row;
}
}
}
impl<T:Float> Mul<Vec4<T>> for Mat4<T> {
type Output = Vec4<T>;
fn mul(self, rhs: Vec4<T>) -> Vec4<T> {
let mut coords = [T::zero(); 4];
for (coord, row) in coords.iter_mut().zip(self.elements.iter()) {
*coord = Vec4 { coords: *row }.dot(&rhs);
}
Vec4 { coords }
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn elements_are_in_expected_locations() {
let target = Mat4::new(
1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0,
);
assert!(target.get_element(0, 0) == 1.0);
assert!(target.get_element(0, 1) == 2.0);
assert!(target.get_element(0, 2) == 3.0);
assert!(target.get_element(0, 3) == 4.0);
assert!(target.get_element(1, 0) == 5.0);
assert!(target.get_element(1, 1) == 6.0);
assert!(target.get_element(1, 2) == 7.0);
assert!(target.get_element(1, 3) == 8.0);
assert!(target.get_element(2, 0) == 9.0);
assert!(target.get_element(2, 1) == 10.0);
assert!(target.get_element(2, 2) == 11.0);
assert!(target.get_element(2, 3) == 12.0);
assert!(target.get_element(3, 0) == 13.0);
assert!(target.get_element(3, 1) == 14.0);
assert!(target.get_element(3, 2) == 15.0);
assert!(target.get_element(3, 3) == 16.0);
}
#[test]
fn from_rows_places_values_in_rows() {
let target = Mat4::from_rows(
&Vec4::new(1.0, 2.0, 3.0, 4.0),
&Vec4::new(5.0, 6.0, 7.0, 8.0),
&Vec4::new(9.0, 10.0, 11.0, 12.0),
&Vec4::new(13.0, 14.0, 15.0, 16.0),
);
assert!(target.get_element(0, 0) == 1.0);
assert!(target.get_element(0, 1) == 2.0);
assert!(target.get_element(0, 2) == 3.0);
assert!(target.get_element(0, 3) == 4.0);
assert!(target.get_element(1, 0) == 5.0);
assert!(target.get_element(1, 1) == 6.0);
assert!(target.get_element(1, 2) == 7.0);
assert!(target.get_element(1, 3) == 8.0);
assert!(target.get_element(2, 0) == 9.0);
assert!(target.get_element(2, 1) == 10.0);
assert!(target.get_element(2, 2) == 11.0);
assert!(target.get_element(2, 3) == 12.0);
assert!(target.get_element(3, 0) == 13.0);
assert!(target.get_element(3, 1) == 14.0);
assert!(target.get_element(3, 2) == 15.0);
assert!(target.get_element(3, 3) == 16.0);
}
#[test]
fn get_column_returns_expected_value() {
let target = Mat4::from_rows(
&Vec4::new(1.0, 2.0, 3.0, 4.0),
&Vec4::new(5.0, 6.0, 7.0, 8.0),
&Vec4::new(9.0, 10.0, 11.0, 12.0),
&Vec4::new(13.0, 14.0, 15.0, 16.0),
);
assert!(target.get_column(0) == Vec4::new(1.0, 5.0, 9.0, 13.0));
assert!(target.get_column(1) == Vec4::new(2.0, 6.0, 10.0, 14.0));
assert!(target.get_column(2) == Vec4::new(3.0, 7.0, 11.0, 15.0));
assert!(target.get_column(3) == Vec4::new(4.0, 8.0, 12.0, 16.0));
}
#[test]
fn mul_with_mat4_returns_expected_result() {
let a = Mat4::from_rows(
&Vec4::new(1.0, 2.0, 3.0, 4.0),
&Vec4::new(5.0, 6.0, 7.0, 8.0),
&Vec4::new(9.0, 10.0, 11.0, 12.0),
&Vec4::new(13.0, 14.0, 15.0, 16.0),
);
let b = Mat4::from_rows(
&Vec4::new(17.0, 18.0, 19.0, 20.0),
&Vec4::new(21.0, 22.0, 23.0, 24.0),
&Vec4::new(25.0, 26.0, 27.0, 28.0),
&Vec4::new(29.0, 30.0, 31.0, 32.0),
);
let c = Mat4::from_rows(
&Vec4::new(250.0, 260.0, 270.0, 280.0),
&Vec4::new(618.0, 644.0, 670.0, 696.0),
&Vec4::new(986.0, 1028.0, 1070.0, 1112.0),
&Vec4::new(1354.0, 1412.0, 1470.0, 1528.0),
);
assert!(a * b == c);
}
#[test]
fn mul_assign_returns_expected_result() {
let mut a = Mat4::from_rows(
&Vec4::new(1.0, 2.0, 3.0, 4.0),
&Vec4::new(5.0, 6.0, 7.0, 8.0),
&Vec4::new(9.0, 10.0, 11.0, 12.0),
&Vec4::new(13.0, 14.0, 15.0, 16.0),
);
let b = Mat4::from_rows(
&Vec4::new(17.0, 18.0, 19.0, 20.0),
&Vec4::new(21.0, 22.0, 23.0, 24.0),
&Vec4::new(25.0, 26.0, 27.0, 28.0),
&Vec4::new(29.0, 30.0, 31.0, 32.0),
);
let c = Mat4::from_rows(
&Vec4::new(250.0, 260.0, 270.0, 280.0),
&Vec4::new(618.0, 644.0, 670.0, 696.0),
&Vec4::new(986.0, 1028.0, 1070.0, 1112.0),
&Vec4::new(1354.0, 1412.0, 1470.0, 1528.0),
);
a *= b;
assert!(a == c);
}
#[test]
fn mul_with_vec4_returns_expected_result() {
let a = Mat4::from_rows(
&Vec4::new(1.0, 2.0, 3.0, 4.0),
&Vec4::new(5.0, 6.0, 7.0, 8.0),
&Vec4::new(9.0, 10.0, 11.0, 12.0),
&Vec4::new(13.0, 14.0, 15.0, 16.0),
);
let b = Vec4::new(17.0, 18.0, 19.0, 20.0);
let c = Vec4::new(190.0, 486.0, 782.0, 1078.0);
assert!(a * b == c);
}
}

23
src/math/mod.rs
View File

@ -1,23 +0,0 @@
mod number;
pub use number::*;
mod vec2;
pub use vec2::*;
mod vec3;
pub use vec3::*;
mod vec4;
pub use vec4::*;
mod mat2;
pub use mat2::*;
mod mat3;
pub use mat3::*;
mod mat4;
pub use mat4::*;
mod affine3;
pub use affine3::*;

72
src/math/number.rs
View File

@ -1,72 +0,0 @@
use std::{
cmp::PartialOrd,
iter::Sum,
ops::{Add, AddAssign, Div, Mul, MulAssign, Neg, Sub, SubAssign},
};
pub trait HasZero {
fn zero() -> Self;
}
pub trait HasOne {
fn one() -> Self;
}
pub trait Float:
Copy
+ HasZero
+ HasOne
+ Add<Output = Self>
+ AddAssign
+ Div<Output = Self>
+ Mul<Output = Self>
+ MulAssign
+ Sub<Output = Self>
+ SubAssign
+ Neg<Output = Self>
+ Sum
+ PartialOrd
+ From<f32>
+ From<f64>
+ From<u8>
+ From<i8>
+ From<u16>
+ From<i16>
+ From<u32>
+ From<i32>
{
fn abs(self) -> Self;
fn sqrt(self) -> Self;
fn sin(self) -> Self;
fn cos(self) -> Self;
}
impl HasZero for f64 {
fn zero() -> Self {
0.0
}
}
impl HasOne for f64 {
fn one() -> Self {
1.0
}
}
impl Float for f64 {
fn abs(self) -> Self {
self.abs()
}
fn sqrt(self) -> Self {
self.sqrt()
}
fn sin(self) -> Self {
self.sin()
}
fn cos(self) -> Self {
self.cos()
}
}

187
src/math/vec2.rs
View File

@ -1,187 +0,0 @@
use super::Float;
use itertools::izip;
use std::ops::{Add, AddAssign, Mul, MulAssign, Sub, SubAssign};
#[derive(PartialEq, Debug, Copy, Clone)]
pub struct Vec2<T: Float> {
coords: [T; 2],
}
impl<T: Float> Vec2<T> {
pub fn new(x: T, y: T) -> Self {
Vec2 { coords: [x, y] }
}
pub fn x(&self) -> T {
self.coords[0]
}
pub fn y(&self) -> T {
self.coords[1]
}
pub fn dot(&self, rhs: &Vec2<T>) -> T {
self.coords
.iter()
.copied()
.zip(rhs.coords.iter().copied())
.map(|(a_elem, b_elem)| a_elem * b_elem)
.sum()
}
pub fn perp(&self, rhs: &Vec2<T>) -> T {
self.x() * rhs.y() - self.y() * rhs.x()
}
}
impl<T: Float> Add for Vec2<T> {
type Output = Self;
fn add(self, rhs: Self) -> Self {
let mut coords = [T::zero(); 2];
for (r, a, b) in izip!(
coords.iter_mut(),
self.coords.iter().copied(),
rhs.coords.iter().copied()
) {
*r = a + b;
}
Vec2 { coords }
}
}
impl<T: Float> AddAssign for Vec2<T> {
fn add_assign(&mut self, rhs: Self) {
for (a, b) in self.coords.iter_mut().zip(rhs.coords.iter().copied()) {
*a += b;
}
}
}
impl<T: Float> Sub for Vec2<T> {
type Output = Self;
fn sub(self, rhs: Self) -> Self {
let mut coords = [T::zero(); 2];
for (r, a, b) in izip!(
coords.iter_mut(),
self.coords.iter().copied(),
rhs.coords.iter().copied()
) {
*r = a - b;
}
Vec2 { coords }
}
}
impl<T: Float> SubAssign for Vec2<T> {
fn sub_assign(&mut self, rhs: Self) {
for (a, b) in self.coords.iter_mut().zip(rhs.coords.iter().copied()) {
*a -= b;
}
}
}
impl<T: Float> Mul<T> for Vec2<T> {
type Output = Self;
fn mul(self, rhs: T) -> Vec2<T> {
let mut coords = [T::zero(); 2];
for (r, a) in coords.iter_mut().zip(self.coords.iter().copied()) {
*r = a * rhs;
}
Vec2 { coords }
}
}
impl<T: Float> MulAssign<T> for Vec2<T> {
fn mul_assign(&mut self, rhs: T) {
for a in self.coords.iter_mut() {
*a *= rhs;
}
}
}
#[cfg(test)]
mod tests {
use super::*;
use quickcheck::{Arbitrary, Gen};
impl Arbitrary for Vec2<f64> {
fn arbitrary<G: Gen>(g: &mut G) -> Vec2<f64> {
Vec2::new(f64::arbitrary(g), f64::arbitrary(g))
}
}
#[test]
fn x_returns_first_element() {
let target = Vec2::new(1.0, 2.0);
assert!(target.x() == 1.0);
}
#[test]
fn y_returns_second_element() {
let target = Vec2::new(1.0, 2.0);
assert!(target.y() == 2.0);
}
#[test]
fn dot_product_returns_correct_result() {
let a = Vec2::new(1.0, 2.0);
let b = Vec2::new(4.0, 5.0);
assert!(a.dot(&b) == 14.0);
}
#[test]
fn add_returns_correct_result() {
let a = Vec2::new(1.0, 2.0);
let b = Vec2::new(4.0, 5.0);
let c = Vec2::new(5.0, 7.0);
assert!(a + b == c);
}
#[test]
fn add_assign_returns_correct_result() {
let mut a = Vec2::new(1.0, 2.0);
let b = Vec2::new(4.0, 5.0);
let c = Vec2::new(5.0, 7.0);
a += b;
assert!(a == c);
}
#[test]
fn sub_returns_correct_result() {
let a = Vec2::new(1.0, 2.0);
let b = Vec2::new(4.0, 6.0);
let c = Vec2::new(-3.0, -4.0);
assert!(a - b == c);
}
#[test]
fn sub_assign_returns_correct_result() {
let mut a = Vec2::new(1.0, 2.0);
let b = Vec2::new(4.0, 6.0);
let c = Vec2::new(-3.0, -4.0);
a -= b;
assert!(a == c);
}
#[test]
fn mul_by_scalar_returns_correct_result() {
let a = Vec2::new(1.0, 2.0);
let b = 0.5;
let c = Vec2::new(0.5, 1.0);
assert!(a * b == c);
}
#[test]
fn mul_assign_by_scalar_returns_correct_result() {
let mut a = Vec2::new(1.0, 2.0);
let b = 0.5;
let c = Vec2::new(0.5, 1.0);
a *= b;
assert!(a == c);
}
}

537
src/math/vec3.rs
View File

@ -1,537 +0,0 @@
use super::{Float, Mat3};
use itertools::izip;
use std::ops::{Add, AddAssign, Div, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign};
#[derive(Copy, Clone, PartialEq, Debug, Default)]
pub struct Vec3<T: Float> {
pub coords: [T; 3],
}
impl<T: Float> Vec3<T> {
pub fn new(x: T, y: T, z: T) -> Self {
Vec3 { coords: [x, y, z] }
}
pub fn from_slice(v: &[T]) -> Self {
let mut coords = [T::zero(); 3];
coords.clone_from_slice(v);
Vec3 { coords }
}
pub fn zeros() -> Vec3<T> {
Vec3 {
coords: [T::zero(), T::zero(), T::zero()],
}
}
pub fn unit_x() -> Vec3<T> {
Vec3 {
coords: [T::one(), T::zero(), T::zero()],
}
}
pub fn unit_y() -> Vec3<T> {
Vec3 {
coords: [T::zero(), T::one(), T::zero()],
}
}
pub fn unit_z() -> Vec3<T> {
Vec3 {
coords: [T::zero(), T::zero(), T::one()],
}
}
pub fn x(&self) -> T {
self.coords[0]
}
pub fn y(&self) -> T {
self.coords[1]
}
pub fn z(&self) -> T {
self.coords[2]
}
pub fn as_slice(&self) -> &[T] {
&self.coords
}
pub fn dot(&self, rhs: &Vec3<T>) -> T {
self.coords
.iter()
.copied()
.zip(rhs.coords.iter().copied())
.map(|(a_elem, b_elem)| a_elem * b_elem)
.sum()
}
pub fn cross(&self, rhs: &Vec3<T>) -> Vec3<T> {
let x = self.y() * rhs.z() - self.z() * rhs.y();
let y = self.z() * rhs.x() - self.x() * rhs.z();
let z = self.x() * rhs.y() - self.y() * rhs.x();
Vec3 { coords: [x, y, z] }
}
pub fn abs(&self) -> Self {
Vec3::new(self.x().abs(), self.y().abs(), self.z().abs())
}
pub fn norm_squared(&self) -> T {
self.dot(self)
}
pub fn norm(&self) -> T {
self.norm_squared().sqrt()
}
pub fn normalize(&self) -> Self {
let mut coords = [T::zero(); 3];
let inverse_norm = T::one() / self.norm();
for (r, a) in coords.iter_mut().zip(self.coords.iter().copied()) {
*r = a * inverse_norm;
}
Vec3 { coords }
}
pub fn smallest_coord(&self) -> usize {
let x = self.x().abs();
let y = self.y().abs();
let z = self.z().abs();
if x < y {
if x < z {
0
} else {
2
}
} else if y < z {
1
} else {
2
}
}
pub fn component_mul(&self, rhs: &Self) -> Self {
let mut coords = [T::zero(); 3];
for (elem, lhs_elem, rhs_elem) in izip!(
coords.iter_mut(),
self.coords.iter().copied(),
rhs.coords.iter().copied()
) {
*elem = lhs_elem * rhs_elem;
}
Vec3 { coords }
}
}
impl<T: Float> Index<usize> for Vec3<T> {
type Output = T;
fn index(&self, i: usize) -> &T {
&self.coords[i]
}
}
impl<T: Float> IndexMut<usize> for Vec3<T> {
fn index_mut(&mut self, i: usize) -> &mut T {
&mut self.coords[i]
}
}
impl<T: Float> Add<Vec3<T>> for &Vec3<T> {
type Output = Vec3<T>;
fn add(self, rhs: Vec3<T>) -> Vec3<T> {
let mut coords = [T::zero(); 3];
for (r, a, b) in izip!(
coords.iter_mut(),
self.coords.iter().copied(),
rhs.coords.iter().copied()
) {
*r = a + b;
}
Vec3 { coords }
}
}
impl<T: Float> Add<&Vec3<T>> for &Vec3<T> {
type Output = Vec3<T>;
fn add(self, rhs: &Vec3<T>) -> Vec3<T> {
let mut coords = [T::zero(); 3];
for (r, a, b) in izip!(
coords.iter_mut(),
self.coords.iter().copied(),
rhs.coords.iter().copied()
) {
*r = a + b;
}
Vec3 { coords }
}
}
impl<T: Float> Add for Vec3<T> {
type Output = Self;
fn add(self, rhs: Self) -> Self {
let mut coords = [T::zero(); 3];
for (r, a, b) in izip!(
coords.iter_mut(),
self.coords.iter().copied(),
rhs.coords.iter().copied()
) {
*r = a + b;
}
Vec3 { coords }
}
}
impl<T: Float> AddAssign for Vec3<T> {
fn add_assign(&mut self, rhs: Self) {
for (a, b) in self.coords.iter_mut().zip(rhs.coords.iter().copied()) {
*a += b;
}
}
}
impl<T: Float> Neg for Vec3<T> {
type Output = Vec3<T>;
fn neg(self) -> Vec3<T> {
Vec3::new(-self.x(), -self.y(), -self.z())
}
}
impl<T: Float> Sub for &Vec3<T> {
type Output = Vec3<T>;
fn sub(self, rhs: Self) -> Vec3<T> {
let mut coords = [T::zero(); 3];
for (r, a, b) in izip!(
coords.iter_mut(),
self.coords.iter().copied(),
rhs.coords.iter().copied()
) {
*r = a - b;
}
Vec3 { coords }
}
}
impl<T: Float> Sub for Vec3<T> {
type Output = Vec3<T>;
fn sub(self, rhs: Self) -> Vec3<T> {
let mut coords = [T::zero(); 3];
for (r, a, b) in izip!(
coords.iter_mut(),
self.coords.iter().copied(),
rhs.coords.iter().copied()
) {
*r = a - b;
}
Vec3 { coords }
}
}
impl<T: Float> SubAssign for Vec3<T> {
fn sub_assign(&mut self, rhs: Self) {
for (a, b) in self.coords.iter_mut().zip(rhs.coords.iter().copied()) {
*a -= b;
}
}
}
impl<T: Float> Mul<T> for &Vec3<T> {
type Output = Vec3<T>;
fn mul(self, rhs: T) -> Vec3<T> {
let mut coords = [T::zero(); 3];
for (r, a) in coords.iter_mut().zip(self.coords.iter().copied()) {
*r = a * rhs;
}
Vec3 { coords }
}
}
impl<T: Float> Mul<T> for Vec3<T> {
type Output = Vec3<T>;
fn mul(self, rhs: T) -> Vec3<T> {
let mut coords = [T::zero(); 3];
for (r, a) in coords.iter_mut().zip(self.coords.iter().copied()) {
*r = a * rhs;
}
Vec3 { coords }
}
}
impl<T: Float> MulAssign<T> for Vec3<T> {
fn mul_assign(&mut self, rhs: T) {
for a in self.coords.iter_mut() {
*a *= rhs;
}
}
}
impl<T:Float> Mul<Mat3<T>> for &Vec3<T> {
type Output = Vec3<T>;
fn mul(self, rhs: Mat3<T>) -> Vec3<T> {
let mut coords = [T::zero(); 3];
for i in 0..3 {
coords[i] = self.dot(&rhs.get_column(i));
}
Vec3 { coords }
}
}
impl<T:Float> Mul<Mat3<T>> for Vec3<T> {
type Output = Self;
fn mul(self, rhs: Mat3<T>) -> Self {
let mut coords = [T::zero(); 3];
for i in 0..3 {
coords[i] = self.dot(&rhs.get_column(i));
}
Vec3 { coords }
}
}
impl<T:Float> MulAssign<Mat3<T>> for Vec3<T> {
fn mul_assign(&mut self, rhs: Mat3<T>) {
let mut coords = [T::zero(); 3];
for i in 0..3 {
coords[i] = self.dot(&rhs.get_column(i));
}
self.coords = coords;
}
}
impl Mul<Vec3<f64>> for f64 {
type Output = Vec3<f64>;
fn mul(self, rhs: Vec3<f64>) -> Vec3<f64> {
rhs * self
}
}
impl Mul<&Vec3<f64>> for f64 {
type Output = Vec3<f64>;
fn mul(self, rhs: &Vec3<f64>) -> Vec3<f64> {
rhs * self
}
}
impl<T: Float> Div<T> for &Vec3<T> {
type Output = Vec3<T>;
fn div(self, rhs: T) -> Vec3<T> {
let mut coords = [T::zero(); 3];
for (r, a) in coords.iter_mut().zip(self.coords.iter().copied()) {
*r = a / rhs;
}
Vec3 { coords }
}
}
impl<T: Float> Div<T> for Vec3<T> {
type Output = Vec3<T>;
fn div(self, rhs: T) -> Vec3<T> {
let mut coords = [T::zero(); 3];
for (r, a) in coords.iter_mut().zip(self.coords.iter().copied()) {
*r = a / rhs;
}
Vec3 { coords }
}
}
#[cfg(test)]
mod tests {
use super::*;
use quickcheck::{Arbitrary, Gen};
impl<T: Arbitrary + Float> Arbitrary for Vec3<T> {
fn arbitrary<G: Gen>(g: &mut G) -> Vec3<T> {
Vec3::new(T::arbitrary(g), T::arbitrary(g), T::arbitrary(g))
}
}
#[test]
fn x_returns_first_element() {
let target = Vec3::new(1.0, 2.0, 3.0);
assert!(target.x() == 1.0);
}
#[test]
fn y_returns_second_element() {
let target = Vec3::new(1.0, 2.0, 3.0);
assert!(target.y() == 2.0);
}
#[test]
fn z_returns_third_element() {
let target = Vec3::new(1.0, 2.0, 3.0);
assert!(target.z() == 3.0);
}
/*#[test]
fn from_iterator_takes_first_three_elements() {
let target = Vec3::from_iterator([1.0, 2.0, 3.0].iter());
assert!(target = Vec3::new(1.0, 2.0, 3.0));
}*/
#[test]
fn dot_product_returns_correct_result() {
let a = Vec3::new(1.0, 2.0, 3.0);
let b = Vec3::new(4.0, 5.0, 6.0);
assert!(a.dot(&b) == 32.0);
}
#[test]
fn cross_product_returns_correct_result() {
let a = Vec3::new(1.0, 2.0, 3.0);
let b = Vec3::new(4.0, 5.0, 6.0);
let c = Vec3::new(-3.0, 6.0, -3.0);
assert!(a.cross(&b) == c);
}
#[test]
fn norm_returns_expected_value() {
let target = Vec3::new(2.0, 3.0, 6.0);
assert!(target.norm() == 7.0);
}
#[test]
fn normalized_vector_times_norm_yields_original() {
let mut target = Vec3::new(2.0, 3.0, 6.0);
let norm = target.norm();
target = target.normalize();
target *= norm;
assert!(target == Vec3::new(2.0, 3.0, 6.0));
}
#[test]
fn smallest_coord_works_for_x_when_positive() {
let target = Vec3::new(1.0, 2.0, 3.0);
assert!(target.smallest_coord() == 0);
}
#[test]
fn smallest_coord_works_for_x_when_negative() {
let target = Vec3::new(-2.0, -3.0, 3.0);
assert!(target.smallest_coord() == 0);
}
#[test]
fn smallest_coord_works_for_y_when_positive() {
let target = Vec3::new(2.0, 1.0, 3.0);
assert!(target.smallest_coord() == 1);
}
#[test]
fn smallest_coord_works_for_y_when_negative() {
let target = Vec3::new(-3.0, -2.0, 3.0);
assert!(target.smallest_coord() == 1);
}
#[test]
fn smallest_coord_works_for_z_when_positive() {
let target = Vec3::new(3.0, 2.0, 1.0);
assert!(target.smallest_coord() == 2);
}
#[test]
fn smallest_coord_works_for_z_when_negative() {
let target = Vec3::new(3.0, -3.0, -2.0);
assert!(target.smallest_coord() == 2);
}
#[test]
fn add_returns_correct_result() {
let a = Vec3::new(1.0, 2.0, 3.0);
let b = Vec3::new(4.0, 5.0, 6.0);
let c = Vec3::new(5.0, 7.0, 9.0);
assert!(a + b == c);
}
#[test]
fn add_assign_returns_correct_result() {
let mut a = Vec3::new(1.0, 2.0, 3.0);
let b = Vec3::new(4.0, 5.0, 6.0);
let c = Vec3::new(5.0, 7.0, 9.0);
a += b;
assert!(a == c);
}
#[test]
fn sub_returns_correct_result() {
let a = Vec3::new(1.0, 2.0, 3.0);
let b = Vec3::new(4.0, 6.0, 8.0);
let c = Vec3::new(-3.0, -4.0, -5.0);
assert!(a - b == c);
}
#[test]
fn sub_assign_returns_correct_result() {
let mut a = Vec3::new(1.0, 2.0, 3.0);
let b = Vec3::new(4.0, 6.0, 8.0);
let c = Vec3::new(-3.0, -4.0, -5.0);
a -= b;
assert!(a == c);
}
#[test]
fn mul_by_scalar_returns_correct_result() {
let a = Vec3::new(1.0, 2.0, 3.0);
let b = 0.5;
let c = Vec3::new(0.5, 1.0, 1.5);
assert!(a * b == c);
}
#[test]
fn div_by_scalar_returns_correct_result() {
let a = Vec3::new(1.0, 2.0, 3.0);
let b = 2.0;
let c = Vec3::new(0.5, 1.0, 1.5);
assert!(dbg!(a / b) == dbg!(c));
}
#[test]
fn mul_assign_by_scalar_returns_correct_result() {
let mut a = Vec3::new(1.0, 2.0, 3.0);
let b = 0.5;
let c = Vec3::new(0.5, 1.0, 1.5);
a *= b;
assert!(a == c);
}
#[test]
fn mul_with_mat3_returns_expected_result() {
let a = Mat3::from_rows(
&Vec3::new(1.0, 2.0, 3.0),
&Vec3::new(4.0, 5.0, 6.0),
&Vec3::new(7.0, 8.0, 9.0),
);
let b = Vec3::new(10.0, 11.0, 12.0);
let c = Vec3::new(138.0, 171.0, 204.0);
assert!(b * a == c);
}
#[test]
fn mul_assign_with_mat3_returns_expected_result() {
let a = Mat3::from_rows(
&Vec3::new(1.0, 2.0, 3.0),
&Vec3::new(4.0, 5.0, 6.0),
&Vec3::new(7.0, 8.0, 9.0),
);
let mut b = Vec3::new(10.0, 11.0, 12.0);
let c = Vec3::new(138.0, 171.0, 204.0);
b *= a;
assert!(b == c);
}
}

212
src/math/vec4.rs
View File

@ -1,212 +0,0 @@
use super::{Float, Vec3};
use itertools::izip;
use std::ops::{Add, AddAssign, Mul, MulAssign, Sub, SubAssign};
#[derive(PartialEq, Debug, Clone, Copy)]
pub struct Vec4<T: Float> {
pub coords: [T; 4],
}
impl<T: Float> Vec4<T> {
pub fn new(x: T, y: T, z: T, w: T) -> Self {
Vec4 {
coords: [x, y, z, w],
}
}
pub fn x(&self) -> T {
self.coords[0]
}
pub fn y(&self) -> T {
self.coords[1]
}
pub fn z(&self) -> T {
self.coords[2]
}
pub fn w(&self) -> T {
self.coords[3]
}
pub fn xyz(&self) -> Vec3<T> {
let mut coords = [T::zero(); 3];
coords.copy_from_slice(&self.coords[0..3]);
Vec3 { coords }
}
pub fn dot(&self, rhs: &Vec4<T>) -> T {
self.coords
.iter()
.copied()
.zip(rhs.coords.iter().copied())
.map(|(a_elem, b_elem)| a_elem * b_elem)
.sum()
}
}
impl<T: Float> Add for Vec4<T> {
type Output = Self;
fn add(self, rhs: Self) -> Self {
let mut coords = [T::zero(); 4];
for (r, a, b) in izip!(
coords.iter_mut(),
self.coords.iter().copied(),
rhs.coords.iter().copied()
) {
*r = a + b;
}
Vec4 { coords }
}
}
impl<T: Float> AddAssign for Vec4<T> {
fn add_assign(&mut self, rhs: Self) {
for (a, b) in self.coords.iter_mut().zip(rhs.coords.iter().copied()) {
*a += b;
}
}
}
impl<T: Float> Sub for Vec4<T> {
type Output = Self;
fn sub(self, rhs: Self) -> Self {
let mut coords = [T::zero(); 4];
for (r, a, b) in izip!(
coords.iter_mut(),
self.coords.iter().copied(),
rhs.coords.iter().copied()
) {
*r = a - b;
}
Vec4 { coords }
}
}
impl<T: Float> SubAssign for Vec4<T> {
fn sub_assign(&mut self, rhs: Self) {
for (a, b) in self.coords.iter_mut().zip(rhs.coords.iter().copied()) {
*a -= b;
}
}
}
impl<T: Float> Mul<T> for Vec4<T> {
type Output = Self;
fn mul(self, rhs: T) -> Vec4<T> {
let mut coords = [T::zero(); 4];
for (r, a) in coords.iter_mut().zip(self.coords.iter().copied()) {
*r = a * rhs;
}
Vec4 { coords }
}
}
impl<T: Float> MulAssign<T> for Vec4<T> {
fn mul_assign(&mut self, rhs: T) {
for a in self.coords.iter_mut() {
*a *= rhs;
}
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn x_returns_first_element() {
let target = Vec4::new(1.0, 2.0, 3.0, 4.0);
assert!(target.x() == 1.0);
}
#[test]
fn y_returns_second_element() {
let target = Vec4::new(1.0, 2.0, 3.0, 4.0);
assert!(target.y() == 2.0);
}
#[test]
fn z_returns_third_element() {
let target = Vec4::new(1.0, 2.0, 3.0, 4.0);
assert!(target.z() == 3.0);
}
#[test]
fn w_returns_third_element() {
let target = Vec4::new(1.0, 2.0, 3.0, 4.0);
assert!(target.w() == 4.0);
}
#[test]
fn xyz_returns_expected_value() {
let target = Vec4::new(1.0, 2.0, 3.0, 4.0).xyz();
assert!(target.x() == 1.0);
assert!(target.y() == 2.0);
assert!(target.z() == 3.0);
}
#[test]
fn dot_product_returns_correct_result() {
let a = Vec4::new(1.0, 2.0, 3.0, 4.0);
let b = Vec4::new(4.0, 5.0, 6.0, 7.0);
assert!(a.dot(&b) == 60.0);
}
#[test]
fn add_returns_correct_result() {
let a = Vec4::new(1.0, 2.0, 3.0, 4.0);
let b = Vec4::new(4.0, 5.0, 6.0, 7.0);
let c = Vec4::new(5.0, 7.0, 9.0, 11.0);
assert!(a + b == c);
}
#[test]
fn add_assign_returns_correct_result() {
let mut a = Vec4::new(1.0, 2.0, 3.0, 4.0);
let b = Vec4::new(4.0, 5.0, 6.0, 7.0);
let c = Vec4::new(5.0, 7.0, 9.0, 11.0);
a += b;
assert!(a == c);
}
#[test]
fn sub_returns_correct_result() {
let a = Vec4::new(1.0, 2.0, 3.0, 4.0);
let b = Vec4::new(4.0, 6.0, 8.0, 10.0);
let c = Vec4::new(-3.0, -4.0, -5.0, -6.0);
assert!(a - b == c);
}
#[test]
fn sub_assign_returns_correct_result() {
let mut a = Vec4::new(1.0, 2.0, 3.0, 4.0);
let b = Vec4::new(4.0, 6.0, 8.0, 10.0);
let c = Vec4::new(-3.0, -4.0, -5.0, -6.0);
a -= b;
assert!(a == c);
}
#[test]
fn mul_by_scalar_returns_correct_result() {
let a = Vec4::new(1.0, 2.0, 3.0, 4.0);
let b = 0.5;
let c = Vec4::new(0.5, 1.0, 1.5, 2.0);
assert!(a * b == c);
}
#[test]
fn mul_assign_by_scalar_returns_correct_result() {
let mut a = Vec4::new(1.0, 2.0, 3.0, 4.0);
let b = 0.5;
let c = Vec4::new(0.5, 1.0, 1.5, 2.0);
a *= b;
assert!(a == c);
}
}

91
src/mesh.rs
View File

@ -1,91 +0,0 @@
/// Load a model from a Wavefront .obj file
mod wavefront_obj {
use crate::materials::Material;
use crate::math::Vec3;
use crate::raycasting::{Primitive, Triangle};
use obj::{IndexTuple, Obj, SimplePolygon};
use std::io::Result;
use std::path::Path;
use std::sync::Arc;
fn get_vertex_and_normal(
index_tuple: &IndexTuple,
vertex_positions: &[[f32; 3]],
normal_positions: &[[f32; 3]],
) -> (Vec3<f64>, Vec3<f64>) {
let &IndexTuple(vertex_index, _, maybe_normal_index) = index_tuple;
(
{
let vertex_coords = &vertex_positions[vertex_index];
Vec3::new(
vertex_coords[0] as f64,
vertex_coords[1] as f64,
vertex_coords[2] as f64,
)
},
match maybe_normal_index {
Some(normal_index) => {
let normal_coords = &normal_positions[normal_index];
Vec3::new(
normal_coords[0] as f64,
normal_coords[1] as f64,
normal_coords[2] as f64,
)
}
None => Vec3::zeros(),
},
)
}
fn get_triangles(
polygon: &SimplePolygon,
vertex_positions: &[[f32; 3]],
normal_positions: &[[f32; 3]],
material: Arc<dyn Material>,
) -> Vec<Triangle> {
if let Some(v0_index) = polygon.iter().next() {
let (v0_vertex, v0_normal) =
get_vertex_and_normal(v0_index, vertex_positions, normal_positions);
polygon
.iter()
.skip(1)
.zip(polygon.iter().skip(2))
.map(|(v1_index, v2_index)| {
let (v1_vertex, v1_normal) =
get_vertex_and_normal(v1_index, vertex_positions, normal_positions);
let (v2_vertex, v2_normal) =
get_vertex_and_normal(v2_index, vertex_positions, normal_positions);
let vertices = [v0_vertex, v1_vertex, v2_vertex];
let normals = [v0_normal, v1_normal, v2_normal];
Triangle {
vertices,
normals,
material: material.clone(),
}
})
.collect()
} else {
vec![]
}
}
pub fn load_obj(
filename: &Path,
material: Arc<dyn Material>,
) -> Result<Vec<Arc<dyn Primitive>>> {
let obj = Obj::<SimplePolygon>::load(filename)?;
Ok(obj
.objects
.iter()
.flat_map(|object| object.groups.iter())
.flat_map(|group| group.polys.iter())
.flat_map(|poly| get_triangles(poly, &obj.position, &obj.normal, material.clone()))
.map(|triangle| Arc::new(triangle) as Arc<dyn Primitive>)
.collect())
}
}
pub use wavefront_obj::load_obj;

60
src/random_distributions/cosine_weighted_hemisphere.rs
View File

@ -1,60 +0,0 @@
use std::f64::consts::PI;
use crate::math::Vec3;
use super::{RandomDistribution, UnitDisc};
#[derive(Default)]
pub struct CosineWeightedHemisphere {
unit_disc: UnitDisc,
}
impl CosineWeightedHemisphere {
pub fn new() -> CosineWeightedHemisphere {
let unit_disc = UnitDisc::new();
CosineWeightedHemisphere { unit_disc }
}
}
impl RandomDistribution<Vec3<f64>> for CosineWeightedHemisphere {
fn value(&self) -> Vec3<f64> {
let point_on_disc = self.unit_disc.value();
let z = 0.0f64
.max(
1.0 - point_on_disc.x() * point_on_disc.x() - point_on_disc.y() * point_on_disc.y(),
)
.sqrt();
Vec3::new(point_on_disc.x(), point_on_disc.y(), z)
}
fn pdf(&self, v: Vec3<f64>) -> f64 {
(v.x() * v.x() + v.y() * v.y()).sqrt() / PI
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
#[ignore]
fn print_values() {
let target = CosineWeightedHemisphere::new();
for _ in 0..1000 {
let value = target.value();
println!("{}, {}, {}", value.x(), value.y(), value.z());
}
}
#[test]
#[ignore]
fn integral_is_near_area() {
let target = CosineWeightedHemisphere::new();
let integral = (0..100000)
.map(|_| target.value())
.map(|value| 1.0 / target.pdf(value))
.sum::<f64>()
/ 100000.0;
println!("Area: {}\nIntegral: {}", 2.0 * PI, integral);
}
}

67
src/random_distributions/linear_weighted.rs
View File

@ -1,67 +0,0 @@
use rand::distributions::Open01;
use rand::{thread_rng, Rng};
use super::RandomDistribution;
pub struct LinearWeighted {
max_value: f64,
}
impl LinearWeighted {
pub fn new(max_value: f64) -> LinearWeighted {
LinearWeighted { max_value }
}
}
impl RandomDistribution<f64> for LinearWeighted {
fn value(&self) -> f64 {
let mut rng = thread_rng();
rng.sample::<f64, _>(Open01).sqrt() * self.max_value
}
fn pdf(&self, value: f64) -> f64 {
2.0 * value / (self.max_value * self.max_value)
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
#[ignore]
fn print_values() {
let target = LinearWeighted::new(2.0);
for _ in 0..1000 {
let value = target.value();
println!("{}", value);
}
}
#[test]
#[ignore]
fn print_buckets() {
let mut buckets = [0; 20];
let target = LinearWeighted::new(20.0);
for _ in 0..10000 {
let value = target.value();
let i = value as usize;
buckets[i] += 1;
}
for count in buckets {
println!("{}", count);
}
}
#[test]
#[ignore]
fn integral_is_near_area() {
let target = LinearWeighted::new(2.0);
let integral = (0..100000)
.map(|_| target.value())
.map(|value| 1.0 / target.pdf(value))
.sum::<f64>()
/ 100000.0;
println!("Area: {}\nIntegral: {}", 2.0, integral);
}
}

22
src/random_distributions/mod.rs
View File

@ -1,22 +0,0 @@
mod uniform_square;
pub use uniform_square::UniformSquare;
mod unit_disc;
pub use unit_disc::UnitDisc;
mod uniform_hemisphere;
pub use uniform_hemisphere::UniformHemisphere;
mod cosine_weighted_hemisphere;
pub use cosine_weighted_hemisphere::CosineWeightedHemisphere;
mod linear_weighted;
pub use linear_weighted::LinearWeighted;
mod sky_light_pdf;
pub use sky_light_pdf::SkyLightPdf;
pub trait RandomDistribution<T> {
fn value(&self) -> T;
fn pdf(&self, value: T) -> f64;
}

71
src/random_distributions/sky_light_pdf.rs
View File

@ -1,71 +0,0 @@
use std::f64::consts::PI;
use rand::distributions::Open01;
use rand::{thread_rng, Rng};
use crate::math::Vec3;
use super::{LinearWeighted, RandomDistribution};
pub struct SkyLightPdf {
z_distribution: LinearWeighted,
}
impl SkyLightPdf {
pub fn new() -> SkyLightPdf {
let z_distribution = LinearWeighted::new(1.0);
SkyLightPdf { z_distribution }
}
}
impl Default for SkyLightPdf {
fn default() -> SkyLightPdf {
SkyLightPdf::new()
}
}
impl RandomDistribution<Vec3<f64>> for SkyLightPdf {
fn value(&self) -> Vec3<f64> {
let mut rng = thread_rng();
let phi = rng.sample::<f64, _>(Open01) * 2.0 * PI;
let z = self.z_distribution.value();
let r = (1.0 - z * z).sqrt();
Vec3::new(r * phi.cos(), r * phi.sin(), z)
}
fn pdf(&self, value: Vec3<f64>) -> f64 {
let z = value.z();
if z < 0.0 {
0.0
} else {
z / PI
}
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
#[ignore]
fn print_values() {
let target = SkyLightPdf::new();
for _ in 0..1000 {
let value = target.value();
println!("{}, {}, {}", value.x(), value.y(), value.z());
}
}
#[test]
#[ignore]
fn integral_is_near_area() {
let target = SkyLightPdf::new();
let integral = (0..100000)
.map(|_| target.value())
.map(|value| 1.0 / target.pdf(value))
.sum::<f64>()
/ 100000.0;
println!("Area: {}\nIntegral: {}", 2.0 * PI, integral);
}
}

67
src/random_distributions/uniform_hemisphere.rs
View File

@ -1,67 +0,0 @@
use std::f64::consts::PI;
use rand::distributions::{Open01, OpenClosed01};
use rand::{thread_rng, Rng};
use crate::math::Vec3;
use super::RandomDistribution;
#[derive(Default)]
pub struct UniformHemisphere {}
impl UniformHemisphere {
pub fn new() -> UniformHemisphere {
UniformHemisphere {}
}
}
impl RandomDistribution<Vec3<f64>> for UniformHemisphere {
fn value(&self) -> Vec3<f64> {
let mut rng = thread_rng();
let mut result = Vec3::new(
2.0 * rng.sample::<f64, _>(Open01) - 1.0,
2.0 * rng.sample::<f64, _>(Open01) - 1.0,
rng.sample::<f64, _>(OpenClosed01),
);
while result.norm_squared() > 1.0 {
result = Vec3::new(
2.0 * rng.sample::<f64, _>(Open01) - 1.0,
2.0 * rng.sample::<f64, _>(Open01) - 1.0,
rng.sample::<f64, _>(OpenClosed01),
);
}
result.normalize()
}
fn pdf(&self, _: Vec3<f64>) -> f64 {
1.0 / (2.0 * PI)
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
#[ignore]
fn print_values() {
let target = UniformHemisphere::new();
for _ in 0..1000 {
let value = target.value();
println!("{}, {}, {}", value.x(), value.y(), value.z());
}
}
#[test]
#[ignore]
fn integral_is_near_area() {
let target = UniformHemisphere::new();
let integral = (0..1000)
.map(|_| target.value())
.map(|value| 1.0 / target.pdf(value))
.sum::<f64>()
/ 1000.0;
println!("Area: {}\nIntegral: {}", 2.0 * PI, integral);
}
}

63
src/random_distributions/uniform_square.rs
View File

@ -1,63 +0,0 @@
use rand::distributions::Open01;
use rand::{thread_rng, Rng};
use crate::math::Vec2;
use super::RandomDistribution;
#[derive(Debug)]
pub struct UniformSquare {
corner: Vec2<f64>,
size: f64,
}
impl UniformSquare {
pub fn new(corner: Vec2<f64>, size: f64) -> UniformSquare {
UniformSquare { corner, size }
}
}
impl RandomDistribution<Vec2<f64>> for UniformSquare {
fn value(&self) -> Vec2<f64> {
let mut rng = thread_rng();
self.corner
+ Vec2::new(rng.sample::<f64, _>(Open01), rng.sample::<f64, _>(Open01)) * self.size
}
fn pdf(&self, _value: Vec2<f64>) -> f64 {
1.0 / (self.size * self.size)
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
#[ignore]
fn print_values() {
let target = UniformSquare {
corner: Vec2::new(1.5, -2.5),
size: 3.0,
};
for _ in 0..1000 {
let value = target.value();
println!("{}, {}", value.x(), value.y());
}
}
#[test]
#[ignore]
fn integral_is_near_area() {
let target = UniformSquare {
corner: Vec2::new(1.5, -2.5),
size: 3.0,
};
let integral = (0..1000)
.map(|_| target.value())
.map(|value| 1.0 / target.pdf(value))
.sum::<f64>()
/ 1000.0;
println!("Area: {}\nIntegral: {}", 3.0 * 3.0, integral);
}
}

72
src/random_distributions/unit_disc.rs
View File

@ -1,72 +0,0 @@
use std::f64::consts::PI;
use crate::math::Vec2;
use super::{RandomDistribution, UniformSquare};
#[derive(Debug)]
pub struct UnitDisc {
square_distribution: UniformSquare,
}
impl Default for UnitDisc {
fn default() -> UnitDisc {
UnitDisc::new()
}
}
impl UnitDisc {
pub fn new() -> UnitDisc {
let square_distribution = UniformSquare::new(Vec2::new(-1.0, -1.0), 2.0);
UnitDisc {
square_distribution,
}
}
}
impl RandomDistribution<Vec2<f64>> for UnitDisc {
fn value(&self) -> Vec2<f64> {
let offset = self.square_distribution.value();
if offset.x() == 0.0 && offset.y() == 0.0 {
offset
} else {
let (radius, angle) = if offset.x().abs() > offset.y().abs() {
(offset.x(), (PI / 4.0) * offset.y() / offset.x())
} else {
(offset.y(), PI / 2.0 - (PI / 4.0) * offset.x() / offset.y())
};
Vec2::new(angle.cos(), angle.sin()) * radius
}
}
fn pdf(&self, _: Vec2<f64>) -> f64 {
1.0 / PI
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
#[ignore]
fn print_values() {
let target = UnitDisc::new();
for _ in 0..1000 {
let value = target.value();
println!("{}, {}", value.x(), value.y());
}
}
#[test]
#[ignore]
fn integral_is_near_area() {
let target = UnitDisc::new();
let integral = (0..1000)
.map(|_| target.value())
.map(|value| 1.0 / target.pdf(value))
.sum::<f64>()
/ 1000.0;
println!("Area: {}\nIntegral: {}", PI, integral);
}
}

375
src/raycasting.rs Normal file
View File

@ -0,0 +1,375 @@
use nalgebra::{convert, RealField, Vector3};
use super::materials::Material;
use std::rc::Rc;
#[derive(Clone, Debug)]
pub struct Ray<T: RealField> {
origin: Vector3<T>,
direction: Vector3<T>,
}
impl<T: RealField> Ray<T> {
pub fn new(origin: Vector3<T>, direction: Vector3<T>) -> Ray<T> {
Ray {
origin,
direction: direction.normalize(),
}
}
pub fn point_at(&self, t: T) -> Vector3<T> {
return self.origin + self.direction * t;
}
pub fn bias(&self, amount: T) -> Ray<T> {
Ray::new(self.origin + self.direction * amount, self.direction)
}
}
#[derive(Debug)]
pub struct IntersectionInfo<T: RealField> {
pub distance: T,
pub location: Vector3<T>,
pub normal: Vector3<T>,
pub tangent: Vector3<T>,
pub cotangent: Vector3<T>,
pub retro: Vector3<T>,
pub material: Rc<dyn Material<T>>,
}
pub trait Intersect<T: RealField> {
fn intersect<'a>(&'a self, ray: &Ray<T>) -> Option<IntersectionInfo<T>>;
}
pub struct Sphere<T: RealField> {
centre: Vector3<T>,
radius: T,
material: Rc<dyn Material<T>>,
}
impl<T: RealField> Sphere<T> {
pub fn new(centre: Vector3<T>, radius: T, material: Rc<dyn Material<T>>) -> Sphere<T> {
Sphere {
centre,
radius,
material,
}
}
}
impl<T: RealField> Intersect<T> for Sphere<T> {
fn intersect<'a>(&'a self, ray: &Ray<T>) -> Option<IntersectionInfo<T>> {
/*let ray_origin_to_sphere_centre = self.centre - ray.origin;
let radius_squared = self.radius * self.radius;
let is_inside_sphere = ray_origin_to_sphere_centre.norm_squared() <= radius_squared;
// t0/p0 is the point on the ray that's closest to the centre of the sphere
// ray.direction is normalized, so it's not necessary to divide by its length.
let t0 = ray_origin_to_sphere_centre.dot(&ray.direction);
if !is_inside_sphere && t0 < T::zero() {
// Sphere is behind ray origin
return None;
}
// Squared distance between ray origin and sphere centre
let d0_squared = (ray_origin_to_sphere_centre).norm_squared();
// p0, ray.origin and sphere.centre form a right triangle, with p0 at the right corner,
// Squared distance petween p0 and sphere centre, using Pythagoras
let p0_dist_from_centre_squared = d0_squared - t0 * t0;
if p0_dist_from_centre_squared > radius_squared {
// Sphere is in front of ray but ray misses
return None;
}
let p0_dist_from_centre =p0_dist_from_centre_squared.sqrt();
// Two more right triangles are formed by p0, the sphere centre, and the two places
// where the ray intersects the sphere. (Or the ray may be a tangent to the sphere
// in which case these triangles are degenerate. Here we use Pythagoras again to find
.// find the distance between p0 and the two intersection points.
let delta = (radius_squared - p0_dist_from_centre_squared).sqrt();
let distance = if is_inside_sphere {
// radius origin is inside sphere
t0 + delta
} else {
t0 - delta
};
let location = ray.point_at(distance);
let normal = (location - self.centre).normalize();
let tangent = normal.cross(&Vector3::z_axis());
let cotangent = normal.cross(&tangent);
let retro = -ray.direction;*/
let two: T = convert(2.0);
let four: T = convert(4.0);
let a = ray
.direction
.component_mul(&ray.direction)
.iter()
.fold(T::zero(), |a, b| a + *b);
let b = ((ray.origin.component_mul(&ray.direction)
- self.centre.component_mul(&ray.direction))
* two)
.iter()
.fold(T::zero(), |a, b| a + *b);
let c = (ray.origin.component_mul(&ray.origin) + self.centre.component_mul(&self.centre)
- self.centre.component_mul(&ray.origin) * two)
.iter()
.fold(T::zero(), |a, b| a + *b)
- self.radius * self.radius;
let delta_squared: T = b * b - four * a * c;
if delta_squared < T::zero() {
None
} else {
let delta = delta_squared.sqrt();
let one_over_2_a = T::one() / (two * a);
let t1 = (-b - delta) * one_over_2_a;
let t2 = (-b + delta) * one_over_2_a;
let distance = if t1 < T::zero() {
t2
} else if t2 < T::zero() {
t1
} else if t1 < t2 {
t1
} else {
t2
};
if distance <= T::zero() {
None
} else {
let location = ray.point_at(distance);
let normal = (location - self.centre).normalize();
let tangent = normal.cross(&Vector3::z_axis());
let cotangent = normal.cross(&tangent);
let retro = -ray.direction;
Some(IntersectionInfo {
distance,
location,
normal,
tangent,
cotangent,
retro,
material: Rc::clone(&self.material),
})
}
}
}
}
pub struct Plane<T: RealField> {
normal: Vector3<T>,
tangent: Vector3<T>,
cotangent: Vector3<T>,
distance_from_origin: T,
material: Rc<dyn Material<T>>,
}
impl<T: RealField> Plane<T> {
pub fn new(
normal: Vector3<T>,
distance_from_origin: T,
material: Rc<dyn Material<T>>,
) -> Plane<T> {
normal.normalize();
let mut axis_closest_to_tangent = Vector3::zeros();
axis_closest_to_tangent[normal.iamin()] = T::one();
let cotangent = normal.cross(&axis_closest_to_tangent);
let tangent = normal.cross(&cotangent);
Plane {
normal,
tangent,
cotangent,
distance_from_origin,
material,
}
}
}
impl<T: RealField> Intersect<T> for Plane<T> {
fn intersect<'a>(&'a self, ray: &Ray<T>) -> Option<IntersectionInfo<T>> {
let ray_direction_dot_plane_normal = ray.direction.dot(&self.normal);
let point_on_plane = self.normal * self.distance_from_origin;
let point_on_plane_minus_ray_origin_dot_normal =
(point_on_plane - ray.origin).dot(&self.normal);
if ray_direction_dot_plane_normal == convert(0.0) {
//Ray is parallel to plane
if point_on_plane_minus_ray_origin_dot_normal != convert(0.0) {
//Ray is not in plane
return None;
}
}
let t = point_on_plane_minus_ray_origin_dot_normal / ray_direction_dot_plane_normal;
if t < convert(0.0) {
return None;
}
Some(IntersectionInfo {
distance: t,
location: ray.point_at(t),
normal: self.normal,
tangent: self.tangent,
cotangent: self.cotangent,
retro: -ray.direction,
material: Rc::clone(&self.material),
})
}
}
#[cfg(test)]
mod tests {
use quickcheck_macros::quickcheck;
macro_rules! assert_matches {
($expression:expr, $($pattern:tt)+) => {
match $expression {
$($pattern)+ => (),
ref e => panic!("assertion failed: `{:?}` does not match `{}`", e,
stringify!($($pattern)+)),
}
}
}
use super::*;
use crate::materials::LambertianMaterial;
use quickcheck::{Arbitrary, Gen, TestResult};
impl<T: Arbitrary + RealField> Arbitrary for Ray<T> {
fn arbitrary<G: Gen>(g: &mut G) -> Ray<T> {
let origin = <Vector3<T> as Arbitrary>::arbitrary(g);
let direction = <Vector3<T> as Arbitrary>::arbitrary(g);
return Ray::new(origin, direction);
}
}
#[quickcheck]
fn t0_is_origin(ray: Ray<f64>) -> bool {
ray.point_at(0.0) == ray.origin
}
#[quickcheck]
fn t1_is_origin_plus_direction(ray: Ray<f64>) -> bool {
ray.point_at(1.0) == ray.origin + ray.direction
}
#[quickcheck]
fn points_are_colinear(ray: Ray<f64>, t1: f64, t2: f64, t3: f64) -> bool {
let p1 = ray.point_at(t1);
let p2 = ray.point_at(t2);
let p3 = ray.point_at(t3);
let epsilon = [t1, t2, t3, ray.origin[0], ray.origin[1], ray.origin[2]]
.iter()
.fold(0.0, |a, &b| a.max(b.abs()))
* std::f64::EPSILON
* 256.0;
(p2 - p1).cross(&(p3 - p2)).norm() < epsilon
}
#[quickcheck]
fn t_is_distance(ray: Ray<f64>, t: f64) -> bool {
(ray.point_at(t) - ray.origin).norm() - t.abs() < 0.0000000001
}
#[test]
fn ray_intersects_sphere() {
let r = Ray::new(Vector3::new(1.0, 2.0, 3.0), Vector3::new(0.0, 0.0, 1.0));
let s = Sphere::new(
Vector3::new(1.5, 1.5, 15.0),
5.0,
Rc::new(LambertianMaterial::new_dummy()),
);
assert_matches!(s.intersect(&r), Some(_));
}
#[test]
fn ray_does_not_intersect_sphere_when_sphere_is_in_front() {
let r = Ray::new(Vector3::new(1.0, 2.0, 3.0), Vector3::new(0.0, 0.0, 1.0));
let s = Sphere::new(
Vector3::new(-5.0, 1.5, 15.0),
5.0,
Rc::new(LambertianMaterial::new_dummy()),
);
assert_matches!(s.intersect(&r), None);
}
#[test]
fn ray_does_not_intersect_sphere_when_sphere_is_behind() {
let r = Ray::new(Vector3::new(1.0, 2.0, 3.0), Vector3::new(0.0, 0.0, 1.0));
let s = Sphere::new(
Vector3::new(1.5, 1.5, -15.0),
5.0,
Rc::new(LambertianMaterial::new_dummy()),
);
assert_matches!(s.intersect(&r), None);
}
#[test]
fn ray_intersects_sphere_when_origin_is_inside() {
let r = Ray::new(Vector3::new(1.0, 2.0, 3.0), Vector3::new(0.0, 0.0, 1.0));
let s = Sphere::new(
Vector3::new(1.5, 1.5, 2.0),
5.0,
Rc::new(LambertianMaterial::new_dummy()),
);
assert_matches!(s.intersect(&r), Some(_));
}
#[quickcheck]
fn ray_intersects_sphere_centre_at_correct_distance(
ray_origin: Vector3<f64>,
sphere_centre: Vector3<f64>,
radius: f64,
) -> TestResult {
if radius <= 0.0 || radius + 0.000001 >= (ray_origin - sphere_centre).norm() {
return TestResult::discard();
};
let sphere = Sphere::new(
sphere_centre,
radius,
Rc::new(LambertianMaterial::new_dummy()),
);
let ray = Ray::new(ray_origin, sphere_centre - ray_origin);
let info = sphere.intersect(&ray).unwrap();
let distance_to_centre = (sphere_centre - ray.origin).norm();
TestResult::from_bool(
(distance_to_centre - (info.distance + sphere.radius)).abs() < 0.00001,
)
}
#[test]
fn ray_intersects_plane() {
let r = Ray::new(Vector3::new(1.0, 2.0, 3.0), Vector3::new(-1.0, 0.0, 1.0));
let p = Plane::new(
Vector3::new(1.0, 0.0, 0.0),
-5.0,
Rc::new(LambertianMaterial::new_dummy()),
);
assert_matches!(p.intersect(&r), Some(_));
}
#[test]
fn ray_does_not_intersect_plane() {
let r = Ray::new(Vector3::new(1.0, 2.0, 3.0), Vector3::new(1.0, 0.0, 1.0));
let p = Plane::new(
Vector3::new(1.0, 0.0, 0.0),
-5.0,
Rc::new(LambertianMaterial::new_dummy()),
);
assert_matches!(p.intersect(&r), None);
}
#[test]
fn intersection_point_is_on_plane() {
let r = Ray::new(Vector3::new(1.0, 2.0, 3.0), Vector3::new(-1.0, 0.0, 1.0));
let p = Plane::new(
Vector3::new(1.0, 0.0, 0.0),
-5.0,
Rc::new(LambertianMaterial::new_dummy()),
);
match p.intersect(&r) {
Some(IntersectionInfo {
distance: _,
location,
normal: _,
tangent: _,
cotangent: _,
retro: _,
material: _,
}) => assert!((location.x - (-5.0f64)).abs() < 0.0000000001),
None => panic!(),
}
}
}

124
src/raycasting/axis_aligned_bounding_box.rs
View File

@ -1,124 +0,0 @@
use crate::util::Interval;
use super::{IntersectP, Ray};
use itertools::izip;
pub use crate::util::axis_aligned_bounding_box::BoundingBox;
impl IntersectP for BoundingBox {
fn intersect(&self, ray: &Ray) -> bool {
let mut t_interval_in_bounds = Interval::infinite();
for (&ray_origin, &ray_direction, bounds) in izip!(
ray.origin.coords.iter(),
ray.direction.coords.iter(),
self.bounds.iter()
) {
t_interval_in_bounds = t_interval_in_bounds.intersection(Interval::new(
(bounds.get_min() - ray_origin) / ray_direction,
(bounds.get_max() - ray_origin) / ray_direction,
));
if t_interval_in_bounds.is_empty() {
return false;
};
}
true
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::math::Vec3;
use quickcheck::TestResult;
use quickcheck_macros::quickcheck;
fn wrap_value_in_interval(value: f64, interval: Interval) -> f64 {
let distance_from_start = (value - interval.get_min()).abs();
let range = interval.get_max() - interval.get_min();
let multiple_of_range = distance_from_start / range;
return interval.get_min() + multiple_of_range.fract() * range;
}
#[quickcheck]
fn wrap_value_in_interval_produces_values_in_interval(v: f64, a: f64, b: f64) -> bool {
let interval = Interval::new(a, b);
interval.contains_value(wrap_value_in_interval(v, interval))
}
fn wrap_point_into_bounding_box(point: Vec3<f64>, bounds: &BoundingBox) -> Vec3<f64> {
let mut coords = [0.0; 3];
for i in 0..3 {
coords[i] = wrap_value_in_interval(point[i], bounds.bounds[i]);
}
Vec3 { coords }
}
#[quickcheck]
fn correctly_detects_intersections(
ray_origin: Vec3<f64>,
corner1: Vec3<f64>,
corner2: Vec3<f64>,
random_point: Vec3<f64>,
) -> bool {
let bounds = BoundingBox::from_corners(corner1, corner2);
let point_in_bounds = wrap_point_into_bounding_box(random_point, &bounds);
let ray = Ray::new(ray_origin, point_in_bounds - ray_origin);
bounds.intersect(&ray)
}
#[quickcheck]
fn intersect_always_true_when_ray_origin_is_inside_bounds(
ray_origin: Vec3<f64>,
corner1: Vec3<f64>,
corner2: Vec3<f64>,
random_point: Vec3<f64>,
) -> TestResult {
let bounds = BoundingBox::from_corners(corner1, corner2);
let ray_origin = wrap_point_into_bounding_box(ray_origin, &bounds);
let ray = Ray::new(ray_origin, ray_origin - random_point);
TestResult::from_bool(bounds.intersect(&ray))
}
#[quickcheck]
fn no_intersection_when_behind_ray(
ray_origin: Vec3<f64>,
corner1: Vec3<f64>,
corner2: Vec3<f64>,
random_point: Vec3<f64>,
) -> TestResult {
let bounds = BoundingBox::from_corners(corner1, corner2);
if bounds.contains_point(ray_origin) {
return TestResult::discard();
}
let point_in_bounds = wrap_point_into_bounding_box(random_point, &bounds);
let ray = Ray::new(ray_origin, ray_origin - point_in_bounds);
TestResult::from_bool(bounds.intersect(&ray))
}
#[test]
fn intersection_detected_when_ray_parallel_to_axis() {
let target =
BoundingBox::from_corners(Vec3::new(1.0f64, 2.0, 3.0), Vec3::new(4.0, 5.0, 6.0));
let x_ray = Ray::new(Vec3::new(0.0, 3.0, 4.0), Vec3::new(1.0, 0.0, 0.0));
assert!(target.intersect(&x_ray));
let y_ray = Ray::new(Vec3::new(2.0, 0.0, 4.0), Vec3::new(0.0, 1.0, 0.0));
assert!(target.intersect(&y_ray));
let z_ray = Ray::new(Vec3::new(2.0, 3.0, 0.0), Vec3::new(0.0, 0.0, 1.0));
assert!(target.intersect(&z_ray));
}
#[test]
fn intersection_missed_when_ray_parallel_to_axis() {
let target =
BoundingBox::from_corners(Vec3::new(1.0f64, 2.0, 3.0), Vec3::new(4.0, 5.0, 6.0));
let x_ray = Ray::new(Vec3::new(0.0, 0.0, 0.0), Vec3::new(1.0, 0.0, 0.0));
assert!(!target.intersect(&x_ray));
let y_ray = Ray::new(Vec3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 1.0, 0.0));
assert!(!target.intersect(&y_ray));
let z_ray = Ray::new(Vec3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0));
assert!(!target.intersect(&z_ray));
}
}

133
src/raycasting/bounding_volume_hierarchy.rs
View File

@ -1,133 +0,0 @@
use crate::math::{Affine3, Vec3};
use super::{BoundingBox, HasBoundingBox, Intersect, IntersectP, IntersectionInfo, Primitive, Ray};
use std::cmp::Ordering;
use std::sync::Arc;
/// Stores a set of [Primitives](Primitive) and accelerates raycasting
///
/// Organizes the primitives into a binary tree based on their bounds, allowing the
/// closest intersection with a ray to be found efficiently.
///
/// Each node knows the overall bounds of all it's children, which means that a ray that
/// doesn't intersect the [BoundingBox](BoundingBox) of the node doesn't intersect any of
/// the primitives stored in it's children.
pub enum BoundingVolumeHierarchy {
Node {
bounds: BoundingBox,
left: Box<BoundingVolumeHierarchy>,
right: Box<BoundingVolumeHierarchy>,
},
Leaf {
bounds: BoundingBox,
primitives: Vec<Arc<dyn Primitive>>,
},
}
fn centre(bounds: &BoundingBox) -> Vec3<f64> {
Vec3::new(
(bounds.bounds[0].get_min() + bounds.bounds[0].get_max()) / 2.00,
(bounds.bounds[1].get_min() + bounds.bounds[1].get_max()) / 2.0,
(bounds.bounds[2].get_min() + bounds.bounds[2].get_max()) / 2.0,
)
}
fn heuristic_split(primitives: &mut [Arc<dyn Primitive>], bounds: &BoundingBox) -> usize {
let largest_dimension = bounds.largest_dimension();
primitives.sort_unstable_by(|a, b| {
centre(&a.bounding_box())[largest_dimension]
.partial_cmp(&centre(&b.bounding_box())[largest_dimension])
.unwrap_or(Ordering::Equal)
});
primitives.len() / 2
}
impl BoundingVolumeHierarchy {
pub fn build(primitives: &mut [Arc<dyn Primitive>]) -> Self {
BoundingVolumeHierarchy::build_from_slice(primitives)
}
pub fn build_from_slice(primitives: &mut [Arc<dyn Primitive>]) -> Self {
let bounds = primitives
.iter()
.fold(BoundingBox::empty(), |acc, p| acc.union(&p.bounding_box()));
if primitives.len() <= 1 {
let primitives = primitives.to_vec();
BoundingVolumeHierarchy::Leaf { bounds, primitives }
} else {
let pivot = heuristic_split(primitives, &bounds);
let left = Box::new(BoundingVolumeHierarchy::build_from_slice(
&mut primitives[0..pivot],
));
let right = Box::new(BoundingVolumeHierarchy::build_from_slice(
&mut primitives[pivot..],
));
BoundingVolumeHierarchy::Node {
bounds,
left,
right,
}
}
}
}
fn closest_intersection(
a: Option<IntersectionInfo>,
b: Option<IntersectionInfo>,
) -> Option<IntersectionInfo> {
match a {
None => b,
Some(a_info) => match b {
None => Some(a_info),
Some(b_info) => Some(if a_info.distance < b_info.distance {
a_info
} else {
b_info
}),
},
}
}
impl Intersect for BoundingVolumeHierarchy {
fn intersect(&self, ray: &Ray) -> Option<IntersectionInfo> {
match self {
BoundingVolumeHierarchy::Node {
bounds,
left,
right,
} => {
if bounds.intersect(ray) {
closest_intersection(left.intersect(ray), right.intersect(ray))
} else {
None
}
}
BoundingVolumeHierarchy::Leaf { bounds, primitives } => {
if bounds.intersect(ray) {
primitives
.iter()
.map(|elem| elem.intersect(ray))
.fold(None, closest_intersection)
} else {
None
}
}
}
}
}
impl HasBoundingBox for BoundingVolumeHierarchy {
fn bounding_box(&self) -> BoundingBox {
BoundingBox::empty()
}
}
impl Primitive for BoundingVolumeHierarchy {
fn transform(&self, transformation: &Affine3<f64>) -> Box<dyn Primitive> {
todo!()
}
}
#[cfg(test)]
mod test {}

171
src/raycasting/mod.rs
View File

@ -1,171 +0,0 @@
use crate::math::{Affine3,Vec3};
use super::materials::Material;
use std::sync::Arc;
pub mod sphere;
pub use sphere::Sphere;
pub mod plane;
pub use plane::Plane;
pub mod triangle;
pub use triangle::Triangle;
pub mod axis_aligned_bounding_box;
pub use axis_aligned_bounding_box::BoundingBox;
pub mod bounding_volume_hierarchy;
pub use bounding_volume_hierarchy::BoundingVolumeHierarchy;
pub mod vec_aggregate;
/// A ray, consisting or a start point and direction
///
/// This is the basic ray struct used to define things like a line-of-sight
/// going out from the camera of a reflection from a surface.
#[derive(Clone, Debug)]
pub struct Ray {
/// The start point of the ray
pub origin: Vec3<f64>,
/// The direction the ray goes in.
///
/// This vector should always be kept normalized
pub direction: Vec3<f64>,
}
impl Ray {
/// Create a new ray
pub fn new(origin: Vec3<f64>, direction: Vec3<f64>) -> Ray {
Ray {
origin,
direction: direction.normalize(),
}
}
/// Return the point on the ray that is `t` units from the start
pub fn point_at(&self, t: f64) -> Vec3<f64> {
self.origin + self.direction * t
}
/// Create a new ray by moving the original ray along it's direction by `amount`
///
/// `amount` is normally a very small number. This function is useful for ensuring
/// that rounding-errors don;t cause a reflection ray doesn't intersect with the point
/// it's reflected from.
pub fn bias(&self, amount: f64) -> Ray {
Ray::new(self.origin + self.direction * amount, self.direction)
}
}
/// Information about a ray-primitive intersection.
///
/// This struct is returned by [intersect()](Intersect::intersect) and contatins all the
/// information needed to evaluate the rendering function for that intersection.
#[derive(Debug)]
pub struct IntersectionInfo {
/// The distance between the ray origin and the intersection point
pub distance: f64,
/// The intersection point
pub location: Vec3<f64>,
/// The surface normal at the intersection point
pub normal: Vec3<f64>,
/// The surface tangent at the intersection point
///
/// Which surface tangent direction returned is dependent on the [Primitive](Primitive)
/// but should generally be smooth over any given surface
pub tangent: Vec3<f64>,
/// Another surface tangent, perpendicular to `tangent`
///
/// The cross product or `normal` and `tangent`
pub cotangent: Vec3<f64>,
/// The direction from the intersection point back towards the ray
///
/// Equal to `-ray.direction`
pub retro: Vec3<f64>,
/// The [Material](crate::materials::Material) which describes the optical
/// properties of the intersected surface
pub material: Arc<dyn Material>,
}
/// A geometric object that has a [Material](crate::materials::Material) and can be
/// intersected with a [Ray](Ray)
pub trait Intersect: Send + Sync {
/// Test if the ray intersects the object, and return information about the object and intersection.
fn intersect(&self, ray: &Ray) -> Option<IntersectionInfo>;
}
/// A geometric object that can be intersected with a ray
///
/// This is useful for objects that don't have materials (such as [BoundingBox](BoundingBox))
/// and as a (possibly) faster alternative to [Intersect](Intersect) when only a simple
/// intersection test is needed.
pub trait IntersectP: Send + Sync {
/// Test if the ray intersects the object, without calculating any extra information.
fn intersect(&self, ray: &Ray) -> bool;
}
/// Any geometric object for which a [BoundingBox](BoundingBox) can be calculated
pub trait HasBoundingBox: Send + Sync {
/// The axis-aligned bounding box of the object
///
/// The object must fit entirely inside this box.
fn bounding_box(&self) -> BoundingBox;
}
/// A basic geometric primitive such as a sphere or a triangle
pub trait Primitive: Intersect + HasBoundingBox {
// / Create a new object by applying the transformation to this object.
fn transform(&self, transformation: &Affine3<f64>) -> Box<dyn Primitive>;
}
#[cfg(test)]
mod tests {
use quickcheck_macros::quickcheck;
use super::*;
use quickcheck::{Arbitrary, Gen};
impl Arbitrary for Ray {
fn arbitrary<G: Gen>(g: &mut G) -> Ray {
let origin = <Vec3<f64> as Arbitrary>::arbitrary(g);
let direction = <Vec3<f64> as Arbitrary>::arbitrary(g);
return Ray::new(origin, direction);
}
}
#[quickcheck]
fn t0_is_origin(ray: Ray) -> bool {
ray.point_at(0.0) == ray.origin
}
#[quickcheck]
fn t1_is_origin_plus_direction(ray: Ray) -> bool {
ray.point_at(1.0) == ray.origin + ray.direction
}
#[quickcheck]
fn points_are_colinear(ray: Ray, t1: f64, t2: f64, t3: f64) -> bool {
let p1 = ray.point_at(t1);
let p2 = ray.point_at(t2);
let p3 = ray.point_at(t3);
let epsilon = [t1, t2, t3, ray.origin[0], ray.origin[1], ray.origin[2]]
.iter()
.fold(0.0f64, |a, &b| a.max(b.abs()))
* std::f64::EPSILON
* 256.0f64;
(p2 - p1).cross(&(p3 - p2)).norm() < epsilon
}
#[quickcheck]
fn t_is_distance(ray: Ray, t: f64) -> bool {
(ray.point_at(t) - ray.origin).norm() - t.abs() < 0.0000000001
}
}

290
src/raycasting/plane.rs
View File

@ -1,290 +0,0 @@
use crate::materials::Material;
use crate::math::Vec3;
use super::{BoundingBox, HasBoundingBox, Intersect, IntersectionInfo, Primitive, Ray};
use std::sync::Arc;
#[derive(Clone)]
pub struct Plane {
normal: Vec3<f64>,
tangent: Vec3<f64>,
cotangent: Vec3<f64>,
distance_from_origin: f64,
material: Arc<dyn Material>,
}
impl Plane {
pub fn new(normal: Vec3<f64>, distance_from_origin: f64, material: Arc<dyn Material>) -> Plane {
let normal = normal.normalize();
let mut axis_closest_to_tangent = Vec3::zeros();
axis_closest_to_tangent[normal.smallest_coord()] = 1.0;
let cotangent = normal.cross(&axis_closest_to_tangent).normalize();
let tangent = normal.cross(&cotangent);
Plane {
normal,
tangent,
cotangent,
distance_from_origin,
material,
}
}
}
/*impl Transform for Plane {
fn transform(&self, transformation: &Affine3<f64>) -> Self {
Plane {
normal: transformation.transform_vector(&self.normal).normalize(),
cotangent: transformation.transform_vector(&self.cotangent).normalize(),
tangent: self.normal.cross(&self.cotangent),
distance_from_origin: transformation
.transform_vector(&(self.normal * self.distance_from_origin))
.norm(),
material: Arc::clone(&self.material),
}
}
}*/
impl Intersect for Plane {
fn intersect(&self, ray: &Ray) -> Option<IntersectionInfo> {
let ray_direction_dot_plane_normal = ray.direction.dot(&self.normal);
let point_on_plane = self.normal * self.distance_from_origin;
let point_on_plane_minus_ray_origin_dot_normal =
(point_on_plane - ray.origin).dot(&self.normal);
if ray_direction_dot_plane_normal == 0.0 {
//Ray is parallel to plane
if point_on_plane_minus_ray_origin_dot_normal != 0.0 {
//Ray is not in plane
return None;
}
}
let t = point_on_plane_minus_ray_origin_dot_normal / ray_direction_dot_plane_normal;
if t < 0.0 {
return None;
}
Some(IntersectionInfo {
distance: t,
location: ray.point_at(t),
normal: self.normal,
tangent: self.tangent,
cotangent: self.cotangent,
retro: -ray.direction,
material: Arc::clone(&self.material),
})
}
}
impl HasBoundingBox for Plane {
fn bounding_box(&self) -> BoundingBox {
let p0 = self.normal * self.distance_from_origin;
let f = |v: Vec3<f64>| {
Vec3::new(
if v.x() == 0.0 { 0.0 } else { f64::INFINITY },
if v.y() == 0.0 { 0.0 } else { f64::INFINITY },
if v.z() == 0.0 { 0.0 } else { f64::INFINITY },
)
};
let tangent = f(self.tangent);
let cotangent = f(self.cotangent);
let p1 = p0 + tangent;
let p2 = p0 - tangent;
let p3 = p0 + cotangent;
let p4 = p0 - cotangent;
BoundingBox::from_points(&[p1, p2, p3, p4])
}
}
impl Primitive for Plane {
fn transform(&self, transformation: &crate::math::Affine3<f64>) -> Box<dyn Primitive> {
todo!()
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::materials::LambertianMaterial;
use crate::math::Vec3;
#[test]
fn ray_intersects_plane() {
let r = Ray::new(Vec3::new(1.0, 2.0, 3.0), Vec3::new(-1.0, 0.0, 1.0));
let p = Plane::new(
Vec3::new(1.0, 0.0, 0.0),
-5.0,
Arc::new(LambertianMaterial::new_dummy()),
);
if let None = p.intersect(&r) {
panic!("Intersection failed.");
}
}
#[test]
fn ray_does_not_intersect_plane() {
let r = Ray::new(Vec3::new(1.0, 2.0, 3.0), Vec3::new(1.0, 0.0, 1.0));
let p = Plane::new(
Vec3::new(1.0, 0.0, 0.0),
-5.0,
Arc::new(LambertianMaterial::new_dummy()),
);
if let Some(_) = p.intersect(&r) {
panic!("Intersection failed.");
}
}
#[test]
fn intersection_point_is_on_plane() {
let r = Ray::new(Vec3::new(1.0, 2.0, 3.0), Vec3::new(-1.0, 0.0, 1.0));
let p = Plane::new(
Vec3::new(1.0, 0.0, 0.0),
-5.0,
Arc::new(LambertianMaterial::new_dummy()),
);
match p.intersect(&r) {
Some(IntersectionInfo {
distance: _,
location,
normal: _,
tangent: _,
cotangent: _,
retro: _,
material: _,
}) => assert!((location.x() - (-5.0f64)).abs() < 0.0000000001),
None => panic!(),
}
}
#[test]
fn bounding_box_is_correct_for_yz_plane() {
let target = Plane::new(
Vec3::new(1.0, 0.0, 0.0),
2.0,
Arc::new(LambertianMaterial::new_dummy()),
);
let bb = target.bounding_box();
assert!(!bb.contains_point(Vec3::new(1.0, 2.0, 3.0)));
assert!(bb.contains_point(Vec3::new(2.0, 2.0, 3.0)));
assert!(bb.contains_point(Vec3::new(2.0, 2000.0, 3.0)));
assert!(bb.contains_point(Vec3::new(2.0, 0.0, 3.0)));
assert!(bb.contains_point(Vec3::new(2.0, -2000.0, 3.0)));
assert!(bb.contains_point(Vec3::new(2.0, 2.0, 3000.0)));
assert!(bb.contains_point(Vec3::new(2.0, 2.0, 0.0)));
assert!(bb.contains_point(Vec3::new(2.0, 2.0, -3000.0)));
assert!(!bb.contains_point(Vec3::new(3.0, 2.0, 3.0)));
}
#[test]
fn bounding_box_is_correct_for_yz_plane_with_negative_normal() {
let target = Plane::new(
Vec3::new(-1.0, 0.0, 0.0),
2.0,
Arc::new(LambertianMaterial::new_dummy()),
);
let bb = target.bounding_box();
assert!(!bb.contains_point(Vec3::new(1.0, 2.0, 3.0)));
assert!(bb.contains_point(Vec3::new(-2.0, 2.0, 3.0)));
assert!(bb.contains_point(Vec3::new(-2.0, 2000.0, 3.0)));
assert!(bb.contains_point(Vec3::new(-2.0, 0.0, 3.0)));
assert!(bb.contains_point(Vec3::new(-2.0, -2000.0, 3.0)));
assert!(bb.contains_point(Vec3::new(-2.0, 2.0, 3000.0)));
assert!(bb.contains_point(Vec3::new(-2.0, 2.0, 0.0)));
assert!(bb.contains_point(Vec3::new(-2.0, 2.0, -3000.0)));
assert!(!bb.contains_point(Vec3::new(-3.0, 2.0, 3.0)));
}
#[test]
fn bounding_box_is_correct_for_xz_plane() {
let target = Plane::new(
Vec3::new(0.0, 1.0, 0.0),
2.0,
Arc::new(LambertianMaterial::new_dummy()),
);
let bb = target.bounding_box();
assert!(!bb.contains_point(Vec3::new(1.0, 1.0, 3.0)));
assert!(bb.contains_point(Vec3::new(1.0, 2.0, 3.0)));
assert!(bb.contains_point(Vec3::new(1000.0, 2.0, 3.0)));
assert!(bb.contains_point(Vec3::new(0.0, 2.0, 3.0)));
assert!(bb.contains_point(Vec3::new(-1000.0, 2.0, 3.0)));
assert!(bb.contains_point(Vec3::new(1.0, 2.0, 3000.0)));
assert!(bb.contains_point(Vec3::new(1.0, 2.0, 0.0)));
assert!(bb.contains_point(Vec3::new(1.0, 2.0, -3000.0)));
assert!(!bb.contains_point(Vec3::new(1.0, 3.0, 3.0)));
}
#[test]
fn bounding_box_is_correct_for_xz_plane_with_negative_normal() {
let target = Plane::new(
Vec3::new(0.0, -1.0, 0.0),
2.0,
Arc::new(LambertianMaterial::new_dummy()),
);
let bb = target.bounding_box();
assert!(!bb.contains_point(Vec3::new(1.0, -1.0, 3.0)));
assert!(bb.contains_point(Vec3::new(1.0, -2.0, 3.0)));
assert!(bb.contains_point(Vec3::new(1000.0, -2.0, 3.0)));
assert!(bb.contains_point(Vec3::new(0.0, -2.0, 3.0)));
assert!(bb.contains_point(Vec3::new(-1000.0, -2.0, 3.0)));
assert!(bb.contains_point(Vec3::new(1.0, -2.0, 3000.0)));
assert!(bb.contains_point(Vec3::new(1.0, -2.0, 0.0)));
assert!(bb.contains_point(Vec3::new(1.0, -2.0, -3000.0)));
assert!(!bb.contains_point(Vec3::new(1.0, 3.0, 3.0)));
}
#[test]
fn bounding_box_is_correct_for_xy_plane() {
let target = Plane::new(
Vec3::new(0.0, 0.0, 1.0),
2.0,
Arc::new(LambertianMaterial::new_dummy()),
);
let bb = target.bounding_box();
assert!(!bb.contains_point(Vec3::new(1.0, 2.0, 1.0)));
assert!(bb.contains_point(Vec3::new(1.0, 2.0, 2.0)));
assert!(bb.contains_point(Vec3::new(1.0, 2000.0, 2.0)));
assert!(bb.contains_point(Vec3::new(1.0, 0.0, 2.0)));
assert!(bb.contains_point(Vec3::new(1.0, -2000.0, 2.0)));
assert!(bb.contains_point(Vec3::new(2000.0, 2.0, 2.0)));
assert!(bb.contains_point(Vec3::new(0.0, 2.0, 2.0)));
assert!(bb.contains_point(Vec3::new(-2000.0, 2.0, 2.0)));
assert!(!bb.contains_point(Vec3::new(3.0, 2.0, 3.0)));
}
#[test]
fn bounding_box_is_correct_for_xy_plane_with_negative_normal() {
let target = Plane::new(
Vec3::new(0.0, 0.0, -1.0),
-2.0,
Arc::new(LambertianMaterial::new_dummy()),
);
let bb = target.bounding_box();
assert!(!bb.contains_point(Vec3::new(1.0, 2.0, 1.0)));
assert!(bb.contains_point(Vec3::new(1.0, 2.0, 2.0)));
assert!(bb.contains_point(Vec3::new(1.0, 2000.0, 2.0)));
assert!(bb.contains_point(Vec3::new(1.0, 0.0, 2.0)));
assert!(bb.contains_point(Vec3::new(1.0, -2000.0, 2.0)));
assert!(bb.contains_point(Vec3::new(2000.0, 2.0, 2.0)));
assert!(bb.contains_point(Vec3::new(0.0, 2.0, 2.0)));
assert!(bb.contains_point(Vec3::new(-2000.0, 2.0, 2.0)));
assert!(!bb.contains_point(Vec3::new(3.0, 2.0, 3.0)));
}
#[test]
fn bounding_box_is_infinite_when_normal_is_not_aligned_with_axis() {
let target = Plane::new(
Vec3::new(0.1, 0.0, -1.0),
-2.0,
Arc::new(LambertianMaterial::new_dummy()),
);
let bb = target.bounding_box();
assert!(bb.contains_point(Vec3::new(1.0, 2.0, 1.0)));
assert!(bb.contains_point(Vec3::new(1.0, 2.0, 2.0)));
assert!(bb.contains_point(Vec3::new(1.0, 2000.0, 2.0)));
assert!(bb.contains_point(Vec3::new(1.0, 0.0, 2.0)));
assert!(bb.contains_point(Vec3::new(1.0, -2000.0, 2.0)));
assert!(bb.contains_point(Vec3::new(2000.0, 2.0, 2.0)));
assert!(bb.contains_point(Vec3::new(0.0, 2.0, 2.0)));
assert!(bb.contains_point(Vec3::new(-2000.0, 2.0, 2.0)));
assert!(bb.contains_point(Vec3::new(3.0, 2.0, 3.0)));
}
}

265
src/raycasting/sphere.rs
View File

@ -1,265 +0,0 @@
use crate::materials::Material;
use crate::math::Vec3;
use super::{BoundingBox, HasBoundingBox, Intersect, IntersectionInfo, Primitive, Ray};
use std::sync::Arc;
#[derive(Clone, Debug)]
pub struct Sphere {
centre: Vec3<f64>,
radius: f64,
material: Arc<dyn Material>,
}
impl Sphere {
pub fn new(centre: Vec3<f64>, radius: f64, material: Arc<dyn Material>) -> Sphere {
Sphere {
centre,
radius,
material,
}
}
}
/*impl Transform for Sphere {
fn transform(&self, transformation: &Affine3<f64>) -> Self {
Sphere {
centre: transformation.transform_point(&self.centre),
// This is not the most efficient way of calculating the radius,
//but will work as long as the resulting shape is still a sphere.
radius: transformation
.transform_vector(&Vector3::new(self.radius, 0.0, 0.0))
.norm(),
material: Arc::clone(&self.material),
}
}
}*/
impl Intersect for Sphere {
fn intersect<'a>(&'_ self, ray: &Ray) -> Option<IntersectionInfo> {
let r_o = ray.origin;
let centre_coords = self.centre;
let a = ray
.direction
.component_mul(&ray.direction)
.coords
.iter()
.fold(0.0, |a, b| a + *b);
let b = ((r_o.component_mul(&ray.direction) - centre_coords.component_mul(&ray.direction))
* 2.0)
.coords
.iter()
.fold(0.0, |a, b| a + *b);
let c = (r_o.component_mul(&r_o) + centre_coords.component_mul(&centre_coords)
- centre_coords.component_mul(&r_o) * 2.0)
.coords
.iter()
.fold(0.0, |a, b| a + *b)
- self.radius * self.radius;
let delta_squared = b * b - 4.0 * a * c;
if delta_squared < 0.0 {
None
} else {
let delta = delta_squared.sqrt();
let one_over_2_a = 1.0 / (2.0 * a);
let t1 = (-b - delta) * one_over_2_a;
let t2 = (-b + delta) * one_over_2_a;
let distance = if t1 < 0.0 || (t2 >= 0.0 && t1 >= t2) {
t2
} else {
t1
};
if distance <= 0.0 {
None
} else {
let location = ray.point_at(distance);
let normal = (location - self.centre).normalize();
let tangent = normal.cross(&Vec3::unit_z()).normalize();
let cotangent = normal.cross(&tangent);
let retro = -ray.direction;
Some(IntersectionInfo {
distance,
location,
normal,
tangent,
cotangent,
retro,
material: Arc::clone(&self.material),
})
}
}
}
}
impl HasBoundingBox for Sphere {
fn bounding_box(&self) -> BoundingBox {
let radius_xyz = Vec3::new(self.radius, self.radius, self.radius);
BoundingBox::from_corners(self.centre + radius_xyz, self.centre - radius_xyz)
}
}
impl Primitive for Sphere {
fn transform(&self, transformation: &crate::math::Affine3<f64>) -> Box<dyn Primitive> {
todo!()
}
}
#[cfg(test)]
mod tests {
use quickcheck::TestResult;
use quickcheck_macros::quickcheck;
use super::*;
use crate::materials::LambertianMaterial;
#[test]
fn ray_intersects_sphere() {
let r = Ray::new(Vec3::new(1.0, 2.0, 3.0), Vec3::new(0.0, 0.0, 1.0));
let s = Sphere::new(
Vec3::new(1.5, 1.5, 15.0),
5.0,
Arc::new(LambertianMaterial::new_dummy()),
);
if let None = s.intersect(&r) {
panic!("Intersection failed");
}
}
#[test]
fn ray_does_not_intersect_sphere_when_sphere_is_in_front() {
let r = Ray::new(Vec3::new(1.0, 2.0, 3.0), Vec3::new(0.0, 0.0, 1.0));
let s = Sphere::new(
Vec3::new(-5.0, 1.5, 15.0),
5.0,
Arc::new(LambertianMaterial::new_dummy()),
);
if let Some(_) = s.intersect(&r) {
panic!("Intersection passed.");
}
}
#[test]
fn ray_does_not_intersect_sphere_when_sphere_is_behind() {
let r = Ray::new(Vec3::new(1.0, 2.0, 3.0), Vec3::new(0.0, 0.0, 1.0));
let s = Sphere::new(
Vec3::new(1.5, 1.5, -15.0),
5.0,
Arc::new(LambertianMaterial::new_dummy()),
);
if let Some(_) = s.intersect(&r) {
panic!("Intersection failed");
}
}
#[test]
fn ray_intersects_sphere_when_origin_is_inside() {
let r = Ray::new(Vec3::new(1.0, 2.0, 3.0), Vec3::new(0.0, 0.0, 1.0));
let s = Sphere::new(
Vec3::new(1.5, 1.5, 2.0),
5.0,
Arc::new(LambertianMaterial::new_dummy()),
);
if let None = s.intersect(&r) {
panic!("Intersection failed");
}
}
#[quickcheck]
fn ray_intersects_sphere_centre_at_correct_distance(
ray_origin: Vec3<f64>,
sphere_centre: Vec3<f64>,
radius: f64,
) -> TestResult {
if radius <= 0.0 || radius + 0.000001 >= (ray_origin - sphere_centre).norm() {
return TestResult::discard();
};
let sphere = Sphere::new(
sphere_centre,
radius,
Arc::new(LambertianMaterial::new_dummy()),
);
let ray = Ray::new(ray_origin, sphere_centre - ray_origin);
let info = sphere.intersect(&ray).unwrap();
let distance_to_centre = (sphere_centre - ray.origin).norm();
TestResult::from_bool(
(distance_to_centre - (info.distance + sphere.radius)).abs() < 0.00001,
)
}
#[quickcheck]
fn all_points_on_sphere_are_in_bounding_box(sphere_centre: Vec3<f64>, radius_vector: Vec3<f64>) -> bool {
let target_sphere = Sphere::new(
sphere_centre,
radius_vector.norm(),
Arc::new(LambertianMaterial::new_dummy()),
);
let bounding_box = target_sphere.bounding_box();
bounding_box.contains_point(sphere_centre + radius_vector)
}
/*#[quickcheck]
fn translation_moves_centre(
sphere_centre: Vec3<f64>,
radius: f64,
translation_vector: Vec3<f64>,
) -> TestResult {
if radius <= 0.0 {
return TestResult::discard();
};
let sphere = Sphere::new(
sphere_centre,
radius,
Arc::new(LambertianMaterial::new_dummy()),
);
let expected_centre = sphere.centre + translation_vector;
let mut transformation = Affine3::identity();
transformation *= Translation3::from(translation_vector);
let sphere = sphere.transform(&transformation);
TestResult::from_bool(expected_centre == sphere.centre)
}
#[quickcheck]
fn translation_does_not_change_radius(
sphere_centre: Vec3<f64>,
radius: f64,
translation_vector: Vec3<f64>,
) -> TestResult {
if radius <= 0.0 {
return TestResult::discard();
};
let sphere = Sphere::new(
sphere_centre,
radius,
Arc::new(LambertianMaterial::new_dummy()),
);
let expected_radius = sphere.radius;
let mut transformation = Affine3::identity();
transformation *= Translation3::from(translation_vector);
let sphere = sphere.transform(&transformation);
TestResult::from_bool(expected_radius == sphere.radius)
}
#[quickcheck]
fn rotation_about_centre_does_not_move_centre(
sphere_centre: Vec3<f64>,
radius: f64,
rotation_vector: Vec3<f64>,
) -> TestResult {
if radius <= 0.0 {
return TestResult::discard();
};
let sphere = Sphere::new(
sphere_centre,
radius,
Arc::new(LambertianMaterial::new_dummy()),
);
let expected_centre = sphere.centre;
let mut transformation = Affine3::identity();
transformation *= Translation3::from(sphere.centre.coords)
* Rotation3::new(rotation_vector)
* Translation3::from(-sphere.centre.coords);
let sphere = sphere.transform(&transformation);
TestResult::from_bool((expected_centre - sphere.centre).norm() < 0.000000001)
}*/
}

922
src/raycasting/triangle.rs
View File

@ -1,922 +0,0 @@
use crate::materials::Material;
use crate::math::{Vec2, Vec3};
use super::{BoundingBox, HasBoundingBox, Intersect, IntersectionInfo, Primitive, Ray};
use std::sync::Arc;
#[derive(Debug, Clone)]
pub struct Triangle {
pub vertices: [Vec3<f64>; 3],
pub normals: [Vec3<f64>; 3],
pub material: Arc<dyn Material>,
}
/*impl Transform for Triangle {
fn transform(&self, transformation: &Affine3<f64>) -> Self {
let normal_transformation =
Affine3::from_matrix_unchecked(transformation.inverse().matrix().transpose());
Triangle {
vertices: [
transformation.transform_point(&self.vertices[0]),
transformation.transform_point(&self.vertices[1]),
transformation.transform_point(&self.vertices[2]),
],
normals: [
normal_transformation.transform_vector(&self.normals[0]),
normal_transformation.transform_vector(&self.normals[1]),
normal_transformation.transform_vector(&self.normals[2]),
],
material: Arc::clone(&self.material),
}
}
}*/
impl Intersect for Triangle {
fn intersect(&self, ray: &Ray) -> Option<IntersectionInfo> {
let translation = -ray.origin;
let indices = indices_with_index_of_largest_element_last(&ray.direction);
let permuted_ray_direction = permute_vector_elements(&ray.direction, &indices);
let shear_slopes = calculate_shear_to_z_axis(&permuted_ray_direction);
let transformed_vertices: Vec<Vec3<f64>> = self
.vertices
.iter()
.map(|elem| {
apply_shear_to_z_axis(
&permute_vector_elements(&(elem + translation), &indices),
&shear_slopes,
)
})
.collect();
let edge_functions = signed_edge_functions(&transformed_vertices);
if edge_functions.coords.iter().all(|e| e.is_sign_positive())
|| edge_functions.coords.iter().all(|e| e.is_sign_negative())
{
let barycentric_coordinates =
barycentric_coordinates_from_signed_edge_functions(edge_functions.abs());
let transformed_z = barycentric_coordinates
.coords
.iter()
.zip(transformed_vertices.iter())
.map(|(&coord, vertex)| vertex.z() * coord)
.fold(0.0, |acc, z| acc + z);
if transformed_z.is_sign_positive() != permuted_ray_direction.z().is_sign_positive() {
return None;
}
let location = barycentric_coordinates
.coords
.iter()
.zip(self.vertices.iter())
.map(|(&barycentric_coord, vertex)| vertex * barycentric_coord)
.fold(Vec3::zeros(), |a, e| a + e);
let distance = (ray.origin - location).norm();
let normal: Vec3<f64> = barycentric_coordinates
.coords
.iter()
.zip(self.normals.iter())
.fold(Vec3::zeros(), |acc, (&coord, vertex)| acc + vertex * coord)
.normalize();
let cotangent = (self.vertices[0] - self.vertices[1])
.cross(&normal)
.normalize();
let tangent = cotangent.cross(&normal).normalize();
let retro = (ray.origin - location).normalize();
let material = Arc::clone(&self.material);
Some(IntersectionInfo {
distance,
location,
normal,
tangent,
cotangent,
retro,
material,
})
} else {
None
}
}
}
impl HasBoundingBox for Triangle {
fn bounding_box(&self) -> BoundingBox {
BoundingBox::from_points(&self.vertices)
}
}
impl Primitive for Triangle {
fn transform(&self, transformation: &crate::math::Affine3<f64>) -> Box<dyn Primitive> {
todo!()
}
}
fn indices_with_index_of_largest_element_last(v: &Vec3<f64>) -> [usize; 3] {
#[allow(clippy::collapsible_else_if)]
if v.x() > v.y() {
if v.z() > v.x() {
[0, 1, 2]
} else {
[1, 2, 0]
}
} else {
if v.z() > v.y() {
[0, 1, 2]
} else {
[2, 0, 1]
}
}
}
fn is_valid_permutation(indices: &[usize; 3]) -> bool {
(0..2).all(|i: usize| indices.contains(&i))
}
fn permute_vector_elements(v: &Vec3<f64>, indices: &[usize; 3]) -> Vec3<f64> {
debug_assert!(is_valid_permutation(indices));
Vec3::new(v[indices[0]], v[indices[1]], v[indices[2]])
}
fn calculate_shear_to_z_axis(v: &Vec3<f64>) -> Vec2<f64> {
Vec2::new(-v.x() / v.z(), -v.y() / v.z())
}
fn apply_shear_to_z_axis(v: &Vec3<f64>, s: &Vec2<f64>) -> Vec3<f64> {
Vec3::new(v.x() + s.x() * v.z(), v.y() + s.y() * v.z(), v.z())
}
fn signed_edge_function(a: &Vec3<f64>, b: &Vec3<f64>) -> f64 {
a.x() * b.y() - b.x() * a.y()
}
fn signed_edge_functions(vertices: &[Vec3<f64>]) -> Vec3<f64> {
// Iterate over the inputs in such a way that each output element is calculated
// from the twoother elements of the input. ( (y,z) -> x, (z,x) -> y, (x,y) -> z )
let coords: Vec<_> = vertices
.iter()
.cycle()
.skip(1)
.zip(vertices.iter().cycle().skip(2))
.take(vertices.len())
.map(|(v1, v2)| signed_edge_function(v1, v2))
.collect();
Vec3::new(coords[0], coords[1], coords[2])
}
fn barycentric_coordinates_from_signed_edge_functions(e: Vec3<f64>) -> Vec3<f64> {
e * (1.0 / e.coords.iter().fold(0.0, |a, &b| a + b))
}
#[cfg(test)]
mod tests {
use super::*;
/*mod triangle_transform {
use super::*;
use quickcheck_macros::quickcheck;
use crate::materials::LambertianMaterial;
#[quickcheck]
fn transform_by_identity_does_not_change_values(
v0: Vec3,
v1: Vec3,
v2: Vec3,
n0: Vec3,
n1: Vec3,
n2: Vec3,
) -> bool {
let n0 = n0.normalize();
let n1 = n1.normalize();
let n2 = n2.normalize();
let target = Triangle {
vertices: [v0, v1, v2],
normals: [n0, n1, n2],
material: Arc::new(LambertianMaterial::new_dummy()),
};
let target = target.transform(&Affine3::identity());
target.vertices[0] == v0
&& target.vertices[1] == v1
&& target.vertices[2] == v2
&& target.normals[0] == n0
&& target.normals[1] == n1
&& target.normals[2] == n2
}
#[quickcheck]
fn translate_does_not_change_normals(
v0: Vec3,
v1: Vec3,
v2: Vec3,
n0: Vec3,
n1: Vec3,
n2: Vec3,
translation: Vec3,
) -> bool {
let n0 = n0.normalize();
let n1 = n1.normalize();
let n2 = n2.normalize();
let target = Triangle {
vertices: [v0, v1, v2],
normals: [n0, n1, n2],
material: Arc::new(LambertianMaterial::new_dummy()),
};
let transformation = Affine3::identity() * Translation3::from(translation);
let target = target.transform(&transformation);
target.normals[0] == n0 && target.normals[1] == n1 && target.normals[2] == n2
}
#[quickcheck]
fn translate_translates_vertices(
v0: Vec3,
v1: Vec3,
v2: Vec3,
n0: Vec3,
n1: Vec3,
n2: Vec3,
translation: Vec3,
) -> bool {
let n0 = n0.normalize();
let n1 = n1.normalize();
let n2 = n2.normalize();
let target = Triangle {
vertices: [v0, v1, v2],
normals: [n0, n1, n2],
material: Arc::new(LambertianMaterial::new_dummy()),
};
let transformation = Affine3::identity() * Translation3::from(translation);
let target = target.transform(&transformation);
target.vertices[0] == v0 + translation
&& target.vertices[1] == v1 + translation
&& target.vertices[2] == v2 + translation
}
}*/
mod index_of_largest_element {
use super::*;
use quickcheck_macros::quickcheck;
#[quickcheck]
fn result_is_valid_permutation(v: Vec3<f64>) -> bool {
let indices = indices_with_index_of_largest_element_last(&v);
is_valid_permutation(&indices)
}
#[quickcheck]
fn result_includes_x(v: Vec3<f64>) -> bool {
let indices = indices_with_index_of_largest_element_last(&v);
indices.iter().any(|&i| i == 0)
}
#[quickcheck]
fn result_includes_y(v: Vec3<f64>) -> bool {
let indices = indices_with_index_of_largest_element_last(&v);
indices.iter().any(|&i| i == 1)
}
#[quickcheck]
fn result_includes_z(v: Vec3<f64>) -> bool {
let indices = indices_with_index_of_largest_element_last(&v);
indices.iter().any(|&i| i == 2)
}
#[quickcheck]
fn last_index_is_greater_than_or_equal_to_x(v: Vec3<f64>) -> bool {
let indices = indices_with_index_of_largest_element_last(&v);
v[indices[2]] >= v.x()
}
#[quickcheck]
fn last_index_is_greater_than_or_equal_to_y(v: Vec3<f64>) -> bool {
let indices = indices_with_index_of_largest_element_last(&v);
v[indices[2]] >= v.y()
}
#[quickcheck]
fn last_index_is_greater_than_or_equal_to_z(v: Vec3<f64>) -> bool {
let indices = indices_with_index_of_largest_element_last(&v);
v[indices[2]] >= v.z()
}
}
mod permute_vector_elements {
use super::*;
use quickcheck_macros::quickcheck;
#[quickcheck]
fn last_index_is_greater_than_or_equal_to_x(v: Vec3<f64>) -> bool {
let p = permute_vector_elements(&v, &indices_with_index_of_largest_element_last(&v));
p.z() >= v.x()
}
#[quickcheck]
fn last_index_is_greater_than_or_equal_to_y(v: Vec3<f64>) -> bool {
let p = permute_vector_elements(&v, &indices_with_index_of_largest_element_last(&v));
p.z() >= v.y()
}
#[quickcheck]
fn last_index_is_greater_than_or_equal_to_z(v: Vec3<f64>) -> bool {
let p = permute_vector_elements(&v, &indices_with_index_of_largest_element_last(&v));
p.z() >= v.z()
}
}
mod shear_to_z_axis {
use super::*;
use quickcheck_macros::quickcheck;
#[quickcheck]
fn shear_to_z_axis_makes_x_zero(v: Vec3<f64>) -> bool {
let s = calculate_shear_to_z_axis(&v);
apply_shear_to_z_axis(&v, &s).x().abs() < 0.00001
}
#[quickcheck]
fn shear_to_z_axis_makes_y_zero(v: Vec3<f64>) -> bool {
let s = calculate_shear_to_z_axis(&v);
apply_shear_to_z_axis(&v, &s).y().abs() < 0.00001
}
#[quickcheck]
fn shear_to_z_axis_leaves_z_unchanged(v: Vec3<f64>) -> bool {
let s = calculate_shear_to_z_axis(&v);
apply_shear_to_z_axis(&v, &s).z() == v.z()
}
}
mod barycentric_coordinates {
use super::*;
use quickcheck::TestResult;
use quickcheck_macros::quickcheck;
#[quickcheck]
fn sign_of_signed_edge_function_matches_winding(a: Vec3<f64>, b: Vec3<f64>) -> TestResult {
let a_2d = Vec2::new(a.x(), a.y());
let b_2d = Vec2::new(b.x(), b.y());
let c_2d = Vec2::new(0.0, 0.0);
let winding = (b_2d - a_2d).perp(&(c_2d - b_2d));
if winding.abs() < 0.00001 {
TestResult::discard()
} else {
let winding = winding.is_sign_positive();
let area_sign = signed_edge_function(&a, &b).is_sign_positive();
TestResult::from_bool(winding == area_sign)
}
}
#[quickcheck]
fn signed_edge_functions_has_same_result_as_signed_edge_function(
a: Vec3<f64>,
b: Vec3<f64>,
c: Vec3<f64>,
) -> bool {
let es = signed_edge_functions(&vec![a, b, c]);
es[0] == signed_edge_function(&b, &c)
&& es[1] == signed_edge_function(&c, &a)
&& es[2] == signed_edge_function(&a, &b)
}
#[quickcheck]
fn barycentric_coordinates_sum_to_one(a: Vec3<f64>, b: Vec3<f64>, c: Vec3<f64>) -> bool {
let barycentric_coordinates =
barycentric_coordinates_from_signed_edge_functions(signed_edge_functions(&vec![
a, b, c,
]));
(barycentric_coordinates
.coords
.iter()
.fold(0.0, |a, b| a + b)
- 1.0)
.abs()
< 0.00000001
}
}
mod triangle_intersect {
use super::*;
use crate::materials::LambertianMaterial;
use quickcheck::{Arbitrary, TestResult};
use quickcheck_macros::quickcheck;
#[test]
fn intersection_passes_with_ray_along_z_axis_ccw_winding() {
let target_triangle = Triangle {
vertices: [
Vec3::new(0.0, 1.0, 1.0),
Vec3::new(1.0, -1.0, 1.0),
Vec3::new(-1.0, -1.0, 1.0),
],
normals: [Vec3::zeros(); 3],
material: Arc::new(LambertianMaterial::new_dummy()),
};
let target_ray = Ray::new(Vec3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0));
if let None = target_triangle.intersect(&target_ray) {
panic!()
}
}
#[test]
fn intersection_passes_with_ray_along_z_axis_cw_winding() {
let target_triangle = Triangle {
vertices: [
Vec3::new(0.0, 1.0, 1.0),
Vec3::new(-1.0, -1.0, 1.0),
Vec3::new(1.0, -1.0, 1.0),
],
normals: [Vec3::zeros(); 3],
material: Arc::new(LambertianMaterial::new_dummy()),
};
let target_ray = Ray::new(Vec3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0));
if let None = target_triangle.intersect(&target_ray) {
panic!()
}
}
#[test]
fn intersection_passes_with_ray_along_nagative_z_axis_ccw_winding() {
let target_triangle = Triangle {
vertices: [
Vec3::new(0.0, 1.0, -1.0),
Vec3::new(1.0, -1.0, -1.0),
Vec3::new(-1.0, -1.0, -1.0),
],
normals: [Vec3::zeros(); 3],
material: Arc::new(LambertianMaterial::new_dummy()),
};
let target_ray = Ray::new(Vec3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, -1.0));
if let None = target_triangle.intersect(&target_ray) {
panic!()
}
}
#[test]
fn intersection_passes_with_ray_along_negativez_axis_cw_winding() {
let target_triangle = Triangle {
vertices: [
Vec3::new(0.0, 1.0, -1.0),
Vec3::new(-1.0, -1.0, -1.0),
Vec3::new(1.0, -1.0, -1.0),
],
normals: [Vec3::zeros(); 3],
material: Arc::new(LambertianMaterial::new_dummy()),
};
let target_ray = Ray::new(Vec3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, -1.0));
if let None = target_triangle.intersect(&target_ray) {
panic!()
}
}
#[test]
fn intersection_passes_with_ray_along_z_axis_but_translated_ccw_winding() {
let target_triangle = Triangle {
vertices: [
Vec3::new(5.0, 6.0, 6.0),
Vec3::new(6.0, 4.0, 6.0),
Vec3::new(4.0, 4.0, 6.0),
],
normals: [Vec3::zeros(); 3],
material: Arc::new(LambertianMaterial::new_dummy()),
};
let target_ray = Ray::new(Vec3::new(5.0, 5.0, 5.0), Vec3::new(0.0, 0.0, 1.0));
if let None = target_triangle.intersect(&target_ray) {
panic!()
}
}
#[test]
fn intersection_passes_with_ray_at_angle_to_z_axisand_translated_ccw_winding() {
let target_triangle = Triangle {
vertices: [
Vec3::new(6.0, 6.5, 6.0),
Vec3::new(7.0, 4.5, 6.0),
Vec3::new(5.0, 4.5, 6.0),
],
normals: [Vec3::zeros(); 3],
material: Arc::new(LambertianMaterial::new_dummy()),
};
let target_ray = Ray::new(Vec3::new(5.0, 5.0, 5.0), Vec3::new(1.0, 0.5, 1.0));
if let None = target_triangle.intersect(&target_ray) {
panic!()
}
}
fn intersect_with_centroid_and_test_result<
F: Fn(Option<IntersectionInfo>, Vec3<f64>) -> bool,
>(
vertex0: Vec3<f64>,
vertex1: Vec3<f64>,
vertex2: Vec3<f64>,
ray_origin: Vec3<f64>,
test: F,
) -> TestResult {
let centroid: Vec3<f64> = [vertex0, vertex1, vertex2]
.iter()
.fold(Vec3::new(0.0, 0.0, 0.0), |acc, &elem| acc + elem)
* (1.0 / 3.0);
let ray_direction = (centroid - ray_origin).normalize();
let normal = (vertex1 - vertex0).cross(&(vertex2 - vertex0)).normalize();
if normal.dot(&ray_direction).abs() < 0.000_000_1 {
//Discard if triangle is too close to edge-on
return TestResult::discard();
}
let target_triangle = Triangle {
vertices: [
Vec3::from(vertex0),
Vec3::from(vertex1),
Vec3::from(vertex2),
],
normals: [normal; 3],
material: Arc::new(LambertianMaterial::new_dummy()),
};
let ray = Ray::new(ray_origin, ray_direction);
TestResult::from_bool(test(target_triangle.intersect(&ray), centroid))
}
#[quickcheck]
fn intersection_with_centroid_hits(
vertex0: Vec3<f64>,
vertex1: Vec3<f64>,
vertex2: Vec3<f64>,
ray_origin: Vec3<f64>,
) -> TestResult {
let centroid: Vec3<f64> = [vertex0, vertex1, vertex2]
.iter()
.fold(Vec3::new(0.0, 0.0, 0.0), |acc, &elem| acc + elem)
* (1.0 / 3.0);
let ray_direction = (centroid - ray_origin).normalize();
let normal = (vertex1 - vertex0).cross(&(vertex2 - vertex0)).normalize();
if normal.dot(&ray_direction).abs() < 0.000_000_1 {
//Discard if triangle is too close to edge-on
return TestResult::discard();
}
let target_triangle = Triangle {
vertices: [
Vec3::from(vertex0),
Vec3::from(vertex1),
Vec3::from(vertex2),
],
normals: [normal; 3],
material: Arc::new(LambertianMaterial::new_dummy()),
};
let ray = Ray::new(ray_origin, ray_direction);
if let Some(_) = target_triangle.intersect(&ray) {
TestResult::passed()
} else {
TestResult::failed()
}
}
#[quickcheck]
fn intersection_with_centroid_hits_centroid(
vertex0: Vec3<f64>,
vertex1: Vec3<f64>,
vertex2: Vec3<f64>,
ray_origin: Vec3<f64>,
) -> TestResult {
intersect_with_centroid_and_test_result(
vertex0,
vertex1,
vertex2,
ray_origin,
|result, centroid| {
if let Some(IntersectionInfo { location, .. }) = result {
(location - centroid).norm() < 0.000_000_1
} else {
false
}
},
)
}
#[quickcheck]
fn intersection_with_centroid_hits_at_expected_distance(
vertex0: Vec3<f64>,
vertex1: Vec3<f64>,
vertex2: Vec3<f64>,
ray_origin: Vec3<f64>,
) -> TestResult {
intersect_with_centroid_and_test_result(
vertex0,
vertex1,
vertex2,
ray_origin,
|result, centroid| {
if let Some(IntersectionInfo { distance, .. }) = result {
((ray_origin - centroid).norm() - distance).abs() < 0.000_000_1
} else {
false
}
},
)
}
#[quickcheck]
fn intersection_with_centroid_has_expected_normal(
vertex0: Vec3<f64>,
vertex1: Vec3<f64>,
vertex2: Vec3<f64>,
ray_origin: Vec3<f64>,
) -> TestResult {
intersect_with_centroid_and_test_result(
vertex0,
vertex1,
vertex2,
ray_origin,
|result, _| {
if let Some(IntersectionInfo { normal, .. }) = result {
(normal - (vertex1 - vertex0).cross(&(vertex2 - vertex0)).normalize())
.norm()
< 0.000_000_1
} else {
false
}
},
)
}
#[quickcheck]
fn intersection_with_centroid_has_expected_retro(
vertex0: Vec3<f64>,
vertex1: Vec3<f64>,
vertex2: Vec3<f64>,
ray_origin: Vec3<f64>,
) -> TestResult {
intersect_with_centroid_and_test_result(
vertex0,
vertex1,
vertex2,
ray_origin,
|result, centroid| {
let expected_retro = (ray_origin - centroid).normalize();
if let Some(IntersectionInfo { retro, .. }) = result {
(expected_retro - retro).norm() < 0.000_000_1
} else {
false
}
},
)
}
#[derive(Clone, Copy, Debug)]
struct BarycentricCoords {
alpha: f64,
beta: f64,
gamma: f64,
}
impl quickcheck::Arbitrary for BarycentricCoords {
fn arbitrary<G: quickcheck::Gen>(g: &mut G) -> Self {
let e = 0.000_000_1;
let alpha = <f64 as Arbitrary>::arbitrary(g).abs().fract() * (1.0 - e) + e;
let beta = <f64 as Arbitrary>::arbitrary(g).abs().fract() * (1.0 - alpha) + e;
let gamma = 1.0 - (alpha + beta);
BarycentricCoords { alpha, beta, gamma }
}
}
fn intersect_with_barycentric_and_test_result<
F: Fn(Option<IntersectionInfo>, Vec3<f64>) -> bool,
>(
vertex0: Vec3<f64>,
vertex1: Vec3<f64>,
vertex2: Vec3<f64>,
ray_origin: Vec3<f64>,
barycentric_coords: BarycentricCoords,
test: F,
) -> TestResult {
let point = vertex0 * barycentric_coords.alpha
+ vertex1 * barycentric_coords.beta
+ vertex2 * barycentric_coords.gamma;
let ray_direction = (point - ray_origin).normalize();
let normal = (vertex1 - vertex0).cross(&(vertex2 - vertex0)).normalize();
if normal.dot(&ray_direction).abs() < 0.000_000_1 {
//Discard if triangle is too close to edge-on
return TestResult::discard();
}
let target_triangle = Triangle {
vertices: [
Vec3::from(vertex0),
Vec3::from(vertex1),
Vec3::from(vertex2),
],
normals: [normal; 3],
material: Arc::new(LambertianMaterial::new_dummy()),
};
let ray = Ray::new(ray_origin, ray_direction);
TestResult::from_bool(test(target_triangle.intersect(&ray), point))
}
#[quickcheck]
fn point_with_arbitrary_barycentric_coords_hits(
vertex0: Vec3<f64>,
vertex1: Vec3<f64>,
vertex2: Vec3<f64>,
ray_origin: Vec3<f64>,
barycentric_coords: BarycentricCoords,
) -> TestResult {
intersect_with_barycentric_and_test_result(
vertex0,
vertex1,
vertex2,
ray_origin,
barycentric_coords,
|result, _point| {
if let Some(_) = result {
true
} else {
false
}
},
)
}
#[quickcheck]
fn point_with_arbitrary_barycentric_coords_has_expected_normal(
vertex0: Vec3<f64>,
vertex1: Vec3<f64>,
vertex2: Vec3<f64>,
ray_origin: Vec3<f64>,
barycentric_coords: BarycentricCoords,
) -> TestResult {
intersect_with_barycentric_and_test_result(
vertex0,
vertex1,
vertex2,
ray_origin,
barycentric_coords,
|result, _point| {
let expected_normal =
(vertex1 - vertex0).cross(&(vertex2 - vertex0)).normalize();
if let Some(IntersectionInfo { normal, .. }) = result {
(normal - expected_normal).norm().abs() < 0.000_01
} else {
false
}
},
)
}
#[quickcheck]
fn point_with_arbitrary_barycentric_coords_has_expected_distance(
vertex0: Vec3<f64>,
vertex1: Vec3<f64>,
vertex2: Vec3<f64>,
ray_origin: Vec3<f64>,
barycentric_coords: BarycentricCoords,
) -> TestResult {
intersect_with_barycentric_and_test_result(
vertex0,
vertex1,
vertex2,
ray_origin,
barycentric_coords,
|result, point| {
let expected_distance = (point - ray_origin).norm();
if let Some(IntersectionInfo { distance, .. }) = result {
(distance - expected_distance).abs() < 0.000_01
} else {
false
}
},
)
}
#[quickcheck]
fn point_with_arbitrary_barycentric_coords_has_expected_retro(
vertex0: Vec3<f64>,
vertex1: Vec3<f64>,
vertex2: Vec3<f64>,
ray_origin: Vec3<f64>,
barycentric_coords: BarycentricCoords,
) -> TestResult {
intersect_with_barycentric_and_test_result(
vertex0,
vertex1,
vertex2,
ray_origin,
barycentric_coords,
|result, point| {
let expected_retro = (ray_origin - point).normalize();
if let Some(IntersectionInfo { retro, .. }) = result {
(retro - expected_retro).norm().abs() < 0.000_01
} else {
false
}
},
)
}
#[quickcheck]
fn intersection_fails_when_ray_outside_first_edge(
vertex0: Vec3<f64>,
vertex1: Vec3<f64>,
vertex2: Vec3<f64>,
ray_origin: Vec3<f64>,
uv: Vec2<f64>,
) -> bool {
let uv_origin = Vec3::from(vertex0);
let u_axis = (vertex1 - vertex0).normalize();
let w_axis = (vertex2 - vertex0).cross(&u_axis).normalize();
let v_axis = w_axis.cross(&u_axis);
let target_point = uv_origin + u_axis * uv.x() + v_axis * uv.y().abs();
let ray = Ray {
origin: ray_origin,
direction: (target_point - ray_origin).normalize(),
};
let triangle = Triangle {
vertices: [vertex0, vertex1, vertex2],
normals: [Vec3::zeros(); 3],
material: Arc::new(LambertianMaterial::new_dummy()),
};
match triangle.intersect(&ray) {
Some(_) => false,
None => true,
}
}
#[quickcheck]
fn intersection_fails_when_ray_outside_second_edge(
vertex0: Vec3<f64>,
vertex1: Vec3<f64>,
vertex2: Vec3<f64>,
ray_origin: Vec3<f64>,
uv: Vec2<f64>,
) -> bool {
let uv_origin = Vec3::from(vertex0);
let u_axis = (vertex2 - vertex1).normalize();
let w_axis = (vertex1 - vertex0).cross(&u_axis).normalize();
let v_axis = w_axis.cross(&u_axis);
let target_point = uv_origin + u_axis * uv.x() + v_axis * uv.y().abs();
let ray = Ray {
origin: ray_origin,
direction: (target_point - ray_origin).normalize(),
};
let triangle = Triangle {
vertices: [vertex0, vertex1, vertex2],
normals: [Vec3::zeros(); 3],
material: Arc::new(LambertianMaterial::new_dummy()),
};
match triangle.intersect(&ray) {
Some(_) => false,
None => true,
}
}
#[quickcheck]
fn intersection_fails_when_ray_outside_third_edge(
vertex0: Vec3<f64>,
vertex1: Vec3<f64>,
vertex2: Vec3<f64>,
ray_origin: Vec3<f64>,
uv: Vec2<f64>,
) -> bool {
let uv_origin = Vec3::from(vertex0);
let u_axis = (vertex0 - vertex2).normalize();
let w_axis = (vertex1 - vertex2).cross(&u_axis).normalize();
let v_axis = w_axis.cross(&u_axis);
let target_point = uv_origin + u_axis * uv.x() + v_axis * uv.y().abs();
let ray = Ray {
origin: ray_origin,
direction: (target_point - ray_origin).normalize(),
};
let triangle = Triangle {
vertices: [vertex0, vertex1, vertex2],
normals: [Vec3::zeros(); 3],
material: Arc::new(LambertianMaterial::new_dummy()),
};
match triangle.intersect(&ray) {
Some(_) => false,
None => true,
}
}
#[quickcheck]
fn intersection_fails_when_triangle_is_behind_ray(
vertex0: Vec3<f64>,
vertex1: Vec3<f64>,
vertex2: Vec3<f64>,
ray_origin: Vec3<f64>,
barycentric_coords: BarycentricCoords,
) -> bool {
let point_behind_ray = vertex0 * barycentric_coords.alpha
+ vertex1 * barycentric_coords.beta
+ vertex2 * barycentric_coords.gamma;
let ray = Ray {
origin: ray_origin,
direction: (ray_origin - point_behind_ray).normalize(),
};
let triangle = Triangle {
vertices: [vertex0, vertex1, vertex2],
normals: [Vec3::zeros(); 3],
material: Arc::new(LambertianMaterial::new_dummy()),
};
match triangle.intersect(&ray) {
Some(_) => false,
None => true,
}
}
}
}

29
src/raycasting/vec_aggregate.rs
View File

@ -1,29 +0,0 @@
use super::{BoundingBox, HasBoundingBox, Intersect, IntersectionInfo, Primitive, Ray};
use crate::math::Affine3;
impl HasBoundingBox for Vec<Box<dyn Primitive>> {
fn bounding_box(&self) -> BoundingBox {
self.iter().fold(BoundingBox::empty(), |acc, elem| {
acc.union(&elem.bounding_box())
})
}
}
impl Intersect for Vec<Box<dyn Primitive>> {
fn intersect(&self, ray: &Ray) -> Option<IntersectionInfo> {
self.iter()
.flat_map(|primitive| primitive.intersect(ray))
.min_by(
|a, b| match PartialOrd::partial_cmp(&a.distance, &b.distance) {
None => std::cmp::Ordering::Less,
Some(ordering) => ordering,
},
)
}
}
impl Primitive for Vec<Box<dyn Primitive>> {
fn transform(&self, transformation: &Affine3<f64>) -> Box<dyn Primitive> {
todo!()
}
}

58
src/realtype.rs
View File

@ -1,58 +0,0 @@
pub trait NormalizedToU32 {
fn normalized_to_u32(self, num_bits: usize) -> u32;
}
impl NormalizedToU32 for f32 {
fn normalized_to_u32(self, num_bits: usize) -> u32 {
let scale = (num_bits as f32).exp2() - 1.0;
(self * scale) as u32
}
}
impl NormalizedToU32 for f64 {
fn normalized_to_u32(self, num_bits: usize) -> u32 {
let scale = (num_bits as f64).exp2() - 1.0;
(self * scale) as u32
}
}
#[cfg(test)]
mod test {
use super::*;
#[test]
fn zero_f32_normalizes_to_zero() {
let target = 0.0f32;
assert!(target.normalized_to_u32(10) == 0);
}
#[test]
fn one_f32_normalizes_to_all_ones() {
let target = 1.0f32;
assert!(target.normalized_to_u32(10) == 0b1111111111);
}
#[test]
fn half_f32_normalizes_to_half_value() {
let target = 0.5f32;
assert!(target.normalized_to_u32(10) == 511);
}
#[test]
fn zero_f64_normalizes_to_zero() {
let target = 0.0f64;
assert!(target.normalized_to_u32(10) == 0);
}
#[test]
fn one_f64_normalizes_to_all_ones() {
let target = 1.0f64;
assert!(target.normalized_to_u32(10) == 0b1111111111);
}
#[test]
fn half_f64_normalizes_to_half_value() {
let target = 0.5f64;
assert!(target.normalized_to_u32(10) == 511);
}
}

12
src/sampler.rs
View File

@ -1,16 +1,18 @@
use super::raycasting::{IntersectionInfo, Ray}; use super::raycasting::{IntersectionInfo, Ray};
use super::scene::Scene; use super::scene::Scene;
pub struct Sampler<'a> { use nalgebra::RealField;
pub scene: &'a Scene,
pub struct Sampler<'a, T: RealField> {
pub scene: &'a Scene<T>,
} }
impl<'a> Sampler<'a> { impl<'a, T: RealField> Sampler<'a, T> {
pub fn sample(&self, ray: &Ray) -> Option<IntersectionInfo> { pub fn sample(&self, ray: &Ray<T>) -> Option<IntersectionInfo<T>> {
self.scene self.scene
.objects .objects
.iter() .iter()
.flat_map(|object| object.intersect(ray)) .flat_map(|object| object.intersect(&ray))
.min_by( .min_by(
|a, b| match PartialOrd::partial_cmp(&a.distance, &b.distance) { |a, b| match PartialOrd::partial_cmp(&a.distance, &b.distance) {
None => std::cmp::Ordering::Less, None => std::cmp::Ordering::Less,

10
src/scene.rs
View File

@ -1,8 +1,8 @@
use crate::math::Vec3; use nalgebra::{RealField, Vector3};
use crate::raycasting::Primitive; use crate::raycasting::Intersect;
pub struct Scene { pub struct Scene<T: RealField> {
pub camera_location: Vec3<f64>, pub camera_location: Vector3<T>,
pub objects: Vec<Box<dyn Primitive>>, pub objects: Vec<Box<dyn Intersect<T>>>,
} }

51
src/util/algebra_utils.rs
View File

@ -1,51 +0,0 @@
use crate::math::{Mat3, Vec3};
pub fn try_change_of_basis_matrix(x: &Vec3<f64>, y: &Vec3<f64>, z: &Vec3<f64>) -> Option<Mat3<f64>> {
Some(Mat3::from_rows(x, y, z))
}
#[cfg(test)]
mod tests {
use super::*;
#[cfg(test)]
mod change_of_basis_matrix {
use super::*;
use quickcheck_macros::quickcheck;
#[test]
fn produces_isentity_when_passed_axes() {
let target: Mat3<f64> =
try_change_of_basis_matrix(&Vec3::unit_x(), &Vec3::unit_y(), &Vec3::unit_z())
.unwrap();
assert!(target == Mat3::identity())
}
#[quickcheck]
fn swap_xy_does_not_change_z(v: Vec3<f64>) {
let target: Mat3<f64> =
try_change_of_basis_matrix(&Vec3::unit_y(), &Vec3::unit_x(), &Vec3::unit_z())
.unwrap();
let v2 = target * v;
assert!(v2.z() == v.z())
}
#[quickcheck]
fn swap_xy_copies_y_to_x(v: Vec3<f64>) {
let target: Mat3<f64> =
try_change_of_basis_matrix(&Vec3::unit_y(), &Vec3::unit_x(), &Vec3::unit_z())
.unwrap();
let v2 = target * v;
assert!(v2.x() == v.y())
}
#[quickcheck]
fn swap_xy_copies_x_to_y(v: Vec3<f64>) {
let target: Mat3<f64> =
try_change_of_basis_matrix(&Vec3::unit_y(), &Vec3::unit_x(), &Vec3::unit_z())
.unwrap();
let v2 = target * v;
assert!(v2.y() == v.x())
}
}
}

130
src/util/array2d.rs
View File

@ -1,130 +0,0 @@
use std::ops::{Index, IndexMut};
/// 3D row-major dynamic array
#[derive(Clone, Debug)]
pub struct Array2D<T> {
data: Vec<T>,
height: usize,
width: usize,
}
impl<T: Default + Clone> Array2D<T> {
/// Create Array2D with `width` and `height` filled with default values.
pub fn new(height: usize, width: usize) -> Array2D<T> {
Array2D {
data: vec![Default::default(); width * height],
height,
width,
}
}
/// Reset contents of array to all default values
pub fn clear(&mut self) {
self.data.clear();
}
}
impl<T> Array2D<T> {
pub fn get_height(&self) -> usize {
self.height
}
pub fn get_width(&self) -> usize {
self.width
}
/// Return single slice containing all elements (row-major)
pub fn as_slice(&self) -> &[T] {
self.data.as_slice()
}
}
impl<T: Copy> Array2D<T> {
/// Copy `source` into the Array2D at the specified location
pub fn update_block(&mut self, start_row: usize, start_column: usize, source: &Array2D<T>) {
let end_row = start_row + source.height;
let end_column = start_column + source.width;
assert!(end_row <= self.height);
for i in 0..source.height {
self[start_row + i][start_column..end_column].copy_from_slice(&source[i])
}
}
}
impl<T> Index<usize> for Array2D<T> {
type Output = [T];
fn index(&self, index: usize) -> &[T] {
assert!(index < self.height);
&self.data[(index * self.width)..((index + 1) * self.width)]
}
}
impl<T> IndexMut<usize> for Array2D<T> {
fn index_mut(&mut self, index: usize) -> &mut [T] {
assert!(index < self.height);
&mut self.data[(index * self.width)..((index + 1) * self.width)]
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn default_u8_is_all_zeros() {
let target: Array2D<u8> = Array2D::new(10, 12);
for i in 0..10 {
for j in 0..12 {
assert!(target[i][j] == 0);
}
}
}
#[test]
#[should_panic]
fn panics_if_row_outside_bounds() {
let target: Array2D<u8> = Array2D::new(10, 12);
assert!(target[10][6] == 0);
}
#[test]
#[should_panic]
fn panics_if_column_outside_bounds() {
let target: Array2D<u8> = Array2D::new(10, 12);
assert!(target[5][12] == 0);
}
#[test]
fn write_and_read_all_preserves_values() {
let mut target: Array2D<u8> = Array2D::new(10, 12);
for i in 0..10 {
for j in 0..12 {
target[i][j] = (i * 10 + j) as u8;
}
}
for i in 0..10 {
for j in 0..12 {
assert!(target[i][j] == (i * 10 + j) as u8);
}
}
}
#[test]
fn update_block_writes_expected_values_in_block() {
let mut target: Array2D<u8> = Array2D::new(10, 12);
let mut source: Array2D<u8> = Array2D::new(2, 3);
for i in 0..2 {
for j in 0..3 {
source[i][j] = (i * 3 + j) as u8;
}
}
target.update_block(4, 5, &source);
dbg!(&target);
for i in 0..2 {
for j in 0..3 {
assert!(dbg!(target[4 + i][5 + j]) == dbg!(source[i][j]));
}
}
}
}

240
src/util/axis_aligned_bounding_box.rs
View File

@ -1,240 +0,0 @@
use crate::math::Vec3;
use crate::util::Interval;
use itertools::izip;
#[derive(Debug, Clone, Copy)]
pub struct BoundingBox {
pub bounds: [Interval; 3],
}
impl BoundingBox {
pub fn from_corners(a: Vec3<f64>, b: Vec3<f64>) -> Self {
let mut result = BoundingBox {
bounds: [Interval::infinite(); 3],
};
for (bounds_elem, a_elem, b_elem) in
izip!(result.bounds.iter_mut(), a.coords.iter(), b.coords.iter())
{
*bounds_elem = Interval::new(*a_elem, *b_elem);
}
result
}
pub fn empty() -> Self {
BoundingBox {
bounds: [Interval::empty(), Interval::empty(), Interval::empty()],
}
}
pub fn from_point(p: Vec3<f64>) -> Self {
BoundingBox {
bounds: [
Interval::degenerate(p.x()),
Interval::degenerate(p.y()),
Interval::degenerate(p.z()),
],
}
}
pub fn from_points<'a, I>(points: I) -> Self
where
I: IntoIterator<Item = &'a Vec3<f64>>,
{
points
.into_iter()
.fold(BoundingBox::empty(), |acc, p| acc.expand_to_point(p))
}
pub fn expand_to_point(&self, p: &Vec3<f64>) -> Self {
BoundingBox {
bounds: [
self.bounds[0].expand_to_value(p.x()),
self.bounds[1].expand_to_value(p.y()),
self.bounds[2].expand_to_value(p.z()),
],
}
}
pub fn contains_point(&self, p: Vec3<f64>) -> bool {
self.bounds
.iter()
.zip(p.coords.iter())
.all(|(interval, &value)| interval.contains_value(value))
}
pub fn union(&self, other: &BoundingBox) -> BoundingBox {
BoundingBox {
bounds: [
self.bounds[0].union(other.bounds[0]),
self.bounds[1].union(other.bounds[1]),
self.bounds[2].union(other.bounds[2]),
],
}
}
pub fn largest_dimension(&self) -> usize {
let (dimension, _) = self
.bounds
.iter()
.enumerate()
.map(|(index, elem)| {
(
index,
if elem.is_degenerate() {
-1.0
} else {
elem.get_max() - elem.get_min()
},
)
})
.fold((0, 0.0), |(acc, acc_size), (elem, elem_size)| {
if elem_size > acc_size {
(elem, elem_size)
} else {
(acc, acc_size)
}
});
dimension
}
}
#[cfg(test)]
mod tests {
use super::*;
use quickcheck_macros::quickcheck;
#[test]
fn from_corners_with_same_point_yields_degenerate_intervals() {
let test_point = Vec3::new(0f64, 1f64, 2f64);
let target = BoundingBox::from_corners(test_point, test_point);
assert!(target.bounds.iter().all(|e| e.is_degenerate()));
}
#[test]
fn from_corners_yields_same_result_with_any_oposite_corners() {
let corner_000 = Vec3::new(0.0, 0.0, 0.0);
let corner_001 = Vec3::new(0.0, 0.0, 1.0);
let corner_010 = Vec3::new(0.0, 1.0, 0.0);
let corner_011 = Vec3::new(0.0, 1.0, 1.0);
let corner_100 = Vec3::new(1.0, 0.0, 0.0);
let corner_101 = Vec3::new(1.0, 0.0, 1.0);
let corner_110 = Vec3::new(1.0, 1.0, 0.0);
let corner_111 = Vec3::new(1.0, 1.0, 1.0);
let test_inputs: Vec<(Vec3<f64>, Vec3<f64>)> = vec![
(corner_000, corner_111),
(corner_001, corner_110),
(corner_010, corner_101),
(corner_011, corner_100),
(corner_100, corner_011),
(corner_101, corner_010),
(corner_110, corner_001),
(corner_111, corner_000),
];
for (a, b) in test_inputs {
let target = BoundingBox::from_corners(a, b);
assert!(target
.bounds
.iter()
.all(|bounds| bounds.get_min() == 0.0 && bounds.get_max() == 1.0));
}
}
#[quickcheck]
fn union_with_self_yields_self(a: Vec3<f64>, b: Vec3<f64>) -> bool {
let target = BoundingBox::from_corners(a, b);
let result = target.union(&target);
target
.bounds
.iter()
.zip(result.bounds.iter())
.all(|(a, b)| a.get_min() == b.get_min() && a.get_max() == b.get_max())
}
#[quickcheck]
fn union_yields_full_ranges(a: Vec3<f64>, b: Vec3<f64>, c: Vec3<f64>, d: Vec3<f64>) -> bool {
let target1 = BoundingBox::from_corners(a, b);
let target2 = BoundingBox::from_corners(c, d);
let result = target1.union(&target2);
izip!(
result.bounds.iter(),
target1.bounds.iter(),
target2.bounds.iter()
)
.all(|(r, t1, t2)| {
r.get_min() <= t1.get_min()
&& r.get_min() <= t2.get_min()
&& r.get_max() >= t1.get_max()
&& r.get_max() >= t2.get_max()
})
}
#[quickcheck]
fn empty_box_contains_no_points(p: Vec3<f64>) -> bool {
let target = BoundingBox::empty();
!target.contains_point(p)
}
#[quickcheck]
fn from_points_produces_box_that_contains_only_points_bounded_by_inputs_on_all_axes(
p: Vec3<f64>,
a: Vec3<f64>,
b: Vec3<f64>,
c: Vec3<f64>,
d: Vec3<f64>,
e: Vec3<f64>,
) -> bool {
let points = vec![a, b, c, d, e];
let target = BoundingBox::from_points(&points);
let is_in_bounds = points.iter().any(|elem| elem.x() >= p.x())
&& points.iter().any(|elem| elem.x() <= p.x())
&& points.iter().any(|elem| elem.y() >= p.y())
&& points.iter().any(|elem| elem.y() <= p.y())
&& points.iter().any(|elem| elem.z() >= p.z())
&& points.iter().any(|elem| elem.z() <= p.z());
target.contains_point(p) == is_in_bounds
}
#[quickcheck]
fn no_dimension_is_larger_than_largest_dimension(
a: f64,
b: f64,
c: f64,
d: f64,
e: f64,
f: f64,
) -> bool {
let target = BoundingBox {
bounds: [
if a > b {
Interval::empty()
} else {
Interval::new(a, b)
},
if c > d {
Interval::empty()
} else {
Interval::new(c, d)
},
if e > f {
Interval::empty()
} else {
Interval::new(e, f)
},
],
};
let largest_dimension = target.largest_dimension();
let largest_bounds = target.bounds[largest_dimension];
if largest_bounds.is_empty() {
target.bounds.iter().all(|elem| elem.is_empty())
} else {
let largest_size = largest_bounds.get_max() - largest_bounds.get_min();
target
.bounds
.iter()
.all(|elem| elem.is_empty() || !(largest_size < elem.get_max() - elem.get_min()))
}
}
}

25
src/util/binary_tree.rs
View File

@ -1,25 +0,0 @@
pub enum BinaryTree<Value, LeafValue> {
Branch {
value: Value,
left: Box<Self>,
right: Box<Self>,
},
Leaf {
value: LeafValue,
},
None,
}
impl<Value, LeafValue> BinaryTree<Value, LeafValue> {
pub fn count_leaves(&self) -> usize {
match self {
Self::Branch {
value: _,
left,
right,
} => right.count_leaves() + left.count_leaves(),
Self::Leaf { value: _ } => 1,
Self::None => 0,
}
}
}

328
src/util/interval.rs
View File

@ -1,328 +0,0 @@
#[derive(Debug, Clone, Copy)]
pub struct Interval {
min: f64,
max: f64,
}
impl Interval {
pub fn new(a: f64, b: f64) -> Self {
if a > b {
Interval { min: b, max: a }
} else {
Interval { min: a, max: b }
}
}
pub fn empty() -> Self {
Interval {
min: f64::INFINITY,
max: f64::NEG_INFINITY,
}
}
pub fn infinite() -> Self {
Interval {
min: f64::NEG_INFINITY,
max: f64::INFINITY,
}
}
pub fn degenerate(value: f64) -> Self {
Interval {
min: value,
max: value,
}
}
pub fn get_min(&self) -> f64 {
self.min
}
pub fn get_max(&self) -> f64 {
self.max
}
pub fn is_degenerate(self) -> bool {
self.min == self.max
}
pub fn is_empty(self) -> bool {
self.min > self.max
}
pub fn contains_value(&self, value: f64) -> bool {
value >= self.min && value <= self.max
}
pub fn intersection(self, b: Self) -> Self {
Interval {
min: self.min.max(b.min),
max: self.max.min(b.max),
}
}
pub fn union(self, b: Self) -> Self {
if self.is_empty() {
b
} else if b.is_empty() {
self
} else {
Interval {
min: self.min.min(b.min),
max: self.max.max(b.max),
}
}
}
pub fn expand_to_value(self, v: f64) -> Self {
if self.is_empty() {
Interval::degenerate(v)
} else {
Interval {
min: self.min.min(v),
max: self.max.max(v),
}
}
}
}
#[cfg(test)]
mod tests {
use super::*;
use itertools::{Itertools, MinMaxResult};
use quickcheck_macros::quickcheck;
#[test]
fn never_constructed_empty() {
let target1 = Interval::new(5f64, 10f64);
assert!(!target1.is_empty());
let target2 = Interval::new(10f64, 5f64);
assert!(!target2.is_empty());
let target1 = Interval::new(5f64, -10f64);
assert!(!target1.is_empty());
let target2 = Interval::new(10f64, -5f64);
assert!(!target2.is_empty());
let target1 = Interval::new(-5f64, 10f64);
assert!(!target1.is_empty());
let target2 = Interval::new(-10f64, 5f64);
assert!(!target2.is_empty());
let target1 = Interval::new(-5f64, -10f64);
assert!(!target1.is_empty());
let target2 = Interval::new(-10f64, -5f64);
assert!(!target2.is_empty());
}
#[test]
fn empty_is_empty() {
let target: Interval = Interval::empty();
assert!(target.is_empty());
}
#[test]
fn empty_when_min_greater_than_max() {
let target = Interval {
min: 10f64,
max: 5f64,
};
assert!(target.is_empty());
}
#[test]
fn empty_when_min_greater_than_max_with_negative_values() {
let target = Interval {
min: -5f64,
max: -10f64,
};
assert!(target.is_empty());
}
#[test]
fn empty_when_min_greater_than_max_with_mixed_signs() {
let target = Interval {
min: 5f64,
max: -10f64,
};
assert!(target.is_empty());
}
#[test]
fn empty_is_not_degenerate() {
let target = Interval {
min: 10f64,
max: 5f64,
};
assert!(!target.is_degenerate());
}
#[test]
fn empty_is_not_degenerate_with_negative_values() {
let target = Interval {
min: -5f64,
max: -10f64,
};
assert!(!target.is_degenerate());
}
#[test]
fn empty_is_not_degenerate_with_mixed_signs() {
let target = Interval {
min: -5f64,
max: 10f64,
};
assert!(!target.is_degenerate());
}
#[quickcheck]
fn no_value_is_in_interval_returned_by_emtpy(value: f64) -> bool {
!Interval::empty().contains_value(value)
}
#[test]
fn identical_min_max_yields_degenerate() {
let target = Interval {
min: 5f64,
max: 5f64,
};
assert!(target.is_degenerate());
}
#[test]
fn degenerate_is_degenerate() {
let target = Interval::degenerate(5f64);
assert!(target.is_degenerate());
}
#[test]
fn degenerate_is_not_empty() {
let target = Interval {
min: 5f64,
max: 5f64,
};
assert!(target.is_degenerate());
}
#[test]
fn degenerate_is_degenerate_with_negative_value() {
let target = Interval {
min: -5f64,
max: -5f64,
};
assert!(target.is_degenerate());
}
#[test]
fn degenerate_is_not_empty_with_negative_value() {
let target = Interval {
min: -5f64,
max: -5f64,
};
assert!(target.is_degenerate());
}
#[test]
fn degenerate_contains_expected_value() {
let target = Interval::degenerate(5f64);
assert!(target.contains_value(5.0));
}
#[quickcheck]
fn degenerate_does_not_contain_any_values_othter_than_expected_value(value: f64) -> bool {
let target_value = if value == 5f64 { 5.5 } else { 5f64 };
let target = Interval::degenerate(target_value);
!target.contains_value(value)
}
#[test]
fn intersection_with_infinite_is_self() {
let target = Interval::new(5.0, 10.0);
let result = target.intersection(Interval::infinite());
assert!(target.min == result.min);
assert!(target.max == result.max);
}
#[quickcheck]
fn union_with_self_yields_self(a: f64, b: f64) -> bool {
let target = Interval::new(a, b);
let result = target.union(target);
result.min == target.min && result.max == target.max
}
#[quickcheck]
fn union_yields_min_and_max(a: f64, b: f64, c: f64, d: f64) -> bool {
let values = vec![a, b, c, d];
if let MinMaxResult::MinMax(&min, &max) =
values.iter().minmax_by(|a, b| a.partial_cmp(b).unwrap())
{
let target1 = Interval::new(a, b);
let target2 = Interval::new(c, d);
let result = target1.union(target2);
result.min == min && result.max == max
} else {
false
}
}
#[test]
fn union_with_empty_interval_is_correct() {
let empty = Interval {
min: 1f64,
max: -1f64,
};
let not_empty = Interval {
min: 5f64,
max: 10f64,
};
let union1 = not_empty.union(empty);
assert!(union1.min == 5.0);
assert!(union1.max == 10.0);
let union2 = empty.union(not_empty);
assert!(union2.min == 5.0);
assert!(union2.max == 10.0);
}
#[test]
fn union_with_empty_interval_is_correct_when_empty_interval_produced_by_intersection() {
let empty = Interval {
min: 1f64,
max: -1f64,
};
let not_empty = Interval {
min: 5f64,
max: 10f64,
};
let union1 = not_empty.union(empty);
assert!(union1.min == 5.0);
assert!(union1.max == 10.0);
let union2 = empty.union(not_empty);
assert!(union2.min == 5.0);
assert!(union2.max == 10.0);
}
#[quickcheck]
pub fn expand_to_value_creates_interval_that_includes_value(
min: f64,
max: f64,
value: f64,
) -> bool {
// Don't check if min <= max, we want to test empty intervals too
let interval1 = Interval { min, max };
let interval2 = interval1.expand_to_value(value);
interval2.contains_value(value)
}
#[quickcheck]
pub fn expand_to_value_creates_interval_that_includes_original_interval(
b: f64,
a: f64,
value: f64,
) -> bool {
let interval1 = Interval::new(a, b);
let interval2 = interval1.expand_to_value(value);
let interval3 = interval2.intersection(interval1);
// If interval2 contains interval1, that the intersection of the two will
// be equal to interval1
interval1.min == interval3.min && interval1.max == interval3.max
}
}

13
src/util/mod.rs
View File

@ -1,13 +0,0 @@
mod interval;
pub use interval::Interval;
pub mod algebra_utils;
pub mod array2d;
pub use array2d::Array2D;
pub mod axis_aligned_bounding_box;
pub mod binary_tree;
pub mod morton;
pub mod normalizer;
mod tile_iterator;
pub use tile_iterator::{Tile, TileIterator};
pub mod polyhedra;

45
src/util/morton.rs
View File

@ -1,45 +0,0 @@
use crate::math::Vec3;
use crate::realtype::NormalizedToU32;
fn spread_bits(v: u32) -> u32 {
let mut result = 0;
for power in 0..9 {
result |= ((1 << power) & v) << (power * 2);
}
result
}
pub fn morton_order_value_3d(p: Vec3<f64>) -> u32 {
let x = p.x().normalized_to_u32(10);
let y = p.y().normalized_to_u32(10);
let z = p.z().normalized_to_u32(10);
(spread_bits(x) << 2) | (spread_bits(y) << 1) | spread_bits(z)
}
#[cfg(test)]
mod tests {
use super::*;
mod spread_bits {
use super::*;
#[test]
fn zero_yields_zero() {
assert!(spread_bits(0) == 0);
}
#[test]
fn one_yields_one() {
assert!(spread_bits(1) == 1);
}
#[test]
fn b1111_yields_b1001001001() {
assert!(spread_bits(0b1111) == 0b1001001001);
}
#[test]
fn b1010_yields_b1000001000() {
assert!(spread_bits(0b1010) == 0b1000001000);
}
}
}

178
src/util/normalizer.rs
View File

@ -1,178 +0,0 @@
use super::axis_aligned_bounding_box::BoundingBox;
use super::Interval;
use crate::math::Vec3;
use itertools::izip;
#[derive(Debug, Copy, Clone)]
pub struct RealNormalizer {
min: f64,
range: f64,
}
impl RealNormalizer {
pub fn new(interval: Interval) -> Self {
let min = interval.get_min();
let range = interval.get_max() - min;
Self { min, range }
}
pub fn normalize(&self, value: f64) -> f64 {
(value - self.min) / self.range
}
pub fn normalize_and_clamp(&self, value: f64) -> f64 {
((value - self.min) / self.range).clamp(0.0, 1.0)
}
}
#[derive(Debug)]
pub struct Point3Normalizer {
dimension_normalizers: [RealNormalizer; 3],
}
impl Point3Normalizer {
pub fn new(bounds: BoundingBox) -> Self {
let mut normalizer = Point3Normalizer {
dimension_normalizers: [RealNormalizer::new(Interval::empty()); 3],
};
for (normalizer, &bounds) in normalizer
.dimension_normalizers
.iter_mut()
.zip(bounds.bounds.iter())
{
*normalizer = RealNormalizer::new(bounds);
}
normalizer
}
pub fn normalize(&self, point: Vec3<f64>) -> Vec3<f64> {
let mut result = Vec3::new(0.0, 0.0, 0.0);
for (value_out, &value_in, normalizer) in izip!(
result.coords.iter_mut(),
point.coords.iter(),
self.dimension_normalizers.iter()
) {
*value_out = normalizer.normalize(value_in);
}
result
}
pub fn normalize_and_clamp(&self, point: Vec3<f64>) -> Vec3<f64> {
let mut result = Vec3::new(0.0, 0.0, 0.0);
for (value_out, &value_in, normalizer) in izip!(
result.coords.iter_mut(),
point.coords.iter(),
self.dimension_normalizers.iter()
) {
*value_out = normalizer.normalize_and_clamp(value_in);
}
result
}
}
#[cfg(test)]
mod test {
use super::*;
use quickcheck_macros::quickcheck;
#[quickcheck]
fn normalize_zero_to_one_yields_input(value: f64) -> bool {
let target = RealNormalizer::new(Interval::new(0.0, 1.0));
target.normalize(value) == value
}
#[quickcheck]
fn normalize_two_to_three_yields_input_minus_two(value: f64) -> bool {
let target = RealNormalizer::new(Interval::new(2.0, 3.0));
target.normalize(value) == value - 2.0
}
#[quickcheck]
fn normalize_negative_three_to_negative_two_yields_input_plus_three(value: f64) -> bool {
let target = RealNormalizer::new(Interval::new(-3.0, -2.0));
target.normalize(value) == value + 3.0
}
#[quickcheck]
fn normalize_zero_to_two_yields_input_divided_by_two(value: f64) -> bool {
let target = RealNormalizer::new(Interval::new(0.0, 2.0));
target.normalize(value) == value / 2.0
}
#[test]
fn normalize_two_to_four_yields_zero_when_input_is_two() {
let target = RealNormalizer::new(Interval::new(2.0, 4.0));
assert!(target.normalize(2.0) == 0.0)
}
#[test]
fn normalize_two_to_four_yields_one_when_input_is_four() {
let target = RealNormalizer::new(Interval::new(2.0, 4.0));
assert!(target.normalize(4.0) == 1.0)
}
#[quickcheck]
fn normalize_two_to_four_yields_input_divided_by_two_minus_one(value: f64) -> bool {
let target = RealNormalizer::new(Interval::new(2.0, 4.0));
target.normalize(value) == (value - 2.0) / 2.0
}
#[quickcheck]
fn normalize_and_clamp_two_to_four_yields_zero_when_input_less_than_or_equal_two(
value: f64,
) -> bool {
let target = RealNormalizer::new(Interval::new(2.0, 4.0));
target.normalize_and_clamp(value) == 0.0 || value > 2.0
}
#[quickcheck]
fn normalize_and_clamp_two_to_four_yields_one_when_input_greater_than_or_equal_four(
value: f64,
) -> bool {
let target = RealNormalizer::new(Interval::new(2.0, 4.0));
target.normalize_and_clamp(value) == 1.0 || value < 4.0
}
#[quickcheck]
fn normalize_and_clamp_two_to_four_yields_same_value_as_normalize_when_in_range(
value: f64,
) -> bool {
let target = RealNormalizer::new(Interval::new(2.0, 4.0));
target.normalize_and_clamp(value) == target.normalize(value) || value < 2.0 || value > 4.0
}
#[quickcheck]
fn normalize_point3_is_the_same_as_normalize_each_dimension(
a: Vec3<f64>,
b: Vec3<f64>,
c: Vec3<f64>,
) -> bool {
let x_normalizer = RealNormalizer::new(Interval::new(a.x().min(b.x()), a.x().max(b.x())));
let y_normalizer = RealNormalizer::new(Interval::new(a.y().min(b.y()), a.y().max(b.y())));
let z_normalizer = RealNormalizer::new(Interval::new(a.z().min(b.z()), a.z().max(b.z())));
let xyz_normalizer = Point3Normalizer::new(BoundingBox::from_corners(a, b));
let normalized_point = xyz_normalizer.normalize(c);
x_normalizer.normalize(c.x()) == normalized_point.x()
&& y_normalizer.normalize(c.y()) == normalized_point.y()
&& z_normalizer.normalize(c.z()) == normalized_point.z()
}
#[quickcheck]
fn normalize_and_clamp_point3_is_the_same_as_normalize_and_clamp_each_dimension(
a: Vec3<f64>,
b: Vec3<f64>,
c: Vec3<f64>,
) -> bool {
let x_normalizer = RealNormalizer::new(Interval::new(a.x().min(b.x()), a.x().max(b.x())));
let y_normalizer = RealNormalizer::new(Interval::new(a.y().min(b.y()), a.y().max(b.y())));
let z_normalizer = RealNormalizer::new(Interval::new(a.z().min(b.z()), a.z().max(b.z())));
let xyz_normalizer = Point3Normalizer::new(BoundingBox::from_corners(a, b));
let normalized_point = xyz_normalizer.normalize_and_clamp(c);
x_normalizer.normalize_and_clamp(c.x()) == normalized_point.x()
&& y_normalizer.normalize_and_clamp(c.y()) == normalized_point.y()
&& z_normalizer.normalize_and_clamp(c.z()) == normalized_point.z()
}
}

131
src/util/polyhedra.rs
View File

@ -1,131 +0,0 @@
use itertools::izip;
use crate::materials::Material;
use crate::math::Vec3;
use crate::raycasting::{Primitive, Triangle};
use std::sync::Arc;
pub fn triangulate_polygon(
vertices: &[Vec3<f64>],
normal: &Vec3<f64>,
material: Arc<dyn Material>,
) -> Vec<Arc<dyn Primitive>> {
assert!(vertices.len() >= 3);
let hinge = vertices[0];
izip!(vertices.iter().skip(1), vertices.iter().skip(2))
.map(|(a, b)| {
Arc::new(Triangle {
vertices: [hinge, *a, *b],
normals: [*normal, *normal, *normal],
material: Arc::clone(&material),
}) as Arc<dyn Primitive>
})
.collect()
}
pub fn generate_dodecahedron(
centre: Vec3<f64>,
size: f64,
material: Arc<dyn Material>,
) -> Vec<Arc<dyn Primitive>> {
let phi = (1.0 + (5.0_f64).sqrt()) / 2.0;
let phi_inv = 1.0 / phi;
let faces = vec![
vec![
Vec3::new(phi_inv, 0.0, phi),
Vec3::new(-phi_inv, 0.0, phi),
Vec3::new(-1.0, -1.0, 1.0),
Vec3::new(0.0, -phi, phi_inv),
Vec3::new(1.0, -1.0, 1.0),
],
vec![
Vec3::new(phi_inv, 0.0, phi),
Vec3::new(-phi_inv, 0.0, phi),
Vec3::new(-1.0, 1.0, 1.0),
Vec3::new(0.0, phi, phi_inv),
Vec3::new(1.0, 1.0, 1.0),
],
vec![
Vec3::new(phi_inv, 0.0, phi),
Vec3::new(1.0, -1.0, 1.0),
Vec3::new(phi, -phi_inv, 0.0),
Vec3::new(phi, phi_inv, 0.0),
Vec3::new(1.0, 1.0, 1.0),
],
vec![
Vec3::new(-phi_inv, 0.0, phi),
Vec3::new(-1.0, -1.0, 1.0),
Vec3::new(-phi, -phi_inv, 0.0),
Vec3::new(-phi, phi_inv, 0.0),
Vec3::new(-1.0, 1.0, 1.0),
],
vec![
Vec3::new(-1.0, -1.0, 1.0),
Vec3::new(-phi, -phi_inv, 0.0),
Vec3::new(-1.0, -1.0, -1.0),
Vec3::new(0.0, -phi, -phi_inv),
Vec3::new(0.0, -phi, phi_inv),
],
vec![
Vec3::new(0.0, -phi, phi_inv),
Vec3::new(0.0, -phi, -phi_inv),
Vec3::new(1.0, -1.0, -1.0),
Vec3::new(phi, -phi_inv, 0.0),
Vec3::new(1.0, -1.0, 1.0),
],
vec![
Vec3::new(0.0, phi, phi_inv),
Vec3::new(0.0, phi, -phi_inv),
Vec3::new(-1.0, 1.0, -1.0),
Vec3::new(-phi, phi_inv, 0.0),
Vec3::new(-1.0, 1.0, 1.0),
],
vec![
Vec3::new(1.0, 1.0, 1.0),
Vec3::new(phi, phi_inv, 0.0),
Vec3::new(1.0, 1.0, -1.0),
Vec3::new(0.0, phi, -phi_inv),
Vec3::new(0.0, phi, phi_inv),
],
vec![
Vec3::new(1.0, -1.0, -1.0),
Vec3::new(0.0, -phi, -phi_inv),
Vec3::new(-1.0, -1.0, -1.0),
Vec3::new(-phi_inv, 0.0, -phi),
Vec3::new(phi_inv, 0.0, -phi),
],
vec![
Vec3::new(1.0, 1.0, -1.0),
Vec3::new(0.0, phi, -phi_inv),
Vec3::new(-1.0, 1.0, -1.0),
Vec3::new(-phi_inv, 0.0, -phi),
Vec3::new(phi_inv, 0.0, -phi),
],
vec![
Vec3::new(1.0, 1.0, -1.0),
Vec3::new(phi, phi_inv, 0.0),
Vec3::new(phi, -phi_inv, 0.0),
Vec3::new(1.0, -1.0, -1.0),
Vec3::new(phi_inv, 0.0, -phi),
],
vec![
Vec3::new(-1.0, 1.0, -1.0),
Vec3::new(-phi, phi_inv, 0.0),
Vec3::new(-phi, -phi_inv, 0.0),
Vec3::new(-1.0, -1.0, -1.0),
Vec3::new(-phi_inv, 0.0, -phi),
],
];
let scale = size * 3f64.sqrt() / 2.0;
faces
.iter()
.flat_map(|face| {
let normal = (face[1] - face[0]).cross(&(face[2] - face[1]));
let transformed_face: Vec<_> = face.iter().map(|v| centre + v * scale).collect();
triangulate_polygon(&transformed_face, &normal, Arc::clone(&material))
})
.collect()
}

155
src/util/tile_iterator.rs
View File

@ -1,155 +0,0 @@
#[derive(Copy, Clone, Debug)]
pub struct Tile {
pub start_column: usize,
pub end_column: usize,
pub start_row: usize,
pub end_row: usize,
}
impl Tile {
pub fn width(&self) -> usize {
self.end_column - self.start_column
}
pub fn height(&self) -> usize {
self.end_row - self.start_row
}
}
#[derive(Clone)]
pub struct TileIterator {
tile_size: usize,
total_height: usize,
total_width: usize,
current_column: usize,
current_row: usize,
}
impl TileIterator {
pub fn new(total_width: usize, total_height: usize, tile_size: usize) -> TileIterator {
// If tile_size*2 is greater than usize::max_value(), increment would overflow
assert!(tile_size > 0 && tile_size * 2 < usize::MAX);
TileIterator {
tile_size,
total_width,
total_height,
current_column: 0,
current_row: 0,
}
}
}
impl Iterator for TileIterator {
type Item = Tile;
fn next(&mut self) -> Option<Tile> {
if self.current_row >= self.total_height {
None
} else {
let start_column = self.current_column;
let end_column = self.total_width.min(start_column + self.tile_size);
let start_row = self.current_row;
let end_row = self.total_height.min(start_row + self.tile_size);
self.current_column += self.tile_size;
if self.current_column >= self.total_width {
self.current_row += self.tile_size;
self.current_column = 0;
}
Some(Tile {
start_column,
end_column,
start_row,
end_row,
})
}
}
}
#[cfg(test)]
mod tests {
use super::*;
use quickcheck::TestResult;
use quickcheck_macros::quickcheck;
#[test]
fn iterator_has_correct_number_of_tiles_when_dimensions_are_multiple_of_tile_size() {
let target = TileIterator::new(20, 15, 5);
assert!(target.count() == 12);
}
#[test]
fn iterator_has_correct_number_of_tiles_when_width_is_one_under_multiple_of_tile_size() {
let target = TileIterator::new(19, 15, 5);
assert!(target.count() == 12);
}
#[test]
fn iterator_has_correct_number_of_tiles_when_width_is_one_over_multiple_of_tile_size() {
let target = TileIterator::new(21, 15, 5);
assert!(target.count() == 15);
}
#[test]
fn iterator_has_correct_number_of_tiles_when_height_is_one_under_multiple_of_tile_size() {
let target = TileIterator::new(20, 14, 5);
assert!(target.count() == 12);
}
#[test]
fn iterator_has_correct_number_of_tiles_when_height_is_one_over_multiple_of_tile_size() {
let target = TileIterator::new(20, 16, 5);
assert!(target.count() == 16);
}
#[quickcheck]
fn tiles_are_expected_size(width: usize, height: usize, tile_size: usize) -> TestResult {
let max_size = 10000;
// Check size of width and height first, since width*height might overflow.
if width > max_size || height > max_size || width * height > max_size {
return TestResult::discard();
}
if tile_size == 0 {
return TestResult::discard();
}
let mut target = TileIterator::new(width, height, tile_size);
TestResult::from_bool(target.all(|tile| {
tile.end_column - tile.start_column <= tile_size
&& tile.end_row - tile.start_row <= tile_size
}))
}
#[quickcheck]
fn iterator_includes_all_coordinates_exactly_once(
width: usize,
height: usize,
tile_size: usize,
) -> TestResult {
let max_size = 10000;
// Check size of width and height first, since width*height might overflow.
if width > max_size || height > max_size || width * height > max_size {
return TestResult::discard();
}
if tile_size == 0 {
return TestResult::discard();
}
let target = TileIterator::new(width, height, tile_size);
let mut index_counts = vec![0; width * height];
let mut total_count = 0;
for tile in target {
for column in tile.start_column..tile.end_column {
for row in tile.start_row..tile.end_row {
index_counts[row * width + column] += 1;
total_count += 1;
if total_count > width * height {
return TestResult::failed();
}
}
}
}
TestResult::from_bool(index_counts.iter().all(|&elem| elem == 1))
}
}

2
test_data/Readme.md
View File

@ -1,2 +0,0 @@
Thanks to the Stanford Computer Graphics Laboratory (https://graphics.stanford.edu/data/3Dscanrep/)
for the famous Stanford Bunny model.

BIN
test_data/stanford_bunny.obj (Stored with Git LFS)
View File

Binary file not shown.