12
$\begingroup$

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:

screenshot convex hull

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

$\endgroup$
4
  • 4
    $\begingroup$ 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. $\endgroup$
    – user21
    Aug 21, 2014 at 7:39
  • 4
    $\begingroup$ 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. $\endgroup$
    – user484
    Aug 21, 2014 at 9:11
  • 1
    $\begingroup$ 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. $\endgroup$ Aug 21, 2014 at 15:49
  • 1
    $\begingroup$ There are routines to calculate MVE ellipsoids in MATLAB. You could just translate one to Mathematica. $\endgroup$
    – user24615
    Jan 9, 2015 at 18:17

2 Answers 2

12
$\begingroup$

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

enter image description here

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

enter image description here

$\endgroup$
10
$\begingroup$

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

minimum-volume ellipse and 2D convex hull

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

minimum-volume ellipsoid and 3D convex hull

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.