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 was born 19966 days, 17 hours, 58 minutes ago.