Re: Much faster ConvexHull implementation
- To: mathgroup at smc.vnet.net
- Subject: [mg55600] Re: Much faster ConvexHull implementation
- From: "Jens-Peer Kuska" <kuska at informatik.uni-leipzig.de>
- Date: Thu, 31 Mar 2005 01:23:52 -0500 (EST)
- Organization: Uni Leipzig
- 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> <d2b547$79l$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Hi, you can have a MathLink Program based on Ken Clarksons convex hull/delauny triangulation code that that work well for dimensions up to 6. It works internal with 32 bit integers and your data should not have a a higher precision. Regards Jens "Carl K. Woll" <carl at woll2woll.com> schrieb im Newsbeitrag news:d2b547$79l$1 at smc.vnet.net... > 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