Use curvature-based tolerance for ellipse (#96)

* Use curvature-based tolerance for ellipse

* Test SolidModel ellipse keywords

* Fix scaling in discretization_grid

* Fix SchematicLogger deepcopy issue

* Make sure endpoints are not too close together

* Deprecate circle

Author: gpeairs

src/DeviceLayout.jl

@@ -372,6 +372,7 @@ import .Polygons:
     StyleDict,
     addstyle!,
     circle,
+    circle_polygon,
     clip,
     cliptree,
     difference2d,
@@ -395,6 +396,7 @@ export Polygons,
     StyleDict,
     addstyle!,
     circle,
+    circle_polygon,
     clip,
     cliptree,
     difference2d,

src/polygons.jl

@@ -65,6 +65,7 @@ import IntervalSets.endpoints
 
 export Polygon, ClippedPolygon, Ellipse, Circle
 export circle,
+    circle_polygon,
     clip,
     cliptree,
     circularapprox,
@@ -169,14 +170,52 @@ function convert(::Type{Ellipse{T}}, e::Ellipse{S}) where {T, S}
     return Ellipse{T}(convert(Point{T}, e.center), convert(NTuple{2, T}, e.radii), e.angle)
 end
 
-function to_polygons(e::Ellipse; Δθ=5°, kwargs...)
-    r =
-        θ ->
-            (e.radii[1] * e.radii[2]) /
-            sqrt((e.radii[2] * cos(θ))^2 + (e.radii[1] * sin(θ))^2)
-    # θ₀ = atand(e.axis[2], e.axis[1])°
-    p = θ -> e.center + r(θ - e.angle) .* Point(cos(θ), sin(θ))
-    return Polygon(p.(0:Δθ:(360° - Δθ)))
+"""
+    ellipse_curvature(e::Ellipse, θ)
+
+Compute the curvature of an ellipse at parameter θ.
+For an ellipse with semi-major axis `a` and semi-minor axis `b`,
+the curvature is κ(φ) = ab / (a²sin²φ + b²cos²φ)^(3/2)
+where φ is the angle measured from the major axis of the ellipse.
+
+Since θ in the ellipse parameterization is measured from the global x-axis,
+we need φ = θ - e.angle to get the angle relative to the ellipse's major axis.
+"""
+function ellipse_curvature(e::Ellipse, θ)
+    a, b = e.radii[1], e.radii[2]  # a is major axis, b is minor axis
+    φ = θ - e.angle  # Convert from global angle θ to ellipse-relative angle φ
+    return (a * b) / (a^2 * sin(φ)^2 + b^2 * cos(φ)^2)^1.5
+end
+
+function to_polygons(
+    e::Ellipse;
+    atol=DeviceLayout.onenanometer(eltype(e.center)),
+    Δθ=nothing,
+    kwargs...
+)
+    if !isnothing(Δθ) # Use Δθ-based discretization
+        r =
+            θ ->
+                (e.radii[1] * e.radii[2]) /
+                sqrt((e.radii[2] * cos(θ))^2 + (e.radii[1] * sin(θ))^2)
+        p = θ -> e.center + r(θ - e.angle) .* Point(cos(θ), sin(θ))
+        return Polygon(p.(0:Δθ:(360° - Δθ)))
+    else # Use tolerance-based discretization
+        curvature_fn = θ -> ellipse_curvature(e, θ)
+        # t_scale is used to approximately convert "t" (θ) to arclength, use the larger radius to be safe
+        θs = DeviceLayout.discretization_grid(
+            curvature_fn,
+            atol,
+            (0.0, 2π);
+            t_scale=e.radii[1]
+        )
+        r =
+            θ ->
+                (e.radii[1] * e.radii[2]) /
+                sqrt((e.radii[2] * cos(θ))^2 + (e.radii[1] * sin(θ))^2)
+        p = θ -> e.center + r(θ - e.angle) .* Point(cos(θ), sin(θ))
+        return Polygon(p.(θs[1:(end - 1)])) # Don't duplicate last point
+    end
 end
 
 DeviceLayout.magnify(e::Ellipse, mag) = Ellipse(mag .* e.center, mag .* e.radii, e.angle)
@@ -481,11 +520,20 @@ function perimeter(e::Ellipse)
 end
 
 """
-    circle(r, α=10°)
+    circle_polygon(r, Δθ=10°)
 
-Return a circular `Polygon` centered about the origin with radius `r` and angular step `α`.
+Return a circular `Polygon` centered about the origin with radius `r` and angular step `Δθ`.
 """
-circle(r, α=10°) = Polygon(r .* (a -> Point(cos(a), sin(a))).(0:α:(360° - α)))
+circle_polygon(r, Δθ=10°) = Polygon(r .* (a -> Point(cos(a), sin(a))).(0:Δθ:(360° - Δθ)))
+function circle(r, α=10°)
+    @warn """"
+        `circle(r, α)` is deprecated. Use `Circle(r)` or `Circle(center, r)` to create an \
+        exact circle that will be discretized at render time according to rendering keyword \
+        `atol` (default 1nm) or `Δθ` (if provided). To construct the polygon directly, use \
+        `circle_polygon(r, α)`.
+    """
+    return circle_polygon(r, α)
+end
 
 """
     struct Rounded{T <: Coordinate} <: GeometryEntityStyle

src/schematics/schematics.jl

@@ -467,22 +467,36 @@ mutable struct SchematicLogger <: AbstractLogger
     stage::Symbol
     logname::String
     logger::TeeLogger
+end
 
-    function SchematicLogger(name; path="build", level=Logging.Info)
-        logname, logger = if isnothing(path)
-            "output", Base.NullLogger()
-        else
-            mkpath(path)
-            filename = joinpath(path, name * ".log")
-            filename, MinLevelLogger(FormatLogger(_format_log, filename), level)
-        end
-        return new(
-            Dict{Symbol, LogLevel}(),
-            :plan,
-            logname,
-            TeeLogger(current_logger(), logger)
-        )
+function SchematicLogger(name; path="build", level=Logging.Info)
+    logname, logger = if isnothing(path)
+        "output", Base.NullLogger()
+    else
+        mkpath(path)
+        filename = joinpath(path, name * ".log")
+        filename, MinLevelLogger(FormatLogger(_format_log, filename), level)
     end
+    return SchematicLogger(
+        Dict{Symbol, LogLevel}(),
+        :plan,
+        logname,
+        TeeLogger(current_logger(), logger)
+    )
+end
+
+Base.copy(x::SchematicLogger) = SchematicLogger(
+    copy(x.max_level_logged),
+    x.stage,
+    x.logname,
+    TeeLogger(x.logger.loggers...)
+)
+# Needed for an edge case where a `current_logger` was not copyable
+function Base.deepcopy_internal(x::SchematicLogger, stackdict::IdDict)
+    haskey(stackdict, x) && return stackdict[x]
+    y = copy(x)
+    stackdict[x] = y
+    return y
 end
 
 # Overrides required for custom logger

src/solidmodels/render.jl

@@ -94,8 +94,21 @@ function to_primitives(
     return flat
 end
 
-to_primitives(::SolidModel, ent::Ellipse; rounded=true, Δθ=360° / 8, kwargs...) =
-    rounded ? ent : to_polygons(ent; Δθ)
+function to_primitives(::SolidModel, ent::Ellipse; rounded=nothing, Δθ=nothing, kwargs...)
+    if !isnothing(rounded)
+        Base.depwarn(
+            "The `rounded` keyword for Ellipse is deprecated. Use `Δθ=nothing` (default) to keep as ellipse primitive, or `Δθ=some_angle` to discretize to polygon. For the same discretization as `rounded=false` with `Δθ` not specified, use `Δθ=360°/8`",
+            :to_primitives
+        )
+        rounded && return ent  # Keep as ellipse primitive
+        # Otherwise, use the old default Δθ for backward compatibility
+        return to_polygons(ent; Δθ=(isnothing(Δθ) ? 360° / 8 : Δθ), kwargs...)
+    else
+        isnothing(Δθ) && return ent  # Keep as ellipse primitive (default code path)
+        # Otherwise, use specified Δθ
+        return to_polygons(ent; Δθ, kwargs...)  # Discretize to polygon
+    end
+end
 
 # Path nodes that can be drawn with native curves (in OCC)
 # Gmsh does have its own native curves but we don't use them (the APIs and particularly

src/utils.jl

@@ -216,7 +216,7 @@ function discretization_grid(
     bnds::Tuple{Float64, Float64}=(0.0, 1.0);
     t_scale=1.0
 )
-    dt = 0.01
+    dt = 0.01 * bnds[2]
     ts = zeros(4000)
     ts[1] = bnds[1]
     t = bnds[1]
@@ -229,7 +229,7 @@ function discretization_grid(
         )
         t = ts[i - 1]
         # Set dt based on distance from chord assuming constant curvature
-        if cc >= 100 * 8 * tolerance / t_scale^2 # Update dt if curvature is not near zero
+        if cc >= 100 * 8 * tolerance / (bnds[2]^2 * t_scale^2) # Update dt if curvature is not near zero
             dt = uconvert(NoUnits, sqrt(8 * tolerance / cc) / t_scale)
         end
         if t + dt >= bnds[2]
@@ -246,6 +246,9 @@ function discretization_grid(
         ts[i] = min(bnds[2], t + dt)
         t = ts[i]
     end
-
+    # Make sure last two points aren't unnecessarily close together
+    if ts[i - 1] > (ts[i] + ts[i - 2]) / 2
+        ts[i - 1] = (ts[i] + ts[i - 2]) / 2
+    end
     return ts[1:i]
 end

test/coverage/Project.toml

@@ -1,3 +1,4 @@
 [deps]
+DeviceLayout = "ebf59a4a-04ec-49d7-8cd4-c9382ceb8e85"
 LocalCoverage = "5f6e1e16-694c-5876-87ef-16b5274f298e"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

test/test_clipping.jl

@@ -163,7 +163,7 @@ import DeviceLayout.Polygons: circularapprox, circularequality
         @test length(ip) == 1
         @test ip[1] == Polygon(Point{Int}[(2.0, 1.0), (1.0, 1.0), (1.0, 0.0), (2.0, 0.0)])
 
-        c = circle(1, 30°)
+        c = circle_polygon(1, 30°)
         f = θ -> Point(cosd(θ), sind(θ))
         ptrue = f.(0:30:330)
         @test circularapprox(points(c), ptrue)
@@ -293,7 +293,7 @@ import DeviceLayout.Polygons: circularapprox, circularequality
             Point{T}[(2.0μm, 1.0μm), (1.0μm, 1.0μm), (1.0μm, 0.0μm), (2.0μm, 0.0μm)]
         )
 
-        c = circle(1.0μm, 30°)
+        c = circle_polygon(1.0μm, 30°)
         f = θ -> Point(T(cosd(θ)), T(sind(θ)))
         ptrue = f.(0:30:330)
         @test circularapprox(points(c), ptrue)
@@ -330,7 +330,7 @@ import DeviceLayout.Polygons: circularapprox, circularequality
             @test u == union2d([u], r1)
             @test u == union2d(r1, [u])
 
-            c = circle(1, 1°)
+            c = circle_polygon(1, 1°)
             u = union2d(c, c)
             @test u == union2d(c, [c])
             @test u == union2d([c], c)
@@ -426,25 +426,25 @@ import DeviceLayout.Polygons: circularapprox, circularequality
         for CASE = 3:7
             if CASE == 3 # works
                 inner = Rectangle(10μm, 10μm)
-                outer = circle(50μm, 1°)
+                outer = circle_polygon(50μm, 1°)
             elseif CASE == 4 # works
                 inner = Rectangle(10μm, 10μm) |> Translation(0μm, 10μm)
-                outer = circle(50μm, 1°)
+                outer = circle_polygon(50μm, 1°)
             elseif CASE == 5 # DOESN'T work
                 inner = Rectangle(10μm, 10μm) |> Translation(10μm, 0μm)
-                outer = circle(50μm, 1°)
+                outer = circle_polygon(50μm, 1°)
             elseif CASE == 6 # works
                 inner = Rectangle(10μm, 10μm) |> Translation(10μm, 0μm)
-                outer = circle(50μm, 60°)
+                outer = circle_polygon(50μm, 60°)
             elseif CASE == 7 # DOESN'T work
                 inner = Rectangle(10μm, 10μm) |> Translation(10μm, 0μm)
-                outer = circle(50μm, 45°)
+                outer = circle_polygon(50μm, 45°)
             end
             diff = to_polygons(difference2d(outer, inner))
             @test length(diff[1].p) > length(outer.p)
         end
 
-        rout = circle(100μm, π / 50)
+        rout = circle_polygon(100μm, π / 50)
         rin = centered(Rectangle(10μm, 10μm))
         pcb_outline = to_polygons(difference2d(rout, rin))[1]
         @test length(pcb_outline.p) > length(rout.p)

test/test_schematic_solidmodel.jl

@@ -132,7 +132,7 @@ function test_component(component, clip_area, mesh=false, gui=false)
     # Should not render to the solid model.
     render!(
         floorplan.coordinate_system,
-        not_solidmodel(circle(5μm, 1°)),
+        not_solidmodel(circle_polygon(5μm, 1°)),
         SemanticMeta(:circle)
     )
 

test/test_shapes.jl

@@ -13,7 +13,7 @@ rng = MersenneTwister(1337)
     @test length(points((elements(c)[1]))) == 9
 
     # check that polygons aren't rounded if min_angle is set
-    poly = Polygons._round_poly(circle(10μm, 5°), 0.5μm, min_angle=10 / 180 * π)
+    poly = Polygons._round_poly(circle_polygon(10μm, 5°), 0.5μm, min_angle=10 / 180 * π)
     @test length(points(poly)) == 72
 
     @test Polygons._round_poly(Rectangle(10.0, 10.0), 1.0) isa Polygon
@@ -480,6 +480,66 @@ end
 
     @test_nowarn place!(cs, er, GDSMeta())
     @test_nowarn render!(c, cs)
+
+    # Test tolerance-based discretization (new default)
+    e_default = to_polygons(e)  # Should use atol by default
+    e_delta_style = to_polygons(e; Δθ=5°)  # Δθ-based approach
+
+    # The new tolerance-based approach should produce more points than 5° steps
+    @test length(points(e_default)) > length(points(e_delta_style))
+
+    # Test that atol parameter works
+    e_coarse = to_polygons(e; atol=0.1μm)  # Very coarse tolerance
+    e_fine = to_polygons(e; atol=1.0nm)    # Very fine tolerance
+
+    # Coarse tolerance should produce fewer points than fine tolerance
+    @test length(points(e_coarse)) < length(points(e_fine))
+
+    # If tolerance is too high the discretization will get too scared to update dt from 1%
+    # It's OK if an improvement to the discretization algorithm renders this test obsolete
+    # But currently that's the correct thing to do
+    e_too_coarse = to_polygons(e; atol=1.0μm)
+    @test length(points(e_too_coarse)) > length(points(e_coarse))
+
+    # Test backward compatibility - Δθ should still work when explicitly provided
+    e_10deg = to_polygons(e; Δθ=10°)
+    e_5deg = to_polygons(e; Δθ=5°)
+
+    # 10° steps should produce fewer points than 5° steps
+    @test length(points(e_10deg)) < length(points(e_5deg))
+
+    # Test that both atol and Δθ can be provided (Δθ should take precedence for backward compatibility)
+    e_mixed = to_polygons(e; atol=1.0nm, Δθ=10°)
+    @test length(points(e_mixed)) == length(points(e_10deg))
+
+    # Test with circles (special case of ellipse)
+    circ = Circle(Point(0.0μm, 0.0μm), 1.0μm)
+    circ_default = to_polygons(circ)
+    circ_delta = to_polygons(circ; Δθ=5.12°) # default atol gives ~5.12° spacing on this circle
+    circ_delta_coarse = to_polygons(circ; Δθ=10°)
+
+    # Delta method also takes away last point if it's closer than Δθ to the end, so n1=n2+1
+    @test length(points(circ_default)) == length(points(circ_delta)) + 1
+    @test length(points(circ_default)) > length(points(circ_delta_coarse))
+    # Discretization should be very similar to circular_arc
+    @test length(points(circ_default)) == length(circular_arc(2pi, 1.0μm, 1.0nm)) - 1 # last pt duplicated
+
+    # Make sure error is as small as tolerance says
+    area(p) =
+        sum(
+            (gety.(p.p) + gety.(circshift(p.p, -1))) .*
+            (getx.(p.p) - getx.(circshift(p.p, -1)))
+        ) / 2
+    e_fine = to_polygons(e; atol=0.1nm)
+    p = to_polygons(difference2d(e_fine, e_default))[1]
+    @test abs(area(p) / perimeter(p)) < 1nm # on average better than 1nm
+
+    # Last two points are not too close together
+    p = points(to_polygons(e, atol=60nm))
+    @test norm(p[1] - p[end]) > norm(p[end] - p[end - 1]) / 2
+
+    # circle is deprecated
+    @test_logs (:warn, r"deprecated") circle(10.0)
 end
 
 @testset "Sweeping" begin

test/test_solidmodel.jl

@@ -466,6 +466,10 @@ import DeviceLayout.SolidModels.STP_UNIT
     place!(cs, e, SemanticMeta(:test))
     sm = SolidModel("test"; overwrite=true)
     render!(sm, cs)
+    @test SolidModels.to_primitives(sm, e) === e
+    @test SolidModels.to_primitives(sm, e; rounded=true) === e
+    @test length(points(SolidModels.to_primitives(sm, e; rounded=false))) == 8
+    @test length(points(SolidModels.to_primitives(sm, e; Δθ=pi / 2))) == 4
 
     # CurvilinearPolygon
     # A basic, noncurved polygon