8000 go-stuff/curve/fit.go at main · dominikh/go-stuff · GitHub
[go: up one dir, main page]
Skip to content

fit.go

Latest commit

 

History

History
821 lines (784 loc) · 24.3 KB

fit.go

File metadata and controls

821 lines (784 loc) · 24.3 KB
Raw
// SPDX-FileCopyrightText: 2018 Raph Levien
// SPDX-FileCopyrightText: 2024 Dominik Honnef and contributors
//
// SPDX-License-Identifier: MIT
// SPDX-FileAttributionText: https://github.com/linebender/kurbo
package curve
import (
"fmt"
"iter"
"math"
"honnef.co/go/stuff/container/maybe"
"honnef.co/go/stuff/math/polyroot"
)
// As described in [Simplifying Bézier paths], strictly optimizing for
// Fréchet distance can create bumps. The problem is curves with long
// control arms (distance from the control point to the corresponding
// endpoint). We mitigate that by applying a penalty as a multiplier to
// the measured error (approximate Fréchet distance). This is ReLU-like,
// with a value of 1.0 below the elbow, and a given slope above it. The
// values here have been determined empirically to give good results.
//
// [Simplifying Bézier paths]: https://raphlinus.github.io/curves/2023/04/18/bezpath-simplify.html
const (
dPenaltyElbow = 0.65
dPenaltySlope = 2.0
)
const numSamples = 20
// FittableCurve describes curves in a way useful for curve fitting. Curves that
// implement this interface can be used as source curves in [FitToBezPath] and
// [FitToBezPathOpt].
//
// The interface can represent source curves with cusps and corners, though if
// the corners are known in advance, it may be better to run curve fitting on
// subcurves bounded by the corners.
//
// The interface primarily works by sampling the source curve and computing the
// position and derivative at each sample. Those derivatives are then used for
// multiple sub-tasks, including ensuring G1 continuity at subdivision points,
// computing the area and moment of the curve for curve fitting, and casting
// rays for evaluation of a distance metric to test accuracy.
//
// A major motivation is computation of offset curves, which often have cusps,
// but the presence and location of those cusps is not generally known. It is
// also intended for conversion between curve types (for example, piecewise
// Euler spiral or NURBS), and distortion effects such as perspective transform.
//
// Note general similarities to [ParametricCurve] but also important
// differences. Instead of separate methods for evaluating the curve and its
// derivative, there is a single SamplePtDeriv method, which can be more
// efficient and also handles cusps more robustly. Also, there is no method for
// subsegmenting, as that is not needed and would be annoying to implement.
type FittableCurve interface {
// SamplePtTangent evaluates the curve and its tangent at parameter t.
//
// For a regular curve (one not containing a cusp or corner), the
// derivative is a good choice for the tangent vector and the sign
// parameter can be ignored. Otherwise, the sign parameter selects which
// side of the discontinuity the tangent will be sampled from.
//
// Generally, t is in the range [0, 1].
SamplePtTangent(t float64, sign float64) CurveFitSample
// SamplePtDeriv evaluates the point and derivative at parameter t.
//
// In curves with cusps, the derivative can go to zero.
SamplePtDeriv(t float64) (Point, Vec2)
// BreakCusp findf a cusp or corner within the given range.
//
// If the range contains a corner or cusp, return it. If there is more
// than one such discontinuity, any can be reported, as the function will
// be called repeatedly after subdivision of the range.
//
// Do not report cusps at the endpoints of the range, as this may cause
// potentially infinite subdivision. In particular, when a cusp is reported
// and this method is called on a subdivided range bounded by the reported
// cusp, then the subsequent call should not report a cusp there.
//
// The definition of what exactly constitutes a cusp is somewhat loose.
// If a cusp is missed, then the curve fitting algorithm will attempt to
// fit the curve with a smooth curve, which is generally not a disaster but
// will usually result in more subdivision. Conversely, it might be useful
// to report near-cusps, specifically points of curvature maxima where the
// curvature is large but still mathematically finite.
BreakCusp(start, end float64) (float64, bool)
}
type MomentIntegraler interface {
MomentIntegrals(start, end float64) (float64, float64, float64)
}
// MomentIntegrals computes moment integrals.
//
// This function computes the integrals of y dx, x y dx, and y^2 dx over the
// length of this curve. From these integrals it is fairly straightforward
// to derive the moments needed for curve fitting.
//
// By default it uses quadrature integration with Green's theorem, in terms of
// samples evaluated with [FittableCurve.SamplePtDeriv]. If pcf implements
// [MomentIntegraler], then this function defers to it.
func MomentIntegrals(pcf FittableCurve, start, end float64) (float64, float64, float64) {
if pcf, ok := pcf.(MomentIntegraler); ok {
return pcf.MomentIntegrals(start, end)
}
t0 := 0.5 * (start + end)
dt := 0.5 * (end - start)
var a, x, y float64
for _, v := range gaussLegendreCoeffs16 {
wi, xi := v[0], v[1]
t := t0 + xi*dt
p, d := pcf.SamplePtDeriv(t)
a_ := wi * d.X * p.Y
x_ := p.X * a_
y_ := p.Y * a_
a += a_
x += x_
y += y_
}
return a * dt, x * dt, y * dt
}
// FitToCubic fits a single cubic to a range of the source curve.
//
// Returns the cubic segment and the square of the error. Returns false if no
// fitting cubic could be computed.
func FitToCubic(
source FittableCurve,
rangeStart float64,
rangeEnd float64,
accuracy float64,
) (CubicBez, float64, bool) {
start := source.SamplePtTangent(rangeStart, 1.0)
end := source.SamplePtTangent(rangeEnd, -1.0)
d := end.Point.Sub(start.Point)
chord2 := d.Hypot2()
acc2 := accuracy * accuracy
if chord2 <= acc2 {
// Special case very short chords; try to fit a line.
return tryFitLine(source, accuracy, rangeStart, rangeEnd, start.Point, end.Point)
}
th := d.Angle()
mod2Pi := func(th float64) float64 {
thScaled := th * (1.0 / math.Pi) * 0.5
return math.Pi * 2.0 * (thScaled - math.Round(thScaled))
}
th0 := mod2Pi(start.Tangent.Angle() - th)
th1 := mod2Pi(th - end.Tangent.Angle())
area, x, y := MomentIntegrals(source, rangeStart, rangeEnd)
x0, y0 := start.Point.X, start.Point.Y
dx, dy := d.X, d.Y
// Subtract off area of chord
area -= dx * (y0 + 0.5*dy)
// area is oriented area of closed curve segment.
// This quantity is invariant to translation and rotation.
// Subtract off moment of chord
dy3 := dy * (1.0 / 3.0)
x -= dx * (x0*y0 + 0.5*(x0*dy+y0*dx) + dy3*dx)
y -= dx * (y0*y0 + y0*dy + dy3*dy)
// Translate start point to origin; convert raw integrals to moments.
x -= x0 * area
y = 0.5*y - y0*area
// Rotate into place (this also scales up by chordlength for efficiency).
moment := d.X*x + d.Y*y
// moment is the chordlength times the x moment of the curve translated
// so its start point is on the origin, and rotated so its end point is on the
// x axis.
chord2Inv := 1.0 / chord2
unitArea := area * chord2Inv
mx := moment * (chord2Inv * chord2Inv)
// unitArea is signed area scaled to unit chord; mx is scaled x moment
chord := math.Sqrt(chord2)
aff := Translate(Vec2(start.Point)).Mul(Rotate(th)).Mul(Scale(chord, chord))
curveDist := curveDistFromCurve(source, rangeStart, rangeEnd)
var bestC maybe.Option[CubicBez]
var bestErr2 maybe.Option[float64]
fits, fitsN := cubicFit(th0, th1, unitArea, mx)
for _, cfit := range fits[:fitsN] {
cand := cfit.cbez
d0 := cfit.d0
d1 := cfit.d1
c := cand.Transform(aff)
if err2, ok := curveDist.evalDist(source, c, acc2); ok {
scaleF := func(d float64) float64 {
return 1.0 + max(d-dPenaltyElbow, 0.0)*dPenaltySlope
}
scale := pow2(max(scaleF(d0), scaleF(d1)))
err2 := err2 * scale
if err2 < acc2 && (!bestErr2.Set() || err2 < bestErr2.Unwrap()) {
bestC = maybe.Some(c)
bestErr2 = maybe.Some(err2)
}
}
}
if bestC.Set() && bestErr2.Set() {
return bestC.Unwrap(), bestErr2.Unwrap(), true
} else {
return CubicBez{}, 0, false
}
}
// cubicFit returns curves matching area and moment, given unit chord.
func cubicFit(th0 float64, th1 float64, area float64, mx float64) ([4]struct {
cbez CubicBez
d0 float64
d1 float64
}, int) {
// Note: maybe we want to take unit vectors instead of angle? Shouldn't
// matter much either way though.
s0, c0 := math.Sincos(th0)
s1, c1 := math.Sincos(th1)
a4 := -9.0 *
c0 *
(((2.0*s1*c1*c0+s0*(2.0*c1*c1-1.0))*c0-2.0*s1*c1)*c0 -
c1*c1*s0)
a3 := 12.0 *
((((c1*(30.0*area*c1-s1)-15.0*area)*c0+2.0*s0-
c1*s0*(c1+30.0*area*s1))*
c0+
c1*(s1-15.0*area*c1))*
c0 -
s0*c1*c1)
a2 := 12.0 *
((((70.0*mx+15.0*area)*s1*s1+c1*(9.0*s1-70.0*c1*mx-5.0*c1*area))*
c0-
5.0*s0*s1*(3.0*s1-4.0*c1*(7.0*mx+area)))*
c0 -
c1*(9.0*s1-70.0*c1*mx-5.0*c1*area))
a1 := 16.0 *
(((12.0*s0-5.0*c0*(42.0*mx-17.0*area))*s1-
70.0*c1*(3.0*mx-area)*s0-
75.0*c0*c1*area*area)*
s1 -
75.0*c1*c1*area*area*s0)
a0 := 80.0 * s1 * (42.0*s1*mx - 25.0*area*(s1-c1*area))
// TODO: "roots" is not a good name for this variable, as it also contains
// the real part of complex conjugate pairs.
roots := make([]float64, 0, 4)
const epsilon = 1e-12
if math.Abs(a4) > epsilon {
a := a3 / a4
b := a2 / a4
c := a1 / a4
d := a0 / a4
if quads, ok := FactorQuarticInner(a, b, c, d, false); ok {
for _, quad := range quads {
qc1, qc0 := quad[0], quad[1]
qroots := polyroot.NewPolynomial(qc0, qc1, 1).Roots(math.Inf(-1), math.Inf(1), 0, nil)
if len(qroots) == 0 {
// Real part of pair of complex roots
roots = append(roots, -0.5*qc1)
} else {
roots = append(roots, qroots...)
}
}
}
} else if math.Abs(a3) > epsilon {
qroots := polyroot.NewPolynomial(a0, a1, a2, a3).Roots(math.Inf(-1), math.Inf(1), 0, nil)
roots = append(roots, qroots...)
} else if math.Abs(a2) > epsilon || math.Abs(a1) > epsilon || math.Abs(a0) > epsilon {
qroots := polyroot.NewPolynomial(a0, a1, a2).Roots(math.Inf(-1), math.Inf(1), 0, nil)
roots = append(roots, qroots...)
} else {
return [4]struct {
cbez CubicBez
d0 float64
d1 float64
}{{
CubicBez{
Pt(0.0, 0.0),
Pt(1.0/3.0, 0.0),
Pt(2.0/3.0, 0.0),
Pt(1.0, 0.0),
},
1.0 / 3.0,
1.0 / 3.0,
}}, 1
}
s01 := s0*c1 + s1*c0
var outN int
var out [4]struct {
cbez CubicBez
d0 float64
d1 float64
}
for _, root := range roots {
var d0, d1 float64
if root > 0.0 {
d := (root*s0 - area*(10.0/3.0)) / (0.5*root*s01 - s1)
if d > 0.0 {
d0, d1 = root, d
} else {
d0, d1 = s1/s01, 0.0
}
} else {
d0, d1 = 0.0, s0/s01
}
// We could implement a maximum d value here.
if d0 >= 0.0 && d1 >= 0.0 {
out[outN] = struct {
cbez CubicBez
d0 float64
d1 float64
}{
CubicBez{
Pt(0.0, 0.0),
Pt(d0*c0, d0*s0),
Pt(1.0-d1*c1, d1*s1),
Pt(1.0, 0.0),
},
d0,
d1,
}
outN++
}
}
return out, outN
}
// CurveFitSample describes a sample point of a curve for fitting.
type CurveFitSample struct {
// A point on the curve at the sample location.
Point Point
// A vector tangent to the curve at the sample location.
Tangent Vec2
}
// Intersect intersects a ray orthogonal to the tangent with the given cubic.
//
// Returns a vector of t values on the cubic.
func (cfs CurveFitSample) Intersect(c CubicBez) ([3]float64, int) {
p1 := c.P1.Sub(c.P0).Mul(3)
p2 := Vec2(c.P2).Mul(3).Sub(Vec2(c.P1).Mul(6)).Add(Vec2(c.P0).Mul(3))
p3 := c.P3.Sub(c.P0).Sub(c.P2.Sub(c.P1).Mul(3))
c0 := c.P0.Sub(cfs.Point).Dot(cfs.Tangent)
c1 := p1.Dot(cfs.Tangent)
c2 := p2.Dot(cfs.Tangent)
c3 := p3.Dot(cfs.Tangent)
roots := polyroot.NewPolynomial(c0, c1, c2, c3).Roots(0, 1, 0, nil)
var out [3]float64
nn := 0
for _, t := range roots {
if t >= 0.0 && t <= 1.0 {
out[nn] = t
nn++
}
}
return out, nn
}
// FitToBezPath generates a Bézier path that fits the source curve.
//
// This function recursively subdivides the curve in half by the parameter when
// the accuracy is not met. That gives a reasonably optimized result but not
// necessarily the minimum number of segments.
//
// In general, the resulting Bézier path should have a [Fréchet distance] less
// than the provided accuracy parameter. However, this is not a rigorous
// guarantee, as the error metric is computed approximately.
//
// This function is intended for use when the source curve is piecewise
// continuous, with the discontinuities reported by the BreakCusp method. In
// applications (such as stroke expansion) where this property may not hold, it
// is up to the client to detect and handle such cases. Even so, best effort is
// made to avoid infinite subdivision.
//
// When a higher degree of optimization is desired (at considerably more runtime
// cost), consider using [FitToBezPathOpt] instead.
//
// See [Fitting cubic Bézier curves] and [Parallel curves of cubic Béziers] for
// an explanation of the approach used.
//
// [Fréchet distance]: https://en.wikipedia.org/wiki/Fr%C3%A9chet_distance
// [Fitting cubic Bézier curves]: https://raphlinus.github.io/curves/2021/03/11/bezier-fitting.html
// [Parallel curves of cubic Béziers]: https://raphlinus.github.io/curves/2022/09/09/parallel-beziers.html
func FitToBezPath(source FittableCurve, accuracy float64) iter.Seq[PathElement] {
return func(yield func(PathElement) bool) {
fitToBezPathRec(source, 0.0, 1.0, accuracy, yield, true)
}
}
func fitToBezPathRec(
source FittableCurve,
start float64,
end float64,
accuracy float64,
yield func(PathElement) bool,
needMove bool,
) (cont, needMoveRet bool) {
// Discussion question: possibly should take endpoint samples, to avoid
// duplication of that work.
startPt := source.SamplePtTangent(start, 1.0).Point
endPt := source.SamplePtTangent(end, -1.0).Point
if startPt.DistanceSquared(endPt) <= accuracy*accuracy {
if c, _, ok := tryFitLine(source, accuracy, start, end, startPt, endPt); ok {
if needMove {
if !yield(MoveTo(c.P0)) {
return false, true
}
}
return yield(CubicTo(c.P1, c.P2, c.P3)), false
}
}
var t float64
if cusp, ok := source.BreakCusp(start, end); ok {
t = cusp
} else if c, _, ok := FitToCubic(source, start, end, accuracy); ok {
if needMove {
if !yield(MoveTo(c.P0)) {
return false, true
}
}
return yield(CubicTo(c.P1, c.P2, c.P3)), false
} else {
// A smarter approach is possible than midpoint subdivision, but would be
// a significant increase in complexity.
t = 0.5 * (start + end)
}
if t == start || t == end {
// infinite recursion, just draw a line
p1 := startPt.Lerp(endPt, 1.0/3.0)
p2 := endPt.Lerp(startPt, 1.0/3.0)
if needMove {
if !yield(MoveTo(startPt)) {
return false, true
}
}
return yield(CubicTo(p1, p2, endPt)), false
}
cont, needMove = fitToBezPathRec(source, start, t, accuracy, yield, needMove)
if !cont {
return false, needMove
}
cont, needMove = fitToBezPathRec(source, t, end, accuracy, yield, needMove)
if !cont {
return false, needMove
}
return true, needMove
}
// tryFitLine tries fitting a line.
//
// This is especially useful for very short chords, in which the standard
// cubic fit is not numerically stable. The tangents are not considered, so
// it's useful in cusp and near-cusp situations where the tangents are not
// reliable, as well.
//
// Returns the line raised to a cubic and the error, if within tolerance.
func tryFitLine(
source FittableCurve,
accuracy float64,
rangeStart float64,
rangeEnd float64,
start Point,
end Point,
) (CubicBez, float64, bool) {
acc2 := accuracy * accuracy
chordLine := Line{start, end}
const shortN = 7
maxErr2 := 0.0
dt := (rangeEnd - rangeStart) / float64(shortN+1)
< 8BD0 div id="LC473" class="react-file-line html-div" data-testid="code-cell" data-line-number="473" style="position:relative"> for i := range shortN {
t := rangeStart + float64(i+1)*dt
p, _ := source.SamplePtDeriv(t)
err2, _ := chordLine.Nearest(p, accuracy)
if err2 > acc2 {
// Not in tolerance; likely subdivision will help.
return CubicBez{}, 0, false
}
maxErr2 = max(err2, maxErr2)
}
p1 := start.Lerp(end, 1.0/3.0)
p2 := end.Lerp(start, 1.0/3.0)
c := CubicBez{start, p1, p2, end}
return c, maxErr2, true
}
// FitToBezPathOpt generates a highly optimized Bézier path that fits the source
// curve.
//
// This function is considerably slower than [FitToBezPath], as it computes
// optimal subdivision points. Its result is expected to be very close to the
// optimum possible Bézier path for the source curve, in that it has a minimal
// number of curve segments, and a minimal error over all paths with that number
// of segments.
//
// In general, the resulting Bézier path should have a [Fréchet distance] less
// than the provided accuracy parameter. However, this is not a rigorous
// guarantee, as the error metric is computed approximately.
//
// [Fréchet distance]: https://en.wikipedia.org/wiki/Fr%C3%A9chet_distance
func FitToBezPathOpt(source FittableCurve, accuracy float64) BezPath {
var cusps []float64
path := BezPath{}
t0 := 0.0
for {
t1 := 1.0
if len(cusps) > 0 {
t1 = cusps[len(cusps)-1]
}
if t, ok := fitToBezPathOptInner(source, accuracy, t0, t1, &path); ok {
cusps = append(cusps, t)
} else {
if len(cusps) > 0 {
t := cusps[len(cusps)-1]
t0 = t
cusps = cusps[:len(cusps)-1]
} else {
break
}
}
}
return path
}
// fitToBezPathOptInner fits a range without cusps.
//
// Returns true and a cusp's location if a cusp was found.
func fitToBezPathOptInner(
source FittableCurve,
accuracy float64,
rangeStart float64,
rangeEnd float64,
path *BezPath,
) (float64, bool) {
if t, ok := source.BreakCusp(rangeStart, rangeEnd); ok {
return t, true
}
var err float64
if c, err2, ok := FitToCubic(source, rangeStart, rangeEnd, accuracy); ok {
err = math.Sqrt(err2)
if err < accuracy {
if rangeStart == 0.0 {
path.MoveTo(c.P0)
}
path.CubicTo(c.P1, c.P2, c.P3)
return 0, false
}
} else {
err = 2.0 * accuracy
}
t0, t1 := rangeStart, rangeEnd
n := uint(0)
var lastErr float64
loop:
for {
n++
switch kind, v := fitOptSegment(source, accuracy, t0, t1); kind {
case FitResultParamVal:
t0 = v
case FitResultSegmentError:
lastErr = v
break loop
case FitResultCuspFound:
return v, true
}
}
t0 = rangeStart
const epsilon = 1e-9
f := func(x float64) result[float64, float64] {
return fitOptErrDelta(source, x, accuracy, t0, t1, n)
}
k1 := 0.2 / accuracy
ya := -err
yb := accuracy - lastErr
var x float64
if k := solveITPFallible(f, 0.0, accuracy, epsilon, 1, k1, ya, yb); k.isOK {
x = k.ok[1]
} else {
return k.err, true
}
pathLen := len(*path)
for i := range n {
var t1 float64
if i < n-1 {
switch kind, v := fitOptSegment(source, x, t0, rangeEnd); kind {
case FitResultParamVal:
t1 = v
case FitResultSegmentError:
t1 = rangeEnd
case FitResultCuspFound:
path.Truncate(pathLen)
return v, true
}
} else {
t1 = rangeEnd
}
c, _, ok := FitToCubic(source, t0, t1, accuracy)
if !ok {
panic("unreachable")
}
if i == 0 && rangeStart == 0.0 {
path.MoveTo(c.P0)
}
path.CubicTo(c.P1, c.P2, c.P3)
t0 = t1
if t0 == rangeEnd {
// This is unlikely but could happen when not monotonic.
break
}
}
return 0, false
}
func measureOneSeg(source FittableCurve, rangeStart, rangeEnd float64, limit float64) (float64, bool) {
_, err2, ok := FitToCubic(source, rangeStart, rangeEnd, limit)
return math.Sqrt(err2), ok
}
type FitResultKind int
const (
// The parameter value that meets the desired accuracy.
FitResultParamVal FitResultKind = iota + 1
// Error of the measured segment.
FitResultSegmentError
// The parameter value where a cusp was found.
FitResultCuspFound
)
func fitOptSegment(source FittableCurve, accuracy float64, rangeStart, rangeEnd float64) (FitResultKind, float64) {
if t, ok := source.BreakCusp(rangeStart, rangeEnd); ok {
return FitResultCuspFound, t
}
missingErr := accuracy * 2.0
var err float64
if v, ok := measureOneSeg(source, rangeStart, rangeEnd, accuracy); ok {
err = v
} else {
err = missingErr
}
if err <= accuracy {
return FitResultSegmentError, err
}
t0, t1 := rangeStart, rangeEnd
f := func(x float64) result[float64, float64] {
if t, ok := source.BreakCusp(rangeStart, rangeEnd); ok {
return result[float64, float64]{err: t}
}
var err float64
f, ok := measureOneSeg(source, t0, x, accuracy)
if ok {
err = f
} else {
err = missingErr
}
return result[float64, float64]{isOK: true, ok: err - accuracy}
}
const epsilon = 1e-9
k1 := 2.0 / (t1 - t0)
if k := solveITPFallible(f, t0, t1, epsilon, 1, k1, -accuracy, err-accuracy); k.isOK {
return FitResultParamVal, k.ok[0]
} else {
return FitResultCuspFound, k.err
}
}
// fitOptErrDelta returns the delta error (accuracy - error of last segment) or
// a cusp.
func fitOptErrDelta(
source FittableCurve,
accuracy float64,
limit float64,
rangeStart, rangeEnd float64,
n uint,
) result[float64, float64] {
t0, t1 := rangeStart, rangeEnd
for range n - 1 {
switch kind, v := fitOptSegment(source, accuracy, t0, t1); kind {
case FitResultParamVal:
t0 = v
// In this case, n - 1 will work, which of course means the error is highly
// non-monotonic. We should probably harvest that solution.
case FitResultSegmentError:
return result[float64, float64]{isOK: true, ok: accuracy}
case FitResultCuspFound:
return result[float64, float64]{err: v}
default:
panic(fmt.Sprintf("invalid kind; %v", kind))
}
}
var err float64
v, ok := measureOneSeg(source, t0, t1, limit)
if ok {
err = v
} else {
err = accuracy * 2
}
return result[float64, float64]{isOK: true, ok: accuracy - err}
}
// curveDist is an acceleration structure for estimating curve distance.
type curveDist struct {
samples [numSamples]CurveFitSample
arcparams []float64
rangeStart float64
rangeEnd float64
// A "spicy" curve is one with potentially extreme curvature variation,
// so use arc length measurement for better accuracy.
spicy bool
}
func curveDistFromCurve(source FittableCurve, rangeStart, rangeEnd float64) curveDist {
step := (rangeEnd - rangeStart) * (1.0 / float64(numSamples+1))
var lastTan maybe.Option[Vec2]
spicy := false
const spicyThresh = 0.2
var samples [numSamples]CurveFitSample
for i := range numSamples + 2 {
sample := source.SamplePtTangent(rangeStart+float64(i)*step, 1.0)
if tan, ok := lastTan.Get(); ok {
cross := sample.Tangent.Cross(tan)
dot := sample.Tangent.Dot(tan)
if math.Abs(cross) > spicyThresh*math.Abs(dot) {
spicy = true
}
}
lastTan = maybe.Some(sample.Tangent)
if i > 0 && i < numSamples+1 {
samples[i-1] = sample
}
}
return curveDist{
samples: samples,
arcparams: nil,
rangeStart: rangeStart,
rangeEnd: rangeEnd,
spicy: spicy,
}
}
func (cd *curveDist) computeArcParams(source FittableCurve) {
const nSubsamples = 10
start, end := cd.rangeStart, cd.rangeEnd
dt := (end - start) * (1.0 / float64((numSamples+1)*nSubsamples))
arclen := 0.0
for i := range numSamples + 1 {
for j := range nSubsamples {
t := start + dt*(float64(i*nSubsamples+j)+0.5)
_, deriv := source.SamplePtDeriv(t)
arclen += deriv.Hypot()
}
if i < numSamples {
cd.arcparams = append(cd.arcparams, arclen)
}
}
arclenInv := 1.0 / arclen
for i := range cd.arcparams {
cd.arcparams[i] *= arclenInv
}
}
// evalArc evaluates distance based on arc length parametrization.
func (cd *curveDist) evalArc(c CubicBez, acc2 float64) (float64, bool) {
// TODO: this could perhaps be tuned.
const epsilon = 1e-9
cArclen := c.PathLength(epsilon)
maxErr2 := 0.0
for i := range cd.samples {
sample := cd.samples[i]
s := cd.arcparams[i]
t := SolveForPathLength(c, cArclen*s, epsilon)
err := sample.Point.DistanceSquared(c.Eval(t))
maxErr2 = max(err, maxErr2)
if maxErr2 > acc2 {
return 0, false
}
}
return maxErr2, true
}
// evalRay evaluates distance to a cubic approximation.
//
// If distance exceeds stated accuracy, can return false. Note that
// acc2 is the square of the target.
//
// Returns the square of the error, which is intended to be a good
// approximation of the Fréchet distance.
func (cd *curveDist) evalRay(c CubicBez, acc2 float64) (float64, bool) {
maxErr2 := 0.0
for _, sample := range cd.samples {
best := acc2 + 1.0
roots, n := sample.Intersect(c)
for _, t := range roots[:n] {
err := sample.Point.DistanceSquared(c.Eval(t))
best = min(best, err)
}
maxErr2 = max(best, maxErr2)
if maxErr2 > acc2 {
return 0, false
}
}
return maxErr2, true
}
func (cd *curveDist) evalDist(source FittableCurve, c CubicBez, acc2 float64) (float64, bool) {
// Always compute cheaper distance, hoping for early-out.
rayDist, ok := cd.evalRay(c, acc2)
if !ok {
return 0, false
}
if !cd.spicy {
return rayDist, true
}
if len(cd.arcparams) == 0 {
cd.computeArcParams(source)
}
return cd.evalArc(c, acc2)
}
0