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

Re: Much faster ConvexHull implementation

  • To: mathgroup at smc.vnet.net
  • Subject: [mg55758] Re: [mg55556] Much faster ConvexHull implementation
  • From: "Wolf, Hartmut" <Hartmut.Wolf at t-systems.com>
  • Date: Tue, 5 Apr 2005 05:45:08 -0400 (EDT)
  • Sender: owner-wri-mathgroup at wolfram.com

Dear Carl,

just saw the discussion, and (without deep thinking) tried another
simple idea: sort the points, identify the maximum (the point with max
x-coordinate), and then, by sort order, incrementally build the "top-"
and "bottomlines" (we have a convex hull for each set of points
considered so far).  This gives "stupid" but compact code.  It is slower
than yours by a factor of 3, but compiles well and then of course is
faster.

A quick test over a range of up to 200000 random points shows that it
scales similar to your function.

Here is the compiled version:

cch = Compile[{{data, _Real, 2}, {ord, _Integer, 1}}, 
    Module[{high, low, i, j, J}, J = {{0, 1}, {-1, 0}}; 
      {high, low} = If[data[[ord[[-2]],2]] >= 
         data[[ord[[-3]],2]], {ord[[{-1, -2}]], 
         ord[[{-1, -3}]]}, {ord[[{-1, -3}]], 
         ord[[{-1, -2}]]}]; i = Length[data] - 2; 
      While[--i > 0, For[j = Length[high], 
         j > 1 && (data[[high[[j]]]] - data[[high[[j - 
                1]]]]) . J . (data[[ord[[i]]]] - 
             data[[high[[j]]]]) <= 0, --j]; 
        high = Append[Take[high, {1, j}], ord[[i]]]; 
        For[j = Length[low], j > 1 && 
          (data[[low[[j]]]] - data[[low[[j - 1]]]]) . J . 
            (data[[ord[[i]]]] - data[[low[[j]]]]) >= 0, --j]; 
        low = Append[Take[low, {1, j}], ord[[i]]]]; 
      Join[Drop[high, -1], Drop[Reverse[low], -1]]]]; 

ch2[data_] := cch[data, Ordering[data]]

The result is given in the same way as by your function and ConvexHull.

This is not meant very seriously (and I won't try to improve on the
Append-s and Take-s), but I think divide and conquer solutions will only
be better at very large sets of data.


--
Hartmut



>-----Original Message-----
>From: Carl K. Woll [mailto:carl at woll2woll.com] 
To: mathgroup at smc.vnet.net
>Sent: Tuesday, March 29, 2005 10:43 AM
>Subject: [mg55758] [mg55556] Much faster ConvexHull implementation
>
>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: NMinimize--problem with a min-max problem
  • Next by Date: Re: file corrupted
  • Previous by thread: Re: Much faster ConvexHull implementation
  • Next by thread: transpose