```module Graphics.Triangulation.Triangulation where
import qualified Graphics.Triangulation.GJPTriangulation as T
import Data.Tuple.Select
import qualified Data.Vector as V
import Data.Vector (Vector, (!))
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

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.++
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)
```