Port to diagrams >1.3

We had to add NoFace/NoEdge/NoVert constructors to our half-edge
data structures, because using "Maybe HeEdge a" as the result value
of 'indirectToDirect' causes an infinite recursion, since the whole
data structure (which is cyclic and infinite) has to be evaluated
in order to know which constructor to use.

Unfortunately this makes the code quite hard to read.

TODO: add pseudo-code

.gitignore

@@ -11,3 +11,10 @@ dist/
 # cabal
 .cabal-sandbox/
 cabal.sandbox.config
+
+# profiling
+*.prof
+
+_darcs/
+
+.liquid/

Algebra/Polygon.hs

@@ -4,15 +4,14 @@ module Algebra.Polygon where
 
 import Algebra.Vector
 import Data.Maybe
-import Diagrams.TwoD.Types
 import MyPrelude
 
 
 -- |Split a polygon by a given segment which must be vertices of the
 -- polygon (returns empty array otherwise).
-splitPoly :: [P2]
+splitPoly :: [PT]
-          -> (P2, P2)
+          -> Segment
-          -> [[P2]]
+          -> [[PT]]
 splitPoly pts (a, b)
   | elem a pts && elem b pts =
     [b : takeWhile (/= b) shiftedPoly, a : dropWhile (/= b) shiftedPoly]
@@ -22,7 +21,7 @@ splitPoly pts (a, b)
 
 
 -- |Get all edges of a polygon.
-polySegments :: [P2] -> [(P2, P2)]
+polySegments :: [PT] -> [Segment]
 polySegments p@(x':_:_:_) = go p ++ [(last p, x')]
   where
     go (x:y:xs) = (x, y) : go (y:xs)
@@ -33,7 +32,7 @@ polySegments _ = []
 -- |Check whether the given segment is inside the polygon.
 -- This doesn't check for segments that are completely outside
 -- of the polygon yet.
-isInsidePoly :: [P2] -> (P2, P2) -> Bool
+isInsidePoly :: [PT] -> Segment -> Bool
 isInsidePoly pts seg =
   null
     . catMaybes
@@ -42,21 +41,21 @@ isInsidePoly pts seg =
 
 
 -- |Check whether two points are adjacent vertices of a polygon.
-adjacent :: P2 -> P2 -> [P2] -> Bool
+adjacent :: PT -> PT -> [PT] -> Bool
 adjacent u v = any (\x -> x == (u, v) || x == (v, u)) . polySegments
 
 
 -- |Check whether the polygon is a triangle polygon.
-isTrianglePoly :: [P2] -> Bool
+isTrianglePoly :: [PT] -> Bool
 isTrianglePoly [_, _, _] = True
 isTrianglePoly _         = False
 
 
 -- |Get all triangle polygons.
-triangleOnly :: [[P2]] -> [[P2]]
+triangleOnly :: [[PT]] -> [[PT]]
 triangleOnly = filter isTrianglePoly
 
 
 -- |Get all non-triangle polygons.
-nonTriangleOnly :: [[P2]] -> [[P2]]
+nonTriangleOnly :: [[PT]] -> [[PT]]
 nonTriangleOnly = filter (not . isTrianglePoly)

Algebra/Vector.hs

@@ -13,6 +13,13 @@ import GHC.Float
 import MyPrelude
 
 
+type Vec     = V2 Double
+type PT      = P2 Double
+type Coord   = (Double, Double)
+type Segment = (PT, PT)
+type Square  = (Coord, Coord)
+
+
 data Alignment = CW
                | CCW
                | CL
@@ -24,13 +31,13 @@ data Alignment = CW
 -- ((xmin, ymin), (xmax, ymax))
 dimToSquare :: (Double, Double) -- ^ x dimension
             -> (Double, Double) -- ^ y dimension
-            -> ((Double, Double), (Double, Double)) -- ^ square describing those dimensions
+            -> Square           -- ^ square describing those dimensions
 dimToSquare (x1, x2) (y1, y2) = ((x1, y1), (x2, y2))
 
 
 -- |Checks whether the Point is in a given Square.
-inRange :: ((Double, Double), (Double, Double)) -- ^ the square: ((xmin, ymin), (xmax, ymax))
+inRange :: Square -- ^ the square: ((xmin, ymin), (xmax, ymax))
-        -> P2     -- ^ Coordinate
+        -> PT     -- ^ Coordinate
         -> Bool   -- ^ result
 inRange ((xmin, ymin), (xmax, ymax)) (coords -> x :& y)
   = x >= min xmin xmax
@@ -40,7 +47,7 @@ inRange ((xmin, ymin), (xmax, ymax)) (coords -> x :& y)
 
 
 -- |Get the angle between two vectors.
-getAngle :: R2 -> R2 -> Double
+getAngle :: Vec -> Vec -> Double
 getAngle a b =
   acos
     . flip (/) (vecLength a * vecLength b)
@@ -49,50 +56,48 @@ getAngle a b =
 
 
 -- |Get the length of a vector.
-vecLength :: R2 -> Double
+vecLength :: Vec -> Double
 vecLength v = sqrt (x^(2 :: Int) + y^(2 :: Int))
   where
     (x, y) = unr2 v
 
 
 -- |Compute the scalar product of two vectors.
-scalarProd :: R2 -> R2 -> Double
+scalarProd :: Vec -> Vec -> Double
-scalarProd (R2 a1 a2) (R2 b1 b2) = a1 * b1 + a2 * b2
+scalarProd (V2 a1 a2) (V2 b1 b2) = a1 * b1 + a2 * b2
 
 
 -- |Multiply a scalar with a vector.
-scalarMul :: Double -> R2 -> R2
+scalarMul :: Double -> Vec -> Vec
-scalarMul d (R2 a b) = R2 (a * d) (b * d)
+scalarMul d (V2 a b) = V2 (a * d) (b * d)
 
 
 -- |Construct a vector that points to a point from the origin.
-pt2Vec :: P2 -> R2
+pt2Vec :: PT -> Vec
 pt2Vec = r2 . unp2
 
 
 -- |Give the point which is at the coordinates the vector
 -- points to from the origin.
-vec2Pt :: R2 -> P2
+vec2Pt :: Vec -> PT
 vec2Pt = p2 . unr2
 
 
 -- |Construct a vector between two points.
-vp2 :: P2  -- ^ vector origin
+vp2 :: PT  -- ^ vector origin
-    -> P2  -- ^ vector points here
+    -> PT  -- ^ vector points here
-    -> R2
+    -> Vec
 vp2 a b = pt2Vec b - pt2Vec a
 
 
 -- |Computes the determinant of 3 points.
-det :: P2 -> P2 -> P2 -> Double
+det :: PT -> PT -> PT -> Double
 det (coords -> ax :& ay) (coords -> bx :& by) (coords -> cx :& cy) =
   (bx - ax) * (cy - ay) - (by - ay) * (cx - ax)
 
 
 -- |Get the point where two lines intesect, if any.
-intersectSeg' :: (P2, P2) -- ^ first segment
+intersectSeg' :: Segment -> Segment -> Maybe PT
-              -> (P2, P2) -- ^ second segment
-              -> Maybe P2
 intersectSeg' (a, b) (c, d) =
   glossToPt <$> intersectSegSeg (ptToGloss a)
                                 (ptToGloss b)
@@ -105,7 +110,7 @@ intersectSeg' (a, b) (c, d) =
 
 -- |Get the point where two lines intesect, if any. Excludes the
 -- case of end-points intersecting.
-intersectSeg'' :: (P2, P2) -> (P2, P2) -> Maybe P2
+intersectSeg'' :: Segment -> Segment -> Maybe PT
 intersectSeg'' (a, b) (c, d) = case intersectSeg' (a, b) (c, d) of
   Just x  -> if x `notElem` [a,b,c,d] then Just a else Nothing
   Nothing -> Nothing
@@ -115,7 +120,7 @@ intersectSeg'' (a, b) (c, d) = case intersectSeg' (a, b) (c, d) of
 -- * clock-wise
 -- * counter-clock-wise
 -- * collinear
-getOrient :: P2 -> P2 -> P2 -> Alignment
+getOrient :: PT -> PT -> PT -> Alignment
 getOrient a b c = case compare (det a b c) 0 of
   LT -> CW
   GT -> CCW
@@ -125,7 +130,7 @@ getOrient a b c = case compare (det a b c) 0 of
 --- |Checks if 3 points a,b,c do not build a clockwise triangle by
 --- connecting a-b-c. This is done by computing the determinant and
 --- checking the algebraic sign.
-notcw :: P2 -> P2 -> P2 -> Bool
+notcw :: PT -> PT -> PT -> Bool
 notcw a b c = case getOrient a b c of
   CW -> False
   _  -> True
@@ -134,22 +139,22 @@ notcw a b c = case getOrient a b c of
 --- |Checks if 3 points a,b,c do build a clockwise triangle by
 --- connecting a-b-c. This is done by computing the determinant and
 --- checking the algebraic sign.
-cw :: P2 -> P2 -> P2 -> Bool
+cw :: PT -> PT -> PT -> Bool
 cw a b c = not . notcw a b $ c
 
 
 -- |Sort X and Y coordinates lexicographically.
-sortedXY :: [P2] -> [P2]
+sortedXY :: [PT] -> [PT]
 sortedXY = fmap p2 . sortLex . fmap unp2
 
 
 -- |Sort Y and X coordinates lexicographically.
-sortedYX :: [P2] -> [P2]
+sortedYX :: [PT] -> [PT]
 sortedYX = fmap p2 . sortLexSwapped . fmap unp2
 
 
 -- |Sort all points according to their X-coordinates only.
-sortedX :: [P2] -> [P2]
+sortedX :: [PT] -> [PT]
 sortedX xs =
   fmap p2
   . sortBy (\(a1, _) (a2, _) -> compare a1 a2)
@@ -157,7 +162,7 @@ sortedX xs =
 
 
 -- |Sort all points according to their Y-coordinates only.
-sortedY :: [P2] -> [P2]
+sortedY :: [PT] -> [PT]
 sortedY xs =
   fmap p2
   . sortBy (\(_, b1) (_, b2) -> compare b1 b2)
@@ -165,25 +170,25 @@ sortedY xs =
 
 
 -- |Apply a function on the coordinates of a point.
-onPT :: ((Double, Double) -> (Double, Double)) -> P2 -> P2
+onPT :: (Coord -> Coord) -> PT -> PT
 onPT f = p2 . f . unp2
 
 
 -- |Compare the y-coordinate of two points.
-ptCmpY :: P2 -> P2 -> Ordering
+ptCmpY :: PT -> PT -> Ordering
 ptCmpY (coords -> _ :& y1) (coords -> _ :& y2) =
   compare y1 y2
 
 
 -- |Compare the x-coordinate of two points.
-ptCmpX :: P2 -> P2 -> Ordering
+ptCmpX :: PT -> PT -> Ordering
 ptCmpX (coords -> x1 :& _) (coords -> x2 :& _) =
   compare x1 x2
 
 
-posInfPT :: P2
+posInfPT :: PT
 posInfPT = p2 (read "Infinity", read "Infinity")
 
 
-negInfPT :: P2
+negInfPT :: PT
 negInfPT = p2 (negate . read $ "Infinity", negate . read $ "Infinity")

Algorithms/GrahamScan.hs

@@ -3,7 +3,6 @@
 module Algorithms.GrahamScan where
 
 import Algebra.Vector
-import Diagrams.TwoD.Types
 import MyPrelude
 
 
@@ -75,18 +74,18 @@ ys = []
 return [(100, 100), (400, 200)]
 =========================================================
 --}
-grahamCH :: [P2] -> [P2]
+grahamCH :: [PT] -> [PT]
 grahamCH vs = grahamUCH vs ++ (tailInit . grahamLCH $ vs)
 
 
 -- |Get the lower part of the convex hull.
-grahamLCH :: [P2] -> [P2]
+grahamLCH :: [PT] -> [PT]
 grahamLCH vs = uncurry (\x y -> last . scanH x $ y)
                        (first reverse . splitAt 3 . sortedXY $ vs)
 
 
 -- |Get the upper part of the convex hull.
-grahamUCH :: [P2] -> [P2]
+grahamUCH :: [PT] -> [PT]
 grahamUCH vs = uncurry (\x y -> last . scanH x $ y)
                        (first reverse . splitAt 3 . reverse . sortedXY $ vs)
 
@@ -96,9 +95,9 @@ grahamUCH vs = uncurry (\x y -> last . scanH x $ y)
 -- If it's the upper or lower half depends on the input.
 -- Also, the first list is expected to be reversed since we only care
 -- about the last 3 elements and want to stay efficient.
-scanH :: [P2]   -- ^ the first 3 starting points in reversed order
+scanH :: [PT]   -- ^ the first 3 starting points in reversed order
-      -> [P2]   -- ^ the rest of the points
+      -> [PT]   -- ^ the rest of the points
-      -> [[P2]] -- ^ all convex hull points iterations for the half
+      -> [[PT]] -- ^ all convex hull points iterations for the half
 scanH hs@(x:y:z:xs) (r':rs')
   |	notcw z y x       = hs : scanH (r':hs) rs'
   | otherwise         = hs : scanH (x:z:xs) (r':rs')
@@ -112,12 +111,12 @@ scanH hs _            = [hs]
 -- |Compute all steps of the graham scan algorithm to allow
 -- visualizing it.
 -- Whether the upper or lower hull is computed depends on the input.
-grahamCHSteps :: Int -> [P2] -> [P2] -> [[P2]]
+grahamCHSteps :: Int -> [PT] -> [PT] -> [[PT]]
 grahamCHSteps c xs' ys' = take c . scanH xs' $ ys'
 
 
 -- |Get all iterations of the upper hull of the graham scan algorithm.
-grahamUHSteps :: [P2] -> [[P2]]
+grahamUHSteps :: [PT] -> [[PT]]
 grahamUHSteps vs =
   (++) [getLastX 2 . sortedXY $ vs]
   . rmdups
@@ -128,7 +127,7 @@ grahamUHSteps vs =
 
 
 -- |Get all iterations of the lower hull of the graham scan algorithm.
-grahamLHSteps :: [P2] -> [[P2]]
+grahamLHSteps :: [PT] -> [[PT]]
 grahamLHSteps vs =
   (++) [take 2 . sortedXY $ vs]
   . rmdups

Algorithms/KDTree.hs

@@ -42,9 +42,9 @@ instance Not Direction where
 
 
 -- |Construct a kd-tree from a list of points in O(n log n).
-kdTree :: [P2]      -- ^ list of points to construct the kd-tree from
+kdTree :: [PT]      -- ^ list of points to construct the kd-tree from
        -> Direction -- ^ initial direction of the root-node
-       -> KDTree P2 -- ^ resulting kd-tree
+       -> KDTree PT -- ^ resulting kd-tree
 kdTree xs' = go (sortedX xs') (sortedY xs')
   where
     go [] _ _ = KTNil
@@ -67,10 +67,10 @@ kdTree xs' = go (sortedX xs') (sortedY xs')
 -- If you want to partition against the pivot of X, then you pass
 --   partition' (pivot xs) (ys, xs)
 -- and get ((y1, y2), (x1, x2)).
-partition' :: P2                           -- ^ the pivot to partition against
+partition' :: PT                           -- ^ the pivot to partition against
-           -> (P2 -> P2 -> Ordering)       -- ^ ptCmpY or ptCmpX
+           -> (PT -> PT -> Ordering)       -- ^ ptCmpY or ptCmpX
-           -> ([P2], [P2])                 -- ^ both lists (X, Y) or (Y, X)
+           -> ([PT], [PT])                 -- ^ both lists (X, Y) or (Y, X)
-           -> (([P2], [P2]), ([P2], [P2])) -- ^ ((x1, x2), (y1, y2)) or
+           -> (([PT], [PT]), ([PT], [PT])) -- ^ ((x1, x2), (y1, y2)) or
                                            --   ((y1, y2), (x1, x2))
 partition' piv cmp' (xs, ys) = ((x1, x2), (y1, y2))
   where
@@ -83,16 +83,16 @@ partition' piv cmp' (xs, ys) = ((x1, x2), (y1, y2))
 -- |Partition two sorted lists of points X and Y against the pivot of
 -- Y. This function is unsafe as it does not check if there is a valid
 -- pivot.
-partitionY :: ([P2], [P2])                 -- ^ both lists (X, Y)
+partitionY :: ([PT], [PT])                 -- ^ both lists (X, Y)
-           -> (([P2], [P2]), ([P2], [P2])) -- ^ ((x1, x2), (y1, y2))
+           -> (([PT], [PT]), ([PT], [PT])) -- ^ ((x1, x2), (y1, y2))
 partitionY (xs, ys) = partition' (fromJust . pivot $ ys) ptCmpY (xs, ys)
 
 
 -- |Partition two sorted lists of points X and Y against the pivot of
 -- X. This function is unsafe as it does not check if there is a valid
 -- pivot.
-partitionX :: ([P2], [P2])                 -- ^ both lists (X, Y)
+partitionX :: ([PT], [PT])                 -- ^ both lists (X, Y)
-           -> (([P2], [P2]), ([P2], [P2])) -- ^ ((x1, x2), (y1, y2))
+           -> (([PT], [PT]), ([PT], [PT])) -- ^ ((x1, x2), (y1, y2))
 partitionX (xs, ys) = (\(x, y) -> (y, x))
                       . partition' (fromJust . pivot $ xs) ptCmpX $ (ys, xs)
 
@@ -100,9 +100,7 @@ partitionX (xs, ys) = (\(x, y) -> (y, x))
 -- |Execute a range search in O(log n). It returns a tuple
 -- of the points found in the range and also gives back a pretty
 -- rose tree suitable for printing.
-rangeSearch :: KDTree P2                            -- ^ tree to search in
+rangeSearch :: KDTree PT -> Square -> ([PT], Tree String)
-            -> ((Double, Double), (Double, Double)) -- ^ square describing the range
-            -> ([P2], Tree String)
 rangeSearch kd' sq' = (goPt kd' sq', goTree kd' sq' True)
   where
     -- either y1 or x1 depending on the orientation
@@ -112,7 +110,7 @@ rangeSearch kd' sq' = (goPt kd' sq', goTree kd' sq' True)
     -- either the second or first of the tuple, depending on the orientation
     cur' dir = if' (dir == Vertical) snd fst
     -- All points in the range.
-    goPt :: KDTree P2 -> ((Double, Double), (Double, Double)) -> [P2]
+    goPt :: KDTree PT -> Square -> [PT]
     goPt KTNil _ = []
     goPt (KTNode ln pt dir rn) sq =
       [pt | inRange sq pt]
@@ -124,7 +122,7 @@ rangeSearch kd' sq' = (goPt kd' sq', goTree kd' sq' True)
               (goPt rn sq)
               [])
     -- A pretty rose tree suitable for printing.
-    goTree :: KDTree P2 -> ((Double, Double), (Double, Double)) -> Bool -> Tree String
+    goTree :: KDTree PT -> Square -> Bool -> Tree String
     goTree KTNil _ _ = Node "nil" []
     goTree (KTNode ln pt dir rn) sq vis
       | ln == KTNil && rn == KTNil = Node treeText []
@@ -181,7 +179,7 @@ getDirection _                  = Nothing
 
 
 -- |Convert a kd-tree to a rose tree, for pretty printing.
-kdTreeToRoseTree :: KDTree P2 -> Tree String
+kdTreeToRoseTree :: KDTree PT -> Tree String
 kdTreeToRoseTree (KTNil) = Node "nil" []
 kdTreeToRoseTree (KTNode ln val _ rn) =
   Node (show . unp2 $ val) [kdTreeToRoseTree ln, kdTreeToRoseTree rn]

Algorithms/PolygonIntersection.hs

@@ -18,14 +18,14 @@ import           QueueEx
 -- successor are saved for convenience.
 data PolyPT =
   PolyA {
-    id'   :: P2
+    id'   :: PT
-    , pre :: P2
+    , pre :: PT
-    , suc :: P2
+    , suc :: PT
     }
   | PolyB {
-    id'   :: P2
+    id'   :: PT
-    , pre :: P2
+    , pre :: PT
-    , suc :: P2
+    , suc :: PT
     }
   deriving (Show, Eq)
 
@@ -42,7 +42,7 @@ isPolyB = not . isPolyA
 -- |Shift a list of sorted convex hull points of a polygon so that
 -- the first element in the list is the one with the highest y-coordinate.
 -- This is done in O(n).
-sortLexPoly :: [P2] -> [P2]
+sortLexPoly :: [PT] -> [PT]
 sortLexPoly ps = maybe [] (`shiftM` ps) (elemIndex (yMax ps) ps)
   where
     yMax       = foldl1 (\x y -> if ptCmpY x y == GT then x else y)
@@ -50,8 +50,8 @@ sortLexPoly ps = maybe [] (`shiftM` ps) (elemIndex (yMax ps) ps)
 
 -- |Make a PolyPT list out of a regular list of points, so
 -- the predecessor and successors are all saved.
-mkPolyPTList :: (P2 -> P2 -> P2 -> PolyPT) -- ^ PolyA or PolyB function
+mkPolyPTList :: (PT -> PT -> PT -> PolyPT) -- ^ PolyA or PolyB function
-             -> [P2]                       -- ^ polygon points
+             -> [PT]                       -- ^ polygon points
              -> [PolyPT]
 mkPolyPTList f' pts@(x':y':_:_) =
   f' x' (last pts) y' : go f' pts
@@ -64,7 +64,7 @@ mkPolyPTList _ _ = []
 
 -- |Sort the points of two polygons according to their y-coordinates,
 -- while saving the origin of that point. This is done in O(n).
-sortLexPolys :: ([P2], [P2]) -> [PolyPT]
+sortLexPolys :: ([PT], [PT]) -> [PolyPT]
 sortLexPolys (pA'@(_:_), pB'@(_:_)) =
   queueToList $ go (Q.fromList . mkPolyPTList PolyA . sortLexPoly $ pA')
                    (Q.fromList . mkPolyPTList PolyB . sortLexPoly $ pB')
@@ -104,7 +104,7 @@ sortLexPolys _ = []
 
 -- |Get all points that intersect between both polygons. This is done
 -- in O(n).
-intersectionPoints :: [PolyPT] -> [P2]
+intersectionPoints :: [PolyPT] -> [PT]
 intersectionPoints xs' = rmdups . go $ xs'
   where
     go [] = []
@@ -113,7 +113,7 @@ intersectionPoints xs' = rmdups . go $ xs'
 
     -- Get the scan line or in other words the
     -- Segment pairs we are going to check for intersection.
-    scanLine :: [PolyPT] -> ([(P2, P2)], [(P2, P2)])
+    scanLine :: [PolyPT] -> ([Segment], [Segment])
     scanLine sp@(_:_) = (,) (getSegment isPolyA) (getSegment isPolyB)
       where
         getSegment f = fromMaybe []
@@ -124,7 +124,7 @@ intersectionPoints xs' = rmdups . go $ xs'
     -- Gets the actual intersections between the segments of
     -- both polygons we currently examine. This is done in O(1)
     -- since we have max 4 segments.
-    segIntersections :: ([(P2, P2)], [(P2, P2)]) -> [P2]
+    segIntersections :: ([Segment], [Segment]) -> [PT]
     segIntersections (a@(_:_), b@(_:_)) =
       catMaybes
       . fmap (\[x, y] -> intersectSeg' x y)

Algorithms/PolygonTriangulation.hs

@@ -6,7 +6,6 @@ import Algebra.Polygon
 import Algebra.Vector
 import qualified Control.Arrow as A
 import Data.Maybe
-import Diagrams.TwoD.Types
 import Safe
 
 
@@ -19,12 +18,12 @@ data VCategory = VStart
 
 
 -- |Classify all vertices on a polygon into five categories (see VCategory).
-classifyList :: [P2] -> [(P2, VCategory)]
+classifyList :: [PT] -> [(PT, VCategory)]
 classifyList p@(x:y:_:_) =
   -- need to handle the first and last element separately
   [classify (last p) x y] ++ go p ++ [classify (last . init $ p) (last p) x]
   where
-    go :: [P2] -> [(P2, VCategory)]
+    go :: [PT] -> [(PT, VCategory)]
     go (x':y':z':xs) = classify x' y' z' : go (y':z':xs)
     go _          = []
 classifyList _ = []
@@ -32,10 +31,10 @@ classifyList _ = []
 
 -- |Classify a vertex on a polygon given it's next and previous vertex
 -- into five categories (see VCategory).
-classify :: P2  -- ^ prev vertex
+classify :: PT  -- ^ prev vertex
-         -> P2  -- ^ classify this one
+         -> PT  -- ^ classify this one
-         -> P2  -- ^ next vertex
+         -> PT  -- ^ next vertex
-         -> (P2, VCategory)
+         -> (PT, VCategory)
 classify prev v next
   | isVStart prev v next = (v, VStart)
   | isVSplit prev v next = (v, VSplit)
@@ -46,9 +45,9 @@ classify prev v next
 
 -- |Whether the vertex, given it's next and previous vertex,
 -- is a start vertex.
-isVStart :: P2  -- ^ previous vertex
+isVStart :: PT  -- ^ previous vertex
-         -> P2  -- ^ vertice to check
+         -> PT  -- ^ vertice to check
-         -> P2  -- ^ next vertex
+         -> PT  -- ^ next vertex
          -> Bool
 isVStart prev v next =
   ptCmpY next v == LT && ptCmpY prev v == LT && cw next v prev
@@ -56,9 +55,9 @@ isVStart prev v next =
 
 -- |Whether the vertex, given it's next and previous vertex,
 -- is a split vertex.
-isVSplit :: P2  -- ^ previous vertex
+isVSplit :: PT  -- ^ previous vertex
-         -> P2  -- ^ vertice to check
+         -> PT  -- ^ vertice to check
-         -> P2  -- ^ next vertex
+         -> PT  -- ^ next vertex
          -> Bool
 isVSplit prev v next =
   ptCmpY prev v == LT && ptCmpY next v == LT && cw prev v next
@@ -66,9 +65,9 @@ isVSplit prev v next =
 
 -- |Whether the vertex, given it's next and previous vertex,
 -- is an end vertex.
-isVEnd :: P2  -- ^ previous vertex
+isVEnd :: PT  -- ^ previous vertex
-       -> P2  -- ^ vertice to check
+       -> PT  -- ^ vertice to check
-       -> P2  -- ^ next vertex
+       -> PT  -- ^ next vertex
        -> Bool
 isVEnd prev v next =
   ptCmpY prev v == GT && ptCmpY next v == GT && cw next v prev
@@ -76,9 +75,9 @@ isVEnd prev v next =
 
 -- |Whether the vertex, given it's next and previous vertex,
 -- is a merge vertex.
-isVMerge :: P2  -- ^ previous vertex
+isVMerge :: PT  -- ^ previous vertex
-         -> P2  -- ^ vertice to check
+         -> PT  -- ^ vertice to check
-         -> P2  -- ^ next vertex
+         -> PT  -- ^ next vertex
          -> Bool
 isVMerge prev v next =
   ptCmpY next v == GT && ptCmpY prev v == GT && cw prev v next
@@ -86,9 +85,9 @@ isVMerge prev v next =
 
 -- |Whether the vertex, given it's next and previous vertex,
 -- is a regular vertex.
-isVRegular :: P2  -- ^ previous vertex
+isVRegular :: PT  -- ^ previous vertex
-           -> P2  -- ^ vertice to check
+           -> PT  -- ^ vertice to check
-           -> P2  -- ^ next vertex
+           -> PT  -- ^ next vertex
            -> Bool
 isVRegular prev v next =
      (not . isVStart prev v $ next)
@@ -99,7 +98,7 @@ isVRegular prev v next =
 
 
 -- |A polygon P is y-monotone, if it has no split and merge vertices.
-isYmonotone :: [P2] -> Bool
+isYmonotone :: [PT] -> Bool
 isYmonotone poly =
   not
   . any (\x -> x == VSplit || x == VMerge)
@@ -108,12 +107,12 @@ isYmonotone poly =
 
 
 -- |Partition P into y-monotone pieces.
-monotonePartitioning :: [P2] -> [[P2]]
+monotonePartitioning :: [PT] -> [[PT]]
 monotonePartitioning pts
   | isYmonotone pts = [pts]
   | otherwise = go (monotoneDiagonals pts) pts
   where
-    go :: [(P2, P2)] -> [P2] -> [[P2]]
+    go :: [Segment] -> [PT] -> [[PT]]
     go (x:xs) pts'@(_:_)
       | isYmonotone a && isYmonotone b = [a, b]
       | isYmonotone b = b : go xs a
@@ -125,37 +124,37 @@ monotonePartitioning pts
 
 -- |Try to eliminate the merge and split vertices by computing the
 -- diagonals we have to use for splitting the polygon.
-monotoneDiagonals :: [P2] -> [(P2, P2)]
+monotoneDiagonals :: [PT] -> [Segment]
 monotoneDiagonals pts = catMaybes . go $ classifyList pts
   where
-    go :: [(P2, VCategory)] -> [Maybe (P2, P2)]
+    go :: [(PT, VCategory)] -> [Maybe Segment]
     go (x:xs) = case snd x of
       VMerge -> getSeg (belowS . fst $ x) (fst x) : go xs
       VSplit -> getSeg (aboveS . fst $ x) (fst x) : go xs
       _      -> [] ++ go xs
     go [] = []
-    getSeg :: [P2] -- all points above/below the current point
+    getSeg :: [PT] -- all points above/below the current point
-           -> P2   -- current point
+           -> PT   -- current point
-           -> Maybe (P2, P2)
+           -> Maybe Segment
     getSeg [] _ = Nothing
     getSeg (z:zs) pt
       | isInsidePoly pts (z, pt) = Just (z, pt)
       | otherwise = getSeg zs pt
-    aboveS :: P2 -> [P2]
+    aboveS :: PT -> [PT]
     aboveS pt = tail . dropWhile (/= pt) $ sortedYX pts
-    belowS :: P2 -> [P2]
+    belowS :: PT -> [PT]
     belowS pt = reverse . takeWhile (/= pt) $ sortedYX pts
 
 
 -- |Triangulate a y-monotone polygon.
-triangulate :: [P2] -> [[P2]]
+triangulate :: [PT] -> [[PT]]
 triangulate pts =
   go pts . A.first reverse . splitAt 3 . reverse . sortedYX $ pts
   where
-    go :: [P2]          -- current polygon
+    go :: [PT]          -- current polygon
-       -> ([P2], [P2])  -- (stack of visited vertices, rest)
+       -> ([PT], [PT])  -- (stack of visited vertices, rest)
                         --  sorted by Y-coordinate
-       -> [[P2]]
+       -> [[PT]]
     go xs (p@[_, _], r:rs) = go xs (r:p, rs)
     go xs (p@(u:vi:vi1:ys), rs)
       -- case 1 and 3

Algorithms/QuadTree.hs

@@ -56,8 +56,7 @@ data Orient = North | South | East | West
 
 
 -- |Get a sub-square of the current square, e.g. nw, ne, sw or se.
-nwSq, neSq, swSq, seSq :: ((Double, Double), (Double, Double)) -- ^ current square
+nwSq, neSq, swSq, seSq :: Square -> Square
-                       -> ((Double, Double), (Double, Double)) -- ^ sub-square
 nwSq ((xl, yl), (xu, yu)) = (,)  (xl, (yl + yu) / 2)  ((xl + xu) / 2, yu)
 neSq ((xl, yl), (xu, yu)) = (,)  ((xl + xu) / 2, (yl + yu) / 2)  (xu, yu)
 swSq ((xl, yl), (xu, yu)) = (,)  (xl, yl)  ((xl + xu) / 2, (yl + yu) / 2)
@@ -80,9 +79,9 @@ isSEchild _                 = False
 -- |Builds a quadtree of a list of points which recursively divides up 2D
 -- space into quadrants, so that every leaf-quadrant stores either zero or one
 -- point.
-quadTree :: [P2]                                 -- ^ the points to divide
+quadTree :: [PT]        -- ^ the points to divide
-         -> ((Double, Double), (Double, Double)) -- ^ the initial square around the points
+         -> Square      -- ^ the initial square around the points
-         -> QuadTree P2                          -- ^ the quad tree
+         -> QuadTree PT -- ^ the quad tree
 quadTree []   _ = TNil
 quadTree [pt] _ = TLeaf pt
 quadTree pts sq = TNode (quadTree nWPT . nwSq $ sq) (quadTree nEPT . neSq $ sq)
@@ -96,9 +95,9 @@ quadTree pts sq = TNode (quadTree nWPT . nwSq $ sq) (quadTree nEPT . neSq $ sq)
 
 
 -- |Get all squares of a quad tree.
-quadTreeSquares :: ((Double, Double), (Double, Double))  -- ^ the initial square around the points
+quadTreeSquares :: Square      -- ^ the initial square around the points
-                -> QuadTree P2                            -- ^ the quad tree
+                -> QuadTree PT -- ^ the quad tree
-                -> [((Double, Double), (Double, Double))] -- ^ all squares of the quad tree
+                -> [Square]    -- ^ all squares of the quad tree
 quadTreeSquares sq (TNil)              = [sq]
 quadTreeSquares sq (TLeaf _)           = [sq]
 quadTreeSquares sq (TNode nw ne sw se) =
@@ -108,9 +107,7 @@ quadTreeSquares sq (TNode nw ne sw se) =
 
 -- |Get the current square of the zipper, relative to the given top
 -- square.
-getSquareByZipper :: ((Double, Double), (Double, Double)) -- ^ top square
+getSquareByZipper :: Square -> QTZipper a -> Square
-                  -> QTZipper a
-                  -> ((Double, Double), (Double, Double)) -- ^ current square
 getSquareByZipper sq z = go sq (reverse . snd $ z)
   where
     go sq' []              = sq'
@@ -203,7 +200,7 @@ lookupByNeighbors :: [Orient] -> QTZipper a -> Maybe (QTZipper a)
 lookupByNeighbors = flip (foldlM (flip findNeighbor))
 
 
-quadTreeToRoseTree :: QTZipper P2 -> Tree String
+quadTreeToRoseTree :: QTZipper PT -> Tree String
 quadTreeToRoseTree z' = go (rootNode z')
   where
     go z = case z of

CG2.cabal

@@ -65,6 +65,7 @@ executable Gtk
                        Graphics.Diagram.Core
                        Graphics.Diagram.Gtk
                        Graphics.Diagram.Plotter
+                       Graphics.HalfEdge
                        GUI.Gtk
                        MyPrelude
                        Parser.Meshparser
@@ -76,21 +77,20 @@ executable Gtk
 
   -- Other library packages from which modules are imported.
   build-depends:       attoparsec >= 0.12.1.1,
-                       base >=4.6 && <4.8,
+                       base >=4.6,
                        bytestring >= 0.10.4.0,
                        containers >= 0.5.0.0,
                        dequeue >= 0.1.5,
-                       diagrams-lib >=1.2 && <1.3,
+                       diagrams-lib >=1.3,
-                       diagrams-cairo >=1.2 && <1.3,
+                       diagrams-cairo >=1.3,
-                       diagrams-contrib >= 1.1.2.1,
+                       diagrams-contrib >= 1.3.0.0,
-                       directory >=1.2 && <1.3,
+                       directory >=1.2,
                        filepath >= 1.3.0.2,
-                       glade >=0.12 && <0.13,
+                       glade >=0.12,
                        gloss >= 1.2.0.1,
-                       gtk >=0.12 && <0.13,
+                       gtk >=0.12,
-                       multiset-comb >= 0.2.1,
                        safe >= 0.3.8,
-                       transformers >=0.4 && <0.5
+                       transformers >=0.4
 
   -- Directories containing source files.
   -- hs-source-dirs:
@@ -115,6 +115,7 @@ executable Gif
                        Graphics.Diagram.Core
                        Graphics.Diagram.Gif
                        Graphics.Diagram.Plotter
+                       Graphics.HalfEdge
                        MyPrelude
                        Parser.Meshparser
                        Parser.PathParser
@@ -126,18 +127,17 @@ executable Gif
 
   -- Other library packages from which modules are imported.
   build-depends:       attoparsec >= 0.12.1.1,
-                       base >=4.6 && <4.8,
+                       base >=4.6,
                        bytestring >= 0.10.4.0,
                        containers >= 0.5.0.0,
                        dequeue >= 0.1.5,
-                       diagrams-lib >=1.2 && <1.3,
+                       diagrams-lib >=1.3,
-                       diagrams-cairo >=1.2 && <1.3,
+                       diagrams-cairo >=1.3,
-                       diagrams-contrib >= 1.1.2.1,
+                       diagrams-contrib >= 1.3.0.0,
                        gloss >= 1.2.0.1,
                        JuicyPixels >= 3.1.7.1,
-                       multiset-comb >= 0.2.1,
+                       safe >= 0.3.8,
-                       transformers >=0.4 && <0.5,
+                       transformers >=0.4
-                       safe >= 0.3.8
 
   -- Directories containing source files.
   -- hs-source-dirs:
@@ -162,6 +162,7 @@ executable Test
                        Graphics.Diagram.Core
                        Graphics.Diagram.Gif
                        Graphics.Diagram.Plotter
+                       Graphics.HalfEdge
                        MyPrelude
                        Parser.Meshparser
                        Parser.PathParser
@@ -175,18 +176,14 @@ executable Test
 
   -- Other library packages from which modules are imported.
   build-depends:       attoparsec >= 0.12.1.1,
-                       base >=4.6 && <4.8,
+                       base >=4.6,
                        bytestring >= 0.10.4.0,
                        containers >= 0.5.0.0,
-                       dequeue >= 0.1.5,
+                       diagrams-lib >=1.3,
-                       diagrams-lib >=1.2 && <1.3,
+                       diagrams-cairo >=1.3,
-                       diagrams-cairo >=1.2 && <1.3,
+                       diagrams-contrib >= 1.3.0.0,
-                       diagrams-contrib >= 1.1.2.1,
                        gloss >= 1.2.0.1,
-                       JuicyPixels >= 3.1.7.1,
-                       multiset-comb >= 0.2.1,
                        QuickCheck >= 2.4.2,
-                       transformers >=0.4 && <0.5,
                        safe >= 0.3.8
 
   -- Directories containing source files.

GUI/Gtk.hs

@@ -63,9 +63,9 @@ data MyGUI = MkMyGUI {
   -- |Path entry widget for the quad tree.
   quadPathEntry :: Entry,
   -- |Horizontal box containing the path entry widget.
-  vbox7 :: Box,
+  vbox7 :: Graphics.UI.Gtk.Box,
   -- |Horizontal box containing the Rang search entry widgets.
-  vbox10 :: Box,
+  vbox10 :: Graphics.UI.Gtk.Box,
   -- |Range entry widget for lower x bound
   rangeXminEntry :: Entry,
   -- |Range entry widget for upper x bound
@@ -299,9 +299,9 @@ saveAndDrawDiag fp fps mygui =
         renderDiag winWidth winHeight buildDiag =
           renderDia Cairo
             (CairoOptions fps
-               (Dims (fromIntegral winWidth) (fromIntegral winHeight))
+               (mkSizeSpec2D (Just $ fromIntegral winWidth) (Just $ fromIntegral winHeight))
                SVG False)
-            (buildDiag (def{
+            (buildDiag (MyPrelude.def{
                 dotSize       = scaleVal,
                 xDimension    = fromMaybe (0, 500) xDim,
                 yDimension    = fromMaybe (0, 500) yDim,

Graphics/Diagram/AlgoDiags.hs

@@ -2,6 +2,7 @@
 
 module Graphics.Diagram.AlgoDiags where
 
+import Algebra.Vector(PT,Square)
 import Algorithms.GrahamScan
 import Algorithms.QuadTree
 import Algorithms.KDTree
@@ -123,9 +124,7 @@ kdSquares = Diag f
       where
         -- Gets all lines that make up the kdSquares. Every line is
         -- described by two points, start and end respectively.
-        kdLines :: KDTree P2
+        kdLines :: KDTree PT -> Square -> [(PT, PT)]
-                -> ((Double, Double), (Double, Double)) -- ^ square
-                -> [(P2, P2)]
         kdLines (KTNode ln pt Horizontal rn) ((xmin, ymin), (xmax, ymax)) =
           (\(x, _) -> [(p2 (x, ymin), p2 (x, ymax))])
             (unp2 pt)
@@ -180,7 +179,7 @@ kdTreeDiag = Diag f
 
 
 -- |Get the quad tree corresponding to the given points and diagram properties.
-qt :: [P2] -> DiagProp -> QuadTree P2
+qt :: [PT] -> DiagProp -> QuadTree PT
 qt vt p = quadTree vt (diagDimSquare p)
 
 
@@ -193,9 +192,7 @@ quadPathSquare = Diag f
       (uncurry rectByDiagonal # lw thin # lc red)
       (getSquare (stringToQuads (quadPath p)) (qt (mconcat vts) p, []))
       where
-        getSquare :: [Either Quad Orient]
+        getSquare :: [Either Quad Orient] -> QTZipper PT -> Square
-                  -> QTZipper P2
-                  -> ((Double, Double), (Double, Double))
         getSquare [] z = getSquareByZipper (diagDimSquare p) z
         getSquare (q:qs) z = case q of
           Right x -> getSquare qs (fromMaybe z (findNeighbor x z))
@@ -211,9 +208,7 @@ gifQuadPath = GifDiag f
       (uncurry rectByDiagonal # lw thick # lc col)
       <$> getSquares (stringToQuads (quadPath p)) (qt vt p, [])
       where
-        getSquares :: [Either Quad Orient]
+        getSquares :: [Either Quad Orient] -> QTZipper PT -> [Square]
-                   -> QTZipper P2
-                   -> [((Double, Double), (Double, Double))]
         getSquares [] z = [getSquareByZipper (diagDimSquare p) z]
         getSquares (q:qs) z = case q of
           Right x -> getSquareByZipper (diagDimSquare p) z :
@@ -233,12 +228,12 @@ treePretty = Diag f
                       . quadPath
                       $ p)
       where
-        getCurQT :: [Either Quad Orient] -> QTZipper P2 -> QTZipper P2
+        getCurQT :: [Either Quad Orient] -> QTZipper PT -> QTZipper PT
         getCurQT [] z = z
         getCurQT (q:qs) z = case q of
           Right x -> getCurQT qs (fromMaybe z (findNeighbor x z))
           Left x  -> getCurQT qs (fromMaybe z (goQuad x z))
-        prettyRoseTree :: Tree String -> Diagram Cairo R2
+        prettyRoseTree :: Tree String -> Diagram Cairo
         prettyRoseTree tree =
         -- HACK: in order to give specific nodes a specific color
           renderTree (\n -> case head n of

Graphics/Diagram/Core.hs

@@ -15,18 +15,18 @@ data Diag =
   Diag
   {
     mkDiag :: DiagProp
-           -> [[P2]]
+           -> [[PT]]
-           -> Diagram Cairo R2
+           -> Diagram Cairo
   }
   | GifDiag
   {
     mkGifDiag :: DiagProp
               -> Colour Double
-              -> ([P2] -> [[P2]])
+              -> ([PT] -> [[PT]])
-              -> [P2]
+              -> [PT]
-              -> [Diagram Cairo R2]
+              -> [Diagram Cairo]
   }
-  | EmptyDiag (Diagram Cairo R2)
+  | EmptyDiag (Diagram Cairo)
 
 
 -- |Holds the properties for a Diagram, like thickness of 2d points etc.
@@ -49,7 +49,7 @@ data DiagProp = MkProp {
   -- |The path to a quad in the quad tree.
   quadPath :: String,
   -- |The square for the kd-tree range search.
-  rangeSquare :: ((Double, Double), (Double, Double))
+  rangeSquare :: Square
 }
 
 
@@ -134,33 +134,33 @@ maybeDiag b d
   | otherwise = mempty
 
 
-filterValidPT :: DiagProp -> [P2] -> [P2]
+filterValidPT :: DiagProp -> [PT] -> [PT]
 filterValidPT =
   filter
   . inRange
   . diagDimSquare
 
 
-diagDimSquare :: DiagProp -> ((Double, Double), (Double, Double))
+diagDimSquare :: DiagProp -> Square
 diagDimSquare p = dimToSquare (xDimension p) $ yDimension p
 
 
 -- |Draw a list of points.
-drawP :: [P2]              -- ^ the points to draw
+drawP :: [PT]              -- ^ the points to draw
       -> Double            -- ^ dot size
-      -> Diagram Cairo R2  -- ^ the resulting diagram
+      -> Diagram Cairo     -- ^ the resulting diagram
 drawP [] _  = mempty
 drawP vt ds =
   position (zip vt (repeat dot))
   where
-    dot = circle ds :: Diagram Cairo R2
+    dot = circle ds :: Diagram Cairo
 
 
 -- |Create a rectangle around a diagonal line, which has sw
 -- as startpoint and nw as endpoint.
 rectByDiagonal :: (Double, Double)  -- ^ sw point
                -> (Double, Double)  -- ^ nw point
-               -> Diagram Cairo R2
+               -> Diagram Cairo
 rectByDiagonal (xmin, ymin) (xmax, ymax) =
   fromVertices [p2 (xmin, ymin)
                 , p2 (xmax, ymin)
@@ -172,7 +172,7 @@ rectByDiagonal (xmin, ymin) (xmax, ymax) =
 
 -- |Creates a Diagram from a point that shows the coordinates
 -- in text format, such as "(1.0, 2.0)".
-pointToTextCoord :: P2 -> Diagram Cairo R2
+pointToTextCoord :: PT -> Diagram Cairo
 pointToTextCoord pt =
   text ("(" ++ (show . trim') x ++ ", " ++ (show . trim') y ++ ")") # scale 10
   where

Graphics/Diagram/Gif.hs

@@ -16,7 +16,7 @@ import Parser.Meshparser
 
 
 -- |Return a list of tuples used by 'gifMain' to generate an animated gif.
-gifDiag :: DiagProp -> [PT] -> [(Diagram Cairo R2, GifDelay)]
+gifDiag :: DiagProp -> [PT] -> [(Diagram Cairo, GifDelay)]
 gifDiag p xs =
   fmap ((\x -> (x, 50)) . (<> nonChDiag))
    (upperHullList
@@ -35,5 +35,5 @@ gifDiag p xs =
 
 -- |Same as gifDiag, except that it takes a string containing the
 -- mesh file content instead of the the points.
-gifDiagS :: DiagProp -> B.ByteString -> [(Diagram Cairo R2, GifDelay)]
+gifDiagS :: DiagProp -> B.ByteString -> [(Diagram Cairo, GifDelay)]
-gifDiagS p = gifDiag p . filterValidPT p . meshToArr
+gifDiagS p = gifDiag p . filterValidPT p . meshVertices

Graphics/Diagram/Gtk.hs

@@ -2,6 +2,7 @@
 
 module Graphics.Diagram.Gtk where
 
+import Algebra.Vector(PT)
 import qualified Data.ByteString.Char8 as B
 import Data.List(find)
 import Diagrams.Backend.Cairo
@@ -45,7 +46,7 @@ diagTreAlgos =
 
 
 -- |Create the Diagram from the points.
-diag :: DiagProp -> [DiagAlgo] -> [[P2]] -> Diagram Cairo R2
+diag :: DiagProp -> [DiagAlgo] -> [[PT]] -> Diagram Cairo
 diag p das vts = maybe mempty (\x -> mkDiag x p vts)
                   $ mconcat
                       -- get the actual [Diag] array
@@ -57,22 +58,22 @@ diag p das vts = maybe mempty (\x -> mkDiag x p vts)
 
 -- |Create the Diagram from a String which is supposed to be the contents
 -- of an obj file.
-diagS :: DiagProp -> B.ByteString -> Diagram Cairo R2
+diagS :: DiagProp -> B.ByteString -> Diagram Cairo
 diagS p mesh =
   diag p diagAlgos
     . fmap (filterValidPT p)
-    . (\x -> if null x then [meshToArr mesh] else x)
+    . (\x -> if null x then [meshVertices mesh] else x)
-    . facesToArr
+    . meshFaceVertices
     $ mesh
 
 
 -- |Create the tree diagram from a String which is supposed to be the contents
 -- of an obj file.
-diagTreeS :: DiagProp -> B.ByteString -> Diagram Cairo R2
+diagTreeS :: DiagProp -> B.ByteString -> Diagram Cairo
 diagTreeS p mesh =
   diag p diagTreAlgos
     . fmap (filterValidPT p)
-    . (\x -> if null x then [meshToArr mesh] else x)
+    . (\x -> if null x then [meshVertices mesh] else x)
-    . facesToArr
+    . meshFaceVertices
     $ mesh
 

Graphics/HalfEdge.hs

@@ -0,0 +1,241 @@
+{-# OPTIONS_HADDOCK ignore-exports #-}
+
+-- |This module provides methods to build a cyclic half-edge data structure
+-- from an already parsed obj mesh file. As such, it depends on details
+-- of the parsed data.
+--
+-- In particular, 'indirectHeFaces', 'indirectHeVerts' and 'indirectToDirect'
+-- assume specific structure of some input lists. Check their respective
+-- documentation.
+--
+-- As the data structure has a lot of cross-references and the knots are
+-- not really known at compile-time, we have to use helper data structures
+-- such as lists and maps under the hood and tie the knots through
+-- index lookups.
+--
+-- For an explanation of the abstract concept of the half-edge data structure,
+-- check <http://www.flipcode.com/archives/The_Half-Edge_Data_Structure.shtml>
+module Graphics.HalfEdge (
+    HeVert(..)
+  , HeFace(..)
+  , HeEdge(..)
+  , buildHeEdge
+  , buildHeEdgeFromStr
+) where
+
+import Algebra.Vector
+import Control.Applicative
+import Control.Monad
+import qualified Data.ByteString.Char8 as B
+import qualified Data.IntMap.Lazy as Map
+import Data.Maybe
+import Parser.Meshparser
+import Safe
+
+
+-- |The vertex data structure for the half-edge.
+data HeVert a = HeVert {
+    vcoord  :: a        -- the coordinates of the vertex
+  , emedge  :: HeEdge a -- one of the half-edges emanating from the vertex
+} | NoVert
+
+
+-- |The face data structure for the half-edge.
+data HeFace a = HeFace {
+  bordedge :: HeEdge a -- one of the half-edges bordering the face
+} | NoFace
+
+-- |The actual half-edge data structure.
+data HeEdge a = HeEdge {
+    startvert :: HeVert a  -- start-vertex of the half-edge
+  , oppedge   :: HeEdge a  -- oppositely oriented adjacent half-edge
+  , edgeface  :: HeFace a  -- face the half-edge borders
+  , nextedge  :: HeEdge a  -- next half-edge around the face
+} | NoEdge
+
+-- This is a helper data structure of half-edge edges
+-- for tying the knots in 'indirectToDirect'.
+data IndirectHeEdge = IndirectHeEdge {
+    edgeindex  :: Int  -- edge index
+  , svindex    :: Int  -- index of start-vertice
+  , nvindex    :: Int  -- index of next-vertice
+  , indexf     :: Int  -- index of face
+  , offsetedge :: Int  -- offset to get the next edge
+} deriving (Show)
+
+-- This is a helper data structure of half-edge vertices
+-- for tying the knots in 'indirectToDirect'.
+data IndirectHeVert = IndirectHeVert {
+    emedgeindex  :: Int    -- emanating edge index (starts at 1)
+  , edgelist     :: [Int]  -- index of edge that points to this vertice
+} deriving (Show)
+
+-- This is a helper data structure of half-edge faces
+-- for tying the knots in 'indirectToDirect'.
+data IndirectHeFace =
+  IndirectHeFace (Int, [Int]) -- (faceIndex, [verticeindex])
+  deriving (Show)
+
+
+-- |Construct the indirect data structure for half-edge faces.
+-- This function assumes that the input faces are parsed exactly like so:
+--
+-- @
+-- f 1 3 4 5
+-- f 4 6 1 3
+-- @
+--
+-- becomes
+--
+-- > [[1,3,4,5],[4,6,1,3]]
+indirectHeFaces :: [[Int]]          -- ^ list of faces with their respective
+                                    --   list of vertice-indices
+                -> [IndirectHeFace]
+indirectHeFaces = fmap IndirectHeFace . zip [0..]
+
+
+-- |Construct the indirect data structure for half-edge edges.
+indirectHeEdges :: [IndirectHeFace] -> [IndirectHeEdge]
+indirectHeEdges = concat . fmap indirectHeEdge
+  where
+    indirectHeEdge :: IndirectHeFace -> [IndirectHeEdge]
+    indirectHeEdge   (IndirectHeFace (_, []))  = []
+    indirectHeEdge p@(IndirectHeFace (_, pv@(v:_))) = go p 0
+      where
+        go (IndirectHeFace (_, [])) _
+          = []
+        -- connect last to first element
+        go (IndirectHeFace (fi, [vlast])) ei
+          = [IndirectHeEdge ei vlast v fi (negate $ length pv - 1)]
+        -- regular non-last element
+        go (IndirectHeFace (fi, vfirst:vnext:vrest)) ei
+          = (:) (IndirectHeEdge ei vfirst vnext fi 1)
+                (go (IndirectHeFace (fi, vnext:vrest)) (ei + 1))
+
+
+-- |Construct the indirect data structure for half-edge vertices.
+-- It is assumed that the list of points is indexed in order of their
+-- appearance in the obj mesh file.
+indirectHeVerts :: [IndirectHeEdge]          -- ^ list of indirect edges
+                -> Map.IntMap IndirectHeVert -- ^ output map, starts at index 1
+indirectHeVerts hes' = go hes' Map.empty 0
+  where
+    go [] map' _ = map'
+    go (IndirectHeEdge _ _ nv _ offset:hes) map' i
+      = go hes
+           (Map.alter updateMap nv map')
+           (i + 1)
+      where
+        updateMap (Just (IndirectHeVert _ xs))
+          = Just (IndirectHeVert (i + offset) (i:xs))
+        updateMap Nothing
+          = Just (IndirectHeVert (i + offset) [i])
+
+
+-- |Tie the knots!
+-- It is assumed that the list of points is indexed in order of their
+-- appearance in the obj mesh file.
+--
+-- pseudo-code:
+--
+-- @
+-- indirectToDirect :: [a] -- parsed vertices, e.g. 2d points (Double, Double)
+--                  -> [IndirectHeEdge]
+--                  -> [IndirectHeFace]
+--                  -> [IndirectHeVert]
+--                  -> HeEdge a
+-- indirectToDirect points edges faces vertices
+--   = thisEdge (head edges)
+--   where
+--     thisEdge edge
+--       = HeEdge (thisVert (vertices !! svindex edge) $ svindex edge)
+--                (thisOppEdge (svindex edge) $ indexf edge)
+--                (thisFace $ faces !! indexf edge)
+--                (thisEdge $ edges !! (edgeindex edge + offsetedge edge))
+--     thisFace face = HeFace $ thisEdge (edges !! (head . snd $ face))
+--     thisVert vertice coordindex
+--       = HeVert (points !! (coordindex - 1))
+--                (thisEdge $ points !! (emedgeindex vertice - 1))
+--     thisOppEdge startverticeindex faceindex
+--       = case headMay
+--           . filter ((/=) faceindex . indexf)
+--           . fmap (edges !!)
+--           . edgelist         -- getter
+--           $ vertices !! startverticeindex
+--         of Just x  -> thisEdge x
+--            Nothing -> NoEdge
+-- @
+indirectToDirect :: [a]                       -- ^ list of points
+                 -> [IndirectHeEdge]
+                 -> [IndirectHeFace]
+                 -> Map.IntMap IndirectHeVert -- ^ assumed to start at index 1
+                 -> HeEdge a
+indirectToDirect pts pe@(e:_) fs vertmap
+  = thisEdge e
+  where
+    thisEdge (IndirectHeEdge ei sv _ fi off)
+      = case (fs `atMay` fi, pe `atMay` (ei + off), Map.lookup sv vertmap) of
+          (Just face,
+           Just edge,
+           Just vert) -> HeEdge (thisVert vert sv)
+                                (getOppEdge sv fi)
+                                (thisFace face)
+                                (thisEdge edge)
+          _           -> NoEdge
+    thisFace (IndirectHeFace (_, vi:_))
+      = case pe `atMay` vi of
+          Just edge -> HeFace (thisEdge edge)
+          Nothing   -> NoFace
+    thisFace (IndirectHeFace _) = NoFace
+    thisVert (IndirectHeVert eedg _) coordi
+      = case (pts `atMay` (coordi - 1), pe `atMay` (eedg - 1)) of
+          (Just vert, Just edge) -> HeVert vert $ thisEdge edge
+          _                      -> NoVert
+    getOppEdge sv fi
+      = case join
+               $ headMay
+               . filter ((/=) fi . indexf)
+               . catMaybes
+               . fmap (pe `atMay`)
+               . edgelist
+               <$> Map.lookup sv vertmap
+        of Just x  -> thisEdge x
+           Nothing -> NoEdge
+indirectToDirect _ _ _ _ = NoEdge
+
+
+-- |Build the half-edge data structure from a list of points
+-- and from a list of faces.
+-- The points are assumed to have been parsed in order of their appearance
+-- in the .obj mesh file, so that the indices match.
+-- The faces are assumed to have been parsed in order of their appearance
+-- in the .obj mesh file as follows:
+--
+-- @
+-- f 1 3 4 5
+-- f 4 6 1 3
+-- @
+--
+-- becomes
+--
+-- > [[1,3,4,5],[4,6,1,3]]
+buildHeEdge :: [a] -> [[Int]] -> Maybe (HeEdge a)
+buildHeEdge [] _   = Nothing
+buildHeEdge _ []   = Nothing
+buildHeEdge pts fs
+  = let faces' = indirectHeFaces fs
+        edges' = indirectHeEdges faces'
+        verts' = indirectHeVerts edges'
+    in Just $ indirectToDirect pts edges' faces' verts'
+
+
+-- |Build the HeEdge data structure from the .obj mesh file contents.
+buildHeEdgeFromStr :: B.ByteString -- ^ contents of an .obj mesh file
+                   -> HeEdge PT
+buildHeEdgeFromStr bmesh =
+  let pts    = meshVertices bmesh
+      faces' = indirectHeFaces . meshFaces $ bmesh
+      edges  = indirectHeEdges faces'
+      verts  = indirectHeVerts edges
+  in  indirectToDirect pts edges faces' verts
+

Parser/Meshparser.hs

@@ -1,7 +1,8 @@
 {-# OPTIONS_HADDOCK ignore-exports #-}
 
-module Parser.Meshparser (meshToArr, facesToArr) where
+module Parser.Meshparser where
 
+import Algebra.Vector(PT)
 import Control.Applicative
 import Data.Attoparsec.ByteString.Char8
 import Data.Either
@@ -11,19 +12,17 @@ import Diagrams.TwoD.Types
 
 -- |Convert a text String with multiple vertices and faces into
 -- a list of vertices, ordered by the faces specification.
-facesToArr :: B.ByteString -> [[P2]]
+meshFaceVertices :: B.ByteString -> [[PT]]
-facesToArr str = fmap (fmap (\y -> meshToArr str !! (fromIntegral y - 1)))
+meshFaceVertices str = fmap (fmap (\y -> meshVertices str !! (y - 1)))
-                      (faces str)
+                            (meshFaces str)
-  where
-    faces = rights . fmap (parseOnly parseFace) . B.lines
 
 
 -- |Convert a text String with multiple vertices into
 -- an array of float tuples.
-meshToArr :: B.ByteString -- ^ the string to convert
+meshVertices :: B.ByteString -- ^ the string to convert
-          -> [P2]         -- ^ the resulting vertice table
+          -> [PT]         -- ^ the resulting vertice table
-meshToArr =
+meshVertices
-  fmap p2
+  = fmap p2
     . rights
     . fmap (parseOnly parseVertice)
     . B.lines
@@ -37,5 +36,12 @@ parseVertice =
     <*> (many' space *> double)
 
 
-parseFace :: Parser [Integer]
+parseFace :: (Integral a) => Parser [a]
 parseFace = char 'f' *> many1' (many' space *> decimal)
+
+
+meshFaces :: B.ByteString -> [[Int]]
+meshFaces
+  = rights
+    . fmap (parseOnly parseFace)
+    . B.lines

Test/Vector.hs

@@ -21,19 +21,19 @@ newtype PosRoundDouble = PosRoundDouble { getPRD :: Double }
   deriving (Eq, Ord, Show, Read)
 
 
-newtype RoundR2 = RoundR2 { getRR2 :: R2 }
+newtype RoundR2 = RoundR2 { getRR2 :: V2 Double }
   deriving (Eq, Ord, Show, Read)
 
 
-newtype PosRoundR2 = PosRoundR2 { getPRR2 :: R2 }
+newtype PosRoundR2 = PosRoundR2 { getPRR2 :: V2 Double }
   deriving (Eq, Ord, Show, Read)
 
 
-newtype RoundP2 = RoundP2 { getRP2 :: P2 }
+newtype RoundP2 = RoundP2 { getRP2 :: P2 Double }
   deriving (Eq, Ord, Show, Read)
 
 
-newtype PosRoundP2 = PosRoundP2 { getPRP2 :: P2 }
+newtype PosRoundP2 = PosRoundP2 { getPRP2 :: P2 Double }
   deriving (Eq, Ord, Show, Read)
 
 
@@ -72,50 +72,50 @@ instance Arbitrary PosRoundP2 where
                 <*> (arbitrary :: Gen PosRoundDouble)
 
 
-instance Arbitrary R2 where
+instance Arbitrary (V2 Double) where
   arbitrary = curry r2 <$> arbitrary <*> arbitrary
 
 
-instance Arbitrary P2 where
+instance Arbitrary (P2 Double) where
   arbitrary = curry p2 <$> arbitrary <*> arbitrary
 
 
 -- the point describing the lower left corner of the square
 -- must be part of the square
-inRangeProp1 :: ((Double, Double), (Double, Double)) -> Bool
+inRangeProp1 :: Square -> Bool
 inRangeProp1 sq@((x1, y1), _) =
   inRange sq (p2 (x1, y1))
 
 
 -- the point describing the upper right corner of the square
 -- must be part of the square
-inRangeProp2 :: ((Double, Double), (Double, Double)) -> Bool
+inRangeProp2 :: Square -> Bool
 inRangeProp2 sq@(_, (x2, y2)) =
   inRange sq (p2 (x2, y2))
 
 
 -- the point describing the upper left corner of the square
 -- must be part of the square
-inRangeProp3 :: ((Double, Double), (Double, Double)) -> Bool
+inRangeProp3 :: Square -> Bool
 inRangeProp3 sq@((x1, _), (_, y2)) =
   inRange sq (p2 (x1, y2))
 
 
 -- the point describing the lower right corner of the square
 -- must be part of the square
-inRangeProp4 :: ((Double, Double), (Double, Double)) -> Bool
+inRangeProp4 :: Square -> Bool
 inRangeProp4 sq@((_, y1), (x2, _)) =
   inRange sq (p2 (x2, y1))
 
 
 -- generating random points within the square
-inRangeProp5 :: ((Double, Double), (Double, Double)) -> Positive Double -> Positive Double -> Bool
+inRangeProp5 :: Square -> Positive Double -> Positive Double -> Bool
 inRangeProp5 sq@((x1, y1), (x2, y2)) (Positive a) (Positive b) =
   inRange sq (p2 (x1 + ((x2 - x1) / (a + 1)), y1 + ((y2 - y1) / (b + 1))))
 
 
 -- generating random points outside of the square
-inRangeProp6 :: ((Double, Double), (Double, Double)) -> Positive Double -> Positive Double -> Bool
+inRangeProp6 :: Square -> Positive Double -> Positive Double -> Bool
 inRangeProp6 sq@((x1, y1), (x2, y2)) (Positive a) (Positive b) =
   (not . inRange sq $ p2 (max x1 x2 + (a + 1), max y1 y2 + (b + 1)))
   && (not . inRange sq $ p2 (max x1 x2 + (a + 1), max y1 y2 - (b + 1)))
@@ -126,51 +126,51 @@ inRangeProp6 sq@((x1, y1), (x2, y2)) (Positive a) (Positive b) =
 
 
 -- apply id function on the point
-onPTProp1 :: P2 -> Bool
+onPTProp1 :: PT -> Bool
 onPTProp1 pt = onPT id pt == pt
 
 
 -- add a random value to the point coordinates
-onPTProp2 :: P2 -> Positive R2 -> Bool
+onPTProp2 :: PT -> Positive (V2 Double) -> Bool
-onPTProp2 pt (Positive (R2 rx ry))
+onPTProp2 pt (Positive (V2 rx ry))
   = onPT (\(x, y) -> (x + rx, y + ry)) pt /= pt
 
 
 -- angle between two vectors both on the x-axis must be 0
-getAngleProp1 :: Positive R2 -> Positive R2 -> Bool
+getAngleProp1 :: Positive Vec -> Positive Vec -> Bool
-getAngleProp1 (Positive (R2 x1 _)) (Positive (R2 x2 _))
+getAngleProp1 (Positive (V2 x1 _)) (Positive (V2 x2 _))
-  = getAngle (R2 x1 0) (R2 x2 0) == 0
+  = getAngle (V2 x1 0) (V2 x2 0) == 0
 
 
 -- angle between two vectors both on the y-axis must be 0
-getAngleProp2 :: Positive R2 -> Positive R2 -> Bool
+getAngleProp2 :: Positive Vec -> Positive Vec -> Bool
-getAngleProp2 (Positive (R2 _ y1)) (Positive (R2 _ y2))
+getAngleProp2 (Positive (V2 _ y1)) (Positive (V2 _ y2))
-  = getAngle (R2 0 y1) (R2 0 y2) == 0
+  = getAngle (V2 0 y1) (V2 0 y2) == 0
 
 
 -- angle between two vectors both on the x-axis but with opposite direction
 -- must be pi
-getAngleProp3 :: Positive R2 -> Positive R2 -> Bool
+getAngleProp3 :: Positive Vec -> Positive Vec -> Bool
-getAngleProp3 (Positive (R2 x1 _)) (Positive (R2 x2 _))
+getAngleProp3 (Positive (V2 x1 _)) (Positive (V2 x2 _))
-  = getAngle (R2 (negate x1) 0) (R2 x2 0) == pi
+  = getAngle (V2 (negate x1) 0) (V2 x2 0) == pi
 
 
 -- angle between two vectors both on the y-axis but with opposite direction
 -- must be pi
-getAngleProp4 :: Positive R2 -> Positive R2 -> Bool
+getAngleProp4 :: Positive Vec -> Positive Vec -> Bool
-getAngleProp4 (Positive (R2 _ y1)) (Positive (R2 _ y2))
+getAngleProp4 (Positive (V2 _ y1)) (Positive (V2 _ y2))
-  = getAngle (R2 0 (negate y1)) (R2 0 y2) == pi
+  = getAngle (V2 0 (negate y1)) (V2 0 y2) == pi
 
 
 -- angle between vector in x-axis direction and y-axis direction must be
 -- p/2
-getAngleProp5 :: Positive R2 -> Positive R2 -> Bool
+getAngleProp5 :: Positive Vec -> Positive Vec -> Bool
-getAngleProp5 (Positive (R2 x1 _)) (Positive (R2 _ y2))
+getAngleProp5 (Positive (V2 x1 _)) (Positive (V2 _ y2))
-  = getAngle (R2 x1 0) (R2 0 y2) == pi / 2
+  = getAngle (V2 x1 0) (V2 0 y2) == pi / 2
 
 
 -- commutative
-getAngleProp6 :: Positive R2 -> Positive R2 -> Bool
+getAngleProp6 :: Positive Vec -> Positive Vec -> Bool
 getAngleProp6 (Positive v1) (Positive v2)
   = getAngle v1 v2 == getAngle v2 v1
 
@@ -183,7 +183,7 @@ getAngleProp7 (PosRoundR2 v)
 
 
 -- commutative
-scalarProdProp1 :: R2 -> R2 -> Bool
+scalarProdProp1 :: Vec -> Vec -> Bool
 scalarProdProp1 v1 v2 = v1 `scalarProd` v2 == v2 `scalarProd` v1
 
 
@@ -212,9 +212,9 @@ scalarProdProp4 (RoundDouble s1) (RoundDouble s2) (RoundR2 v1) (RoundR2 v2)
 
 
 -- orthogonal
-scalarProdProp5 :: Positive R2 -> Positive R2 -> Bool
+scalarProdProp5 :: Positive Vec -> Positive Vec -> Bool
-scalarProdProp5 (Positive (R2 x1 _)) (Positive (R2 _ y2))
+scalarProdProp5 (Positive (V2 x1 _)) (Positive (V2 _ y2))
-  = scalarProd (R2 x1 0) (R2 0 y2) == 0
+  = scalarProd (V2 x1 0) (V2 0 y2) == 0
 
 
 -- this is almost the same as the function definition
@@ -226,49 +226,49 @@ dimToSquareProp1 (x1, x2) (y1, y2) =
 -- multiply scalar with result of vecLength or with the vector itself...
 -- both results must be the same. We can't check against 0
 -- because of sqrt in vecLength.
-vecLengthProp1 :: PosRoundDouble -> R2 -> Bool
+vecLengthProp1 :: PosRoundDouble -> Vec -> Bool
 vecLengthProp1 (PosRoundDouble r) v
   = abs (vecLength v * r - vecLength (scalarMul r v)) < 0.0001
 
 
 -- convert to vector and back again
-pt2VecProp1 :: P2 -> Bool
+pt2VecProp1 :: PT -> Bool
 pt2VecProp1 pt = (vec2Pt . pt2Vec $ pt) == pt
 
 
 -- unbox coordinates and check if equal
-pt2VecProp2 :: P2 -> Bool
+pt2VecProp2 :: PT -> Bool
 pt2VecProp2 pt = (unr2 . pt2Vec $ pt) == unp2 pt
 
 
 -- convert to point and back again
-vec2PtProp1 :: R2 -> Bool
+vec2PtProp1 :: Vec -> Bool
 vec2PtProp1 v = (pt2Vec . vec2Pt $ v) == v
 
 
 -- unbox coordinates and check if equal
-vec2PtProp2 :: R2 -> Bool
+vec2PtProp2 :: Vec -> Bool
 vec2PtProp2 v = (unp2 . vec2Pt $ v) == unr2 v
 
 
 -- vector from a to b must not be the same as b to a
-vp2Prop1 :: P2 -> P2 -> Bool
+vp2Prop1 :: PT -> PT -> Bool
 vp2Prop1 p1' p2'
   | p1' == origin && p2' == origin = True
   | otherwise = vp2 p1' p2' /= vp2 p2' p1'
 
 
 -- negating vector from a to be must be the same as vector b to a
-vp2Prop2 :: P2 -> P2 -> Bool
+vp2Prop2 :: PT -> PT -> Bool
 vp2Prop2 p1' p2'
   | p1' == origin && p2' == origin = True
-  | otherwise = vp2 p1' p2' == (\(R2 x y) -> negate x ^& negate y)
+  | otherwise = vp2 p1' p2' == (\(V2 x y) -> negate x ^& negate y)
                                  (vp2 p2' p1')
                 &&
-                vp2 p2' p1' == (\(R2 x y) -> negate x ^& negate y)
+                vp2 p2' p1' == (\(V2 x y) -> negate x ^& negate y)
                                  (vp2 p1' p2')
 
 
 -- determinant of the 3 same points is always 0
-detProp1 :: P2 -> Bool
+detProp1 :: PT -> Bool
 detProp1 pt' = det pt' pt' pt' == 0

test_objs/testcube_trans.obj

@@ -0,0 +1,19 @@
+v 9.0 10.0
+v 11.0 10.0
+v 9.0 11.0
+v 11.0 11.0
+v 9.0 11.0
+v 11.0 11.0
+v 9.0 10.0
+v 11.0 10.0
+f 1 2 4 3
+f 3 4 6 5
+f 5 6 8 7
+f 7 8 2 1
+f 2 8 6 4
+f 7 1 3 5
+
+cstype bezier
+deg 3
+curv 1 2 3 4
+end