Convexhull


.

It is well known that for every finite set of points on a plane, there exists a minimal convex polygon such that it encompass all the points (by encompass, we allow the possibility of points being vertices and on the edges of the polygon), such a polygon is called the ‘convex hull’. Given the co ordinates of a set of points, there are several algorithms to obtain the hull (by hull, I mean a subset of the points that form the outline of the convex gon). Along with a friend of mine, I wrote a program in mathematica implementing gift wraping algorithm, which is simple but not optimal (It is ~n^2 compared to others which are ~n log(n))

Although the algorithm was implemented, the code can be further optimised.

srikconvexhull2[points_, precision_: 10^-10] :=
  (* set precision smaller than least difference between any numbers \
in Flatten[points] *)
  Module[{orgpt = {}, pt = {}, anglelist = {}, newcollinearpts = {}, 
    newpt = {}, hull = {}, oldpt = {}, selpoints = {}},

   If[Length[points] <= 3, Return[points]];
   (* send off *)

   orgpt = Part[points, Last[Ordering[points]]];
   (* getting right-top point using ordering *)

   pt = orgpt;
   (* initialising pt to use, and to keep orgpt safe *)

   hull = {orgpt};
   (* we start our hull with right-top point *)

   oldpt = {orgpt[[1]], orgpt[[2]] - precision};
   (* setting up initial old pair where the base(old) point is just \
below the orgpt *)

   While[True,

    If[Sort[hull] == Sort[points], Return[hull]];
    (* this is for the case when all points are on a vertical / 
    horizontal line when our while breaking condition fails *)

    selpoints = setminus[points, {oldpt, pt}];
    (* avoids two previous points *)

    anglelist = positivelineangle[{{oldpt, pt}, #}] & /@ selpoints;
    (*  getting the list of all angles *)

    newcollinearpts = 
     selpoints[[#]] & /@ Flatten[Position[anglelist, Min[anglelist]]];
    (* The set of collinear points that have minimum angular \
deviation from the guiding line *)

    newpt = newcollinearpts[[closestpt[pt, newcollinearpts]]];
    (* getting the actual closet point *)

    If[hull[[1]] == newpt , Return[hull], AppendTo[hull, newpt]];
    (* unless we reach the first element of the hull, 
    we keep appending *)

    oldpt = pt;
    pt = newpt;
    (* setting up the things for the next iteration *)
    ];
   ];

(* support functions *)

(* gives the closest point in the list to pt.note that we are sending \
the address of the closet point and not the closest point itself *)
closestpt[pt_, list_] := 
 Ordering[EuclideanDistance[#, pt] & /@ list] // First

setminus[a_List, b_List] := Module[{new = a},
   If[MemberQ[a, #], new = DeleteCases[new, #]] & /@ b;
   Return[new];
   ];

(* ranges from 0 to 2\[Pi] *)
nanarctan[pt_] := If[pt == {0, 0}, Return[0],
  If[ArcTan[pt[[1]], pt[[2]]] >= 0, Return[ArcTan[pt[[1]], pt[[2]]]], 
   Return[2 \[Pi] + ArcTan[pt[[1]], pt[[2]]]]
   ]
  ]

(* pt1 moves to pt2 in anticlockwise direction, thinking both points \
are on the unit circle *)
positiveangle[{pt1_, pt2_}] := 
 If[nanarctan[pt1] <= nanarctan[pt2], nanarctan[pt2] - nanarctan[pt1],
   2 \[Pi] - (nanarctan[pt1] - nanarctan[pt2])]

(* reflect the point pt1 about the point pt2 on the straight line \
connecting them *)
reflect[pair_] := 2 pair[[2, #]] - pair[[1, #]] & /@ Range[2]

(* pair defines a line. how much positive angle do we have cover so \
that a imaginary leading point with the same slope as pair (and \
pair[[2]] as the base point) to reach pt *)
positivelineangle[{pair_, pt_}] := Module[
   {rbasept = {}, refpt = {}, rleadpt = {}, rpt = {}},
   rbasept = {0, 0};
   refpt = reflect[pair];
   rleadpt = refpt - pair[[2]];
   rpt = pt - pair[[2]];
   Return[positiveangle[{rleadpt, rpt}]];
   ];

(note that '\' usually means a line break)

Here is a visual output of a dataset found here on stackoverflow using the program:

bigdata.png

The red line corresponds to the polygon spanned by the output generated by the above program, green to the polygon output given by the mathematica function 'ConvexHull', part of the computational geometry package.

Please feel free to comment and suggest improvements in our code.

Better programs (mostly): One written by Eric W. Weisstein, pascoletti, Imtek mathematica supplement and asynadak.

Here is the nb file of the program.


29 Aug 2012 10:54

Comments: 0
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-NonCommercial-NoDerivs 3.0 License