Problems with booleans

[tvt173]

Hi @heitzmann ,

I gave it a try using gdstk for boolean operations in a large GDS, however I find multiple issues with it

1. Results are incorrect. I found this even with a simple OR operation on two curvilinear layers. One of the layers in the OR was almost totally omitted, aside from some slivers at the edges. I did export both of the input layers out of gdstk to make sure I was capturing them correctly (I was).
2. It is very slow. I found it to be about 2x slower than gdspy for the same operations on the same example. However, in this case the result out of gdspy was correct.

Sorry, I don't have much I can share yet beyond that, but I just wanted to report for now that there appear to be these issues and make sure you are aware of it. And also see if you might know why we would expect performance to be worse in gdstk.

[heitzmann]

Hi @tvt173, that's very strange. The boolean code is the only part that's exactly the same in both projects, because it comes from a 3rd party library. I'll check the wrapping code to see if I can find any bugs there. If you can send me an example at some point, that's very helpful.

[tvt173]

hi @heitzmann ,

here is my script for the booleans. i had trouble replicating the issue on an example i could share. even when copying the isolated troublesome regions of input geometries, and neglecting the rest of the shapes on the layer, the booleans start to behave fine. but something about having all of the shapes together (even though they are disjoint) causes dropping of certain shapes to happen. this problem does not exist in gdspy.

i keep this script very similar to gdspy syntax so i can easily create an equivalent gdspy script... essentially I only need to replace the polydict with gdspy's cell.get_polygons(by_spec=True)
let me know if there is something wrong fundamentally with this script though

i did also run this on an already-flattened version of the layout to confirm it is not an issue with hierarchy.

def main(gds_in, gds_out, layer_a, layer_b, layer_c):
    start = time()
    print('Reading GDS...')
    lib = gdstk.read_gds(gds_in)
    top = lib.top_level()[0]
    print('Flattening...')
    top = top.flatten()
    polys = top.polygons
    print('Sorting polygons...')
    polydict = {}
    for poly in polys:
        ldt = (poly.layer, poly.datatype)
        if ldt not in polydict:
            polydict[ldt] = []
        polydict[ldt].append(poly)
    layer_a_polys = polydict[layer_a]
    layer_b_polys = polydict[layer_b]
    print('Executing booleans...')
    or_polys = gdstk.boolean(layer_a_polys, layer_b_polys, 'or', layer=layer_c[0], datatype=layer_c[1])
    print('Writing output GDS...')
    lib2 = gdstk.Library('out')
    outcell = lib2.new_cell('top')
    outcell.add(*layer_b_polys)
    outcell.add(*layer_a_polys)
    outcell.add(*or_polys)
    lib2.write_gds(gds_out)
    stop = time()
    elapsed = stop - start
    print(f'Done! {elapsed}s elapsed')

[tvt173]

the problem is also wholly dependent on the layer_a_polys. if I set layer_b_polys instead to [], the same region is still dropped in the OR result

[heitzmann]

@tvt173 Thanks for the example. The issue is replicating the bug without your geometry, but I understand that most layouts we work on are confidential. I played around yesterday trying to come up with large geometries where any bugs would happen and I found 2 instances that I'm working on. Probably latter today I'll push the fixes and maybe you can try them on your layout to see if anything changes.

[tvt173]

thank you. yes, sorry that I cannot share the layout. I tried to clip out the region and distill into something I could share, but like I said, doing that seemed to fix the issue. so indeed it seems to have something to do with the layout being "big". will try the bugfix once you have it

[tvt173]

tried with your new fixes. i seem to be stuck in an infinite loop now. the call to gdstk.boolean never exits

[heitzmann]

The only change in boolean that I ended up implementing is to ask the library to produce only simple polygons at the output. Wihtout it, it would generate self-intersecting polygons and holes. Holes, in particular, are tricky to process and connect to their boundaries.

So the issue is probably in the clipper library. I can try to look into it, but I can't promisse anything at this point.

Alternatively, I cant revert the changes and see if there's any other way to process the self-intersecting holes, but that'll take some time too.

I'll keep posting any updates in this issue.

[tvt173]

hmm ok... i'll try to look also when i get the chance. seems odd that it would be an issue with the clipper library though, since my gdspy clone of this script worked perfectly fine. i feel like there must be some difference in how you are processing the polygons before sending to or after receiving from the clipper lib, no?

[svollenweider]

Hi @heitzmann
I am currently running into the same problem as @tvt173. I was able to reduce one problem to a non-working minimal example and stripped it of all confidential information.
Issue_minimal_nonconfindential.zip
My plan is to union the different paths and structures.
First I convert all of them into polygons and then use either

gdstk::boolean(boxarray,polygons,gdstk::Operation::And,1000,combined);

where boxarray is a box around everything, or any of the below

gdstk::boolean(empty,polygons,gdstk::Operation::Or,1000,combined);

gdstk::boolean(polygons,polygons,gdstk::Operation::Or,1000,combined);

gdstk::boolean(polygons,polygons,gdstk::Operation::And,1000,combined);

All of them seem to yield a hole

between the Path and the Polygon.

[heitzmann]

@svollenweider Thanks for the example! The latest commit fixed the issue there; maybe that was it for the large layout.

It would be quite hard for me to find the problem without the example because I believe all internal functions in gdstk create polygons in the ccw orientation, so, no matter what I tried, it always worked…

@tvt173 Indeed the problem was properly handled in gdspy and I missed the polygon reversal condition when porting. Please, let me know how it goes for you large layout.

[heitzmann]

@svollenweider @tvt173 Did you get a chance of trying out the fixes? I was hoping to release a new version if they are working.

[svollenweider]

Hi @heitzmann
We had a few other priorities to deal with, so I am sorry for the delay. For my code I have to change a few things, since I want/need to get rid of the zlib and lapack dependencies. I just did that and put a few examples into the program. So far not a single one showed any unwanted artifacts.

[heitzmann]

That's great news, @svollenweider thanks! I'll run a few more tests on my own and try a new release soon, then.

[tvt173]

hi @heitzmann , sorry for the delay. i still seem to be getting stuck in an infinite loop when i build the latest master locally and run on my old failing example

[heitzmann]

@tvt173 Sorry to hear that… I've just went through the code but I can't find a where there could be an infinite loop. I've just pushed a minor change to prevent errors if the clipperlib returns degenerate paths (with 0 or 1 vertex) which doesn't seem to happen. Either way, if that was the case, it should result in a segfault, not an infinite loop, so I don't believe it affects your case. I'll keep looking into any differences with gdspy…

If you gather any other info, let me know.

[tvt173]

hi @heitzmann, so I was able to confirm that it was definitely commit 339777 (3397773) which changed the behavior for me from correctly terminating process but incorrect boolean result, to having the process get totally stuck. furthermore, when I sent a SIGINT (CTRL+C) to the process it does not stop. I must forcibly kill it. not sure if that is useful... I have been trying to debug a bit, but not making much progress. sorry, it's been awhile since i've worked with c++, so it's a bit slow-going :)

[tvt173]

hi @heitzmann, noticing a similar problem I believe to be related in gdspy. when I upgrade from 1.6.6 to 1.6.7, I experience crippling performance issues. still trying to figure out if these functions are stalling completely or have just become extremely slow. will let you know if i figure it out. i see this issue come up in multiple different scripts, including my gdspy reference script for this gdstk comparison example.

in the case of this gdspy bug, the source of the issue is almost definitely the new lines around line 4785 of clipper.cpp:

else if ((pnext->Y == p->Y && pprev->Y == p->Y) && ((pnext->X <= p->X && p->X <= pprev->X) || (pprev->X <= p->X && p->X <= pnext->X))) {
          xnew = p->X;
          p1 = pnext;
          break;

what are these lines doing?

[heitzmann]

These lines deal with a special case when the connection point of a hole lies on a horizontal segment of the enclosing polygon. I can't see how this would hang, since there is no loop there. It also assumes the polygon has no self-intersections beyond vertices (that is, no other segment should lie on top of that segment).

[tvt173]

hmm... indeed it does though. my script finishes in about 100 seconds in gdspy 1.6.6, but does not finish even in a full day with gdspy 1.6.7. will have to look into it a bit more i guess

[heitzmann]

Yes, that could be the problem. I've just managed to build an example that seems to suffer from this hanging issue you describe with a distorted checkerboard-like pattern. I'll try to work on it in the next couple of days and keep you posted.

I have to make sure the flag I added to clipperlib is the one responsible for the extra long processing (probably it is, since it also appeared in gdspy) and whether it is really necessary.

[flaport]

It seems to me that the incorrect results part of this issue has been resolved.

Furthermore I could resolve the freezing issue by replacing the following two lines in clipper_tools.cpp:

diff --git a/src/clipper_tools.cpp b/src/clipper_tools.cpp
index 9a319ec..d58646b 100644
--- a/src/clipper_tools.cpp
+++ b/src/clipper_tools.cpp
@@ -218,7 +218,7 @@ ErrorCode boolean(const Array<Polygon*>& polys1, const Array<Polygon*>& polys2,
     ClipperLib::Paths paths1 = polygons_to_paths(polys1, scaling);
     ClipperLib::Paths paths2 = polygons_to_paths(polys2, scaling);

-    ClipperLib::Clipper clpr(ClipperLib::ioStrictlySimple);
+    ClipperLib::Clipper clpr;
     clpr.AddPaths(paths1, ClipperLib::ptSubject, true);
     clpr.AddPaths(paths2, ClipperLib::ptClip, true);

@@ -296,7 +296,7 @@ ErrorCode slice(const Polygon& polygon, const Array<double>& positions, bool x_a
             clip[0][2].Y = clip[0][3].Y = pos;
         }

-        ClipperLib::Clipper clpr(ClipperLib::ioStrictlySimple);
+        ClipperLib::Clipper clpr;
         clpr.AddPaths(subj, ClipperLib::ptSubject, true);
         clpr.AddPaths(clip, ClipperLib::ptClip, true);

You can make a similar replacement in gdspy clipper.cpp to solve the freezing issue there as well.

That said, I'm not sure what I'm changing here. What difference does adding this ioStrictlySimple option make?

[svollenweider]

I assume that with this option you assume that your polygons are stricty simple. This means that vertices or edges are not touching other vertices or edges and the polygon is not self crossing. In a checkerboard pattern, this of course happens.
I would never enable this option if it does what I assume.

[flaport]

right, so removing it (like I did) seems to be a sensible decision?

[heitzmann]

The library documentation is not clear, but I believe the flag refers to the output, not the input. The Clipper object accepts self intersecting algorithms without issues with out without the flag. However, without it, and in very rare cases, it outputs shapes where I cannot properly connect the holes (I assume it is because the output polygons are not simple, which I assume when connecting the holes).

My main issue right now is to properly characterize this situation. I believe the original errors @tvt173 reported fall in this case, but I have yet to find a way to reproduce the issue myself. So far I've been able to reproduce the freezing problem and it is definitely linked to the flag, which means (1) clipperlib seems to have a bug and (2) we should not use the flag (for that, I have to reproduce and fix @tvt173's issues). Or we replace the library altogether…

[svollenweider]

I have also investigated the source code now and compared the c++ version to the delphi version, but could not find any obvious differences. I am also not fully understanding why for example SimplifyPolygons, and all the following functions are not declared static even though they should be.

The flag seems to do call the following codes /functions

 DoSimplePolygons()

 //Line 3053
  if (m_StrictSimple)
  {  
    TEdge* ePrev = e->PrevInAEL;
    if ((e->OutIdx >= 0) && (e->WindDelta != 0) && ePrev && (ePrev->OutIdx >= 0) &&
      (ePrev->Curr.X == e->Curr.X) && (ePrev->WindDelta != 0))
    {
      IntPoint pt = e->Curr;
      #ifdef use_xyz
      SetZ(pt, *ePrev, *e);
      #endif
      OutPt* op = AddOutPt(ePrev, pt);
      OutPt* op2 = AddOutPt(e, pt);
      AddJoin(op, op2, pt); //StrictlySimple (type-3) join
    }
  }


//Line 3150&3164
bool preserveCol = m_PreserveCollinear || m_StrictSimple;
(!preserveCol || !Pt2IsBetweenPt1AndPt3(pp->Prev->Pt, pp->Pt, pp->Next->Pt))))

My main suspect is the Line 3053, since you said that the issue happens with a distorted checkerboard pattern, which would certainly lead to the execution of this line.

It could also be that the issue lies with one of the functions FixupFirstLefts 1/2/3 called in DoSimplePolygons().

[svollenweider]

ErrorInducerMaybe.zip
Does the attached file lead to an error?

[heitzmann]

No, this case seems fine.

[flaport]

Sorry I forgot to mention that me and @tvt173 are working together on this.

For us, the dropped polygon problem seems to be resolved at the current commit (currently d2f151d ) of gdstk. And the freezing problem was resolved by applying my 2-line patch on top of that (i.e. dropping ioStrictlySimple).

I also managed to distill a minimal (non-freezing) example which previously lead to dropped polygons (up until the latest gdstk version 0.6.1, but not anymore at the current commit) with the following or-operation

import gdstk
lib = gdstk.read_gds(infile='test.gds')
cell = lib.top_level()[0].flatten()
out_polys = gdstk.boolean(cell.polygons, [], 'or', layer=0, datatype=0)
out_lib = gdstk.Library("out")
out_cell = gdstk.Cell(name="out")
out_cell.add(*out_polys)
out_lib.add(out_cell)
out_lib.write_gds("out.gds")

working on the following gds file
test.gds.zip:

This is the result of the OR-operation for gdstk version 0.6.1:

This is the (expected) result of the OR-operation for gdstk @ d2f151d:

[heitzmann]

In that case, I'll just revert the changes and leave it in the back burner until I can find a reproducible bug.

[flaport]

Thanks! Would it be possible to make the same change in gdspy as well? (I can also open an issue there if you prefer to track it there as well)

[heitzmann]

@flaport No need to open the issue there. I'll revert it and we'll track it in a single place.

[tvt173]

thanks @heitzmann + all for helping! the booleans in gdspy also appear to be more robust now. will let you know if we notice anything else. on a side note, we can confirm that gdstk is also now significantly faster than gdspy. for some reason though, it seems to be a bit slow on the write_gds operation (?) but still a bit faster overall

[heitzmann]

IO is usually a bottleneck, specially if the hardware is older. Regardless of implementation, writing to disk should be the slow in both Python and C++. There is a chance of small improvements if path conversions and polygon slicing is done in parallel threads, but we'd need a few benchmarks to make sure. This discussion can continue in #8, though.

I'll close this issue for now. If need be, we can re-open it latter.