Visibility.hs

Visibility.hs

module Visibility where

import Data.List (sortBy, minimumBy)
import Data.Ord
import Data.Maybe (mapMaybe)
import Diagrams (Diagram, P2, (#), Located)
import qualified Diagrams as Diag
import Diagrams.Backend.Rasterific (Rasterific)
import qualified Diagrams.Backend.Rasterific as Rast
import Data.Colour (Colour)
import qualified Data.Colour.Names as Color
--import Debug.Trace (traceShow)

data Point = Point Double Double deriving (Eq, Ord, Read, Show)
data Segment = Segment Point Point deriving (Eq, Ord, Read, Show)
data SimplePolygon = Simple [Point] deriving (Eq, Ord, Read, Show)
data Polygon = Polygon SimplePolygon [SimplePolygon] deriving (Eq, Ord, Read, Show)

type Vector = Point

instance Num Point where
    (Point x1 y1) + (Point x2 y2) = Point (x1 + x2) (y1 + y2)
    (Point x1 y1) * (Point x2 y2) = Point (x1 * x2) (y1 * y2)
    abs = error "abs not defined for points"
    signum = error "signum not defined for points"
    fromInteger = error "fromInteger not defined for points"
    negate (Point x y) = Point (-x) (-y)

vecTo :: Point -> Point -> Vector
vecTo p1 p2 = p2 - p1

direction :: Vector -> Double
direction (Point x y) = atan2 y x

sqrDist :: Point -> Point -> Double
sqrDist (Point x1 y1) (Point x2 y2) = (x2 - x1) ** 2 + (y2 - y1) ** 2

cross :: Vector -> Vector -> Double
cross (Point x1 y1) (Point x2 y2) = x1 * y2 - y1 * x2

scale :: Double -> Vector -> Vector
scale a (Point x y) = Point (a * x) (a * y)

rotate :: Vector -> Double -> Vector
rotate (Point x y) a = Point (cos a * x - sin a * y) (sin a * x + cos a * y)

jitter :: Vector -> [Vector]
jitter v = [rotate v (-0.0001), v, rotate v 0.0001]

pairwise :: [a] -> [(a, a)]
pairwise [] = []
pairwise xs = zip xs (tail xs)

segments :: SimplePolygon -> [Segment]
segments (Simple pts) = map (uncurry Segment) pairs
    where pairs = pairwise (pts ++ [head pts])

visibilityPolygon :: Point -> Polygon -> SimplePolygon
visibilityPolygon start poly = Simple cont
    where verts        = vertices poly
          cont         = sortBy (comparing vertDir) (concatMap cast verts)
          cast vert    = map (\v -> rayCast start v poly) (jitter $ start `vecTo` vert)
          vertDir vert = direction $ start `vecTo` vert

intersect :: Point -> Vector -> Segment -> Maybe Point
intersect p r (Segment q q')
    | t <= 0                = Nothing
    | u >= 0 && u <= 1 = Just (q + u `scale` s)
    | otherwise             = Nothing
    where s = q `vecTo` q'
          t = ((q - p) `cross` s) / (r `cross` s)
          u = ((p - q) `cross` r) / (s `cross` r)

rayCast :: Point -> Vector -> Polygon -> Point
rayCast start dir (Polygon outer holes) = minimumBy (comparing (sqrDist start)) inters
    where segs   = segments outer ++ concatMap segments holes
          inters = mapMaybe (start `intersect` dir) segs

vertices :: Polygon -> [Point]
vertices (Polygon (Simple outer) holes) = outer ++ concatMap (\(Simple inner) -> inner) holes




-- Diagrams -------------------------------
pointToP2 :: Point -> P2 Double
pointToP2 (Point x y) = Diag.p2 (x, y)

locate :: P2 Double -> Located (P2 Double)
locate p = p `Diag.at` Diag.p2 (0, 0)

drawPolygon :: Polygon -> Diagram Rasterific
drawPolygon (Polygon outer holes) = mconcat $ reverse $ drawSimplePolygon outer Color.black
                                                      : map (`drawSimplePolygon` Color.white) holes

drawSimplePolygon :: SimplePolygon -> Colour Double -> Diagram Rasterific
drawSimplePolygon (Simple pts) col =
    locPath # Diag.closeLine
            # Diag.strokeLoop
            # Diag.fc col
            # Diag.translate first
    where locPath = Diag.fromVertices (map pointToP2 pts)
          first   = head pts # \(Point x y) -> Diag.r2 (x, y)


star :: Int -> Double -> SimplePolygon
star n rad = Simple $ concat $ zipWith (\i o -> [i, o]) inner outer
    where angleDelta = 2 * pi / fromIntegral n
          polar r phi = Point (r * cos phi) (r * sin phi)
          inner = map (\x -> polar rad (angleDelta * fromIntegral x)) [1..n]
          outer = map (\x -> polar (rad * 2) (angleDelta * fromIntegral x + angleDelta / 2)) [1..n]

polyStar :: Polygon
polyStar = Polygon (star 5 10) [star 5 5]

visPoly :: Double -> Double -> SimplePolygon
visPoly x y = visibilityPolygon (Point x y) polyStar

testDiag :: Double -> Double -> Diagram Rasterific
testDiag x y = Diag.circle 1 # Diag.translate (Diag.r2 (x, y)) # Diag.fc Color.red `Diag.atop` visDiag
    where visDiag = drawPolygon polyStar `Diag.beneath` drawSimplePolygon (visPoly x y) Color.yellow

animation :: [Diagram Rasterific]
animation = map (uncurry testDiag . polar) [0, pi / 40..2 * pi]
    where polar a = (rad a * cos a, rad a * sin a)
          rad   a = (sin (a * 5 - pi / 2) + 1) * 3 + 7

main :: IO ()
main = mapM_ (uncurry render) indexedFrames
    where indexedFrames = zip animation [(1 :: Int)..]
          render diag i = Rast.renderRasterific ("out/frame" ++ show i ++ ".png") (Diag.mkHeight 400) diag