loop/src/Geometry/Polygon.hs

{-# LANGUAGE BangPatterns #-}

module Geometry.Polygon where

import Geometry.Intersect
import qualified Control.Foldl as L
import Data.Maybe
import Geometry.Data
import Geometry.LHS
import Geometry.Vector
import ListHelp
import Linear

-- | Draw an anticlockwise rectangle based on maximal N S W E values.
rectNSWE :: Float -> Float -> Float -> Float -> [Point2]
rectNSWE !n !s !w !e = [V2 w n, V2 w s, V2 e s, V2 e n]

-- | Draw an anticlockwise rectangle around the origin with given height and width
rectWH :: Float -> Float -> [Point2]
rectWH w h = rectNSWE h (- h) (- w) w

trapTBH :: Float -> Float -> Float -> [Point2]
trapTBH t b h = 
    [ V2 (-b) (-h)
    , V2 b (-h)
    , V2 t h
    , V2 (-t) h
    ]

isotriBWH :: Point2 -> Float -> Float -> [Point2]
isotriBWH (V2 x y) w h = [V2 (x - w) y, V2 (x + w) y, V2 x (y + h)]

-- trapezion
trapezionBWHW :: Point2 -> Float -> Float -> Float -> [Point2]
trapezionBWHW (V2 x y) w1 h w2 =
    [ V2 (x - w1) y
    , V2 (x + w1) y
    , V2 (x + w2) (y + h)
    , V2 (x - w2) (y + h)
    ]

rectXH :: Float -> Float -> [Point2]
rectXH x h = rectNSWE h (- h) 0 x

rectXY :: Float -> Float -> [Point2]
rectXY x y = rectNSWE y 0 0 x

rectVV :: Point2 -> Point2 -> [Point2]
rectVV (V2 x y) (V2 a b) = rectNSWE n s w e
  where
    (e,w) | x > a = (x,a)
        | otherwise = (a,x)
    (n,s) | y > b = (y,b)
        | otherwise = (b,y)

square :: Float -> [Point2]
square n = rectWH n n

polyOrthDist :: Int -> Float -> [Point2]
polyOrthDist n x = mapMaybe
            (\(ra, rb) -> intersectLineLine (rotateV ra bl) (rotateV ra br) (rotateV rb bl) (rotateV rb br))
            $ loopPairs rots
  where
    rot = 2 * pi / fromIntegral n
    rots = map ((rot *) . fromIntegral) [0 .. n -1]
    bl = V2 x x
    br = V2 (- x) x

polyCornerDist :: Int -> Float -> [Point2]
polyCornerDist n x = map f rots
  where
    rot = 2 * pi / fromIntegral n
    rots = map ((rot *) . fromIntegral) [0 .. n -1]
    f a = rotateV a (V2 x 0)

mirrorXAxis :: [Point2] -> [Point2]
mirrorXAxis ps = orderPolygon $ ps ++ mapMaybe f ps
  where
    f (V2 _ 0) = Nothing
    f (V2 x y) = Just $ V2 x (- y)

{- | Test whether a point is in a polygon or on the polygon border.
 Supposes the points in the
 polygon are listed in anticlockwise order.
-}
pointInOrOnPolygon :: Point2 -> [Point2] -> Bool
pointInOrOnPolygon !p (x : xs) = all (\l -> not (uncurry isRHS l p)) $ zip (x : xs) (xs ++ [x])
pointInOrOnPolygon _ _ = undefined

{- | Test whether a point is strictly inside a polygon.
 Supposes the points in the polygon are listed in anticlockwise order.
Requires that the polygon is convex.
-}
pointInPoly :: Point2 -> [Point2] -> Bool
pointInPoly !p (x : xs) = all (\l -> uncurry isLHS l p) $ zip (x : xs) (xs ++ [x])
pointInPoly _ [] = False

inSimplePoly :: Point2 -> [Point2] -> Bool
inSimplePoly p (x:xs) = foldl' (flip f) True $ zip (x:xs) (xs ++ [x])
  where
    f (a,b) = case intersectSegRay a b p (p + V2 1 0) of
        Nothing -> id
        Just {} -> not
inSimplePoly _ [] = False

---- implement Dan Sunday point in polygon algorithm?
--wnPointPoly :: Point2 -> Point2 -> Point2 -> Int
--wnPointPoly p x y = 0

circInPolygon :: Point2 -> Float -> [Point2] -> Bool
circInPolygon !p !r (x : xs) = all f $ zip (x : xs) (xs ++ [x])
  where
    f l = uncurry isLHS l (p - r *.* vNormal (normalizeV (uncurry (-.-) l)))
circInPolygon _ _ [] = False

orderPolygonAround ::
    -- | point to order around
    Point2 ->
    [Point2] ->
    [Point2]
orderPolygonAround _ [] = []
orderPolygonAround cen ps = sortOn (\p -> argV (p -.- cen)) ps

orderAroundFirstReverse :: [Point2] -> [Point2]
orderAroundFirstReverse [] = []
orderAroundFirstReverse (a : as) = a : reverse (orderPolygonAround a as)

orderAroundFirst :: [Point2] -> [Point2]
orderAroundFirst [] = []
orderAroundFirst (a : as) = a : orderPolygonAround a as

-- | Reorder points to be anticlockwise around their center.
orderPolygon :: [Point2] -> [Point2]
orderPolygon [] = []
--orderPolygon ps = orderPolygonAround (1/ fromIntegral (length ps) *.* foldr1 (+.+) ps) ps
orderPolygon ps = orderPolygonAround (centroid ps) ps

{- | Adds a point to a convex polygon.
 If the point is inside, returns the original.
 Points ordered anticlockwise, input not checked.
-}
addPointPolygon :: Point2 -> [Point2] -> [Point2]
addPointPolygon p ps
    | pointInOrOnPolygon p ps = ps
    | otherwise = orderPolygon $ p : ps

{- | Creates the convex hull of a set of points.
 Need to verify whether or not this is ordered
-}
convexHull :: [Point2] -> [Point2]
convexHull (x : y : z : xs) = grahamScan $ orderAroundFirst $ sortOn (\(V2 a b) -> (b, a)) (x : y : z : xs)
convexHull _ = error "Tried to create the convex hull of two or fewer points"

---- assumes the points go "anticlockwise" around a non-self intersecting shape
--convexPartition :: [Point2] -> [[Point2]]
--convexPartition (x:y:z:[]) = [[x,y,z]]
--convexPartition (x:y:z:xs)
--    | isLHS x y z = [x,y,z] : convexPartition (x:z:xs)
--    | otherwise = convexPartition (y:z:xs <> [x])
--convexPartition _ = error "unexpected shape for convexPartition"

triangulateEarClip :: [Point2] -> [[Point2]]
triangulateEarClip [x,y,z] = [[x,y,z]]
triangulateEarClip (x:y:z:xs)
    | isLHS x y z && not (any (`pointInPoly` [x,y,z]) xs) = [x,y,z] : triangulateEarClip (x:z:xs)
    | otherwise = triangulateEarClip (y:z:xs ++ [x])
triangulateEarClip _ = error "triangulateEarClip: non-simple polygon input?"

{- | Creates the convex hull of a set of points.
 assumes no repetition of points: try nubbing!
-}
convexHullSafe :: [Point2] -> [Point2]
--convexHullSafe (x:y:z:xs) = grahamScan $ orderAroundFirst $ sortOn (\(V2 a b) -> (b,a)) (x:y:z:xs)
convexHullSafe (x : y : z : xs) = grahamScan $ orderAroundFirst $ sortOn (\(V2 a b) -> (b, a)) (x : y : z : xs)
convexHullSafe _ = []

grahamScan :: [Point2] -> [Point2]
grahamScan = foldr push []
  where
    push p stack = grahamEliminate (p : stack)

{- | Remove second element if top three elements are not counterclockwise.
 Repeat if necessary. See
 https://codereview.stackexchange.com/questions/206019/graham-scan-algorithm-in-haskell
-}
grahamEliminate :: [Point2] -> [Point2]
grahamEliminate (x : y : z : xs)
    | not $ isLHS x y z = grahamEliminate (x : z : xs)
grahamEliminate xs = xs

-- not sure what definition of centroid is applicable here
centroid :: (Num (f a),Functor f, Fractional a,Foldable t) => t (f a) -> f a
centroid = L.fold $ (^/) <$> L.Fold (+) 0 id <*> L.genericLength

centroidNum :: (Fractional a, Foldable t) => t a -> a
centroidNum = L.fold $ (/) <$> L.Fold (+) 0 id <*> L.genericLength

shrinkPolyOnEdges :: Float -> [Point2] -> [Point2]
shrinkPolyOnEdges x (p : q : ps) = map (shrinkVert x) . slideWindow 3 $ (p : q : ps) ++ [p, q]
shrinkPolyOnEdges _ _ = error "too few vertices in polygon"

shrinkVert :: Float -> [Point2] -> Point2
shrinkVert d [x, y, z] = x +.+ (d *.* normalizeV (x -.- y)) +.+ (d *.* normalizeV (z -.- y))
shrinkVert _ _ = error "wrong number of vertices"

-- divide a polygon into two along a line
-- assumes the line intersects the polygon exactly twice, nocolinearity
-- this may duplicate points
-- this should be tested: there are many possible points of failure...
cutPoly :: Point2 -> Point2 -> [Point2] -> ([Point2],[Point2])
cutPoly a b (p:ps)
    | not (isLHS a b p) && not (isRHS a b p) = cutPolyL a b ps p p [p] [p]
    | isLHS a b p = cutPolyL a b ps p p [p] []
    | otherwise = cutPolyR a b ps p p [] [p]
cutPoly _ _ [] = error "cutPoly empty poly"

cutPolyL :: Point2 -> Point2 -> [Point2] -> Point2 -> Point2 -> [Point2] -> [Point2] -> ([Point2], [Point2])
cutPolyL a b (p:ps) e x ls rs
    | isLHS a b p = cutPolyL a b ps e p (p\:ls) rs
    | otherwise = case intersectLineLine x p a b of
        Nothing -> error "cutPolyL nonintersecting lines"
        Just p' -> cutPolyR a b ps e p (p'\:ls) (p\:p'\:rs)
cutPolyL a b [] e x ls rs = case intersectSegLine x e a b of
    Nothing -> (reverse ls,reverse rs)
    Just p -> (reverse (p\:ls), reverse (p\:rs))

cutPolyR :: Point2 -> Point2 -> [Point2] -> Point2 -> Point2 -> [Point2] -> [Point2] -> ([Point2], [Point2])
cutPolyR a b (p:ps) e x ls rs
    | isRHS a b p = cutPolyR a b ps e p ls (p\:rs)
    | otherwise = case intersectLineLine x p a b of
        Nothing -> error "cutPolyR nonintersecting lines"
        Just p' -> cutPolyL a b ps e p (p\:p'\:ls) (p'\:rs)
cutPolyR a b [] e x ls rs = case intersectSegLine x e a b of
    Nothing -> (reverse ls,reverse rs)
    Just p -> (reverse (p\:ls), reverse (p\:rs))

polyInPoly :: Point2 -> [Point2] -> [Point2] -> [Point2]
polyInPoly p ps = mapMaybe f
  where
    f q = intersectRayPoly p q ps

infixr 5 \:
(\:) :: Eq a => a -> [a] -> [a]
(\:) x (y:ys) 
    | x /= y = x:y:ys
    | otherwise = y:ys 
(\:) x [] = [x]