Voronoi Diagrams

Voronoi.hs

{-# 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) ]