KETTriangulation.hs

-- |
-- Module    : Graphics.Triangulation.KETTriangulation
-- Copyright :(C) 1997, 1998, 2008 Joern Dinkla, www.dinkla.net
--
-- Triangulation of simple polygons after Kong, Everett, Toussaint 91
-- with some changes by T.Vogt: return indices instead of coordinates of triangles
--
-- see
--     Joern Dinkla, Geometrische Algorithmen in Haskell, Diploma Thesis,
--     University of Bonn, Germany, 1998. 

module Graphics.Triangulation.KETTriangulation (ketTri) where
import List      ( (\\) )
import Data.Array (Array(..), (!), bounds)

type F2 = (Float,Float)
type Points = Array Int (Float,Float)

ketTri :: Points -> [Int] -> [(Int,Int,Int)]
ketTri points poly | (length vertices) > 3 = scan points vs stack rs
                   | otherwise = []
  where (p1:p2:p3:qs) = vertices
        vs            = qs ++ [p1]
        stack         = [p3, p2, p1, last vertices]
        rs            = reflexVertices points vertices
        vertices | polygon_direction points poly = poly -- make vertices of polygon counterclockwise
                 | otherwise                     = reverse poly

scan :: Points -> [Int] -> [Int] -> [Int] -> [(Int,Int,Int)]
scan points [] _ _                   = []
scan points [v] [x_p, x_i, _, _] rs  = [(x_i, x_p, v)]
scan points (v:vs) ss@[_,_,_] rs     = scan points vs (v:ss) rs
scan points vs@(v:vs') ss@(x_p:x_i:ss'@(x_m:x_mm:xs)) rs
  | isEar (map (points!) rs) (points!x_m) (points!x_i) (points!x_p) = (x_p, x_i, x_m) : scan points vs (x_p:ss') rs'
  | otherwise = scan points vs' (v:ss) rs
  where rs'   = rs \\ (isConvex x_m x_p v ++ isConvex x_mm x_m x_p)
	isConvex im i ip      = if isLeftTurn (points!im) (points!i) (points!ip) then [i] else []

isEar :: [F2] -> F2 -> F2 -> F2 -> Bool
isEar [] _ _ _         = True 
isEar rs m x p         = isLeftTurn m x p && not (any ( (m,x,p) `containsBNV`) rs)

reflexVertices                :: Points -> [Int] -> [Int]
reflexVertices points ps      = [ x | (m,x,p) <- angles ps, isRightTurnOrOn (points!m) (points!x) (points!p) ]

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)

containsBNV (s,t,v) p         = (a==b && b==c)
  where a                     = isLeftTurn s t p
        b                     = isLeftTurn t v p
        c                     = isLeftTurn v s p

angles                        :: [a] -> [(a,a,a)]
angles xs                     = zip3 (rotateR xs) xs (rotateL xs)

rotateL xs                    = tail xs ++ [head xs]
rotateR xs                    = [last xs] ++ init xs

-- | the direction (clockwise or counterclockwise) of a polygon can be obtained by looking at a maximal point
polygon_direction :: Points -> [Int] -> Bool
polygon_direction points poly = isLeftTurn (points!lminus) (points!l) (points!lplus)
  where l = maxim (map (points!) poly) 0 0 (0,0)
        lminus | l == fst (bounds points) = snd (bounds points)
               | otherwise = l - 1
        lplus | l == snd (bounds points) = fst (bounds points)
              | otherwise = l + 1
        -- the index of the right-/upmost point
        maxim []     count ml (mx,my) = ml
        maxim ((x,y):xs) count ml (mx,my) | (x > mx) && (y >= my) = maxim xs (count+1) count (x,y)
                                          | otherwise             = maxim xs (count+1) ml (mx,my)