```module Graphics.Triangulation.Triangulation where
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)
```