Skip to main content

computational geometry - 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:


screenshot convex hull



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



Answer



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


Comments

Popular posts from this blog

plotting - How to draw lines between specified dots on ListPlot?

I would like to create a plot where I have unconnected dots and some connected. So far, I have figured out how to draw the dots. My code is the following: ListPlot[{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4,13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full] I have thought using ListLinePlot command, but I don't know how to specify to the command to draw only selected lines between the dots. Do have any suggestions/hints on how to do that? Thank you. Answer One possibility would be to use Epilog with Line : ListPlot[ {{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4, 13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full, Epilog -> { Line[ ...

dynamic - How can I make a clickable ArrayPlot that returns input?

I would like to create a dynamic ArrayPlot so that the rectangles, when clicked, provide the input. Can I use ArrayPlot for this? Or is there something else I should have to use? Answer ArrayPlot is much more than just a simple array like Grid : it represents a ranged 2D dataset, and its visualization can be finetuned by options like DataReversed and DataRange . These features make it quite complicated to reproduce the same layout and order with Grid . Here I offer AnnotatedArrayPlot which comes in handy when your dataset is more than just a flat 2D array. The dynamic interface allows highlighting individual cells and possibly interacting with them. AnnotatedArrayPlot works the same way as ArrayPlot and accepts the same options plus Enabled , HighlightCoordinates , HighlightStyle and HighlightElementFunction . data = {{Missing["HasSomeMoreData"], GrayLevel[ 1], {RGBColor[0, 1, 1], RGBColor[0, 0, 1], GrayLevel[1]}, RGBColor[0, 1, 0]}, {GrayLevel[0], GrayLevel...

list manipulation - Selecting multiple columns from a matrix?

Sample data: data = { {{2013, 1, 1}, 24.13, 167.67, 231.82}, {{2013, 1, 2}, 32.15, 170.92, 225.99}, {{2013, 1, 3}, 35.43, 172.68, 221.67}, {{2013, 1, 4}, 36.73, 173.05, 218.32}, {{2013, 1, 5}, 58.19, 165.96, 197.05}, {{2013, 1, 6}, 69.99, 163.50, 187.52}, {{2013, 1, 7}, 71.37, 154.21, 175.58}, {{2013, 1, 8}, 72.51, 149.66, 163.25}}; I want a DateListPlot with three graphs, so for a matrix formed by columns 1 and 2, one for columns 1 and 3, and 1 for columns 1 and 4. At the moment I'm using this code: data2 = Transpose[{data[[All, 1]], data[[All, 2]]}]; data3 = Transpose[{data[[All, 1]], data[[All, 3]]}]; data4 = Transpose[{data[[All, 1]], data[[All, 4]]}]; DateListPlot[{data2, data3, data4}, Joined -> True, Filling -> {3 -> {1}}] but I have a hunch that this can be done more efficiently. I don't like the Transpose s in particular. Any ideas? edit (for extra credit) What if I need to multiply the second column by 2, which in my solution is simp...