Much faster ConvexHull implementation
- To: mathgroup at smc.vnet.net
- Subject: [mg55556] Much faster ConvexHull implementation
- From: "Carl K. Woll" <carl at woll2woll.com>
- Date: Tue, 29 Mar 2005 03:42:40 -0500 (EST)
- References: <d1tvc0$rli$1@smc.vnet.net> <200503270742.CAA06233@smc.vnet.net> <opsoa9q5xpiz9bcq@monster.ma.dl.cox.net> <011401c53310$74dde680$6400a8c0@Main> <opsob4uhh3iz9bcq@monster.ma.dl.cox.net> <02a501c533ac$f76aa4c0$6400a8c0@Main> <opsoc793kbiz9bcq@monster.ma.dl.cox.net>
- Sender: owner-wri-mathgroup at wolfram.com
Hi all,
Inspired by the "point in a convex hull" discussion, I noticed that the
DiscreteMath`ComputationalGeometry` version of ConvexHull appears to scale
as O(n^2) instead of O(n log n) which is the usual lower bound. So, for
those who are interested, I will present a much faster convex hull
implementation that seems to scale as O(n log n). To get an idea of the
speed difference, on a random sample of 1000 points, my version of convex
hull is 60 times faster. Of course, with larger samples, the disparity is
even worse due to the different scaling behavior.
The algorithm is very simple. Given a line and some points above the line,
the function finds the point that is highest above the line. Then, it forms
a left line and a right line to this point, finds the points above these
lines and repeats. Do the same thing with a line and points below the line
and you have the convex hull.
The output of my version of convex hull gives the points in counterclockwise
order and is identical to the built-in ConvexHull, except for the case where
there are vertical lines at the right edge. In this case, my version of
convex hull starts with the top end of the right vertical line and goes
counterclockwise, whereas the built-in ConvexHull starts at the bottom end
of the right vertical line and goes counterclockwise. Also, if there are
vertical lines at the right or left edge with more than 2 points, then my
version omits these extra points, while the built-in ConvexHull includes
them.
At any rate, my version of convex hull can be found below. Any comments are
appreciated.
Carl Woll
It's best to make sure your default input format is InputForm when you copy
the function below to your notebook. At least on my machine, copying the
following code into a StandardForm input cell introduces invisible
multiplications so that executing the code results in Null^9 and the code
doesn't work.
convex[pts_] := Module[{spts, ss, toppts, bottompts},
spts = Sort[Transpose[{N[pts], Range[Length[pts]]}]];
ss = Drop[Split[spts[[All,1,1]]], {2, -2}];
If[spts[[Length[ss[[1]]],1]] === spts[[1,1]],
topleftindex = {};
topleft = spts[[1,1]]; ,
topleftindex = {spts[[Length[ss[[1]]],2]]};
topleft = spts[[Length[ss[[1]]],1]];
];
If[spts[[-Length[ss[[-1]]],1]] === spts[[-1,1]],
bottomrightindex = {};
bottomright = spts[[-1,1]]; ,
bottomrightindex = {spts[[-Length[ss[[-1]]],2]]};
bottomright = spts[[-Length[ss[[-1]]],1]];
];
topline = Interpolation[{topleft, spts[[-1,1]]}, InterpolationOrder -> 1];
bottomline = Interpolation[{spts[[1,1]], bottomright},
InterpolationOrder -> 1];
toppts = Cases[spts, {{x_, y_}, _} /; y - topline[x] > 0];
bottompts = Cases[spts, {{x_, y_}, _} /; y - bottomline[x] < 0];
Join[
Reverse[toppart[toppts, topline, Null, spts[[-1,2]]]],
topleftindex,
bottompart[bottompts, bottomline, spts[[1,2]], Null],
bottomrightindex
]
]
toppart[pts_, line_, l_, r_] := Module[{newpt, leftline, rightline, leftpts,
rightpts},
newpt = Ordering[pts[[All,1,2]] - line[pts[[All,1,1]]], -1][[1]];
leftline = Interpolation[{leftend[line], pts[[newpt,1]]},
InterpolationOrder -> 1];
rightline = Interpolation[{pts[[newpt,1]], rightend[line]},
InterpolationOrder -> 1];
leftpts = Cases[Take[pts, newpt - 1], {{x_, y_}, _} /; y - leftline[x] >
0];
rightpts = Cases[Drop[pts, newpt], {{x_, y_}, _} /; y - rightline[x] > 0];
Join[ toppart[leftpts, leftline, l, pts[[newpt,2]]],
toppart[rightpts, rightline, pts[[newpt,2]], r]
]
]
toppart[{{pt_, index_Integer}}, line_, l_, r_] := {index, r}
toppart[{}, line_, l_, r_] := {r}
bottompart[pts_, line_, l_, r_] := Module[{newpt, leftline, rightline,
leftpts, rightpts},
newpt = Ordering[pts[[All,1,2]] - line[pts[[All,1,1]]], 1][[1]];
leftline = Interpolation[{leftend[line], pts[[newpt,1]]},
InterpolationOrder -> 1];
rightline = Interpolation[{pts[[newpt,1]], rightend[line]},
InterpolationOrder -> 1];
leftpts = Cases[Take[pts, newpt - 1], {{x_, y_}, _} /; y - leftline[x] <
0];
rightpts = Cases[Drop[pts, newpt], {{x_, y_}, _} /; y - rightline[x] < 0];
Join[ bottompart[leftpts, leftline, l, pts[[newpt,2]]],
bottompart[rightpts, rightline, pts[[newpt,2]], r]
]
]
bottompart[{{pt_, index_Integer}}, line_, l_, r_] := {l, index}
bottompart[{}, line_, l_, r_] := {l}
leftend[interp_] := {#1, interp[#1]}&[interp[[1,1,1]]]
rightend[interp_] := {#1, interp[#1]}&[interp[[1,1,2]]]
- References:
- Re: point in convex hull
- From: "Carl K. Woll" <carlw@u.washington.edu>
- Re: point in convex hull