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