--{-# LANGUAGE TupleSections #-}
module Geometry.Zone
    ( ddaExt
    , ddaExt'
    , ddaSq
    , ddaSqStream
    , ddaStream
    , ddaStreamX
    , ddaStreamY
    , xIntercepts
    , yIntercepts
    , divTo
    , zoneOfPoint
    ) where
import Geometry.Data

import Streaming
import qualified Streaming.Prelude as S
import Data.Foldable
import qualified Data.IntMap.Strict as IM
import qualified Data.IntSet as IS
--import Control.Monad

--foldl2'
--    :: (b -> a -> a -> b) 
--    -> b 
--    -> [a]
--    -> b
--foldl2' f s (t:ts) = fst $ foldl' g (s, t) ts
--  where
--    g (r,x) y = (f r x y,y)
--foldl2' _ s _ = s

--sortArguments 
--    :: Ord a 
--    => (a -> a -> b)
--    -> a -> a -> b
--sortArguments f x y
--    | x < y = f x y
--    | otherwise = f y x
--sortArgumentsReverse
--    :: Ord a 
--    => (a -> a -> [b])
--    -> a -> a -> [b]
--sortArgumentsReverse f x y
--    | x < y = f x y
--    | otherwise = reverse $ f y x
--
--intervalBounds
--    :: Float    -- ^ interval threshold
--    -> Float    -- ^ First endpoint
--    -> Float    -- ^ Second endpoint
--    -> [Float]
--intervalBounds = sortArgumentsReverse . f
--  where 
--    f r a b 
--        | x > b = [a]
--        | otherwise = (a : [x,x+r..b])
--      where
--        x = floorTo r a + r

--floorTo :: Float -> Float -> Float
--floorTo r x = r * (fromIntegral ((floor $ x / r) :: Int))

--ceilingTo :: Float -> Float -> Float
--ceilingTo r x = r * (fromIntegral ((ceiling $ x / r) :: Int))

divTo :: Float -> Float -> Int
{-# INLINE divTo #-}
divTo s = floor . (/s)

modTo :: Float -> Float -> Float
modTo s x = x - s * fromIntegral (divTo s x)

--remTo :: Float -> Float -> Float
--remTo s x = x - s * fromIntegral (quotTo s x)
--
--quotTo :: Float -> Float -> Int
--{-# INLINE quotTo #-}
--quotTo s = truncate . (/s)


--flipV :: Point2 -> Point2
--{-# INLINE flipV #-}
--flipV (V2 a b) = V2 b a

--applyInverted 
--    :: (Point2 -> Point2 -> [Point2])
--    -> Point2 -> Point2 -> [Point2]
--applyInverted f sp@(V2 sx sy) ep@(V2 ex ey)
--    | abs (sx-ex) > abs (sy-ey) = f sp ep
--    | otherwise = map flipV $ f (flipV sp) (flipV ep)

zoneOfPoint :: Float -> Point2 -> V2 Int
{-# INLINE zoneOfPoint #-}
zoneOfPoint s = fmap (divTo s)

--increasingInterval :: Int -> Int -> [Int]
--increasingInterval x y
--    | y > x = [x .. y]
--    | otherwise = [y .. x]

-- | Determines a "square" zone of points for a line
ddaSq :: Float -> V2 Float -> V2 Float -> IM.IntMap IS.IntSet
ddaSq s (V2 sx sy) (V2 ex ey) = IM.fromSet (const ys) xs
  where
    maxMin a b | a >= b = (a,b)
        | otherwise = (b,a)
    (maxx,minx) = maxMin (divTo s sx) (divTo s ex)
    (maxy,miny) = maxMin (divTo s sy) (divTo s ey)
    xs = IS.fromDistinctAscList [minx-1..maxx+1]
    ys = IS.fromDistinctAscList [miny-1..maxy+1]

--ddaSqStream :: Monad m => Float -> V2 Float -> V2 Float -> Stream (Of (V2 Int)) m ()
--ddaSqStream s (V2 sx sy) (V2 ex ey) = S.each [V2 x y | x <- [minx..maxx], y <- [miny..maxy]]
--  where
--    maxMin a b | a >= b = (a,b)
--        | otherwise = (b,a)
--    (maxx,minx) = maxMin (divTo s sx) (divTo s ex)
--    (maxy,miny) = maxMin (divTo s sy) (divTo s ey)

-- | Determines which horizontal and vertical lines on a grid are crossed by a
-- line. For each adds the x-y index of the square to the right or above the
-- crossed grid line. Also adds the index of the square containing the start
-- point.
ddaExt' :: Float -> V2 Float -> V2 Float -> IM.IntMap IS.IntSet
ddaExt' s sp@(V2 sx sy) ep@(V2 ex ey) 
    | x1 <= x2 = addsp . addys . IM.fromDistinctAscList $ zip [x1 .. x2] 
                    $ map (IS.singleton . divTo s) [x1y,x1y+ydx..]
    | otherwise = ddaExt' s ep sp
  where
    addsp im = let V2 x y = zoneOfPoint s sp
        in insertXY im (x,y)
    x1 = divTo s sx
    x2 = divTo s ex
    x1y = fx' sp ep $ s * fromIntegral x1
    ydx = s * ydx' sp ep
    addys m = add2s m ypairs 
    y1 = divTo s sy
    y2 = divTo s ey
    y1x = fy' sp ep $ s * fromIntegral y1
    y2x = fy' sp ep $ s * fromIntegral y2
    xdy = s * xdy' sp ep
    ypairs 
        | y1 <= y2  = zip (map (divTo s) [y1x,y1x+xdy..])
                         [y1 .. y2] 
        | otherwise = zip (map (divTo s) [y2x,y2x+xdy..])
                         [y2-1 .. y1-1] 
ddaSqStream :: Float -> Point2 -> Point2 -> Stream (Of (V2 Int)) Identity ()
ddaSqStream s sp ep = S.each [V2 x y | x <- makeInterval sx ex, y <- makeInterval sy ey]
  where
    V2 sx sy = zoneOfPoint s sp
    V2 ex ey = zoneOfPoint s ep

makeInterval :: Int -> Int -> [Int]
makeInterval x y
    | x < y = [x-1..y+1]
    | otherwise = [y-1..x+1]

ddaStream :: Float -> Point2 -> Point2 -> Stream (Of (V2 Int)) Identity ()
ddaStream s sp ep = S.map (zoneOfPoint s)
    $  S.yield sp
    <> xIntercepts s sp ep
    <> yIntercepts s sp ep

ddaStreamX :: Float -> Point2 -> Point2 -> Stream (Of (V2 Int)) Identity ()
ddaStreamX s sp ep = S.map (zoneOfPoint s) $ xIntercepts s sp ep

ddaStreamY :: Float -> Point2 -> Point2 -> Stream (Of (V2 Int)) Identity ()
ddaStreamY s sp ep = S.map (zoneOfPoint s) $ yIntercepts s sp ep

xIntercepts :: Float -> Point2 -> Point2 -> Stream (Of Point2) Identity ()
{-# INLINE xIntercepts #-}
xIntercepts s (V2 sx sy) (V2 ex ey) 
    | xdx == 0  = mempty
    | xdx > 0   = S.each $ zipWith V2 [sx',sx'+xdx*50..ex] ([sy',sy'+ydx*50..ey] ++ repeat ey)
    | otherwise = S.each $ zipWith V2 [sx'-50,sx'+xdx*50-50..ex-50] ([sy',sy'+ydx*50..ey] ++ repeat ey)
  where
    xdx = signum (ex - sx)
    ydx = (ey - sy) / abs (ex - sx) -- carefull: if this is zero

    sy' = sy + ydx * abs (sx - sx')
    sx' | xdx < 0 = sx - modTo s sx
        | otherwise = s + sx - modTo s sx

yIntercepts :: Float -> Point2 -> Point2 -> Stream (Of Point2) Identity ()
yIntercepts s sp ep = S.map f $ xIntercepts s (f sp) (f ep)
  where
    f (V2 x y) = V2 y x


-- | Determines which horizontal and vertical lines on a grid are crossed by a
-- line. For each adds the x-y index of the square to the right or above the
-- crossed grid line. Also adds the index of the square containing the start
-- point.
-- Not correct, eg ddaExt 10 (V2 40 50) (V2 0 0)
ddaExt :: Float -> V2 Float -> V2 Float -> IM.IntMap IS.IntSet
ddaExt s sp@(V2 sx sy) ep@(V2 ex ey) 
    | x1 <= x2 = addsp . addys . IM.fromDistinctAscList $ zip [x1 .. x2] 
                    $ map (IS.singleton . divTo s) [x1y,x1y+ydx..]
    | otherwise = addsp . addys . IM.fromDistinctAscList $ zip [x2-1 .. x1-1]
                    $ map (IS.singleton . divTo s) [x2y,x2y+ydx..]
  where
    addsp im = let V2 x y = zoneOfPoint s sp
        in insertXY im (x,y)
    x1 = divTo s sx
    x2 = divTo s ex
    x1y = fx' sp ep $ s * fromIntegral x1
    x2y = fx' sp ep $ s * fromIntegral x2
    ydx = s * ydx' sp ep
    addys m = add2s m ypairs 
    y1 = divTo s sy
    y2 = divTo s ey
    y1x = fy' sp ep $ s * fromIntegral y1
    y2x = fy' sp ep $ s * fromIntegral y2
    xdy = s * xdy' sp ep
    ypairs 
        | y1 <= y2  = zip (map (divTo s) [y1x,y1x+xdy..])
                         [y1 .. y2] 
        | otherwise = zip (map (divTo s) [y2x,y2x+xdy..])
                         [y2-1 .. y1-1] 

ydx' :: Point2 -> Point2 -> Float
{-# INLINE ydx' #-}
ydx' (V2 sx sy) (V2 ex ey)
    | sx == ex = 0
    | otherwise = (ey - sy) / (ex - sx)
fx' :: Point2 -> Point2 -> Float -> Float
{-# INLINE fx' #-}
fx' sp@(V2 sx sy) ep@(V2 _ ey) x
    | sy == ey = sy
    | otherwise = sy + ydx' sp ep * (x - sx)

xdy' :: Point2 -> Point2 -> Float
{-# INLINE xdy' #-}
xdy' (V2 sx sy) (V2 ex ey)
    | sy == ey = 0
    | otherwise = (ex - sx) / (ey - sy)
fy' :: Point2 -> Point2 -> Float -> Float
{-# INLINE fy' #-}
fy' sp@(V2 sx sy) ep@(V2 ex _) y
    | sx == ex = sx
    | otherwise = sx + xdy' sp ep * (y - sy)

add2s :: IM.IntMap IS.IntSet -> [(Int,Int)] -> IM.IntMap IS.IntSet
{-# INLINE add2s #-}
add2s = foldl' 
    (\m (k,x) -> IM.insertWith (\_ old -> IS.insert x old) k (IS.singleton x) m)

insertXY :: IM.IntMap IS.IntSet -> (Int,Int) -> IM.IntMap IS.IntSet
{-# INLINE insertXY #-}
insertXY m (k,x) = IM.insertWith (\_ old -> IS.insert x old) k (IS.singleton x) m

--addV2s :: IM.IntMap IS.IntSet -> [V2 Int] -> IM.IntMap IS.IntSet
--{-# INLINE addV2s #-}
--addV2s imis = foldl' 
--    (\m (V2 k x) -> IM.insertWith (\_ old -> IS.insert x old) k (IS.singleton x) m)
--    imis

--pairsToIntMapSet :: [V2 Int] -> IM.IntMap IS.IntSet
--pairsToIntMapSet = foldl' 
--    (\m (V2 k x) -> IM.insertWith (\_ old -> IS.insert x old) k (IS.singleton x) m)
--    IM.empty