diff --git a/Geom2D.hs b/Geom2D.hs new file mode 100644 --- /dev/null +++ b/Geom2D.hs @@ -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 + diff --git a/Geom2D/CubicBezier.hs b/Geom2D/CubicBezier.hs new file mode 100644 --- /dev/null +++ b/Geom2D/CubicBezier.hs @@ -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 + + + diff --git a/Geom2D/CubicBezier/Approximate.hs b/Geom2D/CubicBezier/Approximate.hs new file mode 100644 --- /dev/null +++ b/Geom2D/CubicBezier/Approximate.hs @@ -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) diff --git a/Geom2D/CubicBezier/Basic.hs b/Geom2D/CubicBezier/Basic.hs new file mode 100644 --- /dev/null +++ b/Geom2D/CubicBezier/Basic.hs @@ -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) + diff --git a/Geom2D/CubicBezier/Curvature.hs b/Geom2D/CubicBezier/Curvature.hs new file mode 100644 --- /dev/null +++ b/Geom2D/CubicBezier/Curvature.hs @@ -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 diff --git a/Geom2D/CubicBezier/Intersection.hs b/Geom2D/CubicBezier/Intersection.hs new file mode 100644 --- /dev/null +++ b/Geom2D/CubicBezier/Intersection.hs @@ -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 diff --git a/Geom2D/CubicBezier/Numeric.hs b/Geom2D/CubicBezier/Numeric.hs new file mode 100644 --- /dev/null +++ b/Geom2D/CubicBezier/Numeric.hs @@ -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 + diff --git a/Geom2D/CubicBezier/Outline.hs b/Geom2D/CubicBezier/Outline.hs new file mode 100644 --- /dev/null +++ b/Geom2D/CubicBezier/Outline.hs @@ -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 diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/Math/BernsteinPoly.hs b/Math/BernsteinPoly.hs new file mode 100644 --- /dev/null +++ b/Math/BernsteinPoly.hs @@ -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) diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/cubicbezier.cabal b/cubicbezier.cabal new file mode 100644 --- /dev/null +++ b/cubicbezier.cabal @@ -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