packages feed

cubicbezier 0.2.0 → 0.3.0

raw patch · 8 files changed

+451/−143 lines, 8 filesdep +cubicbezierdep +parsecdep +tastydep −deepseqdep ~basedep ~containersdep ~integration

Dependencies added: cubicbezier, parsec, tasty, tasty-hunit

Dependencies removed: deepseq

Dependency ranges changed: base, containers, integration

Files

Geom2D.hs view
@@ -10,6 +10,7 @@ data Point = Point {   pointX :: {-# UNPACK #-} !Double,   pointY :: {-# UNPACK #-} !Double}+           deriving Eq  instance Show Point where   show (Point x y) =@@ -65,7 +66,9 @@         b' = x2 - x1         d = sqrt(a'*a' + b'*b') --- | Return the the distance from a point to the line.+-- | Return the signed distance from a point to the line.  If the+-- distance is negative, the point lies to the right of the line+         lineDistance :: Line -> Point -> Double lineDistance l = \(Point x y) -> a*x + b*y + c   where (a, b, c) = lineEquation l
Geom2D/CubicBezier/Basic.hs view
@@ -212,9 +212,17 @@   bezierSubsegment c t u :   tail (splitBezierN c $ u:rest) --- | Return True if all the control points are colinear within tolerance.+-- | Return False if some points fall outside a line with a thickness of the given tolerance.++-- fat line calculation taken from the bezier-clipping algorithm (Sederberg) colinear :: CubicBezier -> Double -> Bool-colinear (CubicBezier !a !b !c !d) eps =-  abs (ld b) < eps && abs (ld c) < eps+colinear (CubicBezier !a !b !c !d) eps = dmax - dmin < eps   where ld = lineDistance (Line a d)+        d1 = ld b+        d2 = ld c+        (dmin, dmax) | d1*d2 > 0 = (3/4 * minimum [0, d1, d2],+                                    3/4 * maximum [0, d1, d2])+                     | otherwise = (4/9 * minimum [0, d1, d2],+                                    4/9 * maximum [0, d1, d2])+
Geom2D/CubicBezier/Intersection.hs view
@@ -1,4 +1,4 @@-{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE BangPatterns, MultiWayIf #-} -- | Intersection routines using Bezier Clipping.  Provides also functions for finding the roots of onedimensional bezier curves.  This can be used as a general polynomial root solver by converting from the power basis to the bernstein basis. module Geom2D.CubicBezier.Intersection        (bezierIntersection, bezierLineIntersections, bezierFindRoot, closest)@@ -94,6 +94,14 @@   -- no intersection   | isNothing chop_interval = [] +  -- within tolerance      +  | max (umax - umin) (new_tmax - new_tmin) < eps =+    if revCurves+    then [ (umin + (umax-umin)/2,+            new_tmin + (new_tmax-new_tmin)/2) ]+    else [ (new_tmin + (new_tmax-new_tmin)/2,+            umin + (umax-umin)/2) ]+   -- not enough reduction, so split the curve in case we have   -- multiple intersections   | prevClip > 0.8 && newClip > 0.8 =@@ -109,17 +117,8 @@       in bezierClip ql newP umin half_t new_tmin new_tmax newClip eps (not revCurves) ++          bezierClip qr newP half_t umax new_tmin new_tmax newClip eps (not revCurves) -  -- within tolerance      -  | max (umax - umin) (new_tmax - new_tmin) < eps =-    if revCurves-    then [ (umin + (umax-umin)/2,-            new_tmin + (new_tmax-new_tmin)/2) ]-    else [ (new_tmin + (new_tmax-new_tmin)/2,-            umin + (umax-umin)/2) ]-   -- iterate with the curves reversed.-  | otherwise =-      bezierClip q newP umin umax new_tmin new_tmax newClip eps (not revCurves)+  | otherwise = bezierClip q newP umin umax new_tmin new_tmax newClip eps (not revCurves)    where     d = lineDistance (Line q0 q3)@@ -137,13 +136,10 @@     new_tmin = tmax * chop_tmin + tmin * (1 - chop_tmin)     new_tmax = tmax * chop_tmax + tmin * (1 - chop_tmax) --- | Find the intersections between two Bezier curves within given--- tolerance, using the Bezier Clip algorithm. Returns the parameters--- for both curves.+-- | Find the intersections between two Bezier curves, using the+-- Bezier Clip algorithm. Returns the parameters for both curves. bezierIntersection :: CubicBezier -> CubicBezier -> Double -> [(Double, Double)]-bezierIntersection p q eps = bezierClip p q 0 1 0 1 0 eps' False-  where-    eps' = min (bezierParamTolerance p eps) (bezierParamTolerance q eps)+bezierIntersection p q eps = bezierClip p q 0 1 0 1 0 eps False  ------------------------ Line intersection ------------------------------------- -- Clipping a line uses a simplified version of the Bezier Clip algorithm,@@ -198,11 +194,18 @@  -- | Find the closest value(s) on the bezier to the given point, within tolerance. closest :: CubicBezier -> Point -> Double -> [Double]-closest cb (Point px py) eps = bezierFindRoot poly 0 1 eps+closest cb p@(Point px py) eps =+  map fst $ filter (\(_, d) -> abs (d - closestDist) < eps/2) $+  zip tVals dists   where+    closestDist = minimum dists+    dists = map (vectorDistance p . evalBezier cb) tVals+    tVals = 0:1:bezierFindRoot poly 0 1 eps     (bx, by) = bezierToBernstein cb     bx' = bernsteinDeriv bx     by' = bernsteinDeriv by     poly = (bx ~- listToBernstein [px, px, px, px]) ~* bx' ~+            (by ~- listToBernstein [py, py, py, py]) ~* by' +-- let cb = (CubicBezier (Point 0 0) (Point 3 4) (Point 10 4) (Point 31 2)); cb1 = fst (splitBezier cb 0.83242); cb2 = CubicBezier {bezierC0 = Point 4.542593123258268 2.7028033902052537, bezierC1 = Point 9.036628467934 3.788306467438, bezierC2 = Point 16.832161 3.4493180000000002, bezierC3 = Point 31.0 2.0}+-- bezierIntersection (CubicBezier (Point 0 0) (Point 3 4) (Point 10 4) (Point 31 2)) (CubicBezier (Point 0 0) (Point 6 8) (Point 2 42) (Point 4 15)) 1e-10
Geom2D/CubicBezier/MetaPath.hs view
@@ -1,8 +1,8 @@ {-# LANGUAGE BangPatterns #-} -- | This module implements an extension to paths as used in--- D.E.Knuth's /Metafont/.  Metafont gives a more intuitive method to--- specify paths than bezier curves.  I'll give a brief overview of--- the metafont curves.  For a more in depth explanation look at+-- D.E.Knuth's /Metafont/.  Metafont gives an alternate way+-- to specify paths using bezier curves.  I'll give a brief overview of+-- the metafont curves.  A more in depth explanation can be found in  -- /The MetafontBook/. --  -- Each spline has a tension parameter, which is a relative measure of@@ -75,7 +75,7 @@ data MetaNodeType = Open                   | Curl {curlgamma :: Double}                   | Direction {nodedir :: Point}-                  deriving Show+                  deriving (Eq, Show)  data Tension = Tension {tensionValue :: Double}              | TensionAtLeast {tensionValue :: Double}@@ -87,7 +87,7 @@   show (OpenMetaPath nodes lastpoint) =     showPath nodes ++ showPoint lastpoint -showPath :: [(Point, MetaJoin)] -> [Char]+showPath :: [(Point, MetaJoin)] -> String showPath = concatMap showNodes   where     showNodes (p, Controls u v) =@@ -112,16 +112,29 @@ -- | Create a normal path from a metapath. unmeta :: MetaPath -> Path unmeta (OpenMetaPath nodes endpoint) =-  unmetaOpen (sanitizeOpen nodes) endpoint+  unmetaOpen (flip sanitize endpoint $ removeEmptyDirs nodes) endpoint  unmeta (CyclicMetaPath nodes) =-  case span (bothOpen . snd) nodes of-    (l, []) -> unmetaCyclic l-    (l, (m:n)) ->-      if leftOpen $ snd m+  case spanList bothOpen (removeEmptyDirs nodes) of+    ([], []) -> error "empty metapath"+    (l, []) -> if fst (last l) == fst (head l)+               then unmetaAsOpen l []+               else unmetaCyclic l+    (l, m:n) ->+      if leftOpen (m:n)       then unmetaAsOpen (l++[m]) n       else unmetaAsOpen l (m:n) +-- solve a cyclic metapath as an open path if possible.+-- rotate to the defined node, and rotate back after+-- solving the path.+unmetaAsOpen :: [(Point, MetaJoin)] -> [(Point, MetaJoin)] -> Path+unmetaAsOpen l m = ClosedPath (l'++m') +  where n = length m+        OpenPath o _ =+          unmetaOpen (sanitizeCycle (m++l)) (fst $ head (m ++l))+        (m',l') = splitAt n o+ unmetaOpen :: [(Point, MetaJoin)] -> Point -> Path unmetaOpen nodes endpoint =   let subsegs = openSubSegments nodes endpoint@@ -130,19 +143,16 @@  -- decompose into a list of subsegments that need to be solved. openSubSegments :: [(Point, MetaJoin)] -> Point -> [MetaPath]-openSubSegments l p = openSubSegments' (tails l) p--openSubSegments' :: [[(Point, MetaJoin)]] -> Point -> [MetaPath]-openSubSegments' [[]] _ = []-openSubSegments' [] _   = []-openSubSegments' l lastPoint = case break breakPoint l of-  (m, n:o) ->-    let point = case o of-          (((p,_):_):_) -> p-          _ -> lastPoint-    in OpenMetaPath (map head (m ++ [n])) point :-       openSubSegments' o lastPoint-  _ -> error "openSubSegments': unexpected end of segments"+openSubSegments [] _   = []+openSubSegments l lastPoint =+  case spanList (not . breakPoint) l of+    (m, n:o) ->+      let point = case o of+            ((p,_):_) -> p+            _ -> lastPoint+      in OpenMetaPath (m ++ [n]) point :+         openSubSegments o lastPoint+    _ -> error "openSubSegments': unexpected end of segments"  -- join subsegments into one segment joinSegments :: [Path] -> [(Point, PathJoin)]@@ -154,29 +164,20 @@ unmetaCyclic :: [(Point, MetaJoin)] -> Path unmetaCyclic nodes =   let points = map fst nodes-      chords = zipWith (^-^) points (last points : points)-      tensionsA = (map (tensionL . snd) nodes)-      tensionsB = (map (tensionR . snd) nodes)+      chords = zipWith (^-^) (tail $ cycle points) points+      tensionsA = map (tensionL . snd) nodes+      tensionsB = map (tensionR . snd) nodes       turnAngles = zipWith turnAngle chords (tail $ cycle chords)       thetas = solveCyclicTriD $                eqsCycle tensionsA                points                tensionsB                turnAngles-      phis = zipWith (\x y -> -(x+y)) turnAngles (tail thetas ++ [head thetas])+      phis = zipWith (\x y -> -(x+y)) turnAngles (tail $ cycle thetas)   in ClosedPath $ zip points $-     zipWith6 unmetaJoin points (tail points ++ [head points])+     zipWith6 unmetaJoin points (tail $ cycle points)      thetas phis tensionsA tensionsB --- solve a cyclic metapath as an open path if possible.--- rotate to the defined node, and rotate back after--- solving the path.-unmetaAsOpen :: [(Point, MetaJoin)] -> [(Point, MetaJoin)] -> Path-unmetaAsOpen l m = ClosedPath (l'++m') -  where n = length m-        OpenPath o _ = unmetaOpen (sanitizeCycle (m++l)) (fst $ head m)-        (m',l') = splitAt n o- -- solve a subsegment unmetaSubSegment :: MetaPath -> Path @@ -189,82 +190,106 @@   let points = map fst nodes ++ [lastpoint]       joins = map snd nodes       chords = zipWith (^-^) (tail points) points-      tensionsA = map tensionL  joins+      tensionsA = map tensionL joins       tensionsB = map tensionR joins       turnAngles = zipWith turnAngle chords (tail chords) ++ [0]       thetas = solveTriDiagonal $                eqsOpen points joins chords turnAngles                (map tensionValue tensionsA)                (map tensionValue tensionsB)-      phis = zipWith (\x y -> -x-y) turnAngles (tail thetas)-      pathjoins = zipWith6 unmetaJoin points (tail points) thetas phis tensionsA tensionsB+      phis = zipWith (\x y -> -(x+y)) turnAngles (tail thetas)+      pathjoins =+        zipWith6 unmetaJoin points (tail points) thetas phis tensionsA tensionsB   in OpenPath (zip points pathjoins) lastpoint  unmetaSubSegment _ = error "unmetaSubSegment: subsegment should not be cyclic" -bothOpen :: MetaJoin -> Bool-bothOpen (MetaJoin Open _ _ Open) = True+removeEmptyDirs :: [(Point, MetaJoin)] -> [(Point, MetaJoin)]+removeEmptyDirs = map remove+  where remove (p, MetaJoin (Direction (Point 0 0)) tl tr jr) = remove (p, MetaJoin Open tl tr jr)+        remove (p, MetaJoin jl tl tr (Direction (Point 0 0))) = (p, MetaJoin jl tl tr Open)+        remove j = j++-- if p == q, it will become a control point+bothOpen :: [(Point, MetaJoin)] -> Bool+bothOpen ((p, MetaJoin Open _ _ Open):(q, _):_) = p /= q  +bothOpen [(_, MetaJoin Open _ _ Open)] = True bothOpen _ = False -leftOpen :: MetaJoin -> Bool-leftOpen (MetaJoin Open _ _ _) = True+leftOpen :: [(Point, MetaJoin)] -> Bool+leftOpen ((p, MetaJoin Open _ _ _):(q, _):_) = p /= q  +leftOpen [(_, MetaJoin Open _ _ _)] = True leftOpen _ = False -replaceLast :: [a] -> a -> [a]-replaceLast [] _ = []-replaceLast [_] n = [n]-replaceLast (l:ls) n = l:replaceLast ls n- sanitizeCycle :: [(Point, MetaJoin)] -> [(Point, MetaJoin)]-sanitizeCycle l = replaceLast ls l'-  where-    (l':ls) = sanitizeRest (last l: l)---- replace open nodetypes with more defined nodetypes if possible-sanitizeOpen :: [(Point, MetaJoin)] -> [(Point, MetaJoin)]-sanitizeOpen [] = []+sanitizeCycle [] = []+sanitizeCycle l = take n $ tail $+                  sanitize (drop (n-1) $ cycle l) (fst $ head l)+  where n = length l --- starting open => curl-sanitizeOpen ((p, MetaJoin Open t1 t2 m):rest) =-  sanitizeRest ((p, MetaJoin (Curl 1) t1 t2 m):rest)-sanitizeOpen l = sanitizeRest l-   -sanitizeRest :: [(Point, MetaJoin)] -> [(Point, MetaJoin)]-sanitizeRest [] = []+sanitize :: [(Point, MetaJoin)] -> Point -> [(Point, MetaJoin)]+sanitize [] _ = []  -- ending open => curl-sanitizeRest [(p, MetaJoin m t1 t2 Open)] =-  [(p, MetaJoin m t1 t2 (Curl 1))]+sanitize [(p, MetaJoin m t1 t2 Open)] r =+  if p == r+  then [(p, Controls p p)]+  else [(p, MetaJoin m t1 t2 (Curl 1))] -sanitizeRest (node1@(p, MetaJoin m1 tl tr m2): node2@(q, MetaJoin n1 sl sr n2): rest) =-  case (m2, n1) of-    (Curl g, Open) -> -- curl, open => curl, curl-      node1 : sanitizeRest ((q, MetaJoin (Curl g) sl sr n2):rest)-    (Open, Curl g) -> -- open, curl => curl, curl-      (p, MetaJoin m1 tl tr (Curl g)) : sanitizeRest (node2:rest)-    (Direction dir, Open) ->   -- given, open => given, given-      node1 : sanitizeRest ((q, (MetaJoin (Direction dir) sl sr n2)) : rest)-    (Open, Direction dir) ->   -- open, given => given, given-      (p, MetaJoin m1 tl tr (Direction dir)) : sanitizeRest (node2:rest)-    _ -> node1 : sanitizeRest (node2:rest)+sanitize ((p, MetaJoin m1 tl tr Open): rest@(node2:node3:_)) r+  | (fst node2 == fst node3) && (metaTypeL (snd node2) == Open) =+    (p, MetaJoin m1 tl tr (Curl 1)) : sanitize rest r+    +sanitize (node1@(p, MetaJoin m1 tl tr m2): node2@(q, MetaJoin n1 sl sr n2): rest) r+  | p == q =+    -- if two consecutive points are the same, just make a curve with all control points the same+    -- we still have to propagate a curl or given direction.+    let newnode = (p, Controls p p)+    in case (m2, n1) of+      (Curl g, Open) -> -- curl, open => explicit, curl+        newnode : sanitize ((q, MetaJoin (Curl g) sl sr n2):rest) r+      (Direction dir, Open) ->   -- given, open => explicit, given+        newnode : sanitize ((q, MetaJoin (Direction dir) sl sr n2) : rest) r+      (Open, Open) ->   -- open, open => explicit, curl+        newnode : sanitize ((q, MetaJoin (Curl 1) sl sr n2) : rest) r+      _ -> newnode : sanitize (node2:rest) r+  | otherwise =+    case (m2, n1) of+      (Curl g, Open) -> -- curl, open => curl, curl+        node1 : sanitize ((q, MetaJoin (Curl g) sl sr n2):rest) r+      (Open, Curl g) -> -- open, curl => curl, curl+        (p, MetaJoin m1 tl tr (Curl g)) : sanitize (node2:rest) r+      (Direction dir, Open) ->   -- given, open => given, given+        node1 : sanitize ((q, MetaJoin (Direction dir) sl sr n2) : rest) r+      (Open, Direction dir) ->   -- open, given => given, given+        (p, MetaJoin m1 tl tr (Direction dir)) : sanitize (node2:rest) r+      _ -> node1 : sanitize (node2:rest) r -sanitizeRest ((p, m): (q, n): rest) =+sanitize ((p, m): (q, n): rest) r =   case (m, n) of-    (Controls _u v, MetaJoin Open t1 t2 mt2) ->  -- explicit, open => explicit, given-      (p, m) : sanitizeRest ((q, MetaJoin (Direction (q^-^v)) t1 t2 mt2): rest)-    (MetaJoin mt1 tl tr Open, Controls u _v) ->  -- open, explicit => given, explicit-      (p, MetaJoin mt1 tl tr (Direction (u^-^p))) : sanitizeRest ((q, n): rest)-    _ -> (p, m) : sanitizeRest ((q, n) : rest)+    (Controls _u v, MetaJoin Open t1 t2 mt2) -- explicit, open => explicit, given+      | q == v    -> (p, m) : sanitize ((q, MetaJoin (Curl 1) t1 t2 mt2): rest) r+      | otherwise -> (p, m) : sanitize ((q, MetaJoin (Direction (q^-^v)) t1 t2 mt2): rest) r+    (MetaJoin mt1 tl tr Open, Controls u _v) -- open, explicit => given, explicit+      | u == p    -> (p, MetaJoin mt1 tl tr (Curl 1)) : sanitize ((q, n): rest) r +      | otherwise -> (p, MetaJoin mt1 tl tr (Direction (u^-^p))) : sanitize ((q, n): rest) r+    _ -> (p, m) : sanitize ((q, n) : rest) r -sanitizeRest (n:l) = n:sanitizeRest l+sanitize (n:l) r = n:sanitize l r +spanList :: ([a] -> Bool) -> [a] -> ([a], [a])+spanList _ xs@[] =  (xs, xs)+spanList p xs@(x:xs')+  | p xs =  let (ys,zs) = spanList p xs' in (x:ys,zs)+  | otherwise    =  ([],xs)+ -- break the subsegment if the angle to the left or the right is defined or a curl. breakPoint :: [(Point, MetaJoin)] -> Bool breakPoint ((_,  MetaJoin _ _ _ Open):(_, MetaJoin Open _ _ _):_) = False breakPoint _ = True  -- solve the tridiagonal system for t[i]:--- a[n] t[i-1] + b[i] t[i] + c[b] t[i+1] = d[i]+-- a[n] t[i-1] + b[n] t[i] + c[n] t[i+1] = d[i] -- where a[0] = c[n] = 0 -- by first rewriting it into -- the system t[i] + u[i] t[i+1] = v[i]@@ -282,8 +307,19 @@     solutions = reverse $ scanl nextsol vn twovars     nextsol ti (u, v) = v - u*ti --- 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)])+-- solveTriDiagonal2 :: (Double, Double, Double) -> V.Vector (Double, Double, Double, Double) -> V.Vector Double+-- solveTriDiagonal2 (!b0, !c0, !d0) rows = solutions+--   where+--     solutions = undefined+--     twovars = V.scanl nextrow (c0/b0, d0/b0) rows+--     solutions = scanr V.unsafeInit+--     nextrow (u, v) (ai, bi, ci, di) =+--       (ci/(bi - u*ai), (di - v*ai)/(bi - u*ai)) +-- 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]@@ -302,22 +338,23 @@     nextsol t (!u, !v, !w) = (v + w*t0 - t)/u  turnAngle :: Point -> Point -> Double+turnAngle (Point 0 0) _ = 0 turnAngle (Point x y) q = vectorAngle $ rotateVec p $* q   where p = Point x (-y) -zipPrev :: [a] -> [(a, a)]-zipPrev [] = []-zipPrev l = zip (last l : l) l+zipNext :: [b] -> [(b, b)]+zipNext [] = []+zipNext l = zip l (tail $ cycle l)  -- find the equations for a cycle containing only open points eqsCycle :: [Tension] -> [Point] -> [Tension]          -> [Double] -> [(Double, Double, Double, Double)] eqsCycle tensionsA points tensionsB turnAngles =    zipWith4 eqTension-  (zipPrev (map tensionValue tensionsA))-  (zipPrev dists)-  (zipPrev turnAngles)-  (zipPrev (map tensionValue tensionsB))+  (zipNext (map tensionValue tensionsA))+  (zipNext dists)+  (zipNext turnAngles)+  (zipNext (map tensionValue tensionsB))   where      dists = zipWith vectorDistance points (tail $ cycle points) @@ -326,20 +363,22 @@  eqsOpen :: [Point] -> [MetaJoin] -> [Point] -> [Double]         -> [Double] -> [Double] -> [(Double, Double, Double, Double)]-eqsOpen _ [join] [delta] _ _ _ =-  case join of-    MetaJoin (Curl _) _ _ (Curl _) ->-      [(0, 1, 0, 0), (0, 1, 0, 0)]-    MetaJoin (Curl g) t1 t2 (Direction dir) ->+eqsOpen _ [MetaJoin mt1 t1 t2 mt2] [delta] _ _ _ =+  let replaceType Open = Curl 1+      replaceType t = t+  in case (replaceType mt1, replaceType mt2) of+    (Curl g, Direction dir) ->       [eqCurl0 g (tensionValue t1) (tensionValue t2) 0,        (0, 1, 0, turnAngle delta dir)]-    MetaJoin (Direction dir) t1 t2 (Curl g) ->+    (Direction dir, Curl g) ->       [(0, 1, 0, turnAngle delta dir),        eqCurlN g (tensionValue t1) (tensionValue t2)]-    MetaJoin (Direction dir) _ _ (Direction dir2) ->+    (Direction dir, Direction dir2) ->       [(0, 1, 0, turnAngle delta dir),        (0, 1, 0, turnAngle delta dir2)]-    _ -> error "eqsOpen: illegal nodetype in subsegment"+    (Curl _, Curl _) ->+      [(0, 1, 0, 0), (0, 1, 0, 0)]+    _ -> undefined  eqsOpen points joins chords turnAngles tensionsA tensionsB =   eq0 : restEquations joins tensionsA dists turnAngles tensionsB@@ -348,13 +387,15 @@     eq0 = case head joins of       (MetaJoin (Curl g) _ _ _) -> eqCurl0 g (head tensionsA) (head tensionsB) (head turnAngles)       (MetaJoin (Direction dir) _ _ _) -> (0, 1, 0, turnAngle (head chords) dir)-      _ -> error "eqsOpen: illegal subsegment first nodetype"+      (MetaJoin Open _ _ _) -> eqCurl0 1 (head tensionsA) (head tensionsB) (head turnAngles)+      (Controls _ _) -> error "eqsOpen: illegal join"      restEquations [lastnode] (tensionA:_) _ _ (tensionB:_) =       case lastnode of         MetaJoin _ _ _ (Curl g) -> [eqCurlN g tensionA tensionB]+        MetaJoin _ _ _ Open -> [eqCurlN 1 tensionA tensionB]         MetaJoin _ _ _ (Direction dir) -> [(0, 1, 0, turnAngle (last chords) dir)]-        _  -> error "eqsOpen: illegal subsegment last nodetype"+        (Controls _ _) -> error "eqsOpen: illegal join"      restEquations (_:othernodes) (tensionA:restTA) (d:restD) (turn:restTurn) (tensionB:restTB) =       eqTension (tensionA, head restTA) (d, head restD) (turn, head restTurn) (tensionB, head restTB) :@@ -369,7 +410,7 @@ eqTension (tensionA', tensionA) (dist', dist) (psi', psi) (tensionB', tensionB) =   (a, b+c, d, -b*psi' - d*psi)   where-    a = (tensionB' * tensionB' / (tensionA' * dist'))+    a = tensionB' * tensionB' / (tensionA' * dist')     b = (3 - 1/tensionA') * tensionB' * tensionB' / dist'     c = (3 - 1/tensionB) * tensionA * tensionA / dist     d = tensionA * tensionA / (tensionB * dist)@@ -391,7 +432,7 @@     b = chi/tensionB + 3 - 1/tensionA     chi = gamma*tensionA*tensionA / (tensionB*tensionB) --- magic formula for getting the control points by John Hobby+-- getting the control points unmetaJoin :: Point -> Point -> Double -> Double -> Tension -> Tension -> PathJoin unmetaJoin !z0 !z1 !theta !phi !alpha !beta   | abs phi < 1e-4 && abs theta < 1e-4 = JoinLine@@ -414,17 +455,20 @@           TensionAtLeast _ | bounded ->             min ss' (st/stf)           _ -> ss'-        u = z0 ^+^ rr *^ Point (dx*ct - dy*st) (dy*ct + dx*st)  -- z0 + rr * (rotate theta chord)-        v = z1 ^-^ ss *^ Point (dx*cf + dy*sf) (dy*cf - dx*sf)  -- z1 - ss * (rotate (-phi) chord)+        -- u = z0 + rr * (rotate theta chord)+        u = z0 ^+^ rr *^ Point (dx*ct - dy*st) (dy*ct + dx*st)+        -- v = z1 - ss * (rotate (-phi) chord)+        v = z1 ^-^ ss *^ Point (dx*cf + dy*sf) (dy*cf - dx*sf)  constant1, constant2, sqrt2 :: Double-constant1 = 3/2*(sqrt 5 - 1)-constant2 = 3/2*(3 - sqrt 5)+constant1 = 3*(sqrt 5 - 1)/2+constant2 = 3*(3 - sqrt 5)/2 sqrt2 = sqrt 2  -- another magic formula by John Hobby. velocity :: Double -> Double -> Double          -> Double -> Tension -> Double velocity st sf ct cf t =+  min 4 $    (2 + sqrt2 * (st - sf/16)*(sf - st/16)*(ct - cf)) /   ((3 + constant1*ct + constant2*cf) * tensionValue t)
Geom2D/CubicBezier/Numeric.hs view
@@ -3,15 +3,20 @@  -- | @quadraticRoot a b c@ find the real roots of the quadratic equation -- @a x^2 + b x + c = 0@.  It will return one, two or zero roots.+ quadraticRoot :: Double -> Double -> Double -> [Double]-quadraticRoot a b c = result where-  d = b*b - 4*a*c-  q = - (b + signum b * sqrt d) / 2-  x1 = q/a-  x2 = c/q-  result | d < 0     = []-         | d == 0    = [x1]-         | otherwise = [x1, x2]+quadraticRoot a b c+  | a == 0 && b == 0 = []+  | a == 0 = [-c/b]+  | otherwise = result+  where+    d = b*b - 4*a*c+    q = - (b + signum b * sqrt d) / 2+    x1 = q/a+    x2 = c/q+    result | d < 0     = []+           | d == 0    = [x1]+           | otherwise = [x1, x2]  -- | @solveLinear2x2 a b c d e f@ solves the linear equation with two variables (x and y) and two systems: --
Math/BernsteinPoly.hs view
@@ -72,7 +72,8 @@  -- | Evaluate the bernstein polynomial. bernsteinEval :: BernsteinPoly -> Double -> Double-bernsteinEval (BernsteinPoly lp [b]) _ = b+bernsteinEval (BernsteinPoly _ []) _ = error "illegal bernstein polynomial"+bernsteinEval (BernsteinPoly _ [b]) _ = b bernsteinEval (BernsteinPoly lp (b':bs)) t = go t n (b'*u) 2 bs   where u = 1-t         n = fromIntegral lp@@ -83,6 +84,8 @@           ((tmp + tn*bc*b)*u) -- tmp           (i+1)             -- i           rest+        go _ _ _ _ [] = error "impossible"+          -- | Evaluate the bernstein polynomial and its derivatives. bernsteinEvalDerivs :: BernsteinPoly -> Double -> [Double]
cubicbezier.cabal view
@@ -1,8 +1,8 @@ Name:		cubicbezier-Version: 	0.2.0+Version: 	0.3.0 Synopsis:	Efficient manipulating of 2D cubic bezier curves. Category: 	Graphics, Geometry, Typography-Copyright: 	Kristof Bastiaensen (2013)+Copyright: 	Kristof Bastiaensen (2014) Stability:	Unstable License:	BSD3 License-file:	LICENSE@@ -10,7 +10,7 @@ Maintainer:	Kristof Bastiaensen Bug-Reports: 	https://github.com/kuribas/cubicbezier/issues Build-type:	Simple-Cabal-version:	>=1.6+Cabal-version:	>=1.8 Description:	This library supports efficient manipulating of 2D cubic bezier curves.  The original goal   is to support typography, but it is useful for general graphics.  Supported features are:   .@@ -27,7 +27,7 @@  Library   Ghc-options: -Wall-  Build-depends: base >= 3 && < 5, containers > 0.4, integration >= 0.1.1, deepseq >= 1.3.0+  Build-depends: base >= 3 && < 5, containers > 0.4, integration >= 0.1.1   Exposed-Modules:     Geom2D     Geom2D.CubicBezier@@ -40,3 +40,16 @@     Math.BernsteinPoly   Other-Modules:     Geom2D.CubicBezier.Numeric++test-suite test+  type: exitcode-stdio-1.0+  hs-source-dirs:+    tests+  main-is:+    test.hs+  build-depends:+    base >= 4 && < 5,+    tasty >= 0.8,+    tasty-hunit >= 0.9,+    parsec >= 3.0,+    cubicbezier
+ tests/test.hs view
@@ -0,0 +1,229 @@+{-# Language ViewPatterns #-}++import Test.Tasty+import Test.Tasty.HUnit+import Geom2D.CubicBezier+import Control.Monad+import Text.Parsec+import Text.Parsec.String+import Text.Parsec.Error++tests :: TestTree+tests = testGroup "Tests" [unitTests]++num :: Parser Double+num = +  liftM (read.concat) $ sequence+  [ option "" $ string "-"+  , many1 digit+  , option "" $ string "."+  , option "0" $ many digit]++pointP :: Parser Point+pointP = do+  char '('; spaces+  n <- num; spaces+  char ','; spaces+  m <- num ; spaces+  char ')'+  return (Point n m)++nodeP :: Parser MetaNodeType+nodeP = option Open specialNode+  where specialNode = do+          char '{'; spaces+          node <- choice [+            do string "curl"+               spaces+               n <- num+               return (Curl n),+            do p <- pointP+               return $ Direction p]+          spaces; char '}'+          return node++tensionAmount :: Parser Tension+tensionAmount = do+  cons <- option Tension+          (string "atleast" >>+           return TensionAtLeast)+  spaces+  n <- num+  return $ cons n+         +tensionP :: Parser (Tension, Tension)+tensionP =+  option (Tension 1, Tension 1) $+  do string "tension";+     spaces;+     t1 <- tensionAmount+     spaces+     t2 <- option t1 (do string "and"+                         spaces+                         tensionAmount)+     spaces+     string ".."+     return (t1, t2)++mpRest :: Point -> Parser MetaPath+mpRest 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)++mpP :: Parser MetaPath+mpP =+  do p <- pointP+     spaces+     option (OpenMetaPath [] p) (mpRest p)+  <|> do+    string "cycle"+    return (CyclicMetaPath [])++pathRest :: Point -> Parser Path+pathRest p = do+  string ".."; spaces+  string "controls"; spaces+  n <- pointP; spaces+  string "and"; spaces+  m <- 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 =+  do p <- pointP+     spaces+     option (OpenPath [] p) (pathRest p)+  <|> do+    string "cycle"+    return (ClosedPath [])++tryParse :: Parser a -> String -> a+tryParse p s =+  case parse p "" s of+   Left err -> error $ concatMap messageString $+               errorMessages err+   Right res -> res+  +  +doubleEq :: (Ord a, Fractional a) => a -> a -> Bool+doubleEq a b =+  abs (a - b) < 0.01++pointEq :: Point -> Point -> Bool+pointEq (Point a b) (Point c d) =+  doubleEq a c && doubleEq b d++joinEq :: PathJoin -> PathJoin -> 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) =+  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) =+  and (zipWith+   (\(p1, j1) (p2, j2) ->+     pointEq p1 p2 && joinEq j1 j2)+   joins joins2)++thetas :: Path -> [Double]+thetas (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 (JoinCurve c1 _) =+      vectorAngle (c1^-^q) - vectorAngle (r^-^q)+    +thetas (ClosedPath j) =+  thetas (OpenPath j (fst $ head j))++phis :: Path -> [Double]+phis (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 (JoinCurve _ c2) =+      vectorAngle (q^-^r) - vectorAngle (c2^-^r)++phis (ClosedPath j) =+  phis (OpenPath j (fst $ head j))++testPath :: TestName -> String -> TestTree+testPath p1 p2 =+  testCase p1 $ +  assertBool "Incorrect metapath." $+  unmeta (tryParse mpP p1) `pathEq`+  tryParse pathP p2++-- These tests were created by running mf, typing expr after the+-- prompt, and entering the metapaths.+unitTests :: TestTree+unitTests = testGroup "Metafont" [+  testPath "(0,0)..(4,3)"+  "(0,0)..controls (1.33333,1) and (2.66667,2) ..(4,3)",++  testPath "(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)"+  "(0,0)..controls (-2.97739,2.75954) and (2.18452,6.63095)..(4,3)",++  testPath "(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)"+  "(0,0)..controls (0.95523,1.43285) and (3.21622,1.43243)..(4,3)",++  testPath "(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)"+  "(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)"+  "(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)"+  "(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"+  "(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"+  "(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"+  "(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"+  ]++main :: IO ()+main = defaultMain tests