{-# LANGUAGE MultiWayIf, PatternGuards, TemplateHaskell, BangPatterns #-}

Removing overlap from bezier paths in haskell

This document describes an algorithm for removing overlap and performing set operations on bezier paths, but at the same time it is a working module for the haskell cubicbezier package. This way it can serve two purposes at once: someone who wants to implement this algorithm can use this as an explanation of the algorithm, while at the same time it is a working version of the described algorithm.

Note on porting: Porting this code to another language should present no difficulties. However some care must be taken with regards to lazyness. Often many variables inside the where statement aren't evaluated in all guards, so it's important to evaluate only those which appear in the guards. The main state can be modified using mutation instead of copying without any troubles. For modifying the state, the lens library is used. The lens functions can be interpreted in a mutable language as follows:

Let's begin with declaring the module and library imports:

module Geom2D.CubicBezier.Overlap
       (boolPathOp, union, intersection, difference,
        exclusion, FillRule (..))
import Prelude hiding (mapM)
import Geom2D
import Geom2D.CubicBezier.Basic
import Geom2D.CubicBezier.Intersection
import Math.BernsteinPoly
import Data.Traversable (mapM)
import Data.Functor ((<$>))
import Data.List (sortBy, sort)
import Control.Monad.State hiding (mapM)
import Lens.Micro
import Lens.Micro.TH
import Lens.Micro.Mtl
import qualified Data.Map.Strict as M
import qualified Data.Set as S

So what does it mean to remove overlap? Basicly we want to keep curves where one side is inside the filled region, and the other side is outside, and discard the rest.Since that could be true only of a part of the curve, we also need to split each curve when it intersects another curve. How do you know which side is the inside, and which side the outside? There are two methods which are use the most: the even-odd rule and the nonzero rule. Instead of hardwiring it, I use higher-order functions to determine when a turnratio is inside the region to be filled, and how the turnratio changes with each curve.

Checking each pair of curves for intersections would work, but is rather inefficient. We only need to check for overlap when two curves are adjacent. Fortunately there exist a good method from computational geometry, called a sweep line algorithm. The basic idea is to sweep a vertical line over the input, starting from leftmost point to the right (of course the opposite direction is also possible), and to update the input dynamically. We keep track of each curve that intersects the sweepline by using a balanced tree of curves. When adding a new curve, it's only necessary to check for intersections with the curve above and below. Since searching on the tree takes only O(log n) operations, this will save a lot of computation.

The input is processed in horizontal order, and after splitting curves the order must be preserved, so an ordered structure is needed. The standard map library from Data.Map is ideal, and has all the operations needed. This structure is called the X-structure, since the data is ordered by X coordinate.:

type XStruct = M.Map PointEvent [Curve]

Why PointEvent, and not just Point? We need to have a Ord instance for the map, which much match our horizontal ordering. A newtype is ideal, since it has no extra cost, and allows us to define a Ord instance for defining the relative order. The value from the map is a list, since there can be many curves starting from the same point.

newtype PointEvent = PointEvent DPoint
                   deriving Show

When the x-coordinates are equal, I use the y-coordinate to determine the order.

instance Eq PointEvent where
  (PointEvent (Point x1 y1)) == (PointEvent (Point x2 y2)) =
    (x1, y1) == (x2, y2)

instance Ord PointEvent where
  compare (PointEvent (Point x1 y1)) (PointEvent (Point x2 y2)) =
    compare (x1, y2) (x2, y1)

All curves are kept left to right, so we need to remember the direction for the output:

The curves intersecting the sweepline are kept in another balanced Tree, called the Y-structure. These curves are not allowed to overlap, except in the endpoints, and will be ordered vertically. I'll use the Curve datatype to define the ordering of the curves, and to add additional information. The turnRatio field is the turnRatio of the area to the left for a left to right curve, and to the right for a right to left curve. The changeTurn function determines how the turnRatio will change from up to down. This together with a test for the insideness of a certain turnratio, allows for more flexibility. Using this, it is possible to generalize this algorithm to boolean operations!

The FillRule datatype is used for the exported API:

data FillRule = EvenOdd | NonZero
data Curve = Curve {
  _bezier :: !(CubicBezier Double),
  _turnRatio :: !(Int, Int),
  _changeTurn :: !((Int, Int) -> (Int, Int))}

trOne :: (Int, Int)
trOne = (0,0)

makeLenses ''Curve

instance Show Curve where
  show (Curve b a _) =
    "Curve " ++ show b ++ " " ++ show a

type YStruct = S.Set Curve

The total state for the algorithm consists of the X-structure, the Y-structure, and the output found so far. I use a trick to make access to curves above and below the current pointevent more convenient. I use two sets to represent a focus point into the Y-structure, where the left set are the elements less than the pointEvent (above), and the right set the elements greater (below):

data SweepState = SweepState {
  _output :: !(M.Map PointEvent [CubicBezier Double]),
  _yStructLeft :: !YStruct,
  _yStructRight :: !YStruct,
  _xStruct :: !XStruct}
                  deriving Show
makeLenses ''SweepState

Changing the focus point can be done efficiently in O(log n) by mering and splitting again:

changeFocus :: DPoint -> SweepState -> SweepState
changeFocus p sweep =
  let (lStr, rStr) =
        S.split (Curve (CubicBezier p p p p) trOne id) $
        S.union (view yStructLeft sweep) (view yStructRight sweep)
  in set yStructLeft lStr $
     set yStructRight rStr

This handy helper function will pass the first curve above to the given function, and if it doesn't return Nothing, remove it from the state. It does nothing when there is no curve above.

withAbove :: (Curve -> Maybe a) -> State SweepState (Maybe a)
withAbove f = do
  lStr <- use yStructLeft
  if S.null lStr
    then return Nothing
    else let (c, lStr') = S.deleteFindMax lStr
         in case f c of
             Nothing ->
               return Nothing
             Just x -> do
               yStructLeft .= lStr'
               return $ Just x

The same with the curve below.

withBelow :: (Curve -> Maybe a) -> State SweepState (Maybe a)
withBelow f = do
  rStr <- use yStructRight
  if S.null rStr
    then return Nothing
    else let (c, rStr') = S.deleteFindMin rStr
         in case f c of
             Nothing ->
               return Nothing
             Just x -> do
               yStructRight .= rStr'
               return $ Just x

splitYStruct changes the focus and returns and removes any curves which end in the current pointEvent:

splitYStruct :: DPoint -> State SweepState [Curve]
splitYStruct p = do
  modify $ changeFocus p
  let go = do
        mbC <- withAbove $ \c ->
          -- remove and return c if it ends in point p
          guard (cubicC3 (_bezier c) == p) >> Just c
        case mbC of
         Just c ->
           (c:) <$> go
         Nothing -> return []

Some functions on the Sweep state:

Adding and removing curves from the X structure.

insertX :: PointEvent -> [Curve] -> SweepState -> SweepState
insertX p c =
  over xStruct $ M.insertWith (++) p c

xStructAdd :: Curve -> SweepState -> SweepState
xStructAdd c =
  insertX (PointEvent $ cubicC0 $
                        view bezier c) [c]

xStructRemove :: State SweepState (PointEvent, [Curve])
xStructRemove = zoom xStruct $ state M.deleteFindMin

To compare curves vertically, take the the curve which starts the rightmost, and see if it falls below or above the curve. If the first control points are coincident, test the last control points instead. The curves in the Y-structure shouldn't intersect (except in the endpoints), so these cases don't have to be handled. To lookup a single point, I use a singular bezier curve.

instance Eq Curve where
  Curve c1 t1 ct1 == Curve c2 t2 ct2 =
    c1 == c2 && t1 == t2 && ct1 (ct2 t1) == t1
instance Ord Curve where
  compare (Curve c1@(CubicBezier p0 p1 p2 p3) tr1 _)
    (Curve c2@(CubicBezier q0 q1 q2 q3) tr2 _)
    | p0 == q0 = if
        | p3 == q3 ->
            -- compare the midpoint
            case (compVert (evalBezier c1 0.5) c2) of
             LT -> LT
             GT -> GT
             EQ ->
               -- otherwise arbitrary
               compare (tr1, PointEvent p1, PointEvent p2)
               (tr2, PointEvent q1, PointEvent q2)
        | pointX p3 < pointX q3 ->
            case (compVert p3 c2) of
            LT -> LT
            EQ -> LT
            GT -> GT
        | otherwise ->
            case compVert q3 c1 of
             LT -> GT
             EQ -> GT
             GT -> LT
    | pointX p0 < pointX q0 =
      case compVert q0 c1 of
       LT -> GT
       EQ -> LT
       GT -> LT
    | otherwise =
      case (compVert p0 c2) of
      LT -> LT
      EQ -> GT
      GT -> GT

Compare a point with a curve. See if it falls below or above the hull first. Otherwise find the point on the curve with the same X-coordinate by solving a cubic equation.

compVert :: DPoint -> CubicBezier Double -> Ordering
compVert p c
  | p == cubicC0 c ||
    p == cubicC3 c = EQ
  | compH /= EQ = compH
  | otherwise = comparePointCurve p c
      compH = compareHull p c

Test if the point is above or below the curve

comparePointCurve :: Point Double -> CubicBezier Double -> Ordering
comparePointCurve (Point x1 y1) c1@(CubicBezier p0 p1 p2 p3)
  | pointX p0 == x1 &&
    pointX p0 == pointX p1 &&
    pointX p0 == pointX p2 &&
    pointX p0 == pointX p3 =
    compare (pointY p0) y1
  | otherwise = compare y2 y1
    t = findX x1 c1 $
        maximum (map maxp [p0, p1, p2, p3])*1e-12
    maxp (Point x y) = max (abs x) (abs y)
    y2 = pointY $ evalBezier c1 t

Comparing against the hull

Compare a point against the convex hull of the bezier. EQ means the point is inside the hull, LT below and GT above. I am currently only testing against the control points, some testing needs to be done to see what is faster.

belowLine :: DPoint -> DPoint -> DPoint -> Bool
belowLine (Point px py) (Point lx ly) (Point rx ry)
  | lx == rx = True
  | (px >= lx && px <= rx) ||
    (px <= lx && px >= rx) = py < midY
  | otherwise = True
  where midY = ly + (ry-ly) * (rx-lx) / (px-lx)

aboveLine :: DPoint -> DPoint -> DPoint -> Bool
aboveLine (Point px py) (Point lx ly) (Point rx ry)
  | lx == rx = True
  | (px >= lx && px <= rx) ||
    (px <= lx && px >= rx) = py > midY
  | otherwise = True
  where midY = ly + (ry-ly) * (rx-lx) / (px-lx)

compareHull :: DPoint -> CubicBezier Double -> Ordering
compareHull p (CubicBezier c0 c1 c2 c3)
  | pointY p > pointY c0 &&
    pointY p > pointY c1 &&
    pointY p > pointY c2 &&
    pointY p > pointY c3 = LT
  | pointY p < pointY c0 &&
    pointY p < pointY c1 &&
    pointY p < pointY c2 &&
    pointY p < pointY c3 = GT
  | otherwise = EQ


Since the algorithm assumes curves are increasing in the horizontal direction they have to be preprocessed first. I split each curve where the tangent is vertical. If the resulting subsegment is too small however, I just adjust the control point to make the curve vertical at the endpoint.

I also do snaprounding to prevent points closer than the tolerance.

makeXStruct :: ((Int, Int) -> (Int, Int)) -> ((Int, Int) -> (Int, Int)) -> Double -> [CubicBezier Double] -> XStruct
makeXStruct chTr chTrBack tol =
  M.fromListWith (++) .
  concatMap (toCurve . snapRoundBezier tol) .
  concatMap (splitVert tol)
  where toCurve c@(CubicBezier p0 _ _ p3) =
          case compare (pointX p0) (pointX p3) of
           LT -> [(PointEvent p0, [Curve c trOne chTr])]
           GT -> [(PointEvent p3, [Curve (reorient c) trOne chTrBack]),
                  (PointEvent p0, [])]
           -- vertical curve
           EQ | pointY p0 > pointY p3 ->
                [(PointEvent p0, [Curve c trOne chTr])]
              | otherwise ->
                [(PointEvent p3, [Curve (reorient c) trOne chTrBack]),
                 (PointEvent p0, [])]

splitVert :: Double -> CubicBezier Double -> [CubicBezier Double]
splitVert tol curve@(CubicBezier c0 c1 c2 c3) =
  uncurry splitBezierN $
  adjustLast $
  adjustFirst (curve, vert)
  where vert
          | pointX c0 == pointX c1 &&
            pointX c0 == pointX c2 &&
            pointX c0 == pointX c3 = []
          | otherwise = 
              sort $ bezierVert curve
        -- adjust control points to avoid small curve fragments
        -- near the endpoints
        adjustFirst (c@(CubicBezier p0 p1 p2 p3), t:ts)
          | vectorDistance p0 (evalBezier c t) < tol =
              (CubicBezier p0 (Point (pointX p0) (pointY p1)) p2 p3,
        adjustFirst x = x
        adjustLast (c@(CubicBezier p0 p1 p2 p3), ts@(_:_))
          | vectorDistance p3 (evalBezier c $ last ts) < tol =
              (CubicBezier p0 p1 (Point (pointX p3) (pointY p2)) p3,
               init ts)
        adjustLast x = x

main loop

For the main loop, we remove the leftmost point from the X-structure, and do the following steps:

  1. Split any curves which come near the current pointEvent.

  2. Send all curves to the left of the sweepline to the output, after filtering them based on the turning number.

  3. For each curve starting at the point, split if it intersects with the curve above or the curve below. Sort resulting curves vertically. If there are no curves starting from point, test the curves above and below instead. Adjust the turnRatios for each curve.

  4. Insert the points in the Y structure.

  5. Loop until the X-structure is empty

loopEvents :: ((Int, Int) -> Bool) -> Double -> SweepState -> SweepState
loopEvents isInside tol sweep 
  | M.null $ view xStruct sweep = sweep
  | otherwise =
    loopEvents isInside tol $!
    flip execState sweep $ do
      -- remove leftmost point from X structure
      (PointEvent p, curves) <- xStructRemove
      -- change focus, and remove curves ending at current
      -- pointevent from Y structure
      ending <- splitYStruct p
      -- split near curves
      (ending2, rightSubCurves) <- splitNearPoints p tol
      -- output curves to the left of the sweepline.
      modify $ filterOutput (ending ++ ending2) isInside 
      let allCurves = rightSubCurves ++ curves
      if null allCurves
         -- split surrounding curves
        then splitSurround tol
        else do
        -- sort curves
        sorted <- splitAndOrder tol allCurves
        -- split curve above
        curves2 <- splitAbove sorted tol
        -- add curves to Y structure
        addMidCurves curves2 tol

Send curves to output

outputPaths :: (M.Map PointEvent [CubicBezier Double]) -> [ClosedPath Double]
outputPaths m
  | M.null m = []
  | otherwise = outputNext m
    lookupDelete p m =
      case M.lookup (PointEvent p) m of
       Nothing -> Nothing
       Just (x:xs) -> Just (x, m')
         where m' | null xs = M.delete (PointEvent p) m
                  | otherwise = M.insert (PointEvent p) xs m
       _ -> error "outputPaths: empty list inside map."
    outputNext !m
      | M.null m = []
      | otherwise = 
        let ((PointEvent p0, (c0:cs)), m0) =
              M.deleteFindMin m
            m0' | null cs = m0
                | otherwise = M.insert (PointEvent p0) cs m0
        in go m0' c0 [] p0
    go !m !next !prev !start
      | p == start =
          curvesToPath (reverse $ next:prev):
          outputNext m
      | otherwise =
        case lookupDelete p m of
         Nothing -> outputNext m
         Just (x, m') -> go m' x (next:prev) start
      where p = cubicC3 next

curvesToPath :: [CubicBezier Double] -> ClosedPath Double
curvesToPath =
  ClosedPath .
  map (\(CubicBezier p0 p1 p2 _) ->
        (p0, JoinCurve p1 p2))

Filter and output the given curves. The isInside function determines the insideness of a give turnratio. For example for the nonzero-rule, this would be (> 0). This inserts the curve into the output map.

filterOutput :: [Curve] -> ((Int, Int) -> Bool) -> SweepState -> SweepState
filterOutput curves isInside sweep =
  foldl (flip $ outputCurve isInside) sweep curves

outputCurve :: ((Int, Int) -> Bool) -> Curve -> SweepState -> SweepState
outputCurve isInside (Curve c tr op)
  | isInside (op tr) /= isInside tr =
      let c' | isInside tr = reorient c
             | otherwise = c
      in over output (M.insertWith (++) (PointEvent $ cubicC0 c') [c'])
  | otherwise = id

Test for intersections and split:

Since the curves going out of the current pointEvent in the X-structure are unordered, they need to be ordered first. First they are ordered by first derivative. Since it's easier to compare two curves when they don't overlap, remove overlap, and then sort again by comparing the whole curve.

To do this, I implemented a monadic insertion sort. First the curves are split in the statemonad, then they are compared.

splitAndOrder :: Double -> [Curve] -> State SweepState [Curve]
splitAndOrder tol curves =
  sortSplit tol $
  sortBy compDeriv curves

compDeriv :: Curve -> Curve -> Ordering
compDeriv (Curve (CubicBezier p0 p1 _ _) _ _)
  (Curve (CubicBezier q0 q1 _ _) _ _) =
  compare (vectorCross (p1^-^p0) (q1^-^ q0)) 0

Insertion sort, by splitting and comparing. This should be efficient enough, since ordering by derivative should mostly order the curves.

sortSplit :: Double -> [Curve] -> State SweepState [Curve]
sortSplit _ [] = return []
sortSplit tol (x:xs) =
  insertM x tol =<<
  sortSplit tol xs

insertM :: Curve -> Double -> [Curve] -> State SweepState [Curve]
insertM x _ [] = return [x]
insertM x tol (y:ys) =
  case curveOverlap x y tol of
   Just (c1, c2) -> do
     mapM (modify . xStructAdd) c2
     insertM c1 tol ys
   Nothing -> do
     (x', y') <- splitM x y tol
     if x' < y'
       then return (x':y':ys)
       else (y':) <$> insertM x' tol ys

splitM :: Curve -> Curve -> Double -> State SweepState (Curve, Curve)
splitM x y tol =
  case splitMaybe x y tol of
  (Just (a, b), Just (c, d)) -> do
    modify $ insertX (PointEvent $ cubicC0 $ view bezier b) [b, d]
    return (a, c)
  (Nothing, Just (c, d)) -> do
    modify $ insertX (PointEvent $ cubicC0 $ view bezier d) [d]
    return (x, c)
  (Just (a, b), Nothing) -> do
    modify $ insertX (PointEvent $ cubicC0 $ view bezier b) [b]
    return (a, y)
  (Nothing, Nothing) ->
    return (x, y)

Handle intersections of the first curve at point and the curve above. Return the curves with updated turnratios. Some care is needed when one of the curves is intersected at the endpoints, in order not to create singular curves.

updateTurnRatio :: Curve -> Curve -> Curve
updateTurnRatio (Curve _ tr chTr) =
  set turnRatio (chTr tr)

propagateTurnRatio :: Curve -> [Curve] -> [Curve]
propagateTurnRatio cAbove l =
  tail $ scanl updateTurnRatio cAbove l

splitAbove :: [Curve] -> Double -> State SweepState [Curve]
splitAbove [] _ = return []
splitAbove (c:cs) tol = do
  lStr <- use yStructLeft
  if S.null lStr
    then let c' = set turnRatio trOne c
         in return $ c':propagateTurnRatio c' cs
    else do
    let (cAbove, lStr') = S.deleteFindMax lStr
    case splitMaybe c cAbove tol of
     (Nothing, Nothing) ->
       return $ propagateTurnRatio cAbove $ c:cs
     (Just (c1, c2), Nothing) ->
       if cubicC3 (_bezier c1) == cubicC0 (_bezier cAbove)
       then do
         modify $ xStructAdd cAbove . xStructAdd c2
         yStructLeft .= lStr'
         return $ propagateTurnRatio cAbove $ c1:cs
       else do
         modify $ xStructAdd c2
         return $ propagateTurnRatio cAbove $ c1:cs
     (Nothing, Just (c3, c4)) ->
       if cubicC3 (_bezier c3) == cubicC0 (_bezier c)
       then error "curve intersecting pointevent"
       else do
         modify $ xStructAdd c4
         yStructLeft .= S.insert c3 lStr'
         return $ propagateTurnRatio cAbove $ c:cs
     (Just (c1, c2), Just (c3, c4)) -> do
       modify $ xStructAdd c2 . xStructAdd c4
       yStructLeft .= S.insert c3 lStr'
       return $ propagateTurnRatio cAbove $ c1:cs

Split curves near the point. Return the curves starting from this point, and the index of the last split point

splitNearPoints :: DPoint -> Double -> State SweepState ([Curve], [Curve])
splitNearPoints p tol = do
  curves1 <- splitNearDir withAbove p tol
  curves2 <- splitNearDir withBelow p tol
  return (map fst curves1 ++ map fst curves2,
          map snd curves1 ++ map snd curves2)

splitNearDir  :: ((Curve -> Maybe (Curve, Double))
                  -> State SweepState (Maybe (Curve, Double)))
              -> DPoint -> Double
              -> State SweepState [(Curve, Curve)]
splitNearDir dir p tol = do
  mbSplit <- dir $ \curve ->
    (,) curve <$>
    pointOnCurve tol p
    (view bezier curve)
  case mbSplit of
   Nothing -> return []
   Just (curve, t) -> do
     let (c1, c2) = splitBezier (view bezier curve) t
         c1' = adjust curve $ adjustC3 p $
               snapRound tol <$> c1
         c2' = adjust curve $ adjustC0 p $
               snapRound tol <$> c2
     ((c1', c2'):) <$> splitNearDir dir p tol

Add the sorted curves starting at point to the Y-structure, and test last curve with curve below.

addMidCurves :: [Curve] -> Double -> State SweepState ()
addMidCurves [] _ = return ()
addMidCurves [c] tol =
  splitBelow c tol
addMidCurves (c:cs) tol = do
  yStructLeft %= S.insert c 
  addMidCurves cs tol
splitBelow :: Curve -> Double -> State SweepState ()
splitBelow c tol = do
  rStr <- use yStructRight
  let (cBelow, rStr') = S.deleteFindMin rStr
  if S.null rStr
    then yStructLeft %= S.insert c
    case splitMaybe c cBelow tol of
     (Nothing, Nothing) ->
       yStructLeft %= S.insert c
     (Nothing, Just (c3, c4)) ->
       if cubicC3 (_bezier c3) == cubicC0 (_bezier c)
       then error "internal error: splitBelow: curve starting in future"
       else do
         modify $ xStructAdd c4
         yStructLeft %= S.insert c . S.insert c3
         yStructRight .= rStr'
     (Just (c1, c2), Nothing) ->
       if cubicC3 (_bezier c1) == cubicC0 (_bezier cBelow)
       then error "internal error: splitBelow: curve intersecting pointevent."
       else do
         modify $ xStructAdd c2
         yStructLeft %= S.insert c1
     (Just (c1, c2), Just (c3, c4)) -> do
       modify $ xStructAdd c2 . xStructAdd c4
       yStructLeft %= S.insert c1 . S.insert c3
       yStructRight .= rStr'

If no curves start from the point, we have to check if the surrounding curves overlap.

splitSurround :: Double -> State SweepState ()
splitSurround tol = do
  lStr <- use yStructLeft
  rStr <- use yStructRight
  if S.null lStr || S.null rStr
    then return ()
    else do
    let (cAbove, lStr') = S.deleteFindMax lStr
        (cBelow, rStr') = S.deleteFindMin rStr
    case splitMaybe cAbove cBelow tol of
     (Just (c1, c2), Just (c3, c4)) -> do
       modify $ xStructAdd c2 .
         xStructAdd c4
       yStructLeft .= S.insert c1 lStr'
       yStructRight .= S.insert c3 rStr'
     (Just (c1, c2), Nothing) -> do
       modify $ xStructAdd c2
       yStructLeft .= S.insert c1 lStr'
     (Nothing, Just (c1, c2)) -> do
       modify $ xStructAdd c2
       yStructRight .= S.insert c1 rStr'
     (Nothing, Nothing) ->
       return ()

Find curve intersections

Test if both curves intersect. Split one or both of the curves when they intersect. Also snapround each point, and make sure the point of overlap is the same in both curves.

splitMaybe :: Curve -> Curve -> Double ->
              (Maybe (Curve, Curve),
               Maybe (Curve, Curve))
splitMaybe c1 c2 tol =
  (adjustSplit c1 <$> fst n,
   adjustSplit c2 <$> snd n)
    n = nextIntersection b1 b2 tol $
        bezierIntersection b1 b2 pTol
    pTol = min (bezierParamTolerance b1 tol)
           (bezierParamTolerance b2 tol)
    b1 = view bezier c1
    b2 = view bezier c2

adjustSplit :: Curve -> (CubicBezier Double, CubicBezier Double) -> (Curve, Curve)
adjustSplit curve (b1, b2)   =
  (set bezier b1 curve,
   set bezier b2 curve)

adjust :: Curve -> CubicBezier Double -> Curve
adjust curve curve2 = set bezier curve2 curve

snapRoundBezier :: Double -> CubicBezier Double -> CubicBezier Double
snapRoundBezier tol = fmap (snapRound tol)

Given a list of intersection parameters, split at the next intersection, but don't split at the first or last control point, or when the two curves are (nearly) coincident. Note that list of intersections is read lazily, in order not to evaluate more intersections that necessary.

nextIntersection :: CubicBezier Double -> CubicBezier Double -> Double -> [(Double, Double)]
                 -> (Maybe (CubicBezier Double, CubicBezier Double),
                     Maybe (CubicBezier Double, CubicBezier Double))
nextIntersection _ _ _ [] = (Nothing, Nothing)
nextIntersection b1@(CubicBezier p0 _ _ p3) b2@(CubicBezier q0 _ _ q3) tol ((t1, t2): ts)
  | atStart1 && atStart2 =
      nextIntersection b1 b2 tol ts
  | atStart1 =
     Just (adjustC3 p0 $ snapRoundBezier tol b2l,
           adjustC0 p0 $ snapRoundBezier tol b2r))
  | atStart2 =
    (Just (adjustC3 q0 $ snapRoundBezier tol b1l,
           adjustC0 q0 $ snapRoundBezier tol b1r),
  | atEnd1 && atEnd2 = (Nothing, Nothing)
  | atEnd1 =
     Just (adjustC3 p3 $ snapRoundBezier tol b2l,
           adjustC0 p3 $ snapRoundBezier tol b2r))
  | atEnd2 =
    (Just (adjustC3 q3 $ snapRoundBezier tol b1l,
           adjustC0 q3 $ snapRoundBezier tol b1r),
  | bezierEqual b1l b2l tol =
      nextIntersection b1 b2 tol ts
  | otherwise =
      let pMid = snapRound tol <$> cubicC3 b1l
      in (Just (snapRoundBezier tol b1l,
                snapRoundBezier tol b1r),
          Just (adjustC3 pMid $ snapRoundBezier tol b2l,
                adjustC0 pMid $ snapRoundBezier tol b2r))
     x1 = evalBezier b1 t1
     x2 = evalBezier b2 t2
     atStart1 = vectorDistance (cubicC0 b1) x1 < tol
     atStart2 = vectorDistance (cubicC0 b2) x2 < tol
     atEnd1 = vectorDistance (cubicC3 b1) x1 < tol
     atEnd2 = vectorDistance (cubicC3 b2) x2 < tol
     (b1l, b1r) = splitBezier b1 t1
     (b2l, b2r) = splitBezier b2 t2

adjustC0 :: Point a -> CubicBezier a -> CubicBezier a
adjustC0 p (CubicBezier _ p1 p2 p3) = CubicBezier p p1 p2 p3

adjustC3 :: Point a -> CubicBezier a -> CubicBezier a
adjustC3 p (CubicBezier p0 p1 p2 _) = CubicBezier p0 p1 p2 p

Check if curves overlap.

If the curves overlap, combine the overlapping part into one curve. To compare the curves, I first split the longest curve so that the velocities in the first control point match, then compare those curves for equality.

curveOverlap :: Curve -> Curve -> Double
             -> Maybe (Curve, Maybe Curve)
curveOverlap c1 c2 tol
  -- starting in the same point
  | p0 /= q0 = Nothing
  | colinear (view bezier c1) tol = if
      | not $ colinear (view bezier c2) tol ->
      | vectorDistance (p3^-^p0)
        ((q3^-^q0) ^* (d1/d2)) > tol ->
      | p3 == q3 -> 
          Just (combineCurves c2 c1,
      | d1 > d2 ->
          Just (combineCurves c2 c1,
                Just $ adjust c1 $
                CubicBezier q3
                (snapRound tol <$> interpolateVector q3 p3 (1/3))
                (snapRound tol <$> interpolateVector q3 p3 (2/3))
      | otherwise ->
          Just (combineCurves c1 c2,
                Just $ adjust c2 $
                CubicBezier p3
                (snapRound tol <$> interpolateVector p3 q3 (1/3))
                (snapRound tol <$> interpolateVector p3 q3 (2/3))
  -- equalize velocities, and compare           
  | v1 == 0 ||
    v2 == 0 = Nothing
  | v1 > v2 = if bezierEqual b2 b1l tol
              then Just (combineCurves c2 c1,
                         if checkEmpty b1r tol
                         then Nothing
                         else Just $ adjust c1 $
                              adjustC0 (cubicC3 b2) $
                              snapRoundBezier tol b1r)
              else Nothing
  | otherwise =
      if bezierEqual b1 b2l tol
              then Just (combineCurves c1 c2,
                         if checkEmpty b2r tol
                         then Nothing
                         else Just $ adjust c2 $
                              adjustC0 (cubicC3 b1) $
                              snapRoundBezier tol b2r)
              else Nothing
    (b1l, b1r) = splitBezier b1 (v2/v1)
    (b2l, b2r) = splitBezier b2 (v1/v2)
    b1@(CubicBezier p0 p1 _ p3) = view bezier c1
    b2@(CubicBezier q0 q1 _ q3) = view bezier c2
    d1 = vectorDistance p0 p3
    d2 = vectorDistance q0 q3
    v1 = vectorDistance p0 p1
    v2 = vectorDistance q0 q1

checkEmpty :: CubicBezier Double -> Double -> Bool
checkEmpty (CubicBezier p0 p1 p2 p3) tol = 
  p0 == p3 &&
  vectorDistance p0 p1 < tol &&
  vectorDistance p0 p2 < tol

Curves can be combined if they are equal, just by composing their changeTurn functions.

combineCurves :: Curve -> Curve -> Curve
combineCurves c1 c2 =
  over changeTurn (view changeTurn c2 .) c1


snapRound :: Double -> Double -> Double
snapRound tol v =
  fromInteger (round (v/tol)) * tol

Test if the point is on the curve (within tolerance)

pointOnCurve :: Double -> DPoint -> CubicBezier Double -> Maybe Double
pointOnCurve tol p c1
  | (t:_) <-
    closest c1 p tol,
    p2 <- evalBezier c1 t,
    vectorDistance p p2 < tol = Just t
  | otherwise = Nothing

Testing beziers for approximate equality

If the control points of two bezier curves are within a distance eps from each other, then both curves will all so be at least within distance eps from each other. This can be proven easily: subtracting both curves gives the distance curve. Since each control point of this curve lies within a circle of radius eps, by the convex hull property, the curve will also be inside the circle, so the distances between each point will never exceed eps.

bezierEqual :: CubicBezier Double -> CubicBezier Double -> Double -> Bool
bezierEqual cb1@(CubicBezier a0 a1 a2 a3) cb2@(CubicBezier b0 b1 b2 b3) tol
  -- controlpoints equal within tol
  | vectorDistance a1 b1 < tol &&
    vectorDistance a2 b2 < tol &&
    vectorDistance a3 b3 < tol &&
    vectorDistance a0 b0 < tol = True
  -- compare if both are colinear and close together
  | dist < tol &&
    colinear cb1 ((tol-dist)/2) &&
    colinear cb2 ((tol-dist)/2) = True
  | otherwise = False
  where dist = max (abs $ ld b0) (abs $ ld b3)
        ld = lineDistance (Line a0 a3)

Finding the on curve point at the X coordinate

Solve a cubic equation to find the X coordinate. This should be converted to a closed form solver for efficiency.

findX :: Double -> CubicBezier Double -> Double -> Double
findX x c@(CubicBezier p0 p1 p2 p3) eps =
  head $ bezierFindRoot bez 0 1 $
  bezierParamTolerance c (eps/10)
  where bez = listToBernstein 
              (map pointX [p0, p1, p2, p3]) ~-
              listToBernstein [x, x, x, x]

Higher level functions

fillFunction :: FillRule -> Int -> Bool
fillFunction NonZero = (>0)
fillFunction EvenOdd = odd
-- | Union of paths, removing overlap and rounding to the given
-- tolerance.
union :: [ClosedPath Double] -- ^ Paths
         -> FillRule         -- ^ input fillrule
         -> Double           -- ^ Tolerance
         -> [ClosedPath Double]
union p fill tol =
  outputPaths out
    out = view output $
          loopEvents (fillFunction fill . fst) tol sweep
    sweep = SweepState M.empty S.empty S.empty xStr
    xStr = makeXStruct (over _1 $ subtract 1) (over _1 (+1)) tol $
           concatMap closedPathCurves p

-- | combine paths using the given boolean operation    
boolPathOp :: (Bool -> Bool -> Bool) -- ^ operation
          -> [ClosedPath Double]     -- ^ first path (merged with union)
          -> [ClosedPath Double]     -- ^ second path (merged with union)
          -> FillRule                -- ^ input fillrule
          -> Double                  -- ^ tolerance 
          -> [ClosedPath Double]
boolPathOp op p1 p2 fill tol =
  outputPaths $ view output $
  loopEvents isInside tol sweep
    isInside (a, b) = fillFunction fill a `op`
                      fillFunction fill b
    sweep = SweepState M.empty S.empty S.empty xStr
    xStr = M.unionWith (++)
            (over _1 (subtract 1))
            (over _1 (+1)) tol $
            concatMap closedPathCurves p1)
            (over _2 (subtract 1))
            (over _2 (+1)) tol $
            concatMap closedPathCurves p2)

intersection, difference, exclusion ::
  [ClosedPath Double] -> [ClosedPath Double] ->
  FillRule -> Double -> [ClosedPath Double]

-- | path intersection  
intersection = boolPathOp (&&)

-- | path difference
difference = boolPathOp (\a b -> a && not b)

-- | path exclusion
exclusion = boolPathOp (\a b -> if a then not b else b)