1
$\begingroup$

Some time ago I created a printable STL-file with a certain software of my own and now I want to modify it with Mathematica programmatically, mainly by subtracting certain geometric objects (e.g. spheres), for example to get a net -like appearance of what was formerly a massive surface.

As a "warm-up" I imported the surface, repaired it to be sure by RepairMesh, and tried to subtract a sphere centered in the origin of radius 1 with RegionDifference.

But I could not get Mathematica do it, trying a lot of variations of using DiscretizeRegion, BoundaryDiscretizeRegion, TriangulateMesh, BoundaryMesh on the arguments of RegionDifference and on its result value. Nothing worked.

Has anyone an idea what else I could try?

Follows an exported png from part of my notebook:

enter image description here

At last a "meta-question": To insert a part of a notebook into a question (or answer) here, is it the right/best way to select the part of the notebook intended for publication and then use File>Save Selection As from Mathematica's "File"-Menu?

As requested I include the above as plain text:

myQuartic = Import["http://www.aviduratas.de/wc/model-scaled-x7.stl", "BoundaryMeshRegion"]
myQuartic1 = RegionResize[myQuartic, {16.0}]
myQuartic1Repair = RepairMesh[myQuartic1]
(* *)

dcr[reg_, pg_:"Quality", mcm_:0.01] := DiscretizeRegion[reg, PerformanceGoal->pg, MaxCellMeasure->mcm]
bdcr[reg_, pg_:"Quality", mcm_:0.01] := BoundaryDiscretizeRegion[reg, PerformanceGoal->pg, MaxCellMeasure->mcm]

DiscretizeRegion[RegionDifference[myQuartic1Repair, bdcr[Ball[{0,0,0},1]]],Table[{-8,8},3], PerformanceGoal->"Quality", MaxCellMeasure->0.01]
(* During evaluation of In[…]: BoundaryMeshRegion::bsuncl: The boundary surface is not closed because the edges Line[{{126567,126563},{254562,254563},{124643,126266},{268387,268383},{255259,255243},{265372,265373},{194888,194887},{254318,254319},{126273,125167},{124916,124913},{97259,97260},{126449,126450},{124908,124913},{193478,193608},<<23>>,{97283,97281},{205986,205978},{254808,254822},{265338,265350},{126409,126417},{125270,125271},{124936,124940},{194837,193519},{255481,255485},{268500,268499},{125367,125383},{124639,124640},{255245,255247},<<1195>>}] only come from a single face. *)
(* During evaluation of In[…]: DiscretizeRegion::drf: DiscretizeRegion was unable to discretize the region BooleanRegion[<<2>>]. *)
(* DiscretizeRegion[BooleanRegion[#1&&!#2&,{,}],{{-8,8},{-8,8},{-8,8}},PerformanceGoal->Quality,MaxCellMeasure->0.01] *)

In the meantime I tried to overcome the problem by defining explicitly an indicator function for regions and a rule to compute this indicator function for region differences, but it did not solve the problem:

Clear[ff]

Clear[charfun]

charfun[reg_?RegionQ,{x_,y_,z_}] := If[RegionMember[reg,{x,y,z}] == True, 1, 0]
regdiff /: charfun[regdiff[reg1_, reg2_], {x_,y_,z_}] := charfun[reg1, {x,y,z}]  (1-charfun[reg2,{x,y,z}])
Clear[ff]
ff[{x_?NumericQ, y_?NumericQ, z_?NumericQ}] := 
charfun[regdiff[myQuartic1Repair,Ball[{0,0,0},2]], {x,y,z}]
ff[{0.,0.,0.}]
(* 0 *)

charfun[regdiff[myQuartic1Repair,Ball[{0,0,0},2]],{x,y,z}]
(* (1-If[((x|y|z)∈\[DoubleStruckCapitalR]&&x^2+y^2+z^2<=4)==True,1,0]) If[RegionMember[,{x,y,z}]==True,1,0] *)

data=RandomPoint[Ball[{0,0,0},8], 100000];
Graphics3D[Point[data]]
Graphics3D[Point[Select[data, ff[#] == 1 &]]]
dcr[ImplicitRegion[Abs[ff[{x,y,z}]-1] < 0.2, {x,y,z}],Table[{-8,8},3]]
(* During evaluation of In[…]: DiscretizeRegion::drf: DiscretizeRegion was unable to discretize the region ImplicitRegion[Abs[-1+ff[{x,y,z}]]<0.2,{x,y,z}]. *)
dcr[ImplicitRegion[ff[{x,y,z}] >= 0.5, {x,y,z}]]
(* During evaluation of In[…]: DiscretizeRegion::drf: DiscretizeRegion was unable to discretize the region ImplicitRegion[ff[{x,y,z}]>=0.5,{x,y,z}].
(* DiscretizeRegion[ImplicitRegion[ff[{x,y,z}]>=0.5,{x,y,z}],PerformanceGoal->Quality,MaxCellMeasure->0.01] *)

The same as picture:

using a self-defined indicator function

How can I find out what the source of the problem is? I am now even pondering the possibility to do the intersection with an external function call either to a Python wrapper of the CGAL library or to a C++ program calling CGAL. But I would really like to avoid going so far.

Final(?) Note added: I finally did overcome my fear of complicated docker images and python pip dependency errors, installed docker and then used the package PyMesh

PyMesh

to compute what I tried in vain with Mathematica, namely subtracting a centered sphere from my STL-model. No problem, one second computation time, everything done. Import, export: Easy, all features of numpy, scipy etc. available for working on the mesh.

It seems an old prejudice against Mathematica seems to be confirmed here: Great for demos, but if you go to a specific field there are usually much more powerful systems created by and for specialists in this field.

In commutative algebra, you would use Macaulay 2 or Singular or Magma, in group theory GAP or Magma, in scientific computation in general python and Matlab seem to prevail and now, I must concede, in geometric computation, pymesh, with its access to libraries like CGAL, easily overcomes Mathematica.

Really sad, that such a well thought out system like Mathematica with such a great, solid and easy notebook interface and all in all really good language can not deliver the performance its bold advertising lets the users expect.

Improve this question
$\endgroup$
1
  • 4
    $\begingroup$ Please do not post images of your work, especially when the images display at a size that make them difficult to read. Please post your actual Mathematica code in the form of text that can be copied and pasted into a Mathematica notebook. Without such, it will be difficult to reproduce your problem and to experiment with possible solutions. $\endgroup$
    – m_goldberg
    Commented Nov 8, 2019 at 17:42

2 Answers 2

Reset to default
4
$\begingroup$

I have described integrating Mathematica and the open source 3D modeling tool, Blender 2.79b, in previous answers here and here.

Your geometry does have some small features and concavity that can cause many meshers problems. Blender appears to be able to handle it. You will need to learn some python scripting to facilitate the integration but there are plenty of online resources.

Here is a blender workflow to do the boolean difference (you can Uncompress the string to view the python string template.

blenderworkflow[] := 
 Module[{pre, imgset, nbd, differenceScript, file, fileName, 
   outputfile, stext, diffstl},
  nbd = NotebookDirectory[];
  pre = StringTemplate[
     Uncompress[
      "1:eJylU01r20AQzaE/\
ZKoeJIFQkx4LOiSNTA11YmzTEtLgSNYYbyPtqLvrNv732VmtPnIoBCqQ0M6bjzfzZt+\
XtFrfvTs7E01LykDZnn6eX3zaK2qgLcyhFiV4aGmPjDmcFBAIyf5pVZgipfIX7oz+\
zCjYR+yBUnNqEbIMwkW+/hoCB422fLHc3IVDAD+UaqxtGshgo47YQ1hr/JffrLBgT+\
sDXFu7QSjqGjwjZmkOCHqH0jkyZWq1Z5xWLiKKxxxrNFAJZUFSJzAEkgyWRE+\
jtZtBXaGyFHgwUfgoSws/\
hpNEl1XlSpeWjouQRWN7h4ANqTZ1MOXTjXnboD4wFu1FjSxBpo2KfLGP4HLErghnsdk4fk\
fS4LPxY8Fq65u/P3/oPdO++JU9BCPJeacu8/x9LJQRu1dUG7ID+n+u/\
LpUbyPcVe05LPg0ULZK1BfW3PnYr9gLVDqV+Dfi5cqCq9vbb/\
nlTZBwgixwAUE8xHrls14Wb2xRFUaQzMLr+WyWr/\
KbL3nIuOf6mjkDbqWseWimXzH3kx5bezP8ZnVk+aYMzA25+UXOOXELn0C4XOXf5/\
mPcBI1tqiwoT8YOcLDWG032nOb3sQe6oP4/74T/\
yGBirZHWQv5lHHdeKotPr9F22AjGquL24w4OWrcdlLyAIeU/L4AZAty1w=="]][<|
     "nbdir" -> StringReplace[nbd, "\\" -> "\\\\"]|>];
  differenceScript = pre;
  fileName = "blenderboolediff.py";
  file = OpenWrite[fileName];
  WriteString[file, differenceScript];
  Close[file];
  outputfile = CreateFile[];
  Run["blender --background --python blenderboolediff.py >>" <> 
    outputfile];
  stext = OpenRead[outputfile];
  Close[stext];
  DeleteFile[outputfile];
  diffstl = Import["TimMod.stl"];
  diffstl
  ]

Now, you can export Mathematica objects as STL and difference with Blender using the following code (you will need to save the notebook first and it took 3 seconds to execute on my machine):

nbdir = NotebookDirectory[];
SetDirectory@nbdir;
myQuartic = 
 Import["http://www.aviduratas.de/wc/model-scaled-x7.stl", 
  "BoundaryMeshRegion"]
ball = BoundaryDiscretizeRegion[Ball[{0, 0, 0}, 16]]
Export["model.stl", myQuartic];
Export["ball.stl", ball];
stl = blenderworkflow[]

Differenced Object

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you very much for this answer and for making the effort to prepare the above script - it will be useful for me, even if I should stay using PyMesh instead of Blender (at the moment I need only the boolean intersection features from PyMesh and then next some mesh feature extraction, which can be done either in Mathematica or in Python. But I will look at the Blender documentation to get a feeling for the possibilities.) $\endgroup$
    – Jürgen Böhm
    Commented Nov 14, 2019 at 16:31
3
$\begingroup$

I am a bit late for the party but it can be done in Mathematica. Here is how:

Get the file and find the bounds:

file = URLDownload[
   "http://www.aviduratas.de/wc/model-scaled-x7.stl"];
myQuartic = Import[file, "BoundaryMeshRegion"];
RegionBounds[myQuartic]

(*{{-48.5411, 48.4823}, {-49.0028, 49.0095}, {-49.0087, 49.0065}}*)

Create a set of objects - use something more intelligent here:

s = 20;
balls = RegionUnion[
   Flatten[Table[
     Ball[{i, j, k}, 3], {i, -50, 50, s}, {j, -50, 50, s}, {k, -50, 
      50, s}]]];

To show case this, we will want to find the intersections of these:

Show[myQuartic, Region[balls]]

enter image description here

Load and use OpenCascadeLink:

Needs["OpenCascadeLink`"]
s1 = OpenCascadeShapeImport[file];
s2 = OpenCascadeShape[balls];
s3 = OpenCascadeShapeDifference[s1, s2];
bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[s3];

Visualize:

bmesh["Wireframe"[
  "MeshElementStyle" -> Directive[FaceForm[LightBlue], EdgeForm[]]]]

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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