# How to do boolean operations between a large amount of regions

I want to do boolean operations between a large amount of regions.

For convenience, I consider a simplified problem.

The problem is as follows, digging out a large number of small hemispheres on the surface of a hemisphere.

First, i try the function BoundaryDiscretizeRegion

data = Table[{Sqrt[40^2 - r^2] Cos[5 Pi*r/40],
Sqrt[40^2 - r^2] Sin[5 Pi*r/40], r}, {r, 0, 40, 0.2}];
R[0] = BoundaryDiscretizeRegion@
ImplicitRegion[x^2 + y^2 + z^2 <= 40^2 && z >= 0, {x, y, z}];
R[i_] := BoundaryDiscretizeRegion[Ball[data[[i]], 1]];


It works right when when the number of region objects is small.

RegionDifference[R[0], R[1], ViewPoint -> 10 {1, -1, 0.5}]


But when small balls increase to 3 and above, it does not evaluate and it evaluates very slow.

However, i found if we use nested RegionDifference instead of RegionUnion, things become better.

It can evaluate more objects, but also break down when increase to 6 and above.

When i flip through the help documentation, i found a new function CSGRegion is imported.

(CSGR1 =
CSGRegion[
"Difference", {Ball[{0, 0, 0}, 40],
Cuboid[{-40, -40, -40}, {40, 40, 0}]},
ViewPoint -> 10 {1, -1, 0.5}]) // AbsoluteTiming
(CSGR2 =
CSGRegion["Union",
Ball[#, 1] & /@ data[[;; 100]]]) // AbsoluteTiming
CSGRegion["Difference", {CSGR1, CSGR2}, ViewPoint -> 10 {1, -1, 0.5}]


It is very fast and can evaluate much more objects. However, despite taking up only a small amount of memory and very low cpu resources, when the number increases to say 200, Mathematica will crash.

Since in a real problem my sample points will be dense, I also tried, approximating it as a Tube (one object).

RegionDifference[R[0],
DiscretizeGraphics[Tube[data, 1], MaxCellMeasure -> 0.01]]


However, it does not evaluate.

I guest may be Tube is not solid, so i search on MSE for a way to build solid tube, and i finally found in the following question.enter link description here

Here is the code

(*Pixar method;http://jcgt.org/published/0006/01/01/*)
orthogonalDirections[{p1_?VectorQ, p2_?VectorQ}] :=
Module[{s, w, w1, xx, yy, zz}, {xx, yy, zz} = Normalize[p2 - p1];
s = 2 UnitStep[zz] - 1; w = -1/(s + zz); w1 = xx yy w;
{{1 + s w xx^2, s w1, -s xx}, {w1, s + w yy^2, -yy}}]

orthogonalDirections[{p1_?VectorQ, p2_?VectorQ, p3_?VectorQ}] :=
Module[{d, u, v}, {u, v} = Normalize /@ {p3 - p2, p1 - p2};
If[Chop[Norm[u - v] Norm[u + v]] != 0, d = (u + v)/2;
Normalize /@ {d, Cross[u, d]}, orthogonalDirections[{p1, p2}]]]

extend[cs_, q_, d_, nrms_] :=
cs + Outer[Times,
First[
LinearSolve[Transpose[Prepend[-nrms, d]], q - Transpose[cs]]], d]

(*for custom cross-sections*)

crossSection[pointList_?MatrixQ, r_, csList_?MatrixQ] :=
Module[{p1, p2}, {p1, p2} = Take[pointList, 2];
(p1 + #) & /@ (r csList . orthogonalDirections[{p1, p2}])] /;
Last[Dimensions[pointList]] == 3 && Last[Dimensions[csList]] == 2

(*for circular cross-sections*)

crossSection[pointList_?MatrixQ, r_, n_Integer] :=
crossSection[pointList, r,
Composition[Through, {Cos, Sin}] /@ Range[0, 2 Pi, 2 Pi/n]]

makeCap[type : ("Butt" | "Round" | "Square"), s : (-1 | 1),
path_?MatrixQ, tube_?ArrayQ, r_?NumericQ, h_?NumericQ] :=
Module[{d = Take[path, 2 s], cs, p, t0, t1}, cs = tube[[s]];
t0 = h; t1 = 1 - h Boole[type =!= "Square"];
If[s == -1, {t0, t1} = {t1, t0}];
If[type === "Butt", {d[[s]],
Table[ScalingTransform[{t, t, t}, d[[s]]]@cs, {t, t0, t1, s h}]},
p = (s r/(EuclideanDistance @@ d)) ({1, -1} . d);
{d[[s]] + p,
Switch[type, "Round",
Table[
Composition[TranslationTransform[p Cos[\[Pi] t/2]],
ScalingTransform[{1, 1, 1} Sin[\[Pi] t/2], d[[s]]]][cs], {t,
t0, t1, s h}], "Square",
Table[
Composition[TranslationTransform[p],
ScalingTransform[{t, t, t}, d[[s]]]][cs], {t, t0, t1, s h}]]}]]

Options[TubeMesh] = {"CapForm" -> None, "CirclePoints" -> Automatic,
"MeshType" -> Automatic, Tolerance -> Automatic};

TubeMesh[path_?MatrixQ, r_?NumericQ,
opts : OptionsPattern[{TubeMesh, MeshRegion}]] :=
Module[{c0, c1, cf, mt, dims, h, idx, m, n, p0, p1, t0, t1, tol,
tube}, cf = OptionValue["CapForm"]; mt = OptionValue["MeshType"];
If[mt === Automatic,
mt = If[MatchQ[cf, "Butt" | "Round" | "Square"],
BoundaryMeshRegion, MeshRegion]];
tol = OptionValue[Tolerance] /. Automatic -> 0.0015;
n = OptionValue["CirclePoints"];
If[n === Automatic, n = Round[17 tol^(-1/3)/5 - 57 tol^(1/3)/59];
n += Boole[OddQ[n]]]; h = 2/n;
tube =
FoldList[
Function[{p, t},
extend[p, t[[2]], t[[2]] - t[[1]], orthogonalDirections[t]]],
crossSection[path, r, n], Partition[path, 3, 1, {1, 2}, {}]];
If[MatchQ[cf, "Butt" | "Round" | "Square"], {p0, c0} =
makeCap[cf, 1, path, tube, r, h];
{p1, c1} = makeCap[cf, -1, path, tube, r, h];
tube = Join[c0, tube, c1]];
dims = Most[Dimensions[tube]]; tube = Apply[Join, tube];
m = Times @@ dims; idx = Partition[Range[m], Last[dims]];
t0 = t1 = {};
If[MatchQ[cf, "Butt" | "Round" | "Square"], PrependTo[tube, p0];
AppendTo[tube, p1]; idx += 1;
t0 = PadLeft[Partition[First[idx], 2, 1], {Automatic, 3}, 1];
Reverse /@ Partition[Last[idx], 2, 1], {Automatic, 3}, m + 2]];
mt[tube,
Triangle[
Join[t0,
Flatten[
Apply[{Append[Reverse[#1], Last[#2]], Prepend[#2, First[#1]]} &,
Partition[idx, {2, 2}, {1, 1}], {2}], 2], t1]],
FilterRules[{opts}, Options[MeshRegion]]]]


Now the tube is solid

Well done, it give me some results.

But I seem to be too happy.

When i increase the sample points.It throws error.

data = Table[{Sqrt[40^2 - r^2] Cos[18 Pi*r/40],
Sqrt[40^2 - r^2] Sin[18 Pi*r/40], r}, {r, 0, 40, 0.1}];
RegionDifference[R[0], TubeMesh[data, 1, "CapForm" -> "Round"]]


I've tried everything I can, but still can't find a general and stable solution.

Just use the FEM boundary mesh generator:

RegionDifference[
RegionUnion[{Ball[{0, 0, 0}, 40],
Cuboid[{-40, -40, -40}, {40, 40, 0}]}],
RegionUnion[Ball[#, 1] & /@ data]];

Needs["NDSolveFEM"]
bmesh = ToBoundaryMesh[reg, {{-41, 41}, {-41, 41}, {-1, 41}},
"AccuracyGoal" -> 3];

MeshRegion[bmesh]


If the export does not work, use OpenCasdeLink for all of it.

• Wow, it's awesome. After researching this package(OpenCasdeLink) for two days, i was impressed by its powerful capabilities deeply. I have never thought Mathematica could be so good at this. Thank you very much. Apr 11 at 5:37
• Have fun with it ;-) since this is an open source add on you can look at all the code and even contribute... Apr 11 at 5:39

We subdivide the arc according to its arc length.

Here we use NDSolve to do such subdivide as in my previous answers. https://mathematica.stackexchange.com/a/242188/72111

f[t_] = {Sqrt[40^2 - t^2] Cos[5 Pi*t/40],
Sqrt[40^2 - t^2] Sin[5 Pi*t/40], t};
L = NIntegrate[Sqrt[f'[t] . f'[t]], {t, 0, 40}];
t = NDSolveValue[{t'[s]*Norm[f'[t[s]]] == 1, t[0] == 0}, t, {s, 0, L}];
arcs = t /@ Subdivide[0, L, 180];
balls = (r |-> Ball[f[r], 1]) /@ arcs;
CSGR1 = CSGRegion[
"Difference", {Ball[{0, 0, 0}, 40],
Cuboid[{-40, -40, -40}, {40, 40, 0}]},
ViewPoint -> 10 {1, -1, 0.5}];
CSGR2 = CSGRegion["Union", balls]
CSGRegion["Difference", {CSGR1, CSGR2}]


• Thanks very much, but if i modify the 180 to 200, then Mathematica will return two white graphics, and then it will crash. Apr 9 at 4:17
• @hadesth Yes,I have test the two cases. It must be the limitation of CGRRegion in the current version. I also test the Printout3D, too slow to export to STL. Apr 9 at 4:24