"subtract" gives bad output.

[briansturgill]

Expected Behavior

Actual Behavior

Meshlab reports 12 non manifold vertices 64 faces over non manifold edges. (I've seen this several times, here's a simple example.)

Steps to Reproduce the Problem

Run this: program output to .stl const { cylinder } = require('@jscad/modeling').primitives const { translate } = require('@jscad/modeling').transforms const { subtract, union, intersect } = require('@jscad/modeling').booleans

const main = () => { h = 8 size = 3 const inner_d = size*0.8 const outer_d = 4+size; c1 = cylinder({radius: outer_d/2, height: h}) c2 = cylinder({radius: inner_d/2, height: h}) let shape = union( subtract(c1, c2), //translate([20, 0, 0],union(c1, translate([3, 0, 0],c2))), //translate([0, 20, 0],intersect(c1, translate([3, 0, 0],c2))) ) return shape; }

module.exports = { main }

Run meshlab on the stl file.

Specifications

OpenJSCAD 2.2.21 [Linux:5.15.0-39-generic,linux:x64]

• Version:
• Platform:
• Environment: (browser, local server, node.js etc) CLI

[briansturgill]

OK, put updated code in original issue message. The union also has a similar problem. The intersect works perfectly. Is this by any chance, the "T-junction" problem? Update: geom3.validate(shape) reports it is non-manifold. Prusa slicer doesn't see any errors. I was able to print the faulty subtract object with no issues.

[platypii]

Yea, the issues with boolean operations creating non-manifold geometries is well known. #807 #696 #598 (this comment) and others. Here's the simplest example I've found:

union(
  cube({ size: 8 }),
  cube({ center: [0, 0, 4] }),
)

The main issue appears to be creation of T-junctions as you mention. To fix it will require either:

1. Fix it inside the BSP operations, likely using a data structure which tracks adjacent faces, and anytime one face is split along an edge, the adjacent face must also have a vertex introduced.
2. Fix it after the fact with some kind of "repair" function. We kind of already have this, but it does not seem to work reliably. I suspect this fails due to small rounding errors.

[hrgdavor]

1. sounds faster, and I am not sure if it will be more complicated to make a better fixing tool, or change BSP to track these splits. I can only hope someone with better knowledge of computational geometry than me will find this problem interesting to fix :)

[briansturgill]

@platypii Thanks for that example, I'm sure that will be handy while I work on this!

1. Fix it inside the BSP operations, likely using a data structure which tracks adjacent faces, and anytime one face is split along an edge, the adjacent face must also have a vertex introduced.

I'm starting to question if the right criterion is being use see the union algorithm here.

2. Fix it after the fact with some kind of "repair" function. We kind of already have this, but it does not seem to work reliably. I suspect this fails due to small rounding errors.

I find this very unlikely to work. The boolean algorithms only work on manifold sets of polygons. Anytime you do one of the booleans, there's a good chance you're not manifold. Subsequent boolean ops will have undefined results. Not sure you can ever fix that sort of problem in generalize.

[briansturgill]

I'm busy until next week sometime. I intend to go after the problem full blast when I get free. So if you want to beat me to it, now's the time. :-) I've been hunting down reference materials.

https://www.gamers.org/dhs/helpdocs/bsp_faq.html These folks actually have a Java implementation of 2D and 3D versions in Java... it's not straightforward though. https://commons.apache.org/proper/commons-geometry/tutorials/bsp-tree.html

Here's a paper from 1990 about this... http://www.cs.yorku.ca/~amana/research/bsptSetOp.pdf Alas the scan of the paper is faulty at the edge! EDIT: a better PDF or the article: http://www.vjlove.com/publications/papers/siggraph90.pdf

[briansturgill]

@platypii The union example you gave above actually works, i.e. doesn't fail under either CSCAD or the latest JSCAD.

[platypii]

I should have been more clear. The example I gave DOES produce a non-manifold geometry. But on export, generalize is called and is able to fix this example. There are other examples (like yours) that fail to repair even with generalize. The nice thing about my example was just that it demonstrates the same kind of underlying issue, with nice easy round numbers. This can be shown using geom3.validate:

const jscad = require('@jscad/modeling')
const { geom3 } = jscad.geometries
const { cube } = jscad.primitives
const { union } = jscad.booleans

const main = () => {
  const obs = union(
    cube({ size: 8 }),
    cube({ center: [0, 0, 4] }),
  )
  geom3.validate(obs)
  return obs
}

module.exports = { main }

[platypii]

Looking at @briansturgill's example.

Before generalize, it produces 326 non-manifold edges. After generalize there are zero-area polygons, and 71 non-manifold edges. Even after repair.

nonmanifold.js

const jscad = require('@jscad/modeling')
const { subtract } = jscad.booleans
const { geom3 } = jscad.geometries
const { generalize } = jscad.modifiers
const { cylinder } = jscad.primitives

const main = () => {
  const c1 = cylinder({ radius: 3.5 })
  const c2 = cylinder({ radius: 1.2 })
  let shape = subtract(c1, c2)
  try { geom3.validate(shape) } catch (e) { console.log(e.message) }
  shape = generalize({ triangulate: true }, shape)
  try { geom3.validate(shape) } catch (e) { console.log(e.message) }
  return shape
}

module.exports = { main }

Pre-generalize non-manifold edges:

Post-generalize non-manifold edges:

[briansturgill]

@platypii Thanks!

In a manifold object, every edge has exactly two faces. That edge will get split twice. So I modified SplitPolygonByPlane to use a Dictionary (Map) of edges to returned values and then check for the reverse edge to see if the adjacent was already done, returning the cached intersect point if it was there, so that the points would definitely be the same. It made no difference. So, I just kept number track of the number seen and alternates seen. Over half the alternates were not seen, so clearly something is wrong with the algorithm. I don't think this is strictly a precision problem.

I think I'll try steping through one of your simple examples in the debugger next.

[z3dev]

It's a problem of precision.

@platypii lets investigate those zero area polygons. The generalize functions should be removing those.

[briansturgill]

Ummm... I believe that was what the recently removed repair did. It was generally not turned on anyway.

On Tue, Jun 21, 2022 at 4:36 PM Z3 Development @.***> wrote:

It's a problem of precision.

@platypii https://github.com/platypii lets investigate those zero area polygons. The generalize functions should be removing those.

— Reply to this email directly, view it on GitHub https://github.com/jscad/OpenJSCAD.org/issues/1113#issuecomment-1162430985, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABYCO7ZRTTYHQAQHLH742H3VQI7XXANCNFSM5ZG4SNQA . You are receiving this because you were mentioned.Message ID: @.***>

--

Brian

[z3dev]

The serializers call generalize if necessary to create triangular meshes, etc.

[z3dev]

From the paper on BSP trees...

Note that this code is not robust, since numerical stability may cause errors in the classification of a point. The standard solution is to make the plane "thick" by use of an epsilon value.

Which is the current approach being taken to classify polygons for clipping.

[briansturgill]

I've found a better PDF of the naylor journal article: http://www.vjlove.com/publications/papers/siggraph90.pdf

[briansturgill]

Which is the current approach being taken to classify polygons for clipping.

Yes, but the problem (at least with the Cube, Cube example) is that it is not clipping properly. Or more accurately, edges that should have been split were not. As I said, I think the algorithm implementation is faulty. (It's kind of a myth that there is a BSP union algorithm... generally the academics wave their hand and say it's possible. Note, I only have found only one other package that use BSP's for 3D at all, and that is the Apache Commons project I gave the link to earlier.) None of the big packages (that are open source) use BSP trees.

Certainly the choice of plane for the split is very questionable. It doesn't match at all what was said in: https://www.gamers.org/dhs/helpdocs/bsp_faq.html

[briansturgill]

Anybody at an institution that has access to this: https://www.sciencedirect.com/science/article/abs/pii/S0010448512002412?via%3Dihub It sounds like exactly what we are looking for. Title: "Fast and robust Booleans on polyhedra" Abstract This paper presents efficient methods for supporting robust Boolean operations on pairs of polyhedral manifolds. A new spatial hashing technique is proposed for detecting intersections of pairs of polygon faces quickly and robustly. Robust predicates are used to compute the exact topology of the resulting polyhedron. Degenerate intersections are identified explicitly but handled in the topological routines. Geometric positions of new vertices are then approximated using finite-precision computations. Although vertex locations are rounded off, the exact connectivity graph is kept. The faces of the rounded-off triangulated final model thus could intersect each other and have the potential for an invalid embedding.

[hrgdavor]

Anybody at an institution that has access to this:

I am willing to bu iy for $32 if you are eager to try it out.

[briansturgill]

Anybody at an institution that has access to this:

I am willing to bu iy for $32 if you are eager to try it out.

I too am willing to buy it, but only if it looks like it really does what it is billed to do. Abstracts are unreliable. I was hoping someone could take a look at it first to see if was implementable. Academics are prone to hand waving when it comes to the hard parts.

[briansturgill]

So I think I may have figured out what is wrong with the Booleans algorithm. I think it may be for convex polyhedron. We have convex polygons, but definitely not convex polyhedron [that would be wrong]. Anyway, it would explain why some things are missed. If I am right, we'll have to find a different algorithm. I tried find evanw's email but didn't find it. Perhaps we could ask him about the origins of the algorithm, that might quickly decide the matter.

[hrgdavor]

@briansturgill there are some more in the references this one is may 2022 https://arxiv.org/abs/2205.14151 https://arxiv.org/pdf/2205.14151.pdf

[briansturgill]

OK, I have a "fix", feels more like duct-tape, but here goes. In line 245 of insertTjunctions.js: constants.EPS * constants.EPS, should only be constants.EPS I see this a lot in the code... it makes no sense. First of all, squaring EPS makes it SMALLER. This allows for LESS error. Clearly not what is desired. Didn't check it in JavaScript, but in the CSharp version, everything that was producing bad output, now produces good output. Note, I'm not saying that EPS is the correct value here, just that EPS**2 is too small.

[platypii]

Hm, I'm not seeing any effect from changing that line. On your example with the 2 cylinders, I get exactly the same numbers of non-manifold edges with EPS or with EPS * EPS. I looked at the values we were getting for distancesquared and they were all either like 0 or 0.001 or 1e-32. None were affected by choice of EPS = 0.00001 or EPS * EPS = 1e-10. Maybe I'm missing something?

[briansturgill]

Hm, I'm not seeing any effect from changing that line. On your example with the 2 cylinders, I get exactly the same numbers of non-manifold edges with EPS or with EPS * EPS. I looked at the values we were getting for distancesquared and they were all either like 0 or 0.001 or 1e-32. None were affected by choice of EPS = 0.00001 or EPS * EPS = 1e-10. Maybe I'm missing something? @platypii

Perhaps I wasn't clear enough. You do indeed still get manifold issues, etc. The changes only occur when you output an STL file. It is passing meshlab perfectly. In other words, you need to go through snap/insertTjunctions/triangulate.

I just went back and checked and that one change in insertTjunction fixes this sample: (To be clear, I checked that (EPS*EPS) failed and just EPS works (in meshlab). Git is showing no other changes in code. (Eariler I had many but backed them out.)

I even did a clean and rebuilt from scratch, just to be certain.

var g2 = new Geom2();
Geom3 TapPost(double h, double size, double post_wall = 2, bool add_taper = false, double z_rot = 0)
{
    var inner_d = size * 0.8; //Possibly snug, but with PLA I prefer that
    var outer_d = post_wall * 2 + size;
    var a = ExtrudeLinear(height: h, gobj: Circle(outer_d / 2.0));
    var b = (Geom3)Translate((0, 0, -1), ExtrudeLinear(height: h + 2, gobj: Circle(inner_d / 2.0)));
    a = Cylinder(radius: outer_d / 2.0, height: h);
    b = (Geom3)Translate((0, 0, -1), Cylinder(radius: inner_d / 2.0, height: h + 2));
    //a = Triangulate(a);
    //b = Triangulate(b);
    var gobj = (Geom3)Subtract(a, b);
    //gobj = (Geom3)ExtrudeSimpleOutlines(Circle(outer_d / 2.0), Circle(inner_d / 2.0), height: h);
    //gobj = Triangulate(gobj);
    //gobj.Validate();
    if (add_taper)
    {
        var cout = Circle(outer_d / 2.0);
        var cylout = ExtrudeLinear(height: h * 2, gobj: cout);
        cylout.Validate();
        var cin = Circle(inner_d / 2.0);
        var cylin = ExtrudeLinear(height: h * 2 + 2, gobj: cin);
        cylin.Validate();
        var cb = Cuboid((outer_d, outer_d, h * 3 + 2), center: (0, 0, 0));
        cb.Validate();
        var g3 = (Geom3)Subtract(cylout, Translate((0, 0, -1), cylin));
        //g3.Validate();
        g3 = (Geom3)Subtract(g3, Rotate((45, 0, z_rot), Translate((0, 0, h), cb)));
        //g3.Validate();
        gobj = (Geom3)Union(gobj, (Geom3)Translate((0, 0, -h * 2), g3));
    }
    return gobj;
}

/*
var a = Cylinder(radius: 4.5, height: 8, segments: 10);
a.Validate();
var b = (Geom3)Translate((0, 0, -1), Cylinder(radius: 2, height: 10, segments: 10));
b.Validate();
var g = (Geom3)Subtract(a, b);// TapPost(8, 5, add_taper: false);
*/
var g = TapPost(8, 5, add_taper: true);
Save("/tmp/test.stl", g);
//g.Validate();

[briansturgill]

So, EPS (1e-5) might be too large: (I'm refering to its use in insertTjunctions.js, not elsewhere). For:

g = Union(Sphere(radius: 20), Translate((5, 5, 5), Sphere(radius: 20)));
Save("/tmp/test.stl", g);

I set the value to 1e-6 and get no errors from meshlab. At 1e-5. it causes 1 failure that was different (I thought maybe it "fixed" a T-junction that wasn't a T-junction. At 1e-7 Boundary Faces and Boundary edges exist, though I have no idea what that is or if they are bad. At, 1e-8 and below still sees errors from unfixed T-junctions.

I'm going with 1e-6 at the moment.

The real solutions is to of course not to create T-junctions to begin with. One serious problem with the current algorithm is it doesn't require triangulation. Booleans are much simpler when everything is a triangle.

UPDATE: I did more research on boundary faces and edges, they are bad. None of them should be showing either. https://pages.mtu.edu/~shene/COURSES/cs3621/SLIDES/Mesh.pdf

[briansturgill]

OK, I've tried various values, it looks like 1e-6 works the best. This is all very much a hack though.

[briansturgill]

Ok, after mind numbing searching, I think I've found the general approach that should be used to replace 3D boolean processing. This approach is used in Blender and CGAL and apparently in OpenGL as well. It requires that meshes be triangulated. There are a number of papers: https://cims.nyu.edu/gcl/papers/zhou2016mas.pdf https://arxiv.org/pdf/1601.07953.pdf https://igl.ethz.ch/projects/winding-number/robust-inside-outside-segmentation-using-generalized-winding-numbers-siggraph-2013-jacobson-et-al.pdf https://www.dgp.toronto.edu/projects/fast-winding-numbers/fast-winding-numbers-for-soups-and-clouds-siggraph-2018-barill-et-al.pdf These 4 papers all contain Alec Jacobson (ETH & Univ. Toronto) as a primary/co-author.

[elalish]

Sorry to jump in late; I just found your library. I too am a big fan of OpenSCAD, web development, and computational geometry. I've been building a library that might be helpful on the Mesh Boolean front: https://github.com/elalish/manifold

It's C++, but we now have a WASM build that's <200kb; you can see a really basic test here. I built it from the ground up for topological robustness; I can't promise there are no bugs, but I think it avoids a lot of the problems being talked about here. If anyone wants to kick the tires, I could definitely use more bug reports.

As a bonus, it's also several orders of magnitude faster than OpenSCAD: https://elalish.blogspot.com/2022/03/manifold-performance.html (and it's had some performance improvements since this was written).

[hrgdavor]

@elalish wow, that looks promising indeed. We are definitely interested in making the boolean engine plugabble at some point and also explore what there is. I am also very curious seeing that you also have a gpu variant.

[briansturgill]

Sorry to jump in late; I just found your library. I too am a big fan of OpenSCAD, web development, and computational geometry. I've been building a library that might be helpful on the Mesh Boolean front: https://github.com/elalish/manifold ...

@elalish Interesting, does the algorithm you are using require the meshes to be triangulated. (Our current one doesn't, but I'm having trouble find one to replace that directly.) I'm doing a C# translation of JSCAD, with a Python API Wrapper. C++ is somewhat problematic for me, though I guess I could translate that too. :-)

I don't think you're late to the game at all... currently focused on fixing the 2D problem. But 3D needs done too!

[elalish]

So, yes and no. Manifold's API is based on triangle meshes (I wouldn't trust an input polygonal face to be planar). However, a Boolean operation obviously creates polygonal faces from triangle inputs, so I have internal triangulation and decimation routines to return it to a clean triangle mesh. I invented by own robust triangulator as part of writing this library, since none of the existing ones could guarantee manifoldness and handle rounding errors the way I wanted.

The Boolean itself can operate on polygon meshes and I tried that approach (leave the internal data structure polygonal and only triangulate at the end as needed), but it turns out there are nasty edge cases when input polygons are nearly coplanar and when the polygons have holes. It's much more reliable to triangulate at each step.