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

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>