Added transformations, quaternion logic, spectra, and beginning of camera

This commit is contained in:
pingupingou 2025-11-03 02:14:10 +00:00
parent 87b3f8d6dd
commit 45e961a1a0
15 changed files with 2482 additions and 421 deletions

1
.gitignore vendored
View file

@ -1,3 +1,4 @@
/target
*.lock
*.log
*.bak

View file

@ -5,3 +5,4 @@ edition = "2024"
[dependencies]
num-traits = "0.2.19"
once_cell = "1.21.3"

264
README.md
View file

@ -1,264 +0,0 @@
<!-- Improved compatibility of back to top link: See: https://github.com/othneildrew/Best-README-Template/pull/73 -->
<a id="readme-top"></a>
[![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]
<!-- PROJECT LOGO -->
<br />
<div align="center">
<a href="https://github.com/othneildrew/Best-README-Template">
<img src="images/logo.png" alt="Logo" width="80" height="80">
</a>
<h3 align="center">Best-README-Template</h3>
<p align="center">
An awesome README template to jumpstart your projects!
<br />
<a href="https://github.com/othneildrew/Best-README-Template"><strong>Explore the docs »</strong></a>
<br />
<br />
<a href="https://github.com/othneildrew/Best-README-Template">View Demo</a>
&middot;
<a href="https://github.com/othneildrew/Best-README-Template/issues/new?labels=bug&template=bug-report---.md">Report Bug</a>
&middot;
<a href="https://github.com/othneildrew/Best-README-Template/issues/new?labels=enhancement&template=feature-request---.md">Request Feature</a>
</p>
</div>
<!-- TABLE OF CONTENTS -->
<details>
<summary>Table of Contents</summary>
<ol>
<li>
<a href="#about-the-project">About The Project</a>
<ul>
<li><a href="#built-with">Built With</a></li>
</ul>
</li>
<li>
<a href="#getting-started">Getting Started</a>
<ul>
<li><a href="#prerequisites">Prerequisites</a></li>
<li><a href="#installation">Installation</a></li>
</ul>
</li>
<li><a href="#usage">Usage</a></li>
<li><a href="#roadmap">Roadmap</a></li>
<li><a href="#contributing">Contributing</a></li>
<li><a href="#license">License</a></li>
<li><a href="#contact">Contact</a></li>
<li><a href="#acknowledgments">Acknowledgments</a></li>
</ol>
</details>
<!-- ABOUT THE PROJECT -->
## 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.
<p align="right">(<a href="#readme-top">back to top</a>)</p>
### 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]
<p align="right">(<a href="#readme-top">back to top</a>)</p>
<!-- GETTING STARTED -->
## 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
```
<p align="right">(<a href="#readme-top">back to top</a>)</p>
<!-- USAGE EXAMPLES -->
## 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)_
<p align="right">(<a href="#readme-top">back to top</a>)</p>
<!-- ROADMAP -->
## 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).
<p align="right">(<a href="#readme-top">back to top</a>)</p>
<!-- CONTRIBUTING -->
## 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:
<a href="https://github.com/othneildrew/Best-README-Template/graphs/contributors">
<img src="https://contrib.rocks/image?repo=othneildrew/Best-README-Template" alt="contrib.rocks image" />
</a>
<p align="right">(<a href="#readme-top">back to top</a>)</p>
<!-- LICENSE -->
## License
Distributed under the Unlicense License. See `LICENSE.txt` for more information.
<p align="right">(<a href="#readme-top">back to top</a>)</p>
<!-- CONTACT -->
## 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)
<p align="right">(<a href="#readme-top">back to top</a>)</p>
<!-- ACKNOWLEDGMENTS -->
## 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)
<p align="right">(<a href="#readme-top">back to top</a>)</p>
<!-- MARKDOWN LINKS & IMAGES -->
<!-- https://www.markdownguide.org/basic-syntax/#reference-style-links -->
[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

76
src/core/camera.rs Normal file
View file

@ -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<Float>,
}
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;

496
src/core/color.rs Normal file
View file

@ -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<usize> 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<usize> 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<XYZ> 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<Float> 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<Float> for XYZ {
fn mul_assign(&mut self, rhs: Float) {
self.x *= rhs;
self.y *= rhs;
self.z *= rhs;
}
}
impl Mul<XYZ> 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<Float> 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<Float> 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<usize> 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<usize> 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<RGB> 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<Float> 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<Float> for RGB {
fn mul_assign(&mut self, rhs: Float) {
self.r *= rhs;
self.g *= rhs;
self.b *= rhs;
}
}
impl Mul<RGB> 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<Float> 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<Float> 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] }
}
}

75
src/core/colorspace.rs Normal file
View file

@ -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<RGBToSpectrumTable>,
xyz_from_rgb: SquareMatrix<Float, 3>,
rgb_from_xyz: SquareMatrix<Float, 3>,
}
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<Float, 3> = 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<Float, 3> {
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)
}
}

9
src/core/film.rs Normal file
View file

@ -0,0 +1,9 @@
pub enum Film {
RGB(RGBFilm),
GBuffer(GBufferBFilm),
Spectral(SpectralFilm),
}
pub struct RGBFilm;
pub struct GBufferBFilm;
pub struct SpectralFilm;

View file

@ -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<T, const N: usize>:
@ -37,9 +38,11 @@ macro_rules! impl_tuple_core {
}
impl<T: Default + Copy, const N: usize> Default for $Struct<T, N> {
fn default() -> Self {
Self([T::default(); N])
fn default() -> Self { Self([T::default(); N]) }
}
impl<T, const N: usize> $Struct<T, N> where T: Zero + Copy {
#[inline] pub fn zero() -> Self { Self([T::zero(); N]) }
}
impl<T, const N: usize> Index<usize> for $Struct<T, N> {
@ -58,7 +61,11 @@ macro_rules! impl_tuple_core {
Self(result)
}
}
};
}
macro_rules! impl_scalar_ops {
($Struct:ident) => {
impl<T, const N: usize> Mul<T> for $Struct<T, N> where T: Mul<Output = T> + Copy {
type Output = Self;
fn mul(self, rhs: T) -> Self::Output {
@ -67,6 +74,14 @@ macro_rules! impl_tuple_core {
Self(result)
}
}
impl<const N: usize> Mul<$Struct<Float, N>> for Float {
type Output = $Struct<Float, N>;
fn mul(self, rhs: $Struct<Float, N>) -> Self::Output {
rhs * self
}
}
impl<T, const N: usize> MulAssign<T> for $Struct<T, N> 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<T, const N: usize> Add for $Struct<T, N> where T: Add<Output = T> + 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<T, const N: usize> $Op<$Rhs<T, N>> for $Lhs<T, N>
where
T: $Op<Output = T> + Copy,
{
type Output = $Output<T, N>;
fn $op(self, rhs: $Rhs<T, N>) -> Self::Output {
let mut result = self.0;
for i in 0..N { result[i] = self.0[i] + rhs.0[i]; }
Self(result)
for i in 0..N {
result[i] = self.0[i].$op(rhs.0[i]);
}
}
impl<T, const N: usize> AddAssign for $Struct<T, N> where T: AddAssign + Copy {
fn add_assign(&mut self, rhs: Self) {
for i in 0..N { self.0[i] += rhs.0[i]; }
}
}
impl<T, const N: usize> Sub for $Struct<T, N> where T: Sub<Output = T> + 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<T, const N: usize> SubAssign for $Struct<T, N> where T: SubAssign + Copy {
fn sub_assign(&mut self, rhs: Self) {
for i in 0..N { self.0[i] -= 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<T, const N: usize> $OpAssign<$Rhs<T, N>> for $Lhs<T, N>
where
T: $OpAssign + Copy,
{
fn $op_assign(&mut self, rhs: $Rhs<T, N>) {
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<Rhs = Self> {
type Output;
fn dot(self, rhs: Rhs) -> Self::Output;
}
macro_rules! impl_dot_for {
($Lhs:ident, $Rhs:ident) => {
impl<T: Num + Copy, const N: usize> Dot<$Rhs<T, N>> for $Lhs<T, N> {
type Output = T;
fn dot(self, rhs: $Rhs<T, N>) -> 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<T: Copy, const N: usize> From<Vector<T, N>> for Normal<T, N> {
fn from(v: Vector<T, N>) -> Self { Self(v.0) }
@ -134,45 +199,6 @@ impl<T: Copy, const N: usize> From<Normal<T, N>> for Vector<T, N> {
fn from(n: Normal<T, N>) -> Self { Self(n.0) }
}
impl<T, const N: usize> Sub for Point<T, N> where T: Sub<Output = T> + Copy {
type Output = Vector<T, N>;
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<T, const N: usize> Add<Vector<T, N>> for Point<T, N> where T: Add<Output = T> + Copy {
type Output = Self;
fn add(self, rhs: Vector<T, N>) -> Self::Output {
let mut result = self.0;
for i in 0..N { result[i] = self.0[i] + rhs.0[i]; }
Self(result)
}
}
impl<T, const N: usize> AddAssign<Vector<T, N>> for Point<T, N> where T: AddAssign + Copy {
fn add_assign(&mut self, rhs: Vector<T, N>) {
for i in 0..N { self.0[i] += rhs.0[i]; }
}
}
// Point - Vector -> Point
impl<T, const N: usize> Sub<Vector<T, N>> for Point<T, N> where T: Sub<Output = T> + Copy {
type Output = Self;
fn sub(self, rhs: Vector<T, N>) -> Self::Output {
let mut result = self.0;
for i in 0..N { result[i] = self.0[i] - rhs.0[i]; }
Self(result)
}
}
impl<T, const N: usize> SubAssign<Vector<T, N>> for Point<T, N> where T: SubAssign + Copy {
fn sub_assign(&mut self, rhs: Vector<T, N>) {
for i in 0..N { self.0[i] -= rhs.0[i]; }
}
}
macro_rules! impl_accessors {
($Struct:ident) => {
impl<T: Copy> $Struct<T, 2> {
@ -193,11 +219,6 @@ impl_accessors!(Normal);
// Vector stuff
pub trait Dot<Rhs = Self> {
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, Output=Self>, Self::Scalar: NumFloat;
}
macro_rules! impl_vector_math_for {
macro_rules! impl_normed_for {
($Struct:ident) => {
impl<T: Num + Copy, const N: usize> Dot for $Struct<T, N> {
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<T: Num + Copy, const N: usize> Normed for $Struct<T, N> {
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<T, const N: usize> $Struct<T, N>
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<T, const N: usize> Point<T, N>
where
T: NumFloat + Copy,
Point<T, N>: Sub<Output = Vector<T, N>>, // Point - Point -> Vector
Vector<T, N>: Normed<Scalar = T>, // Vector has a norm
Point<T, N>: Sub<Output = Vector<T, N>>,
Vector<T, N>: Normed<Scalar = T>,
{
pub fn distance(self, other: Self) -> T {
(self - other).norm()
@ -241,19 +274,25 @@ where
}
}
// Utility aliases and functions
pub type Point2<T> = Point<T, 2>;
pub type Point2f = Point2<Float>;
pub type Point2i = Point2<i32>;
pub type Point3<T> = Point<T, 3>;
pub type Point3f = Point3<Float>;
pub type Point3i = Point3<i32>;
pub type Vector2<T> = Vector<T, 2>;
pub type Vector2f = Vector2<Float>;
pub type Vector2i = Vector2<i32>;
pub type Vector3<T> = Vector<T, 3>;
pub type Vector3f = Vector3<Float>;
pub type Vector3i = Vector3<i32>;
pub type Normal3<T> = Normal<T, 3>;
pub type Normal3f = Normal3<Float>;
pub type Normal3i = Normal3<i32>;
pub type Vector3fi = Vector<Interval, 3>;
pub type Point3fi = Point<Interval, 3>;
impl<T: Copy> Vector2<T> { pub fn new(x: T, y: T) -> Self { Self([x, y]) } }
impl<T: Copy> Point2<T> { pub fn new(x: T, y: T) -> Self { Self([x, y]) } }
@ -288,11 +327,97 @@ where T: Num + NumFloat + Copy + Neg<Output = T>
}
}
impl<const N: usize> Point<Interval, N> {
pub fn new_from_point(p: Point<Float, N>) -> 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<Float, N>, e: Vector<Float, N>) -> 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<Float, N> {
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<Float, N> {
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<const N: usize> Vector<Interval, N> {
pub fn new_from_vector(v: Vector<Float, N>) -> 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<Float, N>, e: Vector<Float, N>) -> 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<Float, N> {
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<Float, N> {
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<const N: usize> From<Point<Interval, N>> for Point<Float, N> {
fn from(pi: Point<Interval, N>) -> Self {
let mut arr = [0.0; N];
for i in 0..N {
arr[i] = pi[i].midpoint();
}
Point(arr)
}
}
impl<T> Normal3<T>
where T: Num + PartialOrd + Copy + Neg<Output = T>
{
pub fn face_forward(self, v: Vector3<T>) -> 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<T: NumFloat>(v: Vector3f) -> Float {
}
}
// AABB bouding boxes
// AABB BOUNDING BOXES
#[derive(Debug, Copy, Clone, PartialEq)]
pub struct Bounds<T, const N: usize> {
pub p_min: Point<T, N>,
@ -469,7 +594,7 @@ where
<<Point3<T> 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;

View file

@ -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::<SurfaceInteraction>()
}
fn is_medium_interaction(&self) -> bool {
self.as_any().is::<MediumInteraction>()
}
/// Determines the medium for a ray starting at this interaction.
fn get_medium(&self, w: Vector3f) -> Option<Arc<Medium>>;
/// 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<Arc<MediumInterface>>,
pub face_index: bool,
pub area_light: Option<Arc<Light>>,
pub dpdx: Vector3f,
pub dpdy: Vector3f,
pub dudx: Float,
pub dvdx: Float,
pub dudy: Float,
pub dvdy: Float,
// pub shape: Option<Arc<dyn Shape>>,
// pub bsdf: Option<BSDF>,
}
pub struct PhaseFunction;
pub struct MediumInteraction {
pub common: InteractionData,
pub medium: Arc<Medium>,
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<Arc<Medium>> {
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::<SurfaceInteraction>() {
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<Arc<Medium>> {
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,
}
}
}

102
src/core/interval.rs Normal file
View file

@ -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)
}
}

View file

@ -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;

View file

@ -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<Float> {
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<P>(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,
}

View file

@ -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<usize> 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<usize> 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;
}
}

View file

@ -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<Float>) -> 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::<Float>() / (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<usize> for SampledSpectrum {
type Output = Float;
fn index(&self, i: usize) -> &Self::Output { &self.values[i] }
}
impl IndexMut<usize> 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<SampledSpectrum> 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<Float> 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<Float> 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<Float> 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<Float> 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<usize> for SampledWavelengths {
type Output = Float;
fn index(&self, i: usize) -> &Self::Output { &self.lambda[i] }
}
impl IndexMut<usize> 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<Arc<DenselySampledSpectrum>>,
}
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<Float>,
}
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: 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<Spectrum> = Lazy::new(|| {
let dss = DenselySampledSpectrum::from_piecewise_linear(&CIE_X_DATA);
Spectrum::DenselySampled(dss)
});
pub static Y: Lazy<Spectrum> = Lazy::new(|| {
let dss = DenselySampledSpectrum::from_piecewise_linear(&CIE_Y_DATA);
Spectrum::DenselySampled(dss)
});
pub static Z: Lazy<Spectrum> = Lazy::new(|| {
let dss = DenselySampledSpectrum::from_piecewise_linear(&CIE_Z_DATA);
Spectrum::DenselySampled(dss)
});
}

View file

@ -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<T, const N: usize> {
m: [[T; N]; N],
pub m: [[T; N]; N],
}
impl<T, const N: usize> SquareMatrix<T, N>
{
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<T, const N: usize> SquareMatrix<T, N>
Self { m }
}
pub fn zero() -> Self where T: Copy + Zero {
SquareMatrix { m: [[T::zero(); N]; N] }
}
@ -43,6 +47,21 @@ impl<T, const N: usize> SquareMatrix<T, N>
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<T, const N: usize> SquareMatrix<T, N>
{
pub fn inverse(&self) -> Option<Self> {
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,15 +128,42 @@ impl<T, const N: usize> SquareMatrix<T, N>
}
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];
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 {
// 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() {
@ -125,19 +171,17 @@ impl<T, const N: usize> SquareMatrix<T, N>
}
}
// 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).
// Singular matrix
if lum[i][i] == T::zero() {
return T::zero();
}
// Gaussian elimination to form the upper-triangular matrix
// Gaussian elimination
for j in (i + 1)..N {
let factor = lum[j][i] / lum[i][i];
for k in i..N {
@ -158,6 +202,8 @@ impl<T, const N: usize> SquareMatrix<T, N>
det
}
}
}
}
impl<T, const N: usize> Default for SquareMatrix<T, N>
where
@ -272,6 +318,29 @@ where
}
}
impl<T, const N: usize> Mul<T> for SquareMatrix<T, N>
where
T: Copy + Mul<Output = T>,
{
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<const N: usize> Mul<SquareMatrix<Float, N>> for Float
{
type Output = SquareMatrix<Float, N>;
fn mul(self, rhs: SquareMatrix<Float, N>) -> Self::Output {
rhs * self
}
}
impl<T, const N: usize> Div<T> for SquareMatrix<T, N>
where
T: Copy + Div<Output = T>,
@ -320,6 +389,61 @@ where
}
}
impl SquareMatrix<Float, 3> {
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<XYZ> for SquareMatrix<Float, 3>
{
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<RGB> for SquareMatrix<Float, 3>
{
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<T: num_traits::Float> {
m: SquareMatrix<T, 4>,
@ -327,7 +451,11 @@ pub struct Transform<T: num_traits::Float> {
}
impl<T: num_traits::Float + Signed> Transform<T> {
pub fn new(m: &SquareMatrix<T, 4>) -> Self {
pub fn new(m: &SquareMatrix<T, 4>, m_inv: &SquareMatrix<T, 4>) -> Self {
Self { m: *m, m_inv: *m_inv }
}
pub fn from_matrix(m: &SquareMatrix<T, 4>) -> Self {
let inv = match m.inverse() {
Some(inv) => inv ,
None => SquareMatrix::zero(),
@ -336,6 +464,12 @@ impl<T: num_traits::Float + Signed> Transform<T> {
Self { m: *m, m_inv: inv }
}
pub fn identity() -> Self {
let m: SquareMatrix<T, 4> = 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<T: num_traits::Float + Signed> Transform<T> {
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<T, 3>) -> Point<T, 3> {
*self * p
}
pub fn apply_inverse_vector(&self, v: Vector<T, 3>) -> Vector<T, 3> {
*self * v
}
}
impl Transform<Float> {
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<T: num_traits::Float> PartialEq for Transform<T> {
}
}
impl<T> Mul for Transform<T>
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<Transform<Float>> for Float
{
type Output = Transform<Float>;
fn mul(self, rhs: Transform<Float>) -> Self::Output {
Transform { m: rhs.m * self , m_inv: rhs.m_inv * self }
}
}
impl<T> Mul<Point<T, 3>> for Transform<T>
where
T: NumFloat
{
type Output = Point<T, 3>;
fn mul(self, p: Point<T, 3>) -> 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<T> Mul<Vector<T, 3>> for Transform<T>
where
T: NumFloat
{
type Output = Vector<T, 3>;
fn mul(self, p: Vector<T, 3>) -> 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<Quaternion> for Transform<Float> {
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<Float>,
pub end_transform: Transform<Float>,
pub start_time: Float,
pub end_time: Float,
actually_animated: bool,
t: [Vector3f; 2],
r: [Quaternion; 2],
s: [SquareMatrix<Float, 4>; 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<Float> {
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::*;