Z-map of a MeshRegion

Question

TL;DR

How can I efficiently compute a matrix representing a z-map of a non-convex MeshRegion?

Longer explanation

I'm developing a collision detection algorithm for a 3-axis NC machine. My idea is to exploit the fact that we have 3 axis, to turn the problem from a detection of solids intersection to a sort of convolution between two matrices, representing z-maps (one for the workpiece, one for the tool), which should be faster.

My idea about how to create such z-map so far: importing from an .stl file the BoundaryMeshRegion of the workpiece, superimpose a 2D grid and then, for each cell of such grid, having a little Cuboid falling down. When such Cuboid touches the solid, I store the z at which the contact happened. A weaker voxelization, so to say.

The code below works, but it's awfully slow, as it eats 14 GBs of memory and maxes all the cores of an i7-8700 CPU for two minutes to run over a 400x300 mm solid.

Is there any possible improvement or an entirely different route that can provide what I'm looking for without requiring so much time and computational power? Ideally I should be able to go down to the 100-micron level for a 500x500 mm workpiece and the z-map (this is important) should always approximate in excess, if needed (as does my falling Cuboid algorithm).

Code

tick[{left_, right_}, step_] := Module[{min, max, res},
{min, max} = If[right < left, {-left, -right}, {left, right}];
res = Range[min, max, step];
res = If[res[[-1]] == max, Most[res], res];
(*output*)
res = If[right < left, -res, res]
];

rain[x_, y_, {zz_, zmin_}, res_, solid_] := Block[{z = zz},
While[
RegionDisjoint[solid, Cuboid[{x, y, z}, {x + res, y - res, z - res}]] && z >= zmin,
(*PrintTemporary[z];*)
z -= res;
];
(*output*)
z
];

wp = Import["ExampleData/spikey.stl", "BoundaryMeshRegion"]
dt = 1.;(*example collision detection tolerance*)

lim = RegionBounds[wp];

x = tick[lim[[1]], dt];
y = tick[Reverse@lim[[2]], dt];
z = tick[Reverse@lim[[3]], dt];

Monitor[data = Table[
rain[x[[i]], y[[j]], {z[[1]], z[[-2]]}, dt, wp], {j, Length[y]}, {i, Length[x]}],
{ProgressIndicator[j, {0, Length[y]}], j, ProgressIndicator[i, {0, Length[x]}], i}]

MatrixPlot[data, ColorFunction -> "TemperatureMap"]

@HenrikSchumacher Eventually I changed route, I'm going with conservative rasterization of triangles which seems easier to implement. Nevertheless, thank you for your help, as you can see I accepted your answer. — Marco, Jun 17 '18 at 15:29

score 3 · Accepted Answer · answered May 31 '18 at 23:55

R = Import["ExampleData/spikey.stl", "BoundaryMeshRegion"]
lim = RegionBounds[R];
dt = 0.01;
x = tick[lim[[1]], dt];
y = tick[Reverse@lim[[2]], dt];
z = tick[Reverse@lim[[3]], dt];
zmax = Max[z];

getZValue = Compile[{{a, _Integer, 1}},
  Block[{i = 1},
   While[i < Length[a] && Compile`GetElement[a, i] != 1, i++];
   If[i < Length[a], i+1,i]
   ],
  CompilationTarget -> "C",
  RuntimeAttributes -> {Listable},
  Parallelization -> True,
  RuntimeOptions -> "Speed"
  ]

Testing for intersection is a bit expensive. Thus, we use a test which is a bit cheaper by employing RegionDistance and effectivle testing against balls of radius 0.5 dt Sqrt[3.] that are circumscribed around the cubes.

f = RegionDistance[R];
pts = Developer`ToPackedArray@Flatten[Outer[List, x, y, z, 1], 2]; // 
  AbsoluteTiming // First
data = getZValue[
     ArrayReshape[UnitStep[0.5 dt Sqrt[3.] - f[pts]], 
      Length /@ {x, y, z}],
     dt,
     zmax
     ]; // AbsoluteTiming // First
MatrixPlot[data, ColorFunction -> "TemperatureMap"]

0.684264

8.39285

In general, this z-map should be a bit too high. You can refine it by performing an additional step of your rain method.

Z-map of a MeshRegion

1 Answers1