9
$\begingroup$

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}]

enter image description here

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.

enter image description here

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}]

enter image description here

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.

enter image description here

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

enter image description here

Well done, it give me some results.

enter image description here

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"]]

enter image description here

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

$\endgroup$

2 Answers 2

7
$\begingroup$

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["NDSolve`FEM`"]
bmesh = ToBoundaryMesh[reg, {{-41, 41}, {-41, 41}, {-1, 41}}, 
   "AccuracyGoal" -> 3];

MeshRegion[bmesh]

enter image description here

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

$\endgroup$
2
  • 1
    $\begingroup$ 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. $\endgroup$
    – hadesth
    Apr 11 at 5:37
  • $\begingroup$ Have fun with it ;-) since this is an open source add on you can look at all the code and even contribute... $\endgroup$
    – user21
    Apr 11 at 5:39
4
$\begingroup$

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}]

enter image description here enter image description here

$\endgroup$
2
  • $\begingroup$ Thanks very much, but if i modify the 180 to 200, then Mathematica will return two white graphics, and then it will crash. $\endgroup$
    – hadesth
    Apr 9 at 4:17
  • $\begingroup$ @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. $\endgroup$
    – cvgmt
    Apr 9 at 4:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.