cubicbezier 0.5.0.0 → 0.6.0.0
raw patch · 10 files changed
+1233/−476 lines, 10 filesdep +fast-mathdep +vector-spacedep ~basedep ~mtl
Dependencies added: fast-math, vector-space
Dependency ranges changed: base, mtl
Files
- Geom2D.hs +53/−41
- Geom2D/CubicBezier/Approximate.hs +81/−44
- Geom2D/CubicBezier/Basic.hs +232/−24
- Geom2D/CubicBezier/Intersection.hs +63/−78
- Geom2D/CubicBezier/MetaPath.hs +26/−22
- Geom2D/CubicBezier/Numeric.hs +24/−1
- Geom2D/CubicBezier/Overlap.lhs +415/−259
- Geom2D/CubicBezier/Stroke.hs +323/−0
- Math/BernsteinPoly.hs +3/−3
- cubicbezier.cabal +13/−4
Geom2D.hs view
@@ -1,21 +1,24 @@-{-# LANGUAGE TypeFamilies, MultiParamTypeClasses, FlexibleInstances, DeriveFunctor, FunctionalDependencies #-}+{-# LANGUAGE TypeFamilies, MultiParamTypeClasses, FlexibleInstances, DeriveTraversable, FunctionalDependencies #-} -- | Basic 2 dimensional geometry functions.-module Geom2D where+module Geom2D (+ module Data.VectorSpace,+ module Data.Cross,+ module Geom2D ) where import qualified Data.Vector.Generic.Mutable as M import qualified Data.Vector.Generic as G import qualified Data.Vector.Unboxed as V+import Data.VectorSpace+import Data.Cross import Control.Monad-+import Numeric.FastMath -infixl 6 ^+^, ^-^-infixl 7 *^, ^*, ^/ infixr 5 $* data Point a = Point { pointX :: !a,- pointY :: !a}- deriving (Eq, Functor)+ pointY :: !a+ } deriving (Eq, Ord, Functor, Foldable, Traversable) type DPoint = Point Double @@ -31,19 +34,20 @@ xformD :: !a, xformE :: !a, xformF :: !a }- deriving (Eq, Show, Functor)+ deriving (Eq, Show, Functor, Foldable, Traversable) data Line a = Line (Point a) (Point a)- deriving (Show, Eq, Functor)+ deriving (Show, Eq, Functor, Foldable, Traversable)+ data Polygon a = Polygon [Point a]- deriving (Show, Eq, Functor)+ deriving (Show, Eq, Functor, Foldable, Traversable) class AffineTransform a b | a -> b where transform :: Transform b -> a -> a-+ instance Num a => AffineTransform (Transform a) a where {-# INLINE transform #-}- transform (Transform a' b' c' d' e' f') (Transform a b c d e f) =+ transform (Transform a' b' c' d' e' f') (Transform a b c d e f) = Transform (a*a'+b'*d) (a'*b + b'*e) (a'*c+b'*f +c') (d'*a+e'*d) (d'*b+e'*e) (d'*c+e'*f+f') @@ -72,6 +76,7 @@ {-# INLINE basicSet #-} {-# INLINE basicUnsafeCopy #-} {-# INLINE basicUnsafeGrow #-}+ basicInitialize (MV_Point v) = M.basicInitialize v basicLength (MV_Point v) = M.basicLength v basicUnsafeSlice i n (MV_Point v) = MV_Point $ M.basicUnsafeSlice i n v basicOverlaps (MV_Point v1) (MV_Point v2) = M.basicOverlaps v1 v2@@ -108,7 +113,7 @@ {-# INLINE ($*) #-} -- | Calculate the inverse of a transformation.-inverse :: (Eq a, Num a, Fractional a) => Transform a -> Maybe (Transform a)+inverse :: (Eq a, Fractional a) => Transform a -> Maybe (Transform a) inverse (Transform a b c d e f) = case a*e - b*d of 0 -> Nothing det -> Just $! Transform (a/det) (d/det) (-(a*c + d*f)/det) (b/det) (e/det)@@ -117,8 +122,9 @@ -- | Return the parameters (a, b, c) for the normalised equation -- of the line: @a*x + b*y + c = 0@.-lineEquation :: Floating t => Line t -> (t, t, t)-lineEquation (Line (Point x1 y1) (Point x2 y2)) = (a, b, c)+lineEquation :: Floating t => Line t -> ( t, t, t )+lineEquation (Line (Point x1 y1) (Point x2 y2)) =+ a `seq` b `seq` c `seq` (a, b, c) where a = a' / d b = b' / d c = -(y1*b' + x1*a') / d@@ -130,9 +136,10 @@ -- | 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 :: Floating a => Line a -> Point a -> a-lineDistance l = \(Point x y) -> a*x + b*y + c- where (a, b, c) = lineEquation l-{-# SPECIALIZE lineDistance :: Line Double -> DPoint -> Double #-}+lineDistance l = + case lineEquation l of+ (a, b, c) -> \(Point x y) -> a*x + b*y + c+{-# INLINE lineDistance #-} -- | Return the point on the line closest to the given point. closestPoint :: Fractional a => Line a -> Point a -> Point a@@ -165,6 +172,12 @@ vectorMag (Point x y) = sqrt(x*x + y*y) {-# INLINE vectorMag #-} +-- | The lenght of the vector.+vectorMagSquare :: Floating a => Point a -> a+vectorMagSquare (Point x y) = x*x + y*y+{-# INLINE vectorMagSquare #-}++ -- | The angle of the vector, in the range @(-'pi', 'pi']@. vectorAngle :: RealFloat a => Point a -> a vectorAngle (Point 0.0 0.0) = 0.0@@ -182,30 +195,24 @@ where l = vectorMag p {-# INLINE normVector #-} --- | Scale vector by constant.-(*^) :: Num a => a -> Point a -> Point a-s *^ (Point x y) = Point (s*x) (s*y)-{-# INLINE (*^) #-}---- | Scale vector by reciprocal of constant.-(^/) :: Fractional a => Point a -> a -> Point a-(Point x y) ^/ s = Point (x/s) (y/s)-{-# INLINE (^/) #-}+instance Num e => AdditiveGroup (Point e) where+ zeroV = Point 0 0+ {-# INLINE (^+^) #-}+ (Point x1 y1) ^+^ (Point x2 y2) = Point (x1+x2) (y1+y2)+ {-# INLINE negateV #-}+ negateV (Point a b) = Point (-a) (-b)+ {-# INLINE (^-^) #-}+ (Point x1 y1) ^-^ (Point x2 y2) = Point (x1-x2) (y1-y2) --- | Scale vector by constant, with the arguments swapped.-(^*) :: Num a => Point a -> a -> Point a-p ^* s = s *^ p-{-# INLINE (^*) #-}+instance (Num e) => VectorSpace (Point e) where+ type Scalar (Point e) = e+ s *^ (Point x y) = Point (s*x) (s*y) --- | Add two vectors.-(^+^) :: Num a => Point a -> Point a -> Point a-(Point x1 y1) ^+^ (Point x2 y2) = Point (x1+x2) (y1+y2)-{-# INLINE (^+^) #-}+instance (AdditiveGroup e, Num e) => InnerSpace (Point e) where+ (<.>) = (^.^) --- | Subtract two vectors.-(^-^) :: Num a => Point a -> Point a -> Point a-(Point x1 y1) ^-^ (Point x2 y2) = Point (x1-x2) (y1-y2)-{-# INLINE (^-^) #-}+instance (Floating e) => HasNormal (Point e) where+ normalVec = normVector -- | Dot product of two vectors. (^.^) :: Num a => Point a -> Point a -> a@@ -238,6 +245,10 @@ flipVector (Point x y) = Point x (-y) {-# INLINE flipVector #-} +turnAround :: (Num a) => Point a -> Point a+turnAround = negateV+{-# INLINE turnAround #-}+ -- | Create a transform that rotates by the angle of the given vector -- with the x-axis rotateVec :: Floating a => Point a -> Transform a@@ -253,12 +264,12 @@ -- | Rotate vector 90 degrees left. rotate90L :: Floating s => Transform s-rotate90L = rotateVec (Point 0 1)+rotate90L = Transform 0 (-1) 0 1 0 0 {-# INLINE rotate90L #-} -- | Rotate vector 90 degrees right. rotate90R :: Floating s => Transform s-rotate90R = rotateVec (Point 0 (-1))+rotate90R = Transform 0 1 0 (-1) 0 0 {-# INLINE rotate90R #-} -- | Create a transform that translates by the given vector.@@ -270,3 +281,4 @@ idTrans :: Num a => Transform a idTrans = Transform 1 0 0 0 1 0 {-# INLINE idTrans #-}+
Geom2D/CubicBezier/Approximate.hs view
@@ -1,6 +1,7 @@ {-# LANGUAGE BangPatterns, MultiWayIf #-} module Geom2D.CubicBezier.Approximate- (approximatePath, approximateQuadPath, approximatePathMax, approximateCubic)+ (approximatePath, approximateQuadPath, approximatePathMax, approximateQuadPathMax,+ approximateCubic) where import Geom2D import Geom2D.CubicBezier.Basic@@ -60,28 +61,38 @@ | otherwise = approximateQuad' f tol tmin tmax fast where curve = approx1quad f tmin tmax- err = maxDist f curve tmin tmax+ err = maxDist f curve tmin tmax (if fast then 0 else 5) {-# SPECIALIZE approximateQuadPath :: (Double -> (DPoint, DPoint)) -> Double -> Double -> Double -> Bool -> [QuadBezier Double] #-} -- find the distance between the function at t and the quadratic bezier. -- calculate the value and derivative at t, and improve the closeness of t.-quadDist :: (V.Unbox a, Floating a) =>- (a -> (Point a, Point a)) -> QuadBezier a -> a -> a -> a -> a-quadDist f qb tmin tmax t =- let p = fst (f $ interpolate tmin tmax t)- (b, b') = evalBezierDeriv qb t- -- distance from p to the normal at b(t) / velocity- nd = ((p ^-^ b) ^.^ b') / (b'^.^b')- in vectorDistance p $ evalBezier qb (t + nd)+quadDist :: (V.Unbox a, Ord a, Floating a) =>+ (a -> (Point a, Point a)) -> QuadBezier a -> a+ -> a -> Int -> a -> a+quadDist f qb tmin tmax maxiter t =+ quadDist' f qb tmin tmax t (fst (f $ interpolate tmin tmax t)) maxiter (evalBezierDeriv qb t) +quadDist' :: (V.Unbox a, Ord a, Floating a) =>+ (a -> (Point a, Point a))+ -> QuadBezier a -> a -> a -> a -> Point a -> Int -> (Point a, Point a) -> a+quadDist' f qb tmin tmax t p maxiter (b, b')+ | maxiter <= 1 || abs (err2-err1) <= err1 * (1/8) = err2+ | otherwise = quadDist' f qb tmin tmax (t+ndist) p (maxiter-1) (b2, b2')+ where dp = p ^-^ b+ err1 = vectorMag dp+ -- distance from p to the normal at b(t) / velocity at t+ ndist = (dp ^.^ b') / (b' ^.^ b')+ (b2, b2') = evalBezierDeriv qb (t+ndist)+ err2 = vectorDistance p b2+ -- 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 =- quadDist f qb tmin tmax $- goldSearch (quadDist f qb tmin tmax) 4+ QuadBezier a -> a -> a -> Int -> a+maxDist f qb tmin tmax maxiter =+ quadDist f qb tmin tmax maxiter $+ goldSearch (quadDist f qb tmin tmax maxiter) 4 approxquad :: (Ord a, Floating a) => Point a -> Point a -> Point a -> Point a -> QuadBezier a@@ -104,18 +115,24 @@ splitQuad :: (Show a, V.Unbox a, Ord a, Floating a) => a -> a -> (a -> (Point a, Point a)) -> a -> a -> Int -> (a, a, QuadBezier a, a, QuadBezier a)-splitQuad node offset f tmin tmax maxiter+splitQuad node offset f tmin tmax maxiter =+ splitQuad' node offset f tmin tmax maxiter maxiter+ +splitQuad' :: (Show a, V.Unbox a, Ord a, Floating a) =>+ a -> a -> (a -> (Point a, Point a))+ -> a -> a -> Int -> Int -> (a, a, QuadBezier a, a, QuadBezier a)+splitQuad' node offset f tmin tmax maxiter maxiter2 | maxiter < 1 || (err0 < 2*err1 && err0 > err1/2) = (tmid, err0, curve0, err1, curve1) | otherwise =- splitQuad (if err0 < err1 then node+offset else node-offset)- (offset/2) f tmin tmax (maxiter-1)+ splitQuad' (if err0 < err1 then node+offset else node-offset)+ (offset/2) f tmin tmax (maxiter-1) maxiter2 where tmid = interpolate tmin tmax node curve0 = approx1quad f tmin tmid - err0 = maxDist f curve0 tmin tmid+ err0 = maxDist f curve0 tmin tmid maxiter2 curve1 = approx1quad f tmid tmax - err1 = maxDist f curve1 tmid tmax+ err1 = maxDist f curve1 tmid tmax maxiter2 approximateQuad' :: (Show a, V.Unbox a, Ord a, Floating a) => (a -> (Point a, Point a)) -> @@ -130,7 +147,7 @@ else approximateQuad' f tol tmid tmax fast) where (tmid, err0, curve0, err1, curve1) =- splitQuad 0.5 0.25 f tmin tmax (if fast then 0 else 5)+ splitQuad 0.5 0.25 f tmin tmax (if fast then 0 else 5) approximatePath' :: (V.Unbox a, Ord a, Floating a) => (a -> (Point a, Point a)) -> Int ->@@ -145,7 +162,7 @@ else approximatePath' f n tol tmid tmax fast) where (tmid, err0, curve0, err1, curve1) =- splitCubic 0.5 0.25 n f tmin tmax (if fast then 0 else 5)+ splitCubic n 0.5 0.25 f tmin tmax (if fast then 0 else 5) --{-# SPECIALIZE approximatePath' :: (Double -> (Point Double, Point Double)) -> Int -> Double -> Double -> Double -> [CubicBezier Double] #-} -- | Like approximatePath, but limit the number of subcurves.@@ -162,16 +179,17 @@ -> a -- ^ The upper parameter of the function -> Bool -- ^ Calculate the result faster, but with more- -- subcurves. Runs typically 10 times faster, but- -- generates 50% more subcurves. Useful for interactive use.+ -- subcurves. Runs faster (typically 10 times),+ -- but generates more subcurves (about 50%).+ -- Useful for interactive use. -> [CubicBezier a]-approximatePathMax m f n tol tmin tmax fast =- approxMax f tol m ts fast segments+approximatePathMax m f samples tol tmin tmax fast =+ approxMax f tol m fast (splitCubic samples) segments where segments = M.singleton err (FunctionSegment tmin tmax outline) (p0, p0') = f tmin (p1, p1') = f tmax- ts = V.map (\i -> fromIntegral i/(fromIntegral n+1) `asTypeOf` tmin) $- V.enumFromN (1::Int) n+ ts = V.map (\i -> fromIntegral i/(fromIntegral samples+1) `asTypeOf` tmin) $+ V.enumFromN (1::Int) samples points = V.map (fst . f . interpolate tmin tmax) ts curveCb = CubicBezier p0 (p0^+^p0') (p1^-^p1') p1 (outline, err) =@@ -179,24 +197,47 @@ {-# SPECIALIZE approximatePathMax :: Int -> (Double -> (Point Double, Point Double)) -> Int -> Double -> Double -> Double -> Bool -> [CubicBezier Double] #-}-data FunctionSegment a = FunctionSegment {+data FunctionSegment b a = FunctionSegment { fsTmin :: !a, -- the least t param of the segment in the original curve _fsTmax :: !a, -- the max t param of the segment in the original curve- fsCurve :: CubicBezier a -- the curve segment+ fsCurve :: b -- the curve segment } +-- | Like approximateQuadPath, but limit the number of subcurves.+approximateQuadPathMax :: (V.Unbox a, Show a, Floating a, Ord a) =>+ Int -- ^ The maximum number of subcurves+ -> (a -> (Point a, Point a)) -- ^ The function to approximate and it's derivative+ -> a -- ^ The tolerance+ -> a -- ^ The lower parameter of the function + -> a -- ^ The upper parameter of the function+ -> Bool+ -- ^ Calculate the result faster, but with more+ -- subcurves. Runs faster, but generates more+ -- subcurves. Useful for interactive use.+ -> [QuadBezier a]+approximateQuadPathMax m f tol tmin tmax fast =+ approxMax f tol m fast splitQuad segments+ where segments = M.singleton err (FunctionSegment tmin tmax curveQd)+ (p0, p0') = f tmin+ (p1, p1') = f tmax+ curveQd = approxquad p0 p0' p1' p1+ err = maxDist f curveQd tmin tmax (if fast then 0 else 5)+{-# SPECIALIZE approximateQuadPathMax ::+ Int -> (Double -> (Point Double, Point Double))+ -> Double -> Double -> Double -> Bool -> [QuadBezier Double] #-} -- Keep a map from maxError to FunctionSegment for each subsegment to keep -- track of the segment with the maximum error. This ensures a n -- log(n) execution time, rather than n^2 when a list is used. approxMax :: (V.Unbox a, Ord a, Floating a) =>- (a -> (Point a, Point a)) -> a -> Int- -> V.Vector a -> Bool -> M.Map a (FunctionSegment a) ->- [CubicBezier a]-approxMax f tol n ts fast segments- | (n <= 1) || (err < tol) =+ (a -> (Point a, Point a)) -> a -> Int -> Bool+ -> (a -> a -> (a -> (Point a, Point a)) -> a -> a -> Int -> (a, a, b, a, b))+ -> M.Map a (FunctionSegment b a)+ -> [b]+approxMax f tol maxCurves fast splitBez segments+ | (maxCurves <= 1) || (err < tol) = map fsCurve $ sortBy (compare `on` fsTmin) $ map snd $ M.toList segments- | otherwise = approxMax f tol (n-1) ts fast $+ | otherwise = approxMax f tol (maxCurves-1) fast splitBez $ M.insert err_l (FunctionSegment t_min t_mid curve_l) $ M.insert err_r (FunctionSegment t_mid t_max curve_r) newSegments@@ -204,25 +245,21 @@ ((err, FunctionSegment t_min t_max _), newSegments) = M.deleteFindMax segments (t_mid, err_l, curve_l, err_r, curve_r) =- splitCubic 0.5 0.25 n f t_min t_max (if fast then 0 else 5)-{-# SPECIALIZE approxMax :: (Double -> (Point Double, Point Double)) -> Double -> Int- -> V.Vector Double -> Bool -> M.Map Double (FunctionSegment Double) -> [CubicBezier Double] #-}+ splitBez 0.5 0.25 f t_min t_max (if fast then 0 else 5) splitCubic :: (V.Unbox a, Ord a, Floating a) =>- a -> a -> Int -> (a -> (Point a, Point a))+ Int -> a -> a -> (a -> (Point a, Point a)) -> a -> a -> Int -> (a, a, CubicBezier a, a, CubicBezier a)-splitCubic node offset n f tmin tmax maxiter+splitCubic n node offset f tmin tmax maxiter | maxiter < 1 || (err0 < 2*err1 && err0 > err1/2) = (tmid, err0, curve0, err1, curve1) | otherwise = - splitCubic (if err0 < err1 then node+offset else node-offset)- (offset/2) n f tmin tmax (maxiter-1)+ splitCubic n (if err0 < err1 then node+offset else node-offset)+ (offset/2) f tmin tmax (maxiter-1) where tmid = interpolate tmin tmax node (curve0, err0) = approx1cubic n f tmin tmid maxiter (curve1, err1) = approx1cubic n f tmid tmax maxiter-{-# SPECIALIZE splitCubic :: Double -> Double -> Int -> (Double -> (Point Double, Point Double))- -> Double -> Double -> Int -> (Double, Double, CubicBezier Double, Double, CubicBezier Double) #-} approx1cubic :: (V.Unbox a, Ord a, Floating a) => Int -> (a -> (Point a, Point a)) -> a -> a ->
Geom2D/CubicBezier/Basic.hs view
@@ -1,18 +1,26 @@-{-# LANGUAGE BangPatterns, FlexibleInstances, MultiParamTypeClasses, DeriveFunctor, ViewPatterns #-}+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE BangPatterns, FlexibleInstances, MultiParamTypeClasses, DeriveTraversable, ViewPatterns, PatternSynonyms, MultiWayIf #-} module Geom2D.CubicBezier.Basic (CubicBezier (..), QuadBezier (..), AnyBezier (..), GenericBezier(..), PathJoin (..), ClosedPath(..), OpenPath (..), AffineTransform (..), anyToCubic, anyToQuad, openPathCurves, closedPathCurves, curvesToOpen, curvesToClosed,+ consOpenPath, consClosedPath, openClosedPath, closeOpenPath, bezierParam, bezierParamTolerance, reorient, bezierToBernstein, evalBezierDerivs, evalBezier, evalBezierDeriv, findBezierTangent, quadToCubic, bezierHoriz, bezierVert, findBezierInflection, findBezierCusp,- bezierArc, arcLength, arcLengthParam, splitBezier, bezierSubsegment, splitBezierN,- colinear)+ bezierArc, arcLength, arcLengthParam, splitBezier, bezierSubsegment,+ splitBezierN, colinear, closest, findX) where import Geom2D import Geom2D.CubicBezier.Numeric import Math.BernsteinPoly import Numeric.Integration.TanhSinh+import Data.Monoid ()+import Data.List (minimumBy)+import Data.Function (on)+import Data.VectorSpace+import Debug.Trace+ import qualified Data.Vector.Unboxed as V import qualified Data.Vector.Unboxed.Mutable as MV @@ -22,14 +30,14 @@ cubicC1 :: !(Point a), cubicC2 :: !(Point a), cubicC3 :: !(Point a)}- deriving (Eq, Show, Functor)+ deriving (Eq, Show, Functor, Foldable, Traversable) -- | 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, Foldable, Traversable) -- Use a tuple, because it has 0(1) unzip when using unboxed vectors. -- | A bezier curve of any degree.@@ -79,26 +87,61 @@ data PathJoin a = JoinLine | JoinCurve (Point a) (Point a)- deriving (Show, Functor)+ deriving (Show, Functor, Foldable, Traversable) data OpenPath a = OpenPath [(Point a, PathJoin a)] (Point a) - deriving (Show, Functor)+ deriving (Show, Functor, Foldable, Traversable) data ClosedPath a = ClosedPath [(Point a, PathJoin a)]- deriving (Show, Functor)+ deriving (Show, Functor, Foldable, Traversable) +instance Monoid (OpenPath a) where+ mempty = OpenPath [] (error "empty path")+ mappend p1 (OpenPath [] _) = p1+ mappend (OpenPath [] _) p2 = p2+ mappend (OpenPath joins1 _) (OpenPath joins2 p) =+ OpenPath (joins1 ++ joins2) p++instance (Num a) => AffineTransform (PathJoin a) a where+ transform _ JoinLine = JoinLine+ transform t (JoinCurve p q) = JoinCurve (transform t p) (transform t q)++instance (Num a) => AffineTransform (OpenPath a) a where+ transform t (OpenPath s p) = OpenPath (map (transformSeg t) s) (transform t p)++transformSeg :: (Num a) => Transform a -> (Point a, PathJoin a) -> (Point a, PathJoin a)+transformSeg t (p, jn) = (transform t p, transform t jn)++instance (Num a) => AffineTransform (ClosedPath a) a where+ transform t (ClosedPath s) = ClosedPath (map (transformSeg t) s)+ instance (Num a) => AffineTransform (CubicBezier a) a where {-# SPECIALIZE transform :: Transform Double -> CubicBezier Double -> CubicBezier Double #-} transform t (CubicBezier c0 c1 c2 c3) = CubicBezier (transform t c0) (transform t c1) (transform t c2) (transform t c3) +instance (Num a) => AffineTransform (QuadBezier a) a where+ {-# SPECIALIZE transform :: Transform Double -> QuadBezier Double -> QuadBezier Double #-}+ transform t (QuadBezier c0 c1 c2) =+ QuadBezier (transform t c0) (transform t c1) (transform t c2)++-- | construct an open path+consOpenPath :: Point a -> PathJoin a -> OpenPath a -> OpenPath a+consOpenPath p join (OpenPath joins q) =+ OpenPath ((p, join):joins) q++-- | construct a closed path+consClosedPath :: Point a -> PathJoin a -> ClosedPath a -> ClosedPath a+consClosedPath p join (ClosedPath joins) =+ ClosedPath ((p, join):joins)+ -- | Return the open path as a list of curves. openPathCurves :: Fractional a => OpenPath a -> [CubicBezier a] openPathCurves (OpenPath curves p) = go curves p where go [] _ = []- go [(p0, jn)] p = [makeCB p0 jn p]- go ((p0, jn):rest@((p1,_):_)) p =- makeCB p0 jn p1 : go rest p- makeCB p0 (JoinLine) p1 =+ go [(p0, jn)] q = [makeCB p0 jn q]+ go ((p0, jn):rest@((p1,_):_)) q =+ makeCB p0 jn p1 : go rest q+ makeCB p0 JoinLine p1 = CubicBezier p0 (interpolateVector p0 p1 (1/3)) (interpolateVector p0 p1 (2/3)) p1 makeCB p0 (JoinCurve p1 p2) p3 =@@ -128,7 +171,17 @@ where OpenPath cs2 _ = curvesToOpen cs +-- | close an open path, discarding the last point+closeOpenPath :: OpenPath a -> ClosedPath a+closeOpenPath (OpenPath j p) = ClosedPath j +-- | open a closed path+openClosedPath :: ClosedPath a -> OpenPath a+openClosedPath (ClosedPath []) = OpenPath [] (error "empty path")+openClosedPath (ClosedPath j@((p,_):_)) = OpenPath j p+++ -- | safely convert from `AnyBezier' to `CubicBezier` anyToCubic :: (V.Unbox a) => AnyBezier a -> Maybe (CubicBezier a) anyToCubic b@(AnyBezier v)@@ -141,9 +194,10 @@ | degree b == 2 = Just $ unsafeFromVector v | otherwise = Nothing -evalBezierDerivsCubic :: Fractional a =>+evalBezierDerivsCubic :: Num a => CubicBezier a -> a -> [Point a] evalBezierDerivsCubic (CubicBezier a b c d) t =+ p `seq` p' `seq` p'' `seq` p''' `seq` [p, p', p'', p''', Point 0 0] where u = 1-t@@ -152,23 +206,47 @@ da = 3*^(b^-^a) db = 3*^(c^-^b) dc = 3*^(d^-^c)- p = u*^(u*^(u*^a ^+^ 3*t*^b) ^+^ 3*t2*^c) ^+^ t3*^d- p' = u*^(u*^da ^+^ 2*t*^db) ^+^ t2*^dc- p'' = 2*u*^(db^-^da) ^+^ 2*t*^(dc^-^db)+ p = u*^(u*^(u*^a ^+^ (3*t)*^b) ^+^ (3*t2)*^c) ^+^ t3*^d+ p' = u*^(u*^da ^+^ (2*t)*^db) ^+^ t2*^dc+ p'' = (2*u)*^(db^-^da) ^+^ (2*t)*^(dc^-^db) p''' = 2*^(dc^-^2*^db^+^da) {-# SPECIALIZE evalBezierDerivsCubic :: CubicBezier Double -> Double -> [DPoint] #-} -evalBezierDerivsQuad :: Fractional a =>+evalBezierDerivCubic :: Num a =>+ CubicBezier a -> a -> (Point a, Point a)+evalBezierDerivCubic (CubicBezier a b c d) t = p `seq` p' `seq` (p, p')+ where+ u = 1-t+ t2 = t*t+ t3 = t2*t+ da = 3*^(b^-^a)+ db = 3*^(c^-^b)+ dc = 3*^(d^-^c)+ p = u*^(u*^(u*^a ^+^ (3*t)*^b) ^+^ (3*t2)*^c) ^+^ t3*^d+ p' = u*^(u*^da ^+^ (2*t)*^db) ^+^ t2*^dc+{-# SPECIALIZE evalBezierDerivCubic :: CubicBezier Double -> Double -> (DPoint, DPoint) #-} ++evalBezierDerivsQuad :: Num a => QuadBezier a -> a -> [Point a] evalBezierDerivsQuad (QuadBezier a b c) t = [p, p', p'', Point 0 0] where u = 1-t t2 = t*t- p = u*^(u*^a ^+^ 2*t*^b) ^+^ t2*^c+ p = u*^(u*^a ^+^ (2*t)*^b) ^+^ t2*^c p' = 2*^(u*^(b^-^a) ^+^ t*^(c^-^b)) p'' = 2*^(c^-^ 2*^b ^+^ a) {-# SPECIALIZE evalBezierDerivsQuad :: QuadBezier Double -> Double -> [DPoint] #-} +evalBezierDerivQuad :: Num a =>+ QuadBezier a -> a -> (Point a, Point a)+evalBezierDerivQuad (QuadBezier a b c) t = p `seq` p' `seq` (p, p')+ where+ u = 1-t+ t2 = t*t+ p = u*^(u*^a ^+^ (2*t)*^b) ^+^ t2*^c+ p' = 2*^(u*^(b^-^a) ^+^ t*^(c^-^b))+{-# SPECIALIZE evalBezierDerivQuad :: QuadBezier Double -> Double -> (DPoint, DPoint) #-} + -- | Evaluate the bezier and all its derivatives using the modified horner algorithm. evalBezierDerivs :: (GenericBezier b, V.Unbox a, Fractional a) => b a -> a -> [Point a]@@ -193,7 +271,23 @@ where maxVel = 3 * V.maximum (V.zipWith vectorDistance (V.map (uncurry Point) v) (V.map (uncurry Point) $ V.tail v))+{-# NOINLINE [2] bezierParamTolerance #-} +bezierParamToleranceCubic :: (Ord a, Floating a) => CubicBezier a -> a -> a+bezierParamToleranceCubic (CubicBezier p0 p1 p2 p3) eps = eps / maxVel+ where maxVel = (3 *) $ max (vectorDistance p0 p1) $+ max (vectorDistance p1 p2)+ (vectorDistance p2 p3)+{-# SPECIALIZE bezierParamToleranceCubic :: CubicBezier Double -> Double -> Double #-}+{-# RULES "bezierParamTolerance/cubic" bezierParamTolerance = bezierParamToleranceCubic #-}++bezierParamToleranceQuad :: (Ord a, Floating a) => QuadBezier a -> a -> a+bezierParamToleranceQuad (QuadBezier p0 p1 p2) eps = eps / maxVel+ where maxVel = (3 *) $ max (vectorDistance p0 p1) (vectorDistance p1 p2)++{-# SPECIALIZE bezierParamToleranceQuad :: QuadBezier Double -> Double -> Double #-}+{-# RULES "bezierParamTolerance/quad" bezierParamTolerance = bezierParamToleranceQuad #-}+ -- | Reorient to the curve B(1-t). reorient :: (GenericBezier b, V.Unbox a) => b a -> b a reorient = unsafeFromVector . V.reverse . toVector@@ -224,7 +318,7 @@ evalBezierCubic :: Fractional a => CubicBezier a -> a -> Point a evalBezierCubic (CubicBezier a b c d) t =- u*^(u*^(u*^a ^+^ 3*t*^b) ^+^ 3*t2*^c) ^+^ t3*^d+ u*^(u*^(u*^a ^+^ (3*t)*^b) ^+^ (3*t2)*^c) ^+^ t3*^d where u = 1-t t2 = t*t@@ -234,7 +328,7 @@ evalBezierQuad :: Fractional a => QuadBezier a -> a -> Point a evalBezierQuad (QuadBezier a b c) t = - u*^(u*^a ^+^ 2*t*^b) ^+^ t2*^c+ u*^(u*^a ^+^ (2*t)*^b) ^+^ t2*^c where u = 1-t t2 = t*t@@ -249,6 +343,9 @@ evalBezierDeriv bc t = (b,b') where (b:b':_) = evalBezierDerivs bc t+{-# NOINLINE evalBezierDeriv #-} +{-# RULES "evalBezierDeriv/cubic" evalBezierDeriv = evalBezierDerivCubic #-}+{-# RULES "evalBezierDeriv/quad" evalBezierDeriv = evalBezierDerivQuad #-} -- | @findBezierTangent p b@ finds the parameters where -- the tangent of the bezier curve @b@ has the same direction as vector p.@@ -364,7 +461,7 @@ quadToCubic :: (Fractional a) => QuadBezier a -> CubicBezier a quadToCubic (QuadBezier a b c) =- CubicBezier a (1/3*^(a ^+^ 2*^b)) (1/3*^(2*^b ^+^ c)) c+ CubicBezier a ((1/3)*^(a ^+^ 2*^b)) ((1/3)*^(2*^b ^+^ c)) c -- | Split a bezier curve into two curves. splitBezier :: (V.Unbox a, Fractional a) =>@@ -388,7 +485,7 @@ abbc = interpolateVector ab bc t bccd = interpolateVector bc cd t mid = interpolateVector abbc bccd t- in (CubicBezier a ab abbc mid, CubicBezier mid bccd cd d)+ in mid `seq` (CubicBezier a ab abbc mid, CubicBezier mid bccd cd d) {-# SPECIALIZE splitBezierCubic :: CubicBezier Double -> Double -> (CubicBezier Double, CubicBezier Double) #-} -- | Split a bezier curve into two curves.@@ -402,7 +499,6 @@ {-# RULES "splitBezier/cubic" splitBezier = splitBezierCubic #-} {-# RULES "splitBezier/quad" splitBezier = splitBezierQuad #-} - -- | Return the subsegment between the two parameters. bezierSubsegment :: (Ord a, V.Unbox a, Fractional a) => GenericBezier b => b a -> a -> a -> b a@@ -411,9 +507,20 @@ | t2 == 0 = fst $ splitBezier b t1 | otherwise = snd $ flip splitBezier (t1/t2) $ fst $ splitBezier b t2-{-# SPECIALIZE bezierSubsegment :: CubicBezier Double -> Double -> Double -> CubicBezier Double #-} {-# SPECIALIZE bezierSubsegment :: QuadBezier Double -> Double -> Double -> QuadBezier Double #-}+{-# NOINLINE[2] bezierSubsegment #-} +-- | Return the subsegment between the two parameters.+bezierSubsegmentCubic :: (V.Unbox a, Fractional a) => CubicBezier a -> a -> a -> CubicBezier a+bezierSubsegmentCubic b t1 t2 =+ CubicBezier b1 (b1 ^+^ b1' ^* ((t2-t1)/3))+ (b2 ^-^ (b2' ^* ((t2-t1)/3))) b2+ where (b1, b1') = evalBezierDeriv b t1+ (b2, b2') = evalBezierDeriv b t2+{-# SPECIALIZE bezierSubsegmentCubic :: CubicBezier Double -> Double -> Double -> CubicBezier Double #-}+{-# RULES "bezierSubsegment/cubic" bezierSubsegment = bezierSubsegmentCubic #-}++ -- | Split a bezier curve into a list of beziers -- The parameters should be in ascending order or -- the result is unpredictable.@@ -429,6 +536,107 @@ {-# SPECIALIZE splitBezierN :: CubicBezier Double -> [Double] -> [CubicBezier Double] #-} {-# SPECIALIZE splitBezierN :: QuadBezier Double -> [Double] -> [QuadBezier Double] #-} +evalBezierDerivs2Cubic :: CubicBezier Double -> Double -> (DPoint, DPoint, DPoint)+evalBezierDerivs2Cubic (CubicBezier a b c d) t =+ p `seq` p' `seq` p'' `seq`(p, p', p'')+ where+ u = 1-t+ t2 = t*t+ t3 = t2*t+ da = 3*^(b^-^a)+ db = 3*^(c^-^b)+ dc = 3*^(d^-^c)+ p = u*^(u*^(u*^a ^+^ (3*t)*^b) ^+^ (3*t2)*^c) ^+^ t3*^d+ p' = u*^(u*^da ^+^ (2*t)*^db) ^+^ t2*^dc+ p'' = (2*u)*^(db^-^da) ^+^ (2*t)*^(dc^-^db)++-- estimate the next approximation for the point closest to p by+-- dividing the approximate arclength on the osculating circle at t by the+-- speed at t.+closestBezierCurve :: CubicBezier Double -> DPoint -> Double -> Double+closestBezierCurve cb p t+ | vectorMagSquare (p ^-^ b) < v*v*r_v*r_v*10 =+ closestLinePt b b' p+ | otherwise =+ -- circular arc divided by velocity+ (fastVectorAngle (rotateScaleVec (flipVector (p ^-^ c)) $* (Point (-by') bx'))) * r_v+ where+ closestLinePt :: DPoint -> DPoint -> DPoint -> Double+ closestLinePt (Point bx by) (Point bx' by') (Point px py) =+ ((px - bx)*bx' + (py - by)*by')/v+ (b, b'@(Point bx' by'), Point bx'' by'') = evalBezierDerivs2Cubic cb t+ -- radius of curvature / velocity+ r_v = v/denom+ v = bx'*bx' + by'*by'+ denom = bx''*by' - bx'*by''+ -- center of osculating circle+ c = b ^+^ ((rotate90R $* b') ^* r_v)++-- approximation of vectorAngle, fastVectorAngle θ ≊ vectorAngle θ for+-- small |θ|. This avoids a costly atan2 instruction, and I measured a+-- performance increase upto 33%.+fastVectorAngle :: DPoint -> Double+fastVectorAngle (Point x y)+ | abs y < abs x = y/x + if x < 0 then signum y*pi else 0+ | otherwise = signum y*pi/2-(x/y)+ +-- | Find the closest value on the bezier to the given point, within tolerance. Return the first value found.+closest :: CubicBezier Double -> DPoint -> Double -> Double+closest cb p eps =+ iter (closestBezierCurve cb p) eps 0 1 0.5++-- iterate, fallback to bisection if the approximation doesn't converge. +iter :: (Double -> Double) -> Double -> Double -> Double -> Double -> Double +iter f eps mint maxt curt+ | abs dt <= eps = curt + dt+ | dt < 0 =+ if | curt + dt <= mint ->+ if mint == 0 && f 0 <= 0+ then 0+ -- bisect+ else+ let dT = (curt - mint)/2+ in if dT < eps then mint+dT+ else iter f eps mint curt (mint+dT)+ | otherwise ->+ iter f eps mint curt (curt+dt)+ | otherwise =+ if | curt + dt >= maxt ->+ if maxt == 1 && f 1 >= 0+ then 1+ --bisect+ else+ let dT = (maxt - curt)/2+ in if dT < eps then curt +dT+ else iter f eps curt maxt (curt+dT)+ | otherwise ->+ iter f eps curt maxt (curt+dt)+ where+ dt = f curt++-- | Find the x value of the cubic bezier. The bezier must be+-- monotonically increasing in the X coordinate.+findX :: CubicBezier Double -> Double -> Double -> Double+findX (CubicBezier (Point p0 _) (Point p1 _) (Point p2 _) (Point p3 _)) x =+ find0 (p0-x) (p1-x) (p2-x) (p3-x)++find0 :: Double -> Double -> Double -> Double -> Double -> Double+find0 a b c d eps =+ iter bezierZero eps 0 1 0.5+ where+ bezierZero t+ | bx == 0 = t+ | otherwise = -bx/bx'+ where + u = 1-t+ t2 = t*t+ t3 = t2*t+ da = 3*^(b-a)+ db = 3*^(c-b)+ dc = 3*^(d-c)+ bx = u*(u*(u*a + (3*t)*b) + (3*t2)*c) + t3*d+ bx' = u*(u*da + (2*t)*db) + t2*dc+ -- | 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)
Geom2D/CubicBezier/Intersection.hs view
@@ -1,14 +1,14 @@ {-# 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)+ (bezierIntersection, bezierLineIntersections, bezierFindRoot) where import Geom2D import Geom2D.CubicBezier.Basic import Math.BernsteinPoly import Data.Maybe+import Geom2D.CubicBezier.Numeric import qualified Data.Vector.Unboxed as V-import Debug.Trace -- find the convex hull by comparing the angles of the vectors with -- the cross product and backtracking if necessary.@@ -92,65 +92,65 @@ -> Double -> Double -> Double -> Double -> Bool -> [(Double, Double)] bezierClip p@(CubicBezier !p0 !p1 !p2 !p3) q@(CubicBezier !q0 !q1 !q2 !q3)- tmin tmax umin umax prevClip eps revCurves- -- 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 =- if | abs (dmax - dmin) < eps * vectorDistance p0 p3 ->- -- fat line is smaller than tolerance.- if revCurves- then [(umin, tmin), (umax, tmax)]- else [(tmin, umin), (umin, tmin)]- | new_tmax - new_tmin > umax - umin ->- -- split the longest segment- let (pl, pr) = splitBezier newP 0.5- half_t = new_tmin + (new_tmax - new_tmin) / 2- in bezierClip q pl umin umax new_tmin half_t- newClip eps (not revCurves) ++- bezierClip q pr umin umax half_t new_tmax- newClip eps (not revCurves)- | otherwise ->- let (ql, qr) = splitBezier q 0.5- half_t = umin + (umax - umin) / 2- 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)- -- iterate with the curves reversed.- | otherwise =- bezierClip q newP umin umax new_tmin- new_tmax newClip eps (not revCurves)-- where- q3' | q0 == q3 =- q0 ^+^ (rotate90L $* p3 ^-^ p0)- | otherwise = q3- d = lineDistance (Line q0 q3')- d1 = d q1- d2 = d q2- (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])- chop_interval = chopHull dmin dmax $- map d [p0, p1, p2, p3]- Just (chop_tmin, chop_tmax) = chop_interval- newP = bezierSubsegment p chop_tmin chop_tmax- newClip = chop_tmax - chop_tmin- new_tmin = tmax * chop_tmin + tmin * (1 - chop_tmin)- new_tmax = tmax * chop_tmax + tmin * (1 - chop_tmax)+ tmin tmax umin umax prevClip eps revCurves = either id id $ do+ q3' <- if | vectorDistance q0 q3 > max (vectorMag q0) (vectorMag q3) / (2**30) -> Right q3+ | vectorDistance q0 q1 > max (vectorMag q0) (vectorMag q1) / (2**30) -> Right q1+ | vectorDistance q0 q2 > max (vectorMag q0) (vectorMag q2) / (2**30) -> Right q2+ | otherwise -> Left $+ let t = closest p q0 eps+ newT = tmin * (1-t) + tmax * t+ umid = umin + (umax-umin)/2+ in if revCurves then [(umid, newT)] else [(newT, umid)]+ let d = lineDistance (Line q0 q3')+ d1 = d q1+ d2 = d q2+ (dmin, dmax) | d1*d2 > 0 = (3/4 * min 0 (min d1 d2),+ 3/4 * max 0 (max d1 d2))+ | otherwise = (4/9 * min 0 (min d1 d2),+ 4/9 * max 0 (max d1 d2))+ (chop_tmin, chop_tmax) <- maybe (Left []) Right $+ chopHull dmin dmax $+ map d [p0, p1, p2, p3]+ let newP = bezierSubsegment p chop_tmin chop_tmax+ newClip = chop_tmax - chop_tmin+ new_tmin = tmax * chop_tmin + tmin * (1 - chop_tmin)+ new_tmax = tmax * chop_tmax + tmin * (1 - chop_tmax)+ if | -- within tolerance + max (umax - umin) (new_tmax - new_tmin) < eps ->+ if revCurves+ then Right [ (umin + (umax-umin)/2,+ new_tmin + (new_tmax-new_tmin)/2) ]+ else Right [ (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 ->+ if | abs (dmax - dmin) < eps * vectorDistance p0 p3 ->+ -- fat line is smaller than tolerance.+ if revCurves+ then Right [(umin, new_tmin), (umax, new_tmax)]+ else Right [(new_tmin, umin), (new_tmax, umax)]+ | new_tmax - new_tmin > umax - umin ->+ -- split the longest segment+ let (pl, pr) = splitBezier newP 0.5+ half_t = new_tmin + (new_tmax - new_tmin) / 2+ in Right $ bezierClip q pl umin umax new_tmin half_t+ newClip eps (not revCurves) +++ bezierClip q pr umin umax half_t new_tmax+ newClip eps (not revCurves)+ | otherwise ->+ let (ql, qr) = splitBezier q 0.5+ half_t = umin + (umax - umin) / 2+ in Right $ 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)+ -- iterate with the curves swapped.+ | otherwise ->+ Right $ bezierClip q newP umin umax new_tmin+ new_tmax newClip eps (not revCurves) +maxEps :: Double maxEps = 1e-8 -- | Find the intersections between two Bezier curves, using the@@ -210,25 +210,10 @@ -- X-axis. Then we only need to test the Y-values for a zero. bezierLineIntersections :: CubicBezier Double -> Line Double -> Double -> [Double] bezierLineIntersections b (Line p q) eps =- bezierFindRoot (listToBernstein $ map pointY [p0, p1, p2, p3]) 0 1 $- bezierParamTolerance b eps- where (CubicBezier p0 p1 p2 p3) = + filter (\x -> x > 0 && x < 1) $+ cubicRoot (p3 - 3*p2 + 3*p1 - p0) (3*p2 - 6*p1 + 3*p0) (3*p1 - 3*p0) p0+ where (CubicBezier (Point p0 _) (Point p1 _) (Point p2 _) (Point p3 _)) = fromJust (inverse $ translate p $* rotateVec (q ^-^ p)) $* b---- | Find the closest value(s) on the bezier to the given point, within tolerance.-closest :: CubicBezier Double -> DPoint -> Double -> [Double]-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,3 +1,5 @@+{-# LANGUAGE DeriveFoldable #-}+{-# LANGUAGE DeriveTraversable #-} {-# LANGUAGE BangPatterns, DeriveFunctor #-} -- | This module implements an extension to paths as used in -- D.E.Knuth's /Metafont/. Metafont gives an alternate way@@ -65,20 +67,21 @@ import Geom2D.CubicBezier.Numeric data OpenMetaPath a = OpenMetaPath [(Point a, MetaJoin a)] (Point a)+ deriving (Functor, Traversable, Foldable) -- ^ A metapath with endpoints data ClosedMetaPath a = ClosedMetaPath [(Point a, MetaJoin a)] -- ^ A metapath with cycles. The last join -- joins the last point with the first.- deriving (Eq, Functor)+ deriving (Eq, Functor, Traversable, Foldable) data MetaJoin a = MetaJoin { metaTypeL :: MetaNodeType a -- ^ The nodetype going out of the -- previous point. The metafont default is -- @Open@.- , tensionL :: Tension+ , tensionL :: Tension a -- ^ The tension going out of the previous point. -- The metafont default is 1.- , tensionR :: Tension+ , tensionR :: Tension a -- ^ The tension going into the next point. -- The metafont default is 1. , metaTypeR :: MetaNodeType a@@ -87,7 +90,7 @@ } | Controls (Point a) (Point a) -- ^ Specify the control points explicitly.- deriving (Show, Eq, Functor)+ deriving (Show, Eq, Functor, Traversable, Foldable) data MetaNodeType a = Open -- ^ An open node has no direction specified. If@@ -95,7 +98,7 @@ -- same direction going into and going out from -- the node. If it is an endpoint or corner -- point, it will have curl of 1.- | Curl {curlgamma :: Double}+ | Curl {curlgamma :: a} -- ^ The node becomes and endpoint or a corner -- point. The curl specifies how much the segment -- `curves`. A curl of `gamma` means that the@@ -103,28 +106,28 @@ -- following node. | Direction {nodedir :: Point a} -- ^ The node has a given direction.- deriving (Eq, Show, Functor)+ deriving (Eq, Show, Functor, Foldable, Traversable) -data Tension = Tension {tensionValue :: Double}+data Tension a = Tension {tensionValue :: a} -- ^ The tension value specifies how /tense/ the curve is. -- A higher value means the curve approaches a line -- segment, while a lower value means the curve is more -- round. Metafont doesn't allow values below 3/4.- | TensionAtLeast {tensionValue :: Double}+ | TensionAtLeast {tensionValue :: a} -- ^ Like @Tension@, but keep the segment inside the -- bounding triangle defined by the control points, if -- there is one.- deriving (Eq, Show)+ deriving (Eq, Show, Functor, Foldable, Traversable) -instance Show a => Show (ClosedMetaPath a) where+instance (Show a, Real a) => Show (ClosedMetaPath a) where show (ClosedMetaPath nodes) = showPath nodes ++ "cycle" -instance Show a => Show (OpenMetaPath a) where+instance (Show a, Real a) => Show (OpenMetaPath a) where show (OpenMetaPath nodes lastpoint) = showPath nodes ++ showPoint lastpoint -showPath :: Show a => [(Point a, MetaJoin a)] -> String+showPath :: (Show a, Real a) => [(Point a, MetaJoin a)] -> String showPath = concatMap showNodes where showNodes (p, Controls u v) =@@ -137,10 +140,10 @@ | t1 == t2 = printf "tension %s.." (showTension t1) | otherwise = printf "tension %s and %s.." (showTension t1) (showTension t2)- showTension (TensionAtLeast t) = printf "atleast %.3f" t :: String- showTension (Tension t) = printf "%.3f" t :: String+ showTension (TensionAtLeast t) = printf "atleast %.3f" (realToFrac t :: Double) :: String+ showTension (Tension t) = printf "%.3f" (realToFrac t :: Double) :: String typename Open = ""- typename (Curl g) = printf "{curl %.3f}" g :: String+ typename (Curl g) = printf "{curl %.3f}" (realToFrac g :: Double) :: String typename (Direction dir) = printf "{%s}" (showPoint dir) :: String showPoint :: Show a => Point a -> String@@ -159,8 +162,6 @@ path = joinSegments $ map unmetaSubSegment subsegs in OpenPath path endpoint -- unmetaClosed :: ClosedMetaPath Double -> ClosedPath Double unmetaClosed (ClosedMetaPath nodes) = case spanList bothOpen (removeEmptyDirs nodes) of@@ -359,7 +360,7 @@ zipNext l = zip l (tail $ cycle l) -- find the equations for a cycle containing only open points-eqsCycle :: [Tension] -> [DPoint] -> [Tension]+eqsCycle :: [Tension Double] -> [DPoint] -> [Tension Double] -> [Double] -> [(Double, Double, Double, Double)] eqsCycle tensionsA points tensionsB turnAngles = zipWith4 eqTension@@ -428,7 +429,8 @@ d = tensionA * tensionA / (tensionB * dist) -- the equation for a starting curl-eqCurl0 :: Double -> Double -> Double -> Double -> (Double, Double, Double, Double)+eqCurl0 :: Double -> Double -> Double -> Double+ -> (Double, Double, Double, Double) eqCurl0 gamma tensionA tensionB psi = (0, c, d, r) where c = chi/tensionA + 3 - 1/tensionB@@ -437,7 +439,8 @@ r = -d*psi -- the equation for an ending curl-eqCurlN :: Double -> Double -> Double -> (Double, Double, Double, Double)+eqCurlN :: Double -> Double -> Double+ -> (Double, Double, Double, Double) eqCurlN gamma tensionA tensionB = (a, b, 0, 0) where a = (3 - 1/tensionB)*chi + 1/tensionA@@ -445,7 +448,8 @@ chi = gamma*tensionA*tensionA / (tensionB*tensionB) -- getting the control points-unmetaJoin :: DPoint -> DPoint -> Double -> Double -> Tension -> Tension -> PathJoin Double+unmetaJoin :: DPoint -> DPoint -> Double -> Double -> Tension Double+ -> Tension Double -> PathJoin Double unmetaJoin !z0 !z1 !theta !phi !alpha !beta | abs phi < 1e-4 && abs theta < 1e-4 = JoinLine | otherwise = JoinCurve u v@@ -479,7 +483,7 @@ -- another magic formula by John Hobby. velocity :: Double -> Double -> Double- -> Double -> Tension -> Double+ -> Double -> Tension Double -> Double velocity st sf ct cf t = min 4 $ (2 + sqrt2 * (st - sf/16)*(sf - st/16)*(ct - cf)) /
Geom2D/CubicBezier/Numeric.hs view
@@ -3,7 +3,7 @@ -- contains functions that aren't used anymore, but might be useful on -- its own. module Geom2D.CubicBezier.Numeric- (quadraticRoot, solveLinear2x2, + (quadraticRoot, cubicRoot, cubicRootNorm, solveLinear2x2, goldSearch, makeSparse, SparseMatrix, sparseMulT, sparseMul, addMatrix, addVec, lsqMatrix, lsqSolveDist,@@ -38,6 +38,29 @@ result | d < 0 = [] | d == 0 = [x1] | otherwise = [x1, x2]++cubicRoot :: Double -> Double -> Double -> Double -> [Double]+cubicRoot a b c d+ | a == 0 = quadraticRoot b c d+ | otherwise = cubicRootNorm (b/a) (c/a) (d/a)++-- | Find real roots from the normalized cubic equation.+cubicRootNorm :: Double -> Double -> Double -> [Double]+cubicRootNorm a b c+ | r2 < q3 = [m2sqrtQ * cos(t/3) - a/3,+ m2sqrtQ * cos((t+2*pi)/3) - a/3,+ m2sqrtQ * cos((t - 2*pi)/3) - a/3]+ | otherwise = [d+e-a/3]+ where q = (a*a - 3*b) / 9+ q3 = q*q*q+ r = (2*a*a*a - 9*a*b + 27*c) / 54+ t = acos(r/sqrt q3)+ m2sqrtQ = -2*sqrt q+ r2 = r*r+ d = - sign(r)*((abs r + sqrt(r2-q3))**1/3)+ e | d == 0 = 0+ | otherwise = q/d+ -- | @solveLinear2x2 a b c d e f@ solves the linear equation with two variables (x and y) and two systems: --
Geom2D/CubicBezier/Overlap.lhs view
@@ -1,4 +1,4 @@-> {-# LANGUAGE MultiWayIf, PatternGuards, TemplateHaskell, BangPatterns #-}+> {-# LANGUAGE MultiWayIf, PatternGuards, TemplateHaskell, BangPatterns, CPP #-} Removing overlap from bezier paths in haskell =============================================@@ -30,7 +30,7 @@ * modifying state: - `over field fun struct`: `struct.field = fun (struct.field)` - `modify fun`: `state = fun state`- - `field %= fun`: `state.field = fun (state.field)`+ - `modifying field fun`: `state.field = fun (state.field)` Let's begin with declaring the module and library imports: @@ -38,24 +38,34 @@ > (boolPathOp, union, intersection, difference, > exclusion, FillRule (..)) > where-> import Prelude hiding (mapM)+> import Prelude > import Geom2D > import Geom2D.CubicBezier.Basic > import Geom2D.CubicBezier.Intersection+> import Geom2D.CubicBezier.Numeric > import Math.BernsteinPoly-> import Data.Traversable (mapM)+> import Data.Foldable (traverse_) > import Data.Functor ((<$>))-> import Data.List (sortBy, sort)-> import Control.Monad.State hiding (mapM)+> import Data.List (sortBy, sort, intercalate, intersperse)+> import Control.Monad.State.Strict > import Lens.Micro > import Lens.Micro.TH > import Lens.Micro.Mtl > import qualified Data.Map.Strict as M-> import qualified Data.Set as S+> import qualified Data.Set as S +> import Text.Printf+> import Data.Ratio+> import Data.Maybe (isJust, isNothing, mapMaybe) +#ifdef DEBUG+> import System.IO.Unsafe (unsafePerformIO)+> import System.IO+> import Debug.Trace+#endif+ 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+Since that could be true only of a part of the curve, I 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@@ -66,12 +76,12 @@ turnratio changes with each curve. Checking each pair of curves for intersections would work, but is-rather inefficient. We only need to check for overlap when two curves+rather inefficient. I only need to check for overlap when two curves are adjacent. Fortunately there exist a good method from *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+possible), and to update the input dynamically. I keep track of each curve that intersects the sweepline by using a balanced tree of curves. When adding a new curve, it's only necessary to check for intersections with the curve above and below. Since searching on the@@ -94,7 +104,9 @@ point. > newtype PointEvent = PointEvent DPoint-> deriving Show+>+> instance Show PointEvent where+> show (PointEvent (Point px py)) = printf "(%.5g, %.5g)" px py When the x-coordinates are equal, use the y-coordinate to determine the order.@@ -107,7 +119,7 @@ > compare (PointEvent (Point x1 y1)) (PointEvent (Point x2 y2)) = > compare (x1, y2) (x2, y1) -All curves are kept left to right, so we need to remember the+All curves are kept left to right, so I need to remember the curve direction for the output: The curves intersecting the sweepline are kept in another balanced@@ -122,6 +134,9 @@ allows for more flexibility. Using this, it is possible to generalize this algorithm to boolean operations! +The curveRank parameter is used to memoize the order in the Ystruct.+This will avoid costly comparisons for tree rebalancing etc...+ The FillRule datatype is used for the exported API: > data FillRule = EvenOdd | NonZero@@ -129,16 +144,22 @@ > data Curve = Curve { > _bezier :: !(CubicBezier Double), > _turnRatio :: !(Int, Int),-> _changeTurn :: !((Int, Int) -> (Int, Int))}+> _changeTurn :: !((Int, Int) -> (Int, Int)),+> _curveRank :: Maybe (Ratio Integer)} > > trOne :: (Int, Int) > trOne = (0,0)+>+> between :: Ratio Integer -> Ratio Integer -> Ratio Integer+> between a b = a +(b-a)/2 > > makeLenses ''Curve > > instance Show Curve where-> show (Curve b a _) =-> "Curve " ++ show b ++ " " ++ show a+> show (Curve (CubicBezier (Point p0x p0y) (Point p1x p1y)+> (Point p2x p2y) (Point p3x p3y)) (t1, t2) _ o) =+> printf "Curve (%.5g, %.5g) (%.5g, %.5g) (%.5g, %.5g) (%.5g, %.5g) (%i,%i) %s" +> p0x p0y p1x p1y p2x p2y p3x p3y t1 t2 (show o) > > type YStruct = S.Set Curve @@ -149,111 +170,208 @@ Y-structure, where the left set are the elements less than the pointEvent (above), and the right set the elements greater (below): - > data SweepState = SweepState { > _output :: !(M.Map PointEvent [CubicBezier Double]),-> _yStructLeft :: !YStruct,-> _yStructRight :: !YStruct,+> _yStruct :: !YStruct,+> _focusPoint :: DPoint, > _xStruct :: !XStruct}-> deriving Show > > makeLenses ''SweepState -Changing the focus point can be done efficiently in `O(log n)` by-mering and splitting again:+> singularC :: Point Double -> Curve+> singularC p = Curve (CubicBezier p p p p) trOne id Nothing+> -> changeFocus :: DPoint -> SweepState -> SweepState-> changeFocus p sweep =-> let (lStr, rStr) =-> S.split (Curve (CubicBezier p p p p) trOne id) $-> S.union (view yStructLeft sweep) (view yStructRight sweep)-> in set yStructLeft lStr $-> set yStructRight rStr-> sweep+This handy function will be optimized away when DEBUG is false. +> showCurve (CubicBezier p0 p1 p2 p3) =+> showPt p0 ++ showPt p1 ++ showPt p2 ++ showPt p3+>+> showPt :: DPoint -> String+> showPt (Point x y) = "(" ++ show x ++ "," ++ show y ++ ")"+>+#ifdef DEBUG+> type SweepStateM = StateT SweepState IO+>+> traceMessage :: String -> SweepStateM ()+> traceMessage = liftIO . hPutStrLn stderr+> +> assert :: Bool -> String -> SweepStateM ()+> assert p msg = unless p $ do+> liftIO $ hPutStrLn stderr $ "ASSERT " ++ msg+>+#else+> type SweepStateM = State SweepState+>+> traceMessage _ = return ()+>+> assert :: Bool -> String -> SweepStateM ()+> assert _ _ = return ()+#endif++> activate, deactivate :: [Curve] -> SweepStateM ()+> activate cs = +> traverse_ (traceMessage . ("ACTIVATE " ++) . showCurve . view bezier) cs+> +> deactivate cs = +> traverse_ (traceMessage . ("DEACTIVATE " ++) . showCurve . view bezier) cs++ This handy helper function will pass the first curve above to the given function, and if it doesn't return `Nothing`, remove it from the state. It does nothing when there is no curve above. -> withAbove :: (Curve -> Maybe a) -> State SweepState (Maybe a)+> withAbove :: (Curve -> Maybe a) -> SweepStateM (Maybe a) > withAbove f = do-> lStr <- use yStructLeft-> if S.null lStr+> p <- use focusPoint+> yStr <- use yStruct+> let i = yStructIndex (singularC p) yStr+> if i < 0 > then return Nothing-> else let (c, lStr') = S.deleteFindMax lStr+> else let c = S.elemAt i yStr > in case f c of-> Nothing ->-> return Nothing-> Just x -> do-> yStructLeft .= lStr'-> return $ Just x+> Nothing ->+> return Nothing+> Just x -> do+> yStructDel i+> return $ Just x The same with the curve below. -> withBelow :: (Curve -> Maybe a) -> State SweepState (Maybe a)+> withBelow :: (Curve -> Maybe a) -> SweepStateM (Maybe a) > withBelow f = do-> rStr <- use yStructRight-> if S.null rStr+> p <- use focusPoint+> yStr <- use yStruct+> let i = yStructIndex (singularC p) yStr+> s = S.size yStr+> if i >= s-1 > then return Nothing-> else let (c, rStr') = S.deleteFindMin rStr+> else let c = S.elemAt (i+1) yStr > in case f c of-> Nothing ->-> return Nothing-> Just x -> do-> yStructRight .= rStr'-> return $ Just x+> Nothing ->+> return Nothing+> Just x -> do+> yStructDel (i+1)+> return $ Just x `splitYStruct` changes the focus and returns and removes any curves which end in the current pointEvent: -> splitYStruct :: DPoint -> State SweepState [Curve]+> splitYStruct :: DPoint -> SweepStateM [Curve] > splitYStruct p = do-> modify $ changeFocus p-> let go = do-> mbC <- withAbove $ \c ->-> -- remove and return c if it ends in point p-> -> guard (cubicC3 (_bezier c) == p) >> Just c-> case mbC of-> Just c ->-> (c:) <$> go-> Nothing -> return []-> go-+> yStr <- use yStruct+> focusPoint .= p+> traceMessage $ "CHANGEFOCUS " ++ showPt p+> traceMessage "MSG Remove curves ending at pointevent from Y structure" +> let lStr = fst $ S.split (singularC p) yStr+> rightCurves = takeWhile (\c -> cubicC3 (_bezier c) == p) $+> S.toDescList lStr+> nR = length rightCurves+> i = S.size lStr - nR+> replicateM_ nR (yStructDel i)+> return rightCurves+> === Some functions on the Sweep state: Adding and removing curves from the X structure. -> insertX :: PointEvent -> [Curve] -> SweepState -> SweepState+> insertX :: PointEvent -> [Curve] -> SweepStateM () > insertX p c =-> over xStruct $ M.insertWith (++) p c+> modify $ over xStruct $ M.insertWith (++) p c >-> xStructAdd :: Curve -> SweepState -> SweepState-> xStructAdd c =+> xStructAdd :: Curve -> SweepStateM ()+> xStructAdd c = do+> traceMessage $ "XSTRUCTADD " ++ showCurve (view bezier c) > insertX (PointEvent $ cubicC0 $-> view bezier c) [c]+> view bezier c) [c] >-> xStructRemove :: State SweepState (PointEvent, [Curve])-> xStructRemove = zoom xStruct $ state M.deleteFindMin+> xStructRemove :: SweepStateM (PointEvent, [Curve])+> xStructRemove = do+> str <- use xStruct+> return str+> (p, c) <- zoom xStruct $ state M.deleteFindMin+> traverse_ (traceMessage . ("XSTRUCTREM " ++) .+> showCurve . view bezier) c+> str <- use xStruct+> return str+> return (p, c)+>+> yStructIndex :: Curve -> YStruct -> Int+> yStructIndex c str = S.size (fst $ S.split c str) - 1+>+> yStructDel :: Int -> SweepStateM ()+> yStructDel i = do+> str <- use yStruct+> traceMessage $ "YSTRUCTREM " ++ showCurve (view bezier $ S.elemAt i str)+> yStruct .= S.deleteAt i str+> +Insert the curve into the Y structure. First lookup the position of+the curve, then calculate the rank of the curve, using the surrounding+elements. Insert the curve using the rank. This will avoid repeating+expensive operations.++> +> yStructAdd :: Curve -> SweepStateM ()+> yStructAdd c = do+> str <- use yStruct+> traceMessage $ "YSTRUCTADD " ++ showCurve (view bezier c)+> traceMessage $ "YSTRUCT: " +++> (concat $ intersperse "\n " $ map (showCurve . view bezier) $ S.toList str)+> assert (isNothing (_curveRank c))+> "CURVE ALREADY HAS A RANK IN THE YSTRUCT" +> assert (yStructConsistent c str)+> ("Y STRUCT NOT CONSISTENT WITH CURVE:" +++> show (map (compare c) (S.toAscList str)) )+> let i = yStructIndex c str+> s = S.size str+> newC +> | s == 0 = set curveRank (Just 0) c+> | i >= s-1 = set curveRank ((+1) <$> _curveRank (S.elemAt (s-1) str)) c+> | i < 0 = set curveRank (subtract 1 <$> _curveRank (S.elemAt 0 str)) c+> | otherwise = set curveRank (liftM2 between+> (_curveRank $ S.elemAt i str)+> (_curveRank $ S.elemAt (i+1) str)) c+> assert (not $ S.member c str)+> ("CURVE ALREADY IN YSTRUCT: " ++ showCurve (view bezier c))+> assert (S.size str < S.size (S.insert newC str)) $+> "CURVE NOT ADDED TO YSTRUCT" ++ show newC +> yStruct .= S.insert newC str+> +> yStructConsistent :: Curve -> YStruct -> Bool+> yStructConsistent c str =+> all (> c) $ dropWhile (< c) $+> S.toAscList str+>+> yStructOverlap :: Curve -> YStruct -> [String]+> yStructOverlap c str =+> mapMaybe checkOverlap $ S.toAscList str+> where checkOverlap c2 =+> case splitMaybe c c2 1e-5 of+> (Nothing, Nothing) -> Nothing+> _ -> Just (show c2 ++ "\n")+ To compare curves vertically, take the the curve which starts the-rightmost, and see if it falls below or above the curve. If the first control points are coincident, test the-last control points instead. The curves in the Y-structure shouldn't-intersect (except in the endpoints), so these cases don't have to be-handled. To lookup a single point, I use a singular bezier curve.+rightmost, and see if it falls below or above the curve. If the first+control points are coincident, test the last control points instead,+or the midpoint. This works because if the first point is coincident+the curves shouldn't intersect except in the endpoints (see #splitAndOrder).+To lookup a single point, I use a singular bezier curve. > instance Eq Curve where-> Curve c1 t1 ct1 == Curve c2 t2 ct2 =-> c1 == c2 && t1 == t2 && ct1 (ct2 t1) == t1+> Curve _ _ _ (Just o1) == Curve _ _ _ (Just o2) = o1 == o2+> Curve c1 t1 ct1 _ == Curve c2 t2 ct2 _ =+> c1 == c2 && t1 == t2 && ct1 t1 == ct2 t2 > > instance Ord Curve where-> compare (Curve c1@(CubicBezier p0 p1 p2 p3) tr1 _)-> (Curve c2@(CubicBezier q0 q1 q2 q3) tr2 _)+> compare (Curve _ _ _ (Just o1)) (Curve _ _ _ (Just o2)) = compare o1 o2+> compare (Curve c1@(CubicBezier p0 p1 p2 p3) tr1 _ _)+> (Curve c2@(CubicBezier q0 q1 q2 q3) tr2 _ _) > | p0 == q0 = if > | p3 == q3 -> > -- compare the midpoint-> case (compVert (evalBezier c1 0.5) c2) of+> case compVert (evalBezier c1 0.5) c2 of > LT -> LT > GT -> GT > EQ ->@@ -261,7 +379,7 @@ > compare (tr1, PointEvent p1, PointEvent p2) > (tr2, PointEvent q1, PointEvent q2) > | pointX p3 < pointX q3 ->-> case (compVert p3 c2) of+> case compVert p3 c2 of > LT -> LT > EQ -> LT > GT -> GT@@ -276,14 +394,14 @@ > EQ -> LT > GT -> LT > | otherwise =-> case (compVert p0 c2) of+> case compVert p0 c2 of > LT -> LT > EQ -> GT > GT -> GT Compare a point with a curve. See if it falls below or above the hull first. Otherwise find the point on the curve with the same-X-coordinate by solving a cubic equation.+X-coordinate by iterating. > compVert :: DPoint -> CubicBezier Double -> Ordering > compVert p c@@ -305,9 +423,7 @@ > compare (pointY p0) y1 > | otherwise = compare y2 y1 > where-> t = findX x1 c1 $-> maximum (map maxp [p0, p1, p2, p3])*1e-12-> maxp (Point x y) = max (abs x) (abs y)+> t = findX c1 x1 (maximum (map (abs.pointX) [p0, p1, p2, p3])*1e-14) > y2 = pointY $ evalBezier c1 t === Comparing against the hull {#hull}@@ -363,18 +479,18 @@ > concatMap (splitVert tol) > where toCurve c@(CubicBezier p0 _ _ p3) = > case compare (pointX p0) (pointX p3) of-> LT -> [(PointEvent p0, [Curve c trOne chTr])]-> GT -> [(PointEvent p3, [Curve (reorient c) trOne chTrBack]),+> LT -> [(PointEvent p0, [Curve c trOne chTr Nothing])]+> GT -> [(PointEvent p3, [Curve (reorient c) trOne chTrBack Nothing]), > (PointEvent p0, [])] > -- vertical curve > EQ | pointY p0 > pointY p3 ->-> [(PointEvent p0, [Curve c trOne chTr])]+> [(PointEvent p0, [Curve c trOne chTr Nothing])] > | otherwise ->-> [(PointEvent p3, [Curve (reorient c) trOne chTrBack]),+> [(PointEvent p3, [Curve (reorient c) trOne chTrBack Nothing]), > (PointEvent p0, [])] > > splitVert :: Double -> CubicBezier Double -> [CubicBezier Double]-> splitVert tol curve@(CubicBezier c0 c1 c2 c3) =+> splitVert tol curve@(CubicBezier c0 c1 c2 c3) = > uncurry splitBezierN $ > adjustLast $ > adjustFirst (curve, vert)@@ -400,7 +516,7 @@ main loop --------- -For the main loop, we remove the leftmost point from the+For the main loop, I remove the leftmost point from the X-structure, and do the following steps: 1. Split any curves which come near the current pointEvent.@@ -417,37 +533,61 @@ 5. Loop until the X-structure is empty -> loopEvents :: ((Int, Int) -> Bool) -> Double -> SweepState -> SweepState-> loopEvents isInside tol sweep -> | M.null $ view xStruct sweep = sweep-> | otherwise =-> loopEvents isInside tol $!-> flip execState sweep $ do-> -- remove leftmost point from X structure+> loopEvents :: ((Int, Int) -> Bool) -> Double -> SweepStateM ()+> loopEvents isInside tol = do+> xStr <- use xStruct+> unless (M.null xStr) $ do > (PointEvent p, curves) <- xStructRemove-> -- change focus, and remove curves ending at current-> -- pointevent from Y structure+> activate curves+> > ending <- splitYStruct p+> activate ending+> > -- split near curves+> traceMessage "MSG Split curves near the focuspoint." > (ending2, rightSubCurves) <- splitNearPoints p tol+>+> activate ending2+> activate rightSubCurves+> traceMessage "MSG Output curves"+> > -- output curves to the left of the sweepline.-> modify $ filterOutput (ending ++ ending2) isInside +> deactivate (ending ++ ending2)+> filterOutput (ending ++ ending2) isInside > let allCurves = rightSubCurves ++ curves > if null allCurves+> > -- split surrounding curves-> then splitSurround tol+> then do+> traceMessage "MSG Split curves around pointevent."+> splitSurround tol > else do+> > -- sort curves+> traceMessage "MSG Sort curves."+> > sorted <- splitAndOrder tol allCurves+> deactivate allCurves+> activate sorted+> > -- split curve above+> traceMessage "MSG Split curve above sorted curves."+> deactivate sorted > curves2 <- splitAbove sorted tol+> activate curves2+> > -- add curves to Y structure+> traceMessage "MSG Add curves to Y structure."+> deactivate curves2 > addMidCurves curves2 tol+>+> loopEvents isInside tol + Send curves to output --------------------- -> outputPaths :: (M.Map PointEvent [CubicBezier Double]) -> [ClosedPath Double]+> outputPaths :: M.Map PointEvent [CubicBezier Double] -> [ClosedPath Double] > outputPaths m > | M.null m = [] > | otherwise = outputNext m@@ -462,7 +602,7 @@ > outputNext !m > | M.null m = [] > | otherwise = -> let ((PointEvent p0, (c0:cs)), m0) =+> let ((PointEvent p0, c0:cs), m0) = > M.deleteFindMin m > m0' | null cs = m0 > | otherwise = M.insert (PointEvent p0) cs m0@@ -488,20 +628,22 @@ nonzero-rule, this would be `(> 0)`. This inserts the curve into the output map. -> filterOutput :: [Curve] -> ((Int, Int) -> Bool) -> SweepState -> SweepState-> filterOutput curves isInside sweep =-> foldl (flip $ outputCurve isInside) sweep curves+> filterOutput :: [Curve] -> ((Int, Int) -> Bool) -> SweepStateM ()+> filterOutput curves isInside =+> mapM_ (outputCurve isInside) curves >-> outputCurve :: ((Int, Int) -> Bool) -> Curve -> SweepState -> SweepState-> outputCurve isInside (Curve c tr op)+> outputCurve :: ((Int, Int) -> Bool) -> Curve -> SweepStateM ()+> outputCurve isInside (Curve c tr op _) > | isInside (op tr) /= isInside tr = > let c' | isInside tr = reorient c > | otherwise = c-> in over output (M.insertWith (++) (PointEvent $ cubicC0 c') [c'])-> | otherwise = id+> in do traceMessage $ "OUTPUT " ++ showCurve c+> modifying output (M.insertWith (++) (PointEvent $ cubicC0 c') [c'])+> | otherwise =+> traceMessage $ "DISCARD " ++ showCurve c -Test for intersections and split:----------------------------------+Test for intersections and split: (#splitAndOrder)+-------------------------------------------------- Since the curves going out of the current pointEvent in the X-structure are unordered, they need to be ordered first. First they are ordered by@@ -512,31 +654,35 @@ To do this, I implemented a monadic insertion sort. First the curves are split in the statemonad, then they are compared. -> splitAndOrder :: Double -> [Curve] -> State SweepState [Curve]+> splitAndOrder :: Double -> [Curve] -> SweepStateM [Curve] > splitAndOrder tol curves = > sortSplit tol $ > sortBy compDeriv curves > > compDeriv :: Curve -> Curve -> Ordering-> compDeriv (Curve (CubicBezier p0 p1 _ _) _ _)-> (Curve (CubicBezier q0 q1 _ _) _ _) =-> compare (vectorCross (p1^-^p0) (q1^-^ q0)) 0+> compDeriv (Curve (CubicBezier p0 p1 _ _) _ _ _)+> (Curve (CubicBezier q0 q1 _ _) _ _ _) =+> compare (slope (q1^-^ q0)) (slope (p1^-^p0)) +>+> slope (Point 0 0) = 0+> slope (Point x y) = y/x + Insertion sort, by splitting and comparing. This should be efficient enough, since ordering by derivative should mostly order the curves. -> sortSplit :: Double -> [Curve] -> State SweepState [Curve]+> sortSplit :: Double -> [Curve] -> SweepStateM [Curve] > sortSplit _ [] = return [] > sortSplit tol (x:xs) = > insertM x tol =<< > sortSplit tol xs >-> insertM :: Curve -> Double -> [Curve] -> State SweepState [Curve]+> insertM :: Curve -> Double -> [Curve] -> SweepStateM [Curve] > insertM x _ [] = return [x] > insertM x tol (y:ys) = > case curveOverlap x y tol of > Just (c1, c2) -> do-> mapM (modify . xStructAdd) c2+> traverse_ xStructAdd c2 > insertM c1 tol ys > Nothing -> do > (x', y') <- splitM x y tol@@ -544,17 +690,18 @@ > then return (x':y':ys) > else (y':) <$> insertM x' tol ys >-> splitM :: Curve -> Curve -> Double -> State SweepState (Curve, Curve)+> splitM :: Curve -> Curve -> Double -> SweepStateM (Curve, Curve) > splitM x y tol = > case splitMaybe x y tol of > (Just (a, b), Just (c, d)) -> do-> modify $ insertX (PointEvent $ cubicC0 $ view bezier b) [b, d]+> xStructAdd b+> xStructAdd d > return (a, c) > (Nothing, Just (c, d)) -> do-> modify $ insertX (PointEvent $ cubicC0 $ view bezier d) [d]+> xStructAdd d > return (x, c) > (Just (a, b), Nothing) -> do-> modify $ insertX (PointEvent $ cubicC0 $ view bezier b) [b]+> xStructAdd b > return (a, y) > (Nothing, Nothing) -> > return (x, y)@@ -565,49 +712,50 @@ to create singular curves. > updateTurnRatio :: Curve -> Curve -> Curve-> updateTurnRatio (Curve _ tr chTr) =+> updateTurnRatio (Curve _ tr chTr _) = > set turnRatio (chTr tr) > > propagateTurnRatio :: Curve -> [Curve] -> [Curve] > propagateTurnRatio cAbove l = > tail $ scanl updateTurnRatio cAbove l >-> splitAbove :: [Curve] -> Double -> State SweepState [Curve]+> splitAbove :: [Curve] -> Double -> SweepStateM [Curve] > splitAbove [] _ = return [] > splitAbove (c:cs) tol = do-> lStr <- use yStructLeft-> if S.null lStr+> yStr <- use yStruct+> p <- use focusPoint+> let i = yStructIndex (singularC p) yStr+> if i < 0 > then let c' = set turnRatio trOne c-> in return $ c':propagateTurnRatio c' cs-> else do-> let (cAbove, lStr') = S.deleteFindMax lStr-> case splitMaybe c cAbove tol of-> (Nothing, Nothing) ->-> return $ propagateTurnRatio cAbove $ c:cs-> (Just (c1, c2), Nothing) ->-> if cubicC3 (_bezier c1) == cubicC0 (_bezier cAbove)-> then do-> modify $ xStructAdd cAbove . xStructAdd c2-> yStructLeft .= lStr'-> return $ propagateTurnRatio cAbove $ c1:cs-> else do-> modify $ xStructAdd c2-> return $ propagateTurnRatio cAbove $ c1:cs-> (Nothing, Just (c3, c4)) ->-> if cubicC3 (_bezier c3) == cubicC0 (_bezier c)-> then error "curve intersecting pointevent"-> else do-> modify $ xStructAdd c4-> yStructLeft .= S.insert c3 lStr'-> return $ propagateTurnRatio cAbove $ c:cs-> (Just (c1, c2), Just (c3, c4)) -> do-> modify $ xStructAdd c2 . xStructAdd c4-> yStructLeft .= S.insert c3 lStr'-> return $ propagateTurnRatio cAbove $ c1:cs+> in return $ c':propagateTurnRatio c' cs+> else +> let cAbove = S.elemAt i yStr+> in case splitMaybe c cAbove tol of+> (Nothing, Nothing) ->+> return $ propagateTurnRatio cAbove $ c:cs+> (Just (c1, c2), Nothing)+> | cubicC3 (_bezier c1) == cubicC0 (_bezier cAbove)+> -> do+> xStructAdd cAbove; xStructAdd c2+> yStructDel i+> return $ propagateTurnRatio cAbove $ c1:cs+> | otherwise -> do+> xStructAdd c2+> return $ propagateTurnRatio cAbove $ c1:cs+> (Nothing, Just (c3, c4)) -> do+> assert (cubicC3 (_bezier c3) /= cubicC0 (_bezier c)) $+> "curve intersecting pointevent: cAbove " ++ show cAbove+> xStructAdd c4+> yStructDel i; yStructAdd c3+> return $ propagateTurnRatio cAbove $ c:cs+> (Just (c1, c2), Just (c3, c4)) -> do+> xStructAdd c2; xStructAdd c4+> yStructDel i; yStructAdd c3+> return $ propagateTurnRatio cAbove $ c1:cs -Split curves near the point. Return the curves starting from this point, and the index of the last split point+Split curves near the point. Return the curves starting from this point. -> splitNearPoints :: DPoint -> Double -> State SweepState ([Curve], [Curve])+> splitNearPoints :: DPoint -> Double -> SweepStateM ([Curve], [Curve]) > splitNearPoints p tol = do > curves1 <- splitNearDir withAbove p tol > curves2 <- splitNearDir withBelow p tol@@ -615,9 +763,9 @@ > map snd curves1 ++ map snd curves2) > > splitNearDir :: ((Curve -> Maybe (Curve, Double))-> -> State SweepState (Maybe (Curve, Double)))+> -> SweepStateM (Maybe (Curve, Double))) > -> DPoint -> Double-> -> State SweepState [(Curve, Curve)]+> -> SweepStateM [(Curve, Curve)] > splitNearDir dir p tol = do > mbSplit <- dir $ \curve -> > (,) curve <$>@@ -631,75 +779,75 @@ > snapRound tol <$> c1 > c2' = adjust curve $ adjustC0 p $ > snapRound tol <$> c2+> traceMessage $ "MSG Splitting curve " ++ showCurve (view bezier curve) > ((c1', c2'):) <$> splitNearDir dir p tol Add the sorted curves starting at point to the Y-structure, and test last curve with curve below. -> addMidCurves :: [Curve] -> Double -> State SweepState ()+> addMidCurves :: [Curve] -> Double -> SweepStateM () > addMidCurves [] _ = return () > addMidCurves [c] tol = > splitBelow c tol > addMidCurves (c:cs) tol = do-> yStructLeft %= S.insert c > addMidCurves cs tol+> yStructAdd c > -> splitBelow :: Curve -> Double -> State SweepState ()+> splitBelow :: Curve -> Double -> SweepStateM () > splitBelow c tol = do-> rStr <- use yStructRight-> let (cBelow, rStr') = S.deleteFindMin rStr-> if S.null rStr-> then yStructLeft %= S.insert c+> yStr <- use yStruct+> p <- use focusPoint+> let i = yStructIndex (singularC p) yStr+> if i >= S.size yStr-1+> then yStructAdd c > else-> case splitMaybe c cBelow tol of-> (Nothing, Nothing) ->-> yStructLeft %= S.insert c-> (Nothing, Just (c3, c4)) ->-> if cubicC3 (_bezier c3) == cubicC0 (_bezier c)-> then error "internal error: splitBelow: curve starting in future"-> else do-> modify $ xStructAdd c4-> yStructLeft %= S.insert c . S.insert c3-> yStructRight .= rStr'-> (Just (c1, c2), Nothing) ->-> if cubicC3 (_bezier c1) == cubicC0 (_bezier cBelow)-> then error "internal error: splitBelow: curve intersecting pointevent."-> else do-> modify $ xStructAdd c2-> yStructLeft %= S.insert c1-> (Just (c1, c2), Just (c3, c4)) -> do-> modify $ xStructAdd c2 . xStructAdd c4-> yStructLeft %= S.insert c1 . S.insert c3-> yStructRight .= rStr'+> let cBelow = S.elemAt (i+1) yStr+> in case splitMaybe c cBelow tol of+> (Nothing, Nothing) -> +> yStructAdd c+> (Nothing, Just (c3, c4)) -> do+> assert (cubicC3 (_bezier c3) /= cubicC0 (_bezier c)) $+> "splitBelow: curve starting in future: c3 == " +++> show c3 ++ " c == " ++ show c+> traceMessage "MSG Split lower curve only." +> xStructAdd c4+> yStructDel (i+1); yStructAdd c3; yStructAdd c+> (Just (c1, c2), Nothing) -> do+> traceMessage "MSG split curve above lower curve"+> assert (cubicC3 (_bezier c1) /= cubicC0 (_bezier cBelow))+> "SPLITBELOW: CURVE INTERSECTING POINTEVENT." +> xStructAdd c2+> yStructAdd c1+> (Just (c1, c2), Just (c3, c4)) -> do+> traceMessage "MSG split lower curve and curve above."+> xStructAdd c2; xStructAdd c4+> yStructDel (i+1); yStructAdd c3; yStructAdd c1 -If no curves start from the point, we have to check if the surrounding+If no curves start from the point, check if the surrounding curves overlap. -> splitSurround :: Double -> State SweepState ()+> splitSurround :: Double -> SweepStateM () > splitSurround tol = do-> lStr <- use yStructLeft-> rStr <- use yStructRight-> if S.null lStr || S.null rStr-> then return ()-> else do-> let (cAbove, lStr') = S.deleteFindMax lStr-> (cBelow, rStr') = S.deleteFindMin rStr-> case splitMaybe cAbove cBelow tol of+> p <- use focusPoint+> yStr <- use yStruct+> let i = yStructIndex (singularC p) yStr+> s = S.size yStr+> when (i >= 0 && i < s-1) $+> case splitMaybe (S.elemAt i yStr) (S.elemAt (i+1) yStr) tol of > (Just (c1, c2), Just (c3, c4)) -> do-> modify $ xStructAdd c2 .-> xStructAdd c4-> yStructLeft .= S.insert c1 lStr'-> yStructRight .= S.insert c3 rStr'+> xStructAdd c2; xStructAdd c4+> yStructDel i; yStructDel i+> yStructAdd c3; yStructAdd c1 > (Just (c1, c2), Nothing) -> do-> modify $ xStructAdd c2-> yStructLeft .= S.insert c1 lStr'+> xStructAdd c2+> yStructDel i; yStructAdd c1 > (Nothing, Just (c1, c2)) -> do-> modify $ xStructAdd c2-> yStructRight .= S.insert c1 rStr'+> xStructAdd c2+> yStructDel (i+1); yStructAdd c1 > (Nothing, Nothing) -> > return ()- + === Find curve intersections Test if both curves intersect. Split one or both of the curves when@@ -721,12 +869,12 @@ > b2 = view bezier c2 > > adjustSplit :: Curve -> (CubicBezier Double, CubicBezier Double) -> (Curve, Curve)-> adjustSplit curve (b1, b2) =-> (set bezier b1 curve,-> set bezier b2 curve)+> adjustSplit curve (b1, b2) =+> (set curveRank Nothing $ set bezier b1 curve,+> set curveRank Nothing $ set bezier b2 curve) > > adjust :: Curve -> CubicBezier Double -> Curve-> adjust curve curve2 = set bezier curve2 curve+> adjust curve curve2 = set curveRank Nothing $ set bezier curve2 curve > > snapRoundBezier :: Double -> CubicBezier Double -> CubicBezier Double > snapRoundBezier tol = fmap (snapRound tol)@@ -745,40 +893,30 @@ > nextIntersection b1@(CubicBezier p0 _ _ p3) b2@(CubicBezier q0 _ _ q3) tol ((t1, t2): ts) > | atStart1 && atStart2 = > nextIntersection b1 b2 tol ts-> | atStart1 =-> (Nothing,-> Just (adjustC3 p0 $ snapRoundBezier tol b2l,-> adjustC0 p0 $ snapRoundBezier tol b2r))-> | atStart2 =-> (Just (adjustC3 q0 $ snapRoundBezier tol b1l,-> adjustC0 q0 $ snapRoundBezier tol b1r),-> Nothing)-> | atEnd1 && atEnd2 = (Nothing, Nothing)-> | atEnd1 =-> (Nothing,-> Just (adjustC3 p3 $ snapRoundBezier tol b2l,-> adjustC0 p3 $ snapRoundBezier tol b2r))-> | atEnd2 =-> (Just (adjustC3 q3 $ snapRoundBezier tol b1l,-> adjustC0 q3 $ snapRoundBezier tol b1r),-> Nothing) > | bezierEqual b1l b2l tol = > nextIntersection b1 b2 tol ts > | otherwise =-> let pMid = snapRound tol <$> cubicC3 b1l-> in (Just (snapRoundBezier tol b1l,-> snapRoundBezier tol b1r),-> Just (adjustC3 pMid $ snapRoundBezier tol b2l,-> adjustC0 pMid $ snapRoundBezier tol b2r))-> where-> x1 = evalBezier b1 t1-> x2 = evalBezier b2 t2-> atStart1 = vectorDistance (cubicC0 b1) x1 < tol-> atStart2 = vectorDistance (cubicC0 b2) x2 < tol-> atEnd1 = vectorDistance (cubicC3 b1) x1 < tol-> atEnd2 = vectorDistance (cubicC3 b2) x2 < tol-> (b1l, b1r) = splitBezier b1 t1-> (b2l, b2r) = splitBezier b2 t2+> (bs1, bs2)+> where+> bs1 | atStart1 || atEnd1 = Nothing+> | otherwise = Just (adjustC3 pMid $ snapRoundBezier tol b1l,+> adjustC0 pMid $ snapRoundBezier tol b1r)+> bs2 | atStart2 || atEnd2 = Nothing+> | otherwise = Just (adjustC3 pMid $ snapRoundBezier tol b2l,+> adjustC0 pMid $ snapRoundBezier tol b2r)+> pMid | atStart1 = cubicC0 b1+> | atEnd1 = cubicC3 b1+> | atStart2 = cubicC0 b2+> | atEnd1 = cubicC3 b2+> | otherwise = snapRound tol <$> cubicC3 b1l+> x1 = evalBezier b1 t1+> x2 = evalBezier b2 t2+> atStart1 = vectorDistance (cubicC0 b1) x1 < tol+> atStart2 = vectorDistance (cubicC0 b2) x2 < tol+> atEnd1 = vectorDistance (cubicC3 b1) x1 < tol+> atEnd2 = vectorDistance (cubicC3 b2) x2 < tol+> (b1l, b1r) = splitBezier b1 t1+> (b2l, b2r) = splitBezier b2 t2 > > adjustC0 :: Point a -> CubicBezier a -> CubicBezier a > adjustC0 p (CubicBezier _ p1 p2 p3) = CubicBezier p p1 p2 p3@@ -875,8 +1013,7 @@ > pointOnCurve :: Double -> DPoint -> CubicBezier Double -> Maybe Double > pointOnCurve tol p c1-> | (t:_) <--> closest c1 p tol,+> | t <- closest c1 p tol, > p2 <- evalBezier c1 t, > vectorDistance p p2 < tol = Just t > | otherwise = Nothing@@ -906,25 +1043,29 @@ > where dist = max (abs $ ld b0) (abs $ ld b3) > ld = lineDistance (Line a0 a3) -=== Finding the on curve point at the X coordinate {#findx}--Solve a cubic equation to find the X coordinate. This should be-converted to a closed form solver for efficiency.--> findX :: Double -> CubicBezier Double -> Double -> Double-> findX x c@(CubicBezier p0 p1 p2 p3) eps =-> head $ bezierFindRoot bez 0 1 $-> bezierParamTolerance c (eps/10)-> where bez = listToBernstein -> (map pointX [p0, p1, p2, p3]) ~--> listToBernstein [x, x, x, x]- Higher level functions ---------------------- > fillFunction :: FillRule -> Int -> Bool-> fillFunction NonZero = (>0)+> fillFunction NonZero = (/=0) > fillFunction EvenOdd = odd+>+> newSweep :: XStruct -> SweepState+> newSweep xStr = SweepState M.empty S.empty undefined xStr+>+> runSweep :: SweepState -> SweepStateM () -> SweepState+#if DEBUG+> runSweep sweep m = +> unsafePerformIO $ do+> hPutStrLn stderr "XSTRUCTBEGIN" +> mapM_ (hPutStrLn stderr . showCurve) $ map (view bezier) $+> concat $ M.elems $ view xStruct sweep+> hPutStrLn stderr "XSTRUCTEND" +> execStateT m sweep+#else+> runSweep sweep m =+> execState m sweep+#endif > -- | `O((n+m)*log(n+m))`, for `n` segments and `m` intersections. > -- Union of paths, removing overlap and rounding to the given@@ -934,14 +1075,21 @@ > -> Double -- ^ Tolerance > -> [ClosedPath Double] > union p fill tol =-> outputPaths out+> outputPaths $ view output $ runSweep sweep $ +> loopEvents (fillFunction fill . fst) tol > where-> out = view output $-> loopEvents (fillFunction fill . fst) tol sweep-> sweep = SweepState M.empty S.empty S.empty xStr+> sweep = newSweep xStr > xStr = makeXStruct (over _1 $ subtract 1) (over _1 (+1)) tol $ > concatMap closedPathCurves p >+> union' p fill tol =+> outputPaths $ view output $ runSweep sweep $ +> loopEvents (fillFunction fill . fst) tol+> where+> sweep = newSweep xStr+> xStr = makeXStruct (over _1 $ subtract 1) (over _1 (+1)) tol p+>+> > -- | `O((n+m)*log(n+m))`, for `n` segments and `m` intersections. > -- Combine paths using the given boolean operation > boolPathOp :: (Bool -> Bool -> Bool) -- ^ operation@@ -951,12 +1099,13 @@ > -> Double -- ^ tolerance > -> [ClosedPath Double] > boolPathOp op p1 p2 fill tol =-> outputPaths $ view output $-> loopEvents isInside tol sweep+> outputPaths $ view output $ runSweep sweep $+> loopEvents isInside tol > where > isInside (a, b) = fillFunction fill a `op` > fillFunction fill b-> sweep = SweepState M.empty S.empty S.empty xStr+> sweep = newSweep xStr+> > xStr = M.unionWith (++) > (makeXStruct > (over _1 (subtract 1))@@ -979,3 +1128,10 @@ > > -- | path exclusion > exclusion = boolPathOp (\a b -> if a then not b else b)+>++handy for debugging: ++>+> -- mkBezier (a, b) (c, d) (e, f) (g, h) = CubicBezier (Point a b) (Point c d) (Point e f) (Point g h)+> --x = [mkBezier (1.7945835892524933,2.0547397002191734)(1.7945835892524933,1.4415012999400492)(2.2917115821994845,0.9443733069930581)(2.9049499824786085,0.9443733069930581),mkBezier (1.7945835892524933,2.0547397002191734)(1.7945835892524933,2.667978238848005)(2.2917115821994845,3.1651062317949963)(2.9049499824786085,3.1651062317949963),mkBezier (2.110720695890836,7.885290987712098)(2.110720695890836,7.516067576338687)(2.410035446506936,7.216752964072295)(2.7792588578803468,7.216752964072295),mkBezier (2.110720695890836,7.885290987712098)(2.110720695890836,8.25451439908551)(2.410035446506936,8.553829149701608)(2.7792588578803468,8.553829149701608),mkBezier (2.7323502826483033,3.689897110469715)(2.7323502826483033,3.1926138891505564)(3.1354780884677687,2.789486083331091)(3.632761309786927,2.789486083331091),mkBezier (2.7323502826483033,3.689897110469715)(2.7323502826483033,4.187180470138581)(3.1354780884677687,4.590308275958047)(3.632761309786927,4.590308275958047),mkBezier (2.7792588578803468,8.553829149701608)(3.148482269253758,8.553829149701608)(3.4477968815201496,8.25451439908551)(3.4477968815201496,7.885290987712098),mkBezier (2.7792588578803468,7.216752964072295)(3.148482269253758,7.216752964072295)(3.4477968815201496,7.516067576338687)(3.4477968815201496,7.885290987712098),mkBezier (2.9049499824786085,3.1651062317949963)(3.5181885211074406,3.1651062317949963)(4.015316514054431,2.667978238848005)(4.015316514054431,2.0547397002191734),mkBezier (2.9049499824786085,0.9443733069930581)(3.5181885211074406,0.9443733069930581)(4.015316514054431,1.4415012999400492)(4.015316514054431,2.0547397002191734),mkBezier (3.632761309786927,4.590308275958047)(4.130044669455794,4.590308275958047)(4.533172336925551,4.187180470138581)(4.533172336925551,3.689897110469715),mkBezier (3.632761309786927,2.789486083331091)(4.130044669455794,2.789486083331091)(4.533172336925551,3.1926138891505564)(4.533172336925551,3.689897110469715),mkBezier (4.282618249080003,8.774851108917153)(4.282618249080003,8.216804972625443)(4.735004084869613,7.764419275185541)(5.293050221161321,7.764419275185541),mkBezier (4.282618249080003,8.774851108917153)(4.282618249080003,9.332897245208862)(4.735004084869613,9.78528308099847)(5.293050221161321,9.78528308099847),mkBezier (5.293050221161321,9.78528308099847)(5.851096357453031,9.78528308099847)(6.3034821932426395,9.332897245208862)(6.3034821932426395,8.774851108917153),mkBezier (5.293050221161321,7.764419275185541)(5.851096357453031,7.764419275185541)(6.3034821932426395,8.216804972625443)(6.3034821932426395,8.774851108917153),mkBezier (5.986614600147038,2.477041530394307)(5.986614600147038,1.712957261642156)(6.606027271186491,1.093544452252995)(7.3701116782883505,1.093544452252995),mkBezier (5.986614600147038,2.477041530394307)(5.986614600147038,3.241125937496166)(6.606027271186491,3.8605386085356193)(7.3701116782883505,3.8605386085356193),mkBezier (6.789870375071436,2.6888791442566093)(6.789870375071436,2.0688856811437524)(7.292474387803418,1.5662816684117704)(7.912467850916275,1.5662816684117704),mkBezier (6.789870375071436,2.6888791442566093)(6.789870375071436,3.3088726073694663)(7.292474387803418,3.8114766201014487)(7.912467850916275,3.8114766201014487),mkBezier (7.128223124391741,0.32094724927962476)(7.128223124391741,-0.19279170133607)(7.54469088573806,-0.6092594626823897)(8.058429836353755,-0.6092594626823897),mkBezier (7.128223124391741,0.32094724927962476)(7.128223124391741,0.8346863382450272)(7.54469088573806,1.251154099591347)(8.058429836353755,1.251154099591347),mkBezier (7.3701116782883505,3.8605386085356193)(8.134195947040501,3.8605386085356193)(8.753608756429662,3.241125937496166)(8.753608756429662,2.477041530394307),mkBezier (7.3701116782883505,1.093544452252995)(8.134195947040501,1.093544452252995)(8.753608756429662,1.712957261642156)(8.753608756429662,2.477041530394307),mkBezier (7.912467850916275,3.8114766201014487)(8.532461314029131,3.8114766201014487)(9.035065188411407,3.3088726073694663)(9.035065188411407,2.6888791442566093),mkBezier (7.912467850916275,1.5662816684117704)(8.532461314029131,1.5662816684117704)(9.035065188411407,2.0688856811437524)(9.035065188411407,2.6888791442566093),mkBezier (8.058429836353755,1.251154099591347)(8.572168925319158,1.251154099591347)(8.988636686665478,0.8346863382450272)(8.988636686665478,0.32094724927962476),mkBezier (8.058429836353755,-0.6092594626823897)(8.572168925319158,-0.6092594626823897)(8.988636686665478,-0.19279170133607)(8.988636686665478,0.32094724927962476),mkBezier (9.002227746912014,6.856573334748986)(9.002227746912014,6.402849518935107)(9.370043726707832,6.035033539139289)(9.823767542521711,6.035033539139289),mkBezier (9.002227746912014,6.856573334748986)(9.002227746912014,7.310297288912573)(9.370043726707832,7.678113130358684)(9.823767542521711,7.678113130358684),mkBezier (9.823767542521711,7.678113130358684)(10.277491496685299,7.678113130358684)(10.645307338131408,7.310297288912573)(10.645307338131408,6.856573334748986),mkBezier (9.823767542521711,6.035033539139289)(10.277491496685299,6.035033539139289)(10.645307338131408,6.402849518935107)(10.645307338131408,6.856573334748986)]
+ Geom2D/CubicBezier/Stroke.hs view
@@ -0,0 +1,323 @@+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}+module Geom2D.CubicBezier.Stroke+ (penCircle, pathToPen, penStrokeOpen, penStrokeClosed, Pen,+ bezierOffset)+ where+import Geom2D+import Geom2D.CubicBezier+import Data.Monoid+++data Pen a = PenEllipse (Transform a) (Transform a) (Transform a)+ | PenPath [PenSegment a]++data PenSegment a = PenCorner !(Point a) !(Point a)+ | PenCurve !(Point a) !(CubicBezier a)++-- | A circular pen with unit radius.+penCircle :: (Floating a) => Pen a+penCircle = PenEllipse idTrans rotate90L rotate90R+{-# SPECIALIZE penCircle :: Pen Double #-}++-- | Create a pen from a path. For predictable results the path+-- should be convex.+pathToPen :: (Floating a) => ClosedPath a -> Pen a+pathToPen (ClosedPath []) = PenPath []+pathToPen (ClosedPath nodes) =+ PenPath $ pathToPen' $ nodes ++ take 2 nodes++pathToPen' :: Num a => [(Point a, PathJoin a)] -> [PenSegment a]+pathToPen' [] = []+pathToPen' [_] = []+pathToPen' [_, _] = []+pathToPen' ((p, JoinLine):tl@((q, JoinLine):_)) =+ PenCorner (q ^-^ p) q : pathToPen' tl++pathToPen' ((_, JoinCurve _ _):tl@((_, JoinLine):_)) =+ pathToPen' tl++pathToPen' ((p, JoinLine):tl@((q1, JoinCurve q2 q3):(q4, _):_)) =+ PenCurve (q1 ^-^ p) (CubicBezier q1 q2 q3 q4) :+ pathToPen' tl++pathToPen' ((_, JoinCurve _ p3):tl@((q1, JoinCurve q2 q3):(q4, _):_)) =+ PenCurve (q1 ^-^ p3) (CubicBezier q1 q2 q3 q4) :+ pathToPen' tl++ +noTranslate :: Num a => Transform a -> Transform a+noTranslate (Transform a b _ c d _) =+ Transform a b 0 c d 0++instance (Floating a, Eq a) => AffineTransform (Pen a) a where+ {-# SPECIALIZE transform :: Transform Double -> Pen Double -> Pen Double #-}+ transform t (PenEllipse trans _ _) =+ let t2@(Transform a b c d e f) = transform t trans+ in case inverse $ noTranslate t2 of+ Nothing -> pathToPen $+ ClosedPath [+ (Point c f ^+^ p, JoinLine),+ (Point c f ^-^ p, JoinLine)]+ where+ p | a /= 0 && b /= 0 =+ sqrt(1 + a*a/(b*b)) *^ Point a d+ | d /= 0 && e /= 0 =+ sqrt(1 + d*d/(e*e)) *^ Point a d+ | a /= 0 = Point (a+d) 0+ | b /= 0 = Point 0 (b+e)+ -- singular point: create tiny pen instead of an error+ | otherwise = Point 1e-5 1e-5+ Just inv ->+ PenEllipse t2 (transform rotate90L inv) (transform rotate90R inv)++ transform t (PenPath segments) =+ PenPath $ map (transformSegment t) segments++transformSegment :: Num b => Transform b -> PenSegment b -> PenSegment b+transformSegment t (PenCorner p q) =+ PenCorner (transform t (q^+^p) ^-^ q') q'+ where q' = transform t q++transformSegment t (PenCurve p c) =+ PenCurve (transform t (cubicC0 c^+^p) ^-^ cubicC0 c') c'+ where c' = transform t c++offsetPoint :: (Floating a) => a -> Point a -> Point a -> Point a+offsetPoint dist start tangent =+ start ^+^ (rotate90L $* dist *^ normVector tangent)++bezierOffsetPoint :: CubicBezier Double -> Double -> Double -> (DPoint, DPoint)+bezierOffsetPoint cb dist t = (offsetPoint dist p p', p')+ where (p, p') = evalBezierDeriv cb t++-- | Calculate an offset path from the bezier curve to within+-- tolerance. If the distance is positive offset to the left,+-- otherwise to the right. A smaller tolerance may require more bezier+-- curves in the path to approximate the offset curve+bezierOffset :: CubicBezier Double -- ^ The curve+ -> Double -- ^ Offset distance.+ -> Maybe Int -- ^ maximum subcurves+ -> Double -- ^ Tolerance.+ -> Bool -- ^ Calculate the curve faster but with+ -- more subcurves+ -> [CubicBezier Double] -- ^ The offset curve+bezierOffset cb dist (Just m) tol faster =+ approximatePathMax m (bezierOffsetPoint cb dist) 15 tol 0 1 faster++bezierOffset cb dist Nothing tol faster =+ approximatePath (bezierOffsetPoint cb dist) 15 tol 0 1 faster++penOffset :: Pen Double -> Point Double -> Point Double+penOffset (PenEllipse trans leftInv _) dir =+ transform trans $ normVector $ leftInv $* dir++penOffset (PenPath segments) dir =+ pathOffsetPoint (cycle segments) dir++penOffsetFun :: Pen Double -> (Double -> (DPoint, DPoint)) -> Double -> (Point Double, Point Double)+penOffsetFun pen f t =+ (px ^+^ penOffset pen px', px')+ where+ (px, px') = f t++firstPoint :: PenSegment a -> Point a+firstPoint (PenCorner _ p) = p+firstPoint (PenCurve _ c) = cubicC0 c++pathOffsetPoint :: [PenSegment Double] -> Point Double -> Point Double+pathOffsetPoint (PenCorner c p:b:rest) dir+ | vectorCross dir c > 0 = pathOffsetPoint (b:rest) dir+ | vectorCross dir (firstPoint b ^-^ p) > 0 = p+ | otherwise = pathOffsetPoint (b:rest) dir+ +pathOffsetPoint (PenCurve c curve@(CubicBezier p1 p2 p3 p4):b:rest) dir+ | vectorCross dir c > 0 = pathOffsetPoint (b:rest) dir+ | vectorCross dir (p2 ^-^ p1) > 0 = p1+ | vectorCross dir (p3 ^-^ p4) > 0 =+ case findBezierTangent dir curve of+ (t:_) -> evalBezier curve t+ [] -> p4+ | vectorCross dir (firstPoint b ^-^ p4) > 0 = p4+ | otherwise = pathOffsetPoint (b:rest) dir++pathOffsetPoint _ _ = error "unexpected end of list"++segDirs :: [(DPoint, PathJoin Double)] -> Point Double -> [(DPoint, DPoint)]+segDirs [] _ = []+segDirs [(p, JoinLine)] q = [(dp, dp)]+ where dp = q ^-^ p+segDirs [(p1, JoinCurve p2 p3 )] p4 = [(p2 ^-^ p1, p4 ^-^ p3)]+segDirs ((p, JoinLine):r@((q, _):_)) s = (dp, dp): segDirs r s+ where dp = q ^-^ p+segDirs ((p1, JoinCurve p2 p3 ):r@((p4,_):_)) q = (p2 ^-^ p1, p4 ^-^ p3):segDirs r q++penStrokeOpen :: Int -> Double -> Bool -> Pen Double -> OpenPath Double -> [ClosedPath Double]+penStrokeOpen samples tol fast pen (OpenPath segments p) =+ union [closeOpenPath path] NonZero tol+ where+ dirs = segDirs segments (fst $ head segments)+ fdirs = map fst (tail dirs)+ fd = fst $ head dirs+ ld = snd $ last dirs+ ldirs = map snd dirs + pts = map fst (tail segments) ++ [p]+ leftJoins = zipWith (penJoinLeft pen) ldirs fdirs+ leftStrokes = zipWith (strokeLeft samples tol fast pen) segments pts+ rightJoins = zipWith (penJoinRight pen) ldirs fdirs+ rightStrokes = zipWith (strokeRight samples tol fast pen) segments pts+ path =+ mconcat $+ penJoinLeft pen (turnAround fd) fd :+ interleave leftStrokes leftJoins +++ penJoinLeft pen ld (turnAround ld) :+ reverse (interleave rightStrokes rightJoins)++interleave :: [a] -> [a] -> [a]+interleave [] xs = xs+interleave xs [] = xs+interleave (x:xs) (y:ys) = x:y:interleave xs ys ++--penStrokeClosed :: ClosedPath Double -> Pen Double -> Double -> [ClosedPath Double]+penStrokeClosed :: Int -> Double -> Bool -> Pen Double -> ClosedPath Double+ -> [ClosedPath Double]+penStrokeClosed _ _ _ _ (ClosedPath []) = [ClosedPath []]+penStrokeClosed samples tol fast pen (ClosedPath segments) =+ union [closeOpenPath leftPath, closeOpenPath rightPath] NonZero tol+ where+ dirs = segDirs segments (fst $ head segments)+ fdirs = map fst (tail dirs) ++ [fst (head dirs)]+ ldirs = map snd dirs+ pts = map fst (tail segments) ++ [fst (head segments)]+ leftJoins = zipWith (penJoinLeft pen) ldirs fdirs+ leftStrokes = zipWith (strokeLeft samples tol fast pen) segments pts+ rightJoins = zipWith (penJoinRight pen) ldirs fdirs+ rightStrokes = zipWith (strokeRight samples tol fast pen) segments pts+ leftPath =+ mconcat $ interleave leftStrokes leftJoins+ rightPath =+ mconcat $ reverse $ interleave rightStrokes rightJoins++strokeLeft :: Int -> Double -> Bool -> Pen Double -> (DPoint, PathJoin Double) -> DPoint -> OpenPath Double+strokeLeft _ _ _ pen (p, JoinLine) q =+ OpenPath [(p ^+^ offset, JoinLine)] (q ^+^ offset)+ where offset = penOffset pen (q ^-^ p)++strokeLeft samples tol fast pen (p1, JoinCurve p2 p3) p4 =+ curvesToOpen $ approximatePath+ (penOffsetFun pen (evalBezierDeriv (CubicBezier p1 p2 p3 p4)))+ samples tol 0 1 fast++strokeRight :: Int -> Double -> Bool -> Pen Double -> (DPoint, PathJoin Double) -> DPoint -> OpenPath Double+strokeRight _ _ _ pen (p, JoinLine) q =+ OpenPath [(q ^+^ offset, JoinLine)] (p ^+^ offset)+ where offset = penOffset pen (p ^-^ q)++strokeRight samples tol fast pen (p1, JoinCurve p2 p3) p4 =+ curvesToOpen $ approximatePath+ (penOffsetFun pen (evalBezierDeriv (CubicBezier p4 p3 p2 p1)))+ samples tol 0 1 fast++penJoinLeft :: Pen Double -> DPoint -> DPoint -> OpenPath Double+penJoinLeft = penJoin++penJoinRight :: Pen Double -> DPoint -> DPoint -> OpenPath Double+penJoinRight pen from to = penJoin pen (turnAround to) (turnAround from)++ellipticArc :: Transform Double -> Transform Double+ -> Point Double -> Point Double -> CubicBezier Double+ellipticArc trans leftInv from to =+ trans $* bezierArc+ (vectorAngle $ leftInv $* from)+ (vectorAngle $ leftInv $* to)++segmentsToPath :: (Eq a) => [PenSegment a] -> OpenPath a+segmentsToPath [PenCorner _ q] =+ OpenPath [] q+segmentsToPath [PenCurve _ (CubicBezier p1 p2 p3 p4)] =+ OpenPath [(p1, JoinCurve p2 p3)] p4+ +segmentsToPath (PenCorner _ p:r) =+ consOpenPath p JoinLine (segmentsToPath r)++segmentsToPath (PenCurve _ (CubicBezier p1 p2 p3 p4):r) =+ consOpenPath p1 (JoinCurve p2 p3) $+ case r of+ (PenCurve _ (CubicBezier q1 _ _ _):_)+ | p4 /= q1 -> consOpenPath p4 JoinLine $ segmentsToPath r+ _ -> segmentsToPath r++segmentsToPath [] = emptyOpenPath ++emptyOpenPath :: OpenPath a+emptyOpenPath = OpenPath [] (error "empty path")+ +penJoin :: Pen Double -> Point Double+ -> Point Double -> OpenPath Double+penJoin pen@(PenEllipse trans leftInv _) from to+ | dir == 0 = emptyOpenPath+ | dir > 0 &&+ sameQuadrant from to =+ curvesToOpen [ellipticArc trans leftInv from to]+ | otherwise =+ curvesToOpen [ellipticArc trans leftInv from next] <>+ penJoin pen next to+ where next = nextVector from+ dir = vectorCross from to++penJoin (PenPath segments) from to =+ segmentsToPath $+ nextSegments (firstSegment (cycle segments) from) to++firstSegment :: [PenSegment Double] -> Point Double -> [PenSegment Double]+firstSegment segments@(PenCorner c _:q:rest) from+ | vectorCross from c > 0 =+ firstSegment (q:rest) from+ | otherwise = segments++firstSegment segments@(PenCurve c curve@(CubicBezier p1 p2 p3 p4):q:rest) from+ | vectorCross from c > 0 = firstSegment (q:rest) from+ | vectorCross from (p2 ^-^ p1) > 0 = segments+ | vectorCross from (p4 ^-^ p3) > 0 =+ case findBezierTangent from curve of+ (t:_) -> PenCurve from (snd (splitBezier curve t)):q:rest+ _ -> q:rest+ | vectorCross from (firstPoint q ^-^ p4) > 0 =+ PenCorner (firstPoint q ^-^ p4) p4:q:rest+ | otherwise = firstSegment (q:rest) from++firstSegment _ _ = error "firstsegment: finite list" ++nextSegments :: [PenSegment Double] -> Point Double -> [PenSegment Double]+nextSegments (PenCorner c p:q:rest) to+ | vectorCross to c > 0 =+ PenCorner c p: nextSegments (q:rest) to+ | otherwise = []++nextSegments (pc@(PenCurve c curve@(CubicBezier p1 p2 p3 p4)):q:rest) to + | vectorCross to c > 0 = pc: nextSegments (q:rest) to+ | vectorCross to (p2 ^-^ p1) > 0 = []+ | vectorCross to (p4 ^-^ p3) > 0 =+ case findBezierTangent to curve of+ (t:_) -> [PenCurve c (fst (splitBezier curve t))]+ _ -> []+ | vectorCross to (firstPoint q ^-^ p4) > 0 =+ [PenCorner (firstPoint q ^-^ p4) p4]+ | otherwise = pc:firstSegment (q:rest) to++nextSegments _ _ = error "nextSegments: finite list"++sameQuadrant :: (Num a, Eq a) => Point a -> Point a -> Bool+sameQuadrant v w =+ signum (pointX v) /= -signum (pointX w) &&+ signum (pointY v) /= -signum (pointY w)++nextVector :: (Ord a1, Num a1, Num a) => Point a1 -> Point a+nextVector v+ | pointX v >= 0 &&+ pointY v > 0 = Point 1 0+ | pointX v > 0 &&+ pointY v <= 0 = Point 0 (-1)+ | pointX v <= 0 &&+ pointY v < 0 = Point (-1) 0+ | otherwise = Point 0 1+
Math/BernsteinPoly.hs view
@@ -31,13 +31,13 @@ toScaled (BernsteinPoly v) = ScaledPoly $ V.zipWith (*) v $ binCoeff $ V.length v - 1-{-# NOINLINE[2] toScaled #-}+{-# NOINLINE [2] toScaled #-} fromScaled :: (Unbox a, Fractional a) => ScaledPoly a -> BernsteinPoly a fromScaled (ScaledPoly v) = BernsteinPoly $ V.zipWith (/) v $ binCoeff $ V.length v - 1-{-# NOINLINE[2] fromScaled #-}+{-# NOINLINE [2] fromScaled #-} -- | Create a bernstein polynomail from a list of coëfficients. listToBernstein :: (Unbox a, Num a) => [a] -> BernsteinPoly a@@ -198,7 +198,7 @@ -- maximum of either degrees. (~+) :: (Unbox a, Fractional a) => BernsteinPoly a -> BernsteinPoly a -> BernsteinPoly a-(toScaled -> a) ~+ (toScaled -> b) = fromScaled $ addScaled a b+a ~+ b = fromScaled $ addScaled (toScaled a) (toScaled b) {-# INLINE (~+) #-} subScaled :: (Unbox a, Num a) => ScaledPoly a -> ScaledPoly a -> ScaledPoly a
cubicbezier.cabal view
@@ -1,8 +1,8 @@ Name: cubicbezier-Version: 0.5.0.0+Version: 0.6.0.0 Synopsis: Efficient manipulating of 2D cubic bezier curves. Category: Graphics, Geometry, Typography-Copyright: Kristof Bastiaensen (2014)+Copyright: Kristof Bastiaensen (2017) Stability: Unstable License: BSD3 License-file: LICENSE@@ -20,14 +20,22 @@ The module "Geom2D.CubicBezier" exports all the cubic bezier functions. The module "Geom2D" contains general 2D geometry functions and transformations. -source-repository head+Source-repository head type: git location: https://github.com/kuribas/cubicbezier +Flag Debug+ description: Enable debug messages+ Manual: True+ default: False+ Library+ if flag(Debug)+ cpp-options: -DDEBUG Ghc-options: -Wall Build-depends: base >= 3 && < 5, containers >= 0.5.3, integration >= 0.1.1, vector >= 0.10,- matrices >= 0.4.1, microlens >= 0.1.2, microlens-th >= 0.1.2, microlens-mtl >= 0.1.2, mtl >= 2.1.1+ matrices >= 0.4.1, microlens >= 0.1.2, microlens-th >= 0.1.2, microlens-mtl >= 0.1.2, mtl >= 2.1.1,+ fast-math >= 1.0.0, vector-space >= 0.10.4 Exposed-Modules: Geom2D Geom2D.CubicBezier@@ -38,6 +46,7 @@ Geom2D.CubicBezier.Overlap Geom2D.CubicBezier.Intersection Geom2D.CubicBezier.MetaPath+ Geom2D.CubicBezier.Stroke Math.BernsteinPoly Geom2D.CubicBezier.Numeric