module Graphics.Triangulation.Triangulation where
import Graphics.Formats.Collada.ColladaTypes
import Graphics.Formats.Collada.Transformations (cycleNeighbours,cycleN)
import qualified Graphics.Triangulation.GJPTriangulation as T
import Data.Tuple.Select
import qualified Data.Vector as V
import Data.Vector (Vector, (!))
import Graphics.Formats.Collada.Vector2D3D (V2 (V), V3(V3))
import Debug.Trace
import List

type TriangulationFunction = Vector V2 -> [(Int,Int,Int)]
data Tree = Node Int Int [Tree]

instance Show Tree where
         show (Node c p tree) = "Node " ++ (show c) ++ " " ++ (show p) ++ "[" ++ (concat(map show tree)) ++ "]"

-- | since there are a lot of triangulation algorithms
--   a triangulation function can be passed
triangulate :: TriangulationFunction -> Geometry -> Geometry
triangulate f (Geometry name prims (Vertices vname ps ns)) =
               Geometry name (map triPoly prims) (Vertices vname ps ns)
  where
  triPoly (LP (LinePrimitive pIndices          nIndices                   tex col)) =
           PL (LinePrimitive (tri 0 pIndices) (normals pIndices nIndices) tex col)
  -- TO DO: other patterns
  tri :: Int -> Vector (Vector Int) -> Vector (Vector Int)
  tri i pIndices | V.null pIndices = V.empty
                 | otherwise       = (g ( map (ind (V.head pIndices)) (f (v2s ps (V.head pIndices))))) V.++
                                     (tri (i+(V.length (V.head pIndices))) (V.tail pIndices))
  ind pIndices (i0,i1,i2) = (pIndices V.! i0, pIndices V.! i1, pIndices V.! i2)
  g :: [(Int,Int,Int)] -> Vector (Vector Int)
  g [] = V.empty
  g ((i0,i1,i2):xs) = V.cons (V.cons i0 $ V.cons i1 $ V.singleton i2) (g xs)

  normals pIndices nIndices = V.replicate (V.sum (V.map V.length pIndices)) (V.head nIndices)

v2s :: Vector V3 -> Vector Int -> Vector V2
v2s ps pIndices | V.null pIndices = V.empty
                | otherwise = V.cons (V x z) (v2s ps (V.tail pIndices))
  where (V3 x y z) = ps V.! i
        i = (V.head pIndices)

gjpTri :: Vector V2 -> [(Int,Int,Int)] 
gjpTri = T.triangulation


-- | some triangulation algorithms on't support polygons with holes
-- These polygons with (nested) holes have to be cut so that they consist of only one outline
-- I.e. the chars a,b,d,e,g,o,p,q contain holes tat have to be deleted.

deleteHoles :: Geometry -> Geometry
deleteHoles (Geometry name prims    (Vertices vname ps ns)) =
             Geometry name newPrims (Vertices vname ps ns)
  where
  newPrims = zipWith3 (\pInd tex col -> LP (LinePrimitive pInd pInd tex col)) flattenedTrees (map t prims) (map c prims)
  flattenedTrees = zipWith (flatten vs) trees vertices
  trees = map (generateTrees ps insidePoly) vertices
  pI (LP (LinePrimitive pIndices nIndices tex col)) = pIndices
  t (LP (LinePrimitive pIndices nIndices tex col)) = tex
  c (LP (LinePrimitive pIndices nIndices tex col)) = col
  vertices :: [Vector (Vector Int)]
  vertices = map pI prims
  vs = V.map (\(V3 x y z) -> V x z) ps


flatten :: Vector V2 -> [Tree] -> Vector (Vector Int) -> Vector (Vector Int)
flatten _ []                      is = V.empty
flatten vs ((Node c poly tts):ts) is
          | null tts  = V.cons                               (alternate c (pdir poly) (is V.! poly))  (flatten vs ts is)
          | otherwise = V.cons (embed vs (flatten vs tts is) (alternate c (pdir poly) (is V.! poly))) (flatten vs ts is)
  where
  pdir poly = polygonDirection $ V.map (vs V.!) (is V.! poly)

-- |cut a polygon at a good position and insert the contained hole-polygon with opposite direction
embed :: Vector V2 -> Vector (Vector Int) -> Vector Int -> Vector Int
embed vs sub_polys poly | V.null sub_polys = poly
                        | otherwise = embed vs (V.tail sub_polys) ((V.take (n+1) poly) V.++
                                                                   (V.head sub_polys) V.++
                                                                   (V.cons (V.head (V.head sub_polys)) (V.drop n poly)) )
  where n = fst $ rotatePoly (vs V.! (V.head (V.head sub_polys))) (V.map (vs V.!) poly)

-- |make sure that direction (clockwise or ccw) of polygons alternates depending on the nesting number c of poly
alternate :: Int -> Bool -> Vector Int -> Vector Int
alternate c b poly | (b && (even c)) || (not b && (odd c)) = poly
                   | otherwise                             = V.reverse poly

-- |f should be the funtion to test "contains"
-- the trees then are the hierarchy of containedness of outlines
generateTrees :: Vector V3 -> (Vector V2 -> Vector V2 -> Bool) -> Vector (Vector Int) -> [Tree]
generateTrees vs f ps | V.null ps = []
                      | otherwise = (treesList containedPolys []) ++
                                 (map (\x -> Node 0 x []) separateOutlines)
  where containedPolys = filter (\[p0,p1] -> f (pvs p0) (pvs p1)) (combi ++ (map reverse combi))
        combi = combinationsOf 2 [0..((V.length ps)-1)] -- all 2-subsets i.e. [[0,1],[0,2],[1,2]]
        -- separate outlines don't contain other outlines
        separateOutlines = ([0..((V.length ps)-1)]) \\ (nub $ concat containedPolys)
        pvs p = V.map (\(V3 x y z) -> V x z) $ V.map (vs V.!) (ps V.! p)

treesList :: [[Int]] -> [Tree] -> [Tree]
treesList [] trees = trees
treesList ([x,y]:cs) trees = treesList cs (insertTrees [x,y] trees)

insertTrees :: [Int] -> [Tree] -> [Tree]
insertTrees [x,y] trees | or (map fst ins) = map snd ins
                        | otherwise = (map snd ins) ++ [ Node 0 y [Node 1 x []] ]
  where ins = map (insertTree [x,y]) trees

insertTree :: [Int] -> Tree -> (Bool, Tree)
insertTree [x,y] (Node c i []) | y == i = (True, Node c i [Node (c+1) x []] )
                               | otherwise = (False, Node c i [])
insertTree [x,y] (Node c i trees) | y == i = (True, Node c i ((Node (c+1) x []):trees) )
                                  | otherwise = (b, Node c i (map snd subtrees))
  where subtrees = map (insertTree [x,y]) trees
        b = or (map fst subtrees)

-- subsets of size k
-- borrowed from David Amos' library: HaskellForMaths
combinationsOf 0 _ = [[]]
combinationsOf _ [] = []
combinationsOf k (x:xs) = map (x:) (combinationsOf (k-1) xs) ++ combinationsOf k xs

-- |how many positions to rotate a polygon until the start point is nearest to some other point
-- call i.e. with nearest (3,4) [(0,0),(1,2), ... ] 0 0
rotatePoly :: V2 -> Vector V2  -> (Int,Float)
rotatePoly p poly = nearest p poly (-1) 0 0

nearest :: V2 -> Vector V2 -> Float -> Int -> Int -> (Int,Float)
nearest (V x0 y0) ps dist l ml | V.null ps = (ml,dist)
                               | (newDist < dist) || (dist < 0) = nearest (V x0 y0) (V.tail ps) newDist (l+1) l
                               | otherwise                      = nearest (V x0 y0) (V.tail ps) dist    (l+1) ml
  where newDist = (x0-x1)*(x0-x1)+(y0-y1)*(y0-y1)
        (V x1 y1) = V.head ps

-- | returns True iff the first point of the first polygon is inside the second poylgon
insidePoly :: Vector V2 -> Vector V2 -> Bool
insidePoly poly1 poly2 | V.null poly1 = False
                       | V.null poly2 = False
                       | otherwise    = pointInside (V.head poly1) poly2

-- | A point is inside a polygon if it has an odd number of intersections with the boundary (Jordan Curve theorem)
pointInside :: V2 -> Vector V2 -> Bool
pointInside (V x y) poly = (V.length intersectPairs) `mod` 2 == 1
  where intersectPairs = V.filter (\p -> positiveXAxis p && aboveBelow p) allPairs --, specialCases p]
        allPairs = cycleNeighbours poly
        positiveXAxis p = (x0 p) > x || (x1 p) > x -- intersect with positive x-axis
                                                   -- only lines with one point above + one point below can intersect
        aboveBelow p = (((y0 p)> y && (y1 p)< y) || ((y0 p) < y && (y1 p) > y))
        specialCases p = (((dir1 p) > 0 && (dir2 p) <= 0) || ((dir1 p) <= 0 && (dir2 p) > 0))-- cross product for special cases
        dir1 p = cross ((x1 p)-(x0 p),(y1 p)-(y0 p)) (1,0)
        dir2 p = cross ((x1 p)-(x0 p),(y1 p)-(y0 p)) (x-(x0 p),y-(y0 p))
        cross (a0,b0) (a1,b1) = a0*b1 - a1*b0
        x0 p = (\(V x y) -> x) (V.head p)
        x1 p = (\(V x y) -> x) (V.last p)
        y0 p = (\(V x y) -> y) (V.head p)
        y1 p = (\(V x y) -> y) (V.last p)

-- | the direction of a polygon can be obtained by looking at a maximal point
-- returns True if counterclockwise
--         False if clockwise
polygonDirection :: Vector V2 -> Bool
polygonDirection poly | dir > 0 = True
                      | dir < 0 = False
                      | dir ==0 = (x0 > x1) || (y0 < y1)
  where p = V.fromList $ nub $ V.toList poly
        (V x0 y0) = p V.! lminus
        (V x1 y1) = p V.! lplus
        dir = area2 (p!lminus) (p!l) (p!lplus)
        l = maxim poly 0 0 (-1000000,-1000000)
        lminus = (l-1) `mod` (V.length p)
        lplus = (l+1) `mod` (V.length p)
         -- the index of the right-/upmost point

maxim :: Vector V2 -> Int -> Int -> (Float,Float) -> Int
maxim xs count ml (mx,my) | V.null xs = ml
                          | (x > mx) || (x >= mx && y > my) = maxim (V.tail xs) (count+1) count (x, y)
                          | otherwise                       = maxim (V.tail xs) (count+1) ml (mx,my)
  where (V x y) = V.head xs

isRightTurnOrOn m x p = (area2 m x p) <= 0
isLeftTurn :: V2 -> V2 -> V2 -> Bool
isLeftTurn      m x p = (area2 m x p) > 0
area2 (V x2 y2) (V x0 y0) (V x1 y1) = (x1-x0)*(y2-y0)-(x2-x0)*(y1-y0)