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