Add more comments and tests for bspline approximation

Author: gpeairs

src/paths/segments/bspline_approximation.jl

@@ -1,4 +1,4 @@
-# Arbitrary segment: Discretize according to approximation's curvature
+# Sample points used to estimate approximation error
 function _testvals(f::Paths.Segment{T}, f_approx::Paths.BSpline) where {T}
     l = Paths.pathlength(f)
     tgrid = DeviceLayout.discretization_grid(f, _default_curve_atol(T))
@@ -10,6 +10,7 @@ function _testvals(f::Paths.Segment{T}, f_approx::Paths.BSpline) where {T}
     return tgrid, exact, normal, offsets
 end
 
+# For offset segments, sample non-offset curve, plus normals and offset values
 function _testvals(f::Paths.OffsetSegment{T}, f_approx::Paths.BSpline) where {T}
     l = Paths.pathlength(f)
     tgrid = DeviceLayout.discretization_grid(f, _default_curve_atol(T))
@@ -24,8 +25,7 @@ end
 
 # Offset bsplines use tgrid directly to calculate a bit faster
 function _testvals(f::Paths.GeneralOffset{T, BSpline{T}}, f_approx::Paths.BSpline) where {T}
-    h(t) = Paths.Interpolations.hessian(f.seg.r, t)[1]
-    tgrid = DeviceLayout.discretization_grid(h, _default_curve_atol(T))
+    tgrid = DeviceLayout.discretization_grid(f.seg, _default_curve_atol(T))
     sgrid = t_to_arclength.(Ref(f.seg), tgrid)
     offsets = abs.(getoffset.(Ref(f), sgrid))
     # Don't actually calculate the exact curve, we'll only check that the offsets are correct
@@ -40,8 +40,7 @@ function _testvals(
     f::Paths.ConstantOffset{T, BSpline{T}},
     f_approx::Paths.BSpline
 ) where {T}
-    h(t) = Paths.Interpolations.hessian(f.seg.r, t)[1]
-    tgrid = DeviceLayout.discretization_grid(h, _default_curve_atol(T))
+    tgrid = DeviceLayout.discretization_grid(f.seg, _default_curve_atol(T))
     off = abs(getoffset(f))
     # Don't actually calculate the exact curve, we'll only check that the offset is correct
     exact = f.seg.r.(tgrid)
@@ -58,8 +57,12 @@ function _split_testvals(testvals, seg::Segment{T}) where {T}
     tgrid, exact, normal, offset = testvals
     t_half = _t_halflength(seg)
     idx_h2 = findfirst(t -> t > t_half, tgrid)
-    isnothing(idx_h2) &&
-        error("B-spline approximation of $seg is not converging; try increasing `atol`.")
+    isnothing(idx_h2) && # segment is so short there are no points in test discretization
+        @error """
+        B-spline approximation of $seg failed to converge.
+        Check curve for cusps and self-intersections, which may cause approximation to fail.
+        Otherwise, increase `atol` to relax tolerance.
+        """
     t_h1 = tgrid[1:(idx_h2 - 1)]
     t_h2 = tgrid[idx_h2:end]
     new_tgrid_h1 = t_h1 / t_half
@@ -154,7 +157,11 @@ function bspline_approximation(
     split_tv = [testvals]
     while err > atol
         if refine > maxits
-            @warn "Maximum error $err > tolerance $atol with $(refine-1) refinement iterations. Increase `maxits` to refine further or increase `atol` to relax tolerance."
+            @warn """
+            Maximum error $err > tolerance $atol after $(refine-1) refinement iterations.
+            Check curve $f for cusps and self-intersections, which may cause approximation to fail.
+            Increase `maxits` or manually split path to refine further, or increase `atol` to relax tolerance.
+            """
             break
         end
         err = 0.0 * oneunit(T)

src/paths/segments/compound.jl

@@ -42,7 +42,7 @@ function subsegment(s::CompoundSegment{T}, t) where {T}
         t < l1 && return seg, t - l0
         l0 = l1
     end
-    return last(s.segments), t - l0
+    return last(s.segments), t - l0 + pathlength(last(s.segments))
 end
 
 function curvatureradius(s::CompoundSegment, t)
@@ -66,24 +66,17 @@ function (s::CompoundSegment{T})(t) where {T}
         a0 = p0(c[1])
         x = a0 + Point(D0x * t, D0y * t)
         return x::Point{R}
+    elseif t > L
+        h = c[end]
+        h′ = ForwardDiff.derivative(h, pathlength(h))::Point{Float64}
+        D1x, D1y = getx(h′), gety(h′)
+        a = p1(c[end])
+        x = a + Point(D1x * (t - L), D1y * (t - L))
+        return x::Point{R}
     end
 
-    for i = 1:length(c)
-        seg = c[i]
-        l1 = l0 + pathlength(seg)
-        if t < l1
-            x = (seg)(t - l0)
-            return x::Point{R}
-        end
-        l0 = l1
-    end
-
-    h = c[end]
-    h′ = ForwardDiff.derivative(h, pathlength(h))::Point{Float64}
-    D1x, D1y = getx(h′), gety(h′)
-    a = p1(c[end])
-    x = a + Point(D1x * (t - L), D1y * (t - L))
-    return x::Point{R}
+    seg, dt = subsegment(s, t)
+    return seg(dt)::Point{R}
 end
 
 function _split(seg::CompoundSegment{T}, x, tag1=gensym(), tag2=gensym()) where {T}
@@ -135,10 +128,7 @@ end
 function direction(seg::CompoundSegment{T}, t) where {T}
     l0 = zero(T)
     t < l0 && return α0(seg.segments[1])
-    for se in seg.segments
-        l = l0 + pathlength(se)
-        t < l && return direction(se, t - l0)
-        l0 = l
-    end
-    return α1(seg.segments[end])
+    s, dt = subsegment(seg, t)
+    dt > pathlength(s) && return α1(s)
+    return direction(s, dt)
 end

src/paths/segments/offset.jl

@@ -35,7 +35,6 @@ copy(s::OffsetSegment) = OffsetSegment(copy(s.seg), s.offset)
 getoffset(s::ConstantOffset, l...) = s.offset
 getoffset(s::GeneralOffset{T}) where {T} = l -> uconvert(unit(T), s.offset(l))
 getoffset(s::GeneralOffset, l) = getoffset(s)(l)
-getoffset(s::Segment{T}, l...) where {T} = zero(T)
 offset_derivative(seg::ConstantOffset, s) = 0.0
 offset_derivative(seg::GeneralOffset, s) = ForwardDiff.derivative(seg.offset, s)
 
@@ -90,24 +89,32 @@ function curvatureradius(seg::ConstantOffset, s)
 end
 
 function curvatureradius(seg::GeneralOffset, s)
-    # This is not exactly right, particularly when the curvature is changing?
+    # This one is only approximate for BSpline offsets
+    # let s be seg.seg arclength, let t be seg arclength
+    # tangent(seg, s) defined above is dseg/ds
+    # ||d^2/dt^2 seg(s)|| = ||d/dt tangent(seg, s) ds/dt||
+    # ||d/ds tangent(seg, s)|| ||ds/dt||^2
+    # Where ||ds/dt|| = ||tangent(seg, s)|| = 1 for non-offset segments
+    # But needs a correction for offsets
+    # dt/ds = κ n # Sorry t is tangent for the rest of this comment
+    # dn/ds = -κ t
+    # || κ n + (d2_offset n - κ t d_offset) - d_offset κ t - offset κ^2 n - offset t dκ/ds ||
+
     r = curvatureradius(seg.seg, s)
-    reparam = 1 / (sqrt(1 + offset_derivative(seg, s)^2) * (1 - getoffset(seg, s) / r))
+    reparam = 1 / (sqrt((1 - getoffset(seg, s) / r)^2 + offset_derivative(seg, s)^2))
     base_dir = direction(seg.seg, s)
-    base_tangent = Point(cos(base_dir), sin(base_dir))
+    base_tangent = Point(cos(base_dir), sin(base_dir)) # True for non-offset segments (arclength parameterized)
     base_normal = Point(-base_tangent.y, base_tangent.x)
 
-    dn = -(1 / r) * base_tangent
-    ddn = -(1 / r)^2 * [base_normal.x, base_normal.y]
-
     d_offset = offset_derivative(seg, s)
     d2_offset = ForwardDiff.derivative(s_ -> offset_derivative(seg, s_), s)
     d2_seg =
         (
             (1 / r) * base_normal +
-            ddn * getoffset(seg, s) +
-            2 * dn * d_offset +
-            base_normal * d2_offset
+            d2_offset * base_normal +
+            -2 * (1 / r) * base_tangent * d_offset +
+            -((1 / r) * base_normal) * ((1 / r) * getoffset(seg, s))
+            # - getoffset(seg, s) * base_tangent * dκds # Nonzero for BSplines but basically fine to ignore
         ) * reparam^2
     return sign(r) / norm(d2_seg)
 end

test/test_bspline.jl

@@ -124,27 +124,49 @@ end
         end
     end
     # Relaxed tolerance
-    approx = Paths.bspline_approximation(curv[1].curves[1], atol=100nm)
+    approx = Paths.bspline_approximation(curv[1].curves[1], atol=100.0nm)
     @test length(approx.segments) < 20
     # Non-offset curve approximation
     approx = Paths.bspline_approximation(pa[4].seg)
     @test length(approx.segments) < 9
     # Offset curvatureradius
-    g_fd(c, s, ds=1.0nm) = (c(s + ds / 2) - c(s - ds / 2)) / ds
-    h_fd(c, s, ds=1.0nm) = g_fd(s_ -> g_fd(c, s_, ds), s, ds)
-    curvatureradius_fd(c, s, ds=1.0nm) = begin
+    g_fd(c, s, ds=10.0nm) = (c(s + ds / 2) - c(s - ds / 2)) / ds
+    h_fd(c, s, ds=10.0nm) = g_fd(s_ -> g_fd(c, s_, ds), s, ds)
+    curvatureradius_fd(c, s, ds=10.0nm) = begin
         g = g_fd(c, s, ds)
         h = h_fd(c, s, ds)
         ((g.x^2 + g.y^2)^(3 // 2)) / (g.x * h.y - g.y * h.x)
     end
-    # Arbitrary curvature radius calculation is only approximate...
-    c = curv[1].curves[1]
-    @test abs(curvatureradius_fd(c, 10μm) - Paths.curvatureradius(c, 10μm)) < 100nm
-    c = curv[3].curves[1]
-    @test abs(curvatureradius_fd(c, 10μm) - Paths.curvatureradius(c, 10μm)) < 200nm
-    c = curv[5].curves[1]
+    # For BSplines, curvature radius calculation is only approximate, but not bad
+    c = curv[1].curves[1] # ConstantOffset BSpline
+    @test abs(curvatureradius_fd(c, 10μm) - Paths.curvatureradius(c, 10μm)) < 1nm
+    c = curv[3].curves[1] # GeneralOffset BSpline
+    @test abs(curvatureradius_fd(c, 10μm) - Paths.curvatureradius(c, 10μm)) < 50nm
+    c = curv[5].curves[1] # GeneralOffset Turn
+    # Paths.curvatureradius is exact for Turn offsets, the finite difference is wrong
     @test abs(curvatureradius_fd(c, 10μm) - Paths.curvatureradius(c, 10μm)) < 600nm
-    c = curv[8].curves[1]
+    approx = Paths.bspline_approximation(c, atol=100.0nm)
+    pts = DeviceLayout.discretize_curve(c, 100.0nm)
+    pts_approx = DeviceLayout.discretize_curve(approx, 100.0nm)
+    area(p) =
+        sum(
+            (gety.(p.p) + gety.(circshift(p.p, -1))) .*
+            (getx.(p.p) - getx.(circshift(p.p, -1)))
+        ) / 2
+    p = Polygon([pts; reverse(pts_approx)])
+    @test abs(area(p) / perimeter(p)) < 100nm # It's actually ~25nm but the guarantee is ~< tolerance
+    c = curv[8].curves[1] # ConstantOffset Turn
     @test Paths.curvatureradius(c, 10μm) == sign(c.seg.α) * c.seg.r - c.offset
     @test curvatureradius_fd(c, 10μm) ≈ Paths.curvatureradius(c, 10μm) atol = 1nm
+
+    # Failure due to self-intersection
+    pa2 = Path(Point(0.0, 0.0)nm, α0=90°)
+    bspline!(
+        pa2,
+        [Point(-1000.0μm, 1000.0μm), Point(500.0μm, 500.0μm)],
+        0°,
+        Paths.SimpleCPW(20μm, 10μm)
+    )
+    cps = vcat(pathtopolys(pa2)...)
+    @test_logs (:warn, r"Maximum error") Paths.bspline_approximation(cps[1].curves[1])
 end