Support cylinders

This commit is contained in:
Michael Zhang 2023-02-02 01:43:54 -06:00
parent 0ad8932189
commit adf8fa87fa
4 changed files with 124 additions and 17 deletions

View file

@ -1,7 +1,7 @@
imsize 640 480
eye 0 0 3
eye 0 0 15
viewdir 0 0 -1
hfov 120
hfov 60
updir 0 1 0
bkgcolor 0.1 0.1 0.1
@ -18,4 +18,4 @@ sphere 5 5 -1 1
sphere -6 -4 -8 7
mtlcolor 0.5 1 0.5
cylinder 2 1 -2 1 -2 1 1 2
cylinder 5 1 -2 1 -2 1 1 2

View file

@ -90,7 +90,6 @@ fn main() -> Result<()> {
let pixel_in_space = translate_pixel(px, py);
let ray_start = if scene.parallel_projection {
println!("Using parallel projection.");
// For a parallel projection, we'll just take the view direction and
// subtract it from the target point. This means every single
// ray will be viewed from a point at infinity, rather than a single eye

View file

@ -1,4 +1,5 @@
use nalgebra::Vector3;
use nalgebra::{Matrix3, Vector3};
use ordered_float::NotNan;
use crate::scene_data::{Cylinder, Sphere};
@ -73,19 +74,122 @@ impl Ray {
pub fn intersects_cylinder_at(&self, cylinder: &Cylinder) -> Option<f64> {
// Determine rotation matrix for turning the cylinder upright along the
// Z-axis
let rotation_matrix = {
let target_direction = Vector3::new(0.0, 0.0, 1.0);
let axis = target_direction.cross(&cylinder.direction).normalize();
let angle = (target_direction.dot(&cylinder.direction)
/ (target_direction.norm() * cylinder.direction.norm()))
.acos();
let target_direction = Vector3::new(0.0, 0.0, 1.0);
let rotation_matrix =
compute_rotation_matrix(cylinder.direction, target_direction);
// Transform all parameters according to this rotation matrix
let rotated_cylinder_center = rotation_matrix * cylinder.center;
let rotated_ray_origin = rotation_matrix * self.origin;
let rotated_ray_direction = rotation_matrix * self.direction;
// Now that we know the cylinder is upright, we can start checking against
// the formula:
//
// (ox + t*rx - cx)^2 + (oy + t*ry - cy)^2 = r^2
//
// where o{xy} is the ray origin, r{xy} is the ray direction, and c{xy} is
// the cylinder center. The z will be taken care of after the fact. To
// solve, we must put it into the form At^2 + Bt + c = 0. The variables
// are:
//
// A: rx^2 + ry^2
// B: 2(rx(ox - cx) + ry(oy - cy))
// C: (cx - ox)^2 + (cy - oy)^2 - r^2
let (a, b, c) = {
let o = rotated_ray_origin;
let r = rotated_ray_direction;
let c = rotated_cylinder_center;
(
r.x.powi(2) + r.y.powi(2),
2.0 * (r.x * (o.x - c.x) + r.y * (o.y - c.y)),
(c.x - o.x).powi(2) + (c.y - o.y).powi(2) - cylinder.radius.powi(2),
)
};
// TODO: Implement
None
let discriminant = b * b - 4.0 * a * c;
let mut solutions = match discriminant {
// Discriminant < 0, means the equation has no solutions.
d if d < 0.0 => vec![],
// Discriminant == 0
d if d == 0.0 => vec![-b / 2.0 * a],
// Discriminant > 0, 2 solutions available.
d if d > 0.0 => {
vec![
(-b + discriminant.sqrt()) / (2.0 * a),
(-b - discriminant.sqrt()) / (2.0 * a),
]
}
// Probably hit some NaN or Infinity value due to faulty inputs...
_ => unreachable!("Invalid determinant value: {discriminant}"),
};
// We also need to add solutions for the two ends of the cylinder, which
// uses a similar method except backwards: check intersection points
// with the correct z-plane and then see if the points are within the
// circle.
//
// Luckily, this means we only need to care about one dimension at first,
// and don't need to perform the quadratic equation method above.
//
// oz + t * rz = cz +- (len / 2)
// t = (oz + cz +- (len / 2)) / rz
let possible_z_intersections = {
let o = rotated_ray_origin;
let r = rotated_ray_direction;
let c = rotated_cylinder_center;
vec![
(o.z + c.z + cylinder.length / 2.0) / r.z,
(o.z + c.z - cylinder.length / 2.0) / r.z,
]
};
// Filter out all the solutions where the z does not lie in the circle
solutions.extend(possible_z_intersections.into_iter().filter(|t| {
let ray_point = self.origin + self.direction * (*t);
ray_point.x.powi(2) + ray_point.y.powi(2) <= cylinder.radius.powi(2)
}));
// Filter out solutions that don't have a valid Z position.
let solutions = solutions
.into_iter()
.filter(|t| {
let ray_point = self.origin + self.direction * (*t);
let rotated_ray_point = rotation_matrix * ray_point;
let z = rotated_ray_point.z - rotated_cylinder_center.z;
// Check to see if z is between -len/2 and len/2
z.abs() < cylinder.length / 2.0
})
.filter_map(|t| NotNan::new(t).ok());
// Return the minimum solution
solutions.min().map(|t| t.into_inner())
}
}
/// Calculate the rotation matrix between the 2 given vectors
/// Based on the method here: https://math.stackexchange.com/a/897677
fn compute_rotation_matrix(a: Vector3<f64>, b: Vector3<f64>) -> Matrix3<f64> {
let cos_t = a.dot(&b);
let sin_t = a.cross(&b).norm();
let G = Matrix3::new(cos_t, -sin_t, 0.0, sin_t, cos_t, 0.0, 0.0, 0.0, 1.0);
// Basis
let u = a;
let v = (b - a.dot(&b) * a).normalize();
let w = b.cross(&a);
let F = Matrix3::from_columns(&[u, v, w]).try_inverse().unwrap();
F.try_inverse().unwrap() * G * F
}
#[cfg(test)]
mod tests {
use nalgebra::Vector3;

View file

@ -74,9 +74,9 @@ $\frac{\Delta x + \Delta y}{2}$ to the point to get that)
### Parallel Projection Notes
Because of the way I implemented parallel projection, it's recommended to use
`--distance` to force a much bigger distance from the eye for the raycaster.
This is due to the size of the image.
Because of the way I implemented parallel projection, it's recommended to
either put the eye much farther back, or use `--distance` to force a much bigger
distance from the eye for the raycaster. This is due to the size of the image.
### Cylinder Intersection Notes
@ -85,4 +85,8 @@ cylinder, so that the cylinder location is $(0, 0, 0)$ and the direction vector
is normalized into $(0, 0, 1)$.
Then it's a matter of determining if the $x$ and $y$ coordinates fall into the
space constrained by the equation $x^2 + y^2 = r^2$ and if $z \le L$.
space constrained by the equation $(ox + t\times rx - cx)^2 + (oy + t\times ty -
cy)^2 = r^2$ and if $z \le L$.
See the comments in the code for a more detailed explanation of the equation
used.