# Smooth convex hull of a large data set of 3D points

I recently needed to deal with a large data set of 600,000 points in three dimensions. My task was to find a convex hull for this data. I have Mathematica 10, so I could use the function ConvexHullMesh; I obtained this:

I was wondering if there is some way to find a smooth convex hull (maybe an ellipsoid) for my data.

• What exactly do you mean by 'smooth convex hull', or, how is an ellipsoid a convex hull for your data? Could you explain that a bit. Aug 21, 2014 at 7:39
• I second user21's question; there is no standard notion of a "smooth convex hull". If you really really want an ellipsoid, you could try finding the minimum-volume ellipsoid enclosing your set of points.
– user484
Aug 21, 2014 at 9:11
• Yes, I mean, smooth convex hull as differential function. So, the other comment is probably that I need. I'll read the paper, and If I success, I write hear my trie. Aug 21, 2014 at 15:49
• There are routines to calculate MVE ellipsoids in MATLAB. You could just translate one to Mathematica.
– user24615
Jan 9, 2015 at 18:17

Minimum Volume Ellipsoid

Translated from here, this uses the Khachiyan algorithm, and should work for any dimension.

    MinVolEllipse[P_, tolerance_] :=
Module[{d, n, Q, count, err, u, X, M, maximum, j, stepSize, newu, U, A, c},

{n, d} = Dimensions[P];
Q = Append[1] /@ P;

count = 1;
err = 1;
u = ConstantArray[1./n, n];

While[err > tolerance,
X = Q\[Transpose].DiagonalMatrix[u].Q;
M = Diagonal[Q.Inverse[X].Q\[Transpose]];
maximum = Max[M];
j = Position[M, maximum][[1, 1]];
stepSize = (maximum - d - 1)/((d + 1) (maximum - 1));
newu = (1 - stepSize) u;
newu[[j]] += stepSize;
count += 1;
err = Norm[newu - u];
u = newu;
];

U = DiagonalMatrix[u];

A = (1/d) Inverse[P\[Transpose].U.P - Outer[Times, u.P, u.P]];
c = u.P;

{c, A}
]


Usage:

    pts = RandomVariate[
MultinormalDistribution[RandomReal[{-1, 1}, {2}],
With[{m = RandomReal[{0, 1}, {2, 2}]}, m.m\[Transpose]]], 500];

P = MeshCoordinates[ConvexHullMesh[pts]];
tolerance = 0.0001;

{c, A} = MinVolEllipse[P, tolerance];

X = {x, y};
Show[
ConvexHullMesh[pts],
Graphics[{
Point[pts],
{Red, Point[P]}
}],
ContourPlot[(X - c).A.(X - c) == 1, {x, -4, 4}, {y, -4, 4}]
]


In 3D:

    pts = RandomVariate[
MultinormalDistribution[RandomReal[{-1, 1}, {3}],
With[{m = RandomReal[{0, 1}, {3, 3}]}, m.m\[Transpose]]], 100];

P = MeshCoordinates[ConvexHullMesh[pts]];
tolerance = 0.0001;

{c, A} = MinVolEllipse[P, tolerance];

X = {x, y, z};
Show[
ConvexHullMesh[pts, MeshCellStyle -> {{2, All} -> Opacity[0.5, Green]}],
Graphics3D[{
Point[pts],
{Red, Point[P]}
}],
ContourPlot3D[(X - c).A.(X - c) == 1, {x, -3, 3}, {y, -3, 3}, {z, -3, 3},
ContourStyle -> Directive[Red, Opacity[0.2]]
]
]


Here is a faster implementation of the Khachiyan minimum-volume ellipsoid algorithm, using updating formulae presented in this paper:

mve[pts_?MatrixQ, tol_: 1.*^-8] :=
Module[{prec = Precision[pts], bet, bm, c, d, del, dp, h, j, kap, n,
qj, qm, sc, sig, u, zero},
zero = SetPrecision[0, prec]; {n, d} = Dimensions[pts]; dp = d + 1;
qm = PadRight[pts, {n, dp}, N[1, prec]];
u = N[ConstantArray[1/n, n], prec];
bm = Transpose[qm].DiagonalMatrix[u].qm;
kap = Diagonal[qm.LinearSolve[bm, Transpose[qm]]];
c = FixedPoint[(j = Ordering[kap, -1]; qj = Extract[qm, j];
sc = 1/Extract[kap, j]; sig = (1 - sc dp)/(1 - sc);
bet = sig/dp; h = 1 - bet; del = dp (1 - sc)/d;
kap = del (kap - sc sig (qm.LinearSolve[bm, qj])^2);
bm = h bm + bet Outer[Times, qj, qj];
h # + SparseArray[{j -> bet}, n, zero]) &, u,
SameTest -> (SquaredEuclideanDistance[##] <= tol &)].pts;
Ellipsoid[c, d (Drop[bm, -1, -1] - Outer[Times, c, c])]]


The routine returns an Ellipsoid[] object, suitable for display or any other further manipulations (e.g. one can apply an eigendecomposition to the matrix in the second argument of the Ellipsoid[] object to get the axes).

BlockRandom[SeedRandom[42, Method -> "MersenneTwister"]; (* for reproducibility *)
pts = RandomVariate[MultinormalDistribution[RandomReal[{-1, 1}, {2}],
With[{m = RandomReal[1, {2, 2}]}, m.Transpose[m]]], 1024];]

Graphics[{Red, mve[pts], Blue, MeshPrimitives[ConvexHullMesh[pts], 2], Yellow, Point[pts]}]


BlockRandom[SeedRandom[42, Method -> "Rule30CA"];
pts = RandomVariate[MultinormalDistribution[RandomReal[{-1, 1}, {3}],
With[{m = RandomReal[1, {3, 3}]}, Transpose[m].m]], 1024];]

Graphics3D[{Opacity[1/2, Red], mve[pts], Blue, MeshPrimitives[ConvexHullMesh[pts], 2]}]