Improve B-spline approximations

[gpeairs]

Implements #23.

Using multiple B-splines (and so setting tangents at intermediate points) works much better than creating extra waypoints. Haven't benchmarked time but the goal was to reduce entity counts, and for typical curves this is down to 20 segments from over 100, with typical times around 10ms on my laptop. It's still possible to get very slow (1s) approximations with 50 or 60 segments with tight turns on offsets of long B-splines. That may be because calculations (like approximation error) tend to be easiest and most accurate when the curves are roughly arclength parameterized and have slowly-varying curvature, which isn't so true when the radius of curvature is comparable to the offset.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

📢 Thoughts on this report? Let us know!

[gpeairs]

One tricky thing is getting the curvature of a BSpline offset. Because the underlying BSpline has non-constant curvature, a derivative of that curvature appears in the offset curvature. The current version ignores that term. This calculation is only needed to get the exact curvature of an OffsetSegment, which is only used for choosing discretization points (right now only used to see how often to sample for B-spline approximation). Not to mention that users like to avoid sharp turns, and the default tolerance 1nm is usually overkill anyway. So it seems fine not to be exact.

Of course these are just cubic functions, so we can get a third derivative of the spline itself (constant weights between pairs of waypoints!) and just write down the derivative of the curvature. Not sure it's worth it but we've gotten this far. ForwardDiff would be simple but doesn't play well with Unitful.

(There's also one other missing term still.)

[gpeairs]

BSpline curvature derivative would look like this:

d3_weights(::Interpolations.Cubic, _) = (-1, 3, -3, 1)
function d3r_dt3(r, t)
    n_rescale = (length(r.itp.coefs)-2)-1
    wis = Interpolations.weightedindexes((
        d3_weights,
        ), Interpolations.itpinfo(r)..., (t*n_rescale+1,))
    return Interpolations.symmatrix(
        map(inds -> Interpolations.InterpGetindex(r)[inds...], wis)
    )
end

function dκds(b::BSpline, s)
    # return (signed_curvature(b, s+5nm) - signed_curvature(b, s-5nm))/10nm
    t = arclength_to_t(b, s)
    g = Paths.Interpolations.gradient(b.r, t)[1]
    h = Paths.Interpolations.hessian(b.r, t)[1]
    j = d3r_dt3(b.r, t)[1]

    dκdt = ( # d/dt ((g.x*h.y - g.y*h.x) / ||g||^3)
        (g.x * j.y - g.y * j.x) / norm(g)^3 +
        -3 * (g.x * h.y - g.y * h.x) * (g.x * h.x + g.y * h.y) / norm(g)^5
    )

    return dκdt / norm(g)
end

This agrees very well with finite differences. But it doesn't seem to make the BSpline offset curvature calculation more accurate, so there may be a mistake in how I tried using it (now commented out). As noted above I don't think this is a practical concern, so I will leave it for now.

[gpeairs]

A few other things ended up riding along in the course of implementing/testing the change to BSpline approximation:

• Fix atol and bspline_approximation in SolidModel rendering #46
• Add methods for discretizing segments using curvature, a prerequisite for Use tolerance in Curvilinear and Ellipse rendering #47 and the refactor it mentions. These are only used to test BSpline approximations for now.
• Fix an issue where BSpline discretization tolerance could be violated near inflection points where the curvature is changing too much between samples
• Minor correction to GeneralOffset arclength
• Factor out subsegment method for use with CompoundSegment

[simlapointe]

Code looks good, nice job! I also tested it on more paths from my clipped designs that used to require a lot of points and it worked well for all of them.