```{-# LANGUAGE MultiWayIf #-}
module Reanimate.Math.Compatible where

import           Data.List
import           Data.Maybe
import           Data.Ord
import qualified Data.Vector                   as V
import           Linear.V2
import           Linear.Vector
import           Reanimate.Math.Common
import           Reanimate.Math.Polygon
import           Reanimate.Debug
import           Reanimate.Math.Render
import           Reanimate.Svg
import           Reanimate.Animation

-- import           Debug.Trace

truncateP :: V2 Rational -> V2 Rational
truncateP = fmap (realToFrac . (realToFrac :: Rational -> Double))

mkSteinerPoints :: V2 Rational -> V2 Rational -> Int -> [V2 Rational]
mkSteinerPoints a b s_ = [ lerp (i / (s + 1)) b a | i <- [1 .. s] ]
where s = fromIntegral s_

-- 0..i,j..n-1
-- i..j
split1Link :: Polygon -> Int -> Int -> Int -> (Polygon, Polygon)
split1Link p i j s | j < i = split1Link p j i s
split1Link p i j s =
(mkPolygon \$ V.fromList left, mkPolygon \$ V.fromList right)
where
n     = pSize p
sp    = mkSteinerPoints (pAccess p i) (pAccess p j) s
left  = map (pAccess p) [0 .. i] ++ sp ++ map (pAccess p) [j .. n - 1]
right = map (pAccess p) [i .. j] ++ reverse sp

steiner2Link :: Polygon -> Int -> Int -> V2 Rational
steiner2Link p i j | j < i = steiner2Link p j i
| isNeighbour
= error "steiner2Link: Points are neighbours"
| isParent
= error "steiner2Link: Points can directly see each other."
| not (isStraightLine || isGrandparent || oneBendBetween p i j)
= error
++ show (i, j, pParent p i j, pParent p i (pParent p i j))
| otherwise
= truncateP \$ lerp 0.5 (fst vect) (intersects !! 0)
where
distToV       = approxDist (fst vect)
isNeighbour   = i == pNext p j || i == pPrev p j
isParent      = pParent p i j == i
isGrandparent = pParent p i (pParent p i j) == i
isStraightLine =
direction (pAccess p j) (pAccess p \$ pParent p i j) (pAccess p i) == 0
intersects =
sortOn distToV
\$ snd vect
: [ u
| n <- [0 .. pSize p - 1]
, let edge = (pAccess p n, pAccess p \$ pNext p n)
, u <- case rayIntersect vect edge of
Nothing -> []
Just u  -> [u]
, isBetween u edge
, u /= fst vect
, isForward vect u
]
vect
| isStraightLine
= let p1 = lerp 0.5 (pAccess p i) (pAccess p j)
p2 = case p1 - pAccess p i of
V2 x y -> p1 + V2 (-y) x -- rotate 90 degrees.
in  (p1, p2)
| otherwise
= fromMaybe
(error \$ "No window overlap: " ++ show
(isStraightLine, isGrandparent, oneBendBetween p i j)
)
\$ listToMaybe
\$ [ (p1, p1 + (p2 - p1) + (p3 - p1))
| (a, b)       <- ssspWindows iP
, (c, d)       <- ssspWindows jP
, (p1, p2, p3) <- if
| a == c    -> pure (a, b, d)
| a == d    -> pure (a, b, c)
| b == c    -> pure (b, a, d)
| b == d    -> pure (b, a, c)
| otherwise -> []
]
isForward (a, b) v = not (isBetween a (b, v))

-- 0..i,s,j..n-1
-- i..j,s
split2Link :: Polygon -> Int -> Int -> (Polygon, Polygon)
split2Link p i j | j < i = split2Link p j i
split2Link p i j = (mkPolygon \$ V.fromList left, mkPolygon \$ V.fromList right)
where
s     = steiner2Link p i j
n     = pSize p
left  = map (pAccess p) [0 .. i] ++ [s] ++ map (pAccess p) [j .. n - 1]
right = map (pAccess p) [i .. j] ++ [s]

splitNLink p i js = (mkPolygon \$ V.fromList left, mkPolygon \$ V.fromList right)
where
n        = pSize p
left = map (pAccess p) [0 .. i] ++ steiners ++ map (pAccess p) [j .. n - 1]
right    = map (pAccess p) [i .. j] ++ reverse steiners
j        = snd (last js)
steiners = splitNLinks p i js
splitNLinks _p _i []              = []
let
(l, r) = split2Link p i j
p'     = selectContains l r (pAccess p j')
s      = steiner2Link p i j
sIdx =
fromMaybe (error "missing steiner") \$ V.elemIndex s (polygonPoints p')
in
let
(l, r) = split2Link p i j
p'     = selectContains l r (pAccess p j')
s      = steiner2Link p i j
sIdx =
fromMaybe (error "missing steiner") \$ V.elemIndex s (polygonPoints p')
in
let
(l, r) = split2Link p j j'
p'     = selectContains l r (pAccess p i)
p''    = selectContains l r (pAccess p j')
s      = steiner2Link p j j'
s'     = steiner2Link p' i sIdx
sIdx   = fromMaybe (error "missing steiner sIdx")
\$ V.elemIndex s (polygonPoints p')
sIdx' = fromMaybe (error "missing steiner sIdx'")
\$ V.elemIndex s' (polygonPoints p'')
in

selectContains :: Polygon -> Polygon -> V2 Rational -> Polygon
selectContains p1 p2 elt | V.elem elt (polygonPoints p1) = p1
| V.elem elt (polygonPoints p2) = p2
| otherwise = error "elt not member of either polygons"
-- [(P, Vis)] -> Polygon -> [P]
{-
If there is a 1-link path, connect it.
Recurse into the two new polygon pairs.

If there is a 2-link path, connect it, possibly add bogus steiner point.
Recurse into the two new polygon pairs.

If there is a n-link path, connect it, add bogus steiner points.
Recurse into the two new polygon pairs.
-}

-- a 2 b 1 c
-- ab = 2-link of (a,b)
-- abc = 2-link of (cut-a-ab-b, c)
-- steps: a ab abc c
-- a, [(2,b), (1,c)]
-- [(ab, p1,p2)]

{-
, mkGroup
[ withStrokeColor "green" \$mkGroup
[ mkLine (x1,y1) (x2,y2)
, mkLine (x1',y1') (x2',y2') ]
| (a,b) <- mWins
, (i,j) <- oWins
, a == i || a == j || b == i || b == j
, not (sort [a,b] == sort [i,j])
, let V2 x1 y1 = realToFrac <\$> a
V2 x2 y2 = realToFrac <\$> b
V2 x1' y1' = realToFrac <\$> i
V2 x2' y2' = realToFrac <\$> j
]
-}

type Points = V.Vector (V2 Rational)
type Edges = [(Int, Int, Int)]
data Mesh = Mesh { meshPoints :: Points, meshEdges :: Edges }
data MeshPair = MeshPair Points Points Edges
-- The points in a RelMesh are:
--   relMeshStatic ++ x where Ax = B
-- data RelMesh = RelMesh
--   { relMeshStatic :: Points
--   , relMeshEdges  :: Edges
--   , relMeshA      :: Matrix Double
--   , relMeshB      :: Matrix Double
--   }
-- data RelMeshPair = RelMeshPair Points Edges (Matrix Double) (Matrix Double) (Matrix Double) (Matrix Double)
--
--

-- Use SSSP to find 1-link connections
-- Compute visibility windows for each node
-- Check windows for ovelap, giving 2-link connections.
-- Use floyd to generate all other connections
-- (a,b,1)
-- (a,b,2)
-- floyd
compatiblyTriangulateP :: Polygon -> Polygon -> [(Polygon, Polygon)]
compatiblyTriangulateP a b
| pSize a /= pSize b = error "polygon size mismatch"
| otherwise = traceSVG (showStep a b) \$ compatiblyTriangulateP'
(pSetOffset a 0)
(pSetOffset a 0)
(pSetOffset b 0)

showStep :: Polygon -> Polygon -> SVG
showStep a b = mkGroup
[ translate (-3) 0
\$ mkGroup [withFillColor "grey" \$ polygonShape a, polygonNumDots a]
, translate 3    0
\$ mkGroup [withFillColor "grey" \$ polygonShape b, polygonNumDots b]
]

compatiblyTriangulateP' :: Polygon -> Polygon -> Polygon -> [(Polygon, Polygon)]
compatiblyTriangulateP' aOrigin a b
| n == 3 = traceSVG (showStep a b) \$ {- trace ("Done") \$ -} [(a, b)]
| otherwise =
traceSVG (showStep a b) \$ case bestOneLink of
Nothing -> error \$ "no 2-links"
Just (nodeL, nodeR) ->
{-trace
(show
("two link" :: String, toOriginIndex nodeL, toOriginIndex nodeR)
)
\$ -}
let (aL, aR) = if (nodeL, nodeR) `elem` aOneLink
then split1Link a nodeL nodeR 1
(bL, bR) = if (nodeL, nodeR) `elem` bOneLink
then split1Link b nodeL nodeR 1
in  compatiblyTriangulateP' aOrigin aL bL
++ compatiblyTriangulateP' aOrigin aR bR
Just (nodeL, nodeR) ->
{-trace
(show ("one link" :: String, toOriginIndex nodeL, toOriginIndex nodeR)
)
\$ -}
let (aL, aR) = split1Link a nodeL nodeR 0
(bL, bR) = split1Link b nodeL nodeR 0
in  compatiblyTriangulateP' aOrigin aL bL
++ compatiblyTriangulateP' aOrigin aR bR
where
-- toOriginIndex idx =
--   ( idx
--   , fromMaybe (-1) (V.elemIndex (pAccess a idx) (polygonPoints aOrigin))
--     + polygonOffset aOrigin
--   )
n = pSize a
(sortOn (Down . nodeDist)
)
nodeDist (i, j) = min (j - i) (n - j + i)

oneBendBetween :: Polygon -> Int -> Int -> Bool
oneBendBetween _p _a _b = False
-- oneBendBetween p a b =
--     direction
--       (pAccess p b)
--       (pAccess p (pParent p a b))
--       (pAccess p \$ obstructedBy b) == 0 &&
--     (pParent p a b) /= obstructedBy b
--   where
--     obstructedBy n =
--       case pParent p a n of
--         i -> if i == a then n else obstructedBy i

polygonTwoLinks :: Polygon -> [(Int, Int)]
[ (i, j)
| i <- [0 .. n - 1]
, j <- [i + 2 .. n - 1]
, not (i == 0 && j == n - 1)
, pParent p i j /= i -- check for 1-link
, let isTwoLink      = pParent p i (pParent p i j) == i
isStraightLine = oneBendBetween p i j
-- Points on a straight line should be 2-link even though they are not
-- direct grandparents.