# Finding the equation for the upper frontier of the convex hull of a 2 dimensional set of points

Suppose you have

data2D = {{4.4, 14}, {6.7, 15.25}, {6.9, 12.8}, {2.1, 11.1}, {9.5,
14.9}, {13.2, 11.9}, {10.3, 12.3}, {6.8, 9.5}, {3.3, 7.7}, {0.6,
5.1}, {5.3, 2.4}, {8.45, 4.7}, {11.5, 9.6}, {13.8, 7.3}, {12.9,
3.1}, {11, 1.1}};

g = ConvexHullMesh[data2D];

extremepoints = MeshCoordinates[g];


Is there a built in Mathematica way of generating a function f[x_] where the domain is [.6,13.8] (the min and max of the x-coordinate in data2D) and the function delivers the upper frontier value of ConvexHullMesh[data2D]?

• You might think MaxValue[y && x == #, {x, y} \[Element] g] & is a natural approach, but it does not seem to be implemented. (Not for NMaxValue either.) – Michael E2 Oct 28 '14 at 16:48

You can build one out of region functions, at least for mesh regions. One could always discretize other regions, but the answer will be approximate and probably not highly accurate.

More efficient

This one is based on two assumptions about ConvexHullMesh:

1. The linear mesh elements are order counterclockwise around the region.
2. The first element starts at the point with the minimum x coordinate.

Then we can select the top boundary by selecting the line elements whose x coordinates are in the reverse order. The elements themselves will also be in reverse order, but that does not matter to Interpolation. It may help to know that a line element has the form Line[{i, j}], where i, j are integer indices of the corresponding points in the list of MeshCoordinates. Consecutive elements have the form Line[{i, j}], Line[{j, k}], with the index j repeated. (This leads to the need for DeleteDuplicate below.)

This method has the advantage over my original solution of producing a single interpolating function, which is both more efficient and more convenient.

Clear[topy];
topy[reg_?BoundaryMeshRegionQ] := With[{coords = MeshCoordinates[reg]},
Composition[
Interpolation[#, InterpolationOrder -> 1,          (* 6. interpolate *)
"ExtrapolationHandler" -> {Indeterminate &, "WarningMessage" -> False}] &,
coords[[#]] &,                                     (* 5. extract coordinates *)
DeleteDuplicates,                                  (* 4. delete repeated point indices *)
Flatten,                                           (* 3. flatten list of indices *)
# /. Line -> List &,                               (* 2. change Line to List *)
Select[! OrderedQ[coords[[First@#]][[All, 1]]] &]  (* 1. select mesh elements *)
][MeshCells[reg, 1]]
]


Usage:

f = topy[g]
(* InterpolatingFunction[{{0.6, 13.8}}, <>] *)


The interpolating function f can be used like any function, such as the one below.

Original solution

The basic idea is to select boundary elements that lie over a given value of x. The take the maximum of the y coordinates of the lines.

(* For 2D mesh regions only *)
Clear[topy];
topy[reg_?MeshRegionQ, x0_?NumericQ] := topy[BoundaryMesh[reg], x0];
topy[reg_?BoundaryMeshRegionQ, x0_?NumericQ] :=
Module[{coords, lines, x1, x2, y1, y2, yvals},
{{x1, x2}, {y1, y2}} = RegionBounds[reg];
(coords = MeshCoordinates[reg];
lines = MeshCells[reg, 1];
lines = (* select lines over x0 *)
Cases[lines /. i_Integer :> coords[[i]],
Line[{{xa_, _}, {xb_, _}}] /; xa <= x0 <= xb || xb <= x0 <= xa];
yvals = lines /. Line[p_] :> Interpolation[p, InterpolationOrder -> 1][x0];
Max[yvals]) /; x1 < x0 < x2]

topy[g, 1.4]
(* 8.3 *)


Change Max to Min for the lower boundary.

The top function:

f = topy[g, #] &;

Show[
RegionPlot[g],
Plot[f[x], {x, 0, 16}, PlotStyle -> Red]
] • WRI is working on filling out region functionality. One might find in a month or so, that optimization functions can used to do this more or less directly. – Michael E2 Oct 28 '14 at 17:35
• This is beautiful. Thank you. – cphelan Oct 28 '14 at 20:32
• @cphelan Thanks & you're welcome. I've updated it with a better function. It can be extended to convex mesh regions like my first one, with topy[reg_?MeshRegionQ] := topy[BoundaryMesh[reg]]. (If you don't have V10, replace the Select line with Select[#, ! OrderedQ[coords[[First@#]][[All, 1]]] &] &. – Michael E2 Oct 29 '14 at 20:27
• This is working great and very useful. I didn't think I would need a 3D version, but am finding that would be very useful. The code above produces answers, even with 3D, but not the correct answers. Is there a small fix? – cphelan Oct 31 '14 at 19:36
• @cphelan A similar approach with new-in-V10 utilities is easier to extend to any dimension > 1 that is implemented in Mathematica's mesh function capabilities (i.e. 2 or 3 currently). See my new answer. – Michael E2 Nov 1 '14 at 17:54

Here's a pretty neat way, based on the new FEM utilities for dealing with meshes, that works automatically for dimension 2 and 3, positively oriented, convex, simplex-element meshes. (More precisely, instead of being convex, the intersection of the region only with an arbitrary vertical line needs to consist of at most one connected component).

Initialization:

Needs["NDSolveFEM"]

data2D = {{4.4, 14}, {6.7, 15.25}, {6.9, 12.8}, {2.1, 11.1}, {9.5,
14.9}, {13.2, 11.9}, {10.3, 12.3}, {6.8, 9.5}, {3.3, 7.7}, {0.6,
5.1}, {5.3, 2.4}, {8.45, 4.7}, {11.5, 9.6}, {13.8, 7.3}, {12.9, 3.1}, {11, 1.1}};
g = ConvexHullMesh[data2D];

SeedRandom;
chm = ConvexHullMesh[RandomReal[{0, 5}, {20, 3}]];


The function top produces an InterpolatingFunction. The top mesh elements are selected and projected onto the domain. The height coordinates (y or z) are used as the values for the interpolation over the projected mesh coordinates, basecoords. The top elements are determined by their orientation in the domain. ConvexHullMesh produces (or seems to) a positively oriented boundary in the sense of integration theorems of Green/Stokes/Gauss. Here orientation means something slightly different. It yields the vertical component of an outward normal to the element, which will be positive, if the element is on top.

ClearAll[top];
top[reg_?BoundaryMeshRegionQ] /; reg["SimplexMeshQ"] :=
Module[{
dim, orientation,
allsimplices, topsimplices, topindices, basecoords,
values, element, mesh},
dim = reg["Dimension"];
orientation = dim /. {
2 -> (-First[Differences[#[[All, 1]]]] &),
3 -> (Last[Cross @@ Differences[#]] &)
};
allsimplices = MeshCells[reg, dim - 1] /. h_[list : {__Integer}] :> list;
topsimplices = Pick[
allsimplices,
Positive[orientation /@ (reg["Coordinates"][[#]] &) /@ allsimplices]];
topindices = AssociationThread[# -> Range@Length@#] &[Union @@ topsimplices];
With[{topcoords = reg["Coordinates"] ~Part~ Keys[topindices]},
basecoords = topcoords[[All, 1 ;; -2]];
values     = topcoords[[All, dim]]
];
element = dim /. {2 -> LineElement, 3 -> TriangleElement};
mesh = ToElementMesh[
"Coordinates" -> basecoords,
"MeshElements" -> {element[topsimplices /. topindices]}];
ElementMeshInterpolation[{mesh}, values]
]


OP's example:

meshifn = top[g]

Show[
g,
Plot[meshifn[x], Evaluate@Flatten[{x, meshifn["Domain"]}], PlotStyle -> Red],
Frame -> True
] Random convex hull:

meshifn = top[chm]

plot = Plot3D[meshifn[x, y], {x, y} ∈ meshifn["ElementMesh"], PlotPoints -> 100]

Show[
RegionPlot3D[chm, PlotStyle -> LightBlue],
plot,
Axes -> True, AxesLabel -> {x, y, z}
] The funny color of the top surface in the plot on the right is due to the surfaces approximately coinciding.

It wouldn't be too hard to extend this to other sides.

• thank you for for a very instructive look into, esp ElementMeshIntetpolation +1, of course – ubpdqn Nov 2 '14 at 4:01