# Smooth region plot from list of points

I have function that is costly to evaluate, fun[x,y], that takes two real arguments and returns True or False. I am trying to visualize in the x,y-plane the region where the function returns True.

I took a sample of points on a 2-dimensional grid on the x,y-plane, so I generated a list of the form:

{{x1,y1,True}, {x2,y2,True}, {x3,y3,False}, ...}

consisting of the coordinates of each point and the return value of the function.

Using ArrayPlot I can turn this into a visualization of the region where the function is true, but it is ugly because it has very straight edges. Is there a way to obtain a smoothed boundary?

Here is an example

fun[x_, y_] := (Boole[x^2 <= y])
data = Flatten[
Table[{x, y, fun[x, y]}, {x, 0, 1, 0.1}, {y, 0, 1, 0.1}], {1, 2}];
ListPlot[{Cases[data, {_, _, 1}], Cases[data, {_, _, 0}]}[[All,
All, {1, 2}]]]


Of course, fun[x,y] is a dummy example, in reality it is expensive to evaluate. It would be nice if could get a smooth region boundary.

• Can you post an example? – Henrik Schumacher Oct 22 '17 at 14:05
• @HenrikSchumacher I added an example. – becko Oct 22 '17 at 14:20
• Sigh As code? – Henrik Schumacher Oct 22 '17 at 14:21
• @HenrikSchumacher Done – becko Oct 22 '17 at 14:22
• Hm. Can we assume that the value of fun depends on the fulfillment of few rather smooth inequalities? – Henrik Schumacher Oct 22 '17 at 14:27

Then what about this: I interpolate with a BezierFunction and plot it's superlevel set:

dx = 0.1;
fun[x_, y_] := (Boole[x^2 <= y])
data = Table[{x, y, fun[x, y]}, {x, 0, 1, dx}, {y, 0, 1, dx}];
f = BezierFunction[data];
Show[
Quiet@RegionPlot[f[x, y][[3]] >= .5, {x, 0, 1}, {y, 0, 1}],
ListPlot[{
Cases[Flatten[data, 1], {_, _, 1}],
Cases[Flatten[data, 1], {_, _, 0}]}[[All, All, {1, 2}]]
]
]


BezierFunction only accepts arguments in the unit interval/square/cube. At least for rectangular domains, we can rescale and shift coordinates to obtain a slightly more flexible tool:

maxx = 3;
maxy = 2;
minx = miny = -1;
dx = 0.1;
dy = 0.2;
fun[x_, y_] := (Boole[x^2 <= y])
data = Table[{x, y, fun[x, y]}, {x, minx, maxx, dx}, {y, miny, maxy, dy}];
f = BezierFunction[data.DiagonalMatrix[{1./(maxx - minx), 1./(maxy - minx), 1.}]];
Show[
Quiet@RegionPlot[ f[(x - minx)/(maxx - minx) , (y - miny)/(maxy - miny) ][[3]] >= .5, {x, minx, maxx}, {y, miny, maxy}],
ListPlot[{Cases[Flatten[data, 1], {_, _, 1}], Cases[Flatten[data, 1], {_, _, 0}]}[[All, All, {1, 2}]]]
]


However, I do not know how to get rid of the requirement of even spacings of data points... Maybe one can also experiment with Interpolation to obtain suitable functions f. My experience with Interpolation was that the solutions were a bit too tight to data. Using smoothing splines instead may also help.

• That looks nice. I'll try it on my data and let you know. Thanks! – becko Oct 22 '17 at 15:13
• How do the $u,v$ parameters on the Bezier relate to my original $x,y$? In this example the $x,y$ go from 0 to 1, but in general this is not the case. – becko Oct 26 '17 at 7:49
• Thanks for the edit. Yes, the even spacing is too tight a requirement. – becko Oct 26 '17 at 12:13
• Okay, I posted a new approach... – Henrik Schumacher Oct 26 '17 at 12:30
• Yes! I was looking at it. Massive! Thanks a lot for your help! – becko Oct 26 '17 at 12:32

One way is to use the function as RegionFunction in ContourPlot. In the example below, I desined expensivefunction such that it keeps track about the evaluation points. You may use the options PlotPoints and MaxRecursion to control the amount of function evaluations and the accuracy of the sampling of the domain.

expensivefunction[x_?NumericQ, y_?NumericQ] := (
InternalStuffBag[pts, {x, y}];
x^2 + 0.1 Sin[3 Pi x] + y^2 <= 1.
);
f = {x, y} \[Function] expensivefunction[x, y];

pts = InternalBag[{}];
g = ContourPlot[1, {x, -1.2, 1.2}, {y, -1.2, 1.2},
RegionFunction -> f,
Mesh -> All,
BoundaryStyle -> Black,
PlotPoints -> 2,
MaxRecursion -> 2
];
Show[{g, Graphics[Point[InternalBagPart[pts, All]]]}]


With higher values for MaxRecursion you get better resolution, but of course, also mor points to evaluate your function on. The following is the same as above but with MaxRecursion -> 5:

For your example, one could try this

expensivefunction[x_?NumericQ, y_?NumericQ] := (InternalStuffBag[pts, {x, y}];x^2 - y <= 0.);
f = {x, y} \[Function] expensivefunction[x, y];
pts = InternalBag[{}];
g = ContourPlot[1., {x, 0, 1}, {y, 0, 1},
Mesh -> All,
PlotPoints -> {10, 5},
MaxRecursion -> 0,
RegionFunction -> f
];
Show[{g, Graphics[Point[InternalBagPart[pts, All]]]}]
InternalBagLength[pts]


Unfortunately, it seems to be impossible to tell ContourPlot exactly from which points to sample...

• This is a nice approach. But I would like to control the sampled points manually, if possible. So I generate the points myself, then the question is how I can plot a nice region out of it. – becko Oct 22 '17 at 14:22

And yet another approach. First, I try to detect the region of the interface between 0- and 1-valued points. Afterwards I use the minimimizer of the Dirichlet energy with respect 0- and 1-valued boundary conditions in this region; the final interface will be the 1/2-levelset.

First, the data points:

maxx = 3;
maxy = 2;
minx = miny = -1;
fun[x_, y_] := (Boole[x^2 <= y]);
nn = 600;
data0 = Transpose[{RandomReal[{minx, maxx}, nn], RandomReal[{miny, maxy}, nn]}];
data1 = Join[data0, Transpose[{fun @@@ data0}], 2];
listplot = Show[Graphics[], ListPlot[{Cases[data1, {_, _, 1}], Cases[data1, {_, _, 0}]}[[All, All, {1, 2}]],
PlotStyle -> {ColorData[97][3], ColorData[97][4]}]
]


Next, we coarsly detext the interface region:

Needs["NDSolveFEM"];
R = TriangleCreate[];
pts = data1[[All, 1 ;; 2]];
TriangleSetPoints[R, pts];
S = TriangleTriangulate[R, "a1.0"];
pts = TriangleGetPoints[S];
faces = TriangleGetElements[S];
pos = Flatten[Position[Equal @@@ Partition[data1[[Flatten[faces], 3]], 3], False, 1]];
nfaces = faces[[pos]];
plist = Union @@ nfaces;
npts = pts[[plist]];
nfaces = Partition[Lookup[lookuptable, Flatten[nfaces]], 3];
nullpos = Flatten[Position[data1[[plist, 3]], 0, 1]];
onepos = Flatten[Position[data1[[plist, 3]], 1, 1]];
R = MeshRegion[npts, Polygon[nfaces]];
With[{edges = DeveloperToPackedArray[MeshCells[R, 1][[All, 1]]]},
R0 = MeshRegion[MeshCoordinates[R], Line[Select[edges, SubsetQ[nullpos, #] &]]];
R1 = MeshRegion[MeshCoordinates[R], Line[Select[edges, SubsetQ[onepos, #] &]]];
];
GraphicsRow[{Show[listplot, R], Show[listplot, R0, R1]}, ImageSize -> Large]


Finally, we solve a boundary value problem. This can be done easier with the NDSolve facilities, but as I am not that much into these details, I do it a bit more verbatim.

S = DiscretizeRegion[R, MaxCellMeasure -> 0.0001];
n = Length[MeshCoordinates[S]];
nullpos = Flatten[Position[RegionMember[R0]@MeshCoordinates[S], True, 1]];
onepos = Flatten[Position[RegionMember[R1]@MeshCoordinates[S], True, 1]];
plist = Join[onepos, nullpos];
Module[{vd, sd, cdata, mdata, dpde, dbc, load, damping, bcdata, y, x, u},
Rdiscr = ToElementMesh[
"Coordinates" -> MeshCoordinates[S],
"MeshElements" -> {TriangleElement[MeshCells[S, 2][[All, 1]]]},
"MeshOrder" -> 1,
"NodeReordering" -> False];
vd = NDSolveVariableData[{"DependentVariables","Space"} -> {{u}, {x, y}}];
sd = NDSolveSolutionData[{"Space"} -> {Rdiscr}];
cdata = InitializePDECoefficients[vd, sd,
"DiffusionCoefficients" -> {{-IdentityMatrix[2]}},
"MassCoefficients" -> {{1}}, "LoadCoefficients" -> {{0}}
];
mdata = InitializePDEMethodData[vd, sd];
dpde = DiscretizePDE[cdata, mdata, sd];
dbc = DiscretizeBoundaryConditions[bcdata, mdata, sd]; {load,
stiffness, damping, mass} = dpde["All"];
];
costraints = SparseArray[Transpose[{Range[Length[plist]], plist}] -> 1., {Length[plist], n}, 0.];
L = ArrayFlatten[{{stiffness, Transpose[costraints]}, {costraints, 0.}}];
b = Join[ConstantArray[0., n], ConstantArray[1., Length[onepos]], ConstantArray[0., Length[nullpos]]];
x = LinearSolve[L, b][[1 ;; n]];
f = ElementMeshInterpolation[{Rdiscr}, x];
GraphicsGrid[{
{HighlightMesh[S, {0, nullpos}], HighlightMesh[S, {0, onepos}]},
{Plot3D[f[u, v], Element[{u, v}, S]],
Show[ContourPlot[f[u, v] == 1/2, Element[{u, v}, S]], listplot]
}}
, ImageSize -> Large]


This should also work for quite general domains.

• For me the last box does not work and gives an error about things not having the same shape. "12.0.0 for Linux" (Also on 11.3.0 for Linux) – Kvothe Sep 25 '19 at 16:10
• @Bvothe Hm. Apparently, Mathematica does not like if I evaluate  Needs["NDSolveFEM"]; in the Module`. Putting it outside seems to work, though. Thanks for pointing this out! – Henrik Schumacher Sep 25 '19 at 18:58