# How to Interpolate a 3D MeshRegion?

If you were given some discrete MeshRegion called r (you don't know R):

R = ImplicitRegion[
x^6 - 5 x^4 y z + 3 x^4 y^2 + 10 x^2 y^3 z + 3 x^2 y^4 - y^5 z +
y^6 + z^6 <= 1, {x, y, z}];
r = DiscretizeGraphics @ RegionPlot3D[R,
PlotRange -> {{-1.5, 1.5}, {-1.5, 1.5}, {-1.6, 1.6}},
PlotPoints -> 50, Boxed -> False]


What would be a good way to interpolate R, the underlying continuous equational or functional representation? That is, how can you convert a given MeshRegion to an ImplicitRegion or ParametricRegion?

• Do you know anything about R? Anything about its form? Or just have r? – MikeY Dec 28 '18 at 15:38

A very simple method is by interpolating the signed distance function (SignedRegionDistance) of a BoundaryMeshRegion. Therefore, I use BoundaryDiscretizeGraphics instead of DiscretizeGraphics to produce a BoundaryMeshRegion called S. Afterwards, I sample SignedRegionDistance on a regular grid surrounding S and apply Interpolation. This uses bi-cubic splines per default. Of course any other interpolation/approximation technique can be used, for example, trigonometric polynomials, radial basis functions etc.

S = BoundaryDiscretizeGraphics@
RegionPlot3D[R,
PlotRange -> {{-1.5, 1.5}, {-1.5, 1.5}, {-1.6, 1.6}},
PlotPoints -> 50, Boxed -> False];

box = List @@ BoundingRegion[S];
center = Mean[box];
scale = 1.5;
B = Cuboid @@ ({center, center} + scale (box - {center, center}));
n = 20;
nodes = Tuples[Subdivide[#[[1]], #[[2]], n] & /@ Transpose[List @@ B]];
df = SignedRegionDistance[S];
f = Interpolation[Transpose[{nodes, df[nodes]}]];

Show[
S,
SliceContourPlot3D[f[x, y, z], "CenterPlanes", {x, y, z} \[Element] B]
]


And here is a contourplot of the interpolated signed distance function:

ContourPlot3D[
f[x, y, z] == 0,
{x, B[[1, 1]], B[[2, 1]]}, {y, B[[1, 2]], B[[2, 2]]}, {z, B[[1, 3]], B[[2, 3]]}
]


A rediscretized region looks like this:

R1 = ImplicitRegion[
f[x, y, z] <= 0,
{{x, B[[1, 1]], B[[2, 1]]}, {y, B[[1, 2]], B[[2, 2]]}, {z, B[[1, 3]], B[[2, 3]]}}
];
S1 = BoundaryDiscretizeRegion[R1, MaxCellMeasure -> {1 -> 0.1}];
GraphicsRow[{Show[S], Show[S1]}]


(left: original region; right: rediscretized region)

• Can you compute a closed form for the normal to the surface at a point? – M.R. Dec 28 '18 at 10:28
• Well, the normal at a point of a levelset is just the normalized gradient (times $1$ or $-1$, depending on inwards/outwards pointing). I guess normal = {x, y, z} \[Function] Evaluate[ D[f[x, y, z], {{x, y, z}, 1}]/ Sqrt[D[f[x, y, z], {{x, y, z}, 1}].D[ f[x, y, z], {{x, y, z}, 1}]] ]; would do. – Henrik Schumacher Dec 28 '18 at 10:35
• @HenrikSchumacher Only knowing r I Tried RegionPlot3D[RegionDistance[r, {x, y, z}] == 0 , PlotRange -> {{-1.5, 1.5}, {-1.5, 1.5}, {-1.6, 1.6}}] to approximate R directly, but Mathematica doesn't evaluate output. What's wrong with this idea? Thanks! – Ulrich Neumann Dec 28 '18 at 20:50
• @Ulrich There are several problems with this. I) RegionPlot3D is better at plotting volumetric regions. Thus it should better applied to inequalities. E.g., df = RegionDistance[r]; RegionPlot3D[ df[{x, y, z}] <= .1, {x, -1.5, 1.5}, {y, -1.5, 1.5}, {z, -1.6, 1.6}, PlotPoints -> 50] works. II) RegionPlot3D[ df[{x, y, z}] <= 0., {x, -1.5, 1.5}, {y, -1.5, 1.5}, {z, -1.6, 1.6}, PlotPoints -> 50] does not work because the volume is actually a surface and hence has 3-dim Hausdorff measure 0. Hence it is hard to detect numerically. – Henrik Schumacher Dec 28 '18 at 22:43
• The latter tends to happen with all minimizing sets of functions. Another issue is that RegionDistance[r] is not differentiable at the levelset in question. With a signed-distance function, Newton's method is available because the signed distance function (of a smooth embedded hyper-surface) i) is differentiable and ii) has surjective differential (both in a neighborhood of the surface). – Henrik Schumacher Dec 28 '18 at 22:48

In my attempt I try to realize a very simple signed distance function which only needs information about the meshregion r(assumption: triangle mesh in space)

points = MeshCoordinates[ r]  ;
elements = MeshCells[r, 2]  ; (* all triangle elements *)
triangles = elements /. Polygon -> List // Flatten[#, 1] &;(* node indices*)


Furthermore I use/presume a function, which detects the nearest element of a given point.

mnc = RegionMeshMeshNearestCellIndex[r ]
(*mnc[x,y,z] evaluates the index of nearest element*)


The following interpolation idea doesn't need an embedding grid, uses only the given nodes and calculates the signed distance in normal element direction. Because it only considers the nearest element it's some kind of one dimensional interpolation

fUN[x_?NumericQ, y_?NumericQ, z_?NumericQ] := Block[{n\[CapitalDelta], no, p1, p2,p3},
n\[CapitalDelta] = mnc[{x, y, z}][[2]];  (*nearest triangle element *)
{p1, p2, p3} = punkte[[dreiecke[[n\[CapitalDelta] ]]]];(* triangle points*)
no = Cross[p2 - p1, p3 - p1 ];(* normal vector triangle *)
no.({x, y, z} - p1)(* normal distance*)


]

If fUN returns zero the point x,y,z is element of the triangle plane.

RUN = ImplicitRegion[fUN[x, y, z] == 0, {{x, -1.5, 1.5}, {y, -1.5, 1.5}, {z, -1.6,1.6}}] (* interpolated region R*)

DiscretizeRegion[RUN]


The program seems to work quite well (no idea where the sliver comes from...)

Function fUN might be improved checking x,y,z inside triangle and using compile!