Fix Turn GeneralOffset curvature and rare BSpline discretization tolerance violation

Author: gpeairs

src/paths/segments/offset.jl

@@ -24,7 +24,7 @@ struct ConstantOffset{T, S} <: OffsetSegment{T, S}
 end
 
 """
-    struct ConstantOffset{T,S} <: OffsetSegment{T,S}
+    struct GeneralOffset{T,S} <: OffsetSegment{T,S}
 """
 struct GeneralOffset{T, S} <: OffsetSegment{T, S}
     seg::S
@@ -67,7 +67,7 @@ end
 
 function tangent(s::OffsetSegment, t)
     # calculate derivative w.r.t. s.seg arclength
-    # d/dt (s(t) + offset(t) * normal(t))
+    # d/dt (s.seg(t) + offset(s, t) * normal(s.seg, t))
     base_dir = direction(s.seg, t)
     base_tangent = Point(cos(base_dir), sin(base_dir))
     base_normal = Point(-base_tangent.y, base_tangent.x)
@@ -84,38 +84,46 @@ function signed_curvature(seg::Segment, s)
 end
 
 function curvatureradius(seg::ConstantOffset, s)
-    r = curvatureradius(seg.seg, s)
+    r = curvatureradius(seg.seg, s) # Ignore offset * dκds term
     return r - getoffset(seg, s)
 end
 
 function curvatureradius(seg::GeneralOffset, s)
-    # 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
+    # Only approximate for BSpline offsets (ignores offset * dκds term)
+    # let s (above argument) be seg.seg arclength, let l be seg arclength
+    # g is dseg/ds, defined above as `tangent(seg, s)`
+    # ||d^2/dl^2 seg(s)|| = ||d/dl (g ds/dl)||
+    # ||(d/ds g) (ds/dl)^2 + g d^2s/dl^2||
+    # Where ||ds/dl|| = 1/||g|| = 1 for non-offset segments
     # But needs a correction for offsets
-    # dt/ds = κ n # Sorry t is tangent for the rest of this comment
+    # dt/ds = κ n # t is tangent, n normal, κ curvature of base curve
     # dn/ds = -κ t
-    # || κ n + (d2_offset n - κ t d_offset) - d_offset κ t - offset κ^2 n - offset t dκ/ds ||
+    # ||(κ n + (d2_offset n - κ t d_offset) - d_offset κ t - offset κ^2 n - offset t dκ/ds)(ds/dl)^2 + tangent(seg, s) d^2s/dl^2||
 
     r = curvatureradius(seg.seg, s)
-    reparam = 1 / (sqrt((1 - getoffset(seg, s) / r)^2 + offset_derivative(seg, s)^2))
+    ds_dl = 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)) # True for non-offset segments (arclength parameterized)
     base_normal = Point(-base_tangent.y, base_tangent.x)
 
     d_offset = offset_derivative(seg, s)
     d2_offset = ForwardDiff.derivative(s_ -> offset_derivative(seg, s_), s)
+    d2s_dl2 =
+        (-1 / 2) *
+        ds_dl^3 *
+        (
+            -2 * (offset_derivative(seg, s) / r) * (1 - getoffset(seg, s) / r) +
+            2 * offset_derivative(seg, s) * d2_offset
+        )
+
     d2_seg =
         (
             (1 / r) * base_normal +
             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
+            #- getoffset(seg, s) * base_tangent * dκds(seg.seg, s) # Nonzero for BSplines but basically fine to ignore
+        ) * ds_dl^2 + tangent(seg, s) * d2s_dl2
     return sign(r) / norm(d2_seg)
 end
 
@@ -174,6 +182,8 @@ end
 function arclength(f::ConstantOffset{T, BSpline{T}}, t1::Real=1.0; t0::Real=0.0) where {T}
     check_interval(f, t0, t1)
     l_orig = arclength(f.seg, t1; t0=t0)
+    # Constant offset just changes tangent component of gradient
+    # So we can just add the extra contribution due to curvature
     G = StaticArrays.@MVector [zero(Point{T})]
     H = StaticArrays.@MVector [zero(Point{T})]
     l_extra = f.offset * quadgk(t -> _curvature_arclength!(f.seg, G, H, t), t0, t1)[1]
@@ -183,23 +193,23 @@ end
 # For general offsets we also have to add the contribution from varying offset
 function arclength(f::GeneralOffset{T, BSpline{T}}, t1::Real=1.0; t0::Real=0.0) where {T}
     check_interval(f, t0, t1)
-    # If we unwind the curvature, each offset segment will have length (sqrt(1+x^2) ds)
-    # where x is the derivative of the offset
-    l_orig = quadgk(
-        s -> sqrt(1 + (ForwardDiff.derivative(f.offset, s))^2),
-        t_to_arclength(f.seg, t0),
-        t_to_arclength(f.seg, t1)
-    )[1]
-    # Now add curvature contribution as usual
     G = StaticArrays.@MVector [zero(Point{T})]
     H = StaticArrays.@MVector [zero(Point{T})]
-    l_extra = quadgk(
+    # Each offset segment will have length (sqrt(s′^2 +x^2) ds)
+    # where x is the derivative of the offset (normal component of arclength)
+    # and s′ is the curvature*offset contribution (tangent arclength)
+    l = quadgk(
         t ->
-            _curvature_arclength!(f.seg, G, H, t) * getoffset(f, t_to_arclength(f.seg, t)),
+            sqrt(
+                (
+                    1 -
+                    _curvature!(f.seg, G, H, t) * getoffset(f, t_to_arclength(f.seg, t))
+                )^2 + (offset_derivative(f, t_to_arclength(f.seg, t)))^2
+            ) * dsdt(t, f.seg.r),
         t0,
         t1
     )[1]
-    return l_orig + l_extra # + because of _curvature_arclength! sign convention
+    return l
 end
 
 function _curvature_arclength!(b::BSpline, G, H, t)
@@ -211,6 +221,13 @@ function _curvature_arclength!(b::BSpline, G, H, t)
     return uconvert(NoUnits, (G[1].y * H[1].x - G[1].x * H[1].y) / norm(G[1])^2)
 end
 
+function _curvature!(b::BSpline, G, H, t)
+    Paths.Interpolations.gradient!(G, b.r, t)
+    Paths.Interpolations.hessian!(H, b.r, t)
+    return uconvert(NoUnits, (G[1].x * H[1].y - G[1].y * H[1].x) / norm(G[1])^2) /
+           norm(G[1])
+end
+
 function _split(seg::ConstantOffset, x)
     segs = split(seg.seg, x)
     return offset(segs[1], seg.offset), offset(segs[2], seg.offset)

src/utils.jl

@@ -204,10 +204,12 @@ function discretization_grid(s::Paths.Segment, tolerance)
 end
 
 function discretization_grid(s::Paths.BSpline, tolerance)
+    # Assume ds/dt ≈ 1 to avoid converting between arclength and t
     h(t) = Paths.Interpolations.hessian(s.r, t)[1]
     return discretization_grid(h, tolerance)
 end
 
+# Discretize using marching algorithm based on Hessian or curvature
 function discretization_grid(
     ddf,
     tolerance,
@@ -219,20 +221,28 @@ function discretization_grid(
     ts[1] = bnds[1]
     t = bnds[1]
     i = 1
+    cc = norm(ddf(t))
     while t < bnds[2]
         i = i + 1
         i > length(ts) && error(
             "Too many points in curve for GDS. Increase discretization tolerance or split the curve."
         )
         t = ts[i - 1]
-        cc = norm(ddf(t))
-        if cc >= 1e-9 * oneunit(typeof(cc))
+        # Set dt based on distance from chord assuming constant curvature
+        if cc >= 1e-9 * oneunit(typeof(cc)) # Update dt if curvature is not near zero
             dt = uconvert(NoUnits, sqrt(8 * tolerance / cc) / t_scale)
-            # Also assumes third derivative is small
         end
         if t + dt >= bnds[2]
             dt = bnds[2] - t
         end
+        # Check that curvature didn't increase too much (decrease is fine)
+        # Rare but may happen near inflection points
+        cc_next = norm(ddf(t + dt))
+        if (t_scale * dt)^2 * (cc_next - cc) / 24 > tolerance
+            dt = uconvert(NoUnits, sqrt(8 * tolerance / cc_next) / t_scale)
+            cc_next = norm(ddf(t + dt))
+        end
+        cc = cc_next
         ts[i] = min(bnds[2], t + dt)
         t = ts[i]
     end

test/test_bspline.jl

@@ -136,15 +136,15 @@ end
         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
+    end # assumes constant d(arclength)/ds
     # 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
+    # Paths.curvatureradius is exact for Turn offsets
+    @test abs(curvatureradius_fd(c, 10μm) - Paths.curvatureradius(c, 10μm)) < 1nm
     approx = Paths.bspline_approximation(c, atol=100.0nm)
     pts = DeviceLayout.discretize_curve(c, 100.0nm)
     pts_approx = DeviceLayout.discretize_curve(approx, 100.0nm)