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}
209