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

Question

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]?

Answer 1

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]
 ]

Answer 2

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["NDSolve`FEM`"]

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[0];
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.