How to generate higher dimensional convex hull?
New implementation (EDIT:2017-May-08)
I just wanted to share a slightly different implementation (based on the old one, see below). The following function (almost the same as TestCandidate from the old implementation) tests if a candidate point cand can be described by a convex combination of the given corner points corners
ConvexDepenentQ[corners_, cand_] := Module[
{w, ws},
w = Array[ws, Length@corners];
1 == Length@
FindInstance[
w.corners == cand && Total[w] == 1 &&
And @@ Table[w[[i]] >= 0, {i, Length@w}], w]
];
The following function reduces a given data set data to its convex corners, which describe the convex hull. This is done by picking one data point out and testing if it is convex dependent of the rest.
ConvexReduce[data_] := Module[
{corners, ncorners, test},
corners = data;
Do[
ncorners = Delete[corners, Position[corners, data[[i]]]];
test = ConvexDepenentQ[ncorners, data[[i]]];
If[test, corners = ncorners];
, {i, Length@data}
];
corners
];
Naturally, you could get the convex corners just with ConvexReduce but we can speed this up a little by starting with the data points on the bounding rectangle (as in the old implementation), going through the data and only append points which are not convex dependent. At the end we can convex reduce the smaller set.
ConvexHullCorners[data_] := Module[
{corners, rd},
corners = {};
Do[
corners =
Join[corners,
Select[data,
Min[data[[;; , i]]] == #[[i]] ||
Max[data[[;; , i]]] == #[[i]] &]];
, {i, Length@data[[1]]}
];
corners = DeleteDuplicates@corners;
rd = Delete[data, First@Position[data, #] & /@ corners];
Do[
If[ConvexDepenentQ[corners, rd[[i]]], , AppendTo[corners, rd[[i]]]]
, {i, Length@rd}
];
ConvexReduce@DeleteDuplicates@corners
];
At least compared to the algorithms implemented by Loren Petrich, see http://lpetrich.org/Science/ and Jason's answer in this question, my new implementation seems to perform ok. But please, keep in mind that I obly look for the convex corners of the data (this is the easy part), while qhull also delivers the highdimensional facettes of the data (which is the veeeery hard and time expensive part, which I at least dont need). For increasing number of points, this seems to be ok.
d = 6;
np = 40;
ts = {};
Do[
data = RandomReal[{-10, 10}, {i, d}];
AppendTo[ts, {
AbsoluteTiming[ConvexHullCorners@data;][[1]]
, AbsoluteTiming[CHNChanShattering@data;][[1]]
, AbsoluteTiming[CHNGiftWrapping@data;][[1]]
, AbsoluteTiming[CHNIncremental@data;][[1]]
, AbsoluteTiming[CHNQuickHull@data;][[1]]
}];
, {i, 10, np, 10}
]
ListPlot[Transpose@ts, Joined -> True, PlotLegends -> {
"ConvexHullCorners"
, "CHNChan"
, "CHNGift"
, "CHNInc"
, "CHNQuickH"
}]
But specially for increasing dimension, the bounding rectangle approach seems to be a good start.
d = 7;
np = 30;
ts = {};
Do[
data = RandomReal[{-10, 10}, {np, i}];
AppendTo[ts, {
AbsoluteTiming[ConvexHullCorners@data;][[1]]
, AbsoluteTiming[CHNChanShattering@data;][[1]]
, AbsoluteTiming[CHNGiftWrapping@data;][[1]]
, AbsoluteTiming[CHNIncremental@data;][[1]]
, AbsoluteTiming[CHNQuickHull@data;][[1]]
}];
, {i, 2, d}
]
ListPlot[Transpose@ts, Joined -> True, PlotLegends -> {
"ConvexHullCorners"
, "CHNChan"
, "CHNGift"
, "CHNInc"
, "CHNQuickH"
}]
If you want to check, that the corners delivered by the algorithms are identical, use
Sort@ConvexHullCorners@data == Sort@data[[CHNQuickHull[data][[1]]]]
True
If you want to check, that all data points are convex dependent of the corners, use
And @@ (ConvexDepenentQ[ConvexHullCorners@data, #] & /@ data)
True
I am running these computations in Windows 10, Mathematica 11.0.1. Sadly, I was never able to connect to Matlab, so I dont know, if my implementation is slower (probably). Nevertheless, this might help somebody.
Old implementation
Thanks @physicien for the comment on Quickull and Qhull, see his comment on the question. Sadly, I did not understand the paper and so I wanted to give this a shot in Mathematica with some very rude implementation (sorry, I am extremly bad at programming). There is a lot of room for improvement (my solution is extremly slow in high dimensions and for large data), but may be this help also others like me, at least for a start.
For a candidate vector $v \in \mathbb{R}^n$ (cand) I will first test if the candidate is represtable by a convex combination of a given set of coners $\{c_1,c_2,...,c_m\}, c_i \in \mathbb{R}^n$ (corners)
\begin{equation} \sum_{i=1}^m w_i c_i = v \ , \qquad \sum_{i=1}^m w_i = 1 \ , \qquad w_i \geq 0 \ \forall \ i \end{equation}
with the following function.
TestCandidate[corners_, cand_] := Block[
{ws, w, eqs, fi},
ws = Array[w, {Length[corners]}];
eqs = Sum[ws[[i]]*corners[[i]], {i, Length@corners}] - cand;
eqs = Flatten@Join[
Table[eqs[[i]] == 0, {i, Length@eqs}]
, {Total[ws] == 1}
, Table[ws[[i]] >= 0, {i, Length@ws}]
];
eqs = And @@ eqs;
fi = FindInstance[eqs, ws];
Length@fi != 0
];
This functions returns True if the candidate is representable, since if so, an instance solving this can be found, False otherwise, i.e., the candidate is not represtable as a convex combination.
Now I will search the set of corners, which is able to represent any point of a given data set (data) by convex combinations as follows:
Pick initial corner points from data belonging to a $n$-dimensional estimate for a bounding box
Delete the picked points from 1. from the data
Check if any points picked in 1. are representable by a smaller set (this happens sometimes if all points are actually in a $(n-1)$-dimensional plane.
Run through the rest of the data points and add a candidate if it is NOT representable by the updated corners. If a new corner is added, check all subsets of corners in case a new corner allows for the representation of old corners and reduce momentary set of corners.
Code
ConvexHullCorners[data_] := Block[
{datatemp, dim, center, low, pos, upp, corners, cand, cornerstemp},
(***********************)
(*1*)
(*Generate initial bounding rectangle estimate points*)
dim = Length[data[[1]]];
center = Mean[data];
corners = {};
Do[
low = Max[center[[i]] - data[[;; , i]]];
pos = Position[data[[;; , i]], center[[i]] - low][[1, 1]];
AppendTo[corners, data[[pos]]];
upp = Max[data[[;; , i]] - center[[i]]];
pos = Position[data[[;; , i]], center[[i]] + upp][[1, 1]];
AppendTo[corners, data[[pos]]];
, {i, dim}
];
corners = DeleteDuplicates@corners;
(***********************)
(*2*)
(*Delete points from bounding box estimate from data and \
duplicates*)
pos = Flatten[
Table[Position[data, corners[[i]]], {i, Length@corners}], 1];
datatemp = DeleteDuplicates@Delete[data, pos];
(***********************)
(*3*)
(*Check for representable points,
and delete them from corners if so*)
cornerstemp = corners;
Do[
cand = corners[[j]];
pos = Position[cornerstemp, cand][[1]];
cornerstemp = Delete[cornerstemp, pos];
If[TestCandidate[cornerstemp, cand] == False,
PrependTo[cornerstemp, cand]];
, {j, Length[corners]}
];
corners = cornerstemp;
(***********************)
(*4*)
(*Find corners of convex hull*)
Do[
(*Test candidate from rest of data*)
cand = datatemp[[i]];
If[
(*If not representable, add to corners*)
TestCandidate[corners, cand] == False,
AppendTo[corners, cand];
(*Check if just added corner allows for representation of some \
former old corners, and delete them from corners if so*)
cornerstemp = corners;
Do[
cand = corners[[j]];
pos = Position[cornerstemp, cand][[1]];
cornerstemp = Delete[cornerstemp, pos];
If[TestCandidate[cornerstemp, cand] == False,
PrependTo[cornerstemp, cand]];
, {j, Length[corners] - 1}
];
corners = cornerstemp;
];
, {i, Length@datatemp}
];
(***********************)
(*OUT*)
(*Return corners of convex hull*)
corners
];
Some examples:
1) 2D example
Generate some random points in 2D (red points), compute the normal convex hull with Mathematica (blue region), compute the corners with the program given above (green points) and visualize results.
(*Dimension and number of points*)
d = 2;
np = 10^1*2;
(*Generate data*)
data = RandomReal[{-10, 10}, {np, d}];
ch = ConvexHullMesh@data;
datac = ConvexHullCorners@data;
pic = Show[
ch,
Graphics[{PointSize -> 0.05, Green, Point[datac],
Directive[Red, Opacity[0.4]], PointSize -> Large, Point[data]},
ImageSize -> Small]
]
2) 6D example:
Analogously, generate data and extract corners of convex hull with program. Check if all points of the data are really representable by the given corners of the program and check a false point.
(*Dimension and number of points*)
d = 6;
np = 10^1*3;
(*Generate data*)
int = {5, 10};
dis = 7;
data = RandomInteger[int, {np, d}];
(*Extract corner of convex hull*)
datac = ConvexHullCorners@data;
(*Test original data and point outside of convex hull*)
And @@ Table[TestCandidate[datac, data[[i]]], {i, Length@data}]
TestCandidate[datac, RandomReal[int + dis, d]]
True
False
As others mentioned, Qhull can do this. There are multiple ways to access Qhull from Mathematica.
One way is through the mPower package, now part of xCellerator. An answer discussing how to install the package is here: https://mathematica.stackexchange.com/a/18909/12
Another way is the qh-math package, described in this answer: https://mathematica.stackexchange.com/a/13445/12
Yet another way is through MATLAB and MATLink (MATLAB has Qhull built in).
<<MATLink`
OpenMATLAB[]
convhulln = MFunction["convhulln"]/*Round
pts = RandomReal[1, {20,5}] (* 20 points in 5D *)
convhulln[pts]
This returns a list of facets encoded as point-index-lists.
Loren Petrich has done a wonderful job of implementing quite a few convex hull algorithms in Mathematica, and you can get all of his original code from his website. Since it is written under the MIT license I took the liberty of wrapping the code relevant to an N-dimensional convex hull into a package.
Included are the functions CHNGiftWrapping, CHNChanShattering, CHNQuickHull, and CHNIncremental, which use the gift wrapping algorithm, Chan's algorithm, quick hull, and the incremental algorithms, respectively.
All these functions work the same, they take a list of N-dimensional vertices and return a sorted list of vertices which make up the hull and a list of simplexes.
(* load the package *)
<<"https://gist.githubusercontent.com/jasondbiggs/c3d9410af3195da514a442be5b563ab8/raw/80ac5074077f0d4b1366540c8c710e35cb530ddd/NDConvexHull.m"
Here I apply the quick hull function to a list of 4-dimensional points:
SeedRandom[4];
list = RandomReal[5, {50, 4}];
CHNQuickHull[list]
(* {{2, 3, 7, 8, 11, 12, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
26, 27, 28, 29, 31, 32, 33, 34, 35, 36, 37, 38, 40, 41, 42, 45, 46,
47, 48, 49,
50}, {{2, 15, 21, 47}, {2, 15, 47, 50}, {2, 15, 50, 21}, {2, 21, 35,
47}, {2, 21, 50, 35}, {2, 22, 35, 50}, {2, 22, 47, 35}, {2, 22,
50, 47}, {3, 16, 17, 22}, {3, 16, 22, 49}, {3, 16, 49, 17}, {3, 17,
49, 22}, {7, 22, 29, 35}, {7, 22, 34, 49}, {7, 22, 35, 34}, {7,
22, 36, 29}, {7, 22, 49, 36}, {7, 29, 34, 35}, {7, 29, 36, 49}, {7,
29, 49, 34}, {8, 16, 17, 49}, {8, 16, 18, 24}, {8, 16, 24,
17}, {8, 16, 49, 18}, {8, 17, 20, 32}, {8, 17, 24, 20}, {8, 17, 32,
49}, {8, 18, 20, 24}, {8, 18, 23, 20}, {8, 18, 29, 37}, {8, 18,
36, 23}, {8, 18, 37, 36}, {8, 18, 49, 29}, {8, 20, 23, 33}, {8, 20,
33, 41}, {8, 20, 41, 32}, {8, 23, 36, 33}, {8, 29, 49, 37}, {8,
32, 33, 36}, {8, 32, 36, 49}, {8, 32, 41, 33}, {8, 36, 37,
49}, {11, 12, 15, 23}, {11, 12, 17, 40}, {11, 12, 18, 42}, {11, 12,
21, 45}, {11, 12, 23, 18}, {11, 12, 40, 21}, {11, 12, 42,
17}, {11, 12, 45, 15}, {11, 15, 18, 23}, {11, 15, 21, 26}, {11, 15,
26, 27}, {11, 15, 27, 18}, {11, 15, 45, 21}, {11, 16, 17,
24}, {11, 16, 19, 28}, {11, 16, 21, 40}, {11, 16, 24, 19}, {11, 16,
28, 21}, {11, 16, 40, 17}, {11, 17, 42, 24}, {11, 18, 19,
24}, {11, 18, 24, 42}, {11, 18, 27, 29}, {11, 18, 29, 19}, {11, 19,
29, 46}, {11, 19, 35, 28}, {11, 19, 46, 35}, {11, 21, 28,
35}, {11, 21, 35, 38}, {11, 21, 38, 26}, {11, 26, 38, 27}, {11, 27,
38, 46}, {11, 27, 46, 29}, {11, 35, 46, 38}, {12, 15, 23,
45}, {12, 17, 20, 42}, {12, 17, 40, 45}, {12, 17, 45, 20}, {12, 18,
42, 23}, {12, 20, 23, 42}, {12, 20, 45, 23}, {12, 21, 45,
40}, {15, 18, 23, 27}, {15, 21, 26, 50}, {15, 21, 47, 45}, {15, 23,
36, 27}, {15, 23, 45, 48}, {15, 23, 48, 36}, {15, 26, 27,
38}, {15, 26, 38, 50}, {15, 27, 36, 50}, {15, 27, 50, 38}, {15, 31,
36, 47}, {15, 31, 47, 50}, {15, 31, 50, 36}, {15, 36, 48,
47}, {15, 45, 47, 48}, {16, 17, 21, 40}, {16, 17, 22, 47}, {16, 17,
47, 21}, {16, 18, 19, 29}, {16, 18, 24, 19}, {16, 18, 29,
49}, {16, 19, 21, 35}, {16, 19, 28, 21}, {16, 19, 34, 29}, {16, 19,
35, 34}, {16, 21, 47, 35}, {16, 22, 34, 35}, {16, 22, 35,
47}, {16, 22, 49, 34}, {16, 29, 34, 49}, {17, 20, 32, 47}, {17, 20,
42, 24}, {17, 20, 47, 45}, {17, 21, 40, 47}, {17, 22, 32,
49}, {17, 22, 47, 32}, {17, 40, 45, 47}, {18, 20, 24, 42}, {18, 20,
42, 23}, {18, 23, 27, 37}, {18, 23, 37, 36}, {18, 27, 29,
37}, {19, 21, 35, 28}, {19, 29, 35, 34}, {19, 29, 46, 35}, {20, 23,
33, 41}, {20, 23, 41, 45}, {20, 32, 47, 41}, {20, 41, 47,
45}, {21, 26, 50, 38}, {21, 35, 38, 50}, {21, 40, 47, 45}, {22, 29,
35, 36}, {22, 31, 36, 50}, {22, 31, 47, 36}, {22, 31, 50,
47}, {22, 32, 36, 47}, {22, 32, 49, 36}, {22, 35, 50, 36}, {23, 27,
37, 36}, {23, 33, 41, 36}, {23, 36, 41, 48}, {23, 41, 45,
48}, {27, 29, 36, 38}, {27, 29, 37, 36}, {27, 29, 38, 46}, {27, 36,
50, 38}, {29, 35, 36, 38}, {29, 35, 38, 46}, {29, 36, 49,
37}, {32, 33, 36, 41}, {32, 36, 47, 48}, {32, 36, 48, 41}, {32, 41,
48, 47}, {35, 36, 38, 50}, {41, 45, 48, 47}}} *)
You can verify that they all produce the same results:
CHNAlgorithms = {CHNGiftWrapping, CHNChanShattering,
CHNQuickHull, CHNIncremental};
SameQ[Through[CHNAlgorithms][list]]
(* True *)