diff --git a/src/Graphics/Triangulation/KETTriangulation.hs b/src/Graphics/Triangulation/KETTriangulation.hs --- a/src/Graphics/Triangulation/KETTriangulation.hs +++ b/src/Graphics/Triangulation/KETTriangulation.hs @@ -10,6 +10,7 @@ -- University of Bonn, Germany, 1998. module Graphics.Triangulation.KETTriangulation (ketTri) where +import Graphics.Triangulation.Triangulation (polygonDirection, isLeftTurn, isRightTurnOrOn) import List ( (\\) ) import Data.Array (Array(..), (!), bounds) @@ -23,8 +24,8 @@ 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 + vertices | polygonDirection (map (points!) poly) = poly -- make vertices of polygon counterclockwise + | otherwise = reverse poly scan :: Points -> [Int] -> [Int] -> [Int] -> [(Int,Int,Int)] scan points [] _ _ = [] @@ -43,11 +44,6 @@ 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 @@ -58,16 +54,3 @@ 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) diff --git a/src/Graphics/Triangulation/Triangulation.hs b/src/Graphics/Triangulation/Triangulation.hs --- a/src/Graphics/Triangulation/Triangulation.hs +++ b/src/Graphics/Triangulation/Triangulation.hs @@ -1,17 +1,150 @@ module Graphics.Triangulation.Triangulation where import Graphics.Formats.Collada.ColladaTypes -import Data.Array (Array(..),listArray) +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) + 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) + 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) diff --git a/triangulation.cabal b/triangulation.cabal --- a/triangulation.cabal +++ b/triangulation.cabal @@ -1,7 +1,7 @@ Name: triangulation -Version: 0.1 +Version: 0.2 Synopsis: triangulation of polygons -Description: An implementation of a simple triangulation algorithm for polygons without holes, crossings (and maybe other anomalies that I am not aware of). The code is explained in this diploma thesis: . The original author made a very big library that needs a long time to compile. Thats why only one algorithm was extracted and freed from a big net of inner dependencies and types. +Description: An implementation of a simple triangulation algorithm for polygons without crossings. The code is explained in this diploma thesis: . The original author made a very big library that needs a long time to compile. Thats why only one algorithm was extracted and freed from a big net of inner dependencies and types. category: Graphics License: GPL License-file: LICENSE