diff --git a/Geom2D.hs b/Geom2D.hs --- a/Geom2D.hs +++ b/Geom2D.hs @@ -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 diff --git a/Geom2D/CubicBezier/Basic.hs b/Geom2D/CubicBezier/Basic.hs --- a/Geom2D/CubicBezier/Basic.hs +++ b/Geom2D/CubicBezier/Basic.hs @@ -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]) + diff --git a/Geom2D/CubicBezier/Intersection.hs b/Geom2D/CubicBezier/Intersection.hs --- a/Geom2D/CubicBezier/Intersection.hs +++ b/Geom2D/CubicBezier/Intersection.hs @@ -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 diff --git a/Geom2D/CubicBezier/MetaPath.hs b/Geom2D/CubicBezier/MetaPath.hs --- a/Geom2D/CubicBezier/MetaPath.hs +++ b/Geom2D/CubicBezier/MetaPath.hs @@ -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) diff --git a/Geom2D/CubicBezier/Numeric.hs b/Geom2D/CubicBezier/Numeric.hs --- a/Geom2D/CubicBezier/Numeric.hs +++ b/Geom2D/CubicBezier/Numeric.hs @@ -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: -- diff --git a/Math/BernsteinPoly.hs b/Math/BernsteinPoly.hs --- a/Math/BernsteinPoly.hs +++ b/Math/BernsteinPoly.hs @@ -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] diff --git a/cubicbezier.cabal b/cubicbezier.cabal --- a/cubicbezier.cabal +++ b/cubicbezier.cabal @@ -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 diff --git a/tests/test.hs b/tests/test.hs new file mode 100644 --- /dev/null +++ b/tests/test.hs @@ -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