Add layerwise and lazy tiled clipping

[gpeairs]

Update of #41 using a lazy iterator over tiles (giving up on healing). Adds layerwise booleans xor2d_layerwise etc applied across pairs of GeometryStructures.

Closes #100 and #30.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.

📢 Thoughts on this report? Let us know!

[gpeairs]

Example with DemoQPU17, taking the xor2d_layerwise of the artwork before saving vs after a save/load roundtrip that snaps cell contents to 1nm grids:

using FileIO
include("DemoQPU17.jl")
schematic, artwork = DemoQPU17.qpu17_demo(savegds=true)
roundtrip_artwork = load("qpu17_demo.gds")["qpu17_demo"] # Everything rounded to 1nm grid
@time result = xor2d_layerwise(artwork, roundtrip_artwork)

123.436378 seconds (17.01 M allocations: 1.358 GiB, 0.61% gc time)
Dict{GDSMeta, ClippedPolygon{Unitful.Quantity{Float64, 𝐋, Unitful.ContextUnits{(nm,), 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}, nothing}}}} with 4 entries:
  GDSMeta(10, 0) => ClippedPolygon{Quantity{Float64, 𝐋, ContextUnits{(nm,), 𝐋, FreeUnits{(nm,), 𝐋, nothing}, nothing}}}(Top-level PolyNode with 0 immediate children.)
  GDSMeta(21, 0) => ClippedPolygon{Quantity{Float64, 𝐋, ContextUnits{(nm,), 𝐋, FreeUnits{(nm,), 𝐋, nothing}, nothing}}}(Top-level PolyNode with 1449 immediate children.)
  GDSMeta(20, 0) => ClippedPolygon{Quantity{Float64, 𝐋, ContextUnits{(nm,), 𝐋, FreeUnits{(nm,), 𝐋, nothing}, nothing}}}(Top-level PolyNode with 1449 immediate children.)
  GDSMeta(1, 2)  => ClippedPolygon{Quantity{Float64, 𝐋, ContextUnits{(nm,), 𝐋, FreeUnits{(nm,), 𝐋, nothing}, nothing}}}(Top-level PolyNode with 15401 immediate children.)

It's actually a pretty cool effect, effectively leaving sliver polygons only on the boundaries between off-grid points:

The junction pattern (layer 10) has a clean XOR because those rectangles are already on an integer grid in artwork. Everything else is affected by rounding, but we can see that it's all within 1nm tolerance:

to_polygons(result[GDSMeta(1,2)]) |> filter(p -> Polygons.area(p) / perimeter(p) > 1nm) # !is_sliver(p)

Polygon{Unitful.Quantity{Float64, 𝐋, Unitful.ContextUnits{(nm,), 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}, nothing}}}[]

We can also do this with tiling:

@time result = xor2d_layerwise(artwork, roundtrip_artwork, max_tile_size=1mm)
@time collected_result = collect(result[GDSMeta(1, 2)])

 0.130062 seconds (604.86 k allocations: 113.544 MiB)
 7.059415 seconds (18.67 M allocations: 1.411 GiB, 6.21% gc time)

The first step is fast because it hasn't done any clipping yet, but the second step (which does the clipping) is also much faster than the non-tiled version. (This is only the metal negative layer, but the results from the other layers are very fast to collect in comparison.)

Let's count up total polygons and non-sliver polygons to see if it's consistent:

total = 0
non_sliver = 0
for tile_clip in collected_result
    polys = to_polygons(tile_clip)
    total += length(polys)
    non_sliver += count(.!(Polygons.is_sliver.(polys)))
end
@show total
@show non_sliver

total = 17779
non_sliver = 0

We get more total polygons because polygons that touch a tile edge are double-counted, but the result is basically equivalent.

[gpeairs]

I don't like how the return type is Dict{Meta, ClippedPolygon} if not tiled and Dict{Meta, Iterator{ClippedPolygon}} if tiled. Should these just be separate functions? Should the non-tiled one be wrapped as a single-element iterable? Should the iterator be collected and combined into one ClippedPolygon (but then you lose some of the advantage of working lazily tile-by-tile)?

Switched it to wrapping the non-tiled result as a single-element vector. The main use case should be comparing full layouts, and that should really be done tiled anyway for performance reasons, but at least this way you can directly swap to the non-tiled version for testing if you're worried about edge artifacts.

Also just use the exact tile size in the keyword argument, don't try to be clever and fit it to the layout.

[parrangoiz]

I'm having trouble identifying what exactly changed in the code that was moved to clipping.jl. Giving it a quick skim, it looks like it wasn't just a copy-paste from polygons.jl, for example some union2d methods look new or changed. Would you mind summarizing what changed exactly as there is no diff to look at?

[parrangoiz]

Is this a heuristic, or does this formula come from somewhere? Also it looks like it isn't used anywhere other than for testing, would it be cleaner to just define this in the test files and leave it out of the API docs?

[gpeairs]

It's a heuristic, but pretty common, I think. What we really want is something like "you could move the boundary by less than atol and zero out the area", and this is basically that with an "on average". I'm missing a factor of 2 compared with the version here, which I've now added in. The idea is that a rectangle with length l and width atol has 2*area/perimeter approximately atol. I think this would be of practical use in the API for filtering geometry diffs.

The other sliver definition I've found is in ArcGIS, which uses a "thinness parameter" T = 4pi A / P^2 (1 for a circle, 0 for a line), but I like having a tolerance with an interpretable dimension.

[parrangoiz]

I see ok, yeah I can see how this holds for very small rectangles, but my guess is that it's a bad approximation for arbitrary polygons with side lengths >> atol. For "round trip XORs" like the ones you did here all residual polygons after the XOR will have this property (side lengths < atol) so this seems like a useful filtering function. But how often are users going to be doing these kinds of round-trip checks? This seems more like something we'd only use during internal testing. This package's API is already huge so IMO we should start raising the bar for what gets added to it.

[gpeairs]

I think it’s surprisingly hard to get false positives, but it’s possible for “flag” polygons with a long sub-nm sliver sticking out from one corner to be counted as “on average” slivers this way. I’ve seen those generated by layout programs (not ours) that don’t make keyhole polygons correctly, so it’s not totally exotic, although a millimeter-long sliver on a square micron of real area would be pretty extreme.

I’ll move it to the test setup block, since it’s used in multiple test items across the suite.

[gpeairs]

I'm having trouble identifying what exactly changed in the code that was moved to clipping.jl. Giving it a quick skim, it looks like it wasn't just a copy-paste from polygons.jl, for example some union2d methods look new or changed. Would you mind summarizing what changed exactly as there is no diff to look at?

Sorry about that. I tried moving clipping back inside polygons.jl to improve the diff but it didn't help, too much had been moved.

The first commit is a clean move of code from polygons.jl to clipping.jl, so if you look at this diff
or git blame for the clipping file it will show what actually changed.

[gpeairs]

Rebased and added handling of empty entity lists to clip_tiled. The diff after the code move is now here.

[parrangoiz]

Why was it necessary to add this for union2d but not the other boolean ops?

[gpeairs]

Only union2d has single-argument methods. (You can call the others with an empty second argument but you'd only do it on purpose for union.)

[parrangoiz]

I'm not convinced that we need to have a brand new interface for layerwise clipping. IIUC current booleans union2d, difference2d, etc. can currently take arguments that are any combination of

1. GeometryEntity
2. Array{<:GeometryEntity}
3. GeometryStructure
4. GeometryReference
5. Pair{<:Union{GeometryStructure, GeometryReference}}

where options (3, 4) do a flattening operation and then take the elements, and option (5) grabs a specific layer and flattens. So they already have rich behavior, and therefore adding an additional layerwise keyword would seem quite reasonable to me. Basically a boolean flag that affects the behavior of options (3, 4). That way, tiled clipping could also be an option there; right now it's only applied to the layerwise methods and I'm not sure why if the performance benefits can be significant.

[gpeairs]

The main reason I didn't do that is that the return type would then depend on the value of a keyword argument, which seems questionable for both interface and type-stability reasons. The current booleans always return a single ClippedPolygon. Tiling should return an iterable over ClippedPolygon (lazy iterator if it's going to be useful), and layerwise operations have to return a Dict{Meta}. I already have non-tiled layerwise operations wrapping Dict values as a vector with a single ClippedPolygon per layer for interface consistency with the tiled version.

I didn't add union2d_tiled etc because I think in practice big (>1000 polygon?) boolean operations tend to be layerwise anyway, but that would be as simple as adding union2d_tiled(...) = clip_tiled(Clipper.ClipTypeUnion, ...) etc if you think it's worth it. I'm reluctant to add it as a keyword option for the same reason as above.

[parrangoiz]

Ok I guess I can see the argument both ways when it comes to complicating the interface so happy to concede on this one. Regarding union2d_tiled etc, that makes sense, right now I can't think of any situations when I've asked for massive bools that are not layerwise so let's just leave that out of scope here and add it later if the need arises?

[parrangoiz]

Would it be possible to automatically compute the tile size and hide this from the user (instead only offering a tiled::Bool keyword to turn it on or off)? This would be a huge simplification to the interface, I suspect the laziness barrier will keep a lot of people from ever accessing this functionality if they have to manually figure out and pass a tile size that is good.

[gpeairs]

Sure, added a tiled keyword that heuristically aims for 100 entities per tile. Left tile_size in as a manual override.

[parrangoiz]

Can we keep these off the public API for now? I don't see the direct value to users, and I also don't like that you have to manually pass an index to findbox if you don't want it to re-index every single time. This feels like a use case for caching (which we can implement later) since based on my cursory reading of the Wikipedia entry on R-Trees, building the trees can be quite expensive.

[gpeairs]

Yeah, makes sense, removed those from exports and HTML docs. Agree that we should avoid giving users a way to accidentally re-index in a hot loop, e.g. looks like findbox for all tiles combined has a cost comparable to building one tree in the tests here.

[parrangoiz]

I noticed that both this and clip_tiled default to PolyFillTypePositive, but the original clip defaults to PolyFillTypeEvenOdd. This is usually masked by the fact that the user-facing methods like union2d actually pass PolyFillTypePositive to clip by default, but if anyone were to ever use the bare clip or clip_* methods they'd get different behavior.

[gpeairs]

Hmm, looks like the clip default is a holdover from Devices.jl. It's basically never what we want, and internally clip is used directly in one place, where it does use the default but the polyfill type doesn't matter. I'd lean towards changing it to Positive, but if you think that's too risky (clip is exported) I could change the clip_tiled default to EvenOdd for consistency, or leave it as is.

[parrangoiz]

I know that for our internal users it wouldn't be a risk (99% sure at least) since I've never seen code that uses clip directly; for external users I don't know. Your call but if it were up to me I'd just change the clip default to Positive (implicitly assuming that that's already the behavior that pretty much everyone's been using through union2d, difference2d, etc.)

By the way, related to this - this is why I have a bias towards not exporting functions unless there's a pretty good reason to. It makes it harder to introduce non-breaking changes down the line.

[gpeairs]

OK, default changed. Fair point on exporting functions. Really clip shouldn't be public at all now that we have union2d etc, that's just another Devices.jl legacy.

[parrangoiz]

I missed this on my first pass, but this sounds quite concerning. I agree this sounds OK for XOR but for difference2d you could easily get errors. Why aren't we clipping entities to their containing tiles to avoid this altogether?

[gpeairs]

I think that should work, it’s just an extra intersection operation now that #196 lets us skip the intermediate interiorcuts. I’ll benchmark to see how expensive it is. If it’s significant I can try only intersecting entities found touching tile edges with findbox, or add a keyword option.

[gpeairs]

OK, intersecting everything with the whole tile doubles the time, but if we only intersect entities touching a tile edge, there's no measurable difference in the QPU17 case. Added logic to do that (always, no keyword option) and added a test covering it.

[parrangoiz]

Nice! Thanks. Yeah this would have been a correctness problem in some cases but I'm pretty sure with this fix it should be fine now.

[parrangoiz]

It seems like there can be a lot of repeated indices here, which will translate to a more expensive membership check in line 2536 below. I tried wrapping the output of mapreduce with Set and it halved the runtime of xor2d_layerwise using a contrived example I constructed (both in the case of auto tile size and manual tile size). It also reduced memory allocation significantly (>2x). Try it on your end and see if you can see a similar improvement.

[gpeairs]

Good call, change applied.