MATHEMATICA IN ACTION for issue 9-2 of Mathematica in Education and Research

Random Polygons

A robust algorithm to connect the dots.

by Stan Wagon

. Introduction

The problem of generating random polygons has generated a lot of interest in the computational geometry community. In this column we will describe and implement the Auer–Held steady-growth heuristic for generating random-looking polygons. The project is a nice vehicle for demonstrating some of the basic tools of computational geometry, and also for discussing techniques to speed up a numerically heavy program. While the algorithm is a little slow for large numbers of points, it has several virtues:

1. It is acceptably fast for generating n-gons where n ≤ 50.

2. It is a method for getting a polygon through any set of points. Thus, to get an n-gon, one can just generate n random points and then apply the steady-growth method.

3. It is intuitively very easy to understand: start with two points and splice in the unused points in an order that guarantees one will never get stuck.

One approach to the problem of defining a random polygon is to fix the vertex-set S and define a method to be random if the probability of generating any particular polygon on S is (1/k), where k is the number of polygons on S. This should not be taken as a definitive definition, but it is useful. A method that respects this definition is the one that simply generates random permutations of S, stopping as soon as a good permutation — one that yields a simple polygon — is found. This will yield the desired uniform distribution, but it is too slow, since, when S has 20 or more points, the search for a good permutation can take a very long time. The steady-growth heuristic fails to have this uniformity, since there are some polygons (the smallest has 11 vertices) that steady growth will never generate [Auer 1996].

My interest in this subject was inspired by Will Carver and Peter Wiringa, students in an algorithms class at Macalester College who successfully implemented the steady-growth algorithm in Mathematica.

. The Steady-Growth Algorithm

Here are the steps of the Auer–Held steady-growth algorithm [Auer and Held 1996]. Assume that n points in the plane are given and let CH denote convex hull. The following procedure will connect the points so that they form a polygon 𝒫.

1. Let 𝒫 consist of two random points in the set; call the rest unused.

2. Do n - 2 times:

   2a. Find an unused point x so that CH(𝒫 ⋃ {x}) contains none of the other unused points.

   2b. With x as in step 2a, find an edge in 𝒫 so that x can see the entire edge. Extend 𝒫 to include x by splicing it in: delete the visible edge and introduce two new edges from the ends of the visible edge to x.

These ideas are quite simple, but there is an important subtlety: Given a polygon 𝒫 and a point x outside it, it is not necessarily true that x can see an edge of 𝒫. A counterexample is shown in Figure . However, if x is outside CH(𝒫), then there must be a visible edge. The proof is not difficult [Auer and Held 1996]: given x, look at the chain facing x from the leftmost to the rightmost vertex of 𝒫 that it can see (the existence of such a chain uses the convex hull hypothesis). If the leftmost edge in this chain is visible we are done; if not, consider the blocking edges and look at the first blocking edge whose leftmost vertex is visible to x. If the edge is entirely visible, we are done. If not, look at the blocking edges for this one, and continue until a visible edge is found.

The basic ideas of the steady-growth algorithm are simple, but getting an acceptably fast implementation requires several refinements. For example, it is tempting to use the ConvexHull routine in the DiscreteMath`ComputationalGeometry package, but that would be wasteful since it is possible to implement step 2a without using a convex hull routine.

. Computational Geometry Basics

Almost any problem in computational geometry requires fast routines to do certain basic things. The real workhorse is a routine that determines the orientation of a triangle ABC (which is taken to be +1 if the order is counterclockwise, -1 if clockwise, and 0 if the points are collinear). We will not prove it here (see [O’Rourke 1998]), but there is a very simple formula for the orientation in terms of the coordinates of the points: it is the sign of the determinant of (({{(Subscript[x, 1]), (Subscript[y, 1]), 1}, {(Subscript[x, 2]), (Subscript[y, 2]), 1}, {(Subscript[x, 3]), (Subscript[y, 3]), 1}})).  In fact, half of this determinant gives the signed area of the triangle. We use Chop so that three points that are collinear are properly interpreted as being so.

Orientation[{{x1_, y1_}, {x2_, y2_}, {x3_, y3_}}] :=  
      Sign[Chop[-x2 y1 + x3 y1 + x1 y2 - x3 y2 - x1 y3 + 
            x2 y3]]

This routine is a real workhorse, since it underlies so many others. For example, to decide if point R is left of the directed line form P to Q, one need only look at the orientation of PQR.

LeftOf[p_, q_, r_]  :=  Orientation[{p, q, r}] > 0

And now we can decide if two line segments intersect in their interior by checking various LeftOfs as follows.

\(IntersectInteriorQ[{a_, b_}, {c_, d_}] := 
	
      Xor[LeftOf[a, b, c], LeftOf[a, b, d]] && 
	
        Xor[LeftOf[c, d, a], LeftOf[c, d, b]] && 
	
        Intersection[{a, b}, {c, d}] == {}\)

These ideas also lead to a quick way of determining whether a point x lies inside Δpqr. It suffices to check that the three orientations from x to the edges are the same.

\(InsideTriangleQ[{p_, q_, r_}, x_] := \ 
      3 == 
        Abs[
	Orientation[{p, q, x}] + 
            Orientation[{q, r, x}] + 
\ \ \ \ \ \ \ \ \ \ \ Orientation[{r, 
                p, x}]]\)

All these routines can be made faster by compiling, and I have done so in the RandomPolygon package that accompanies this column. A final interesting general problem is to find the maximum and/or minimum of a set S of points based on some function on S. The Sort command allows a second argument that specifies a sorting function: Sort[s, f[#1] < f[#2] &]. This works for the problem at hand but is a little slow since the function is evaluated at each sorting step. It is quite a bit faster to form the f-values of the elements of S, take the max (or min) of the result, find its position in the set of f-values, and then get the corresponding element of S. Here is a routine that finds both the min and max. A cute trick is the use of Position[s, m, 1, 1]; the fourth 1 indicates that the search should stop when 1 is found. The third argument, 1, is a level specification that is necessary here only as a placeholder, so that we can use a fourth argument.

\(\(MinAndMax[s_, 
          f_] := \((s[[{
\ \ \ \ \ \ \ \(Position[#, 
                      Min[#], 1, 1]\)[[1, 
                    1]], 
\ \ \ \ \ \ \ \(Position[#, 
                      Max[#], 1, 1]\)[[1, 
                    1]]}]] &)\)[
          f /@ s];\)\)

\(MinAndMax[{\(-2\), \(-3\), 1, \(-5\), 0, 8, \(-6\)}, #\^2 &]\)

A final geometric routine we need is one that computes the angle formed at x by the points p and q. It is convenient to think in terms of two vectors v and w; we want signed angles and use the sign of an orientation to get it. Of course, the formula underlying this is the familiar \(TraditionalForm\`v·
          w = \(\(|\)\(v\)\(|\)\(|\)\(w\)\(|\)\(\ \)\(cos\ θ\)\)\).

angle[v_, w_] := 
        Sign[v[[1]] w[[2]] - 
              v[[2]] w[[1]]] ArcCos[
            v . (w/(Sqrt[v . v  w . w]))];

And we need a utility to add the first element of a set to the end.

extend[s_] := 
      Append[s, s[[1]]]

. The Main Routines

Now we can turn our attention to efficient ways of implementing the main steps (2a and 2b) of the steady-growth algorithm. Here is one observation that allows us to find the new point quickly, without computing any convex hulls. If 𝒫 is the current polygon, and x one of the unused points, we first identify the two supporting vertices of 𝒫 (the two vertices at the extreme ends of x’s vision). This is done by choosing an arbitrary point of 𝒫 (poly[[1]] in the code that follows) and taking the max and min of the directed angles formed by x, the special point, and the other points of the polygon.

\(\(AdditionalTriangle[poly_, x_] := 
        MinAndMax[
\ \ \ \ \ Range[Length[poly]], 
          angle[poly[[1]] - x, 
              poly[[#]] - 
                x] &];\)\)

samplepoly = {{0, 0}, {10, 0}, {5, 3}, {20, 6}, {5, 7}, {10, 9}, {0, 
            10}, {2, 2}};
x = {15, 14};
AdditionalTriangle[samplepoly, x]

This output indicates that the fourth and second vertices are the two supporting vertices (see Figure ). The reason this is called AdditionalTriangle is that the triangle formed by these two vertices and the point x is the new part of the convex hull of 𝒫 ⋃ {x}. Since no unused points lie in the convex hull of 𝒫, it is sufficient to check only this new triangle to see if it contains any other unused points. If it does not, then x is the new point. But if it contains only one of the unused elements, then we are also done, since we may choose that one element to be the new point! This case is handled by the If statement in the code. The routine that follows accepts a pair {𝒫, unused} and returns the pair: {𝒫  ⋃  {x},  
        unused points less x}. Note how a third argument to Select is always used — it is 1 or 2 — so that the search stops as soon as the necessary information is obtained. This setup allows us to quickly get the entire order of the selection of the unused points by simply iterating the NewPoint function n - 2 times, where n is the size of the initial point set.

\(\(NewPoint[{poly_, unused_}] := 
        Module[{
	p, tri, inside}, 
	
          p = First[
              Select[unused, 
		\((tri = 
                      Append[poly[[
                          AdditionalTriangle[
                            poly, #]]], #]; 
		Length[
                        inside = 
                          Select[
                            DeleteCases[
                              unused, #], 
			InsideTriangleQ[tri, #] &, 
                            2]] < 2)\) &, 1]]; 
	
          If[Length[inside] == 1, 
            p = inside[[1]]]; 
	{Append[poly, p], DeleteCases[unused, p]}];\)\)

n = 12;
pts = Array[Random[] &,  {n, 2}];
Flatten[Position[pts, #] &  /@  
First[
          Nest[NewPoint, {Take[pts, 2],  Drop[pts, 2]}, n - 2]]]

Now we turn to the problem of finding a visible edge. The VisibleEdge routine relies on the visibleQ subroutine. The former looks at the edges of the polygon and sorts them in order of their angle around x. Then it looks for the first edge that is visible. This sorting, suggested by Auer and Held, speeds up the search for a visible edge. The visibleQ routine accepts blockers (the set of potential blocking edges), {p,q} (the edge in question), and x (the new point), and returns True if the edge is not blocked by any of the blockers. Here I use some knowledge about the blockers to speed up the check; this is a bottleneck of the program, so such optimization is a good idea. First note that the blockers are defined in the same order around the polygon as the edge pq. So if the edge pq is blocked by a blocker, then one of the following must hold:

(1) the second vertex of the blocker is inside Δpqx;
(2) the blocker intersects the segment xp;
(3) the blocker has the form rp and the blocker intersects xq.

VisibleEdge[pol_, x_] := 
        Module[{edges, epol = extend[pol]}, 
          If[Length[pol] == 2, pol, 
            edges = Map[
                epol[[{#, # + 1}]] &, 
                Last /@ Sort[
                    Transpose[{Apply[ArcTan, (# - x &) /@ pol, {1}], 
                        Range[Length[pol]]}]]]; 
Select[
                edges, visibleQ[DeleteCases[edges, #], #, x] &, 
                1][[1]]]];

\(visibleQ[blockers_, {p_, q_}, x_] := {} == 
	
        Select[blockers, 
		InsideTriangleQ[{p, q, 
                  x}, #[[2]]] || 
		IntersectInteriorQ[#, {x, 
                  p}] || 
		\((p == #[[2]] && IntersectInteriorQ[#, {x, q}])\) &, 1]\)

Here is a sample output (see Figure ), where we replace the points by their indices in the sample polygon.

Flatten[Position[samplepoly, #] &  /@  
        VisibleEdge[samplepoly, x]]

Now we can define the main routine. Using Fold is efficient. We use it by first getting the order of the points to be spliced in, and then splicing them in one at a time, using Fold. Note the use of the intricate (but fast) pure function

 \(\((a = \(VisibleEdge[#1, #2]\)[[1]]; 
\ \ #1 /. {x___, a, 
                y___} :> {x, a, #2, y})\) &\). 
 
Here #1 will be a polygon and #2 the point to be spliced in. The first vertex of a visible edge is stored in a, and then a replacement is used to splice in #2 just after a.

\(\(RandomPolygon[n_Integer]\  := \ 
        RandomPolygon[Array[Random[] &, \ {n, 2}]];\)
\)
\(\(RandomPolygon[pts_] := 
        With[{order = 
	\(Nest[NewPoint, {Take[pts, 2], Drop[pts, 2]}, 
                  Length[pts] - 
                    2]\)[[1]]}, 
	
          Fold[\((
		a = \(VisibleEdge[#1, #2]\)[[1]]; 
		#1 /. {x___, a, y___} :> {x, a, #2, 
                      y})\) &, 
		\ \ \ \ \ \ Take[order, 2], 
            Drop[order, 2]]];\)\)

The output is just the list of points. The timing here is not too impressive, but the use of compilation, discussed briefly in section 5 and used several times in the RandomPolygon package, cuts the running time in half.

RandomPolygon[30] // Short // Timing

Figure  shows several of the outputs of RandomPolygon.

. A Comprehensive Package

The RandomPolygon package in the electronic supplement is more sophisticated than the basic routines presented here. Here are some of its features:

• A ShowSteps option creates a movie showing the addition of each successive point.

• A smoothing option causes the polygon to be smoothed: each point is replaced by the two points that are 25% of the way along the two edges emanating from it. This doubles the number of points, but eliminates the sharp edges, whichi leads to more visually pleasing polygons (see Figure ). The constant, 0.25, can be changed by the SmoothingParameter option. It is possible that such smoothing will lead to a self-intersection, and an error trap (which can be suppressed) checks for that. The implementation of this option is via SmoothingIterations, which can be 1, for a single smoothing, or a larger integer, for repeated smoothing.

• A Method option allows two methods to be used: SteadyGrowth (the default) and PermuteAndReject. The latter runs through permutations of the points until a good permutation is found. This works only for small sets of points, but it has the benefit that each polygon on the points is equally likely to be generated.

• A SmoothCurve option causes interpolation to be used to generate smooth function for the x- and y-coordinates. This works better if Smoothing is used as well, but there is no check to see if the result is a simple curve and the results are generally not very satisfactory.

• A ShowPolygon option to toggle the display of the polygon; PolygonColor, ShowPoints, PointStyle, BorderStyle, and Background options to control the style of various aspects of the displayed polygon.

• It includes a triangulation routine (using diagonals only) using the most basic basic triangulation method as discussed in chapter 23 of [Wagon 1999].

• It includes an Angles function that returns the list of angles in a polygon. This allows one to do experiments on the distribution of angles. And there is an AngleInterval option that causes the algorithm to generate polygons until it finds one whose angles lie in the specified interval.

• When the input to RandomPolygon is a set of points, they are permuted right away by a random permutation. This guarantees that different runs on the same point set will yield a different polygon.

We now present several applications of these options. First, here is how triple smoothing can yield a more visually pleasing, if less random, polygon (see Figure ). The check for self-intersection is slow, so we turn it off.

SeedRandom[2];  RandomPolygon[25];
SeedRandom[2];  
    RandomPolygon[
      25, 
SmoothingIterations -> 
        3,  
CorrectnessCheckWhenSmoothing -> 
        False];

Here is an example of the triangulation routine, with the triangles colored different shades. The polygon and triangulation are shown in Figure , where a simple algorithm (the triangles of a triangulation always form a tree) was used to 2-color the triangles.

poly = RandomPolygon[15];
Triangulate[poly]

The PermuteAndReject method simply tries random permutations until a good one — one that yields a polygon — is found. This is too slow for any serious work, but it does allow one to look at truly random n-gons when n is small. We can compare the two methods by looking at the angles generated in a large set of trials. If n is small, say n = 6, then, not surprisingly, the angle distributions generated are quite similar; in other words, for such a small value of n it seems as if seady growth is coming close to generating a truly random polygon. But for larger n, things are different. The results of an experiment involving 1000 different 12-gons (but with a careful use of SeedRandom to guarantee that the same sets of points are used for each pair of trials) are shown in Figure , with the angles gathered into bins of size 10°. There is a statistically significant difference, with the steady-growth method showing fewer angles at the extremes, and many more in the flat range (120°–180°).

However, things change as n grows. Figure  shows the result of an experiment that fixes five sets of 50 points and then generates 40 polygons on each set and bins the resulting 10000 angles with a bin-size of 10°. The resulting histograms show definite spikes near 0  °, 180  °, and 360  °. The results are consistent with a similar experiment carried out by Auer and Held.

d = 40; n = 50;
pts = Table[{Random[], Random[]}, {5}, {n}];

angledata = 
        Flatten[
Table[
            Angles[RandomPolygon[pts[[j]], 
                ShowPolygon -> False]], {j, 5}, {i, d}]];

Needs["Statistics`DataManipulation`"];

binsize = 10;
Histogram[angledata, 
        Frame -> True, 
HistogramScale -> 
          100 binsize, 
HistogramCategories -> 
          360/binsize];

Finally, the reader interested in visually appealing polygons might want to restrict the possible angles. This can be done with the AngleInterval option, which can be used to restrict the angles to, say, the interval {10, 350}. Of course, it then takes more time to generate polygons; further, one must be careful to match the angle interval to the number of vertices to avoid a too-small probability of success. Figure  shows an angle-restricted 25-gon.

RandomPolygon[25, 
        AngleInterval -> {10, 
            350}, 
ShowPoints -> False];

Several of the routines in the package were compiled for speed. Here is an example involving LeftOf: the three arguments, p, q, r, are specified to be real of depth 1, so p[[i]] can then be used in the compiled definition.

LeftOf = 
        Compile[{{p, _Real, 1}, {q, _Real, 1}, {r, _Real, 1}}, 
          p[[2]] (r[[1]] - 
                    q[[1]]) + 
              q[[2]] (p[[1]] - 
                    r[[1]]) + 
              r[[2]] (q[[1]] - 
                    p[[1]]) > 
            0];

LeftOfSlow[p_, q_, r_] := 
        p[[2]] (r[[1]] - 
                  q[[1]]) + 
            q[[2]] (p[[1]] - 
                  r[[1]]) + 
            r[[2]] (q[[1]] - 
                  p[[1]]) > 
          0;

Here are some timings to show how significant the speedup is.

{p, q, r} = {{0.1, 0.2}, {0.3, 0.4}, {0.9, 0.6}};
First[Timing[Do[LeftOf[p, q, r], {(10^4)}]]]

First[Timing[Do[LeftOfSlow[p, q, r], {(10^4)}]]]

There are other algorithms in the literature that the serious creator of random polygons will want to look at (the web site in [Auer and Held 1996] is a good place to start). But for rudimentary testing of polygon algorithms on n-gons with 100 or fewer vertices, the steady-growth algorithm seems adequately robust and fast.

References

[Auer 1996] T. Auer, Heuristics for the generation of random polygons, Master’s Thesis, Univ. Salzburg, Salzburg, Austria, 1996 .

[Auer and Held 1996] T. Auer and M. Held, Heuristics for the generation of random polygons, in Proceedings of the Eighth Canadian Conference on Computational Geometry, Ottawa, Carleton Univ. Press, 1996, 38–44 (additional material at www.cosy.sbg.ac.at/~held/projects/rpg/rpg.html).

[O’Rourke 1998] J. O’Rourke, Computational Geometry in C, 2nd ed., Cambridge Univ. Pr., New York, 1998.

[Wagon 1999] S. Wagon, Mathematica in Action, 2nd ed., Springer/TELOS, New York, 1999.

About the Author

Stan Wagon uses Mathematica extensively in his research, teaching, and exposition, and is especially appreciative of the new ways Mathematica allows us to look at mathematical objects. He is author or coauthor of several books, including Animating Calculus, VisualDSolve, Which Way Did the Bicycle Go?, the just-released and much expanded second edition of Mathematica in Action, and the recently released A Course in Computational Number Theory. He teaches a high-altitude course, Rocky Mountain Mathematica, every summer in the mountains of Colorado, and has used Mathematica to help design geometric entries in snow sculpture contests.

Stan Wagon
Department of Mathematics and Computer Science
Macalester College
St. Paul, MN 55105
wagon@macalester.edu

The RandomPolygon` Package

If[Names[#]=!={}, Remove[#]] &  /@  {"Global`extend","Global`InitialStep",
      "Global`PermuteAndReject","Global`PolygonQ","Global`RandomPolygon",
      "Global`ShowPoints","Global`ShowPolygon","Global`ShowSteps",
      "Global`Smooth","Global`SmoothCurve",
      "Global`SmoothingParameter","Global`SteadyGrowth","Global`Triangulate"};

BeginPackage["RandomPolygon`", 
 {"Utilities`FilterOptions`", "DiscreteMath`ComputationalGeometry`",
"DiscreteMath`Permutations`"}]; 
 
Off[General::spell1];

CorrectnessCheckWhenSmoothing::usage = "CorrectnessCheckWhenSmoothing is an option to RandomPolygon that, when False, turns off the time-consuming check for self-intersections.";

AngleInterval::usage = "AngleInterval is an option to RandomPolygon[n] that, when set to {a1, a2}, causes the algorithm to generate polygons until it has found one whose angles lie between a1 and a2. If smoothing is used, that is done before the angle computation.";

BorderStyle::usage = "BorderStyle is an option to RandomPolygon that sets the style of the border. A setting of None causes no border to appear.";

PolygonColor::usage = "PolygonColor is an option to RandomPolygon that sets the color of the displayed polygon.";

SmoothingIterations::usage = "SmoothingIterations is an option to RandomPolygon that sets the number of times the smoothing operation should be carried out. Each iteration doubles the number of vertices by replacing each vertex with the two vertices that are 1/4 of the way to their neighbors. The parameter, 0.25, can be changed by the Smoothing Parameter option. The effect is to eliminating some of the sharp edges.";

Triangulate::usage = "Triangulate[poly] returns a triangulation of the polygon, using indices for the triangles.";

Angles::usage = "Angles[p] returns the list of angles in the polygon p, measured in degrees.";

SteadyGrowth::usage = "SteadyGrowth is a possible setting to the Method option for RandomPolygon, and calls upon the Auer-Held steady-growth algorithm.";

PermuteAndReject::usage = "PermuteAndReject is a possible setting to the Method option for RandomPolygon. This generates random permutations until a polygon is reached.";

RandomPolygon::usage = "RandomPolygon[n] gives a random polygon on n points in the unit square, via the Auer-Held steady-growth algorithm. RandomPolygon[pts] gives a random polygon using the points pts."; 

extend::usage = "extend[a] appends the first element of a to the end of a."; 

Smooth::usage = "Smooth[polygon, c] gives a polygon with twice as many points that is smoother than the original, c is the smoothing parameter, which defaults to 4. There is a chance this will introduce self-intersections."; 

InitialStep::usage = "InitialStep is an option to ShowPolygon that, when ShowSteps is True, specifies the starting point of the movie.";

ShowPoints::usage = "ShowPoints is an option to RandomPolygon that, when ShowPolygon is True, determines whether the vertices are shown.";

ShowPolygon::usage = "ShowPolygon is an option to RandomPolygon that causes the polygon to be shown graphically."; 

SmoothingParameter::usage = "SmoothingParameter is an option to RandomPolygon that sets the constant used for the smoothing process. The default, 0.25, replaces each point by the points 1/4 of the way along to the neighbors. The value must be less than 0.5."; 

PolygonQ::usage = "PolygonQ[p] returns True if p is a simple polygon. It is assumed that p is a list of distinct points with no superposition among its edges."; 

RandomPolygon::badwarn = "The smoothing process has led to a nonpolygon; a different random polygon is being generated.";

SmoothCurve::usage = "SmoothCurve is an option to RandomPolygon that causes the polygon to be interpolated into (fx,fy) that defines a smooth curve. However, there is a good chance the curve will have self-intersections."; 

ShowSteps::usage = "ShowSteps is an option to RandomPolygon that, when True, causes each step of the construction to be shown.";

RandomPolygon::notrac = "The ShowSteps option can be used only when the method chosen is SteadyGrowth.";

PointStyle::usage = "PointStyle is an option to RandomPolygon that specifies the style of the vertices.";

Begin["`Private`"];

If[$VersionNumber < 4, Unprotect[Mod];
   Mod[a_, n_, 1] := Mod[a, n] /. 0 -> n; Protect[Mod]]; 

Angles[poly_] := (epol=Prepend[extend[poly], poly[[-1]]];
  1/Degree Table[
  angle[epol[[i]]-epol[[i+1]], epol[[i]]- epol[[i-1]]]/. 
          xx_?Negative :> 2 Pi + xx, {i, 2, Length[epol]-1}]);

MinAndMax[x_, f_] := (
   x[[{Position[#, Min[#], 1, 1][[1,1]], 
       Position[#, Max[#], 1, 1][[1,1]]}]]&)[f /@ x]; 
       
LeftOf = Compile[{{p,_Real,1}, {q,_Real,1}, {r, _Real, 1}},
  p[[2]](r[[1]]-q[[1]]) + q[[2]](p[[1]]-r[[1]]) + r[[2]](q[[1]]-p[[1]]) > 0];

extend[s_] := Append[s, s[[1]]]; 

norm = Compile[{{v, _Real, 1}}, Sqrt[v . v]]; 

InsideTriangleQ = Compile[{{q, _Real, 2}, {p, _Real, 1}}, 
 3 == Abs[
  Sign[p[[2]](q[[1,1]] - q[[3,1]]) + 
       q[[3,2]](p[[1]] - q[[1,1]]) + 
       q[[1,2]](q[[3,1]] - p[[1]])] + 
  Sign[p[[2]](q[[2,1]] - q[[1,1]]) +
       q[[1,2]](p[[1]] - q[[2,1]]) +
       q[[2,2]](q[[1,1]] - p[[1]])] + 
  Sign[p[[2]](q[[3,1]] - q[[2,1]]) +
       q[[2,2]](p[[1]] - q[[3,1]]) +
       q[[3,2]](q[[2,1]]- p[[1]])]] ];          

IntersectInteriorQ[{a_, b_}, {c_, d_}] := 
  Xor[LeftOf[a, b, c], LeftOf[a, b, d]] && 
  Xor[LeftOf[c, d, a], LeftOf[c, d, b]] && a != c && 
   a != b && b != c && b != d

angle = Compile[{{v, _Real, 1}, {w, _Real, 1}}, 
    Sign[v[[1]] * w[[2]] - v[[2]] * w[[1]]]*
    ArcCos[v . w / Sqrt[v.v * w.w]]]; 
    
AdditionalTriangle[poly_, p_] := 
  MinAndMax[Range[Length[poly]], angle[poly[[1]] - p, poly[[#]] - p] &]; 
 
VisibleEdge[pol_, pt_] := Module[{edges, epol = extend[pol]},
  If[Length[pol] == 2, Return[pol],
    edges = Map[epol[[{#, # + 1}]]&, Last /@ 
       Sort[Transpose[{
          Apply[ArcTan, (# - pt &) /@ pol, {1}], Range[Length[pol]]}]]];
  Select[edges, visibleQ[DeleteCases[edges, #], #, pt]&, 1][[1]]]];

visibleQ[blockers_, {p_,q_}, pt_] := {} == 
    Select[blockers,
    (* the next conditions use in several places the fact that the
       edges in blockers are in the same order as {p, q} *)
      InsideTriangleQ[{p, q, pt}, #[[2]]] ||
      IntersectInteriorQ[#, {pt, p}] || 
      (p == #[[2]] && IntersectInteriorQ[#, {pt, q}])&, 1]
    


NewPoint[{poly_, unused_}] := Module[{p, tri, s, n = Length[unused]},
    p = Select[Range[n],
    (tri = Append[poly[[AdditionalTriangle[poly, unused[[#]]]]],
     unused[[#]]];
    Length[s = Select[Delete[Range[n],#],
   InsideTriangleQ[tri, unused[[#]]]&, 2]] < 2)&, 1][[1]]; 
      If[Length[s] == 1, p = s[[1]]];
      {Append[poly, unused[[p]]], Delete[unused, p]}];

RandomPolygonPermuteAndReject[pts_] := (ans = pts[[DiscreteMath`Permutations`RandomPermutation[Length[pts]]]];
 While[!PolygonQ[ans], ans = pts[[DiscreteMath`Permutations`RandomPermutation[Length[pts]]]]];
 If[Orientation[ans]==-1, Reverse[ans], ans]);

RandomPolygonRaw[pts_] := (ans=Fold[(
   a = VisibleEdge[#1, #2][[1]]; 
   #1 /.{x___, a, y___} :> {x, a, #2, y}) &, 
    {#[[1]], #[[2]]}, Drop[#, 2]]& [
   (* order is here for use with ShowSteps option *)
order = Nest[NewPoint,{Take[pts,2], Drop[pts,2]}, Length[pts]-2][[1]]];
If[Orientation[ans]==-1, Reverse[ans], ans]);
  
Options[RandomPolygon] = {ShowSteps -> False, 
   ShowPoints -> False, InitialStep -> 1,
    SmoothingParameter -> 0.25, SmoothCurve -> False,
    SmoothingIterations->0, Background->GrayLevel[1],
    ShowPolygon -> True, Method -> SteadyGrowth, 
    PolygonColor->GrayLevel[0.6],BorderStyle->{},
    AngleInterval -> {0,360},
    CorrectnessCheckWhenSmoothing->True,
    PointStyle ->{PointSize[0.01], RGBColor[1,0,0]} }; 

RandomPolygon[n_Integer /; n >= 3, opts___] := Module[{m,mm,s},
  ai = AngleInterval /. {opts} /. Options[RandomPolygon];
  s = Random[Integer, 10^6];
  If[ai != {0, 360},
    {m, mm} = Sort[Angles[
       RandomPolygon[Array[Random[]&, {n, 2}], ShowPolygon->False]]][[{1,-1}]]; 
    While[!(ai[[1]] <= m < mm <= ai[[2]]), 
      s++;SeedRandom[s];
       {m, mm} = Sort[Angles[
   RandomPolygon[pts=Array[Random[]&, {n, 2}], 
        ShowPolygon->False]]][[{1,-1}]] ]]; 
SeedRandom[s];RandomPolygon[pts=Array[Random[]&, {n, 2}],opts]];

RandomPolygon[pts_List, opts___] := Module[
  {ans, ans1,smQ, m, smc, scQ, n = Length[pts]}, 
    {trQ, smc, scQ,ist, showQ,ptQ,met,smoothn, psty,bg,
pcol, borsty, chQ} =
    {ShowSteps, SmoothingParameter, 
     SmoothCurve, InitialStep, ShowPolygon, ShowPoints, Method, SmoothingIterations, PointStyle,
Background, PolygonColor, BorderStyle, CorrectnessCheckWhenSmoothing} /.
        {opts} /. Options[RandomPolygon];  

pts1 = pts[[RandomPermutation[n]]];

If[trQ && met =!= SteadyGrowth, Message[RandomPolygon::notrac];    
     Return[HoldForm[RandomPolygon[pts1]]]];
smQ = (smoothn > 0);
If[!ListQ[psty], psty = {psty}];
If[!ListQ[borsty] && borsty =!= None, borsty = {borsty}];


rawpoly = Which[met===SteadyGrowth, RandomPolygonRaw,
        met === PermuteAndReject, RandomPolygonPermuteAndReject][pts1];
  ans = If[smQ, Nest[Smooth[#, smc]&, rawpoly, smoothn],rawpoly]; 
  If[smQ && chQ && !PolygonQ[ans], Message[RandomPolygon::badwarn]; 
          Return[RandomPolygon[n, opts]]]; 
  If[scQ, Interpolation /@ 
    {Transpose[{Range[0, 1, 1/Length[ans]], First /@ extend[ans]}], 
     Transpose[{Range[0, 1, 1/Length[ans]], Last /@ extend[ans]}]},
  If[trQ && met === SteadyGrowth, Do[
     pp = Take[order,j-1]; np = order[[j]]; 
     ans1=Select[ans, #==np||MemberQ[pp, #]&];
     pos = Position[ans1, np][[1,1]];
     pol = Sort[pp, Position[ans,#1][[1,1]]<Position[ans,#2][[1,1]]&];
Show[Graphics[{
      {RGBColor[1,1,.5], Polygon[pol]},
      PointSize[0.02], RGBColor[0,0,1], Point/@ pp,
      {RGBColor[0,0,1],Polygon[ans1[[Mod[pos+{-1,0,1},Length[ans1],1]]]] }, Hue[0], Point /@ Take[order, {j,n}],
      {RGBColor[0,1,0],PointSize[.04],Point[np],
      Hue[0],
        PointSize[.02],Point[order[[j]]]},GrayLevel[0],Thickness[0.01],
      Line[extend[pp[[ConvexHull[pp]]]]]}],
  PlotRange -> {{-.10, 1.1}, {-.1, 1.1}}, AspectRatio->1,
  Frame -> True, FrameTicks -> False],
     {j, Max[3, ist], Length[order]}];
  Show[Graphics[{
     {RGBColor[1,1,.5], Polygon[ans],
     GrayLevel[0], Line[extend[ans[[ConvexHull[ans]]]]],
     PointSize[0.02], RGBColor[0,0,1], Point/@ ans}}],
       PlotRange->{{-0.1,1.1},{-0.1,1.1}},
       AspectRatio -> 1, Frame -> True, FrameTicks->False]
   ];

    If[showQ, Show[Graphics[{{pcol,Polygon[ans]},
     Append[If[borsty=!= None,borsty,{}], Line[extend[ans]]],
      If[ptQ, Prepend[Point /@ ans,psty],{}]}],
  FilterOptions[Graphics, opts],
  PlotRange -> All,
  AspectRatio -> Automatic,Frame -> True,FrameTicks ->None,
  DefaultColor -> RGBColor[0, 0, .2], Background -> bg]];
   ans]]; 

Smooth[data_, c_:0.25] := Flatten[Transpose[{
  (1/c - 1) data + RotateRight[data], 
  (1/c - 1) data + RotateLeft[data]} * c], 1] /; c < 0.5;

PolygonQ[x_] := Module[{ans = True, 
  edges = Partition[extend[x], 2, 1]}, 
  Do[If[IntersectInteriorQ[edges[[i]], edges[[j]]], 
    ans = False; Break[]], {i, 3, Length[x]}, {j, i - 2}]; 
  ans];

(* the following routines come from my computational geometry package, so there is some duplication. *)

SACom = Compile[{x1, y1, x2, y2, x3, y3}, 
    0.5*(-x2*y1 + x3*y1 + x1*y2 - x3*y2 - x1*y3 + x2*y3)]; 

SignedArea[triangle_] := SACom @@ Flatten[triangle];
Area[p_] := Abs[SignedArea[p]];

ConvexVertex[p_] := First[First[Position[p, First[Sort[p]]]]];

PolygonDiagonal[poly_] := Module[{v = ConvexVertex[poly], invertices, 
    u, w}, u = v - 1 /. 0 -> Length[poly]; 
    w = v + 1 /. 1 + Length[poly] -> 1;
	invertices = Select[Delete[poly, {{u}, {v}, {w}}], 
      InsideTriangleQ[poly[[{u, v, w}]], #1] &]; 
    Sort[If[invertices == {}, {u, w}, 
      Identity[{Position[poly, Last[Sort[invertices, 
            Area[{poly[[u]], #1, poly[[w]]}] <= Area[{poly[[u]], #2, 
                poly[[w]]}] & ]]][[1,1]], v}]]]];

Triangulate[poly_List] := If[Orientation[poly] == 1, 
   Triangulate[poly, Range[Length[poly]]], 
   Reverse /@ (Length[poly] + 1 - Triangulate[Reverse[poly], 
      Range[Length[poly]]])];

Triangulate[poly_, inds_] := 
  Module[{d = PolygonDiagonal[poly[[inds]]]}, 
   Join[Triangulate[poly, Take[inds, d]], Triangulate[poly, 
     Drop[inds, d + {1, -1}]]]];

Triangulate[_, {i1_, i2_, i3_}] := {{i1, i2, i3}};

Orientation[triangle_] := Sign[SignedArea[triangle]] /; 
   Length[triangle] == 3;

Orientation[poly_] := Orientation[poly[[
ConvexVertex[poly] + {-1, 0, 1} /. 
      {0 -> Length[poly], Length[poly] + 1 -> 1}]]] /; Length[poly] > 3;

End[]; 
EndPackage[];