module Graphics.Triangulation.Triangulation where
import Graphics.Formats.Collada.ColladaTypes
import Graphics.Formats.Collada.Transformations (cycleNeighbours,cycleN)
import Data.Array (Array(..), listArray, (!))
import Debug.Trace
import List

type Points = Array Int (Float,Float)
type TriangulationFunction = Points -> [Int] -> [(Int,Int,Int)]
data Tree = Node Int Int [Tree]
type F2 = (Float,Float)

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 pIndices) (normals pIndices nIndices) tex col)
  -- TO DO: other patterns
  tri pIndices = map (\(x,y,z) -> [x,y,z]) (concat (map (f arr) pIndices) )
  normals pIndices nIndices = replicate (length (concat pIndices)) (head nIndices) -- TO DO: Why not (tri pIndices)
  arr = listArray (0,l-1) $ map (\(x,y,z) -> (x,z)) ps
  l = length ps

-- | 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 trees indices
  arr = listArray (0,l-1) $ map (\(x,y,z) -> (x,z)) ps
  l = length ps
  trees = map (generateTrees arr insidePoly) indices
  pI (LP (LinePrimitive pIndices nIndices tex col)) = pIndices
  t (LP (LinePrimitive pIndices nIndices tex col)) = tex
  c (LP (LinePrimitive pIndices nIndices tex col)) = col
  indices = map pI prims

  flatten :: [Tree] -> [[Int]] -> [[Int]]
  flatten []                    is = []
  flatten ((Node c poly []):ts) is =                            (alternate c (pdir (is!!poly)) (is!!poly)) : (flatten ts is)
  flatten ((Node c poly ps):ts) is = (embed arr (flatten ps is) (alternate c (pdir (is!!poly)) (is!!poly))): (flatten ts is)
  pdir poly = polygonDirection (map (arr!) poly)

-- |cut a polygon at a good position and insert the contained hole-polygon with opposite direction
embed :: Points -> [[Int]] -> [Int] -> [Int]
embed _      []            poly = poly
embed points (s:sub_polys) poly = embed points sub_polys ((take (n+1) poly) ++ s ++ [head s] ++ (drop n poly))
  where n = fst (rotatePoly (head s) points poly)

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

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

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

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

insertTree :: Points -> [Int] -> Tree -> (Bool, Tree)
insertTree points [x,y] (Node c i []) | y == i = (True, Node c i [Node (c+1) x []] )
                                      | otherwise = (False, Node c i [])
insertTree points [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 points [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 :: Int -> Points -> [Int] -> (Int,Float)
rotatePoly p points poly = (fst tup, snd tup)
  where tup = nearest (points!p) (map (points!) poly) (-1) 0 0

nearest :: F2 -> [F2] -> Float -> Int -> Int -> (Int,Float)
nearest _       []           dist l ml = (ml,dist)
nearest (x0,y0) ((x1,y1):ps) dist l ml | (newDist < dist) || (dist < 0) = nearest (x0,y0) ps newDist (l+1) l
                                       | otherwise                      = nearest (x0,y0) ps dist    (l+1) ml
  where newDist = (x0-x1)*(x0-x1)+(y0-y1)*(y0-y1)

-- | returns True iff the first point of the first polygon is inside the second poylgon
insidePoly :: Points -> [Int] -> [Int] -> Bool
insidePoly _ [] _ = False
insidePoly _ _ [] = False
insidePoly points poly1 poly2 = pointInside (points!(head poly1)) (map (points!) poly2)

-- | A point is inside a polygon if it has an odd number of intersections with the boundary (Jordan Curve theorem)
pointInside :: F2 -> [F2] -> Bool
pointInside (x,y) poly = (length intersectPairs) `mod` 2 == 1
  where intersectPairs = [ p | p <- allPairs, positiveXAxis p, aboveBelow p] --, 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 = fst (head p)
        x1 p = fst (last p)
        y0 p = snd (head p)
        y1 p = snd (last p)

-- | the direction of a polygon can be obtained by looking at a maximal point
-- returns True if counterclockwise
--         False if clockwise
polygonDirection :: [F2] -> Bool
polygonDirection poly | dir > 0 = True
                      | dir < 0 = False
                      | dir ==0 = (fst (p!!lminus) > fst (p!!lplus)) || (snd (p!!lminus) < snd (p!!lplus))
 where p = nub poly
       dir = area2 (p!!lminus) (p!!l) (p!!lplus)
       l = maxim p 0 0 (0,0)
       lminus = (l-1) `mod` (length p)
       lplus = (l+1) `mod` (length p)
        -- the index of the right-/upmost point
       maxim []     count ml (mx,my) = ml
       maxim ((x,y):xs) count ml (mx,my) | (x > mx) || (x >= mx && (y > my)) = maxim xs (count+1) count (x,y)
                                         | otherwise                          = maxim xs (count+1) ml (mx,my)
isRightTurnOrOn m x p = (area2 m x p) <= 0
isLeftTurn :: F2 -> F2 -> F2 -> Bool
isLeftTurn      m x p = (area2 m x p) > 0
area2 (x2,y2) (x0,y0) (x1,y1) = (x1-x0)*(y2-y0)-(x2-x0)*(y1-y0)