Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2005
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2005

[Date Index] [Thread Index] [Author Index]

Search the Archive

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



  • Prev by Date: Re: Re: point in convex hull
  • Next by Date: Re: front end complaint (ui design flaw?)
  • Previous by thread: Re: Re: point in convex hull
  • Next by thread: Re: Much faster ConvexHull implementation