How to Interpolate a 3D MeshRegion?

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)

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 = Region`Mesh`MeshNearestCellIndex[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!