kurbo/
offset.rs

Help
1// Copyright 2022 the Kurbo Authors
2// SPDX-License-Identifier: Apache-2.0 OR MIT
3
4//! Computation of offset curves of cubic Béziers, based on a curve fitting
5//! approach.
6//!
7//! See the [Parallel curves of cubic Béziers] blog post for a discussion of how
8//! this algorithm works and what kind of results can be expected. In general, it
9//! is expected to perform much better than most published algorithms. The number
10//! of curve segments needed to attain a given accuracy scales as O(n^6) with
11//! accuracy.
12//!
13//! In general, to compute the offset curve (also known as parallel curve) of
14//! a cubic Bézier segment, create a [`CubicOffset`] struct with the curve
15//! segment and offset, then use [`fit_to_bezpath`] or [`fit_to_bezpath_opt`]
16//! depending on how much time to spend optimizing the resulting path.
17//!
18//! [`fit_to_bezpath`]: crate::fit_to_bezpath
19//! [`fit_to_bezpath_opt`]: crate::fit_to_bezpath_opt
20//! [Parallel curves of cubic Béziers]: https://raphlinus.github.io/curves/2022/09/09/parallel-beziers.html
21use core::ops::Range;
22
23#[cfg(not(feature = "std"))]
24use crate::common::FloatFuncs;
25
26use crate::{
27    common::solve_itp, CubicBez, CurveFitSample, ParamCurve, ParamCurveDeriv, ParamCurveFit, Point,
28    QuadBez, Vec2,
29};
30
31/// The offset curve of a cubic Bézier.
32///
33/// This is a representation of the offset curve of a cubic Bézier segment, for
34/// purposes of curve fitting.
35///
36/// See the [module-level documentation] for a bit more discussion of the approach,
37/// and how this struct is to be used.
38///
39/// [module-level documentation]: crate::offset
40pub struct CubicOffset {
41    /// Source curve.
42    c: CubicBez,
43    /// Derivative of source curve.
44    q: QuadBez,
45    /// Offset.
46    d: f64,
47    // c0 + c1 t + c2 t^2 is the cross product of second and first
48    // derivatives of the underlying cubic, multiplied by offset (for
49    // computing cusp).
50    c0: f64,
51    c1: f64,
52    c2: f64,
53}
54
55impl CubicOffset {
56    /// Create a new curve from Bézier segment and offset.
57    ///
58    /// This method should only be used if the Bézier is smooth. Use
59    /// [`new_regularized`] instead to deal with a wider range of inputs.
60    ///
61    /// [`new_regularized`]: Self::new_regularized
62    pub fn new(c: CubicBez, d: f64) -> Self {
63        let q = c.deriv();
64        let d0 = q.p0.to_vec2();
65        let d1 = 2.0 * (q.p1 - q.p0);
66        let d2 = q.p0.to_vec2() - 2.0 * q.p1.to_vec2() + q.p2.to_vec2();
67        CubicOffset {
68            c,
69            q,
70            d,
71            c0: d * d1.cross(d0),
72            c1: d * 2.0 * d2.cross(d0),
73            c2: d * d2.cross(d1),
74        }
75    }
76
77    /// Create a new curve from Bézier segment and offset, with numerical robustness tweaks.
78    ///
79    /// The dimension represents a minimum feature size; the regularization is allowed to
80    /// perturb the curve by this amount in order to improve the robustness.
81    pub fn new_regularized(c: CubicBez, d: f64, dimension: f64) -> Self {
82        Self::new(c.regularize(dimension), d)
83    }
84
85    fn eval_offset(&self, t: f64) -> Vec2 {
86        let dp = self.q.eval(t).to_vec2();
87        let norm = Vec2::new(-dp.y, dp.x);
88        // TODO: deal with hypot = 0
89        norm * self.d / dp.hypot()
90    }
91
92    fn eval(&self, t: f64) -> Point {
93        // Point on source curve.
94        self.c.eval(t) + self.eval_offset(t)
95    }
96
97    /// Evaluate derivative of curve.
98    fn eval_deriv(&self, t: f64) -> Vec2 {
99        self.cusp_sign(t) * self.q.eval(t).to_vec2()
100    }
101
102    // Compute a function which has a zero-crossing at cusps, and is
103    // positive at low curvatures on the source curve.
104    fn cusp_sign(&self, t: f64) -> f64 {
105        let ds2 = self.q.eval(t).to_vec2().hypot2();
106        ((self.c2 * t + self.c1) * t + self.c0) / (ds2 * ds2.sqrt()) + 1.0
107    }
108}
109
110impl ParamCurveFit for CubicOffset {
111    fn sample_pt_tangent(&self, t: f64, sign: f64) -> CurveFitSample {
112        let p = self.eval(t);
113        const CUSP_EPS: f64 = 1e-8;
114        let mut cusp = self.cusp_sign(t);
115        if cusp.abs() < CUSP_EPS {
116            // This is a numerical derivative, which is probably good enough
117            // for all practical purposes, but an analytical derivative would
118            // be more elegant.
119            //
120            // Also, we're not dealing with second or higher order cusps.
121            cusp = sign * (self.cusp_sign(t + CUSP_EPS) - self.cusp_sign(t - CUSP_EPS));
122        }
123        let tangent = self.q.eval(t).to_vec2() * cusp.signum();
124        CurveFitSample { p, tangent }
125    }
126
127    fn sample_pt_deriv(&self, t: f64) -> (Point, Vec2) {
128        (self.eval(t), self.eval_deriv(t))
129    }
130
131    fn break_cusp(&self, range: Range<f64>) -> Option<f64> {
132        const CUSP_EPS: f64 = 1e-8;
133        // When an endpoint is on (or very near) a cusp, move just far enough
134        // away from the cusp that we're confident we have the right sign.
135        let break_cusp_help = |mut x, mut d| {
136            let mut cusp = self.cusp_sign(x);
137            while cusp.abs() < CUSP_EPS && d < 1.0 {
138                x += d;
139                let old_cusp = cusp;
140                cusp = self.cusp_sign(x);
141                if cusp.abs() > old_cusp.abs() {
142                    break;
143                }
144                d *= 2.0;
145            }
146            (x, cusp)
147        };
148        let (a, cusp0) = break_cusp_help(range.start, 1e-12);
149        let (b, cusp1) = break_cusp_help(range.end, -1e-12);
150        if a >= b || cusp0 * cusp1 >= 0.0 {
151            // Discussion point: maybe we should search for double cusps in the interior
152            // of the range.
153            return None;
154        }
155        let s = cusp1.signum();
156        let f = |t| s * self.cusp_sign(t);
157        let k1 = 0.2 / (b - a);
158        const ITP_EPS: f64 = 1e-12;
159        let x = solve_itp(f, a, b, ITP_EPS, 1, k1, s * cusp0, s * cusp1);
160        Some(x)
161    }
162}
163
164#[cfg(test)]
165mod tests {
166    use super::CubicOffset;
167    use crate::{fit_to_bezpath, fit_to_bezpath_opt, CubicBez, PathEl};
168
169    // This test tries combinations of parameters that have caused problems in the past.
170    #[test]
171    fn pathological_curves() {
172        let curve = CubicBez {
173            p0: (-1236.3746269978635, 152.17981429574826).into(),
174            p1: (-1175.18662093517, 108.04721798590596).into(),
175            p2: (-1152.142883879584, 105.76260301083356).into(),
176            p3: (-1151.842639804639, 105.73040758939104).into(),
177        };
178        let offset = 3603.7267536453924;
179        let accuracy = 0.1;
180        let offset_path = CubicOffset::new(curve, offset);
181        let path = fit_to_bezpath_opt(&offset_path, accuracy);
182        assert!(matches!(path.iter().next(), Some(PathEl::MoveTo(_))));
183        let path_opt = fit_to_bezpath(&offset_path, accuracy);
184        assert!(matches!(path_opt.iter().next(), Some(PathEl::MoveTo(_))));
185    }
186
187    /// Cubic offset that used to trigger infinite recursion.
188    #[test]
189    fn infinite_recursion() {
190        const DIM_TUNE: f64 = 0.25;
191        const TOLERANCE: f64 = 0.1;
192        let c = CubicBez::new(
193            (1096.2962962962963, 593.90243902439033),
194            (1043.6213991769548, 593.90243902439033),
195            (1030.4526748971193, 593.90243902439033),
196            (1056.7901234567901, 593.90243902439033),
197        );
198        let co = CubicOffset::new_regularized(c, -0.5, DIM_TUNE * TOLERANCE);
199        fit_to_bezpath(&co, TOLERANCE);
200    }
201
202    #[test]
203    fn test_cubic_offset_simple_line() {
204        let cubic = CubicBez::new((0., 0.), (10., 0.), (20., 0.), (30., 0.));
205        let offset = CubicOffset::new(cubic, 5.);
206        let _optimized = fit_to_bezpath(&offset, 1e-6);
207    }
208}
kurbo/offset.rs

kurbo/
offset.rs
