Not quite Boolean operations

[ramcdona]

I have a CSG engine in my program that was developed long before modern ideas of robustness were developed. My code mostly works, but will not suffice for some planned growth. I am considering switching to Geogram as a replacement. This is attractive because it seems that Geogram's approach to CSG is very similar to the one I currently use. This is important because I also use these routines in some non-conventional ways -- and I will need to modify Geogram to behave similarly.

My program (and as I understand Geogram too) approaches CSG operations with the following steps.

1. Two meshes (A, B) are intersected. Some acceleration algorithm (say an octree) is used, but this ends up being a triangle-to-triangle operation.
2. The line of intersection between intersecting triangles is formed (tA, tB).
3. Intersecting triangles are re-triangulated into sub-triangles that obey their parent triangles boundaries and the new intersection line.
4. After intersection is complete, all the triangles are classified -- AinB, AoutB, BinA, BoutA.
5. Depending on the desired operation, different sets of triangles are chosen (say AoutB & BoutA for union).

Sometimes I need to perform Boolean operations where one of my two meshes is open (i.e. it does not enclose a watertight volume). This is not a bad or sloppy surface that is accidentally non-watertight. Instead, this is on purpose -- For example, one mesh is a sphere and the second mesh represents a thin surface. I will always know when I'm dealing with a thin surface, so I can follow a different code path (I don't need Geogram to auto detect this).

The first important difference is in classification -- since one surface is thin, it has no inside and outside, so it can not be used to classify triangles on the other surface. If we were operating on (A,C) where C is a thin surface, the end result of classification would be (A, CinA, CoutA).

There are different operations I might want to perform when doing this intersection -- sometimes I want to keep (A,CinA) or (A,CoutA), or just (CinA) or just (CoutA).

Somewhat similarly, there are times when I want to intersect two meshes (one or both may or may not be open or closed) and I wish to perform steps 1,2,3 -- but then I want to keep all of A or all of B. I.e. I want to 'impress' the intersection line of mesh B on A (and re-triangulate A to satisfy this intersection curve), but then only keep triangles from A.

From poking around a little, the Geogram API is set up for traditional Boolean operations, and it looks like all the pieces I need exist, but I'm not sure how to achieve what I need.

How would I go about achieving these goals with Geogram?

[BrunoLevy]

Hi,

Yes, the algorithm works mostly as described in your outline:

• 1: get all pairs of intersecting triangles using an axis-aligned bounding box tree
• 2: compute the exact intersection of each triangle pair, determine where intersection points are (on vertex, edge, or inside triangle), using this function
• 3: gather all the segments to be inserted in each triangle, and re-triangulate each triangle using a 2D constrained Delaunay triangulation
• 4: duplicate all the surfaces, sort the triangles around the non-manifold edges radial sort
• 5: classify each surface patch. Keep the ones that are on the boundary of the region that corresponds to your boundary expression (that is, the ones for which the result of the boolean expression changes when crossing the patch)
• 6: simplify coplanar facets by determining their boundaries and retriangulate them using (again) a constrained Delaunay triangulation

Everything is done in arbitrary precision. I have two kernels for that, one based on arithmetic expansions (shipped with geogram), but it is a bit slow, and can reach its precision limits with complicated inputs, so I developped a more efficient / more general one based on multi-precision Geogramplus, marketed by the Tessael company.

(I'm currently writing an article about the algorithm, will post it on arXiv when it is ready).

I think that what you suggest is possible, it is mostly a modification in the MeshSurfaceIntersection::classify() function: suppose you have a set of closed surfaces, a boolean expression E (that combines the closed surfaces) and a thin open surface S , and what you want is the portions of the thin open surface for which the boolean expression evaluates as true. Then you can probably use the classification algorithm nearly unchanged, by duplicating the triangles of the thin surface, and creating a new boolean expression E and S, and classifying with respect to it. You will also need to assemble good non-regression tests: believe me, it is super important, it is the first thing to do ! To test my CSG operations, I'm using a large catalog of OpenSCAD models that I downloaded from Thingiverse, plus thingi10K to test the routine that detects and solves self-intersections in meshes. If your test database is large enough, then everything that can happen will happen, and you'll make your code robust, else it is nearly impossible to figure out in advance all the corner cases (at least I'm unable to do so !).

[ramcdona]

Thanks for this great reply. I've spent much of the past day studying classify() -- so far I don't understand the magic of how it does what it does. I'll keep trying.

I appreciate the difference between Geogram and Geogram+ -- though I'm surprised to hear that anything using GMP could be the faster solution (not a dig, just that an arbitrary precision approach comes with a lot of overhead to overcome). I'm interested in your work and eventual paper from an understanding point of view. However, my project is open source and I don't know how I would incorporate a commercial library.

[BrunoLevy]

Everything will become clearer with the article and the figures (it is the type of things that are well explained with a couple of drawings). I also have a powerpoint presentation here that explains some generalities (but it does not have explanations on classify yet unfortunately, but it has a couple of drawings that may help understand this radial sort thing...)

[ramcdona]

I hope you will continue to tolerate some questions while you work on your paper.

I see some images (say on p. 98 of your presentation) that illustrate some thin-surface intersection work. It appears to be related to geophysical fault modeling. While the application is different than my own, it seems closely related. Very encouraging.

One thing that was confusing -- you mentioned complex Boolean expressions -- and BooleanExpression includes a sophisticated boolean expression parser with multiple variables. However, the API appears to only handle simple boolean expressions between two variables.

mesh_union, mesh_intersection, mesh_difference, and mesh_boolean_operation All take two meshes as input and seem to perform a single operation.

This approach would be reasonable (and many programs work this way). It would leave the ordering of repeated evaluations up to the program using the API in order to evaluate more complex Boolean trees.

However, I now see that the lower-level class MeshSurfaceIntersection looks like it can actually compute problems with complex Boolean expressions -- including input of more than two meshes.

First, you would copy multiple meshes into an empty result mesh

copy_operand(result,A,0);
copy_operand(result,B,1);
copy_operand(result,C,2);
copy_operand(result,D,3);
// ...

Then you could continue to use MeshSurfaceIntersection, passing the appropriate expression to classify().

I haven't tested this (yet), but this gives me the idea that one of my needs (to 'impress' one mesh on another, but keep only the first) could be satisfied with a call like....

Mesh A; // Thing to be kept.
Mesh B; // Object whose imprint should be made on A.
Mesh result;
mesh_boolean_operation( result, A, B, string("A"), true);

My hope/expectation is that this would compute the intersections between A and B and then keep all of A, discarding all of B.

[ramcdona]

It seems to me that the operand inclusion bits include all the information to perform my alternative operations. Right now, the only defined operations are (and, or, xor, diff) which correspond to

and: AinB,BinA
or: AoutB,BoutA
xor: (AoutB,-BinA) (BoutA,-AinB)
A-B: AoutB,-BinA

(where - means to flip the normal direction)

My desire is to be able to handle operations like (where A is watertight and B is not)

A,BinA
A,BoutA
BinA
BoutA

One immediate challenge is one of nomenclature / symbology. Currently, the following syntax is accepted for BooleanExpression:

• and: '&' or '*'
• or: '|' or '+'
• xor: '^'
• difference: '-'
• special: '*' for union
• one can use '(' and ')' to group subexpression

What symbols should be used for my 'new' operations? !@#$%<>?:;\ /

Although an expression parser for traditional Boolean ops is pretty concise, I'm afraid that handling all the permutations of my possible operations will lead to a hot mess of a syntax (handled bidirectionally, that is 8 new symbols, yuck).

Instead, I would propose that we use a more abstract syntax based on the in/out operator.

• in: '<'
• out: '>'

Which would allow us to write things like

A, BinA as A+(B<A)

I believe that all possible operations could actually be defined using just +,<,>

Upon further thought, we will also need a symbol to flip the orientation of a surface. Although unary - fits the bill, I think it would be too confusing with the appropriately defined difference operator. I would propose something like / to flip the orientation.

So my revised minimal set of operations would be +, /, <, >

[BrunoLevy]

you got it right, the "low-level" API can do more general n-ary boolean operation than what is exposed in the "high-level" API (union, intersection, difference, with two meshes only), and it works exactly as you have guessed. To test the functions, I have a simplified openSCAD clone here. It has a CSG API that implements (most of) OpenSCAD's operations. For now, to evaluate an OpenSCAD expression, I simply apply the CSG operations in the order I read them from the files. Once it fully works (I still have a few glitches to fix), I plan to implement some tree transformations, then it will use the more general boolean ops that fully use the boolean operation parser. It uses a bit n-ary operations, for union (but union does not use classify()), and for difference (OpenSCAD's difference(A,B,C,D...) corresponds to
differente(A,union(B,C,D....)). More elaborate tree transforms are for the future.

about what you want to do that is, if I understood well, get the part of a possibly open surface that falls into a region defined by a general boolean operation combining several closed surface, my intuition is that the current version of classify() will not do it, because classify() works as follows:

1. assign an operand bit to each input surface.
2. compute the intersections between the surfaces. Sometimes there can be degeneracies with coplanar surfaces, this is why we have an "operand bit" rather than "operand id": whenever a triangle is on two surface, its "operand bitS" (plural) will be a logical OR between the operand bits of each input surface.
3. double all the surfaces. The idea of doubling the surfaces is that each volume has two "skins": an external one, considered to be outside the volume, and an internal one, considered to be inside the volume. It is important for the subsequent steps
4. now what we want is computing the "inclusion bits" for each facet, that is, inclusion bit b of facet f is set if facet f is inside operand b. We do that by using geometry when we have no choice, and combinatorics when we can (combinatorics costs much less than geometry). We need two types of geometric operations: sorting the triangles around the non-manifold edges, and for the "isolated bubbles", determining whether they float in space (new connected component) or whether they represent an internal void within another component. Then "inclusion bits" are propagated using the combinatorial representation (connections between the triangles, that are surfacic connections (alpha2) and volumetric connections (alpha3). This part will be clearer with all the little drawings that I need to create.
5. once you have the inclusion bits, you can find the outer skin of the region determined by the boolean operation as follows: for each pair of triangles t1,t2 that correspond to a duplicated skin, evaluate the boolean operation on the inclusion bits of t1 and t2, if it gives a different result, then it means t1 and t2 are on the boundary of the result of your boolean operation (the boundary of something corresponds to points where you fall inside by moving a little bit in a direction and you fall outside by moving a little bit in the opposite direction, this is what we compute, but combinatorially, using the pairs of triangles on the duplicated surfaces.

Now in your case, you have your set of closed surfaces and the additional open surface. The idea is to proceed as follows:

• insert all the surfaces into a single mesh, using add_operand(). The additional open surface does not need an operand bit
• duplicate all surfaces, as before, except the open surface
• evaluate the inclusion bits everywhere, as before
• evaluate the boolean operation on the triangles of the open surface, keep only the ones that give true

Now can we do the same without writing a new classify() function ? Maybe yes, by cheating a bit: I think that if we insert all the operands and the open surface, assign an operand bit to everybody, and for the open surface, duplicate it as everybody else, but keep the border edges as border edges (so that combinatorially it does not become a "flat volume", it stays a surface), then the algorithm will give you the expected result. We still need to modify classify() a bit for that, but it is a smaller modification.

[ramcdona]

Thanks again for the helpful response.

I have tested a little with a mix of thick and thin surfaces. So far, it gave promising results -- but I wasn't able to explain or understand those results.

Without any changes to Geogram, what do you expect to have happen with open surfaces?

Apparently all surfaces are duplicated to form the outside and inside skins. If there are border edges, they are stitched together to construct a thin volume.

A boolean operation is performed.

What happens with an intentionally open (thin, flat) mesh? What happens with an unintentionally open (say a sphere with a small hole in it) mesh?

It may be that I can achieve what I need with no modifications -- or at least enough to keep moving forward for a while.

[ramcdona]

I need to return to this issue and to make progress. You had previously mentioned an article that was in the works that would help someone like me figure out how to make some of these changes. Has that paper been made available yet?

What do you think is the best way to proceed? Is it better for you to give me some help, guiding me along the way? Or is it better for me to file and Issue report with clearly stated requirements and let you and your team work on it without me?

I am very thankful for your time and effort that you provide so generously -- I do not want to be a burden. At the same time, I need to come to a way to make steady progress on my problem.

[ramcdona]

I've been studying the code quite diligently and am slowly gaining some understanding.

It seems that for non-manifold surfaces, the inner and outer shell get merged together. This is mentioned here, but I see this happen for large holes and for intentionally planar objects.

I don't necessarily see this as a problem, but we need to handle it properly in classify() (actually tentatively_classify_component_vertex_fast()).

Right now, classification gives erratic results when the classifier is non-manifold. This is unsurprising because a point can't actually be inside or outside a non-watertight thing. In the extreme, think of the thing as a flat plane -- in 3D, you can't be inside or outside a plane.

I propose that we either detect non-manifold shapes (perhaps test for number of shells != 2) or let the user identify an object as non-manifold when it is input.

When classifying A in B, for non-manifold B, I propose we always return false (or true, but consistency is the key). We might want to pair this with returning true for B in B, but I'm not sure.