8
$\begingroup$

I would like to create a cubic or 3d grid mesh with a constant size for a surface mesh that represents a hollow torus saved in an STL file. mesh

I am import the vertex using:

pts = Import["Torus.stl", "VertexData"]
size = Length[pts]

 xmin = Min[pts[[All, 1]]];   ymin = Min[pts[[All, 2]]];   zmin =  Min[pts[[All, 3]]];
 xmax = Max[pts[[All, 1]]];    ymax = Min[pts[[All, 2]]];   zmax = Min[pts[[All, 3]]];


 (*positive points*)

   pts[[All,1]] =  pts[[All,1]]+ Abs[xmin]
   pts[[All,2]] =  pts[[All,2]]+ Abs[ymin]
   pts[[All,3]] =  pts[[All,3]]+ Abs[zmin]

 (*Create the points of the space, then I need to evaluate if these points are near to the vertex*)

    xspace = Table[i, {i, xmin , xmax, 0.4}]; 
    yspace = Table[i, {i, ymin , ymax, 0.4}];  
    zspace = Table[i, {i, zmin , zmax, 0.4}];


  (*Use a for loop to determine If the points of the background mesh are near to the vertex or not, and save the points near to the vertex points*)
  dataMesh = {};
  For[k = 1, k < size + 1, k++,
     For[j = 1, j < size + 1, j++,
       For[i = 1, i < size + 1, i++,
  tol = 0.4;
  Posx = xspace[[i]]; Posy = yspace[[j]]; Posz = zspace[[k]];
   .
   .
   .
  ]]]

In evaluating this loop, Mathematica takes several hours. Is there another way to generate this type of mesh from an STL File? . Since it is a curved surface, it is to be expected that the meshing is not perfect.

Improve this question
$\endgroup$
4
  • $\begingroup$ What do you mean by 'background mesh' and 'cubic mesh' ? Maybe an illustration of what you're trying to do would help? $\endgroup$
    – flinty
    Dec 25, 2020 at 16:07
  • $\begingroup$ I have uploaded the STL file, I mean a cubic grid. The idea is to generate a grid mesh, where the lenght of each cube is 0.4 or less. $\endgroup$
    – F.Mark
    Dec 25, 2020 at 16:29
  • $\begingroup$ Right, you are trying to voxelize a torus. $\endgroup$
    – flinty
    Dec 25, 2020 at 16:30
  • $\begingroup$ yes but, the idea is that each cube or grid element has a constant size. In this way is necessary to define a cut-off at the edge of the torus (tol in my program). $\endgroup$
    – F.Mark
    Dec 25, 2020 at 17:03

1 Answer 1

Reset to default
14
$\begingroup$

Update 2: workflow to create perfectly cubical voxels

In update 1, I discovered that although MaxCellMeasurewill allow you to control the resolution of the base mesh, ToElementMesh makes some internal choices to refine the mesh. Unfortunately, this refinement makes it virtually impossible to guarantee that the voxels are perfect cubes. Therefore, I created a workflow that builds the base mesh by hand to create a cubic voxel.

Helper functions

For simplicity, I wrapped the workflow into a function called voxelize as shown below:

(*Import required FEM package*)
Needs["NDSolve`FEM`"];
(*Tensor product mesh from:https://wolfram.com/xid/0rs5ccudm-eqv31q*)


pointsToMesh[data_] :=
  MeshRegion[Transpose[{data}], 
   Line@Table[{i, i + 1}, {i, Length[data] - 1}]];
Clear[voxelize]
voxelize[fil_, nelm_ : 10] := 
 Module[{stl, rb, extents, mcm, paddedextents, rp, crd,
   inc, em, vol, rmfv, reg, regmarkerfn, mean, regmarkers, mesh,
   subset, torusinc, gc, full, cut, unique, rules, torusem},
  (*Import required STL file*)
  stl = Import[fil];
  (*Center region*)
  stl = TransformedRegion[stl, 
    TranslationTransform[-Mean@Transpose@RegionBounds[stl]]];
  (*Get bounding box of the shape*)
  rb = RegionBounds[stl];
  extents = Flatten@Differences[rb, {0, 1}];
  (*Calculate maximum cell measure*)
  mcm = Min@Differences[rb, {0, 1}]/nelm;
  (*Pad the extents so we can make the voxels perfect cubes*)
  paddedextents = mcm Round[#/mcm] & /@ extents;
  rp = RegionProduct @@ ((pointsToMesh@
         Subdivide[#1, #2, #3]) & @@@ ({-#/2, #/2, Round[#/mcm]} & /@ 
        paddedextents));
  crd = MeshCoordinates[rp];
  (*grab hexa element incidents RegionProduct mesh*)
  inc = Delete[0] /@ MeshCells[rp, 3];
  (*Get bounding cuboid*)
   em = ToElementMesh["Coordinates" -> crd, 
    "MeshElements" -> {HexahedronElement[inc]}];
  (*Create volume mesh of torus*)
  vol = MeshRegion@ToElementMesh[stl];
  (*Get region member function of torus*)
  rmfv = RegionMember[vol];
  (*Association for Clearer Region Assignment*)
  reg = <|"main" -> 1, "incl" -> 2|>;
  regmarkerfn = If[rmfv[#], reg["incl"], reg["main"]] &;
  (*Get mean coordinate of each hexa for region marker assignment*)
  mean = Mean /@ GetElementCoordinates[em["Coordinates"], #] & /@ 
     ElementIncidents[em["MeshElements"]] // First;
  regmarkers = regmarkerfn /@ mean;
  (*Create marked region element mesh*)
       mesh = 
   ToElementMesh["Coordinates" -> em["Coordinates"], 
    "MeshElements" -> {HexahedronElement[inc, regmarkers]}];
  (*Find hexa subset that belongs to torus*)
  subset = 
   Flatten@Position [
     ElementMarkers[
      First@mesh["MeshElements"]], _?(# == reg["incl"] &), 1];
  torusinc = First[ ElementIncidents[mesh["MeshElements"]]][[subset]];
  (*ToElementMesh does not like extra coordinates*)
  (*Therefore,the element incidents need to be renumbered*)
  unique = Union@Flatten@torusinc;
       rules = AssociationThread[# -> Range@Length@#] &@unique;
  torusem = 
   ToElementMesh["Coordinates" -> crd[[unique]], 
    "MeshElements" -> {HexahedronElement[torusinc /. rules]}];
  (*GraphicsComplex does not care about the extra coordinates*)
  gc = GraphicsComplex[mesh["Coordinates"], {Hexahedron[torusinc]}];
  full = Graphics3D[gc, Boxed -> False];
  cut = Graphics3D[gc, Boxed -> False,
      PlotRange -> (mesh["Bounds"] /. {x_, y_} -> {0, y})];
  <|"stl" -> stl, "nelm" -> nelm, "mesh" -> mesh, "gc" -> gc, 
   "full" -> full, "cut" -> cut, "mean" -> mean,
   "crd" -> torusem["Coordinates"], "inc" -> torusinc, 
   "torusmesh" -> torusem|>
  ]
stlfile = "torus.stl";

Background mesh

With voxelize, it is easy to demonstrate that your input resolution is your output resolution. For example, if I want 4 elements across the minimum dimension, I simply execute:

vox = voxelize[stlfile, 4];
Show[vox["stl"], vox["mesh"]["Wireframe"]]

Background mesh with voxelize

As you can see, there are 4 elements across the vertical/minimum dimension.

Resolution study

It is easy to set up a resolution study like so:

Grid[Transpose@Table[With[{v = voxelize[stlfile, #]},
      {StringTemplate["Full extracted torus (nelm=``)"][v["nelm"]],
       v["full"], StringTemplate["Cutaway (nelm=``)"][v["nelm"]],
       v["cut"]}
      ] &[i], {i, 2, 10, 3}], Frame -> All]

Resolution study

Accessing the torus mesh info

In response to @F.Mark's comment, I added additional keys to the output of the voxelize module function to access the torus voxelized mesh. For example:

vox = voxelize[stlfile, 2];
vox["torusmesh"]["Wireframe"]
Short[vox["crd"], 5]
Short[vox["inc"], 5]

Torus data access

More complex objects

Here's an example to show that it can handle more complex STL objects.

vox = voxelize["http://exampledata.wolfram.com/gear.1", 10];
Show[{Graphics3D[{Opacity[1], vox["gc"]}, Boxed -> False]}]

Gear

Update 1: workflow to create volumetric region

The following uses ToElementMesh to create a bounding hexa mesh around the torus. Then it uses a region member function to identify the hexahedra that lie within the torus to create a separate region. These hexahedra can then be extracted by the region marker.

(*Import required FEM package*)
Needs["NDSolve`FEM`"];
(*Import required STL file*)
stl = Import["torus.stl"]
(*Get bounding box of the shape*)
rb = RegionBounds[stl];
(*Desired number of elements across minimum dimension*)
nelm = 10;
(*Calculate maximum cell measure*)
mcm = Min@Differences[rb, {0, 1}]/nelm;
(*Get bounding cuboid*)
em = ToElementMesh[Cuboid[Delete[0]@Transpose@rb], 
   MaxCellMeasure -> {"Length" -> mcm}, "MeshOrder" -> 1];
inc = ElementIncidents[em["MeshElements"]][[1]];
(*Create volume mesh of torus*)
vol = MeshRegion@ToElementMesh[stl];
(*Get region member function of torus*)
rmfv = RegionMember[vol];
(*Association for Clearer Region Assignment*)
reg = <|"main" -> 1, "incl" -> 2|>;
regmarkerfn = If[rmfv[#], reg["incl"], reg["main"]] &;
(*Get mean coordinate of each hexa for region marker assignment*)
mean = Mean /@ GetElementCoordinates[em["Coordinates"], #] & /@ 
    ElementIncidents[em["MeshElements"]] // First;
regmarkers = regmarkerfn /@ mean;
(*Create and view element mesh*)
Print["Converted Hexa Element Mesh Cutaway Drawing"]
mesh = ToElementMesh["Coordinates" -> em["Coordinates"], 
   "MeshElements" -> {HexahedronElement[inc, regmarkers]}];
mesh[
  "Wireframe"["MeshElement" -> "MeshElements", 
    "MeshElementStyle" -> (Directive[Opacity[0.5], FaceForm[#](*, 
              EdgeForm[]*)] &  /@ {Orange, Blue}),
    PlotRange -> (rb /. {x_, y_} -> {0, y})]]

Volume mesh of bounding region and torus

To extract the torus is a GraphicsComplex, you could do the following:

subset = Flatten@
   Position [
    ElementMarkers[
     First@mesh["MeshElements"]], _?(# == reg["incl"] &), 1];
torusinc = First[ ElementIncidents[mesh["MeshElements"]]][[subset]];
gc = GraphicsComplex[mesh["Coordinates"], {Hexahedron[torusinc]}];
Print["Full extracted torus"]
Graphics3D[gc, Boxed -> False]
Print["Cutaway to show that it is a volume mesh"]
Graphics3D[gc, Boxed -> False,
   PlotRange -> (rb /. {x_, y_} -> {0, y})]

Extracted torus

Previous answer

Note: I may delete this answer since, after closer inspection, I found out that it is only returning a boundary mesh and not a volume mesh.

You could use RegionImage and ImageMesh (note that Windows 10 gave me a virus warning on the rar file, so I made my own torus.stl):

stl = Import["torus.stl"]
ImageMesh@RegionImage[stl, RasterSize -> 100]

enter image description here

Improve this answer
$\endgroup$
9
  • $\begingroup$ Thanks so much!!, do you know if is possible to save the position (x,y,z) of the centers of each cube that contain the torus? For example, If I save the mesh using, Export["mesh.txt", mesh["Coordinates"], "Table"], I save the box and the torus, and I do not see the difference between them. $\endgroup$
    – F.Mark
    Dec 26, 2020 at 2:25
  • $\begingroup$ I have obtained the position of the centers, thank you very much. Could you explain the code a little more, how does the variable em and mcm work? I don't really understand how they work $\endgroup$
    – F.Mark
    Dec 26, 2020 at 14:34
  • $\begingroup$ @F.Mark for posterity, to select the centroids to the elements, you need only issue the command Select[mean, rmfv[#] &]. The element mesh, em, is a simple hex mesh of the bounding box of the structure of interest. The structure of interest will be contained in this bounding box. To control the element size, we use the variable, mcm , or maximum cell measure. The assignment calculates the extents of the bounding box and uses the minimum extent divided by the number of elements desired across at minimum extent to determine the maximum cell measure. $\endgroup$
    – Tim Laska
    Dec 26, 2020 at 17:21
  • $\begingroup$ Thank you very much, finally, do you know if you can define the size of the grid, for example, 2x2x2 or 4x4x4? If I modify the nelm parameter I can decrease the size $\endgroup$
    – F.Mark
    Dec 26, 2020 at 18:26
  • $\begingroup$ @F.Mark, the way I designed the workflow was to define the number of elements across the minor diameter. So, nelm =10 means there were 10 elements across the minor diameter. To double the resolution, simply make nelm=20. For efficiency, I snapped the bounding box to each dimension. The problem with three-dimensional modeling is that the mesh sizes grow quickly. It would be quite inefficient to use an NxNxN for a high aspect ratio model. $\endgroup$
    – Tim Laska
    Dec 26, 2020 at 21:06

Your Answer

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

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