triangulation 0.1 → 0.2
raw patch · 3 files changed
+141/−25 lines, 3 files
Files
- src/Graphics/Triangulation/KETTriangulation.hs +3/−20
- src/Graphics/Triangulation/Triangulation.hs +136/−3
- triangulation.cabal +2/−2
src/Graphics/Triangulation/KETTriangulation.hs view
@@ -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)
src/Graphics/Triangulation/Triangulation.hs view
@@ -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)
triangulation.cabal view
@@ -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: <www.dinkla.net/download/GeomAlgHaskell.pdf>. 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: <www.dinkla.net/download/GeomAlgHaskell.pdf>. 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