packages feed

cubicbezier (empty) → 0.1.0

raw patch · 12 files changed

+1397/−0 lines, 12 filesdep +basedep +containersdep +integrationsetup-changed

Dependencies added: base, containers, integration

Files

+ Geom2D.hs view
@@ -0,0 +1,147 @@+-- | Basic 2 dimensional geometry functions.+module Geom2D where++infixl 6 ^+^, ^-^+infixl 7 *^, ^*, ^/+infixr 5 $*++data Point = Point {+  pointX :: Double,+  pointY :: Double}++instance Show Point where+  show (Point x y) =+    "Point " ++ show x ++ " " ++ show y++-- | A transformation (x, y) -> (ax + by + c, dx + ey + d)+data Transform = Transform {+  xformA :: Double,+  xformB :: Double,+  xformC :: Double,+  xformD :: Double,+  xformE :: Double,+  xformF :: Double }+               deriving Show++data Line = Line Point Point+data Polygon = Polygon [Point]++class AffineTransform a where+  transform :: Transform -> a -> a++instance AffineTransform Transform where+  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')+    +instance AffineTransform Point where+  transform (Transform a b c d e f) (Point x y) =+    Point (a*x + b*y + c) (d*x + e*y + f)++instance AffineTransform Polygon where+  transform t (Polygon p) = Polygon $ map (transform t) p++-- | Operator for applying a transformation.+($*) :: AffineTransform a => Transform -> a -> a+t $* p = transform t p++-- | Calculate the inverse of a transformation.+inverse :: Transform -> Maybe Transform+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)+         (-(b*c + e*f)/det)++-- | Return the parameters (a, b, c) for the normalised equation+-- of the line: @a*x + b*y + c = 0@.+lineEquation :: Line -> (Double, Double, Double)+lineEquation (Line (Point x1 y1) (Point x2 y2)) = (a, b, c)+  where a = a' / d+        b = b' / d+        c = -(y1*b' + x1*a') / d+        a' = y1 - y2+        b' = x2 - x1+        d = sqrt(a'*a' + b'*b')++-- | Return the the distance from a point to the line.+lineDistance :: Line -> Point -> Double+lineDistance l = \(Point x y) -> a*x + b*y + c+  where (a, b, c) = lineEquation l++-- | The lenght of the vector.+vectorMag :: Point -> Double+vectorMag (Point x y) = sqrt(x*x + y*y)++-- | The angle of the vector, in the range @(-'pi', 'pi']@.+vectorAngle :: Point -> Double+vectorAngle (Point 0.0 0.0) = 0.0+vectorAngle (Point x y) = atan2 y x++-- | The unitvector with the given angle.+dirVector :: Double -> Point+dirVector angle = Point (cos angle) (sin angle)++-- | The unit vector with the same direction.+normVector :: Point -> Point+normVector p@(Point x y) = Point (x/l) (y/l)+  where l = vectorMag p++-- | Scale vector by constant.+(*^) :: Double -> Point -> Point+s *^ (Point x y) = Point (s*x) (s*y)++-- | Scale vector by reciprocal of constant.+(^/) :: Point -> Double -> Point+(Point x y) ^/ s = Point (x/s) (y/s)++-- | Scale vector by constant, with the arguments swapped.+(^*) :: Point -> Double -> Point+p ^* s = s *^ p++-- | Add two vectors.+(^+^) :: Point -> Point -> Point+(Point x1 y1) ^+^ (Point x2 y2) = Point (x1+x2) (y1+y2)++-- | Subtract two vectors.+(^-^) :: Point -> Point -> Point+(Point x1 y1) ^-^ (Point x2 y2) = Point (x1-x2) (y1-y2)++-- | Dot product of two vectors.+(^.^) :: Point -> Point -> Double+(Point x1 y1) ^.^ (Point x2 y2) = x1*x2 + y1*y2++-- | Cross product of two vectors.+vectorCross :: Point -> Point -> Double+vectorCross (Point x1 y1) (Point x2 y2) = x1*y2 - y1*x2++-- | Distance between two vectors.+vectorDistance :: Point -> Point -> Double+vectorDistance p q = vectorMag (p^-^q)++-- | Interpolate between two vectors.+interpolateVector :: Point -> Point -> Double -> Point+interpolateVector a b t = t*^b ^+^ (1-t)*^a++-- | Create a transform that rotates by the angle of the given vector+-- with the x-axis+rotateVec :: Point -> Transform+rotateVec v = Transform x (-y) 0 y x 0+  where Point x y = normVector v++-- | Create a transform that rotates by the given angle (radians).+rotate :: Double -> Transform+rotate a = Transform (cos a) (negate $ sin a) 0+           (sin a) (cos a) 0++-- | Rotate vector 90 degrees left.+rotate90L :: Transform+rotate90L = rotateVec (Point 0 1)++-- | Rotate vector 90 degrees right.+rotate90R :: Transform+rotate90R = rotateVec (Point 0 (-1))++-- | Create a transform that translates by the given vector.+translate :: Point -> Transform+translate (Point x y) = Transform 1 0 x 0 1 y+
+ Geom2D/CubicBezier.hs view
@@ -0,0 +1,18 @@+-- | Export all the cubic bezier functions, but not the numeric helper functions++module Geom2D.CubicBezier+       (module Geom2D.CubicBezier.Basic,+        module Geom2D.CubicBezier.Approximate,+        module Geom2D.CubicBezier.Outline,+        module Geom2D.CubicBezier.Curvature,+        module Geom2D.CubicBezier.Intersection+       ) where++import Geom2D.CubicBezier.Basic+import Geom2D.CubicBezier.Approximate+import Geom2D.CubicBezier.Outline+import Geom2D.CubicBezier.Curvature+import Geom2D.CubicBezier.Intersection+       +       +
+ Geom2D/CubicBezier/Approximate.hs view
@@ -0,0 +1,128 @@+module Geom2D.CubicBezier.Approximate (+  approximateCurve, approximateCurveWithParams)+       where+import Geom2D+import Geom2D.CubicBezier.Numeric+import Geom2D.CubicBezier.Basic+import Data.Function+import Data.List+import Data.Maybe++-- | @approximateCurve b pts eps@ finds the least squares fit of a bezier+-- curve to the points @pts@.  The resulting bezier has the same first+-- and last control point as the curve @b@, and have tangents colinear with @b@.+-- return the curve, the parameter with maximum error, and maximum error.+-- Calculate to withing eps tolerance.++approximateCurve :: CubicBezier -> [Point] -> Double -> (CubicBezier, Double, Double)+approximateCurve curve@(CubicBezier p1 _ _ p4) pts eps =+  approximateCurveWithParams curve pts (approximateParams curve p1 p4 pts) eps++-- | Like approximateCurve, but also takes an initial guess of the+-- parameters closest to the points.  This might be faster if a good+-- guess can be made.++approximateCurveWithParams :: CubicBezier -> [Point] -> [Double] -> Double -> (CubicBezier, Double, Double)+approximateCurveWithParams curve pts ts eps =+  let (c, newTs) = fromMaybe (curve, ts) $+                   approximateCurve' curve pts ts 40 (bezierParamTolerance curve eps) 1+      curvePts   = map (evalBezier c) newTs+      distances  = zipWith vectorDistance pts curvePts+      (t, maxError) = maximumBy (compare `on` snd) (zip ts distances)+  in (c, t, maxError)++add6 (a, b, c, d, e, f) (a', b', c', d', e', f') =+  (a+a', b+b', c+c', d+d', e+e', f+f')+++-- find the least squares between the points p_i and B(t_i) for+-- bezier curve B, where pts contains the points p_i and ts+-- the values of t_i .+-- The tangent at the beginning and end is maintained.+-- Since the start and end point remains the same,+-- we need to find the new value of p2' = p1 + alpha1 * (p2 - p1)+-- and p3' = p4 + alpha2 * (p3 - p4)+-- minimizing (sum |B(t_i) - p_i|^2) gives a linear equation+-- with two unknown values (alpha1 and alpha2), which can be+-- solved easily+leastSquares :: CubicBezier -> [Point] -> [Double] -> Maybe CubicBezier+leastSquares (CubicBezier (Point p1x p1y) (Point p2x p2y) (Point p3x p3y) (Point p4x p4y)) pts ts = let+  calcParams t (Point px py)  = let+    t2 = t * t; t3 = t2 * t+    ax = 3 * (p2x - p1x) * (t3 - 2 * t2 + t)+    ay = 3 * (p2y - p1y) * (t3 - 2 * t2 + t)+    bx = 3 * (p3x - p4x) * (t2 - t3)+    by = 3 * (p3y - p4y) * (t2 - t3)+    cx = (p4x - p1x) * (3 * t2 - 2 * t3) + p1x - px+    cy = (p4y - p1y) * (3 * t2 - 2 * t3) + p1y - py+    in (ax * ax + ay * ay,+        ax * bx + ay * by,+        ax * cx + ay * cy,+        bx * ax + by * ay,+        bx * bx + by * by,+        bx * cx + by * cy)+  (a, b, c, d, e, f) = foldl1' add6 $ zipWith calcParams ts pts+  in do (alpha1, alpha2) <- solveLinear2x2 a b c d e f+        let cp1 = Point (alpha1 * (p2x - p1x) + p1x) (alpha1 * (p2y - p1y) + p1y)+            cp2 = Point (alpha2 * (p3x - p4x) + p4x) (alpha2 * (p3y - p4y) + p4y)+        Just $ CubicBezier (Point p1x p1y) cp1 cp2 (Point p4x p4y)++-- calculate the least Squares bezier curve by choosing approximate values+-- of t, and iterating again with an improved estimate of t, by taking the+-- the values of t for which the points are closest to the curve++approximateCurve' :: CubicBezier -> [Point] -> [Double] -> Int -> Double -> Double -> Maybe (CubicBezier, [Double])+approximateCurve' curve pts ts maxiter eps prevDeltaT = do+  newCurve <- leastSquares curve pts ts+  let deltaTs = zipWith (calcDeltaT newCurve) pts ts+      ts' = map (max 0 . min 1) $ zipWith (-) ts deltaTs+  newCurve <- leastSquares curve pts ts'+  let deltaTs' = zipWith (calcDeltaT newCurve) pts ts'+      newTs = interpolateTs ts ts' deltaTs deltaTs'+      thisDeltaT = maximum $ map abs $ zipWith (-) newTs ts+  if maxiter < 1 ||+     -- Because convergence may be slow initially, make sure it is converging:+     (prevDeltaT < eps/2  && thisDeltaT < prevDeltaT / 2)+    then do c <- leastSquares curve pts newTs+            return (c, newTs)+    else approximateCurve' curve pts newTs (maxiter - 1) eps thisDeltaT++-- improve convergence by making a better estimate for t+-- it is based on the observation that the ratio  +-- r = dt_2 / dt_1, with dt_2 = t_2 - t_1 and dt_1 = t_1 - t_0+-- for successive approximations of t changes little.+-- The infinite sum (dt_1 + dt_1 * r + dt_1 * r^2 + dt_1 * r^3 ...)+-- can easily be calculated by dt_1 * (1 / (1 - r))+-- which becomes dt_1^2 / (dt_1 - dt_2)+-- Only do this if it appears to converge for all values of t+-- If the value of t changes too much keep the old value.+-- This improves the convergence by a factor of about 10+interpolateTs :: [Double] -> [Double] -> [Double] -> [Double] -> [Double]+interpolateTs ts ts' deltaTs deltaTs' =+  map (max 0 . min 1) (+    if all id $ zipWith (\dT dT' -> dT * dT' > 0 && dT' / dT < 1) deltaTs deltaTs'+    then zipWith3 (\t dT dT' -> let+                      newDt = (dT * dT / (dT - dT'))+                      in t - (if abs newDt > 0.2 then dT' else newDt)) ts deltaTs deltaTs'+    else zipWith (-) ts' deltaTs')++-- approximate t by calculating the distances between all points+-- and dividing by the total sum+approximateParams :: CubicBezier -> Point -> Point -> [Point] -> [Double]+approximateParams cb start end pts = let+  segments = start : (pts ++ [end])+  dists = zipWith vectorDistance segments (tail segments)+  total = sum dists+  improve p t = t - calcDeltaT cb p t+  in zipWith improve pts $ map (/ total) $ scanl1 (+) dists++-- find a value of t where B(t) is closer between the bezier curve and+-- the point (ptx, pty), by solving f' = 0 where+-- f(t) = (X(t) - x)^2 + (Y(t) - y)^2, the square of the distance between the bezier and the point+-- the reduction of t is one iteration of Newton Raphson:  f'(t)/f''(t)+-- using more iterations doesn't appear to give an improvement+-- See Curve Fitting with Piecewise Parametric Cubics by Stone & Plass+calcDeltaT curve (Point ptx pty) t = let+  [Point bezx bezy, Point dbezx dbezy, Point ddbezx ddbezy, _] = evalBezierDerivs curve t+  in ((bezx - ptx) * dbezx + (bezy - pty) * dbezy) /+     (dbezx * dbezx + dbezy * dbezy + (bezx - ptx) * ddbezx + (bezy - pty) * ddbezy)
+ Geom2D/CubicBezier/Basic.hs view
@@ -0,0 +1,208 @@+{-# LANGUAGE BangPatterns #-}+module Geom2D.CubicBezier.Basic+       (CubicBezier (..), PathJoin (..), Path (..), AffineTransform (..), +        bezierParam, bezierParamTolerance, reorient, bezierToBernstein,+        evalBezier, evalBezierDeriv, evalBezierDerivs, findBezierTangent,+        bezierHoriz, bezierVert, findBezierInflection, findBezierCusp,+        arcLength, arcLengthParam, splitBezier, bezierSubsegment, splitBezierN,+        colinear)+       where+import Geom2D+import Geom2D.CubicBezier.Numeric+import Math.BernsteinPoly+import Numeric.Integration.TanhSinh++data CubicBezier = CubicBezier {+  bezierC0 :: Point,+  bezierC1 :: Point,+  bezierC2 :: Point,+  bezierC3 :: Point} deriving Show++data PathJoin = JoinLine | JoinCurve Point Point+data Path = Path Point [(PathJoin, Point)]++instance AffineTransform CubicBezier where+  transform t (CubicBezier c0 c1 c2 c3) =+    CubicBezier (transform t c0) (transform t c1) (transform t c2) (transform t c3)++++-- | Return True if the param lies on the curve, iff it's in the interval @[0, 1]@.+bezierParam :: Double -> Bool+bezierParam t = t >= 0 && t <= 1++-- | Convert a tolerance from the codomain to the domain of the bezier curve.+-- Should be good enough, but may not hold for high very tolerance values.++-- The magnification of error from the domain to the codomain of the+-- curve approaches the length of the tangent for small errors.  We+-- can use the maximum of the convex hull of the derivative, and double it to+-- have some margin for larger values.+bezierParamTolerance :: CubicBezier -> Double -> Double+bezierParamTolerance (CubicBezier p1 p2 p3 p4) eps = eps / maxDist+  where +    maxDist = 6 * maximum [vectorDistance p1 p2,+                           vectorDistance p2 p3,+                           vectorDistance p3 p4]++-- | Reorient to the curve B(1-t).+reorient :: CubicBezier -> CubicBezier+reorient (CubicBezier p0 p1 p2 p3) = CubicBezier p3 p2 p1 p0 ++-- | Give the bernstein polynomial for each coordinate.+bezierToBernstein :: CubicBezier -> (BernsteinPoly, BernsteinPoly)+bezierToBernstein (CubicBezier a b c d) = (listToBernstein $ map pointX coeffs,+                                           listToBernstein $ map pointY coeffs)+  where coeffs = [a, b, c, d]++-- | Calculate a value on the curve.+evalBezier :: CubicBezier -> Double -> Point+evalBezier b t = Point (bernsteinEval x t) (bernsteinEval y t)+  where (x, y) = bezierToBernstein b++-- | Calculate a value and the first derivative on the curve.+evalBezierDeriv :: CubicBezier -> Double -> (Point, Point)+evalBezierDeriv b =+  let (px, py) = bezierToBernstein b+      px' = bernsteinDeriv px+      py' = bernsteinDeriv py+  in \t -> (Point (bernsteinEval px t) (bernsteinEval py t),+            Point (bernsteinEval px' t) (bernsteinEval py' t))++-- | Calculate a value and all derivatives on the curve.+evalBezierDerivs :: CubicBezier -> Double -> [Point]+evalBezierDerivs b t = zipWith Point (bernsteinEvalDerivs px t)+                       (bernsteinEvalDerivs py t)+  where (px, py) = bezierToBernstein b++-- | @findBezierTangent p b@ finds the parameters where+-- the tangent of the bezier curve @b@ has the same direction as vector p.++-- Use the formula tx * B'y(t) - ty * B'x(t) = 0 where+-- B'x is the x value of the derivative of the Bezier curve.+findBezierTangent :: Point -> CubicBezier -> [Double]+findBezierTangent (Point tx ty) (CubicBezier (Point x0 y0) (Point x1 y1) (Point x2 y2) (Point x3 y3)) = +  filter bezierParam $ quadraticRoot a b c+    where+      a = tx*((y3 - y0) + 3*(y1 - y2)) - ty*((x3 - x0) + 3*(x1 - x2))+      b = 2*(tx*((y2 + y0) - 2*y1) - ty*((x2 + x0) - 2*x1))+      c = tx*(y1 - y0) - ty*(x1 - x0)++-- | Find the parameter where the bezier curve is horizontal.+bezierHoriz :: CubicBezier -> [Double]+bezierHoriz = findBezierTangent (Point 1 0)++-- | Find the parameter where the bezier curve is vertical.+bezierVert :: CubicBezier -> [Double]+bezierVert = findBezierTangent (Point 0 1)++-- | Find inflection points on the curve.++-- Use the formula B''x(t) * B'y(t) - B''y(t) * B'x(t) = 0+-- with B'x(t) the x value of the first derivative at t,+-- B''y(t) the y value of the second derivative at t+findBezierInflection :: CubicBezier -> [Double]+findBezierInflection (CubicBezier (Point x0 y0) (Point x1 y1) (Point x2 y2) (Point x3 y3)) =+  filter bezierParam $ quadraticRoot a b c+    where+      ax = x1 - x0+      bx = x3 - x1 - ax+      cx = x3 - x2 - ax - 2*bx+      ay = y1 - y0+      by = y2 - y1 - ay+      cy = y3 - y2 - ay - 2*by+      a = bx*cy - by*cx+      b = ax*cy - ay*cx+      c = ax*by - ay*bx++-- | Find the cusps of a bezier.++-- find a cusp.  We look for points where the tangent is both horizontal+-- and vertical, which is only true for the zero vector.+findBezierCusp :: CubicBezier -> [Double]+findBezierCusp b = filter vertical $ bezierHoriz b+  where vertical = (== 0) . pointY . snd . evalBezierDeriv b++-- | @arcLength c t tol finds the arclength of the bezier c at t, within given tolerance tol.++arcLength :: CubicBezier -> Double -> Double -> Double+arcLength b@(CubicBezier c0 c1 c2 c3) t eps =+  if eps / maximum [vectorDistance c0 c1,+                    vectorDistance c1 c2,+                    vectorDistance c2 c3] > 1e-10+  then (signum t *) $ fst $+       arcLengthEstimate (fst $ splitBezier b t) eps+  else arcLengthQuad b t eps++arcLengthQuad :: CubicBezier -> Double -> Double -> Double+arcLengthQuad b t eps = result $ absolute eps $+                        trap distDeriv 0 t+  where distDeriv t' = vectorMag $ snd $ evalD t'+        evalD = evalBezierDeriv b ++outline (CubicBezier c0 c1 c2 c3) =+  sum [vectorDistance c0 c1,+       vectorDistance c1 c2,+       vectorDistance c2 c3]++arcLengthEstimate :: CubicBezier -> Double -> (Double, (Double, Double))+arcLengthEstimate b eps = (arclen, (estimate, ol))+  where+    estimate = (4*(olL+olR) - ol) / 3+    (bl, br) = splitBezier b 0.5+    ol = outline b+    (arcL, (estL, olL)) = arcLengthEstimate bl eps+    (arcR, (estR, olR)) = arcLengthEstimate br eps+    arclen | (abs(estL + estR - estimate) < eps) = estL + estR+           | otherwise = arcL + arcR++-- | arcLengthParam c len tol finds the parameter where the curve c has the arclength len,+-- within tolerance tol.+arcLengthParam b len eps =+  arcLengthP b len ol (len/ol) 1 eps+  where ol = outline b++-- Use the Newton rootfinding method.  Start with large tolerance+-- values, and decrease tolerance as we go closer to the root.+arcLengthP !b !len !tot !t !dt !eps+  | abs diff < eps = t - newDt+  | otherwise = arcLengthP b len tot (t - newDt) newDt eps+  where diff = arcLength b t (max (abs (dt*tot/50)) (eps/2)) - len+        newDt = diff / vectorMag (snd $ evalBezierDeriv b t)++-- | Split a bezier curve into two curves.+splitBezier :: CubicBezier -> Double -> (CubicBezier, CubicBezier)+splitBezier (CubicBezier a b c d) t =+  let ab = interpolateVector a b t+      bc = interpolateVector b c t+      cd = interpolateVector c d t+      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)++-- | Return the subsegment between the two parameters.+bezierSubsegment :: CubicBezier -> Double -> Double -> CubicBezier+bezierSubsegment b t1 t2 +  | t1 > t2   = bezierSubsegment b t2 t1+  | otherwise = snd $ flip splitBezier (t1/t2) $+                fst $ splitBezier b t2++-- | Split a bezier curve into a list of beziers+-- The parameters should be in ascending order or+-- the result is unpredictable.+splitBezierN :: CubicBezier -> [Double] -> [CubicBezier]+splitBezierN c [] = [c]+splitBezierN c [t] = [a, b] where+  (a, b) = splitBezier c t+splitBezierN c (t:u:rest) =+  bezierSubsegment c 0 t :+  bezierSubsegment c t u :+  tail (splitBezierN c $ u:rest)++-- | Return True if all the control points are colinear within tolerance.+colinear :: CubicBezier -> Double -> Bool+colinear (CubicBezier a b c d) eps =+  abs (ld b) < eps && abs (ld c) < eps+  where ld = lineDistance (Line a d)+
+ Geom2D/CubicBezier/Curvature.hs view
@@ -0,0 +1,63 @@+module Geom2D.CubicBezier.Curvature+       (curvature, radiusOfCurvature, curvatureExtrema, findRadius)+where+import Geom2D+import Geom2D.CubicBezier.Basic+import Geom2D.CubicBezier.Intersection+import Math.BernsteinPoly++-- | Curvature of the Bezier curve.+curvature :: CubicBezier -> Double -> Double+curvature b t+  | t == 0 = curvature' b+  | t == 1 = negate $ curvature' $ reorient b+  | t < 0.5 = curvature' $ snd $ splitBezier b t+  | otherwise = negate $ curvature' $ reorient $ fst $ splitBezier b t++curvature' (CubicBezier c0 c1 c2 c3) = 2/3 * b/a^3+  where +    a = vectorDistance c1 c0+    b = (c1^-^c0) `vectorCross` (c2^-^c1)++-- | Radius of curvature of the Bezier curve.  This+-- is the reciprocal of the curvature.+radiusOfCurvature :: CubicBezier -> Double -> Double+radiusOfCurvature b t = 1 / curvature b t++extrema :: CubicBezier -> BernsteinPoly+extrema (CubicBezier p0 p1 p2 p3) =+  let bez = [p0, p1, p2, p3]+      x' = bernsteinDeriv $ listToBernstein $ map pointX bez+      y' = bernsteinDeriv $ listToBernstein $ map pointY bez+      x'' = bernsteinDeriv x'+      y'' = bernsteinDeriv y'+      x''' = bernsteinDeriv x''+      y''' = bernsteinDeriv y''+  in -- (y'^2 + x'^2) * (x'*y''' - y'*x''') -+     -- 3 * (x'*y'' - y'*x'') * (y'*y'' + x'*x'')+   (y'~*y' ~+ x'~*x') ~* (x'~*y''' ~- y'~*x''') ~-+   3 *~ (x'~*y'' ~- y'~*x'') ~* (y'~*y'' ~+ x'~*x'')++-- | Find extrema of the curvature, but not inflection points.+curvatureExtrema :: CubicBezier -> Double -> [Double]+curvatureExtrema b eps = bezierFindRoot (extrema b) 0 1 $+                         bezierParamTolerance b eps++radiusSquareEq :: CubicBezier -> Double -> BernsteinPoly+radiusSquareEq (CubicBezier p0 p1 p2 p3) d =+  let bez = [p0, p1, p2, p3]+      x' = bernsteinDeriv $ listToBernstein $ map pointX bez+      y' = bernsteinDeriv $ listToBernstein $ map pointY bez+      x'' = bernsteinDeriv x'+      y'' = bernsteinDeriv y'+      a =  x'~*x' ~+  y'~*y'+      b =  x'~*y'' ~-  x''~*y'+  in (a~*a~*a) ~- (d*d) *~ b~*b++-- | Find points on the curve that have a certain radius of curvature.+findRadius :: CubicBezier  -- ^ the curve+           -> Double       -- ^ distance+           -> Double       -- ^ tolerance+           -> [Double]     -- ^ result+findRadius b d eps = bezierFindRoot (radiusSquareEq b d) 0 1 $+                     bezierParamTolerance b eps
+ Geom2D/CubicBezier/Intersection.hs view
@@ -0,0 +1,190 @@+{-# LANGUAGE BangPatterns #-}+-- | 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)+       where+import Geom2D+import Geom2D.CubicBezier.Basic+import Math.BernsteinPoly+import Data.Maybe+++-- find the convex hull by comparing the angles of the vectors with+-- the cross product and backtracking if necessary.+findOuter' upper !dir !p1 l@(p2:rest)+  -- backtrack if the direction is outward+  | if upper+    then dir `vectorCross` (p2^-^p1) > 0 -- left turn+    else dir `vectorCross` (p2^-^p1) < 0 = Left l+  -- succeed+  | otherwise = case findOuter' upper (p2^-^p1) p2 rest of+    Left m -> findOuter' upper dir p1 m+    Right m -> Right (p1:m)++findOuter' _ _ p1 p = Right (p1:p)++-- find the outermost point.  It doesn't look at the x values.+findOuter upper (p1:p2:rest) =+  case findOuter' upper (p2^-^p1) p2 rest of+    Right l -> p1:l+    Left l -> findOuter upper (p1:l)+findOuter _ l = l    ++-- take the y values and turn it in into a convex hull with upper en+-- lower points separated.+makeHull :: [Double] -> ([Point], [Point])+makeHull ds =+  let n      = fromIntegral $ length ds - 1+      points = zipWith Point [i/n | i <- [0..n]] ds+  in (findOuter True points,+      findOuter False points)++-- test if the chords cross the fat line+-- use continuation passing style+testBelow :: Double -> [Point] -> Maybe Double -> Maybe Double+testBelow dmin [] _ = Nothing+testBelow dmin [_] _ = Nothing+testBelow dmin (p:q:rest) cont+  | pointY p >= dmin = cont+  | pointY p > pointY q = Nothing+  | pointY q < dmin = testBelow dmin (q:rest) cont+  | otherwise = Just $ intersectPt dmin p q++testBetween :: Double -> Point -> Maybe Double -> Maybe Double+testBetween dmax (Point x y) cont+  | y <= dmax = Just x+  | otherwise = cont++-- test if the chords cross the line y=dmax somewhere+testAbove :: Double -> [Point] -> Maybe Double+testAbove dmax [] = Nothing+testAbove dmax [_] = Nothing+testAbove dmax (p:q:rest)+  | pointY p < pointY q = Nothing+  | pointY q > dmax = testAbove dmax (q:rest)+  | otherwise = Just $ intersectPt dmax p q++-- find the x value where the line through the two points+-- intersect the line y=d+intersectPt d (Point x1 y1) (Point x2 y2) =+  x1 + (d  - y1) * (x2 - x1) / (y2 - y1)++-- make a hull and test over which interval the+-- curve is garuanteed to lie inside the fat line+chopHull dmin dmax ds = do+  let (upper, lower) = makeHull ds+  left_t <- testBelow dmin upper $+            testBetween dmax (head upper) $+            testAbove dmax lower+  right_t <- testBelow dmin (reverse upper) $+             testBetween dmax (last upper) $+             testAbove dmax (reverse lower)+  Just (left_t, right_t)++bezierClip p@(CubicBezier !p0 !p1 !p2 !p3) q@(CubicBezier !q0 !q1 !q2 !q3)+  tmin tmax umin umax prevClip eps reverse++  -- no intersection+  | isNothing chop_interval = []++  -- not enough reduction, so split the curve in case we have+  -- multiple intersections+  | prevClip > 0.8 && newClip > 0.8 =+    if new_tmax - new_tmin > umax - umin -- split the longest segment+    then 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 reverse) +++         bezierClip q pr umin umax half_t new_tmax newClip eps (not reverse)+    else 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 reverse) +++         bezierClip qr newP half_t umax new_tmin new_tmax newClip eps (not reverse)++  -- within tolerance      +  | max (umax - umin) (new_tmax - new_tmin) < eps =+    if reverse+    then [ (umin + (umax-umin)/2,+            new_tmin + (new_tmax-new_tmin)/2) ]+    else [ (new_tmin + (new_tmax-new_tmin)/2,+            umin + (umax-umin)/2) ]++  -- iterate with the curves reversed.+  | otherwise =+      bezierClip q newP umin umax new_tmin new_tmax newClip eps (not reverse)++  where+    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)++-- | Find the intersections between two Bezier curves within given+-- tolerance, using the Bezier Clip algorithm. Returns the parameters+-- for both curves.+bezierIntersection :: CubicBezier -> CubicBezier -> Double -> [(Double, Double)]+bezierIntersection p q eps = bezierClip p q 0 1 0 1 0 eps' False+  where+    eps' = min (bezierParamTolerance p eps) (bezierParamTolerance q eps)++------------------------ Line intersection -------------------------------------+-- Clipping a line uses a simplified version of the Bezier Clip algorithm,+-- and uses the (thin) line itself instead of the fat line.++-- | Find the zero of a 1D bezier curve of any degree.  Note that this+-- can be used as a bernstein polynomial root solver by converting from+-- the power basis to the bernstein basis.+bezierFindRoot :: BernsteinPoly -- ^ the bernstein coefficients of the polynomial+               -> Double  -- ^ The lower bound of the interval +               -> Double  -- ^ The upper bound of the interval+               -> Double  -- ^ The accuracy+               -> [Double] -- ^ The roots found+bezierFindRoot p tmin tmax eps+  -- no intersection+  | chop_interval == Nothing = []++  -- not enough reduction, so split the curve in case we have+  -- multiple intersections+  | clip > 0.8 =+    let (p1, p2) = bernsteinSplit newP 0.5+        half_t = new_tmin + (new_tmax - new_tmin) / 2+    in bezierFindRoot p1 new_tmin half_t eps +++       bezierFindRoot p2 half_t new_tmax eps++  -- within tolerance+  | new_tmax - new_tmin < eps =+      [new_tmin + (new_tmax-new_tmin)/2]++      -- iterate+  | otherwise =+        bezierFindRoot newP new_tmin new_tmax eps++  where+    chop_interval = chopHull 0 0 (bernsteinCoeffs p)+    Just (chop_tmin, chop_tmax) = chop_interval+    newP = bernsteinSubsegment p chop_tmin chop_tmax+    clip = chop_tmax - chop_tmin+    new_tmin = tmax * chop_tmin + tmin * (1 - chop_tmin)+    new_tmax = tmax * chop_tmax + tmin * (1 - chop_tmax)++-- | Find the intersections of the curve with a line.++-- Apply a transformation to the bezier that maps the line onto the+-- X-axis.  Then we only need to test the Y-values for a zero.+bezierLineIntersections :: CubicBezier -> Line -> 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) = +          fromJust (inverse $ translate p $* rotateVec (q ^-^ p)) $* b
+ Geom2D/CubicBezier/Numeric.hs view
@@ -0,0 +1,27 @@+-- | Some numerical computations used by the cubic bezier functions+module Geom2D.CubicBezier.Numeric where++-- | @quadraticRoot a b c@ find the real roots of the quadratic equation+-- @a x^2 + b x + c = 0@.  It will return one, two or zero roots.+quadraticRoot :: Double -> Double -> Double -> [Double]+quadraticRoot a b c = result where+  d = b*b - 4*a*c+  q = - (b + signum b * sqrt d) / 2+  x1 = q/a+  x2 = c/q+  result | d < 0     = []+         | d == 0    = [x1]+         | otherwise = [x1, x2]++-- | @solveLinear2x2 a b c d e f@ solves the linear equation with two variables (x and y) and two systems:+-- +-- >a x + b y + c = 0+-- >d x + e y + f = 0+-- +-- Returns @Nothing@ if no solution is found.+solveLinear2x2 :: Double -> Double -> Double -> Double -> Double -> Double -> Maybe (Double, Double)+solveLinear2x2 a b c d e f =+  case det of 0 -> Nothing+              _ -> Just ((c * e - b * f) / det, (a * f - c * d)  / det)+  where det = d * b - a * e+
+ Geom2D/CubicBezier/Outline.hs view
@@ -0,0 +1,105 @@+-- | Offsetting bezier curves and stroking curves.++module Geom2D.CubicBezier.Outline+       (bezierOffset, bezierOffsetMax)+       where+import Geom2D+import Geom2D.CubicBezier.Basic+import Geom2D.CubicBezier.Approximate+import Geom2D.CubicBezier.Curvature+import qualified Data.Map as M+import Data.Function+import Data.List++offsetPoint :: Double -> Point -> Point -> Point+offsetPoint dist start tangent =+  start ^+^ (rotate90L $* dist *^ normVector tangent)++bezierOffsetPoint :: CubicBezier -> Double -> Double -> Point+bezierOffsetPoint cb dist t =+  uncurry (offsetPoint dist) $+  evalBezierDeriv cb t++-- Approximate the bezier curve offset by dist.  A positive value+-- means to the left, a negative to the right.+approximateOffset :: CubicBezier -> Double -> Double -> (CubicBezier, Double, Double)+approximateOffset cb@(CubicBezier p1 p2 p3 p4) dist tol =+  approximateCurveWithParams offsetCb points ts tol+  where tan1     = p2 ^-^ p1+        tan2     = p4 ^-^ p3+        offsetCb = CubicBezier+                   (offsetPoint dist p1 tan1)+                   (offsetPoint dist p2 tan1)+                   (offsetPoint dist p3 tan2)+                   (offsetPoint dist p4 tan2)+        points   = map (bezierOffsetPoint cb dist) ts+        ts = [i/16 | i <- [1..15]]++-- subdivide the original curve and approximate the offset until+-- the maximum error is below tolerance+offsetSegment :: Double -> Double -> CubicBezier -> [CubicBezier]+offsetSegment dist tol cb+  | err <= tol = [cb_out]+  | otherwise     = offsetSegment dist tol cb_l +++                    offsetSegment dist tol cb_r +  where+    (cb_out, t, err) = approximateOffset cb dist tol+    (cb_l, cb_r) = splitBezier cb t++data OutlineSegment = OutlineSegment {+  os_t_min :: Double,  -- the least t param of the segment in the original curve+  os_t_err :: Double,  -- the param where the error is maximal+  os_curve :: CubicBezier, -- the segment on the original curve+  os_outline :: CubicBezier } -- the outline of the segment++-- Keep a map from maxError to OutlineSegment 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.+offsetMax :: Double -> Double -> Int ->+             M.Map Double OutlineSegment ->+             [CubicBezier]+offsetMax dist tol n segments+  | n <= 1 = error "minimum segments to offset is 1"+  | (n == 1) || (err < tol) = map os_outline $+                              sortBy (compare `on` os_t_min) $+                              M.elems segments++    -- split the maximum curve in two and add the two segments to the map+  | otherwise = offsetMax dist tol (n-1) $+                M.insert err_l (OutlineSegment t_min t_err_l cb_l outline_l) $+                M.insert err_r (OutlineSegment t_err t_err_r cb_r outline_r) $+                newSegments+  where+    ((err, OutlineSegment t_min t_err curve _), newSegments) = M.deleteFindMax segments+    (cb_l, cb_r) = splitBezier curve t_err+    (outline_l, t_err_l, err_l)  = approximateOffset cb_l dist tol+    (outline_r, t_err_r, err_r)  = approximateOffset cb_r dist tol+    +offsetSegmentMax :: Int -> Double -> Double -> CubicBezier -> [CubicBezier]+offsetSegmentMax n dist tol cb =+  offsetMax dist tol n segments+  where segments              = M.singleton err (OutlineSegment 0 t_err cb outline)+        (outline, t_err, err) = approximateOffset cb dist tol++-- | 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 -- ^ The curve+             -> Double      -- ^ Offset distance.+             -> Double      -- ^ Tolerance.+             -> [CubicBezier]        -- ^ The offset curve+bezierOffset cb dist tol =+  --Path $ map BezierSegment $+  concatMap (offsetSegment dist tol) $+  splitBezierN cb $+  findRadius cb dist tol++-- | Like bezierOffset, but limit the number of subpaths for each+-- smooth subsegment.  The number should not be smaller than one.+bezierOffsetMax :: Int -> CubicBezier -> Double -> Double -> [CubicBezier]+bezierOffsetMax n cb dist tol =+  -- Path $ map BezierSegment $+  concatMap (offsetSegmentMax n dist tol) $+  splitBezierN cb $+  findRadius cb dist tol
+ LICENSE view
@@ -0,0 +1,339 @@+GNU GENERAL PUBLIC LICENSE+                       Version 2, June 1991++ Copyright (C) 1989, 1991 Free Software Foundation, Inc.,+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA+ Everyone is permitted to copy and distribute verbatim copies+ of this license document, but changing it is not allowed.++                            Preamble++  The licenses for most software are designed to take away your+freedom to share and change it.  By contrast, the GNU General Public+License is intended to guarantee your freedom to share and change free+software--to make sure the software is free for all its users.  This+General Public License applies to most of the Free Software+Foundation's software and to any other program whose authors commit to+using it.  (Some other Free Software Foundation software is covered by+the GNU Lesser General Public License instead.)  You can apply it to+your programs, too.++  When we speak of free software, we are referring to freedom, not+price.  Our General Public Licenses are designed to make sure that you+have the freedom to distribute copies of free software (and charge for+this service if you wish), that you receive source code or can get it+if you want it, that you can change the software or use pieces of it+in new free programs; and that you know you can do these things.++  To protect your rights, we need to make restrictions that forbid+anyone to deny you these rights or to ask you to surrender the rights.+These restrictions translate to certain responsibilities for you if you+distribute copies of the software, or if you modify it.++  For example, if you distribute copies of such a program, whether+gratis or for a fee, you must give the recipients all the rights that+you have.  You must make sure that they, too, receive or can get the+source code.  And you must show them these terms so they know their+rights.++  We protect your rights with two steps: (1) copyright the software, and+(2) offer you this license which gives you legal permission to copy,+distribute and/or modify the software.++  Also, for each author's protection and ours, we want to make certain+that everyone understands that there is no warranty for this free+software.  If the software is modified by someone else and passed on, we+want its recipients to know that what they have is not the original, so+that any problems introduced by others will not reflect on the original+authors' reputations.++  Finally, any free program is threatened constantly by software+patents.  We wish to avoid the danger that redistributors of a free+program will individually obtain patent licenses, in effect making the+program proprietary.  To prevent this, we have made it clear that any+patent must be licensed for everyone's free use or not licensed at all.++  The precise terms and conditions for copying, distribution and+modification follow.++                    GNU GENERAL PUBLIC LICENSE+   TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION++  0. This License applies to any program or other work which contains+a notice placed by the copyright holder saying it may be distributed+under the terms of this General Public License.  The "Program", below,+refers to any such program or work, and a "work based on the Program"+means either the Program or any derivative work under copyright law:+that is to say, a work containing the Program or a portion of it,+either verbatim or with modifications and/or translated into another+language.  (Hereinafter, translation is included without limitation in+the term "modification".)  Each licensee is addressed as "you".++Activities other than copying, distribution and modification are not+covered by this License; they are outside its scope.  The act of+running the Program is not restricted, and the output from the Program+is covered only if its contents constitute a work based on the+Program (independent of having been made by running the Program).+Whether that is true depends on what the Program does.++  1. You may copy and distribute verbatim copies of the Program's+source code as you receive it, in any medium, provided that you+conspicuously and appropriately publish on each copy an appropriate+copyright notice and disclaimer of warranty; keep intact all the+notices that refer to this License and to the absence of any warranty;+and give any other recipients of the Program a copy of this License+along with the Program.++You may charge a fee for the physical act of transferring a copy, and+you may at your option offer warranty protection in exchange for a fee.++  2. You may modify your copy or copies of the Program or any portion+of it, thus forming a work based on the Program, and copy and+distribute such modifications or work under the terms of Section 1+above, provided that you also meet all of these conditions:++    a) You must cause the modified files to carry prominent notices+    stating that you changed the files and the date of any change.++    b) You must cause any work that you distribute or publish, that in+    whole or in part contains or is derived from the Program or any+    part thereof, to be licensed as a whole at no charge to all third+    parties under the terms of this License.++    c) If the modified program normally reads commands interactively+    when run, you must cause it, when started running for such+    interactive use in the most ordinary way, to print or display an+    announcement including an appropriate copyright notice and a+    notice that there is no warranty (or else, saying that you provide+    a warranty) and that users may redistribute the program under+    these conditions, and telling the user how to view a copy of this+    License.  (Exception: if the Program itself is interactive but+    does not normally print such an announcement, your work based on+    the Program is not required to print an announcement.)++These requirements apply to the modified work as a whole.  If+identifiable sections of that work are not derived from the Program,+and can be reasonably considered independent and separate works in+themselves, then this License, and its terms, do not apply to those+sections when you distribute them as separate works.  But when you+distribute the same sections as part of a whole which is a work based+on the Program, the distribution of the whole must be on the terms of+this License, whose permissions for other licensees extend to the+entire whole, and thus to each and every part regardless of who wrote it.++Thus, it is not the intent of this section to claim rights or contest+your rights to work written entirely by you; rather, the intent is to+exercise the right to control the distribution of derivative or+collective works based on the Program.++In addition, mere aggregation of another work not based on the Program+with the Program (or with a work based on the Program) on a volume of+a storage or distribution medium does not bring the other work under+the scope of this License.++  3. You may copy and distribute the Program (or a work based on it,+under Section 2) in object code or executable form under the terms of+Sections 1 and 2 above provided that you also do one of the following:++    a) Accompany it with the complete corresponding machine-readable+    source code, which must be distributed under the terms of Sections+    1 and 2 above on a medium customarily used for software interchange; or,++    b) Accompany it with a written offer, valid for at least three+    years, to give any third party, for a charge no more than your+    cost of physically performing source distribution, a complete+    machine-readable copy of the corresponding source code, to be+    distributed under the terms of Sections 1 and 2 above on a medium+    customarily used for software interchange; or,++    c) Accompany it with the information you received as to the offer+    to distribute corresponding source code.  (This alternative is+    allowed only for noncommercial distribution and only if you+    received the program in object code or executable form with such+    an offer, in accord with Subsection b above.)++The source code for a work means the preferred form of the work for+making modifications to it.  For an executable work, complete source+code means all the source code for all modules it contains, plus any+associated interface definition files, plus the scripts used to+control compilation and installation of the executable.  However, as a+special exception, the source code distributed need not include+anything that is normally distributed (in either source or binary+form) with the major components (compiler, kernel, and so on) of the+operating system on which the executable runs, unless that component+itself accompanies the executable.++If distribution of executable or object code is made by offering+access to copy from a designated place, then offering equivalent+access to copy the source code from the same place counts as+distribution of the source code, even though third parties are not+compelled to copy the source along with the object code.++  4. You may not copy, modify, sublicense, or distribute the Program+except as expressly provided under this License.  Any attempt+otherwise to copy, modify, sublicense or distribute the Program is+void, and will automatically terminate your rights under this License.+However, parties who have received copies, or rights, from you under+this License will not have their licenses terminated so long as such+parties remain in full compliance.++  5. You are not required to accept this License, since you have not+signed it.  However, nothing else grants you permission to modify or+distribute the Program or its derivative works.  These actions are+prohibited by law if you do not accept this License.  Therefore, by+modifying or distributing the Program (or any work based on the+Program), you indicate your acceptance of this License to do so, and+all its terms and conditions for copying, distributing or modifying+the Program or works based on it.++  6. Each time you redistribute the Program (or any work based on the+Program), the recipient automatically receives a license from the+original licensor to copy, distribute or modify the Program subject to+these terms and conditions.  You may not impose any further+restrictions on the recipients' exercise of the rights granted herein.+You are not responsible for enforcing compliance by third parties to+this License.++  7. If, as a consequence of a court judgment or allegation of patent+infringement or for any other reason (not limited to patent issues),+conditions are imposed on you (whether by court order, agreement or+otherwise) that contradict the conditions of this License, they do not+excuse you from the conditions of this License.  If you cannot+distribute so as to satisfy simultaneously your obligations under this+License and any other pertinent obligations, then as a consequence you+may not distribute the Program at all.  For example, if a patent+license would not permit royalty-free redistribution of the Program by+all those who receive copies directly or indirectly through you, then+the only way you could satisfy both it and this License would be to+refrain entirely from distribution of the Program.++If any portion of this section is held invalid or unenforceable under+any particular circumstance, the balance of the section is intended to+apply and the section as a whole is intended to apply in other+circumstances.++It is not the purpose of this section to induce you to infringe any+patents or other property right claims or to contest validity of any+such claims; this section has the sole purpose of protecting the+integrity of the free software distribution system, which is+implemented by public license practices.  Many people have made+generous contributions to the wide range of software distributed+through that system in reliance on consistent application of that+system; it is up to the author/donor to decide if he or she is willing+to distribute software through any other system and a licensee cannot+impose that choice.++This section is intended to make thoroughly clear what is believed to+be a consequence of the rest of this License.++  8. If the distribution and/or use of the Program is restricted in+certain countries either by patents or by copyrighted interfaces, the+original copyright holder who places the Program under this License+may add an explicit geographical distribution limitation excluding+those countries, so that distribution is permitted only in or among+countries not thus excluded.  In such case, this License incorporates+the limitation as if written in the body of this License.++  9. The Free Software Foundation may publish revised and/or new versions+of the General Public License from time to time.  Such new versions will+be similar in spirit to the present version, but may differ in detail to+address new problems or concerns.++Each version is given a distinguishing version number.  If the Program+specifies a version number of this License which applies to it and "any+later version", you have the option of following the terms and conditions+either of that version or of any later version published by the Free+Software Foundation.  If the Program does not specify a version number of+this License, you may choose any version ever published by the Free Software+Foundation.++  10. If you wish to incorporate parts of the Program into other free+programs whose distribution conditions are different, write to the author+to ask for permission.  For software which is copyrighted by the Free+Software Foundation, write to the Free Software Foundation; we sometimes+make exceptions for this.  Our decision will be guided by the two goals+of preserving the free status of all derivatives of our free software and+of promoting the sharing and reuse of software generally.++                            NO WARRANTY++  11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES+PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,+REPAIR OR CORRECTION.++  12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE+POSSIBILITY OF SUCH DAMAGES.++                     END OF TERMS AND CONDITIONS++            How to Apply These Terms to Your New Programs++  If you develop a new program, and you want it to be of the greatest+possible use to the public, the best way to achieve this is to make it+free software which everyone can redistribute and change under these terms.++  To do so, attach the following notices to the program.  It is safest+to attach them to the start of each source file to most effectively+convey the exclusion of warranty; and each file should have at least+the "copyright" line and a pointer to where the full notice is found.++    Haskell library for manipulating cubic bezier curves+    Copyright (C) 2013  Kristof Bastiaensen++    This program is free software; you can redistribute it and/or modify+    it under the terms of the GNU General Public License as published by+    the Free Software Foundation; either version 2 of the License, or+    (at your option) any later version.++    This program is distributed in the hope that it will be useful,+    but WITHOUT ANY WARRANTY; without even the implied warranty of+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the+    GNU General Public License for more details.++    You should have received a copy of the GNU General Public License along+    with this program; if not, write to the Free Software Foundation, Inc.,+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.++Also add information on how to contact you by electronic and paper mail.++If the program is interactive, make it output a short notice like this+when it starts in an interactive mode:++    Gnomovision version 69, Copyright (C) year name of author+    Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.+    This is free software, and you are welcome to redistribute it+    under certain conditions; type `show c' for details.++The hypothetical commands `show w' and `show c' should show the appropriate+parts of the General Public License.  Of course, the commands you use may+be called something other than `show w' and `show c'; they could even be+mouse-clicks or menu items--whatever suits your program.++You should also get your employer (if you work as a programmer) or your+school, if any, to sign a "copyright disclaimer" for the program, if+necessary.  Here is a sample; alter the names:++  Yoyodyne, Inc., hereby disclaims all copyright interest in the program+  `Gnomovision' (which makes passes at compilers) written by James Hacker.++  {signature of Ty Coon}, 1 April 1989+  Ty Coon, President of Vice++This General Public License does not permit incorporating your program into+proprietary programs.  If your program is a subroutine library, you may+consider it more useful to permit linking proprietary applications with the+library.  If this is what you want to do, use the GNU Lesser General+Public License instead of this License.
+ Math/BernsteinPoly.hs view
@@ -0,0 +1,130 @@+module Math.BernsteinPoly+       (BernsteinPoly(..), bernsteinSubsegment, listToBernstein, zeroPoly, (~*), (*~), (~+),+        (~-), degreeElevate, bernsteinSplit, bernsteinEval,+        bernsteinEvalDerivs, bernsteinDeriv)+       where++import Data.List++data BernsteinPoly = BernsteinPoly {+  bernsteinDegree :: Int,+  bernsteinCoeffs :: [Double] }+                   deriving Show++infixl 7 ~*, *~+infixl 6 ~+, ~-++-- | Create a bernstein polynomail from a list of coëfficients.+listToBernstein :: [Double] -> BernsteinPoly+listToBernstein [] = zeroPoly+listToBernstein l = BernsteinPoly (length l - 1) l++-- | The constant zero.+zeroPoly :: BernsteinPoly+zeroPoly = BernsteinPoly 0 [0]++-- | Return the subsegment between the two parameters.+bernsteinSubsegment :: BernsteinPoly -> Double -> Double -> BernsteinPoly+bernsteinSubsegment b t1 t2 +  | t1 > t2   = bernsteinSubsegment b t2 t1+  | otherwise = snd $ flip bernsteinSplit (t1/t2) $+                fst $ bernsteinSplit b t2++-- multiply two bezier curves+-- control point i from the product of beziers P * Q+-- is sum (P_j * Q_k) where j + k = i+1++-- | Multiply two bernstein polynomials.  The final degree+-- will be the sum of either degrees.  This operation takes O((n+m)^2)+-- with n and m the degree of the beziers.++(~*) :: BernsteinPoly -> BernsteinPoly -> BernsteinPoly+(BernsteinPoly la a) ~* (BernsteinPoly lb b) =+  BernsteinPoly (la+lb) $+  zipWith (flip (/)) (binCoeff (la + lb)) $+                 init $ map sum $+                 zipWith (zipWith (*)) (repeat a') (down b') +++                 zipWith (zipWith (*)) (tail $ tails a') (repeat $ reverse b')+  where down l = tail $ scanl (flip (:)) [] l -- [[1], [2, 1], [3, 2, 1], ...+        a' = zipWith (*) a (binCoeff la)+        b' = zipWith (*) b (binCoeff lb)++degreeElevate' :: BernsteinPoly -> Int -> BernsteinPoly+degreeElevate' b 0 = b+degreeElevate' (BernsteinPoly lp p) times =+  degreeElevate' (BernsteinPoly (lp+1) (head p:inner p 1)) (times-1)+  where+    inner [a] _ = [a]+    inner (a:b:rest) i =+      (i*a/fromIntegral lp + b*(1 - i/fromIntegral lp))+      : inner (b:rest) (i+1)++-- find the binomial coefficients of degree n.+binCoeff :: Int -> [Double]+binCoeff n = map fromIntegral $+             scanl (\x m -> x * (n-m+1) `quot` m) 1 [1..n]++-- | Degree elevate a bernstein polynomail.+degreeElevate :: BernsteinPoly -> Int -> BernsteinPoly+degreeElevate l times = degreeElevate' l times++-- | Evaluate the bernstein polynomial.+bernsteinEval :: BernsteinPoly -> Double -> Double+bernsteinEval (BernsteinPoly lp p) t = foldl' addcoeff 0 $+                                       zip3 ts (binCoeff lp) p+  where ts = iterate (*t) 1+        u = 1-t+        addcoeff a (s, d, b) = (a*u + b*s*d)++-- | Evaluate the bernstein polynomial and its derivatives.+bernsteinEvalDerivs :: BernsteinPoly -> Double -> [Double]+bernsteinEvalDerivs b t+  | bernsteinDegree b == 0 = [bernsteinEval b t]+  | otherwise = bernsteinEval b t :+                bernsteinEvalDerivs (bernsteinDeriv b) t++-- | Find the derivative of a bernstein polynomial.+bernsteinDeriv :: BernsteinPoly -> BernsteinPoly+bernsteinDeriv (BernsteinPoly 0 _) = zeroPoly+bernsteinDeriv (BernsteinPoly lp p) =+  BernsteinPoly (lp-1) $+  map (* fromIntegral lp) $+  zipWith (-) (tail p) p++-- | Split a bernstein polynomial+bernsteinSplit :: BernsteinPoly -> Double -> (BernsteinPoly, BernsteinPoly)+bernsteinSplit (BernsteinPoly lp p) t =+  (BernsteinPoly lp $ map head controls,+   BernsteinPoly lp $ reverse $ map last controls)+  where+    interp a b = (1-t)*a + t*b+    terp [_] = []+    terp l = let ctrs = zipWith interp l (tail l)+             in ctrs : terp ctrs+    controls = p:terp p++-- | Sum two bernstein polynomials.  The final degree will be the maximum of either+-- degrees.+(~+) :: BernsteinPoly -> BernsteinPoly -> BernsteinPoly+ba@(BernsteinPoly la a) ~+ bb@(BernsteinPoly lb b)+  | la < lb = BernsteinPoly lb $+              zipWith (+) (bernsteinCoeffs $ degreeElevate ba $ lb-la) b+  | la > lb = BernsteinPoly la $+              zipWith (+) a (bernsteinCoeffs $ degreeElevate bb $ la-lb)+  | otherwise = BernsteinPoly la $+                zipWith (+) a b++-- | Subtract two bernstein polynomials.  The final degree will be the maximum of either+-- degrees.+(~-) :: BernsteinPoly -> BernsteinPoly -> BernsteinPoly+ba@(BernsteinPoly la a) ~- bb@(BernsteinPoly lb b)+  | la < lb = BernsteinPoly lb $+              zipWith (-) (bernsteinCoeffs $ degreeElevate ba (lb-la)) b+  | la > lb = BernsteinPoly la $+              zipWith (-) a (bernsteinCoeffs $ degreeElevate bb (la-lb))+  | otherwise = BernsteinPoly la $+                zipWith (-) a b++-- | Scale a bernstein polynomial by a constant.+(*~) :: Double -> BernsteinPoly -> BernsteinPoly+a *~ (BernsteinPoly lb b) = BernsteinPoly lb (map (*a) b)
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ cubicbezier.cabal view
@@ -0,0 +1,40 @@+Name:		cubicbezier+Version: 	0.1.0+Synopsis:	Efficient manipulating of 2D cubic bezier curves.+Category: 	Graphics, Geometry, Typography+Copyright: 	Kristof Bastiaensen (2013)+Stability:	Unstable+License:	GPL-2+License-file:	LICENSE+Author:		Kristof Bastiaensen+Maintainer:	Kristof Bastiaensen+Bug-Reports: 	https://github.com/kuribas/cubicbezier/issues+Build-type:	Simple+Cabal-version:	>=1.6+Description:	This library supports efficient manipulating of 2D cubic bezier curves.  The original goal+  is to support typography, but it is useful for general graphics.  Supported features are:+  .+  Evaluating bezier curves and derivatives, affine transformations on bezier curves, arclength and inverse arclength, intersections between two curves, intersection between a curve and a line, curvature and radius of curvature, finding tangents parallel to a vector, finding inflection points and cusps.+  .+  It also supports polynomial root finding with Bernstein polynomials.+  .+  The module "Geom2D.CubicBezier" exports all the cubic bezier functions.  The module "Geom2D"+  contains general 2D geometry functions and transformations.+ +source-repository head+  type:		git+  location:	https://github.com/kuribas/cubicbezier++Library+  Build-depends: base >= 3 && < 5, containers > 0.4, integration >= 0.1.1+  Exposed-Modules:+    Geom2D+    Geom2D.CubicBezier+    Geom2D.CubicBezier.Basic+    Geom2D.CubicBezier.Approximate+    Geom2D.CubicBezier.Outline+    Geom2D.CubicBezier.Curvature+    Geom2D.CubicBezier.Intersection+    Math.BernsteinPoly+  Other-Modules:+    Geom2D.CubicBezier.Numeric