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