diff --git a/.gitignore b/.gitignore index bc2b56a..b6f4691 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ /target *.lock *.log +*.bak diff --git a/Cargo.toml b/Cargo.toml index ee1a3f7..20457b5 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -5,3 +5,4 @@ edition = "2024" [dependencies] num-traits = "0.2.19" +once_cell = "1.21.3" diff --git a/README.md b/README.md deleted file mode 100644 index c21ec0d..0000000 --- a/README.md +++ /dev/null @@ -1,264 +0,0 @@ - - - -[![Contributors][contributors-shield]][contributors-url] -[![Forks][forks-shield]][forks-url] -[![Stargazers][stars-shield]][stars-url] -[![Issues][issues-shield]][issues-url] -[![Unlicense License][license-shield]][license-url] -[![LinkedIn][linkedin-shield]][linkedin-url] - - - - -
-
- - Logo - - -

Best-README-Template

- -

- An awesome README template to jumpstart your projects! -
- Explore the docs ยป -
-
- View Demo - · - Report Bug - · - Request Feature -

-
- - - - -
- Table of Contents -
    -
  1. - About The Project - -
    2. -
  2. - Getting Started - -
    3. -
  3. Usage
    4. -
  4. Roadmap
    5. -
  5. Contributing
    6. -
  6. License
    7. -
  7. Contact
    8. -
  8. Acknowledgments
    9. -
-
- - - - -## About The Project - -[![Product Name Screen Shot][product-screenshot]](https://example.com) - -There are many great README templates available on GitHub; however, I didn't find one that really suited my needs so I created this enhanced one. I want to create a README template so amazing that it'll be the last one you ever need -- I think this is it. - -Here's why: -* Your time should be focused on creating something amazing. A project that solves a problem and helps others -* You shouldn't be doing the same tasks over and over like creating a README from scratch -* You should implement DRY principles to the rest of your life :smile: - -Of course, no one template will serve all projects since your needs may be different. So I'll be adding more in the near future. You may also suggest changes by forking this repo and creating a pull request or opening an issue. Thanks to all the people have contributed to expanding this template! - -Use the `BLANK_README.md` to get started. - -

(back to top)

- - - -### Built With - -This section should list any major frameworks/libraries used to bootstrap your project. Leave any add-ons/plugins for the acknowledgements section. Here are a few examples. - -* [![Next][Next.js]][Next-url] -* [![React][React.js]][React-url] -* [![Vue][Vue.js]][Vue-url] -* [![Angular][Angular.io]][Angular-url] -* [![Svelte][Svelte.dev]][Svelte-url] -* [![Laravel][Laravel.com]][Laravel-url] -* [![Bootstrap][Bootstrap.com]][Bootstrap-url] -* [![JQuery][JQuery.com]][JQuery-url] - -

(back to top)

- - - - -## Getting Started - -This is an example of how you may give instructions on setting up your project locally. -To get a local copy up and running follow these simple example steps. - -### Prerequisites - -This is an example of how to list things you need to use the software and how to install them. -* npm - ```sh - npm install npm@latest -g - ``` - -### Installation - -_Below is an example of how you can instruct your audience on installing and setting up your app. This template doesn't rely on any external dependencies or services._ - -1. Get a free API Key at [https://example.com](https://example.com) -2. Clone the repo - ```sh - git clone https://github.com/github_username/repo_name.git - ``` -3. Install NPM packages - ```sh - npm install - ``` -4. Enter your API in `config.js` - ```js - const API_KEY = 'ENTER YOUR API'; - ``` -5. Change git remote url to avoid accidental pushes to base project - ```sh - git remote set-url origin github_username/repo_name - git remote -v # confirm the changes - ``` - -

(back to top)

- - - - -## Usage - -Use this space to show useful examples of how a project can be used. Additional screenshots, code examples and demos work well in this space. You may also link to more resources. - -_For more examples, please refer to the [Documentation](https://example.com)_ - -

(back to top)

- - - - -## Roadmap - -- [x] Add Changelog -- [x] Add back to top links -- [ ] Add Additional Templates w/ Examples -- [ ] Add "components" document to easily copy & paste sections of the readme -- [ ] Multi-language Support - - [ ] Chinese - - [ ] Spanish - -See the [open issues](https://github.com/othneildrew/Best-README-Template/issues) for a full list of proposed features (and known issues). - -

(back to top)

- - - - -## Contributing - -Contributions are what make the open source community such an amazing place to learn, inspire, and create. Any contributions you make are **greatly appreciated**. - -If you have a suggestion that would make this better, please fork the repo and create a pull request. You can also simply open an issue with the tag "enhancement". -Don't forget to give the project a star! Thanks again! - -1. Fork the Project -2. Create your Feature Branch (`git checkout -b feature/AmazingFeature`) -3. Commit your Changes (`git commit -m 'Add some AmazingFeature'`) -4. Push to the Branch (`git push origin feature/AmazingFeature`) -5. Open a Pull Request - -### Top contributors: - - - contrib.rocks image - - -

(back to top)

- - - - -## License - -Distributed under the Unlicense License. See `LICENSE.txt` for more information. - -

(back to top)

- - - - -## Contact - -Your Name - [@your_twitter](https://twitter.com/your_username) - email@example.com - -Project Link: [https://github.com/your_username/repo_name](https://github.com/your_username/repo_name) - -

(back to top)

- - - - -## Acknowledgments - -Use this space to list resources you find helpful and would like to give credit to. I've included a few of my favorites to kick things off! - -* [Choose an Open Source License](https://choosealicense.com) -* [GitHub Emoji Cheat Sheet](https://www.webpagefx.com/tools/emoji-cheat-sheet) -* [Malven's Flexbox Cheatsheet](https://flexbox.malven.co/) -* [Malven's Grid Cheatsheet](https://grid.malven.co/) -* [Img Shields](https://shields.io) -* [GitHub Pages](https://pages.github.com) -* [Font Awesome](https://fontawesome.com) -* [React Icons](https://react-icons.github.io/react-icons/search) - -

(back to top)

- - - - - -[contributors-shield]: https://img.shields.io/github/contributors/othneildrew/Best-README-Template.svg?style=for-the-badge -[contributors-url]: https://github.com/othneildrew/Best-README-Template/graphs/contributors -[forks-shield]: https://img.shields.io/github/forks/othneildrew/Best-README-Template.svg?style=for-the-badge -[forks-url]: https://github.com/othneildrew/Best-README-Template/network/members -[stars-shield]: https://img.shields.io/github/stars/othneildrew/Best-README-Template.svg?style=for-the-badge -[stars-url]: https://github.com/othneildrew/Best-README-Template/stargazers -[issues-shield]: https://img.shields.io/github/issues/othneildrew/Best-README-Template.svg?style=for-the-badge -[issues-url]: https://github.com/othneildrew/Best-README-Template/issues -[license-shield]: https://img.shields.io/github/license/othneildrew/Best-README-Template.svg?style=for-the-badge -[license-url]: https://github.com/othneildrew/Best-README-Template/blob/master/LICENSE.txt -[linkedin-shield]: https://img.shields.io/badge/-LinkedIn-black.svg?style=for-the-badge&logo=linkedin&colorB=555 -[linkedin-url]: https://linkedin.com/in/othneildrew -[product-screenshot]: images/screenshot.png -[Next.js]: https://img.shields.io/badge/next.js-000000?style=for-the-badge&logo=nextdotjs&logoColor=white -[Next-url]: https://nextjs.org/ -[React.js]: https://img.shields.io/badge/React-20232A?style=for-the-badge&logo=react&logoColor=61DAFB -[React-url]: https://reactjs.org/ -[Vue.js]: https://img.shields.io/badge/Vue.js-35495E?style=for-the-badge&logo=vuedotjs&logoColor=4FC08D -[Vue-url]: https://vuejs.org/ -[Angular.io]: https://img.shields.io/badge/Angular-DD0031?style=for-the-badge&logo=angular&logoColor=white -[Angular-url]: https://angular.io/ -[Svelte.dev]: https://img.shields.io/badge/Svelte-4A4A55?style=for-the-badge&logo=svelte&logoColor=FF3E00 -[Svelte-url]: https://svelte.dev/ -[Laravel.com]: https://img.shields.io/badge/Laravel-FF2D20?style=for-the-badge&logo=laravel&logoColor=white -[Laravel-url]: https://laravel.com -[Bootstrap.com]: https://img.shields.io/badge/Bootstrap-563D7C?style=for-the-badge&logo=bootstrap&logoColor=white -[Bootstrap-url]: https://getbootstrap.com -[JQuery.com]: https://img.shields.io/badge/jQuery-0769AD?style=for-the-badge&logo=jquery&logoColor=white -[JQuery-url]: https://jquery.com diff --git a/src/core/camera.rs b/src/core/camera.rs new file mode 100644 index 0000000..eddd38c --- /dev/null +++ b/src/core/camera.rs @@ -0,0 +1,76 @@ +use crate::core::film::Film; +use crate::core::geometry::{Point2f, Point3f, Ray, RayDifferential}; +use crate::core::medium::Medium; +use crate::core::pbrt::{Float, RenderingCoordinateSystem}; +use crate::core::spectrum::SampledSpectrum; +use crate::core::transform::{Transform, AnimatedTransform}; + +pub enum Camera { + Perspective(PerspectiveCamera), + Ortographic(OrtographicCamera), + Spherical(SphericalCamera), + Realistic(RealisticCamera), +} + +impl Camera { + pub fn create(&self, name: str, medium: Medium, camera_transform: CameraTransform, film: Film) -> Float { + match self { + Camera::Perspective(s) => s.create(name, medium, camera_transform, film), + Camera::Ortographic(s) => s.create(name, medium, camera_transform, film), + Camera::Spherical(s) => s.create(name, medium, camera_transform, film), + Camera::Realistic(s) => s.create(name, medium, camera_transform, film), + } + } +} + +pub struct CameraSample { + p_film: Point2f, + p_lens: Point2f, + time: Float, + filter_weight: Float, +} + +pub struct CameraRay { + ray: Ray, + weight: SampledSpectrum, +} + +pub struct CameraDifferential { + ray: RayDifferential, + weight: SampledSpectrum, +} + +pub struct CameraTransform { + render_from_camera: AnimatedTransform, + world_from_render: Transform, +} + +impl CameraTransform { + pub fn from_world(world_from_camera: AnimatedTransform, rendering_space: RenderingCoordinateSystem) -> Self { + let world_from_render = match rendering_space { + RenderingCoordinateSystem::Camera => { + let t_mid = (world_from_camera.start_time + world_from_camera.end_time) / 2.; + world_from_camera.interpolate(t_mid); + } + RenderingCoordinateSystem::CameraWorld => { + let t_mid = (world_from_camera.start_time + world_from_camera.end_time) / 2.; + let p_camera = world_from_camera.apply_point(Point3f::default(), t_mid); + Transform::translate(p_camera.to_vec()); + } + RenderingCoordinateSystem::World => Transform::identity(), + }; + } + + pub fn render_from_camera(&self, p: Point3f, time: Float) -> Point3f { + self.render_from_camera.apply_point(p, time) + } + + pub fn camera_from_render(&self, p: Point3f, time: Float) -> Point3f { + self.render_from_camera.apply_inverse_point(p, time) + } +} + +pub struct PerspectiveCamera; +pub struct OrtographicCamera; +pub struct SphericalCamera; +pub struct RealisticCamera; diff --git a/src/core/color.rs b/src/core/color.rs new file mode 100644 index 0000000..d2f69ac --- /dev/null +++ b/src/core/color.rs @@ -0,0 +1,496 @@ +use std::ops::{Index, IndexMut, Add, AddAssign, Sub, SubAssign, Mul, MulAssign, Div, DivAssign, Neg}; +use std::fmt; + +use crate::core::pbrt::{Float, lerp}; +use crate::core::geometry::Point2f; +use crate::core::spectrum::Spectrum; + +#[derive(Debug, Clone)] +pub struct XYZ { + x: Float, + y: Float, + z: Float, +} + +impl XYZ { + pub fn new(x: Float, y: Float, z: Float) -> Self { + Self { x, y, z } + } + + pub fn x(&self) -> Float { + self.x + } + + pub fn y(&self) -> Float { + self.y + } + + pub fn z(&self) -> Float { + self.z + } + + pub fn average(&self) -> Float { + (self.x + self.y + self.z) / 3.0 + } + + pub fn xy(&self) -> Point2f { + let sum = self.x + self.y + self.z; + Point2f::new(self.x / sum, self.y / sum) + } + + pub fn from_xyy(xy: Point2f, y: Float) -> Self { + if xy.y() == 0.0 { + return Self { x: 0.0, y: 0.0, z: 0.0} + } + Self::new(xy.x() * y / xy.y(), y, (1.0 - xy.x() - xy.y()) * y / xy.y()) + } +} + +impl Index for XYZ { + type Output = Float; + fn index(&self, index: usize) -> &Self::Output { + debug_assert!(index < 3); + match index { + 0 => &self.x, + 1 => &self.y, + _ => &self.z, + } + } +} + +impl IndexMut for XYZ { + fn index_mut(&mut self, index: usize) -> &mut Self::Output { + debug_assert!(index < 3); + match index { + 0 => &mut self.x, + 1 => &mut self.y, + _ => &mut self.z, + } + } +} + +impl Neg for XYZ { + type Output = Self; + fn neg(self) -> Self { + Self { + x: -self.x, + y: -self.y, + z: -self.z, + } + } +} + +impl Add for XYZ { + type Output = Self; + fn add(self, rhs: Self) -> Self { + Self { + x: self.x + rhs.x, + y: self.y + rhs.y, + z: self.z + rhs.z, + } + } +} + +impl AddAssign for XYZ { + fn add_assign(&mut self, rhs: Self) { + self.x += rhs.x; + self.y += rhs.y; + self.z += rhs.z; + } +} + +impl Sub for XYZ { + type Output = Self; + fn sub(self, rhs: Self) -> Self { + Self { + x: self.x - rhs.x, + y: self.y - rhs.y, + z: self.z - rhs.z, + } + } +} + +impl SubAssign for XYZ { + fn sub_assign(&mut self, rhs: Self) { + self.x -= rhs.x; + self.y -= rhs.y; + self.z -= rhs.z; + } +} + +impl Sub for Float { + type Output = XYZ; + fn sub(self, rhs: XYZ) -> XYZ { + XYZ::new(self - rhs.x, self - rhs.y, self - rhs.z) + } +} + +impl Mul for XYZ { + type Output = Self; + fn mul(self, rhs: Self) -> Self { + Self { + x: self.x * rhs.x, + y: self.y * rhs.y, + z: self.z * rhs.z, + } + } +} + +impl MulAssign for XYZ { + fn mul_assign(&mut self, rhs: Self) { + self.x *= rhs.x; + self.y *= rhs.y; + self.z *= rhs.z; + } +} + +// Scalar multiplication (xyz * float) +impl Mul for XYZ { + type Output = Self; + fn mul(self, rhs: Float) -> Self { + Self { + x: self.x * rhs, + y: self.y * rhs, + z: self.z * rhs, + } + } +} + +impl MulAssign for XYZ { + fn mul_assign(&mut self, rhs: Float) { + self.x *= rhs; + self.y *= rhs; + self.z *= rhs; + } +} + +impl Mul for Float { + type Output = XYZ; + fn mul(self, rhs: XYZ) -> XYZ { + rhs * self + } +} + +impl Div for XYZ { + type Output = Self; + fn div(self, rhs: Self) -> Self { + Self { + x: self.x / rhs.x, + y: self.y / rhs.y, + z: self.z / rhs.z, + } + } +} + +impl DivAssign for XYZ { + fn div_assign(&mut self, rhs: Self) { + self.x /= rhs.x; + self.y /= rhs.y; + self.z /= rhs.z; + } +} + +impl Div for XYZ { + type Output = Self; + fn div(self, rhs: Float) -> Self { + debug_assert_ne!(rhs, 0.0, "Division by zero."); + let inv = 1.0 / rhs; + self * inv + } +} + +impl DivAssign for XYZ { + fn div_assign(&mut self, rhs: Float) { + debug_assert_ne!(rhs, 0.0, "Division by zero."); + let inv = 1.0 / rhs; + *self *= inv; + } +} + +impl fmt::Display for XYZ { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "[ {}, {}, {} ]", self.x, self.y, self.z) + } +} + + +#[derive(Debug, Clone)] +pub struct RGB { + pub r: Float, + pub g: Float, + pub b: Float, +} + + +impl RGB { + pub fn new(r: Float, g: Float, b: Float) -> Self { + Self { r, g, b } + } + + pub fn average(&self) -> Float { + (self.r + self.g + self.b) / 3.0 + } +} + +impl Index for RGB { + type Output = Float; + fn index(&self, index: usize) -> &Self::Output { + debug_assert!(index < 3); + match index { + 0 => &self.r, + 1 => &self.g, + _ => &self.b, + } + } +} + +impl IndexMut for RGB { + fn index_mut(&mut self, index: usize) -> &mut Self::Output { + debug_assert!(index < 3); + match index { + 0 => &mut self.r, + 1 => &mut self.g, + _ => &mut self.b, + } + } +} + +impl Neg for RGB { + type Output = Self; + fn neg(self) -> Self { + Self { + r: -self.r, + g: -self.g, + b: -self.b, + } + } +} + +impl Add for RGB { + type Output = Self; + fn add(self, rhs: Self) -> Self { + Self { + r: self.r + rhs.r, + g: self.g + rhs.g, + b: self.b + rhs.b, + } + } +} + +impl AddAssign for RGB { + fn add_assign(&mut self, rhs: Self) { + self.r += rhs.r; + self.g += rhs.g; + self.b += rhs.b; + } +} + +impl Sub for RGB { + type Output = Self; + fn sub(self, rhs: Self) -> Self { + Self { + r: self.r - rhs.r, + g: self.g - rhs.g, + b: self.b - rhs.b, + } + } +} + +impl SubAssign for RGB { + fn sub_assign(&mut self, rhs: Self) { + self.r -= rhs.r; + self.g -= rhs.g; + self.b -= rhs.b; + } +} + +impl Sub for Float { + type Output = RGB; + fn sub(self, rhs: RGB) -> RGB { + RGB::new(self - rhs.r, self - rhs.g, self - rhs.b) + } +} + +impl Mul for RGB { + type Output = Self; + fn mul(self, rhs: Self) -> Self { + Self { + r: self.r * rhs.r, + g: self.g * rhs.g, + b: self.b * rhs.b, + } + } +} + +impl MulAssign for RGB { + fn mul_assign(&mut self, rhs: Self) { + self.r *= rhs.r; + self.g *= rhs.g; + self.b *= rhs.b; + } +} + +// Scalar multiplication (ryz * float) +impl Mul for RGB { + type Output = Self; + fn mul(self, rhs: Float) -> Self { + Self { + r: self.r * rhs, + g: self.g * rhs, + b: self.b * rhs, + } + } +} + +impl MulAssign for RGB { + fn mul_assign(&mut self, rhs: Float) { + self.r *= rhs; + self.g *= rhs; + self.b *= rhs; + } +} + +impl Mul for Float { + type Output = RGB; + fn mul(self, rhs: RGB) -> RGB { + rhs * self + } +} + +impl Div for RGB { + type Output = Self; + fn div(self, rhs: Self) -> Self { + Self { + r: self.r / rhs.r, + g: self.g / rhs.g, + b: self.b / rhs.b, + } + } +} + +impl DivAssign for RGB { + fn div_assign(&mut self, rhs: Self) { + self.r /= rhs.r; + self.g /= rhs.g; + self.b /= rhs.b; + } +} + +impl Div for RGB { + type Output = Self; + fn div(self, rhs: Float) -> Self { + debug_assert_ne!(rhs, 0.0, "Division by zero."); + let inv = 1.0 / rhs; + self * inv + } +} + +impl DivAssign for RGB { + fn div_assign(&mut self, rhs: Float) { + debug_assert_ne!(rhs, 0.0, "Division by zero."); + let inv = 1.0 / rhs; + *self *= inv; + } +} + +impl fmt::Display for RGB { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "[ {}, {}, {} ]", self.r, self.g, self.b) + } +} + + +pub const RES: usize = 64; +pub type CoefficientArray = [[[[[Float; 3]; RES]; RES]; RES]; 3]; + +pub struct RGBToSpectrumTable { + z_nodes: &'static [f32], + coeffs: &'static CoefficientArray, +} + +#[derive(Debug, Default, Copy, Clone)] +pub struct RGBSigmoidPolynomial { + c0: Float, + c1: Float, + c2: Float, +} + +impl RGBSigmoidPolynomial { + pub fn new(c0: Float, c1: Float, c2: Float) -> Self { + Self { c0, c1, c2 } + } + + pub fn evaluate(&self, lambda: Float) -> Float { + let eval = match crate::core::pbrt::evaluate_polynomial(lambda, &[self.c0, self.c1, self.c2]) { + Some(value) => value, + None => { panic!("evaluate_polynomial returned None with non-empty coefficients, this should not happen.") }, + }; + + Self::s(eval) + } + + pub fn max_value(&self) -> Float { + let lambda = -self.c1 / (2.0 * self.c0); + let result = self.evaluate(360.0).max(self.evaluate(830.0)); + if lambda >= 360.0 && lambda <= 830.0 { + return result.max(self.evaluate(lambda)) + } + result + } + + fn s(x: Float) -> Float { + if x.is_infinite() { + if x > 0.0 { + return 1.0 + } else { + return 0.0 + } + } + 0.5 + x / (2.0 * (1.0 + (x * x)).sqrt()) + } +} + +impl RGBToSpectrumTable { + pub fn new(z_nodes: &'static [f32], coeffs: &'static CoefficientArray) -> Self { + Self { z_nodes, coeffs } + } + + pub fn to_polynomial(&self, rgb: RGB) -> RGBSigmoidPolynomial { + if rgb[0] == rgb[1] && rgb[1] == rgb[2] { + return RGBSigmoidPolynomial::new(0.0, 0.0, (rgb[0] - 0.5)/(rgb[0] * (1.0 - rgb[0])).sqrt()) + } + let maxc; + if rgb[0] > rgb[1] { + if rgb[0] > rgb[2] { + maxc = 0; + } else { + maxc = 2; + } + } else { + if rgb[1] > rgb[2] { + maxc = 1; + } else { + maxc = 2; + } + } + + let z = rgb[maxc]; + let x = rgb[(maxc + 1) % 3] * (RES - 1) as Float / z; + let y = rgb[(maxc + 2) % 3] * (RES - 1) as Float / z; + + let xi = x.min(RES as Float - 2.0); + let yi = y.min(RES as Float - 2.0); + let zi = crate::core::pbrt::find_interval(RES, |i: usize| self.z_nodes[i] < z); + let dx = (x - xi) as usize; + let dy = (y - yi) as usize; + let dz = (z - self.z_nodes[zi]) / (self.z_nodes[zi + 1] - self.z_nodes[zi]); + let mut c = [0.0; 3]; + for i in 0..3 { + let co = |dx: usize, dy: usize, dz: usize| self.coeffs[maxc][zi as usize + dz][yi as usize + dy][xi as usize + dx][i]; + c[i] = lerp(dz, + lerp(dy as Float, lerp(dx as Float, co(0, 0, 0) as Float, co(1, 0, 0)) as Float, lerp(dx as Float, co(0, 1, 0) as Float, co(1, 1, 0) as Float)), + lerp(dy as Float, lerp(dx as Float, co(0, 0, 1) as Float, co(1, 0, 1)) as Float, lerp(dx as Float, co(0, 1, 1) as Float, co(1, 1, 1) as Float))); + } + RGBSigmoidPolynomial { c0: c[0], c1: c[1], c2: c[2] } + } +} diff --git a/src/core/colorspace.rs b/src/core/colorspace.rs new file mode 100644 index 0000000..05e095e --- /dev/null +++ b/src/core/colorspace.rs @@ -0,0 +1,75 @@ +use crate::core::pbrt::Float; +use crate::core::geometry::Point2f; +use crate::core::transform::SquareMatrix; +use crate::core::color::{RGBSigmoidPolynomial, RGBToSpectrumTable, RGB, XYZ}; +use crate::core::spectrum::{DenselySampledSpectrum, Spectrum, SampledSpectrum}; + +use std::cmp::{Eq, PartialEq}; +use std::sync::Arc; + +#[derive(Clone)] +pub struct RGBColorspace { + r: Point2f, + g: Point2f, + b: Point2f, + w: Point2f, + illuminant: Spectrum, + rgb_to_spectrum_table: Arc, + xyz_from_rgb: SquareMatrix, + rgb_from_xyz: SquareMatrix, + +} + +impl RGBColorspace { + pub fn new(r: Point2f, + g: Point2f, + b: Point2f, + illuminant: Spectrum, + rgb_to_spectrum_table: RGBToSpectrumTable) -> Self { + let w_xyz = illuminant.to_xyz(); + let w = w_xyz.xy(); + let r_xyz = XYZ::from_xyy(r, 1.0); + let g_xyz = XYZ::from_xyy(g, 1.0); + let b_xyz = XYZ::from_xyy(b, 1.0); + let rgb_values = [[r_xyz.x(), g_xyz.x(), b_xyz.x()], [r_xyz.y(), g_xyz.y(), b_xyz.y()], [r_xyz.z(), g_xyz.z(), g_xyz.z()]]; + let rgb: SquareMatrix = SquareMatrix { m: rgb_values }; + let c = match rgb.inverse() { + Some(inv_matrix) => { inv_matrix * w_xyz }, + None => { panic!("Cannot create RGBColorspace: The RGB primaries form a singular matrix."); } + }; + let xyz_from_rgb_m = [[c[0], 0.0, 0.0], [0.0, c[1], 0.0], [0.0, 0.0, c[2]]]; + let xyz_from_rgb = rgb * SquareMatrix { m: xyz_from_rgb_m }; + let rgb_from_xyz = xyz_from_rgb.inverse().expect("Failed to invert the XYZfromRGB matrix. Is it singular?"); + + Self { r, g, b, w, illuminant, rgb_to_spectrum_table: Arc::new(rgb_to_spectrum_table), xyz_from_rgb, rgb_from_xyz } + } + + pub fn to_xyz(&self, rgb: RGB) -> XYZ { + self.xyz_from_rgb.transform_to_xyz(rgb) + } + + pub fn to_rgb(&self, xyz: XYZ) -> RGB { + self.rgb_from_xyz.transform_to_rgb(xyz) + } + + pub fn to_rgb_coeffs(&self, rgb: RGB) -> RGBSigmoidPolynomial { + self.rgb_to_spectrum_table.to_polynomial(rgb) + + } + + pub fn convert_colorspace(&self, other: &RGBColorspace) -> SquareMatrix { + if self == other { + return SquareMatrix::default() + } + + self.rgb_from_xyz * other.xyz_from_rgb + } +} + +impl PartialEq for RGBColorspace { + fn eq(&self, other: &Self) -> bool { + self.r == other.r && self.g == other.g && self.b == other.b && self.w == other.w && + Arc::ptr_eq(&self.rgb_to_spectrum_table, &other.rgb_to_spectrum_table) + } +} + diff --git a/src/core/film.rs b/src/core/film.rs new file mode 100644 index 0000000..49a0215 --- /dev/null +++ b/src/core/film.rs @@ -0,0 +1,9 @@ +pub enum Film { + RGB(RGBFilm), + GBuffer(GBufferBFilm), + Spectral(SpectralFilm), +} + +pub struct RGBFilm; +pub struct GBufferBFilm; +pub struct SpectralFilm; diff --git a/src/core/geometry.rs b/src/core/geometry.rs index 9ceba3b..35ce9cf 100644 --- a/src/core/geometry.rs +++ b/src/core/geometry.rs @@ -1,10 +1,11 @@ -use num_traits::{Num, Bounded, Float as NumFloat}; +use num_traits::{Num, Bounded, Float as NumFloat, Zero, Signed}; use std::ops::{Sub, SubAssign, Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Index, IndexMut}; use std::sync::Arc; use std::f32::consts::PI; -use crate::core::pbrt::Float; +use crate::core::pbrt::{Float, next_float_up, next_float_down}; use crate::core::pbrt; +use crate::core::interval::Interval; use crate::core::medium::Medium; pub trait Tuple: @@ -37,9 +38,11 @@ macro_rules! impl_tuple_core { } impl Default for $Struct { - fn default() -> Self { - Self([T::default(); N]) - } + fn default() -> Self { Self([T::default(); N]) } + } + + impl $Struct where T: Zero + Copy { + #[inline] pub fn zero() -> Self { Self([T::zero(); N]) } } impl Index for $Struct { @@ -58,7 +61,11 @@ macro_rules! impl_tuple_core { Self(result) } } + }; +} +macro_rules! impl_scalar_ops { + ($Struct:ident) => { impl Mul for $Struct where T: Mul + Copy { type Output = Self; fn mul(self, rhs: T) -> Self::Output { @@ -67,6 +74,14 @@ macro_rules! impl_tuple_core { Self(result) } } + + impl Mul<$Struct> for Float { + type Output = $Struct; + fn mul(self, rhs: $Struct) -> Self::Output { + rhs * self + } + } + impl MulAssign for $Struct where T: MulAssign + Copy { fn mul_assign(&mut self, rhs: T) { for i in 0..N { self.0[i] *= rhs; } @@ -93,39 +108,89 @@ impl_tuple_core!(Vector); impl_tuple_core!(Point); impl_tuple_core!(Normal); -macro_rules! impl_tuple_ops { - ($Struct:ident) => { - impl Add for $Struct where T: Add + Copy { - type Output = Self; - fn add(self, rhs: Self) -> Self::Output { +impl_scalar_ops!(Vector); +impl_scalar_ops!(Normal); + + +macro_rules! impl_op { + ($Op:ident, $op:ident, $Lhs:ident, $Rhs:ident, $Output:ident) => { + impl $Op<$Rhs> for $Lhs + where + T: $Op + Copy, + { + type Output = $Output; + fn $op(self, rhs: $Rhs) -> Self::Output { let mut result = self.0; - for i in 0..N { result[i] = self.0[i] + rhs.0[i]; } - Self(result) - } - } - impl AddAssign for $Struct where T: AddAssign + Copy { - fn add_assign(&mut self, rhs: Self) { - for i in 0..N { self.0[i] += rhs.0[i]; } - } - } - impl Sub for $Struct where T: Sub + Copy { - type Output = Self; - fn sub(self, rhs: Self) -> Self::Output { - let mut result = self.0; - for i in 0..N { result[i] = self.0[i] - rhs.0[i]; } - Self(result) - } - } - impl SubAssign for $Struct where T: SubAssign + Copy { - fn sub_assign(&mut self, rhs: Self) { - for i in 0..N { self.0[i] -= rhs.0[i]; } + for i in 0..N { + result[i] = self.0[i].$op(rhs.0[i]); + } + $Output(result) } } }; } -impl_tuple_ops!(Vector); -impl_tuple_ops!(Normal); +macro_rules! impl_op_assign { + ($OpAssign:ident, $op_assign:ident, $Lhs:ident, $Rhs:ident) => { + impl $OpAssign<$Rhs> for $Lhs + where + T: $OpAssign + Copy, + { + fn $op_assign(&mut self, rhs: $Rhs) { + for i in 0..N { + self.0[i].$op_assign(rhs.0[i]); + } + } + } + }; +} + +// Addition +impl_op!(Add, add, Vector, Vector, Vector); +impl_op!(Add, add, Point, Vector, Point); +impl_op!(Add, add, Vector, Point, Point); +impl_op!(Add, add, Normal, Normal, Normal); + +// Subtraction +impl_op!(Sub, sub, Vector, Vector, Vector); +impl_op!(Sub, sub, Point, Vector, Point); +impl_op!(Sub, sub, Point, Point, Vector); +impl_op!(Sub, sub, Normal, Normal, Normal); + +// AddAssign +impl_op_assign!(AddAssign, add_assign, Vector, Vector); +impl_op_assign!(AddAssign, add_assign, Point, Vector); +impl_op_assign!(AddAssign, add_assign, Normal, Normal); + +// SubAssign +impl_op_assign!(SubAssign, sub_assign, Vector, Vector); +impl_op_assign!(SubAssign, sub_assign, Point, Vector); +impl_op_assign!(SubAssign, sub_assign, Normal, Normal); + +pub trait Dot { + type Output; + fn dot(self, rhs: Rhs) -> Self::Output; +} + +macro_rules! impl_dot_for { + ($Lhs:ident, $Rhs:ident) => { + impl Dot<$Rhs> for $Lhs { + type Output = T; + fn dot(self, rhs: $Rhs) -> T { + let mut sum = T::zero(); + for i in 0..N { + sum = sum + self[i] * rhs[i]; + } + sum + } + } + }; +} + +impl_dot_for!(Normal, Vector); +impl_dot_for!(Vector, Normal); +impl_dot_for!(Vector, Vector); +impl_dot_for!(Normal, Normal); impl From> for Normal { fn from(v: Vector) -> Self { Self(v.0) } @@ -134,45 +199,6 @@ impl From> for Vector { fn from(n: Normal) -> Self { Self(n.0) } } -impl Sub for Point where T: Sub + Copy { - type Output = Vector; - fn sub(self, rhs: Self) -> Self::Output { - let mut result = self.0; - for i in 0..N { result[i] = self.0[i] - rhs.0[i]; } - Vector(result) - } -} - -// Point + Vector -> Point -impl Add> for Point where T: Add + Copy { - type Output = Self; - fn add(self, rhs: Vector) -> Self::Output { - let mut result = self.0; - for i in 0..N { result[i] = self.0[i] + rhs.0[i]; } - Self(result) - } -} -impl AddAssign> for Point where T: AddAssign + Copy { - fn add_assign(&mut self, rhs: Vector) { - for i in 0..N { self.0[i] += rhs.0[i]; } - } -} - -// Point - Vector -> Point -impl Sub> for Point where T: Sub + Copy { - type Output = Self; - fn sub(self, rhs: Vector) -> Self::Output { - let mut result = self.0; - for i in 0..N { result[i] = self.0[i] - rhs.0[i]; } - Self(result) - } -} -impl SubAssign> for Point where T: SubAssign + Copy { - fn sub_assign(&mut self, rhs: Vector) { - for i in 0..N { self.0[i] -= rhs.0[i]; } - } -} - macro_rules! impl_accessors { ($Struct:ident) => { impl $Struct { @@ -193,11 +219,6 @@ impl_accessors!(Normal); // Vector stuff -pub trait Dot { - type Output; - fn dot(self, rhs: Rhs) -> Self::Output; -} - pub trait Normed { type Scalar; fn norm_squared(&self) -> Self::Scalar; @@ -205,16 +226,8 @@ pub trait Normed { fn normalize(self) -> Self where Self: Sized, Self: Div, Self::Scalar: NumFloat; } -macro_rules! impl_vector_math_for { +macro_rules! impl_normed_for { ($Struct:ident) => { - impl Dot for $Struct { - type Output = T; - fn dot(self, rhs: Self) -> T { - let mut sum = T::zero(); - for i in 0..N { sum = sum + self[i] * rhs[i]; } - sum - } - } impl Normed for $Struct { type Scalar = T; fn norm_squared(&self) -> T { self.dot(*self) } @@ -223,14 +236,34 @@ macro_rules! impl_vector_math_for { }; } -impl_vector_math_for!(Vector); -impl_vector_math_for!(Normal); +impl_normed_for!(Vector); +impl_normed_for!(Normal); + +macro_rules! impl_abs { + ($Struct:ident) => { + impl $Struct + where + T: Signed + Copy, + { + pub fn abs(self) -> Self { + let mut result = self.0; + for i in 0..N { + result[i] = result[i].abs(); + } + Self(result) + } + } + }; +} + +impl_abs!(Vector); +impl_abs!(Normal); impl Point where T: NumFloat + Copy, - Point: Sub>, // Point - Point -> Vector - Vector: Normed, // Vector has a norm + Point: Sub>, + Vector: Normed, { pub fn distance(self, other: Self) -> T { (self - other).norm() @@ -241,19 +274,25 @@ where } } - // Utility aliases and functions pub type Point2 = Point; +pub type Point2f = Point2; +pub type Point2i = Point2; pub type Point3 = Point; pub type Point3f = Point3; pub type Point3i = Point3; pub type Vector2 = Vector; +pub type Vector2f = Vector2; +pub type Vector2i = Vector2; pub type Vector3 = Vector; pub type Vector3f = Vector3; pub type Vector3i = Vector3; pub type Normal3 = Normal; pub type Normal3f = Normal3; pub type Normal3i = Normal3; +pub type Vector3fi = Vector; +pub type Point3fi = Point; + impl Vector2 { pub fn new(x: T, y: T) -> Self { Self([x, y]) } } impl Point2 { pub fn new(x: T, y: T) -> Self { Self([x, y]) } } @@ -288,11 +327,97 @@ where T: Num + NumFloat + Copy + Neg } } +impl Point { + pub fn new_from_point(p: Point) -> Self { + let mut arr = [Interval::default(); N]; + for i in 0..N { + arr[i] = Interval::new(p[i]); + } + Self(arr) + } + + pub fn new_with_error(p: Point, e: Vector) -> Self { + let mut arr = [Interval::default(); N]; + for i in 0..N { + arr[i] = Interval::new_from_value_and_error(p[i], e[i]); + } + Self(arr) + } + + pub fn error(&self) -> Vector { + let mut arr = [0.0; N]; + for i in 0..N { + arr[i] = self[i].width() / 2.0; + } + Vector(arr) + } + + pub fn midpoint(&self) -> Point { + let mut arr = [0.0; N]; + for i in 0..N { + arr[i] = self[i].midpoint(); + } + Point(arr) + } + + pub fn is_exact(&self) -> bool { + self.0.iter().all(|interval| interval.width() == 0.0) + } +} + +impl Vector { + pub fn new_from_vector(v: Vector) -> Self { + let mut arr = [Interval::default(); N]; + for i in 0..N { + arr[i] = Interval::new(v[i]); + } + Self(arr) + } + + pub fn new_with_error(v: Vector, e: Vector) -> Self { + let mut arr = [Interval::default(); N]; + for i in 0..N { + arr[i] = Interval::new_from_value_and_error(v[i], e[i]); + } + Self(arr) + } + + pub fn error(&self) -> Vector { + let mut arr = [0.0; N]; + for i in 0..N { + arr[i] = self[i].width() / 2.0; + } + Vector(arr) + } + + pub fn midpoint(&self) -> Vector { + let mut arr = [0.0; N]; + for i in 0..N { + arr[i] = self[i].midpoint(); + } + Vector(arr) + } + + pub fn is_exact(&self) -> bool { + self.0.iter().all(|interval| interval.width() == 0.0) + } +} + +impl From> for Point { + fn from(pi: Point) -> Self { + let mut arr = [0.0; N]; + for i in 0..N { + arr[i] = pi[i].midpoint(); + } + Point(arr) + } +} + impl Normal3 where T: Num + PartialOrd + Copy + Neg { pub fn face_forward(self, v: Vector3) -> Self { - if self.dot(v.into()) < T::zero() { -self } else { self } + if self.dot(v) < T::zero() { -self } else { self } } } @@ -318,7 +443,7 @@ pub fn spherical_phi(v: Vector3f) -> Float { } } -// AABB bouding boxes +// AABB BOUNDING BOXES #[derive(Debug, Copy, Clone, PartialEq)] pub struct Bounds { pub p_min: Point, @@ -469,7 +594,7 @@ where < as Sub>::Output as Normed>::Scalar: NumFloat, { let two = T::one() + T::one(); - let center = (self.p_min + self.diagonal()) / two; + let center = self.p_min + self.diagonal() / two; let radius = if self.contains(center) { center.distance(self.p_max) } else { T::zero() }; (center, radius) } @@ -483,6 +608,16 @@ pub struct Frame { } impl Frame { + pub fn from_x(x: Vector3f) -> Self { + let (y, z) = x.normalize().coordinate_system(); + Self { x: x.normalize(), y, z } + } + + pub fn from_y(y: Vector3f) -> Self { + let (z, x) = y.normalize().coordinate_system(); + Self { x, y: y.normalize(), z } + } + pub fn from_z(z: Vector3f) -> Self { let (x, y) = z.normalize().coordinate_system(); Self { x, y, z: z.normalize() } @@ -512,6 +647,59 @@ impl Ray { self.o + self.d * t } + pub fn offset_origin(p: &Point3fi, n: &Normal3f, w: &Vector3f) -> Point3f { + let d: Float = n.abs().dot(p.error()); + let normal: Vector3f = Vector3f::from(*n); + + let mut offset = p.midpoint(); + if w.dot(normal) < 0.0 { + offset -= normal * d; + } else { + offset += normal * d; + } + + for i in 0..3 { + if n[i] > 0.0 { + offset[i] = next_float_up(offset[i]); + } else if n[i] < 0.0 { + offset[i] = next_float_down(offset[i]); + } + } + offset + } + + pub fn spawn(pi: &Point3fi, n: &Normal3f, time: Float, d: Vector3f) -> Ray { + let origin = Self::offset_origin(pi, n, &d); + Ray { + o: origin, + d, + time, + medium: None, + differential: None, + } + } + + pub fn spawn_to_point(p_from: &Point3fi, n: &Normal3f, time: Float, p_to: Point3f) -> Ray { + let d = p_to - p_from.midpoint(); + Self::spawn(p_from, n, time, d) + } + + pub fn spawn_to_interaction(p_from: &Point3fi, n_from: &Normal3f, time: Float, p_to: &Point3fi, n_to: &Normal3f) -> Ray { + let dir_for_offset = p_to.midpoint() - p_from.midpoint(); + let pf = Self::offset_origin(p_from, n_from, &dir_for_offset); + let pt = Self::offset_origin(p_to, n_to, &(pf - p_to.midpoint())); + + let d = pt - pf; + + Ray { + o: pf, + d, + time, + medium: None, + differential: None, + } + } + pub fn scale_differentials(&mut self, s: Float) { if let Some(differential) = &mut self.differential { differential.rx_origin = self.o + (differential.rx_origin - self.o) * s; diff --git a/src/core/interaction.rs b/src/core/interaction.rs index e69de29..dbeb73c 100644 --- a/src/core/interaction.rs +++ b/src/core/interaction.rs @@ -0,0 +1,195 @@ +use crate::core::geometry::{Vector3f, Normal3f, Point2f, Point3f, Point3fi, Normed, Ray, Dot}; +use crate::core::pbrt::Float; +use crate::core::medium::{Medium, MediumInterface}; +use crate::core::light::Light; + +use std::any::Any; +use std::sync::Arc; + +pub struct InteractionData { + pub pi: Point3fi, + pub time: Float, + pub wo: Vector3f, // Outgoing direction +} + +pub trait Interaction: Send + Sync { + fn get_common(&self) -> &InteractionData; + fn as_any(&self) -> &dyn Any; + + fn p(&self) -> Point3f { self.get_common().pi.into() } + fn pi(&self) -> Point3fi { self.get_common().pi } + fn time(&self) -> Float { self.get_common().time } + fn wo(&self) -> Vector3f { self.get_common().wo } + + fn is_surface_interaction(&self) -> bool { + self.as_any().is::() + } + + fn is_medium_interaction(&self) -> bool { + self.as_any().is::() + } + + /// Determines the medium for a ray starting at this interaction. + fn get_medium(&self, w: Vector3f) -> Option>; + + /// Spawns a new ray starting from this interaction point. + fn spawn_ray(&self, d: Vector3f) -> Ray; + + /// Spawns a ray from this interaction to a specified point. + fn spawn_ray_to_point(&self, p2: Point3f) -> Ray { + let origin = self.p(); + let direction = p2 - origin; + Ray { + o: origin, + d: direction, + time: self.time(), + medium: self.get_medium(direction), + differential: None, + } + } + + /// Spawns a ray from this interaction to another one. + /// The default implementation simply targets the other interaction's point. + fn spawn_ray_to_interaction(&self, other: &dyn Interaction) -> Ray { + self.spawn_ray_to_point(other.p()) + } +} + +pub struct ShadingGeometry { + pub n: Normal3f, + pub dpdu: Vector3f, + pub dpdv: Vector3f, + pub dndu: Normal3f, + pub dndv: Normal3f, +} + +pub struct SurfaceInteraction { + pub common: InteractionData, + pub n: Normal3f, + pub uv: Point2f, + pub dpdu: Vector3f, + pub dpdv: Vector3f, + pub dndu: Normal3f, + pub dndv: Normal3f, + pub shading: ShadingGeometry, + pub medium_interface: Option>, + pub face_index: bool, + pub area_light: Option>, + pub dpdx: Vector3f, + pub dpdy: Vector3f, + pub dudx: Float, + pub dvdx: Float, + pub dudy: Float, + pub dvdy: Float, + // pub shape: Option>, + // pub bsdf: Option, +} + +pub struct PhaseFunction; + +pub struct MediumInteraction { + pub common: InteractionData, + pub medium: Arc, + pub phase: PhaseFunction, +} + +impl Interaction for SurfaceInteraction { + fn get_common(&self) -> &InteractionData { &self.common } + fn as_any(&self) -> &dyn Any { self } + + fn get_medium(&self, w: Vector3f) -> Option> { + self.medium_interface.as_ref().and_then(|interface| { + if self.n.dot(w) > 0.0 { + interface.outside.clone() + } else { + interface.inside.clone() + } + }) + } + + fn spawn_ray(&self, d: Vector3f) -> Ray { + let mut ray = Ray::spawn(&self.pi(), &self.n, self.time(), d); + ray.medium = self.get_medium(d); + ray + } + + fn spawn_ray_to_point(&self, p2: Point3f) -> Ray { + let mut ray = Ray::spawn_to_point(&self.pi(), &self.n, self.time(), p2); + ray.medium = self.get_medium(ray.d); + ray + } + + fn spawn_ray_to_interaction(&self, other: &dyn Interaction) -> Ray { + // Check if the other interaction is a surface to use the robust spawn method + if let Some(si_to) = other.as_any().downcast_ref::() { + let mut ray = Ray::spawn_to_interaction(&self.pi(), &self.n, self.time(), &si_to.pi(), &si_to.n); + ray.medium = self.get_medium(ray.d); + ray + } else { + // Fallback for non-surface interactions + self.spawn_ray_to_point(other.p()) + } + } +} + +impl SurfaceInteraction { + pub fn new(pi: Point3fi, uv: Point2f, wo: Vector3f, dpdu: Vector3f, dpdv: Vector3f, dndu: Normal3f, dndv: Normal3f, time: Float, flip: bool) -> Self { + let mut n = Normal3f::from(dpdu.cross(dpdv).normalize()); + let mut shading_n = n; + + if flip { + n *= -1.0; + shading_n *= -1.0; + } + + Self { + common: InteractionData { pi, time, wo }, + n, + uv, + dpdu, + dpdv, + dndu, + dndv, + shading: ShadingGeometry { n: shading_n, dpdu, dpdv, dndu, dndv }, + medium_interface: None, + face_index: false, + area_light: None, + dpdx: Vector3f::zero(), + dpdy: Vector3f::zero(), + dudx: 0.0, + dudy: 0.0, + dvdx: 0.0, + dvdy: 0.0, + } + } + + pub fn set_shading_geometry(&mut self, ns: Normal3f, dpdus: Vector3f, dpdvs: Vector3f, dndus: Normal3f, dndvs: Normal3f, orientation: bool) { + self.shading.n = ns; + if orientation { + self.n = self.n.face_forward(self.shading.n.into()); + } + self.shading.dpdu = dpdus; + self.shading.dpdv = dpdvs; + self.shading.dndu = dndus.into(); + self.shading.dndv = dndvs.into(); + } +} + +impl Interaction for MediumInteraction { + fn get_common(&self) -> &InteractionData { &self.common } + fn as_any(&self) -> &dyn Any { self } + + fn get_medium(&self, _w: Vector3f) -> Option> { + Some(self.medium.clone()) + } + + fn spawn_ray(&self, d: Vector3f) -> Ray { + Ray { + o: self.p(), + d, + time: self.time(), + medium: Some(self.medium.clone()), + differential: None, + } + } +} diff --git a/src/core/interval.rs b/src/core/interval.rs new file mode 100644 index 0000000..6976bd6 --- /dev/null +++ b/src/core/interval.rs @@ -0,0 +1,102 @@ +use crate::core::pbrt::{Float, next_float_up, next_float_down}; +use std::ops::{Add, Sub, Mul, Div}; + +#[derive(Debug, Copy, Clone, PartialEq)] +pub struct Interval { + pub low: Float, + pub high: Float, +} + +impl Interval { + pub fn new(v: Float) -> Self { + Self { low: v, high: v } + } + + pub fn new_from_endpoints(v1: Float, v2: Float) -> Self { + if v1 < v2 { + Self { low: v1, high: v2 } + } else { + Self { low: v2, high: v1 } + } + } + + pub fn new_from_value_and_error(val: Float, err: Float) -> Self { + Self { + low: val - err, + high: val + err, + } + } + + pub fn width(&self) -> Float { + self.high - self.low + } + + pub fn midpoint(&self) -> Float { + (self.low + self.high) / 2.0 + } + + pub fn contains(&self, v: Float) -> bool { + v >= self.low && v <= self.high + } + + pub fn is_empty(&self) -> bool { + self.low > self.high + } +} + +impl Default for Interval { + fn default() -> Self { + Self { + low: Float::INFINITY, + high: Float::NEG_INFINITY, + } + } +} + +impl Add for Interval { + type Output = Self; + fn add(self, rhs: Self) -> Self::Output { + Self { + low: next_float_down(self.low + rhs.low), + high: next_float_up(self.high + rhs.high), + } + } +} + +impl Sub for Interval { + type Output = Self; + fn sub(self, rhs: Self) -> Self::Output { + Self { + low: next_float_down(self.low - rhs.high), + high: next_float_up(self.high - rhs.low), + } + } +} + +impl Mul for Interval { + type Output = Self; + fn mul(self, rhs: Self) -> Self::Output { + let prods = [ + self.low * rhs.low, + self.high * rhs.low, + self.low * rhs.high, + self.high * rhs.high, + ]; + + Self { + low: next_float_down(prods[0].min(prods[1]).min(prods[2]).min(prods[3])), + high: next_float_up(prods[0].max(prods[1]).max(prods[2]).max(prods[3])), + } + } +} + +impl Div for Interval { + type Output = Self; + fn div(self, rhs: Self) -> Self::Output { + if rhs.contains(0.0) { + return Interval { low: Float::NEG_INFINITY, high: Float::INFINITY }; + } + + self * Interval::new_from_endpoints(1.0 / rhs.low, 1.0 / rhs.high) + } +} diff --git a/src/core/mod.rs b/src/core/mod.rs index c78e8a3..0681321 100644 --- a/src/core/mod.rs +++ b/src/core/mod.rs @@ -1,7 +1,12 @@ +pub mod camera; +pub mod color; +pub mod colorspace; +pub mod film; pub mod geometry; pub mod integrator; pub mod light; pub mod interaction; +pub mod interval; pub mod material; pub mod medium; pub mod pbrt; diff --git a/src/core/pbrt.rs b/src/core/pbrt.rs index 4856305..b0b687b 100644 --- a/src/core/pbrt.rs +++ b/src/core/pbrt.rs @@ -1,4 +1,5 @@ use num_traits::Num; +use std::ops::{Add, Mul}; pub type Float = f32; @@ -54,3 +55,122 @@ where } r } + +#[inline] +pub fn float_to_bits(f: Float) -> u32 { + f.to_bits() +} + +#[inline] +pub fn bits_to_float(ui: u32) -> Float { + f32::from_bits(ui) +} + +// Corresponding types for 64-bit floats +#[inline] +pub fn f64_to_bits(f: f64) -> u64 { + f.to_bits() +} + +#[inline] +pub fn next_float_up(v: Float) -> Float { + if v.is_infinite() && v > 0.0 { + return v; + } + let v = if v == -0.0 { 0.0 } else { v }; + + let mut ui = float_to_bits(v); + if v >= 0.0 { + ui += 1; + } else { + ui -= 1; + } + bits_to_float(ui) +} + +#[inline] +pub fn next_float_down(v: Float) -> Float { + if v.is_infinite() && v < 0.0 { + return v; + } + let v = if v == 0.0 { -0.0 } else { v }; + + let mut ui = float_to_bits(v); + if v > 0.0 { + ui -= 1; + } else { + ui += 1; + } + bits_to_float(ui) +} + +#[inline] +pub fn evaluate_polynomial(t: Float, coeffs: &[Float]) -> Option { + if coeffs.is_empty() { + return None; + } + + let mut result = coeffs[0]; + for &c in &coeffs[1..] { + result = t.mul_add(result, c); + } + Some(result) +} + +#[inline] +pub fn find_interval

(sz: usize, pred: P) -> usize +where + P: Fn(usize) -> bool, +{ + if sz <= 2 { + return 0; + } + + let mut low = 1; + let mut high = sz - 1; + + while low < high { + let mid = low + (high - low) / 2; + if pred(mid) { + low = mid + 1; + } else { + high = mid; + } + } + + let result = low - 1; + clamp_t(result, 0, sz - 2) +} + + + +#[inline] +pub fn safe_asin(x: Float) -> Float { + if x >= -1.0001 && x <= 1.0001 { + clamp_t(x, -1., 1.).asin() + } else { + panic!("Not valid value for asin") + } +} + + +#[inline] +pub fn sinx_over_x(x: Float) -> Float { + if 1. - x * x == 1. { + return 1.; + } + return x.sin() / x; +} + +#[inline] +pub fn gamma(n: i32) -> Float { + return (n as Float * MACHINE_EPSILON) / (1. - n as Float * MACHINE_EPSILON); +} + + +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum RenderingCoordinateSystem { + Camera, + CameraWorld, + World, +} diff --git a/src/core/quaternion.rs b/src/core/quaternion.rs index 58769e8..fe7652d 100644 --- a/src/core/quaternion.rs +++ b/src/core/quaternion.rs @@ -1,6 +1,8 @@ -use crate::core::geometry::Vector3f; -use crate::core::pbrt::Float; +use crate::core::geometry::{Vector3f, Dot, Normed}; +use crate::core::pbrt::{safe_asin, sinx_over_x, Float}; +use std::ops::{Index, IndexMut}; use std::ops::{Add, AddAssign, Sub, SubAssign, Mul, MulAssign, Div, DivAssign, Neg}; +use std::f32::consts::PI; #[derive(Copy, Clone, PartialEq)] pub struct Quaternion { @@ -56,3 +58,52 @@ impl Neg for Quaternion { type Output = Self; fn neg(self) -> Self { Self { v: -self.v, w: -self.w }} } + +impl Index for Quaternion { + type Output = Float; + #[inline] + fn index(&self, index: usize) -> &Self::Output { + match index { + 0 => &self.v[0], + 1 => &self.v[1], + 2 => &self.v[2], + 3 => &self.w, + _ => panic!("Quaternion index out of bounds: {}", index), + } + } +} + +impl IndexMut for Quaternion { + #[inline] + fn index_mut(&mut self, index: usize) -> &mut Self::Output { + match index { + 0 => &mut self.v[0], + 1 => &mut self.v[1], + 2 => &mut self.v[2], + 3 => &mut self.w, + _ => panic!("Quaternion index out of bounds: {}", index), + } + } +} + +impl Quaternion { + pub fn dot(&self, rhs: Quaternion) -> Float { + self.v.dot(rhs.v) + self.w * rhs.w + } + + #[inline] + pub fn angle_between(&self, rhs: Quaternion) -> Float { + if self.dot(rhs) < 0.0 { + return PI - 2. * safe_asin((self.v + rhs.v).norm() / 2.) + } else { + return 2. * safe_asin((rhs.v - self.v).norm() / 2.) + } + } + + pub fn slerp(t: Float, q1: Quaternion, q2: Quaternion) -> Quaternion { + let theta = q1.angle_between(q2); + let sin_theta_over_theta = sinx_over_x(theta); + return q1 * (1. - t) * sinx_over_x((1. - t) * theta) / sin_theta_over_theta + + q2 * t * sinx_over_x(t * theta) / sin_theta_over_theta; + } +} diff --git a/src/core/spectrum.rs b/src/core/spectrum.rs index 62ab452..31a9c31 100644 --- a/src/core/spectrum.rs +++ b/src/core/spectrum.rs @@ -1,8 +1,591 @@ -use crate::core::pbrt::Float; +use std::fmt; +use std::ops::{Add, AddAssign, Sub, SubAssign, Mul, MulAssign, Div, DivAssign, Neg, Index, IndexMut}; +use once_cell::sync::Lazy; +use std::sync::Arc; -#[derive(Debug, Default, Copy, Clone)] -pub struct RGBSpectrum { - pub c: [Float; 3], +use crate::core::colorspace::RGBColorspace; +use crate::core::pbrt::Float; +use crate::core::color::{RGB, XYZ, RGBSigmoidPolynomial}; + +pub const CIE_Y_INTEGRAL: Float = 106.856895; +pub const N_SPECTRUM_SAMPLES: usize = 1200; +pub const LAMBDA_MIN: i32 = 360; +pub const LAMBDA_MAX: i32 = 830; + +#[derive(Debug, Copy, Clone)] +pub struct SampledSpectrum { + values: [Float; N_SPECTRUM_SAMPLES], } -pub struct Spectrum; +impl Default for SampledSpectrum { + fn default() -> Self { Self { values: [0.0; N_SPECTRUM_SAMPLES] } } +} + +impl SampledSpectrum { + pub fn new(c: Float) -> Self { + Self { + values: [c; N_SPECTRUM_SAMPLES], + } + } + + pub fn from_vector(v: Vec) -> Self { + let mut s = Self::new(0.0); + for i in 0..N_SPECTRUM_SAMPLES { + s.values[i] = v[i]; + } + s + } + + pub fn has_nans(&self) -> bool { + self.values.iter().any(|&v| v.is_nan()) + } + + pub fn min_component_value(&self) -> Float { + self.values.iter().fold(Float::INFINITY, |a, &b| a.min(b)) + } + + pub fn max_component_value(&self) -> Float { + self.values.iter().fold(Float::NEG_INFINITY, |a, &b| a.max(b)) + } + + pub fn average(&self) -> Float { + self.values.iter().sum::() / (N_SPECTRUM_SAMPLES as Float) + } + + pub fn safe_div(&self, rhs: SampledSpectrum) -> Self { + let mut r = SampledSpectrum::new(0.0); + for i in 0..N_SPECTRUM_SAMPLES { + r.values[i] = if rhs[i] != 0.0 {self.values[i] / rhs.values[i]} else { 0.0 } + } + r + } + + pub fn to_xyz(&self, lambda: &SampledWavelengths) -> XYZ { + let x = spectra::X.sample(lambda); + let y = spectra::Y.sample(lambda); + let z = spectra::Z.sample(lambda); + let pdf = lambda.pdf(); + XYZ::new( + (*self * x).safe_div(pdf).average(), + (*self * y).safe_div(pdf).average(), + (*self * z).safe_div(pdf).average(), + ) / CIE_Y_INTEGRAL + } + + pub fn to_rgb(&self, lambda: &SampledWavelengths, c: &RGBColorspace) -> RGB { + let xyz = self.to_xyz(lambda); + c.to_rgb(xyz) + } +} + + +impl Index for SampledSpectrum { + type Output = Float; + fn index(&self, i: usize) -> &Self::Output { &self.values[i] } +} + +impl IndexMut for SampledSpectrum { + fn index_mut(&mut self, i: usize) -> &mut Self::Output { &mut self.values[i] } +} + +impl Add for SampledSpectrum { + type Output = Self; + fn add(self, rhs: Self) -> Self { + let mut ret = self; + for i in 0..N_SPECTRUM_SAMPLES { + ret.values[i] += rhs.values[i]; + } + ret + } +} + +impl AddAssign for SampledSpectrum { + fn add_assign(&mut self, rhs: Self) { + for i in 0..N_SPECTRUM_SAMPLES { + self.values[i] += rhs.values[i]; + } + } +} + +impl Sub for SampledSpectrum { + type Output = Self; + fn sub(self, rhs: Self) -> Self { + let mut ret = self; + for i in 0..N_SPECTRUM_SAMPLES { + ret.values[i] -= rhs.values[i]; + } + ret + } +} + +impl SubAssign for SampledSpectrum { + fn sub_assign(&mut self, rhs: Self) { + for i in 0..N_SPECTRUM_SAMPLES { + self.values[i] -= rhs.values[i]; + } + } +} + +impl Sub for Float { + type Output = SampledSpectrum; + fn sub(self, rhs: SampledSpectrum) -> SampledSpectrum { + let mut ret = SampledSpectrum::new(0.0); + for i in 0..N_SPECTRUM_SAMPLES { + ret.values[i] = self - rhs.values[i]; + } + ret + } +} + +impl Mul for SampledSpectrum { + type Output = Self; + fn mul(self, rhs: Self) -> Self { + let mut ret = self; + for i in 0..N_SPECTRUM_SAMPLES { + ret.values[i] *= rhs.values[i]; + } + ret + } +} + +impl MulAssign for SampledSpectrum { + fn mul_assign(&mut self, rhs: Self) { + for i in 0..N_SPECTRUM_SAMPLES { + self.values[i] *= rhs.values[i]; + } + } +} + +impl Mul for SampledSpectrum { + type Output = Self; + fn mul(self, rhs: Float) -> Self { + let mut ret = self; + for i in 0..N_SPECTRUM_SAMPLES { + ret.values[i] *= rhs; + } + ret + } +} + +impl MulAssign for SampledSpectrum { + fn mul_assign(&mut self, rhs: Float) { + for i in 0..N_SPECTRUM_SAMPLES { + self.values[i] *= rhs; + } + } +} + +impl DivAssign for SampledSpectrum { + fn div_assign(&mut self, rhs: Self) { + for i in 0..N_SPECTRUM_SAMPLES { + debug_assert_ne!(0.0, rhs.values[i]); + self.values[i] /= rhs.values[i]; + } + } +} + +impl Div for SampledSpectrum { + type Output = Self; + + fn div(self, rhs: Float) -> Self::Output { + debug_assert_ne!(rhs, 0.0); + let mut ret = self; + for i in 0..N_SPECTRUM_SAMPLES { + ret.values[i] /= rhs; + } + ret + } +} + +impl DivAssign for SampledSpectrum { + fn div_assign(&mut self, rhs: Float) { + debug_assert_ne!(rhs, 0.0); + for i in 0..N_SPECTRUM_SAMPLES { + self.values[i] /= rhs; + } + } +} + +impl Neg for SampledSpectrum { + type Output = Self; + + fn neg(self) -> Self::Output { + let mut ret = SampledSpectrum::new(0.0); + for i in 0..N_SPECTRUM_SAMPLES { + ret.values[i] = -self.values[i]; + } + ret + } +} + +#[derive(Debug, Copy, Clone, PartialEq)] +pub struct SampledWavelengths { + pub lambda: [Float; N_SPECTRUM_SAMPLES], + pub pdf: [Float; N_SPECTRUM_SAMPLES], +} + +impl SampledWavelengths { + pub fn pdf(&self) -> SampledSpectrum { + SampledSpectrum::from_vector(self.pdf.to_vec()) + } + + pub fn secondary_terminated(&self) -> bool { + for i in 1..N_SPECTRUM_SAMPLES { + if self.pdf[i] != 0.0 { + return false + } + } + true + } + + pub fn terminate_secondary(&mut self) { + if !self.secondary_terminated() { + for i in 1..N_SPECTRUM_SAMPLES { + self.pdf[i] = 0.0; + } + self.pdf[0] /= N_SPECTRUM_SAMPLES as Float; + } + } + + pub fn sample_uniform(u: Float, lambda_min: Float, lambda_max: Float) -> Self { + let mut lambda = [0.0; N_SPECTRUM_SAMPLES]; + lambda[0] = crate::core::pbrt::lerp(u, lambda_min, lambda_min); + let delta = (lambda_max - lambda_min) / N_SPECTRUM_SAMPLES as Float; + for i in 1..N_SPECTRUM_SAMPLES { + lambda[i] = lambda[i - 1] + delta; + if lambda[i] > lambda_max { + lambda[i] = lambda_min + (lambda[i] - lambda_max); + } + } + + let mut pdf = [0.0; N_SPECTRUM_SAMPLES]; + for i in 0..N_SPECTRUM_SAMPLES { + pdf[i] = 1.0 / (lambda_max - lambda_min); + } + + Self { lambda, pdf } + } + + pub fn sample_visible_wavelengths(u: Float) -> Float { + 538.0 - 138.888889 * Float::atanh(0.85691062 - 1.82750197 * u) + } + + pub fn visible_wavelengths_pdf(lambda: Float) -> Float { + if lambda < 360.0 || lambda > 830.0 { return 0.0; } + 0.0039398042 / (Float::cosh(0.0072 * (lambda - 538.0))).sqrt() + } + + pub fn sample_visible(u: Float) -> Self { + let mut lambda = [0.0; N_SPECTRUM_SAMPLES]; + let mut pdf = [0.0; N_SPECTRUM_SAMPLES]; + + for i in 0..N_SPECTRUM_SAMPLES { + let mut up = u + i as Float / N_SPECTRUM_SAMPLES as Float; + if up > 1.0 { up -= 1.0; } + lambda[i] = Self::sample_visible_wavelengths(up); + pdf[i] = Self::visible_wavelengths_pdf(lambda[i]); + } + Self { lambda, pdf } + } +} + +impl Index for SampledWavelengths { + type Output = Float; + fn index(&self, i: usize) -> &Self::Output { &self.lambda[i] } +} + +impl IndexMut for SampledWavelengths { + fn index_mut(&mut self, i: usize) -> &mut Self::Output { &mut self.lambda[i] } +} + +#[derive(Debug, Clone)] +pub enum Spectrum { + Constant(ConstantSpectrum), + DenselySampled(DenselySampledSpectrum), + PiecewiseLinear(PiecewiseLinearSpectrum), + Blackbody(BlackbodySpectrum), + RGBAlbedo(RGBAlbedoSpectrum), +} + +impl Spectrum { + pub fn sample_at(&self, lambda: Float) -> Float { + match self { + Spectrum::Constant(s) => s.sample_at(lambda), + Spectrum::DenselySampled(s) => s.sample_at(lambda), + Spectrum::PiecewiseLinear(s) => s.sample_at(lambda), + Spectrum::Blackbody(s) => s.sample_at(lambda), + Spectrum::RGBAlbedo(s) => s.sample_at(lambda), + } + } + + pub fn max_value(&self) -> Float { + match self { + Spectrum::Constant(_s) => 0.0, + Spectrum::DenselySampled(_s) => 0.0, + Spectrum::PiecewiseLinear(_s) => 0.0, + Spectrum::Blackbody(_s) => 0.0, + Spectrum::RGBAlbedo(s) => s.max_value(), + } + } + + pub fn sample(&self, wavelengths: &SampledWavelengths) -> SampledSpectrum { + let mut s = SampledSpectrum::default(); + for i in 0..N_SPECTRUM_SAMPLES { + s[i] = self.sample_at(wavelengths[i]); + } + s + } + + pub fn to_xyz(&self) -> XYZ { + XYZ::new(inner_product(&spectra::X, self), + inner_product(&spectra::Y, self), + inner_product(&spectra::Z, self))/CIE_Y_INTEGRAL + } +} + +pub fn inner_product(f: &Spectrum, g: &Spectrum) -> Float { + let mut integral = 0.0; + for lambda in LAMBDA_MIN..=LAMBDA_MAX { + integral += f.sample_at(lambda as f32) * g.sample_at(lambda as f32); + } + integral +} + + +#[derive(Debug, Clone, Copy)] +pub struct ConstantSpectrum { + c: Float, +} + +impl ConstantSpectrum { + pub fn new(c: Float) -> Self { Self { c } } + + pub fn sample_at(&self, _lambda: Float) -> Float { + self.c + } + + fn max_value(&self) -> Float { self.c } +} + +#[derive(Debug, Clone, Copy)] +pub struct RGBAlbedoSpectrum { + rsp: RGBSigmoidPolynomial, +} + +impl RGBAlbedoSpectrum { + pub fn new(cs: &RGBColorspace, rgb: RGB) -> Self { + Self { rsp: cs.to_rgb_coeffs(rgb) } + } + + pub fn sample_at(&self, lambda: Float) -> Float { + self.rsp.evaluate(lambda) + } + + pub fn max_value(&self) -> Float { + self.rsp.max_value() + } + + fn sample(&self, lambda: &SampledWavelengths) -> SampledSpectrum { + let mut s = SampledSpectrum::default(); + for i in 0..N_SPECTRUM_SAMPLES { + s[i] = self.rsp.evaluate(lambda[i]); + } + s + } +} + +#[derive(Debug, Clone, Copy)] +pub struct UnboundedRGBSpectrum { + scale: Float, + rsp: RGBSigmoidPolynomial, +} + +impl UnboundedRGBSpectrum { + pub fn new(cs: RGBColorspace, rgb: RGB) -> Self { + let m = rgb.r.max(rgb.g).max(rgb.b); + let scale = 2.0 * m; + let lambda; + if scale == 1.0 { + lambda = rgb / scale; + } else { + lambda = RGB::new(0.0, 0.0, 0.0); + } + Self { scale, rsp: cs.to_rgb_coeffs(lambda) } + } + pub fn sample_at(&self, lambda: Float) -> Float { + self.scale * self.rsp.evaluate(lambda) + } + + pub fn max_value(&self) -> Float { + self.scale * self.rsp.max_value() + } +} + +#[derive(Debug, Clone, Default)] +pub struct RGBIlluminantSpectrum { + scale: Float, + rsp: RGBSigmoidPolynomial, + illuminant: Option>, +} + +impl RGBIlluminantSpectrum { + pub fn sample_at(&self, lambda: Float) -> Float { + match &self.illuminant { + Some(illuminant) => self.scale * self.rsp.evaluate(lambda) * illuminant.sample_at(lambda), + None => 0.0, + } + } + + pub fn max_value(&self) -> Float { + match &self.illuminant { + Some(illuminant) => self.scale * self.rsp.max_value() * illuminant.max_value(), + None => 0.0, + } + } +} + + +#[derive(Debug, Clone)] +pub struct DenselySampledSpectrum { + lambda_min: i32, + lambda_max: i32, + values: Vec, +} + +impl DenselySampledSpectrum { + pub fn new(lambda_min: i32, lambda_max: i32) -> Self { + let n_values = (lambda_max - lambda_min + 1).max(0) as usize; + Self { lambda_min, lambda_max, values: vec![0.0; n_values] } + } + + pub fn from_spectrum(spec: &Spectrum, lambda_min: i32, lambda_max: i32) -> Self { + let mut s = Self::new(lambda_min, lambda_max); + if s.values.is_empty() { return s; } + for lambda in lambda_min..=lambda_max { + let index = (lambda - lambda_min) as usize; + s.values[index] = spec.sample_at(lambda as Float); + } + s + } + + pub fn from_function(f: F, lambda_min: i32, lambda_max: i32) -> Self where F: Fn(Float) -> Float { + let mut s = Self::new(lambda_min, lambda_max); + if s.values.is_empty() { return s; } + for lambda in lambda_min..=lambda_max { + let index = (lambda - lambda_min) as usize; + s.values[index] = f(lambda as Float); + } + s + } + + pub fn sample_at(&self, lambda: Float) -> Float { + let offset = (lambda.round() as i32) - self.lambda_min; + if offset < 0 || offset as usize >= self.values.len() { + 0.0 + } else { + self.values[offset as usize] + } + } + + pub fn max_value(&self) -> Float { + self.values.iter().fold(Float::MIN, |a,b| a.max(*b)) + } +} + +#[derive(Debug, Clone)] +pub struct PiecewiseLinearSpectrum { + samples: Vec<(Float, Float)>, +} + +impl PiecewiseLinearSpectrum { + pub fn from_interleaved(data: &[Float]) -> Self { + assert!(data.len() % 2 == 0, "Interleaved data must have an even number of elements"); + let mut samples = Vec::new(); + for pair in data.chunks(2) { + samples.push((pair[0], pair[1])); + } + // PBRT requires the samples to be sorted by wavelength for interpolation. + samples.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap()); + Self { samples } + } + + fn sample_at(&self, lambda: Float) -> Float { + if self.samples.is_empty() { + return 0.0; + } + // Handle boundary conditions + if lambda <= self.samples[0].0 { + return self.samples[0].1; + } + if lambda >= self.samples.last().unwrap().0 { + return self.samples.last().unwrap().1; + } + + let i = self.samples.partition_point(|s| s.0 < lambda); + let s1 = self.samples[i - 1]; + let s2 = self.samples[i]; + + let t = (lambda - s1.0) / (s2.0 - s1.0); + (1.0 - t) * s1.1 + t * s2.1 + } +} + +#[derive(Debug, Clone, Copy)] +pub struct BlackbodySpectrum { + temperature: Float, + normalization_factor: Float, +} + +// Planck's Law +impl BlackbodySpectrum { + const C: Float = 299792458.0; + const H: Float = 6.62606957e-34; + const KB: Float = 1.3806488e-23; + + pub fn new(temperature: Float) -> Self { + let lambda_max = 2.8977721e-3 / temperature * 1e9; + let max_val = Self::planck_law(lambda_max, temperature); + Self { + temperature, + normalization_factor: if max_val > 0.0 { 1.0 / max_val } else { 0.0 }, + } + } + + fn planck_law(lambda_nm: Float, temp: Float) -> Float { + if temp <= 0.0 { return 0.0; } + let lambda_m = lambda_nm * 1e-9; // Convert nm to meters + let c1 = 2.0 * Self::H * Self::C * Self::C; + let c2 = (Self::H * Self::C) / Self::KB; + + let numerator = c1 / lambda_m.powi(5); + let denominator = (c2 / (lambda_m * temp)).exp() - 1.0; + + if denominator.is_infinite() { 0.0 } else { numerator / denominator } + } + + fn sample_at(&self, lambda: Float) -> Float { + Self::planck_law(lambda, self.temperature) * self.normalization_factor + } +} + +#[derive(Debug, Clone, Copy)] +pub struct RGBSpectrum { pub c: [Float; 3] } + +#[derive(Debug, Clone, Copy)] +pub struct RGBUnboundedSpectrum(pub RGBSpectrum); + +pub mod spectra { + use super::*; + pub static X: Lazy = Lazy::new(|| { + let dss = DenselySampledSpectrum::from_piecewise_linear(&CIE_X_DATA); + Spectrum::DenselySampled(dss) + }); + pub static Y: Lazy = Lazy::new(|| { + let dss = DenselySampledSpectrum::from_piecewise_linear(&CIE_Y_DATA); + Spectrum::DenselySampled(dss) + }); + pub static Z: Lazy = Lazy::new(|| { + let dss = DenselySampledSpectrum::from_piecewise_linear(&CIE_Z_DATA); + Spectrum::DenselySampled(dss) + }); +} diff --git a/src/core/transform.rs b/src/core/transform.rs index 96cf857..4f2dc22 100644 --- a/src/core/transform.rs +++ b/src/core/transform.rs @@ -3,17 +3,20 @@ use std::ops::{Add, Div, Index, IndexMut, Mul}; use std::fmt; use std::fmt::{Display, Debug}; -use crate::core::geometry::{Vector, Point, Vector3f, Point3f}; +use crate::core::geometry::{Vector, Point, Vector3f, Point3f, Point3fi}; +use crate::core::color::{RGB, XYZ}; +use crate::core::pbrt::{gamma, Float}; +use crate::core::quaternion::Quaternion; #[derive(Copy, Clone)] pub struct SquareMatrix { - m: [[T; N]; N], + pub m: [[T; N]; N], } impl SquareMatrix { - fn new() -> Self where T: Copy + Zero { - Self::zero() + pub fn new(data: [[T; N]; N]) -> Self { + Self { m: data } } fn identity() -> Self where T: Copy + Zero + One { @@ -24,6 +27,7 @@ impl SquareMatrix Self { m } } + pub fn zero() -> Self where T: Copy + Zero { SquareMatrix { m: [[T::zero(); N]; N] } } @@ -43,6 +47,21 @@ impl SquareMatrix true } + fn trace(&self) -> T + where + T: Zero + Copy + { + let mut sum = T::zero(); + for i in 0..N { + for j in 0..N { + if i == j { + sum = sum + self.m[i][j]; + } + } + } + sum + } + pub fn transpose(&self) -> Self where T: Copy + Zero { let mut result = SquareMatrix::zero(); for i in 0..N { @@ -61,11 +80,11 @@ impl SquareMatrix { pub fn inverse(&self) -> Option { if N == 0 { - return Some(Self::new()); + return Some(Self::identity()); } let mut mat = self.m; - let mut inv = Self::new(); + let mut inv = Self::identity(); for i in 0..N { let mut pivot_row = i; @@ -109,53 +128,80 @@ impl SquareMatrix } pub fn determinant(&self) -> T { - if N == 0 { - return T::one(); - } + let m = &self.m; - let mut lum = self.m; // Make a mutable copy of the matrix data - let mut parity = 0; // Number of row swaps + match N { + 0 => T::one(), + 1 => m[0][0], + 2 => m[0][0] * m[1][1] - m[0][1] * m[1][0], + 3 => { + let a = m[0][0]; let b = m[0][1]; let c = m[0][2]; + let d = m[1][0]; let e = m[1][1]; let f = m[1][2]; + let g = m[2][0]; let h = m[2][1]; let i = m[2][2]; - for i in 0..N { - // Pivoting: find the row with the largest element in the current column - let mut max_row = i; - for k in (i + 1)..N { - if lum[k][i].abs() > lum[max_row][i].abs() { - max_row = k; + a * (e * i - f * h) - + b * (d * i - f * g) + + c * (d * h - e * g) + } + 4 => { + let det3 = |m11: T, m12: T, m13: T, + m21: T, m22: T, m23: T, + m31: T, m32: T, m33: T| -> T { + m11 * (m22 * m33 - m23 * m32) - + m12 * (m21 * m33 - m23 * m31) + + m13 * (m21 * m32 - m22 * m31) + }; + + let c0 = det3(m[1][1], m[1][2], m[1][3], m[2][1], m[2][2], m[2][3], m[3][1], m[3][2], m[3][3]); + let c1 = det3(m[1][0], m[1][2], m[1][3], m[2][0], m[2][2], m[2][3], m[3][0], m[3][2], m[3][3]); + let c2 = det3(m[1][0], m[1][1], m[1][3], m[2][0], m[2][1], m[2][3], m[3][0], m[3][1], m[3][3]); + let c3 = det3(m[1][0], m[1][1], m[1][2], m[2][0], m[2][1], m[2][2], m[3][0], m[3][1], m[3][2]); + + m[0][0] * c0 - m[0][1] * c1 + m[0][2] * c2 - m[0][3] * c3 + } + _ => { // Fallback to LU decomposition for N > 4 + let mut lum = self.m; + let mut parity = 0; + + for i in 0..N { + let mut max_row = i; + for k in (i + 1)..N { + if lum[k][i].abs() > lum[max_row][i].abs() { + max_row = k; + } + } + + if max_row != i { + lum.swap(i, max_row); + parity += 1; + } + + // Singular matrix + if lum[i][i] == T::zero() { + return T::zero(); + } + + // Gaussian elimination + for j in (i + 1)..N { + let factor = lum[j][i] / lum[i][i]; + for k in i..N { + lum[j][k] = lum[j][k] - factor * lum[i][k]; + } + } } - } - // If the pivot is on a different row, swap them - if max_row != i { - lum.swap(i, max_row); - parity += 1; - } - - // If the pivot element is zero, the determinant is zero - // (the matrix is singular). - if lum[i][i] == T::zero() { - return T::zero(); - } - - // Gaussian elimination to form the upper-triangular matrix - for j in (i + 1)..N { - let factor = lum[j][i] / lum[i][i]; - for k in i..N { - lum[j][k] = lum[j][k] - factor * lum[i][k]; + let mut det = T::one(); + for i in 0..N { + det = det * lum[i][i]; } + + if parity % 2 != 0 { + det = -det; + } + + det } } - - let mut det = T::one(); - for i in 0..N { - det = det * lum[i][i]; - } - - if parity % 2 != 0 { - det = -det; - } - - det } } @@ -272,6 +318,29 @@ where } } +impl Mul for SquareMatrix +where + T: Copy + Mul, +{ + type Output = Self; + fn mul(mut self, rhs: T) -> Self::Output { + for row in self.m.iter_mut() { + for element in row.iter_mut() { + *element = *element * rhs; + } + } + self + } +} + +impl Mul> for Float +{ + type Output = SquareMatrix; + fn mul(self, rhs: SquareMatrix) -> Self::Output { + rhs * self + } +} + impl Div for SquareMatrix where T: Copy + Div, @@ -320,6 +389,61 @@ where } } +impl SquareMatrix { + pub fn transform_to_xyz(&self, rgb: RGB) -> XYZ { + let r = rgb[0]; + let g = rgb[1]; + let b = rgb[2]; + + let x = self.m[0][0] * r + self.m[0][1] * g + self.m[0][2] * b; + let y = self.m[1][0] * r + self.m[1][1] * g + self.m[1][2] * b; + let z = self.m[2][0] * r + self.m[2][1] * g + self.m[2][2] * b; + + XYZ::new(x, y, z) + } + + // This name is also perfectly clear + pub fn transform_to_rgb(&self, xyz: XYZ) -> RGB { + let x = xyz[0]; + let y = xyz[1]; + let z = xyz[2]; + + let r = self.m[0][0] * x + self.m[0][1] * y + self.m[0][2] * z; + let g = self.m[1][0] * x + self.m[1][1] * y + self.m[1][2] * z; + let b = self.m[2][0] * x + self.m[2][1] * y + self.m[2][2] * z; + + RGB::new(r, g, b) + } +} + +impl Mul for SquareMatrix +{ + type Output = XYZ; + fn mul(self, rhs: XYZ) -> Self::Output { + let mut result_arr = [0.0; 3]; + for i in 0..3 { + for j in 0..3 { + result_arr[i] = result_arr[i] + self.m[i][j] * rhs[j]; + } + } + XYZ::new(result_arr[0], result_arr[1], result_arr[2]) + } +} + +impl Mul for SquareMatrix +{ + type Output = RGB; + fn mul(self, rhs: RGB) -> Self::Output { + let mut result_arr = [0.0; 3]; + for i in 0..3 { + for j in 0..3 { + result_arr[i] = result_arr[i] + self.m[i][j] * rhs[j]; + } + } + RGB::new(result_arr[0], result_arr[1], result_arr[2]) + } +} + #[derive(Copy, Clone)] pub struct Transform { m: SquareMatrix, @@ -327,7 +451,11 @@ pub struct Transform { } impl Transform { - pub fn new(m: &SquareMatrix) -> Self { + pub fn new(m: &SquareMatrix, m_inv: &SquareMatrix) -> Self { + Self { m: *m, m_inv: *m_inv } + } + + pub fn from_matrix(m: &SquareMatrix) -> Self { let inv = match m.inverse() { Some(inv) => inv , None => SquareMatrix::zero(), @@ -336,6 +464,12 @@ impl Transform { Self { m: *m, m_inv: inv } } + pub fn identity() -> Self { + let m: SquareMatrix = SquareMatrix::identity(); + let m_inv = m.inverse().expect("Singular matrix"); + Self { m, m_inv } + } + pub fn is_identity(&self) -> bool { self.m.is_identity() } @@ -344,17 +478,154 @@ impl Transform { Self { m: self.m_inv, m_inv: self.m} } - pub fn translate(delta: Vector3f) -> Transform { - Transform { - m: Matrix::new( - 1.0, 0.0, 0.0, delta.x, 0.0, 1.0, 0.0, delta.y, 0.0, 0.0, 1.0, delta.z, 0.0, 0.0, - 0.0, 1.0, - ), - m_inv: Matrix4x4::new( - 1.0, 0.0, 0.0, -delta.x(), 0.0, 1.0, 0.0, -delta.y(), 0.0, 0.0, 1.0, -delta.z(), 0.0, - 0.0, 0.0, 1.0, - ), + pub fn apply_inverse(&self, p: Point) -> Point { + *self * p + } + + pub fn apply_inverse_vector(&self, v: Vector) -> Vector { + *self * v + } + +} + +impl Transform { + pub fn apply(&self, pi: &Point3fi) -> Point3fi { + let p = pi.midpoint(); + let x = p.x(); + let y = p.y(); + let z = p.z(); + let xp = self.m[0][0] * x + self.m[0][1] * y + self.m[0][2] * z + self.m[0][3]; + let yp = self.m[1][0] * x + self.m[1][1] * y + self.m[1][2] * z + self.m[1][3]; + let zp = self.m[2][0] * x + self.m[2][1] * y + self.m[2][2] * z + self.m[2][3]; + let wp = self.m[3][0] * x + self.m[3][1] * y + self.m[3][2] * z + self.m[3][3]; + let mut p_error = Vector3f::default(); + + if pi.is_exact() { + p_error[0] = gamma(3) * (self.m[0][0] * x).abs() + (self.m[0][1] * y).abs() + + (self.m[0][2] + z).abs() + self.m[0][3].abs(); + p_error[1] = gamma(3) * (self.m[1][0] * x).abs() + (self.m[1][1] * y).abs() + + (self.m[1][2] + z).abs() + self.m[1][3].abs(); + p_error[2] = gamma(3) * (self.m[2][0] * x).abs() + (self.m[2][1] * y).abs() + + (self.m[2][2] + z).abs() + self.m[2][3].abs(); } + + if wp == 1. { + return Point3fi::new_with_error(Point([xp, yp, zp]), p_error); + } else { + return Point3fi::new_with_error(Point([xp/wp, yp/wp, zp/wp]), p_error) + } + } + + pub fn to_quaternion(&self) -> Quaternion { + let trace = self.m.trace(); + let mut quat = Quaternion::default(); + if trace > 0. { + let mut s = (trace + 1.).sqrt(); + quat.w = 0.5 * s; + s = 0.5 / s; + quat.v[0] = (self.m[2][1] - self.m[1][2]) * s; + quat.v[1] = (self.m[2][1] - self.m[1][2]) * s; + quat.v[2] = (self.m[2][1] - self.m[1][2]) * s; + } else { + let nxt = [1, 2, 0]; + let mut q = [0.; 3]; + let mut i = 0; + if self.m[1][1] > self.m[0][0] { + i = 1; + } + if self.m[2][2] > self.m[i][i] { + i = 2; + } + let j = nxt[i]; + let k = nxt[j]; + let mut s = (self.m[i][i] - (self.m[j][j] + self.m[k][k]) + 1.).sqrt(); + q[i] = 0.5 * s; + if s != 0. { + s = 0.5 / s; + } + quat.w = (self.m[k][j] - self.m[j][k]) * s; + q[j] = (self.m[j][i] + self.m[i][j]) * s; + q[k] = (self.m[k][i] + self.m[i][k]) * s; + quat.v[0] = q[0]; + quat.v[1] = q[1]; + quat.v[2] = q[2]; + } + quat + } + + pub fn translate(delta: Vector3f) -> Self { + Transform { + m: SquareMatrix { + m: [ + [1.0, 0.0, 0.0, delta.x()], + [0.0, 1.0, 0.0, delta.y()], + [0.0, 0.0, 1.0, delta.z()], + [0.0, 0.0, 0.0, 1.0], + ], + }, + m_inv: SquareMatrix { + m: [ + [1.0, 0.0, 0.0, -delta.x()], + [0.0, 1.0, 0.0, -delta.y()], + [0.0, 0.0, 1.0, -delta.z()], + [0.0, 0.0, 0.0, 1.0], + ], + }, + } + } + + + pub fn scale(x: Float, y: Float, z: Float) -> Self { + let m = SquareMatrix::new([ + [x, 0., 0., 0.], + [0., y, 0., 0.], + [0., 0., z, 0.], + [0., 0., 0., 1.], + ]); + + let m_inv = SquareMatrix::new([ + [1./x, 0., 0., 0.], + [0., 1./y, 0., 0.], + [0., 0., 1./z, 0.], + [0., 0., 0., 1.], + ]); + Self { m, m_inv } + } + + pub fn rotate_x(theta: Float) -> Self { + let sin_theta = theta.to_radians().sin(); + let cos_theta = theta.to_radians().cos(); + let m = SquareMatrix::new([ + [1., 0., 0., 0.], + [0., cos_theta, -sin_theta, 0.], + [0., sin_theta, cos_theta, 0.], + [0., 0., 0., 1.], + ]); + Self { m, m_inv: m.transpose() } + } + + pub fn rotate_y(theta: Float) -> Self { + let sin_theta = theta.to_radians().sin(); + let cos_theta = theta.to_radians().cos(); + let m = SquareMatrix::new([ + [cos_theta, 0., sin_theta, 0.], + [0., 1., 0., 0.], + [-sin_theta, 0., cos_theta, 0.], + [0., 0., 0., 1.], + ]); + Self { m, m_inv: m.transpose() } + } + + pub fn rotate_z(theta: Float) -> Self { + let sin_theta = theta.to_radians().sin(); + let cos_theta = theta.to_radians().cos(); + let m = SquareMatrix::new([ + [cos_theta, -sin_theta, 0., 0.], + [sin_theta, cos_theta, 0., 0.], + [0., 0., 1., 0.], + [0., 0., 0., 1.], + ]); + Self { m, m_inv: m.transpose() } } } @@ -364,6 +635,158 @@ impl PartialEq for Transform { } } +impl Mul for Transform +where T: NumFloat +{ + type Output = Self; + fn mul(self, rhs: Self) -> Self::Output { + Transform { m: self.m * rhs.m, m_inv: rhs.m_inv * self.m_inv } + } +} + +impl Mul> for Float +{ + type Output = Transform; + fn mul(self, rhs: Transform) -> Self::Output { + Transform { m: rhs.m * self , m_inv: rhs.m_inv * self } + } +} + +impl Mul> for Transform +where + T: NumFloat +{ + type Output = Point; + fn mul(self, p: Point) -> Self::Output { + let xp = self.m[0][0] * p.x() + self.m[0][1] * p.y() + self.m[0][2] * p.z() + self.m[0][3]; + let yp = self.m[1][0] * p.x() + self.m[1][1] * p.y() + self.m[1][2] * p.z() + self.m[1][3]; + let zp = self.m[2][0] * p.x() + self.m[2][1] * p.y() + self.m[2][2] * p.z() + self.m[2][3]; + let wp = self.m[3][0] * p.x() + self.m[3][1] * p.y() + self.m[3][2] * p.z() + self.m[3][3]; + + if wp == T::one() { + Point([xp, yp, zp]) + } else { + Point([xp/wp, yp/wp, zp/wp]) + } + } +} + +impl Mul> for Transform +where + T: NumFloat +{ + type Output = Vector; + fn mul(self, p: Vector) -> Self::Output { + let xp = self.m[0][0] * p.x() + self.m[0][1] * p.y() + self.m[0][2] * p.z() + self.m[0][3]; + let yp = self.m[1][0] * p.x() + self.m[1][1] * p.y() + self.m[1][2] * p.z() + self.m[1][3]; + let zp = self.m[2][0] * p.x() + self.m[2][1] * p.y() + self.m[2][2] * p.z() + self.m[2][3]; + let wp = self.m[3][0] * p.x() + self.m[3][1] * p.y() + self.m[3][2] * p.z() + self.m[3][3]; + + if wp == T::one() { + Vector([xp, yp, zp]) + } else { + Vector([xp/wp, yp/wp, zp/wp]) + } + } +} + +impl From for Transform { + fn from(q: Quaternion) -> Self { + let xx = q.v.x() * q.v.x(); + let yy = q.v.y() * q.v.y(); + let zz = q.v.z() * q.v.z(); + let xy = q.v.x() * q.v.y(); + let xz = q.v.x() * q.v.z(); + let yz = q.v.y() * q.v.z(); + let wx = q.v.x() * q.w; + let wy = q.v.y() * q.w; + let wz = q.v.z() * q.w; + + let mut m = [[0.0; 4]; 4]; + + m[0][0] = 1. - 2. * (yy + zz); + m[0][1] = 2. * (xy - wz); + m[0][2] = 2. * (xz + wy); + + m[1][0] = 2. * (xy + wz); + m[1][1] = 1. - 2. * (xx + zz); + m[1][2] = 2. * (yz - wx); + + m[2][0] = 2. * (xz - wy); + m[2][1] = 2. * (yz + wx); + m[2][2] = 1. - 2. * (xx + yy); + + m[3][3] = 1.; + + let m_sq = SquareMatrix::new(m); + // For a pure rotation, the inverse is the transpose. + let m_inv_sq = m_sq.transpose(); + + Transform::new(&m_sq, &m_inv_sq) + } +} + +pub struct DerivativeTerm { + kc: Float, + kx: Float, + ky: Float, + kz: Float, +} + +impl DerivativeTerm { + pub fn new(kc: Float, kx: Float, ky: Float, kz: Float) -> Self { + Self { kc, kx, ky, kz } + } +} + +pub struct AnimatedTransform { + pub start_transform: Transform, + pub end_transform: Transform, + pub start_time: Float, + pub end_time: Float, + actually_animated: bool, + t: [Vector3f; 2], + r: [Quaternion; 2], + s: [SquareMatrix; 2], + has_rotation: bool, + c1: [DerivativeTerm; 3], + c2: [DerivativeTerm; 3], + c3: [DerivativeTerm; 3], + c4: [DerivativeTerm; 3], + c5: [DerivativeTerm; 3], +} + +impl AnimatedTransform { + pub fn interpolate(&self, time: Float) -> Transform { + if !self.actually_animated || time <= self.start_time { return self.start_transform } + if time >= self.end_time { return self.end_transform } + + let dt = (time - self.start_time) / (self.end_time - self.start_time); + let trans = (1.0 - dt) * self.t[0] + dt * self.t[1]; + + let rotate = Quaternion::slerp(dt, self.r[0], self.r[1]); + + let scale = (1.0 - dt) * self.s[0] + dt * self.s[1]; + + // Return interpolated matrix as product of interpolated components + return Transform::translate(trans) * Transform::from(rotate) * Transform::from_matrix(&scale); + } + + pub fn apply_inverse_point(&self, p: Point3f, time: Float) -> Point3f { + if !self.actually_animated { + return self.start_transform.apply_inverse(p); + } + return self.interpolate(time).apply_inverse(p) + } + + pub fn apply_inverse_vector(&self, v: Vector3f, time: Float) -> Vector3f { + if !self.actually_animated { + return self.start_transform.apply_inverse_vector(v); + } + return self.interpolate(time).apply_inverse_vector(v) + } +} + #[cfg(test)] mod tests { use super::*;