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 pub fn new(c: CubicBez, d: f64) -> Self {
69 let q = c.deriv();
70 let d0 = q.p0.to_vec2();
71 let d1 = 2.0 * (q.p1 - q.p0);
72 let d2 = q.p0.to_vec2() - 2.0 * q.p1.to_vec2() + q.p2.to_vec2();
73 CubicOffset {
74 c,
75 q,
76 d,
77 c0: d * d1.cross(d0),
78 c1: d * 2.0 * d2.cross(d0),
79 c2: d * d2.cross(d1),
80 }
81 }
82
83 fn eval_offset(&self, t: f64) -> Vec2 {
84 let dp = self.q.eval(t).to_vec2();
85 let norm = Vec2::new(-dp.y, dp.x);
86 // TODO: deal with hypot = 0
87 norm * self.d / dp.hypot()
88 }
89
90 fn eval(&self, t: f64) -> Point {
91 // Point on source curve.
92 self.c.eval(t) + self.eval_offset(t)
93 }
94
95 /// Evaluate derivative of curve.
96 fn eval_deriv(&self, t: f64) -> Vec2 {
97 self.cusp_sign(t) * self.q.eval(t).to_vec2()
98 }
99
100 // Compute a function which has a zero-crossing at cusps, and is
101 // positive at low curvatures on the source curve.
102 fn cusp_sign(&self, t: f64) -> f64 {
103 let ds2 = self.q.eval(t).to_vec2().hypot2();
104 ((self.c2 * t + self.c1) * t + self.c0) / (ds2 * ds2.sqrt()) + 1.0
105 }
106}
107
108impl ParamCurveFit for CubicOffset {
109 fn sample_pt_tangent(&self, t: f64, sign: f64) -> CurveFitSample {
110 let p = self.eval(t);
111 const CUSP_EPS: f64 = 1e-8;
112 let mut cusp = self.cusp_sign(t);
113 if cusp.abs() < CUSP_EPS {
114 // This is a numerical derivative, which is probably good enough
115 // for all practical purposes, but an analytical derivative would
116 // be more elegant.
117 //
118 // Also, we're not dealing with second or higher order cusps.
119 cusp = sign * (self.cusp_sign(t + CUSP_EPS) - self.cusp_sign(t - CUSP_EPS));
120 }
121 let tangent = self.q.eval(t).to_vec2() * cusp.signum();
122 CurveFitSample { p, tangent }
123 }
124
125 fn sample_pt_deriv(&self, t: f64) -> (Point, Vec2) {
126 (self.eval(t), self.eval_deriv(t))
127 }
128
129 fn break_cusp(&self, range: Range<f64>) -> Option<f64> {
130 const CUSP_EPS: f64 = 1e-8;
131 // When an endpoint is on (or very near) a cusp, move just far enough
132 // away from the cusp that we're confident we have the right sign.
133 let break_cusp_help = |mut x, mut d| {
134 let mut cusp = self.cusp_sign(x);
135 while cusp.abs() < CUSP_EPS && d < 1.0 {
136 x += d;
137 let old_cusp = cusp;
138 cusp = self.cusp_sign(x);
139 if cusp.abs() > old_cusp.abs() {
140 break;
141 }
142 d *= 2.0;
143 }
144 (x, cusp)
145 };
146 let (a, cusp0) = break_cusp_help(range.start, 1e-12);
147 let (b, cusp1) = break_cusp_help(range.end, -1e-12);
148 if a >= b || cusp0 * cusp1 >= 0.0 {
149 // Discussion point: maybe we should search for double cusps in the interior
150 // of the range.
151 return None;
152 }
153 let s = cusp1.signum();
154 let f = |t| s * self.cusp_sign(t);
155 let k1 = 0.2 / (b - a);
156 const ITP_EPS: f64 = 1e-12;
157 let x = solve_itp(f, a, b, ITP_EPS, 1, k1, s * cusp0, s * cusp1);
158 Some(x)
159 }
160}
161