pbrt/shared/src/shapes/curves.rs
Wito Wiala 731a37abae More fixes to allow for spirv compilation. Namely,
issues with array initialization with functions ([object::default()]),
initializing arrays from a function (implemented a simple gpu version
using manual MaybeUninit pointers),
changing enums with distinct types to structs or changing the selection
logic, changing pointer subtraction in light samplers to a scan (this
will come back to bite me in the ass), and ignoring the data module,
since SPIR-V cant use pointers in statics.
2026-02-20 21:11:17 +00:00

391 lines
13 KiB
Rust

use crate::Float;
use crate::core::geometry::{
Bounds3f, DirectionCone, Normal3f, Point2f, Point3f, Point3fi, Ray, Vector2f, Vector3f,
VectorLike,
};
use crate::core::interaction::{Interaction, InteractionTrait, SurfaceInteraction};
use crate::core::shape::{ShapeIntersection, ShapeSample, ShapeSampleContext, ShapeTrait};
use crate::utils::gpu_array_from_fn;
use crate::utils::math::{clamp, lerp, square};
use crate::utils::splines::{
bound_cubic_bezier, cubic_bezier_control_points, evaluate_cubic_bezier, subdivide_cubic_bezier,
};
use crate::utils::transform::{Transform, look_at};
use num_traits::Float as NumFloat;
#[repr(C)]
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum CurveType {
Flat,
Cylinder,
Ribbon,
}
#[repr(C)]
#[derive(Debug, Clone, Copy)]
pub struct CurveCommon {
pub curve_type: CurveType,
pub cp_obj: [Point3f; 4],
pub width: [Float; 2],
pub n: [Normal3f; 2],
pub normal_angle: Float,
pub inv_sin_normal_angle: Float,
pub render_from_object: Transform,
pub object_from_render: Transform,
pub reverse_orientation: bool,
pub transform_swap_handedness: bool,
}
impl CurveCommon {
#[allow(clippy::too_many_arguments)]
pub fn new(
c: &[Point3f],
w0: Float,
w1: Float,
curve_type: CurveType,
norm: &[Normal3f],
render_from_object: Transform,
object_from_render: Transform,
reverse_orientation: bool,
) -> Self {
let transform_swap_handedness = render_from_object.swaps_handedness();
let width = [w0, w1];
assert_eq!(c.len(), 4);
let cp_obj: [Point3f; 4] = c[..4].try_into().unwrap();
let mut n: [Normal3f; 2] = gpu_array_from_fn(|_| Normal3f::default());
let mut normal_angle: Float = 0.;
let mut inv_sin_normal_angle: Float = 0.;
if norm.len() == 2 {
n[0] = norm[0].normalize();
n[1] = norm[1].normalize();
normal_angle = n[0].angle_between(n[1]);
inv_sin_normal_angle = 1. / normal_angle.sin();
}
Self {
curve_type,
cp_obj,
width,
n,
normal_angle,
inv_sin_normal_angle,
render_from_object,
object_from_render,
reverse_orientation,
transform_swap_handedness,
}
}
}
#[repr(C)]
#[derive(Debug, Clone, Copy)]
pub struct CurveShape {
pub common: CurveCommon,
pub u_min: Float,
pub u_max: Float,
}
#[repr(C)]
#[derive(Clone, Copy, Debug)]
struct IntersectionContext {
pub ray: Ray,
pub object_from_ray: Transform,
pub common: CurveCommon,
}
impl CurveShape {
pub fn new(common: CurveCommon, u_min: Float, u_max: Float) -> Self {
Self {
common,
u_min,
u_max,
}
}
fn intersect_ray(&self, r: &Ray, t_max: Float) -> Option<ShapeIntersection> {
let ray = self
.common
.object_from_render
.apply_to_ray(r, &mut Some(t_max));
let cp_obj = cubic_bezier_control_points(&self.common.cp_obj, self.u_min, self.u_max);
// Project curve control points to plane perpendicular to ray
let mut dx = ray.d.cross(cp_obj[3] - cp_obj[0]);
if dx.norm_squared() == 0. {
(dx, _) = ray.d.coordinate_system();
}
let ray_from_object = look_at(ray.o, ray.o + ray.d, dx).expect("Inversion error");
let cp: [Point3f; 4] = gpu_array_from_fn(|i| ray_from_object.apply_to_point(cp_obj[i]));
let max_width = lerp(self.u_min, self.common.width[0], self.common.width[1]).max(lerp(
self.u_max,
self.common.width[0],
self.common.width[1],
));
let curve_bounds = Bounds3f::from_points(cp[0], cp[1])
.union(Bounds3f::from_points(cp[2], cp[3]))
.expand(0.5 * max_width);
let ray_bounds =
Bounds3f::from_points(Point3f::zero(), Point3f::new(0., 0., ray.d.norm() * t_max));
if !ray_bounds.overlaps(&curve_bounds) {
return None;
}
let l0 = (0..2)
.map(|i| {
(cp[i].x() - 2.0 * cp[i + 1].x() + cp[i + 2].x())
.abs()
.max((cp[i].y() - 2.0 * cp[i + 1].y() + cp[i + 2].y()).abs())
.max((cp[i].z() - 2.0 * cp[i + 1].z() + cp[i + 2].z()).abs())
})
.fold(0.0, Float::max);
let max_depth = if l0 > 0. {
let eps = self.common.width[0].max(self.common.width[1]) * 0.05;
let r0: i32 = (1.41421356237 * 6. * l0 / (8. * eps)).log2() as i32 / 2;
clamp(r0, 0, 10)
} else {
0
};
let context = IntersectionContext {
ray,
object_from_ray: ray_from_object.inverse(),
common: self.common.clone(),
};
self.recursive_intersect(&context, t_max, &cp, self.u_min, self.u_max, max_depth)
}
fn recursive_intersect(
&self,
context: &IntersectionContext,
mut t_max: Float,
cp: &[Point3f],
u0: Float,
u1: Float,
depth: i32,
) -> Option<ShapeIntersection> {
if depth > 0 {
let cp_split = subdivide_cubic_bezier(cp);
let u = [u0, (u0 + u1) / 2., u1];
let mut best_hit: Option<ShapeIntersection> = None;
for seg in 0..2 {
let cps: &[Point3f] = &cp_split[3 * seg..3 * seg + 4];
let max_width = lerp(u[seg], self.common.width[0], self.common.width[1]).max(lerp(
u[seg + 1],
self.common.width[0],
self.common.width[1],
));
let curve_bounds = Bounds3f::from_points(cps[0], cps[1])
.union(Bounds3f::from_points(cps[2], cps[3]))
.expand(0.5 * max_width);
let ray_bounds = Bounds3f::from_points(
Point3f::zero(),
Point3f::new(0., 0., context.ray.d.norm() * t_max),
);
if !ray_bounds.overlaps(&curve_bounds) {
continue;
}
if let Some(hit) =
self.recursive_intersect(context, t_max, cps, u[seg], u[seg + 1], depth - 1)
{
best_hit = Some(hit);
t_max = best_hit.as_ref().unwrap().t_hit;
}
}
best_hit
} else {
self.intersect_segment(context, t_max, cp, u0, u1)
}
}
fn intersect_segment(
&self,
context: &IntersectionContext,
t_max: Float,
cp: &[Point3f],
u0: Float,
u1: Float,
) -> Option<ShapeIntersection> {
let edge1 = (cp[1].y() - cp[0].y()) * -cp[0].y() + cp[0].x() * (cp[0].x() - cp[1].x());
let edge2 = (cp[2].y() - cp[3].y()) * -cp[3].y() + cp[3].x() * (cp[3].x() - cp[2].x());
if edge1 <= 0.0 || edge2 <= 0.0 {
return None;
}
let segment_dir = Point2f::new(cp[3].x(), cp[3].y()) - Point2f::new(cp[0].x(), cp[0].y());
let denom = segment_dir.norm_squared();
if denom == 0. {
return None;
}
let w = Vector2f::new(cp[0].x(), cp[0].y()).dot(-segment_dir) / denom;
let u = clamp(lerp(w, u0, u1), u0, u1);
let ray_length = context.ray.d.norm();
let mut hit_width = lerp(u, self.common.width[0], self.common.width[1]);
let mut n_hit = Normal3f::zero();
if let CurveType::Ribbon = context.common.curve_type {
n_hit = if context.common.normal_angle == 0. {
context.common.n[0]
} else {
let sin0 =
((1. - u) * self.common.normal_angle).sin() * self.common.inv_sin_normal_angle;
let sin1 = (u * self.common.normal_angle).sin() * self.common.inv_sin_normal_angle;
sin0 * self.common.n[0] + sin1 * self.common.n[1]
};
hit_width = hit_width * n_hit.dot(context.ray.d.into()).abs() / ray_length;
}
let (pc, dpcdw) = evaluate_cubic_bezier(cp, clamp(w, 0., 1.));
let ray_length = context.ray.d.norm();
if !self.valid_hit(pc, hit_width, t_max, ray_length) {
return None;
}
// Hit is valid, obtain normals, and interaction differentials
self.intersection_result(context, pc, dpcdw, u, hit_width, n_hit, ray_length)
}
fn hit_ribbon(
&self,
u: Float,
hit_width: Float,
context: &IntersectionContext,
) -> (Float, Normal3f) {
let n_hit = if context.common.normal_angle == 0. {
context.common.n[0]
} else {
let sin0 =
((1. - u) * self.common.normal_angle).sin() * self.common.inv_sin_normal_angle;
let sin1 = (u * self.common.normal_angle).sin() * self.common.inv_sin_normal_angle;
sin0 * self.common.n[0] + sin1 * self.common.n[1]
};
let new_hit_width =
hit_width * n_hit.dot(context.ray.d.into()).abs() / context.ray.d.norm();
(new_hit_width, n_hit)
}
fn valid_hit(&self, pc: Point3f, hit_width: Float, t_max: Float, ray_length: Float) -> bool {
let pt_curve_dist_sq = square(pc.x()) + square(pc.y());
if pt_curve_dist_sq > square(hit_width * 0.5) || pc.z() < 0.0 || pc.z() > ray_length * t_max
{
return false;
}
true
}
#[allow(clippy::too_many_arguments)]
fn intersection_result(
&self,
context: &IntersectionContext,
pc: Point3f,
dpcdw: Vector3f,
u: Float,
hit_width: Float,
n_hit: Normal3f,
ray_length: Float,
) -> Option<ShapeIntersection> {
let t_hit = pc.z() / ray_length;
let pt_curve_dist = (square(pc.x()) + square(pc.y())).sqrt();
let edge_func = dpcdw.x() * -pc.y() + dpcdw.y() * pc.x();
let v = if edge_func > 0. {
0.5 + pt_curve_dist / hit_width
} else {
0.5 - pt_curve_dist / hit_width
};
let (_, dpdu) = evaluate_cubic_bezier(&self.common.cp_obj, u);
let dpdv = match context.common.curve_type {
CurveType::Ribbon => Vector3f::from(n_hit).cross(dpdu).normalize() * hit_width,
_ => {
let dpdu_plane = context.object_from_ray.apply_inverse_vector(dpdu);
let mut dpdv_plane =
Vector3f::new(-dpdu_plane.y(), dpdu_plane.x(), 0.).normalize() * hit_width;
if context.common.curve_type == CurveType::Cylinder {
let theta = lerp(v, -90., 90.);
let rot = Transform::rotate_around_axis(-theta, dpdu_plane);
dpdv_plane = rot.apply_to_vector(dpdv_plane);
}
context.object_from_ray.apply_to_vector(dpdv_plane)
}
};
let p_error = Vector3f::fill(hit_width);
let flip_normal = self.common.reverse_orientation ^ self.common.transform_swap_handedness;
let pi = Point3fi::new_with_error(context.ray.at(t_hit), p_error);
let intr = SurfaceInteraction::new(
pi,
Point2f::new(u, v),
-context.ray.d,
dpdu,
dpdv,
Normal3f::default(),
Normal3f::default(),
context.ray.time,
flip_normal,
);
Some(ShapeIntersection { intr, t_hit })
}
}
impl ShapeTrait for CurveShape {
fn bounds(&self) -> Bounds3f {
let cs_span = self.common.cp_obj;
let obj_bounds = bound_cubic_bezier(&cs_span, self.u_min, self.u_max);
let width0 = lerp(self.u_min, self.common.width[0], self.common.width[1]);
let width1 = lerp(self.u_max, self.common.width[0], self.common.width[1]);
let obj_bounds_expand = obj_bounds.expand(width0.max(width1) * 0.5);
self.common
.render_from_object
.apply_to_bounds(obj_bounds_expand)
}
fn normal_bounds(&self) -> DirectionCone {
DirectionCone::entire_sphere()
}
fn area(&self) -> Float {
let cp_obj = cubic_bezier_control_points(&self.common.cp_obj, self.u_min, self.u_max);
let width0 = lerp(self.u_min, self.common.width[0], self.common.width[1]);
let width1 = lerp(self.u_max, self.common.width[0], self.common.width[1]);
let avg_width = (width0 + width1) / 2.;
let mut approx_length = 0.;
for i in 0..3 {
approx_length += cp_obj[i].distance(cp_obj[i + 1]);
}
approx_length * avg_width
}
fn intersect_p(&self, ray: &Ray, t_max: Option<Float>) -> bool {
self.intersect_ray(ray, t_max.unwrap_or(Float::INFINITY))
.is_some()
}
fn intersect(&self, ray: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
self.intersect_ray(ray, t_max.unwrap_or(Float::INFINITY))
}
fn pdf(&self, _interaction: &Interaction) -> Float {
unimplemented!()
}
fn pdf_from_context(&self, _ctx: &ShapeSampleContext, _wi: Vector3f) -> Float {
unimplemented!()
}
fn sample(&self, _u: Point2f) -> Option<ShapeSample> {
unimplemented!()
}
fn sample_from_context(&self, _ctx: &ShapeSampleContext, _u: Point2f) -> Option<ShapeSample> {
unimplemented!()
}
}