Skip to content

Instantly share code, notes, and snippets.

@masterdezign
Last active January 21, 2022 00:28
Show Gist options
  • Star 4 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save masterdezign/ea46c5a3dc46d9fa873eb3030a54ad8d to your computer and use it in GitHub Desktop.
Save masterdezign/ea46c5a3dc46d9fa873eb3030a54ad8d to your computer and use it in GitHub Desktop.
Convex hull algorithm in Haskell
-- Worth reading http://www.geeksforgeeks.org/convex-hull-set-2-graham-scan/
import Text.Printf
import Data.List
type Point = (Double, Double)
-- Euclidean distance
dist :: (Double, Double) -> (Double, Double) -> Double
dist (x1, y1) (x2, y2) = sqrt (f x1 x2 + f y1 y2)
where f a b = (a - b)^2
-- Use Graham scan, https://en.wikipedia.org/wiki/Graham_scan
-- Three points are a counter-clockwise turn if ccw > 0, clockwise if
-- ccw < 0, and collinear if ccw = 0 because ccw is a determinant that
-- gives twice the signed area of the triangle formed by p1, p2 and p3.
ccv :: Num a => (a, a) -> (a, a) -> (a, a) -> a
ccv (x1, y1) (x2, y2) (x3, y3) = (x2 - x1) * (y3 - y1) - (y2 - y1) * (x3 - x1)
-- Find bottom-most point taking into account the x coordinate
findMinYPoint :: [Point] -> Point
findMinYPoint (p:points) = f p points
where
f prev [] = prev
f p0@(x0,y0) (p1@(x1,y1):points) | y1 < y0 = f p1 points
-- NB: select point with smaller x if y's are equal
| (y1 == y0) && (x1 < x0) = f p1 points
| otherwise = f p0 points
filterSameAngle :: [(Point, Double)] -> [Point]
filterSameAngle lst = map fst . filter snd $ r
where
r = zipWith (\(pa, ra) (pb, rb) -> (pa, ra /= rb)) lst ((drop 1 lst) ++ [head lst])
pretty :: [Point] -> String
pretty = unlines . map h
where h (x, y) = printf "\t%0.f %.0f" x y
build :: [Point] -> [Point] -> [Point]
build hull [] = hull
build hs@(h1:[]) ps@(p:points) = build (p:h1:[]) points
build hs@(h2:h1:hull) ps@(p:points) = hull'
where
rightTurn = ccv h1 h2 p < 0
collinear = ccv h2 h1 p == 0
hull' | rightTurn = build (h1:hull) ps -- Remove head and retry
| collinear = build (p:h1:hull) points -- Replace the head with p
| otherwise = build (p:hs) points -- Add p
convexHull :: [(Double, Double)] -> [(Double, Double)]
convexHull points = hull
where
-- 1) Find the bottom-most point
p0@(x0, y0) = findMinYPoint points
-- 2) Sort in increasing order of the angle they and the point P make with the x-axis
sorted' = let
o (a, ra) (b, rb) | ra > rb = GT
-- NB: Nearest points first, further points GT
| (ra == rb) && ((dist a p0) > (dist b p0)) = GT
| otherwise = LT
in sortBy o hullP
-- For performance, instead of atan, do not calculate the
-- angle. Use e.g. cos for ordering intead
f (x, y) = r'
where r = atan $ (y - y0) / (x - x0)
-- NB: Avoid negative angles values
r' | r < 0 = r + 2*pi
| otherwise = r
hullP = map (\p -> (p, f p)) $ delete p0 points
-- 3) NB: Remove points with same angles that are closer from p0
sorted = filterSameAngle sorted'
-- 4) Recursively build the convex hull
hull = build [p0] sorted
main :: IO ()
main = do
print $ convexHull [(1,1), (2,5), (3,3), (5,3), (3,2), (2,2)]
-- [(2.0,5.0),(5.0,3.0),(1.0,1.0)]
print $ convexHull [(1,1), (3,3), (3,1), (1, 3), (2, 2)]
-- [(1.0,3.0),(3.0,3.0),(3.0,1.0),(1.0,1.0)]
print $ convexHull [(3,2), (2,5), (4,5)]
-- [(2.0,5.0),(4.0,5.0),(3.0,2.0)]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment