refactor: use 64-bit types instead of 32-bit

I don't know why I wasn't using 64 bit floats from the beginning,
honestly. I had weird priorities back then
This commit is contained in:
Kiana Sheibani 2024-10-14 18:06:50 -04:00
parent 9879184d47
commit 763a4ff923
Signed by: toki
GPG key ID: 6CB106C25E86A9F7
9 changed files with 469 additions and 243 deletions

View file

@ -2,11 +2,11 @@ extern crate nalgebra as na;
use std::cmp::Ordering;
use na::*;
use na::geometry::Point3;
use na::*;
use super::{bound::*, Surface};
use crate::types::*;
use super::{Surface, bound::*};
pub struct Triangle {
pub v1: usize, // Handles to 3 vertices.
@ -14,35 +14,45 @@ pub struct Triangle {
pub v3: usize,
normal: Unit3f, // Precalculated normal vector.
area: f32, // Precalculated area for barycentric calculations.
area: f64, // Precalculated area for barycentric calculations.
texture: Box<dyn Fn(f32, f32, f32) -> Texture> // Texture map.
// Uses barycentric coordinates as input.
texture: Box<dyn Fn(f64, f64, f64) -> Texture>, // Texture map.
// Uses barycentric coordinates as input.
}
pub struct TriangleMesh {
pub vertices: Vec<Point3f>,
pub triangles: Vec<Triangle>
pub triangles: Vec<Triangle>,
}
fn tri_area(a: &Point3f, b: &Point3f, c: &Point3f) -> f32 {
let prlg_area: f32 = (b - a).cross(&(c - a)).norm();
fn tri_area(a: &Point3f, b: &Point3f, c: &Point3f) -> f64 {
let prlg_area: f64 = (b - a).cross(&(c - a)).norm();
prlg_area / 2.0
}
impl Triangle {
fn vertex1<'a>(&self, vertices: &'a Vec<Point3f>) -> &'a Point3f { &vertices[self.v1] }
fn vertex2<'a>(&self, vertices: &'a Vec<Point3f>) -> &'a Point3f { &vertices[self.v2] }
fn vertex3<'a>(&self, vertices: &'a Vec<Point3f>) -> &'a Point3f { &vertices[self.v3] }
fn vertex1<'a>(&self, vertices: &'a Vec<Point3f>) -> &'a Point3f {
&vertices[self.v1]
}
fn vertex2<'a>(&self, vertices: &'a Vec<Point3f>) -> &'a Point3f {
&vertices[self.v2]
}
fn vertex3<'a>(&self, vertices: &'a Vec<Point3f>) -> &'a Point3f {
&vertices[self.v3]
}
// Conversion of barycentric coordinates to
// a point on the triangle.
fn from_bary(&self, vertices: &Vec<Point3f>, t: f32, u: f32, v: f32) -> Point3f {
Point::from(t * self.vertex1(vertices).coords + u * self.vertex2(vertices).coords + v * self.vertex3(vertices).coords)
fn from_bary(&self, vertices: &Vec<Point3f>, t: f64, u: f64, v: f64) -> Point3f {
Point::from(
t * self.vertex1(vertices).coords
+ u * self.vertex2(vertices).coords
+ v * self.vertex3(vertices).coords,
)
}
// Conversion of a point to barycentric coordinates.
fn to_bary(&self, vertices: &Vec<Point3f>, point: Point3f) -> (f32, f32, f32) {
fn to_bary(&self, vertices: &Vec<Point3f>, point: Point3f) -> (f64, f64, f64) {
let t = tri_area(self.vertex2(vertices), self.vertex3(vertices), &point) / self.area;
let u = tri_area(self.vertex1(vertices), self.vertex3(vertices), &point) / self.area;
let v = tri_area(self.vertex1(vertices), self.vertex2(vertices), &point) / self.area;
@ -50,32 +60,39 @@ impl Triangle {
(t, u, v)
}
fn intersect_(&self, vertices: &Vec<Point3f>, ray: Ray) -> Option<(f32, f32, f32)> {
fn intersect_(&self, vertices: &Vec<Point3f>, ray: Ray) -> Option<(f64, f64, f64)> {
let vect2_1 = self.vertex2(vertices) - self.vertex1(vertices);
let vect3_1 = self.vertex3(vertices) - self.vertex1(vertices);
let p_vect = ray.direction.cross(&vect3_1);
let det = p_vect.dot(&vect2_1);
if det.abs() < 1e-3 { return None; }
if det.abs() < 1e-3 {
return None;
}
let t_vect = ray.origin - self.vertex1(vertices);
let u = t_vect.dot(&p_vect) / det;
if u < 0.0 || u > 1.0 { return None; }
if u < 0.0 || u > 1.0 {
return None;
}
let q_vect = t_vect.cross(&vect2_1);
let v = ray.direction.dot(&q_vect) / det;
if v < 0.0 || (u + v) > 1.0 { return None; }
if v < 0.0 || (u + v) > 1.0 {
return None;
}
let t = 1.0 - u - v;
Some((t, u, v))
}
fn intersect(&self, vertices: &Vec<Point3f>, ray: Ray) -> Option<f32> {
self.intersect_(vertices, ray).map(|(t, u, v)| distance(&ray.origin, &self.from_bary(vertices, t, u, v)))
fn intersect(&self, vertices: &Vec<Point3f>, ray: Ray) -> Option<f64> {
self.intersect_(vertices, ray)
.map(|(t, u, v)| distance(&ray.origin, &self.from_bary(vertices, t, u, v)))
}
fn get_texture(&self, vertices: &Vec<Point3f>, point: Point3f) -> Texture {
@ -86,45 +103,81 @@ impl Triangle {
#[allow(dead_code)]
impl TriangleMesh {
pub fn new(vertices: Vec<Point3f>, tris: Vec<(usize, usize, usize, Box<dyn Fn(f32, f32, f32) -> Texture>)>) -> Self {
let triangles = tris.into_iter()
.map(|(v1, v2, v3, f)| Triangle {
v1, v2, v3,
normal: Unit::new_normalize((&vertices[v2] - &vertices[v1]).cross(&(&vertices[v3] - &vertices[v1]))),
area: tri_area(&vertices[v1], &vertices[v2], &vertices[v3]),
texture: f
}).collect();
TriangleMesh {
vertices, triangles
}
}
pub fn new_solid(vertices: Vec<Point3f>, tris: Vec<(usize, usize, usize)>, texture: Texture) -> Self {
let triangles = tris.into_iter()
.map(|(v1, v2, v3)| Triangle {
v1, v2, v3,
normal: Unit::new_normalize((&vertices[v2] - &vertices[v1]).cross(&(&vertices[v3] - &vertices[v1]))),
area: tri_area(&vertices[v1], &vertices[v2], &vertices[v3]),
texture: Box::new(move |_, _, _| texture)
}).collect();
pub fn new(
vertices: Vec<Point3f>,
tris: Vec<(usize, usize, usize, Box<dyn Fn(f64, f64, f64) -> Texture>)>,
) -> Self {
let triangles = tris
.into_iter()
.map(|(v1, v2, v3, f)| Triangle {
v1,
v2,
v3,
normal: Unit::new_normalize(
(&vertices[v2] - &vertices[v1]).cross(&(&vertices[v3] - &vertices[v1])),
),
area: tri_area(&vertices[v1], &vertices[v2], &vertices[v3]),
texture: f,
})
.collect();
TriangleMesh {
vertices,
triangles
triangles,
}
}
pub fn singleton<F: 'static>(vertex1: Point3f, vertex2: Point3f, vertex3: Point3f, texture: F) -> Self
where F: Fn(f32, f32, f32) -> Texture
{ TriangleMesh::new(vec![vertex1, vertex2, vertex3], vec![(0, 1, 2, Box::new(texture))]) }
pub fn new_solid(
vertices: Vec<Point3f>,
tris: Vec<(usize, usize, usize)>,
texture: Texture,
) -> Self {
let triangles = tris
.into_iter()
.map(|(v1, v2, v3)| Triangle {
v1,
v2,
v3,
normal: Unit::new_normalize(
(&vertices[v2] - &vertices[v1]).cross(&(&vertices[v3] - &vertices[v1])),
),
area: tri_area(&vertices[v1], &vertices[v2], &vertices[v3]),
texture: Box::new(move |_, _, _| texture),
})
.collect();
TriangleMesh {
vertices,
triangles,
}
}
pub fn singleton_solid(vertex1: Point3f, vertex2: Point3f, vertex3: Point3f, texture: Texture) -> Self
{ TriangleMesh::singleton(vertex1, vertex2, vertex3, move |_, _, _| texture) }
pub fn singleton<F: 'static>(
vertex1: Point3f,
vertex2: Point3f,
vertex3: Point3f,
texture: F,
) -> Self
where
F: Fn(f64, f64, f64) -> Texture,
{
TriangleMesh::new(
vec![vertex1, vertex2, vertex3],
vec![(0, 1, 2, Box::new(texture))],
)
}
pub fn singleton_solid(
vertex1: Point3f,
vertex2: Point3f,
vertex3: Point3f,
texture: Texture,
) -> Self {
TriangleMesh::singleton(vertex1, vertex2, vertex3, move |_, _, _| texture)
}
fn closest_tri(&self, point: Point3f) -> &Triangle {
self.triangles.iter()
self.triangles
.iter()
.map(move |tri| {
let rel_pos = point - tri.vertex1(&self.vertices);
let proj_point3 = rel_pos - (*tri.normal * tri.normal.dot(&rel_pos));
@ -139,15 +192,17 @@ impl TriangleMesh {
(tri, distance(&point, &point_new))
})
.min_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(Ordering::Equal))
.unwrap().0
.unwrap()
.0
}
}
impl Surface for TriangleMesh {
fn intersect(&self, ray: Ray) -> Option<f32> {
self.triangles.iter()
.filter_map(|tri| tri.intersect(&self.vertices, ray))
.min_by(|a, b| a.partial_cmp(&b).unwrap_or(Ordering::Equal))
fn intersect(&self, ray: Ray) -> Option<f64> {
self.triangles
.iter()
.filter_map(|tri| tri.intersect(&self.vertices, ray))
.min_by(|a, b| a.partial_cmp(&b).unwrap_or(Ordering::Equal))
}
fn normal(&self, point: Point3f) -> Unit3f {
@ -160,22 +215,26 @@ impl Surface for TriangleMesh {
// Uses Welzl's algorithm to solve the bounding sphere problem
fn bound(&self) -> Bound {
fn smallest_sphere_plane(points: Vec<&Point3f>, boundary: Vec<&Point3f>) -> (Point3f, f32) {
fn smallest_sphere_plane(points: Vec<&Point3f>, boundary: Vec<&Point3f>) -> (Point3f, f64) {
if points.len() == 0 || boundary.len() == 3 {
match boundary.len() {
0 => (Point3::new(0.0, 0.0, 0.0), 0.0),
1 => (*boundary[0], 0.0),
2 => { let half_span = 0.5 * (boundary[1] - boundary[0]);
(*boundary[0] + half_span, half_span.norm()) },
2 => {
let half_span = 0.5 * (boundary[1] - boundary[0]);
(*boundary[0] + half_span, half_span.norm())
}
3 => triangle_sphere(boundary[0], boundary[1], boundary[2]),
_ => unreachable!()
_ => unreachable!(),
}
} else {
let removed = points[0];
let points = Vec::from(&points[1..]);
let bound = smallest_sphere(points.clone(), boundary.clone());
if distance(&bound.0, removed) < bound.1 { return bound; }
if distance(&bound.0, removed) < bound.1 {
return bound;
}
let mut boundary = boundary.clone();
boundary.push(removed);
@ -184,32 +243,44 @@ impl Surface for TriangleMesh {
}
}
fn triangle_sphere(point1: &Point3f, point2: &Point3f, point3: &Point3f) -> (Point3f, f32) {
fn triangle_sphere(point1: &Point3f, point2: &Point3f, point3: &Point3f) -> (Point3f, f64) {
let a = point3 - point1;
let b = point2 - point1;
let crs = b.cross(&a);
let to_center = (crs.cross(&b) * a.norm_squared() + a.cross(&crs) * b.norm_squared())
/ (2.0 * crs.norm_squared());
/ (2.0 * crs.norm_squared());
let radius = to_center.norm();
(point1 + to_center, radius)
}
fn tetrahedron_sphere(point1: &Point3f, point2: &Point3f, point3: &Point3f, point4: &Point3f) -> (Point3f, f32) {
let matrix = Matrix4::from_rows(&[point1.to_homogeneous().transpose(),
point2.to_homogeneous().transpose(),
point3.to_homogeneous().transpose(),
point4.to_homogeneous().transpose()]);
fn tetrahedron_sphere(
point1: &Point3f,
point2: &Point3f,
point3: &Point3f,
point4: &Point3f,
) -> (Point3f, f64) {
let matrix = Matrix4::from_rows(&[
point1.to_homogeneous().transpose(),
point2.to_homogeneous().transpose(),
point3.to_homogeneous().transpose(),
point4.to_homogeneous().transpose(),
]);
let a = matrix.determinant() * 2.0;
if (a != 0.0) {
let mut matrix_mut = matrix.clone();
let squares = Vector4::new(point1.coords.norm_squared(), point2.coords.norm_squared(), point3.coords.norm_squared(), point4.coords.norm_squared());
let squares = Vector4::new(
point1.coords.norm_squared(),
point2.coords.norm_squared(),
point3.coords.norm_squared(),
point4.coords.norm_squared(),
);
matrix_mut.set_column(0, &squares);
let center_x = matrix_mut.determinant();
@ -231,23 +302,27 @@ impl Surface for TriangleMesh {
}
}
fn smallest_sphere(points: Vec<&Point3f>, boundary: Vec<&Point3f>) -> (Point3f, f32) {
fn smallest_sphere(points: Vec<&Point3f>, boundary: Vec<&Point3f>) -> (Point3f, f64) {
if points.len() == 0 || boundary.len() == 4 {
match boundary.len() {
0 => (Point3::new(0.0, 0.0, 0.0), 0.0),
1 => (*boundary[0], 0.0),
2 => { let half_span = 0.5 * (boundary[1] - boundary[0]);
(*boundary[0] + half_span, half_span.norm()) },
2 => {
let half_span = 0.5 * (boundary[1] - boundary[0]);
(*boundary[0] + half_span, half_span.norm())
}
3 => triangle_sphere(boundary[0], boundary[1], boundary[2]),
4 => tetrahedron_sphere(boundary[0], boundary[1], boundary[2], boundary[3]),
_ => unreachable!()
_ => unreachable!(),
}
} else {
let removed = points[0];
let points = Vec::from(&points[1..]);
let bound = smallest_sphere(points.clone(), boundary.clone());
if distance(&bound.0, removed) < bound.1 { return bound; }
if distance(&bound.0, removed) < bound.1 {
return bound;
}
let mut boundary = boundary.clone();
boundary.push(removed);
@ -257,14 +332,18 @@ impl Surface for TriangleMesh {
}
extern crate rand;
use rand::thread_rng;
use rand::seq::SliceRandom;
use rand::thread_rng;
let mut points: Vec<&Point3f> = self.vertices.iter().collect();
points.shuffle(&mut thread_rng());
let (center, radius) = smallest_sphere(points, Vec::new());
Bound { center, radius: radius + 1e-3, bypass: false }
Bound {
center,
radius: radius + 1e-3,
bypass: false,
}
}
}