Download or view ConvexHull.frink in plain text format
/** This program calculates the convex hull of a list of points.
It uses Graham's scan algorithm which is approximately O(n log n) with
respect to the number of points.
You wil generally call this with ConvexHull.hull[points] where points is
an array of [x,y] points.
These points can be obtained or passed to Frink's graphics polygon class:
https://frinklang.org/#PolygonMethods
To graph your solution, call:
ConvexHull.plotHull[points].show[]
where points is an array of [x,y] points.
This somewhat follows the discussion at:
http://www.geeksforgeeks.org/convex-hull-set-2-graham-scan/
It should be noted that the sample implementation at that site is broken
in several ways.
*/
class ConvexHull
{
// This functions calculates the convex hull of an array of points specified
// as [x,y] coordinates. It returns the result as an array of [x,y] points.
// These can be obtained or passed to Frink's graphics polygon class:
// https://frinklang.org/#PolygonMethods
class hull[points] :=
{
points = toArray[points] // Points can come in as a set or enumeration
size = length[points]
// Find the smallest y value, smallest right if there's a tie
ymin = points@0@1
minIndex = 0
for i = 1 to size-1
{
[x, y] = points@i
if ((y < ymin) or (ymin == y and x < points@minIndex@0))
{
ymin = y
minIndex = i
}
}
// Place the smallest y value at index 0 by swapping
temp = points@0
points@0 = points@minIndex
points@minIndex = temp
// Sort items 1 through size-1 with respect to point 0.
// A point p1 comes betwen p2 in sorted output if p2 has a larger polar
// angle (in counterclockwise direction) than p1.
// TODO: Implement sort with sub-region begin and end specified
result = sort[rest[points], getFunction["compare", 3], points@0]
result.pushFirst[points@0]
points = result
// If two or more points make the same angle with p0, remove all but
// the farthest. The sort above made the farthest last in the list.
m = 1
for i=1 to size-1
{
// Keep removing element i while i and i+1 are colinear with p0
while i < size-1 and orientation[points@0, points@i, points@(i+1)] == 0
i = i + 1
points@m = points@i
m = m + 1 // Update size of modified array
}
// If modified array has fewer than 3 points, no convex hull is possible
if m < 3
return undef
// Create an empty stack and push first three points onto it.
// TODO: Should first[array] return an array? It currently returns an
// EnumeratingExpression which led to confusion
s = toArray[first[points, 3]]
// Process remaining n-3 points
for i = 2 to m-1
{
// Keep removing top while the angle formed by points next-to-top,
// top, and points@i makes a non-left turn. The first clause is
// apparently necessary to fix a bug in the original
// TODO: Factor out multiple calls to length[s]
while length[s] >= 2 and orientation[s@(length[s]-2), s@(length[s]-1), points@i] != 2
s.pop[]
s.push[points@i]
}
return s
}
// This is a comparison function that compares points with respect to the
// first point p0. (Which is passed as an extra "data" argument to the sort
// function)
class compare[p1, p2, p0] :=
{
o = orientation[p0, p1, p2]
// If points are colinear, put nearest point first
if o == 0
return -(distanceSquared[p0, p2] <=> distanceSquared[p0, p1])
return o==2 ? -1 : 1
}
/** Returns the squared distance between two points. */
class distanceSquared[p1, p2] :=
{
[x1,y1] = p1
[x2,y2] = p2
// TODO: This doesn't make sense with disparate units
return (x1-x2)^2 + (y1-y2)^2
}
/** Finds the orientation of the ordered triplet of points [p,q,r].
Returns:
0 if p, q, and r are colinear
1 if clockwise
2 if counterclockwise
*/
class orientation[p, q, r] :=
{
[px, py] = p
[qx, qy] = q
[rx, ry] = r
v = (qy - py) * (rx - qx) - (qx - px) * (ry - qy)
units = px py
if v == 0 units // This makes units come out right
return 0 // Colinear
else
return v > 0 units ? 1 : 2 // Clockwise or counterclockwise
}
/** This is a function to calculate and plot a hull of points. It returns
a graphics object which you can display or print with its show[] or
print[] methods. */
class plotHull[points] :=
{
g = new graphics
for [x,y] = points
g.fillEllipseCenter[x,y,.1, .1]
hull = hull[points]
ph = new polygon[hull]
g.color[0,0,1]
g.add[ph]
g.color[1,0,0]
// Draw the actually used hull points in red
for [x,y] = hull
g.fillEllipseCenter[x,y,.2, .2]
return g
}
}
/*
// Testing stuff.
p = [ [0,3], [1,1], [2,2], [4,4], [0,0], [1,2], [3,1], [3,3] ]
ConvexHull.plotHull[p].show[]
println[ConvexHull.hull[p]] // Sample from website above (in reverse order)
p = [ [0,3], [1,1], [2,2], [4,4], [0,0], [1,2], [3,1], [3,3], [2,4] ]
ConvexHull.plotHull[p].show[]
println[ConvexHull.hull[p]]
println["Colinear3"]
p = [ [0,0], [1,1], [2,2]]
ConvexHull.plotHull[p].show[]
println[ConvexHull.hull[p]]
println["Colinear4"]
p = [ [0,0], [1,1], [2,2], [3,3]]
ConvexHull.plotHull[p].show[]
println[ConvexHull.hull[p]]
mean = 5
sd = 1
d = new array
for p = 1 to 2000
d.push[ [randomGaussian[mean, sd], randomGaussian[mean, sd]] ]
g = ConvexHull.plotHull[d]
g.show[]
*/
// Testing with units of measure. Not sure if this can be made to make sense
//p2 = [ [0,3s], [1,1s], [2,2s], [4,4s], [0,0s], [1,2s], [3,1s], [3,3s], [2,4s] ]
//println[ConvexHull.hull[p2]]
Download or view ConvexHull.frink in plain text format
This is a program written in the programming language Frink.
For more information, view the Frink
Documentation or see More Sample Frink Programs.
Alan Eliasen, eliasen@mindspring.com