pbrt_r3/shapes/
curve.rs

1use crate::core::prelude::*;
2
3use std::sync::Arc;
4
5thread_local!(static TESTS: StatPercent = StatPercent::new("Intersections/Ray-curve intersection tests"));
6thread_local!(static REFINEMENT_LEVEL: StatIntDistribution = StatIntDistribution::new("Intersections/Curve refinement level"));
7thread_local!(static N_CURVES: StatCounter = StatCounter::new("Scene/Curves"));
8thread_local!(static N_SPLIT_CURVES: StatCounter = StatCounter::new("Scene/Split curves"));
9
10// CurveType Declarations
11#[derive(Debug, Clone, Copy, PartialEq)]
12pub enum CurveType {
13    Flat,
14    Cylinder,
15    Ribbon,
16}
17
18#[derive(Debug, Clone, Copy)]
19pub struct CurveNormal {
20    pub n: [Normal3f; 2],
21    pub normal_angle: Float,
22    pub inv_sin_normal_angle: Float,
23}
24
25// CurveCommon Declarations
26#[derive(Debug, Clone)]
27pub struct CurveCommon {
28    pub base: BaseShape,
29    pub t: CurveType,
30    pub cp_obj: [Point3f; 4],
31    pub width: [Float; 2],
32    pub normals: Option<CurveNormal>,
33}
34
35// Curve Declarations
36#[derive(Debug, Clone)]
37pub struct Curve {
38    common: Arc<CurveCommon>,
39    u_min: Float,
40    u_max: Float,
41}
42
43// Curve Utility Functions
44#[inline]
45fn lerp3(t: Float, p0: &Point3f, p1: &Point3f) -> Point3f {
46    return Point3f::new(
47        lerp(t, p0.x, p1.x),
48        lerp(t, p0.y, p1.y),
49        lerp(t, p0.z, p1.z),
50    );
51}
52
53#[inline]
54fn blossom_bezier(cp: &[Point3f], u0: Float, u1: Float, u2: Float) -> Point3f {
55    let a = [
56        lerp3(u0, &cp[0], &cp[1]),
57        lerp3(u0, &cp[1], &cp[2]),
58        lerp3(u0, &cp[2], &cp[3]),
59    ];
60    let b = [lerp3(u1, &a[0], &a[1]), lerp3(u1, &a[1], &a[2])];
61    return lerp3(u2, &b[0], &b[1]);
62}
63
64#[inline]
65fn subdivide_bezier(cp: &[Point3f]) -> [Point3f; 7] {
66    let cp_split = [
67        cp[0],
68        (cp[0] + cp[1]) * 0.5,
69        (cp[0] + 2.0 * cp[1] + cp[2]) * 0.25,
70        (cp[0] + 3.0 * cp[1] + 3.0 * cp[2] + cp[3]) * 0.125,
71        (cp[1] + 2.0 * cp[2] + cp[3]) * 0.25,
72        (cp[2] + cp[3]) * 0.5,
73        cp[3],
74    ];
75    return cp_split;
76}
77
78#[inline]
79fn eval_bezier(cp: &[Point3f], u: Float) -> (Point3f, Vector3f) {
80    let cp1 = [
81        lerp3(u, &cp[0], &cp[1]),
82        lerp3(u, &cp[1], &cp[2]),
83        lerp3(u, &cp[2], &cp[3]),
84    ];
85    let cp2 = [lerp3(u, &cp1[0], &cp1[1]), lerp3(u, &cp1[1], &cp1[2])];
86    let cp3 = lerp3(u, &cp2[0], &cp2[1]);
87    let delta = cp2[1] - cp2[0];
88    let deriv = if delta.length_squared() > 0.0 {
89        3.0 * delta
90    } else {
91        // For a cubic Bezier, if the first three control points (say) are
92        // coincident, then the derivative of the curve is legitimately (0,0,0)
93        // at u=0.  This is problematic for us, though, since we'd like to be
94        // able to compute a surface normal there.  In that case, just punt and
95        // take the difference between the first and last control points, which
96        // ain't great, but will hopefully do.
97        cp[3] - cp[0]
98    };
99    return (cp3, deriv);
100}
101
102#[inline]
103fn transform_lookat(e: &Point3f, l: &Point3f, u: &Vector3f) -> Transform {
104    return Transform::look_at(e.x, e.y, e.z, l.x, l.y, l.z, u.x, u.y, u.z);
105}
106
107impl CurveCommon {
108    pub fn new(
109        o2w: &Transform,
110        w2o: &Transform,
111        reverse_orientation: bool,
112        c: &[Point3f],
113        w0: Float,
114        w1: Float,
115        t: CurveType,
116        norm: &Option<Vec<Normal3f>>,
117    ) -> Self {
118        let cp_obj: Vec<Point3f> = c.to_vec();
119        let cp_obj: [Point3f; 4] = cp_obj.try_into().unwrap();
120        let width = [w0, w1];
121        let normals: Option<CurveNormal> = match norm {
122            Some(n) => {
123                let n0 = n[0].normalize();
124                let n1 = n[1].normalize();
125                let normal_angle = Float::acos(Float::clamp(Normal3f::dot(&n0, &n1), 0.0, 1.0));
126                let inv_sin_normal_angle = 1.0 / Float::sin(normal_angle);
127                Some(CurveNormal {
128                    n: [n0, n1],
129                    normal_angle,
130                    inv_sin_normal_angle,
131                })
132            }
133            None => None,
134        };
135        N_CURVES.with(|n| n.inc());
136        CurveCommon {
137            base: BaseShape::new(o2w, w2o, reverse_orientation),
138            t,
139            cp_obj,
140            width,
141            normals,
142        }
143    }
144}
145
146fn get_bounds(cp: &[Point3f]) -> Bounds3f {
147    let min = cp.iter().fold(cp[0], |a, b| -> Vector3f {
148        return Vector3f::new(
149            Float::min(a[0], b[0]),
150            Float::min(a[1], b[1]),
151            Float::min(a[2], b[2]),
152        );
153    });
154    let max = cp.iter().fold(cp[0], |a, b| -> Vector3f {
155        return Vector3f::new(
156            Float::max(a[0], b[0]),
157            Float::max(a[1], b[1]),
158            Float::max(a[2], b[2]),
159        );
160    });
161    return Bounds3f::from(((min[0], min[1], min[2]), (max[0], max[1], max[2])));
162}
163
164fn recursive_intersect(
165    common: &CurveCommon,
166    is_shadow: bool,
167    ray: &Ray,
168    cp: &[Point3f],
169    ray_to_object: &Transform,
170    u0: Float,
171    u1: Float,
172    t_max: Float,
173    depth: i32,
174) -> Option<(Float, Option<SurfaceInteraction>)> {
175    let mut t_max = t_max;
176    let width = &common.width;
177    let t = common.t;
178    let normals = &common.normals;
179    let ray_length = ray.d.length();
180
181    if depth > 0 {
182        // Split curve segment into sub-segments and test for intersection
183        let cp_split = subdivide_bezier(cp);
184        // For each of the two segments, see if the ray's bounding box
185        // overlaps the segment before recursively checking for
186        // intersection with it.
187        let u = [u0, (u0 + u1) / 2.0, u1];
188        let mut t_hit = None;
189        // Pointer to the 4 control poitns for the current segment.
190        for seg in 0..2 {
191            let z_max = ray_length * t_max;
192            //0:0,1,2,3
193            //1:3,4,5,6
194            let cps = &cp_split[(3 * seg)..];
195            let max_width = Float::max(
196                lerp(u[seg], width[0], width[1]),
197                lerp(u[seg + 1], width[0], width[1]),
198            );
199            let max_radius = 0.5 * max_width;
200
201            // As above, check y first, since it most commonly lets us exit
202            // out early.
203            let max_y = Float::max(
204                Float::max(cps[0].y, cps[1].y),
205                Float::max(cps[2].y, cps[3].y),
206            );
207            let min_y = Float::min(
208                Float::min(cps[0].y, cps[1].y),
209                Float::min(cps[2].y, cps[3].y),
210            );
211            if ((max_y + max_radius) < 0.0) || ((min_y - max_radius) > 0.0) {
212                continue;
213            }
214
215            let max_x = Float::max(
216                Float::max(cps[0].x, cps[1].x),
217                Float::max(cps[2].x, cps[3].x),
218            );
219            let min_x = Float::min(
220                Float::min(cps[0].x, cps[1].x),
221                Float::min(cps[2].x, cps[3].x),
222            );
223            if ((max_x + max_radius) < 0.0) || ((min_x - max_radius) > 0.0) {
224                continue;
225            }
226
227            let max_z = Float::max(
228                Float::max(cps[0].z, cps[1].z),
229                Float::max(cps[2].z, cps[3].z),
230            );
231            let min_z = Float::min(
232                Float::min(cps[0].z, cps[1].z),
233                Float::min(cps[2].z, cps[3].z),
234            );
235            if ((max_z + max_radius) < 0.0) || ((min_z - max_radius) > z_max) {
236                continue;
237            }
238
239            {
240                let t_t_hit = recursive_intersect(
241                    common,
242                    is_shadow,
243                    ray,
244                    cps,
245                    ray_to_object,
246                    u[seg],
247                    u[seg + 1],
248                    t_max,
249                    depth - 1,
250                );
251                if let Some(t_t_hit) = t_t_hit {
252                    // If we found an intersection and this is a shadow ray,
253                    // we can exit out immediately.
254                    if is_shadow {
255                        return Some(t_t_hit);
256                    }
257                    //if t_t_hit.0 < t_max {
258                    t_max = Float::min(t_max, t_t_hit.0);
259                    t_hit = Some(t_t_hit);
260                    //}
261                }
262            }
263        }
264        return t_hit;
265    } else {
266        // Intersect ray with curve segment
267
268        // Test ray against segment endpoint boundaries
269
270        // Test sample point against tangent perpendicular at curve start
271        let edge = (cp[1].y - cp[0].y) * -cp[0].y + cp[0].x * (cp[0].x - cp[1].x);
272        if edge < 0.0 {
273            return None;
274        }
275
276        // Test sample point against tangent perpendicular at curve end
277        let edge = (cp[2].y - cp[3].y) * -cp[3].y + cp[3].x * (cp[3].x - cp[2].x);
278        if edge < 0.0 {
279            return None;
280        }
281
282        // Compute line $w$ that gives minimum distance to sample point
283        let segment_direction = Point2f::new(cp[3][0], cp[3][1]) - Point2f::new(cp[0][0], cp[0][1]);
284        let denom = segment_direction.length_squared();
285        if denom == 0.0 {
286            return None;
287        }
288        let w = Vector2f::dot(&-Vector2f::new(cp[0][0], cp[0][1]), &segment_direction) / denom;
289        // Compute $u$ coordinate of curve intersection point and _hitWidth_
290        let u = Float::clamp(lerp(w, u0, u1), u0, u1);
291
292        let mut hit_width = lerp(u, width[0], width[1]);
293        let mut n_hit = Normal3f::new(1.0, 0.0, 0.0);
294        if t == CurveType::Ribbon {
295            if let Some(n) = normals {
296                // Scale _hitWidth_ based on ribbon orientation
297                let sin0 = Float::sin((1.0 - u) * n.normal_angle) * n.inv_sin_normal_angle;
298                let sin1 = Float::sin(u * n.normal_angle) * n.inv_sin_normal_angle;
299                n_hit = (sin0 * n.n[0] + sin1 * n.n[1]).normalize();
300                hit_width *= Vector3f::abs_dot(&n_hit, &ray.d) / ray_length;
301            }
302        }
303
304        // Test intersection point against curve width
305        let z_max = ray_length * t_max;
306        let (pc, dpcdw) = eval_bezier(cp, Float::clamp(w, 0.0, 1.0));
307        let pt_curve_dist2 = pc.x * pc.x + pc.y * pc.y;
308        if pt_curve_dist2 > (hit_width * hit_width * 0.25) {
309            return None;
310        }
311        if pc.z < 0.0 || pc.z > z_max {
312            return None;
313        }
314
315        // Compute $v$ coordinate of curve intersection point
316        let pt_curve_dist = Float::sqrt(pt_curve_dist2);
317        let edge_func = dpcdw.x * -pc.y + pc.x * dpcdw.y;
318        let v = if edge_func > 0.0 {
319            0.5 + pt_curve_dist / hit_width
320        } else {
321            0.5 - pt_curve_dist / hit_width
322        };
323
324        // Compute hit _t_ and partial derivatives for curve intersection
325
326        // FIXME: this tHit isn't quite right for ribbons...
327        let t_hit = pc.z / ray_length;
328        //println!("{}", t_hit);
329
330        TESTS.with(|stat| stat.add_num(1));
331
332        if !is_shadow {
333            // Compute error bounds for curve intersection
334            let p_error_coef = 2.0;
335            let p_error = Vector3f::new(
336                p_error_coef * hit_width,
337                p_error_coef * hit_width,
338                p_error_coef * hit_width,
339            );
340            //let p_error = ray_to_object.transform_vector(&p_error);//object space;
341
342            // Compute $\dpdu$ and $\dpdv$ for curve intersection
343            let (_, dpdu) = eval_bezier(&common.cp_obj, u); //object space
344
345            let dpdv; // = Vector3f::zero();
346            if t == CurveType::Ribbon {
347                dpdv = Vector3f::normalize(&Vector3f::cross(&n_hit, &dpdu)) * hit_width;
348            } else {
349                let object_to_ray = ray_to_object.inverse();
350                // Compute curve $\dpdv$ for flat and cylinder curves
351                let dpdu_plane = object_to_ray.transform_vector(&dpdu).normalize(); //ray space
352                let p_plane = Vector3f::new(0.0, 0.0, 1.0);
353                let mut dpdv_plane = Vector3f::cross(&p_plane, &dpdu_plane);
354                //let mut dpdv_plane =
355                //    Vector3f::normalize(&Vector3f::new(-dpdu_plane.y, dpdu_plane.x, 0.0))
356                //        * hit_width;
357                if t == CurveType::Cylinder {
358                    // Rotate _dpdvPlane_ to give cylindrical appearance
359                    let theta = lerp(v, -90.0, 90.0);
360                    let rot = Transform::rotate(-theta, dpdu_plane.x, dpdu_plane.y, dpdu_plane.z);
361                    dpdv_plane = rot.transform_vector(&dpdv_plane);
362                }
363                dpdv = ray_to_object.transform_vector(&dpdv_plane).normalize(); //object space
364            }
365            let p = ray.o + t_hit * ray.d;
366            let uv = Point2f::new(u, v); //u, v
367            let wo = -ray.d;
368            let mut n = common.base.calc_normal(&dpdu, &dpdv);
369            //if t == CurveType::Cylinder {
370            //    n = (p - p_center).normalize();
371            //} else {
372            if Vector3f::dot(&ray.d, &n) > 0.0 {
373                n *= -1.0;
374            }
375            //}
376            let si = SurfaceInteraction::new(
377                &p,
378                &p_error,
379                &uv,
380                &wo,
381                &n,
382                &dpdu,
383                &dpdv,
384                &Vector3f::zero(),
385                &Vector3f::zero(),
386                ray.time,
387                0,
388            );
389            let si = common
390                .base
391                .object_to_world
392                .transform_surface_interaction(&si);
393            return Some((t_hit, Some(si)));
394        } else {
395            return Some((t_hit, None));
396        }
397    }
398}
399
400fn log2(x: Float) -> i32 {
401    if x < 1.0 {
402        return 0;
403    } else {
404        return ((Float::log2(x) + 0.5) * 10.0) as i32;
405    }
406}
407
408impl Curve {
409    pub fn new(common: &Arc<CurveCommon>, u_min: Float, u_max: Float) -> Self {
410        Curve {
411            common: common.clone(),
412            u_min,
413            u_max,
414        }
415    }
416
417    fn intersect_common(
418        &self,
419        r: &Ray,
420        is_shadow: bool,
421    ) -> Option<(Float, Option<SurfaceInteraction>)> {
422        TESTS.with(|stat| stat.add_denom(1));
423
424        let u_min = self.u_min;
425        let u_max = self.u_max;
426        let common = self.common.as_ref();
427        // Transform _Ray_ to object space
428        let (ray, _o_err, _d_err) = common.base.world_to_object.transform_ray(r);
429        //let ray_length = ray.d.length();
430        let t_max = ray.t_max.get();
431        //println!("rl:{}, tmax:{}", ray_length, t_max);
432
433        // Compute object-space control points for curve segment, _cpObj_
434        let cp_obj = [
435            blossom_bezier(&common.cp_obj, u_min, u_min, u_min),
436            blossom_bezier(&common.cp_obj, u_min, u_min, u_max),
437            blossom_bezier(&common.cp_obj, u_min, u_max, u_max),
438            blossom_bezier(&common.cp_obj, u_max, u_max, u_max),
439        ];
440
441        // Project curve control points to plane perpendicular to ray
442
443        // Be careful to set the "up" direction passed to LookAt() to equal the
444        // vector from the first to the last control points.  In turn, this
445        // helps orient the curve to be roughly parallel to the x axis in the
446        // ray coordinate system.
447        //
448        // In turn (especially for curves that are approaching stright lines),
449        // we get curve bounds with minimal extent in y, which in turn lets us
450        // early out more quickly in recursiveIntersect().
451        let delta = cp_obj[3] - cp_obj[0];
452        let mut dx = Vector3f::cross(&ray.d, &delta);
453        if dx.length_squared() == 0.0 {
454            // If the ray and the vector between the first and last control
455            // points are parallel, dx will be zero.  Generate an arbitrary xy
456            // orientation for the ray coordinate system so that intersection
457            // tests can proceeed in this unusual case.
458            let (_dx, _) = coordinate_system(&ray.d);
459            dx = _dx;
460        }
461        let object_to_ray = transform_lookat(&ray.o, &(ray.o + ray.d), &dx);
462        let cp = [
463            object_to_ray.transform_point(&cp_obj[0]),
464            object_to_ray.transform_point(&cp_obj[1]),
465            object_to_ray.transform_point(&cp_obj[2]),
466            object_to_ray.transform_point(&cp_obj[3]),
467        ];
468
469        // Before going any further, see if the ray's bounding box intersects
470        // the curve's bounding box. We start with the y dimension, since the y
471        // extent is generally the smallest (and is often tiny) due to our
472        // careful orientation of the ray coordinate ysstem above.
473        let max_width = Float::max(
474            lerp(u_min, common.width[0], common.width[1]),
475            lerp(u_max, common.width[0], common.width[1]),
476        );
477        let max_radius = 0.5 * max_width;
478
479        // Check for non-overlap in y.
480        {
481            let max_y = Float::max(Float::max(cp[0].y, cp[1].y), Float::max(cp[2].y, cp[3].y));
482            let min_y = Float::min(Float::min(cp[0].y, cp[1].y), Float::min(cp[2].y, cp[3].y));
483            if ((max_y + max_radius) < 0.0) || ((min_y - max_radius) > 0.0) {
484                return None;
485            }
486        }
487
488        // Check for non-overlap in x.
489        {
490            let max_x = Float::max(Float::max(cp[0].x, cp[1].x), Float::max(cp[2].x, cp[3].x));
491            let min_x = Float::min(Float::min(cp[0].x, cp[1].x), Float::min(cp[2].x, cp[3].x));
492            if ((max_x + max_radius) < 0.0) || ((min_x - max_radius) > 0.0) {
493                return None;
494            }
495        }
496
497        // Check for non-overlap in z.
498        {
499            let ray_length = ray.d.length();
500            let z_max = ray_length * t_max;
501
502            let max_z = Float::max(Float::max(cp[0].z, cp[1].z), Float::max(cp[2].z, cp[3].z));
503            let min_z = Float::min(Float::min(cp[0].z, cp[1].z), Float::min(cp[2].z, cp[3].z));
504            if ((max_z + max_radius) < 0.0) || ((min_z - max_radius) > z_max) {
505                return None;
506            }
507        }
508
509        // Compute refinement depth for curve, _maxDepth_
510        let mut l0 = 0.0;
511        for i in 0..2 {
512            l0 = Float::max(
513                l0,
514                Float::max(
515                    Float::max(
516                        Float::abs(cp[i].x - 2.0 * cp[i + 1].x + cp[i + 2].x),
517                        Float::abs(cp[i].y - 2.0 * cp[i + 1].y + cp[i + 2].y),
518                    ),
519                    Float::abs(cp[i].z - 2.0 * cp[i + 1].z + cp[i + 2].z),
520                ),
521            );
522        }
523        let eps = Float::max(common.width[0], common.width[1]) * 0.05; // width / 20
524
525        // Compute log base 4 by dividing log2 in half.
526        //let r0 = log2(1.41421356237 * 6.0 * l0 / (8.0 * eps)) / 2;
527        let r0 = log2(SQRT_2 * 6.0 * l0 / (8.0 * eps)) / 2;
528        //let r0 = log2(1.41421356237 * 6.0 * l0 / (8.0 * eps) * 16.0) / 2;//x16
529
530        let max_depth = i32::clamp(r0, 0, 10); //
531
532        REFINEMENT_LEVEL.with(|stat| stat.add(max_depth as u64));
533
534        return recursive_intersect(
535            &common,
536            is_shadow,
537            &ray,
538            &cp,
539            &Transform::inverse(&object_to_ray),
540            u_min,
541            u_max,
542            t_max,
543            max_depth,
544        );
545    }
546}
547
548impl Shape for Curve {
549    fn object_bound(&self) -> Bounds3f {
550        let u_min = self.u_min;
551        let u_max = self.u_max;
552        let common = self.common.as_ref();
553        let cp_obj = [
554            blossom_bezier(&common.cp_obj, u_min, u_min, u_min),
555            blossom_bezier(&common.cp_obj, u_min, u_min, u_max),
556            blossom_bezier(&common.cp_obj, u_min, u_max, u_max),
557            blossom_bezier(&common.cp_obj, u_max, u_max, u_max),
558        ];
559        let b = get_bounds(&cp_obj);
560        let width = [
561            lerp(u_min, common.width[0], common.width[1]),
562            lerp(u_max, common.width[0], common.width[1]),
563        ];
564        let radius = Float::max(width[0], width[1]) * 0.5;
565        let b = b.expand(radius);
566        return b;
567    }
568
569    fn world_bound(&self) -> Bounds3f {
570        let common = self.common.as_ref();
571        let b = self.object_bound();
572        let b = common.base.object_to_world.transform_bounds(&b);
573        return b;
574    }
575
576    fn intersect(&self, r: &Ray) -> Option<(Float, SurfaceInteraction)> {
577        if let Some((t_hit, isect)) = self.intersect_common(r, false) {
578            return Some((t_hit, isect.unwrap()));
579        } else {
580            return None;
581        }
582    }
583
584    fn intersect_p(&self, r: &Ray) -> bool {
585        return self.intersect_common(r, true).is_some();
586    }
587
588    fn area(&self) -> Float {
589        // Compute object-space control points for curve segment, _cpObj_
590        let u_min = self.u_min;
591        let u_max = self.u_max;
592        let common = self.common.as_ref();
593        let cp_obj = [
594            blossom_bezier(&common.cp_obj, u_min, u_min, u_min),
595            blossom_bezier(&common.cp_obj, u_min, u_min, u_max),
596            blossom_bezier(&common.cp_obj, u_min, u_max, u_max),
597            blossom_bezier(&common.cp_obj, u_max, u_max, u_max),
598        ];
599        let width0 = lerp(u_min, common.width[0], common.width[1]);
600        let width1 = lerp(u_max, common.width[0], common.width[1]);
601        let avg_width = (width0 + width1) * 0.5;
602        let mut approx_length = 0.0;
603        for i in 0..3 {
604            approx_length += Vector3f::distance(&cp_obj[i], &cp_obj[i + 1]);
605        }
606        return approx_length * avg_width;
607    }
608
609    fn sample(&self, _u: &Point2f) -> Option<(Interaction, Float)> {
610        //LOG(FATAL) << "Curve::Sample not implemented.";
611        return None;
612    }
613}
614
615fn create_curve(
616    o2w: &Transform,
617    w2o: &Transform,
618    reverse_orientation: bool,
619    c: &[Point3f],
620    w0: Float,
621    w1: Float,
622    t: CurveType,
623    norm: &Option<Vec<Normal3f>>,
624    split_depth: i32,
625) -> Result<Vec<Arc<dyn Shape>>, PbrtError> {
626    let mut segments: Vec<Arc<dyn Shape>> = Vec::new();
627    let common = Arc::new(CurveCommon::new(
628        o2w,
629        w2o,
630        reverse_orientation,
631        c,
632        w0,
633        w1,
634        t,
635        norm,
636    ));
637    let n_segments = (1 << split_depth) as usize;
638    segments.reserve(n_segments);
639    for i in 0..n_segments {
640        let u_min = (i as Float) / (n_segments as Float);
641        let u_max = ((i + 1) as Float) / (n_segments as Float);
642        let curve = Arc::new(Curve::new(&common, u_min, u_max));
643        segments.push(curve);
644    }
645    N_SPLIT_CURVES.with(|c| c.add(n_segments as u64));
646    return Ok(segments);
647}
648
649fn get_curve_type(s: &str) -> Result<CurveType, PbrtError> {
650    return match s {
651        "flat" => Ok(CurveType::Flat),
652        "ribbon" => Ok(CurveType::Ribbon),
653        "cylinder" => Ok(CurveType::Cylinder),
654        _ => {
655            let msg = format!("Unknown curve type \"{}\".  Using \"cylinder\".", s);
656            return Err(PbrtError::error(&msg));
657        }
658    };
659}
660
661fn convert_to_point(f: &[Float]) -> Vec<Point3f> {
662    let l = f.len() / 3;
663    let mut v = Vec::new();
664    for i in 0..l {
665        v.push(Vector3f::new(f[3 * i + 0], f[3 * i + 1], f[3 * i + 2]));
666    }
667    return v;
668}
669
670pub fn create_curve_shape(
671    o2w: &Transform,
672    w2o: &Transform,
673    reverse_orientation: bool,
674    params: &ParamSet,
675) -> Result<Vec<Arc<dyn Shape>>, PbrtError> {
676    let width = params.find_one_float("width", 1.0);
677    let width0 = params.find_one_float("width0", width);
678    let width1 = params.find_one_float("width1", width);
679
680    let degree = params.find_one_int("degree", 3);
681    if degree != 2 && degree != 3 {
682        let msg = format!(
683            "Invalid degree {}: only degree 2 and 3 curves are supported.",
684            degree
685        );
686        return Err(PbrtError::error(&msg));
687    }
688    let degree = degree as usize;
689    let basis = params.find_one_string("basis", "bezier");
690    if basis != "bezier" && basis != "bspline" {
691        let msg = format!(
692            "Invalid basis \"{}\": only \"bezier\" and \"bspline\" are supported.",
693            degree
694        );
695        return Err(PbrtError::error(&msg));
696    }
697    let cp = params.get_points_ref("P").ok_or("\"P\"")?;
698    let cp = convert_to_point(&cp);
699    let ncp = cp.len();
700    let n_segments: usize = match &basis as &str {
701        "bezier" => {
702            // After the first segment, which uses degree+1 control points,
703            // subsequent segments reuse the last control point of the previous
704            // one and then use degree more control points.
705            if ((ncp - 1 - degree) % degree) != 0 {
706                let msg = format!(
707                    "Invalid number of control points {}: for the degree {} Bezier basis {} + n * {} are required, for n >= 0.",
708                    ncp, degree, degree +1, degree
709                );
710                return Err(PbrtError::error(&msg));
711            }
712            (ncp - 1) / degree
713        }
714        "bspline" => ncp - degree,
715        _ => {
716            panic!();
717        }
718    };
719
720    let curve_type = params.find_one_string("type", "flat");
721    let curve_type = get_curve_type(&curve_type)?;
722
723    let n = params.get_points_ref("N");
724    /*
725    let mut n = if let Some(nn) = n.as_ref() {
726        Some(convert_to_point(nn))
727    } else {
728        None
729    };
730    */
731    let mut n = n.as_ref().map(|nn| convert_to_point(nn));
732    if let Some(nn) = n.as_ref() {
733        let nnorm = nn.len();
734        if curve_type != CurveType::Ribbon {
735            //Warning("Curve normals are only used with \"ribbon\" type curves.");
736            n = None;
737        } else if nnorm != (n_segments + 1) {
738            let msg = format!("Invalid number of normals {}: must provide {} normals for ribbon curves with {} segments.", nnorm, n_segments + 1, n_segments);
739            return Err(PbrtError::error(&msg));
740        }
741    } else {
742        if curve_type == CurveType::Ribbon {
743            return Err(PbrtError::error(
744                "Must provide normals \"N\" at curve endpoints with ribbon curves.",
745            ));
746        }
747    }
748    let sd = params.find_one_float("splitdepth", 3.0) as i32;
749    let sd = params.find_one_int("splitdepth", sd);
750
751    let mut curves: Vec<Arc<dyn Shape>> = Vec::new();
752    for seg in 0..n_segments {
753        // First, compute the cubic Bezier control points for the current
754        // segment and store them in segCpBezier. (It is admittedly
755        // wasteful storage-wise to turn b-splines into Bezier segments and
756        // wasteful computationally to turn quadratic curves into cubics,
757        // but yolo.)
758        let cp_base = &cp[(seg * degree)..];
759        let mut seg_cp_bezier = vec![
760            Point3f::zero(),
761            Point3f::zero(),
762            Point3f::zero(),
763            Point3f::zero(),
764        ];
765        if basis == "bezier" {
766            if degree == 2 {
767                // Elevate to degree 3.
768                seg_cp_bezier[0] = cp_base[0];
769                seg_cp_bezier[1] = lerp3(2.0 / 3.0, &cp_base[0], &cp_base[1]);
770                seg_cp_bezier[2] = lerp3(1.0 / 3.0, &cp_base[1], &cp_base[2]);
771                seg_cp_bezier[3] = cp_base[2];
772            } else {
773                //for i in 0..4 {
774                //    seg_cp_bezier[i] = cp_base[i];
775                //}
776                seg_cp_bezier.copy_from_slice(&cp_base[..4]);
777            }
778        } else {
779            // Uniform b-spline.
780            if degree == 2 {
781                // First compute equivalent Bezier control points via some
782                // blossiming.  We have three control points and a uniform
783                // knot vector; we'll label the points p01, p12, and p23.
784                // We want the Bezier control points of the equivalent
785                // curve, which are p11, p12, and p22.
786                let p01 = cp_base[0];
787                let p12 = cp_base[1];
788                let p23 = cp_base[2];
789
790                // We already have p12.
791                let p11 = lerp3(0.5, &p01, &p12);
792                let p22 = lerp3(0.5, &p12, &p23);
793
794                // Now elevate to degree 3.
795                seg_cp_bezier[0] = p11;
796                seg_cp_bezier[1] = lerp3(2.0 / 3.0, &p11, &p12);
797                seg_cp_bezier[2] = lerp3(1.0 / 3.0, &p12, &p22);
798                seg_cp_bezier[3] = p22;
799            } else {
800                // Otherwise we will blossom from p012, p123, p234, and p345
801                // to the Bezier control points p222, p223, p233, and p333.
802                // https://people.eecs.berkeley.edu/~sequin/CS284/IMGS/cubicbsplinepoints.gif
803                let p012 = cp_base[0];
804                let p123 = cp_base[1];
805                let p234 = cp_base[2];
806                let p345 = cp_base[3];
807
808                let p122 = lerp3(2.0 / 3.0, &p012, &p123);
809                let p223 = lerp3(1.0 / 3.0, &p123, &p234);
810                let p233 = lerp3(2.0 / 3.0, &p123, &p234);
811                let p334 = lerp3(1.0 / 3.0, &p234, &p345);
812
813                let p222 = lerp3(0.5, &p122, &p223);
814                let p333 = lerp3(0.5, &p233, &p334);
815
816                seg_cp_bezier[0] = p222;
817                seg_cp_bezier[1] = p223;
818                seg_cp_bezier[2] = p233;
819                seg_cp_bezier[3] = p333;
820            }
821        }
822        let w0 = lerp((seg as Float) / (n_segments as Float), width0, width1);
823        let w1 = lerp(((seg + 1) as Float) / (n_segments as Float), width0, width1);
824        let mut c = create_curve(
825            o2w,
826            w2o,
827            reverse_orientation,
828            &seg_cp_bezier,
829            w0,
830            w1,
831            curve_type,
832            &n,
833            sd,
834        )?;
835        curves.append(&mut c);
836    }
837
838    return Ok(curves);
839}
pbrt_r3/shapes/curve.rs

pbrt_r3/shapes/
curve.rs
