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];
t1 = PadRight[
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.
