Created
November 12, 2018 08:32
-
-
Save fmthoma/20568da67230ab5be3a07bf4b4b37f5f to your computer and use it in GitHub Desktop.
Voronoi Diagrams
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{-# LANGUAGE RecordWildCards #-} | |
import Control.Applicative (liftA2) | |
import Control.Monad (replicateM) | |
import Data.Foldable (for_) | |
import Data.List (find, foldl') | |
import Data.Maybe (mapMaybe) | |
import Prelude hiding ((**)) | |
import System.Environment (getArgs) | |
import System.Random (randomRIO) | |
import Graphics.Rendering.Cairo | |
data Point = Point Rational Rational | |
deriving (Eq, Show) | |
data Vector = Vector Rational Rational | |
deriving (Eq, Show) | |
newtype Polygon | |
= Polygon [Point] | |
deriving (Eq, Show) | |
data LineSegment = LineSegment Point Point | |
deriving (Eq, Show) | |
data Line = Line Point Vector | |
deriving (Eq, Show) | |
data Intersection | |
= Intersection Point | |
| VirtualIntersectionL Point | |
| VirtualIntersectionR Point | |
| VirtualIntersection Point | |
| Parallel | |
| Identical | |
deriving (Eq, Show) | |
intersectLS :: LineSegment -> LineSegment -> Intersection | |
intersectLS (LineSegment p1@(Point p1x p1y) q1@(Point q1x q1y)) (LineSegment p2 q2) | |
| p1 == p2 && q1 == q2 = Identical | |
| det d1 d2 == 0 = Parallel | |
| isVirtualL && isVirtualR = VirtualIntersection intersection | |
| isVirtualL = VirtualIntersectionL intersection | |
| isVirtualR = VirtualIntersectionR intersection | |
| otherwise = Intersection intersection | |
where | |
d1 = delta p1 q1 | |
d2 = delta p2 q2 | |
u = det (delta p1 p2) d2 / det d1 d2 | |
v = det (delta p1 p2) d1 / det d1 d2 | |
isVirtualL = u < 0 || u > 1 | |
isVirtualR = v < 0 || v > 1 | |
intersection = Point (p1x + u * (q1x - p1x)) (p1y + u * (q1y - p1y)) | |
intersectL :: LineSegment -> Line -> Intersection | |
intersectL ls (Line p d) = case intersectLS ls (LineSegment p (p `vplus` d)) of | |
VirtualIntersection intersection -> VirtualIntersectionL intersection | |
VirtualIntersectionR intersection -> Intersection intersection | |
other -> other | |
det :: Vector -> Vector -> Rational | |
det (Vector x1 y1) (Vector x2 y2) = x1*y2 - x2*y1 | |
delta :: Point -> Point -> Vector | |
delta (Point x1 y1) (Point x2 y2) = Vector (x2-x1) (y2-y1) | |
vplus :: Point -> Vector -> Point | |
vplus (Point x y) (Vector dx dy) = Point (x + dx) (y + dy) | |
angle :: Vector -> Vector -> Double | |
angle (Vector x1 y1) (Vector x2 y2) = atan2 (fromRational (x2 - x1)) (fromRational (y2 - y1)) | |
angle' :: Point -> Point -> Point -> Double | |
angle' pivot p1 p2 = angle (delta pivot p1) (delta pivot p2) | |
data VoronoiFace = VF { center :: Point, face :: Polygon } | |
deriving (Eq, Show) | |
data Voronoi = Voronoi { bounds :: Polygon, faces :: [VoronoiFace] } | |
deriving (Eq, Show) | |
addPoint :: Voronoi -> Point -> Voronoi | |
addPoint Voronoi{..} p = Voronoi bounds (newFace : faces') | |
where | |
newFace :: VoronoiFace | |
newFace = foldl' (\nf f -> updateFace nf (center f)) (VF p bounds) faces | |
faces' :: [VoronoiFace] | |
faces' = fmap (\f -> updateFace f (center newFace)) faces | |
updateFace :: VoronoiFace -> Point -> VoronoiFace | |
updateFace f p = VF (center f) (clipFace (face f) (Line q d)) | |
where | |
q = let Point px py = p | |
Point fx fy = center f | |
in Point ((px + fx) / 2) ((py + fy) / 2) | |
d = let Vector dx dy = delta (center f) p | |
in Vector (-dy) dx | |
-- | Keep everything that's *left* of the line | |
clipFace :: Polygon -> Line -> Polygon | |
clipFace (Polygon ps) l@(Line p d) = Polygon $ do | |
(p1, p2) <- zip ps (tail (cycle ps)) | |
case intersectL (LineSegment p1 p2) l of | |
Intersection intersection -> case det (delta p1 p2) d of | |
discriminant | discriminant > 0 -> [p1, intersection] -- line leaving | |
| otherwise -> [intersection] -- line returning | |
_otherwise -> case det d (delta p p1) of | |
discriminant | discriminant < 0 -> [] -- p1 right of line | |
| otherwise -> [p1] | |
drawFace :: VoronoiFace -> RGB -> Render () | |
drawFace VF{..} = drawPoly face | |
drawPoly :: Polygon -> RGB -> Render () | |
drawPoly (Polygon []) _ = pure () | |
drawPoly (Polygon (p:ps)) RGB{..} = do | |
newPath | |
let Point x y = p in moveTo (fromRational x) (fromRational y) | |
for_ ps $ \p -> let Point x y = p in lineTo (fromRational x) (fromRational y) | |
closePath | |
setSourceRGB r g b | |
fillPreserve | |
setSourceRGB (r + 0.1) (g + 0.1) (b + 0.1) | |
stroke | |
main :: IO () | |
main = do | |
[count, file] <- getArgs | |
let w, h :: Num a => a | |
w = 2560 | |
h = 1440 | |
initialPolygon = Polygon [ Point 1 0, Point w 0, Point w h, Point 0 h ] | |
points <- replicateM (read count) (randomPoint w h) | |
withSVGSurface file w h $ \surface -> do | |
let allFaces = faces $ foldl' addPoint (Voronoi initialPolygon []) points | |
backgroundFaces = filter (\face -> not $ any (pointInPolygon (center face) . snd) haskellLogo) allFaces | |
pictureFaces = mapMaybe (\face -> fmap (\(color, _) -> (color, face)) $ find (pointInPolygon (center face) . snd) haskellLogo) allFaces | |
for_ backgroundFaces $ \face -> renderWith surface (drawFace face (RGB 0.1 0.1 0.1)) | |
for_ pictureFaces $ \(color, face) -> renderWith surface (drawFace face color) | |
where | |
randomPoint :: Int -> Int -> IO Point | |
randomPoint width height = liftA2 Point (randomCoordinate width) (randomCoordinate height) | |
randomCoordinate :: Int -> IO Rational | |
randomCoordinate mx = fmap fromIntegral $ randomRIO (0, mx) | |
data RGB = RGB { r :: Double, g :: Double, b :: Double } | |
-- Ray-casting algorithm. | |
pointInPolygon :: Point -> Polygon -> Bool | |
pointInPolygon p (Polygon ps) = odd (length intersections) | |
where | |
-- The test ray comes from outside the polygon from the left, and ends at | |
-- the point to be tested. | |
testRay = LineSegment (Point (leftmostPolyX - 1) pointY) p | |
where | |
leftmostPolyX = minimum (edges >>= \(LineSegment (Point x1 _) (Point x2 _)) -> [x1,x2]) | |
Point _ pointY = p | |
intersections = flip filter edges $ \edge -> | |
case intersectLS testRay edge of | |
Intersection _ -> True | |
_other -> False | |
edges = zipWith LineSegment ps (tail (cycle ps)) | |
haskellLogo :: [(RGB, Polygon)] | |
haskellLogo = zip colors $ fmap translateAndResize [left, lambda, upper, lower] | |
where | |
left = Polygon [Point 0 340.15625, Point 113.386719 170.078125, Point 0 0, Point 85.039062 0, Point 198.425781 170.078125, Point 85.039062 340.15625, Point 0 340.15625] | |
lambda = Polygon [Point 113.386719 340.15625, Point 226.773438 170.078125, Point 113.386719 0, Point 198.425781 0, Point 425.195312 340.15625, Point 340.15625 340.15625, Point 269.292969 233.859375, Point 198.425781 340.15625, Point 113.386719 340.15625] | |
upper = Polygon [Point 387.402344 240.945312, Point 349.609375 184.253906, Point 481.890625 184.25, Point 481.890625 240.945312, Point 387.402344 240.945312] | |
lower = Polygon [Point 330.710938 155.90625, Point 292.914062 99.214844, Point 481.890625 99.210938, Point 481.890625 155.90625, Point 330.710938 155.90625] | |
translateAndResize (Polygon ps) = Polygon (fmap (\(Point x y) -> Point (3*x + 600) (3*y + 250)) ps) | |
colors = | |
[ RGB (0x45/255) (0x3a/255) (0x62/255) | |
, RGB (0x5e/255) (0x50/255) (0x86/255) | |
, RGB (0x8f/255) (0x4e/255) (0x8b/255) | |
, RGB (0x8f/255) (0x4e/255) (0x8b/255) ] |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment