cubicbezier 0.4.0.2 → 0.5.0.0

raw patch · 11 files changed

+271/−182 lines, 11 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• Geom2D.hs +19/−0
• Geom2D/CubicBezier.hs +0/−2
• Geom2D/CubicBezier/Approximate.hs +3/−25
• Geom2D/CubicBezier/Basic.hs +15/−3
• Geom2D/CubicBezier/MetaPath.hs +11/−32
• Geom2D/CubicBezier/Numeric.hs +78/−9
• Geom2D/CubicBezier/Outline.hs +1/−1
• Geom2D/CubicBezier/Overlap.lhs +25/−26
• Math/BernsteinPoly.hs +18/−12
• cubicbezier.cabal +1/−1
• tests/test.hs +100/−71

@@ -146,6 +146,20 @@     py = (dy*u - dx*v) / m {-# specialize closestPoint :: Line Double -> Point Double -> Point Double #-}   +-- | Calculate the intersection of two lines.  If the determinant is+-- less than tolerance (parallel or coincident lines), return Nothing.+lineIntersect :: (Ord a, Floating a) => Line a -> Line a -> a -> Maybe (Point a)+lineIntersect (Line p1 p2) (Line p3 p4) eps+  | abs det <= eps = Nothing+  | otherwise = Just $ (a*^d2 ^-^ b*^d1) ^/ det+  where+    d1 = p1 ^-^ p2+    d2 = p3 ^-^ p4+    det = vectorCross d1 d2+    a = vectorCross p1 p2 +    b = vectorCross p3 p4+{-# SPECIALIZE lineIntersect :: Line Double -> Line Double -> Double -> Maybe DPoint #-}    + -- | The lenght of the vector. vectorMag :: Floating a => Point a -> a vectorMag (Point x y) = sqrt(x*x + y*y)@@ -251,3 +265,8 @@ translate :: Num a => Point a -> Transform a translate (Point x y) = Transform 1 0 x 0 1 y {-# INLINE translate #-}++-- | The identity transformation.+idTrans :: Num a => Transform a+idTrans = Transform 1 0 0 0 1 0+{-# INLINE idTrans #-}

@@ -4,7 +4,6 @@        (module Geom2D.CubicBezier.Basic,         module Geom2D.CubicBezier.Approximate,         module Geom2D.CubicBezier.Overlap,-        module Geom2D.CubicBezier.Outline,         module Geom2D.CubicBezier.Curvature,         module Geom2D.CubicBezier.Intersection,         module Geom2D.CubicBezier.MetaPath,@@ -14,7 +13,6 @@ import Geom2D import Geom2D.CubicBezier.Basic import Geom2D.CubicBezier.Approximate-import Geom2D.CubicBezier.Outline import Geom2D.CubicBezier.Curvature import Geom2D.CubicBezier.Overlap import Geom2D.CubicBezier.Intersection

@@ -75,35 +75,13 @@       nd = ((p ^-^ b) ^.^ b') / (b'^.^b')   in vectorDistance p $ evalBezier qb (t + nd) -phi :: (Floating a) => a-phi = (-1 + sqrt 5) / 2--goldSearch :: (Ord a, Floating a) => (a -> a) -> a-goldSearch f =-  goldSearch' f 0 x1 x2 1 (f 0)-  (f x1) (f x2) (f 1) 4-    where x1 = 1 - phi-          x2 = phi--goldSearch' :: (Ord a, Floating a) =>-               (a -> a) -> a -> a -> a ->-               a -> a -> a -> a -> a -> Int -> a-goldSearch' f x0 x1 x2 x3 y0 y1 y2 y3 maxiter-  | maxiter < 1 = maximum [y0, y1, y2, y3]-  | y1 < y2 =-    let x25 = x1 + phi*(x3-x1)-        y25 = f x25-    in goldSearch' f x1 x2 x25 x3 y1 y2 y25 y3 (maxiter-1)-  | otherwise =-    let x05 = x2 + phi*(x0-x2)-        y05 = f x05-    in goldSearch' f x0 x05 x1 x2 y0 y05 y1 y2 (maxiter-1)- -- find maximum distance using golden section search maxDist :: (V.Unbox a, Ord a, Floating a) =>            (a -> (Point a, Point a)) ->            QuadBezier a -> a -> a -> a-maxDist f qb tmin tmax = goldSearch (quadDist f qb tmin tmax)+maxDist f qb tmin tmax =+  quadDist f qb tmin tmax $+  goldSearch (quadDist f qb tmin tmax) 4  approxquad :: (Ord a, Floating a) =>               Point a -> Point a -> Point a -> Point a -> QuadBezier a

@@ -6,7 +6,7 @@         bezierParam, bezierParamTolerance, reorient, bezierToBernstein,         evalBezierDerivs, evalBezier, evalBezierDeriv, findBezierTangent, quadToCubic,         bezierHoriz, bezierVert, findBezierInflection, findBezierCusp,-        arcLength, arcLengthParam, splitBezier, bezierSubsegment, splitBezierN,+        bezierArc, arcLength, arcLengthParam, splitBezier, bezierSubsegment, splitBezierN,         colinear)        where import Geom2D@@ -22,14 +22,14 @@   cubicC1 :: !(Point a),   cubicC2 :: !(Point a),   cubicC3 :: !(Point a)}-                   deriving (Eq, Show, Functor)+  deriving (Eq, Show, Functor)  -- | A quadratic bezier curve. data QuadBezier a = QuadBezier {   quadC0 :: !(Point a),   quadC1 :: !(Point a),   quadC2 :: !(Point a)}-                  deriving (Eq, Show, Functor)+  deriving (Eq, Show, Functor)  -- Use a tuple, because it has 0(1) unzip when using unboxed vectors. -- | A bezier curve of any degree.@@ -296,6 +296,18 @@ findBezierCusp :: CubicBezier Double -> [Double] findBezierCusp b = filter vertical $ bezierHoriz b   where vertical = (== 0) . pointY . snd . evalBezierDeriv b++-- | @bezierArc startAngle endAngle@ approximates an arc on the unit circle with+-- a single cubic béziér curve.  Maximum deviation is <0.03% for arcs+-- 90° degrees or less.+bezierArc :: Double -> Double -> CubicBezier Double+bezierArc start end = CubicBezier p0 p1 p2 p3+  where+    p0 = dirVector start+    p3 = dirVector end+    p1 = p0 ^+^ k *^ (rotate90L $* p0)+    p2 = p1 ^+^ k *^ (rotate90R $* p3)+    k = 4/3 * tan((end-start)/4)  -- | @arcLength c t tol finds the arclength of the bezier c at t, within given tolerance tol. 

@@ -61,7 +61,8 @@ import Geom2D.CubicBezier.Basic import Data.List import Text.Printf-import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as V+import Geom2D.CubicBezier.Numeric  data OpenMetaPath a = OpenMetaPath [(Point a, MetaJoin a)] (Point a)                       -- ^ A metapath with endpoints@@ -209,7 +210,7 @@       tensionsA = map (tensionL . snd) nodes       tensionsB = map (tensionR . snd) nodes       turnAngles = zipWith turnAngle chords (tail $ cycle chords)-      thetas = solveCyclicTriD $+      thetas = solveCyclicTriD2 $                eqsCycle tensionsA                points                tensionsB@@ -234,7 +235,7 @@       tensionsA = map tensionL joins       tensionsB = map tensionR joins       turnAngles = zipWith turnAngle chords (tail chords) ++ [0]-      thetas = solveTriDiagonal $+      thetas = solveTriDiagonal2 $                eqsOpen points joins chords turnAngles                (map tensionValue tensionsA)                (map tensionValue tensionsB)@@ -335,40 +336,18 @@ -- where u[n] = 0 -- then solving for t[n] -- see metafont the program: ¶ 283-solveTriDiagonal :: [(Double, Double, Double, Double)] -> [Double]-solveTriDiagonal [] = error "solveTriDiagonal: not enough equations"-solveTriDiagonal ((_, b0, c0, d0): rows) =-  V.toList $ solveTriDiagonal2 (b0, c0, d0) (V.fromList rows)--solveTriDiagonal2 :: (Double, Double, Double) -> V.Vector (Double, Double, Double, Double) -> V.Vector Double-solveTriDiagonal2 (!b0, !c0, !d0) rows = solutions-  where-    twovars = V.scanl nextrow (c0/b0, d0/b0) rows-    solutions = V.scanr nextsol vn (V.unsafeInit twovars)-    vn = snd $ V.unsafeLast twovars-    nextsol (u, v) ti = v - u*ti-    nextrow (u, v) (ai, bi, ci, di) =-      (ci/(bi - u*ai), (di - v*ai)/(bi - u*ai))+solveTriDiagonal2 :: [(Double, Double, Double, Double)] -> [Double]+solveTriDiagonal2 [] = error "solveTriDiagonal: not enough equations"+solveTriDiagonal2 ((_, b0, c0, d0): rows) =+  V.toList $ solveTriDiagonal (b0, c0, d0) (V.fromList rows)  -- test = ((80.0,58.0,51.0),[(-432.0,78.0,102.0,503.0),(71.0,-82.0,20.0,2130.0),(52.39,-10.43,4.0,56.0),(34.0,38.0,0.0,257.0)]) -- [-15.726940528143576,22.571642107784243,-78.93751365259996,-297.27313545829384,272.74438435742667]        -- solve the cyclic tridiagonal system. -- see metafont the program: ¶ 286-solveCyclicTriD :: [(Double, Double, Double, Double)] -> [Double]-solveCyclicTriD rows = solutions-  where-    (!un, !vn, !wn): threevars =-      reverse $ tail $ scanl nextrow (0, 0, 1) rows-    nextrow (!u, !v, !w) (!ai, !bi, !ci, !di) =-      (ci/(bi - ai*u), (di - ai*v)/(bi - ai*u), -ai*w/(bi - ai*u))-    (totvn, totwn) = foldl (\(v', w') (u, v, w) ->-                             (v - u*v', w - u*w'))-                     (0, 1) threevars-    t0 = (vn - un*totvn) / (1 - (wn - un*totwn))-    solutions = scanl nextsol t0-                ((un, vn, wn) : reverse (tail threevars))-    nextsol t (!u, !v, !w) = (v + w*t0 - t)/u+solveCyclicTriD2 :: [(Double, Double, Double, Double)] -> [Double]+solveCyclicTriD2 = V.toList . solveCyclicTriD . V.fromList  turnAngle :: DPoint -> DPoint -> Double turnAngle (Point 0 0) _ = 0@@ -411,7 +390,7 @@        (0, 1, 0, turnAngle delta dir2)]     (Curl _, Curl _) ->       [(0, 1, 0, 0), (0, 1, 0, 0)]-    _ -> undefined+    _ -> error "illegal end of open path"  eqsOpen points joins chords turnAngles tensionsA tensionsB =   eq0 : restEquations joins tensionsA dists turnAngles tensionsB

@@ -1,6 +1,16 @@--- | Some numerical computations used by the cubic bezier functions-module Geom2D.CubicBezier.Numeric where-import Data.Vector.Unboxed as V+{-# LANGUAGE BangPatterns #-}+-- | Numerical computations used by the cubic bezier functions.  Also+-- contains functions that aren't used anymore, but might be useful on+-- its own.+module Geom2D.CubicBezier.Numeric+       (quadraticRoot, solveLinear2x2, +        goldSearch, +        makeSparse, SparseMatrix, sparseMulT, sparseMul,+        addMatrix, addVec, lsqMatrix, lsqSolveDist,+        decompLDL, lsqSolve,+        solveTriDiagonal, solveCyclicTriD)+       where+import qualified Data.Vector.Unboxed as V import Data.Vector.Unboxed.Mutable as MV import Data.Matrix.Unboxed as M import qualified Data.Matrix.Generic as G@@ -8,7 +18,7 @@ import Control.Monad.ST import Control.Monad -+sign :: (Ord a, Num a, Num t) => a -> t sign x | x < 0 = -1        | otherwise = 1 @@ -42,11 +52,39 @@   where det = d * b - a * e {-# SPECIALIZE solveLinear2x2 :: Double -> Double -> Double -> Double -> Double -> Double -> Maybe (Double, Double) #-} +phi :: (Floating a) => a+phi = (-1 + sqrt 5) / 2+++goldSearch :: (Ord a, Floating a) => (a -> a) -> Int -> a+goldSearch f maxiter =+  goldSearch' f 0 x1 x2 1 (f 0) +  (f x1) (f x2) (f 1) maxiter+    where x1 = 1 - phi+          x2 = phi+{-# SPECIALIZE goldSearch :: (Double -> Double) -> Int -> Double #-}          ++goldSearch' :: (Ord a, Floating a) =>+               (a -> a) -> a -> a -> a ->+               a -> a -> a -> a -> a -> Int -> a+goldSearch' f x0 x1 x2 x3 y0 y1 y2 y3 maxiter+  | maxiter < 1 = snd $ maximum [(y0, x0), (y1, x1),+                                 (y2, x2), (y3, x3)]+  | y1 < y2 =+    let x25 = x1 + phi*(x3-x1)+        y25 = f x25+    in goldSearch' f x1 x2 x25 x3 y1 y2 y25 y3 (maxiter-1)+  | otherwise =+    let x05 = x2 + phi*(x0-x2)+        y05 = f x05+    in goldSearch' f x0 x05 x1 x2 y0 y05 y1 y2 (maxiter-1)++ data SparseMatrix a =   SparseMatrix (V.Vector Int)   (V.Vector (Int, Int)) (M.Matrix a)                       -makeSparse :: Unbox a => Vector Int+makeSparse :: Unbox a => V.Vector Int               -- ^ The column index of the first element of each row.               -- Should be ascending in order.               -> M.Matrix a@@ -78,6 +116,7 @@       | v `V.unsafeIndex` e > i = e       | otherwise = nextEnd (e+1) i + -- | Given a rectangular matrix M, calculate the symmetric square -- matrix MᵀM which can be used to find a least squares solution to -- the overconstrained system.@@ -195,7 +234,7 @@         \i -> do           let vi = v `V.unsafeIndex` i           s <- liftM (V.foldl' (-) vi) $-               V.forM (enumFromN 0 i) $+               V.forM (V.enumFromN 0 i) $                \j -> liftM ((m `M.unsafeIndex` (j, i-j)) *)                      (MV.unsafeRead sol1 j)           MV.unsafeWrite sol1 i s@@ -205,7 +244,7 @@         \i -> do           let vi = v `V.unsafeIndex` i           s <- liftM (V.foldl' (-) vi) $-               V.forM (enumFromN 1 (width-1)) $+               V.forM (V.enumFromN 1 (width-1)) $                \j -> liftM ((m `M.unsafeIndex` (i-j, j)) *)                      (MV.unsafeRead sol1 $ i-j)           MV.unsafeWrite sol1 i s@@ -216,7 +255,7 @@           solI <- MV.unsafeRead sol1 (vars-i-1)           let d = m `M.unsafeIndex` (vars-i-1, 0)           s <- liftM (V.foldl' (-) (solI/d)) $-               V.forM (enumFromN 0 i) $+               V.forM (V.enumFromN 0 i) $                \j -> liftM ((m `M.unsafeIndex` (vars-i-1, j+1)) *)                      (MV.unsafeRead sol1 $ vars-i+j)           MV.unsafeWrite sol1 (vars-i-1) s@@ -227,7 +266,7 @@           solI <- MV.unsafeRead sol1 (vars-i-1)           let d = m `M.unsafeIndex` (vars-i-1, 0)           s <- liftM (V.foldl' (-) (solI/d)) $-               V.forM (enumFromN 0 (width-1)) $+               V.forM (V.enumFromN 0 (width-1)) $                \j -> liftM ((m `M.unsafeIndex` (vars-i-1, j+1)) *)                      (MV.unsafeRead sol1 $ vars-i+j)           MV.unsafeWrite sol1 (vars-i-1) s@@ -265,3 +304,33 @@     m1 = SparseMatrix r s m1'     m2 = SparseMatrix r s m2' {-# SPECIALIZE lsqSolveDist :: SparseMatrix (Double, Double) -> V.Vector (Double, Double) -> V.Vector Double #-}++-- | solve a tridiagonal system.  see metafont the program: ¶ 283+solveTriDiagonal :: (Unbox a, Fractional a) =>+                    (a, a, a) -> V.Vector (a, a, a, a) -> V.Vector a+solveTriDiagonal (!b0, !c0, !d0) rows = solutions+  where+    twovars = V.scanl nextrow (c0/b0, d0/b0) rows+    solutions = V.scanr nextsol vn (V.unsafeInit twovars)+    vn = snd $ V.unsafeLast twovars+    nextsol (u, v) ti = v - u*ti+    nextrow (u, v) (ai, bi, ci, di) =+      (ci/(bi - u*ai), (di - v*ai)/(bi - u*ai))+{-# SPECIALIZE solveTriDiagonal :: (Double, Double, Double) -> V.Vector (Double, Double, Double, Double) -> V.Vector Double #-}++-- | solve a cyclic tridiagonal system.  see metafont the program: ¶ 286+solveCyclicTriD :: (Unbox a, Fractional a) => V.Vector (a, a, a, a) -> V.Vector a+solveCyclicTriD rows = solutions+  where+    threevars = V.tail $ V.scanl nextrow (0, 0, 1) rows+    nextrow (!u, !v, !w) (!ai, !bi, !ci, !di) =+      (ci/(bi - ai*u), (di - ai*v)/(bi - ai*u), -ai*w/(bi - ai*u))+    (totvn, totwn) = V.foldr (\(u, v, w) (v', w') ->+                               (v - u*v', w - u*w'))+                     (0, 1) (V.init threevars)+    t0 = (vn - un*totvn) / (1 - (wn - un*totwn))+    (!un, !vn, !wn) = V.last threevars+    solutions = V.scanl nextsol t0 $ V.cons (un, vn, wn) $+                V.init $ V.init $ threevars+    nextsol t (!u, !v, !w) = (v + w*t0 - t)/u+{-# SPECIALIZE solveCyclicTriD :: V.Vector (Double, Double, Double, Double) -> V.Vector Double #-}

@@ -1,7 +1,7 @@ -- | Offsetting bezier curves and stroking curves.  module Geom2D.CubicBezier.Outline-       (bezierOffset)+       (bezierOffset, bezierOffsetPoint)        where import Geom2D import Geom2D.CubicBezier.Basic

@@ -3,21 +3,20 @@ Removing overlap from bezier paths in haskell ============================================= -This document describes an algorithm for removing overlap and-performing set operations on bezier paths, but at the same time it is-a working module for the haskell `cubicbezier` package.  This way it-can serve two purposes at once: someone who wants to implement this-algorithm can use this as an explanation of the algorithm, while at-the same time it is a working version of the described algorithm.+This literate programming document serves as both the description +and code for the algorithm for removing overlap and+performing set operations on bezier paths.  It can be used both for+understanding the code, and for porting the used algorithm to other +implementations. -**Note on porting**: Porting this code to another language should-present no difficulties.  However some care must be taken with regards-to lazyness.  Often many variables inside the `where` statement aren't+**Note on porting**: Porting should be fairly straightforward. +However some care must be taken with regards to lazyness.  Often+many variables inside the `where` statement aren't evaluated in all guards, so it's important to evaluate only those-which appear in the guards.  The main state can be modified using-mutation instead of copying without any troubles.  For modifying the-state, the `lens` library is used.  The lens functions can be-interpreted in a mutable language as follows:+which appear in the guards.  This library uses copying instead of +mutation, but for the most part mutation can be used as well.+It uses the `lens` library for modifying state.  The following translation+could be made to a mutable language:     * reading state:     - `view field struct`: `struct.field`@@ -54,11 +53,11 @@ > import qualified Data.Map.Strict as M > import qualified Data.Set as S -So what does it mean to remove overlap?  Basicly we want to keep-curves where one side is inside the filled region, and the other side-is outside, and discard the rest.Since that could be true only of a-part of the curve, we also need to split each curve when it intersects-another curve.  How do you know which side is the inside, and which+The basic idea is to keep curves where one side is inside the filled+region, and the other side is outside, and discard the rest. +Since that could be true only of a part of the curve, we also need to+split each curve when it intersects+another curve.  How to know which side is the inside, and which side the outside?  There are two methods which are use the most: the [*even-odd rule*](https://en.wikipedia.org/wiki/Even%E2%80%93odd_rule) and the [*nonzero rule*](https://en.wikipedia.org/wiki/Nonzero-rule).@@ -69,7 +68,7 @@ Checking each pair of curves for intersections would work, but is rather inefficient.  We only need to check for overlap when two curves are adjacent.  Fortunately there exist a good method from-*computational geometry*, called a *sweep line algorithm*.  The basic+*computational geometry*, called the *sweep line algorithm*.  The  idea is to sweep a vertical line over the input, starting from leftmost point to the right (of course the opposite direction is also possible), and to update the input dynamically.  We keep track of each@@ -87,9 +86,9 @@  > type XStruct = M.Map PointEvent [Curve] -Why `PointEvent`, and not just `Point`?  We need to have a `Ord`-instance for the map, which much match our horizontal ordering.  A-newtype is ideal, since it has no extra cost, and allows us to define+I use `PointEvent` instead of just `Point`.  This way I can have a `Ord`+instance for the map, which must match the horizontal ordering.  A+newtype is ideal, since it has no extra cost, and allows me to define a Ord instance for defining the relative order.  The value from the map is a list, since there can be many curves starting from the same point.@@ -97,7 +96,7 @@ > newtype PointEvent = PointEvent DPoint >                    deriving Show -When the x-coordinates are equal, I use the y-coordinate to determine+When the x-coordinates are equal, use the y-coordinate to determine the order.  > instance Eq PointEvent where@@ -109,13 +108,13 @@ >     compare (x1, y2) (x2, y1)  All curves are kept left to right, so we need to remember the-direction for the output:+curve direction for the output:  The curves intersecting the sweepline are kept in another balanced Tree, called the *Y-structure*.  *These curves are not allowed to overlap*, except in the endpoints, and will be ordered vertically.-I'll use the `Curve` datatype to define the ordering of the curves,-and to add additional information.  The `turnRatio` field is the+The `Curve` datatype defines the ordering of the curves,+and adds additional information.  The `turnRatio` field is the turnRatio of the area to the left for a left to right curve, and to the right for a right to left curve.  The `changeTurn` function determines how the turnRatio will change from up to down.  This

@@ -11,6 +11,7 @@         bernsteinEvalDeriv, binCoeff, convolve, bernsteinEvalDerivs, bernsteinDeriv)        where import Data.Vector.Unboxed as V+import Data.Vector.Unboxed.Mutable as M import qualified Data.Vector as B  data BernsteinPoly a = BernsteinPoly {@@ -60,17 +61,23 @@  -- | Calculate the convolution of two vectors. convolve :: (Unbox a, Num a) => Vector a -> Vector a -> Vector a-convolve a b =-  V.map (\i -> V.sum $-               V.zipWith (*) a $-               V.reverse $-               V.unsafeTake i b)-  (V.enumFromN 1 $ V.length b)-  V.++ V.map (\i -> V.sum $-                    V.zipWith (*)-                    (V.unsafeDrop i a)-                    (V.reverse b))-  (V.enumFromN 1 $ V.length a-1)+convolve x h = V.create $ do+  let xN = V.length x+  let hN = V.length h+  let xIndices = V.enumFromN 0 xN+  let hIndices = V.enumFromN 0 hN++  xM <- V.unsafeThaw x+  hM <- V.unsafeThaw h+  yM <- M.replicate (xN + hN - 1) 0++  V.forM_ xIndices $ \i -> do+    a <- M.unsafeRead xM i+    V.forM_ hIndices $ \j -> do+      b <- M.unsafeRead hM j+      M.unsafeModify yM (+ a * b) (i + j)++  return yM {-# SPECIALIZE convolve :: Vector Double -> Vector Double -> Vector Double #-}  -- | Multiply two bernstein polynomials using convolution.  The final@@ -114,7 +121,6 @@  -- | Evaluate the bernstein polynomial using the horner rule adapted -- for bernstein polynomials.- bernsteinEval :: (Unbox a, Fractional a)                  => BernsteinPoly a -> a -> a bernsteinEval (BernsteinPoly v) _

@@ -1,5 +1,5 @@ Name:		cubicbezier-Version: 	0.4.0.2+Version: 	0.5.0.0 Synopsis:	Efficient manipulating of 2D cubic bezier curves. Category: 	Graphics, Geometry, Typography Copyright: 	Kristof Bastiaensen (2014)

@@ -5,12 +5,7 @@ import Text.Parsec import Text.Parsec.String import Text.Parsec.Error-import MPTest-import NumTest -tests :: TestTree-tests = testGroup "Tests" [mfTests, numTests]- num :: Parser Double num =    liftM (read.concat) $ sequence@@ -19,7 +14,7 @@   , option "" $ string "."   , option "0" $ many digit] -pointP :: Parser Point+pointP :: Parser DPoint pointP = do   char '('; spaces   n <- num; spaces@@ -28,7 +23,7 @@   char ')'   return (Point n m) -nodeP :: Parser MetaNodeType+nodeP :: Parser (MetaNodeType Double) nodeP = option Open specialNode   where specialNode = do           char '{'; spaces@@ -65,54 +60,78 @@      string ".."      return (t1, t2) -mpRest :: Point -> Parser MetaPath-mpRest p = do+openMetaRest :: DPoint -> Parser (OpenMetaPath Double)+openMetaRest p = do   leftNode <- nodeP; spaces   string ".."; spaces   (tl, tr) <- tensionP; spaces   rightNode <- nodeP; spaces-  mp <- mpP-  return $ case mp of-    OpenMetaPath joins q ->-      (OpenMetaPath ((p, MetaJoin leftNode tl-                         tr rightNode):joins) q)-    CyclicMetaPath joins ->-      CyclicMetaPath ((p, MetaJoin leftNode tl-                         tr rightNode):joins)+  OpenMetaPath joins q <- openMetaP+  return $ OpenMetaPath ((p, MetaJoin leftNode tl+                          tr rightNode):joins) q -mpP :: Parser MetaPath-mpP =-  do p <- pointP-     spaces-     option (OpenMetaPath [] p) (mpRest p)-  <|> do-    string "cycle"-    return (CyclicMetaPath [])+closedMetaRest :: DPoint -> Parser (ClosedMetaPath Double)+closedMetaRest p = do+  leftNode <- nodeP; spaces+  string ".."; spaces+  (tl, tr) <- tensionP; spaces+  rightNode <- nodeP; spaces+  ClosedMetaPath joins <- closedMetaP+  return $ ClosedMetaPath ((p, MetaJoin leftNode tl+                             tr rightNode):joins) -pathRest :: Point -> Parser Path-pathRest p = do+openRest :: DPoint -> Parser (OpenPath Double)+openRest p = do   string ".."; spaces   string "controls"; spaces-  n <- pointP; spaces+  p1 <- pointP; spaces   string "and"; spaces-  m <- pointP; spaces+  p2 <- pointP; spaces   string ".."; spaces-  path <- pathP-  return $ case path of-    OpenPath joins q ->-      (OpenPath ((p, JoinCurve n m):joins) q)-    ClosedPath joins ->-      ClosedPath ((p, JoinCurve n m):joins)-  -pathP :: Parser Path-pathP =+  OpenPath joins q <- openP+  return $ OpenPath ((p, JoinCurve p1 p2):joins) q++closedRest :: DPoint -> Parser (ClosedPath Double)+closedRest p = do+  string ".."; spaces+  string "controls"; spaces+  p1 <- pointP; spaces+  string "and"; spaces+  p2 <- pointP; spaces+  string ".."; spaces+  ClosedPath joins <- closedP+  return $ ClosedPath ((p, JoinCurve p1 p2):joins)++openP :: Parser (OpenPath Double)+openP =    do p <- pointP      spaces-     option (OpenPath [] p) (pathRest p)+     option (OpenPath [] p) (openRest p)++closedP :: Parser (ClosedPath Double)+closedP =+  do p <- pointP+     spaces+     closedRest p   <|> do     string "cycle"     return (ClosedPath []) +openMetaP :: Parser (OpenMetaPath Double)+openMetaP =+  do p <- pointP+     spaces+     option (OpenMetaPath [] p) (openMetaRest p)++closedMetaP :: Parser (ClosedMetaPath Double)+closedMetaP =+  do p <- pointP+     spaces+     closedMetaRest p+  <|> do+    string "cycle"+    return (ClosedMetaPath [])+ tryParse :: Parser a -> String -> a tryParse p s =   case parse p "" s of@@ -125,105 +144,115 @@ doubleEq a b =   abs (a - b) < 0.01 -pointEq :: Point -> Point -> Bool+pointEq :: DPoint -> DPoint -> Bool pointEq (Point a b) (Point c d) =   doubleEq a c && doubleEq b d -joinEq :: PathJoin -> PathJoin -> Bool-joinEq JoinLine JoinLine = True+joinEq :: PathJoin Double -> PathJoin Double -> Bool+joinEq JoinLine JoinLine  = True joinEq (JoinCurve a b) (JoinCurve c d) =   pointEq a c && pointEq b d joinEq _ _ = True -pathEq :: Path -> Path -> Bool-pathEq (OpenPath joins p) (OpenPath joins2 q) =+openPathEq (OpenPath joins p) (OpenPath joins2 q) =   pointEq p q && length joins == length joins2 &&   and (zipWith    (\(p1, j1) (p2, j2) ->      pointEq p1 p2 && joinEq j1 j2)    joins joins2) -pathEq (ClosedPath joins) (ClosedPath joins2) =+closedPathEq (ClosedPath joins) (ClosedPath joins2) =   and (zipWith    (\(p1, j1) (p2, j2) ->      pointEq p1 p2 && joinEq j1 j2)    joins joins2) -thetas :: Path -> [Double]-thetas (OpenPath j p) =+openThetas :: OpenPath Double -> [Double]+openThetas (OpenPath j p) =   zipWith3 theta   (map fst j)   (tail (map fst j) ++ [p])   (map snd j)   where-    theta q r (JoinLine) = 0+    theta q r JoinLine = 0     theta q r (JoinCurve c1 _) =       vectorAngle (c1^-^q) - vectorAngle (r^-^q)     -thetas (ClosedPath j) =-  thetas (OpenPath j (fst $ head j))+closedThetas (ClosedPath j) =+  openThetas (OpenPath j (fst $ head j)) -phis :: Path -> [Double]-phis (OpenPath j p) =+openPhis :: OpenPath Double -> [Double]+openPhis (OpenPath j p) =   zipWith3 phi   (map fst j)   (tail (map fst j) ++ [p])   (map snd j)   where-    phi q r (JoinLine) = 0+    phi q r JoinLine = 0     phi q r (JoinCurve _ c2) =       vectorAngle (q^-^r) - vectorAngle (c2^-^r) -phis (ClosedPath j) =-  phis (OpenPath j (fst $ head j))+closedPhis (ClosedPath j) =+  openPhis (OpenPath j (fst $ head j)) -testPath :: TestName -> String -> TestTree-testPath p1 p2 =+testOpen :: TestName -> String -> TestTree+testOpen p1 p2 =   testCase p1 $    assertBool "Incorrect metapath." $-  unmeta (tryParse mpP p1) `pathEq`-  tryParse pathP p2+  unmetaOpen (tryParse openMetaP p1) `openPathEq`+  tryParse openP p2 +testClosed :: TestName -> String -> TestTree+testClosed p1 p2 =+  testCase p1 $ +  assertBool "Incorrect metapath." $+  unmetaClosed (tryParse closedMetaP p1) `closedPathEq`+  tryParse closedP p2+ -- These tests were created by running mf, typing expr after the -- prompt, and entering the metapaths. mfTests :: TestTree mfTests = testGroup "Metafont" [-  testPath "(0,0)..(4,3)"+  testOpen "(0,0)..(4,3)"   "(0,0)..controls (1.33333,1) and (2.66667,2) ..(4,3)", -  testPath "(0,0){(1,-2)}..(4,3)"+  testOpen "(0,0){(1,-2)}..(4,3)"   "(0,0)..controls (1.81548,-3.63095) and (6.97739,0.24046)..(4,3)", -  testPath "(0,0)..{(1,-2)}(4,3)"+  testOpen "(0,0)..{(1,-2)}(4,3)"   "(0,0)..controls (-2.97739,2.75954) and (2.18452,6.63095)..(4,3)", -  testPath "(0,0){curl 2}..(4,3)"+  testOpen "(0,0){curl 2}..(4,3)"   "(0,0)..controls (1.33333,1) and (2.66667,2)..(4,3)", -  testPath "(0,0){(2, 3)}..{(1, 2)}(4,3)"+  testOpen "(0,0){(2, 3)}..{(1, 2)}(4,3)"   "(0,0)..controls (0.95523,1.43285) and (3.21622,1.43243)..(4,3)", -  testPath "(0,0)..(4,3)..(-2, 1)"+  testOpen "(0,0)..(4,3)..(-2, 1)"   "(0,0)..controls (2.08194,-1.42896) and (4.78885,0.60123)..(4,3)..controls (2.67747,7.02158) and (-3.35492,5.01077)..(-2,1)", -  testPath "(1,1)..tension 0.8 and 1.2..(3,4)..tension 10 ..(-10,-10)"+  testOpen "(1,1)..tension 0.8 and 1.2..(3,4)..tension 10 ..(-10,-10)"   "(1,1)..controls (-2.7088,-12.93713) and (13.27118,14.12433)..(3,4)..controls (2.54623,3.55272) and (-9.58751,-9.5144)..(-10,-10)",   -  testPath "(0,0){curl 2}..(4,3)..(-2, 1)"+  testOpen "(0,0){curl 2}..(4,3)..(-2, 1)"   "(0,0)..controls (1.14464,-2.66646) and (6.04007,-0.56508)..(4,3)..controls (2.2501,6.05801) and (-2.43489,4.49635)..(-2,1)", -  testPath "(0,0){(-3, -2)}..(4,3)..(-2, 1)"+  testOpen "(0,0){(-3, -2)}..(4,3)..(-2, 1)"   "(0,0)..controls (-3.65675,-2.43784) and (1.35551,2.07506)..(4,3)..controls (27.8797,11.35223) and (-26.11505,-6.64606)..(-2,1)", -  testPath "(0,0)..(2,3)..(4,4)..cycle"+  testClosed "(0,0)..(2,3)..(4,4)..cycle"   "(0,0)..controls (-0.27211,1.267) and (0.9676,2.15346)..(2,3)..controls (2.60509,3.49615) and (3.2241,4.08679)..(4,4)..controls (12.90535,3.00386) and (1.91997,-8.93997)..cycle", -  testPath "(0,0)..tension 0.9 and 1.1 ..(2,3)..(4,4)..cycle"+  testClosed "(0,0)..tension 0.9 and 1.1 ..(2,3)..(4,4)..cycle"   "(0,0)..controls (-0.39941,1.39384) and (0.99234,2.26094)..(2,3)..controls (2.62666,3.45963) and (3.22433,4.07909)..(4,4)..controls (12.2955,3.15413) and (2.40324,-8.38663)..cycle", -    testPath "(0,0)..(2,3){(1,1)}..(4,4)..cycle"+  testClosed "(0,0)..(2,3){(1,1)}..(4,4)..cycle"   "(0,0)..controls (-0.24208,1.27483) and (1.07744,2.07744)..(2,3)..controls (2.56248,3.56248) and (3.22197,4.11229)..(4,4)..controls (12.86206,2.72092) and (1.68616,-8.87949)..cycle"   ]++tests :: TestTree+tests = testGroup "Tests" [mfTests]+  main :: IO () main = defaultMain tests