Fun with Curves: Splitting and Combining

Let's start with a pair of Bézier splines running roughly between the same starting and ending points. Let's suppose we want to turn all the regions of intersection into filled polygons. Like this:

How might we go about this?

Let's break it down. We need to:

Find the points of intersection
Break the edges at the points of intersection
Find the relevant edges (and subedges) of both curves
Assemble the polygons from the right edges in the right order

Intersection Points

These curves consist of a sequence of quadratic and cubic Bézier curves. To find the intersection points between two edges, we run a recursive divide-and-conquer algorithm:

Find bounding boxes of two edges
If the bounding boxes do not overlap, there is no intersection
If the bounding boxes do overlap, break each edge in two and get intersection points from each subpair (each subpart of edge 1 against each subpart of edge 2).
When the size of the bounding boxes gets down to some threshold level, pick one end of one of the edges to act as the intersection point

In generally this behaves well and converges reasonably, although if you set the threshold too small and the point of intersection is right on a split boundary it can get bad. If you're operating in canvas space, a threshold of 1 is fine.

Bounding Boxes

How to we find bounding boxes in the first place?

I'm using axis-aligned bounding boxes here, because the overlap test is stupid-simple (just check min/max coordinate overlap).

Each curve is defined by a parametric function over t (from 0 to 1). Where the derivative of that function is 0, the curve reaches an extremum. If that occurs where the parameter t is in range (i.e. between 0 and 1) then the point at that value is relevant. If we take the minimum and maximum coordinates at those points as well as at the end points, we have the coordinates of the bounding box.

A quadratic Bézier is defined by the parametric equation B(t) = (1-t)²P0 + 2(1-t)tP1 + t²P2, where P0 is the starting point, P1 is the control point, and P2 is the end point.

Derivatives of Bézier curves of degree N are easy to compute: they are Bézier curves of degree N-1. In particular, the derivative of B(t) above is B'(t) = (1-t)a + tb where a = 2(P0-P1) and b = 2(P2-P1). Solving a linear equation is t is easy: B'(t) = 0 where t = -a/(b - a).

Bear in mind these are all vector equations, we we end up with solutions for x and for y. (Find the point value by substituting the t we solved for into the parametric equation.)

B(t) = (1-t)²P0 + 2(1-t)tP1 + t²P2
B'(t) = (1-t)(2(P0-p1)) + t(2(P0-P1))
t = -2(P0-P1) / (2(P0-P1) - 2(P0-P1))

A cubic Bézier is defined by the parametric equation B(t) = (1-t)³P0 + 3(1-t)²tP1 + 3(1-t)t²P2 + t³P3 where P0 is the starting point, P2 and P3 are control points, and P4 is the end point. The derivative is B'(t) = (1-t)²A + 2(1-t)tB + t²C where A = 3(P1-P0), B = 3(P2-P1), and C = 3(P3-P2). We then have a quadratic equation B'(t) = at² + bt + c = 0 with a = A - 2B + C = 3(-P0 + 3P1 - 3P2 + P3), b = 2(B - A) = 6(P0 - 2P1 + P2), and c = A = 3(P1 - P0). The solution to this is t = (-b ± √(b² - 4ac))/2a.

B(t) = B(t) = (1-t)³P0 + 3(1-t)²tP1 + 3(1-t)t²P2 + t³P3
B'(t) = B'(t) = (1-t)²3(P1-P0) + 2(1-t)t3(P2-P1) + t²3(P3-P2)
t = -6(P0 - 2P1 + P2) ±
    √((6(P0 - 2P1 + P2))² - 4(3(-P0 + 3P1 - 3P2 + P3))(3(P1 - P0))) /
    2(3(-P0 + 3P1 - 3P2 + P3))

Splitting the Curve

We want to split a Bézier quadratic or cubic edge in half, maintaining the overall shape. More generally, when we get to step 2 of the algorithm given at the start, we want to split the curve at an arbitrary point. It is easiest to frame this in turns of the parameter t of the parametric equations for the edges. The procedure involves computing intermediate points that are linear combinations of existing points, weighted by t.

For a quadratic edge, the two subedges are defined by:

P01 = P0 + (P1 - P0)*t
P12 = P1 + (P2 - P0)*t
P012 = P01 + (P12 - P01)*t

The resulting subedges:

quad: start=P0, control=P01, end=P012
quad: start=P012, control=P12, end=P2

For a cubic edge edge, the two subedges are defined by:

P01 = P0 + (P1 - P0)*t
P12 = P1 + (P2 - P1)*t
P23 = P2 + (P3 - P2)*t
P012 = P01 + (P12 - P01)*t
P123 = P12 + (P23 - P12)*t
P0123 = P012 + (P123 - P012)*t

The resulting subedges:

cubic: start=P0, controls=P01 and P012, end=P0123, 
cubic: start=P0123, controls=P123 and P23, end=P3

Certain operations are less likely to run into edge case issues if we try to avoid creating quadratic edges that are really straight lines or cubic edges that are really quadratic, so it is helpful to check the closeness of the various points within some epsilon, and downgrade edges as appropriate. For example, if P123 and P23 are essentially the same, the second part of the cubic split can be a quadratic with just P123 as the control point, and so forth.

Paths, Multiple intersections

So far we have dealt only with edges, but we have complete paths, so a little adjustment is needed. This is easy enough: check edge by edge, and to slice a path at a point, make the first slice of the path consist of the edges before the slice point plus the first part of the sliced edge, and make the second slice of the path consist of the second part of the sliced edge plus all edges after the sliced edge.

A further adjustment is needed when we have multiple intersections of the same edge. Say we have two intersection points and we want to create a three-way slice of the edge. Once we make the first cut, the t value for the second intersection point is no longer valid: we need to adjust it to be relative to the new subedge.

If we have t1 and t2, where t1 is less than t2, to adjust t2 to be appropriate relative to the second part of the slice at t1, then use

t2' = (t2 - t1) / (1 - t1)

One way to think of this is in terms of t as a fraction of length along the edge. This isn't quite right, but works well enough. Call the length of the edge l and lengths of the subedges at the t1 slice l1 and l2. That is, we want t2 * l - l1 = t2' * l2. With l1 = t1 * l and l2 = l - l1 we have t2 * l - t1 * l = t2' * (l - t1 * l). The lengths factor out and solving for t2' gives the equation above.

Apply the adjusted t2' to the sliced edge, and we get the new slice at the appropriate point. The same process can be applied repeatedly to t values for subsequent slices for an arbitrary number of intersection point.

Forming Polygons

The general case of this is very hard, but this is not the general case: we have two curves that more or less start and end in the same place, so the way I do this is to assume that contiguous spans of edges in the curves go together, and partition them (with edge slices) into polygon chunks. This is bad for the general case and does fail for particularly wiggly paths in this situation as well, but usually does pretty well.

I first run through one curve edge by edge looking for intersections with any edges on the other curve. If I find one I record the edges and the t values involved. Then I'll run though the edges of the curve again, collecting formed polygons and subsequent edges with no intersections (these will become edges in the new polygon). The edges from the second path to use to finish the polygon are those following the last one up to the latest path2 edge we see intersecting with this path1 edge. Edges get sliced appropriately, and path2 edges are reversed and used in reverse order so that a loop is formed.

Dividing Space

We know how to chop up a Bézier paths at arbitrary points. We know how to find intersections. So now I can have fun dividing up a Bézier polygon into pieces and filling the pieces with quadrangles created by chopping up the path into points, and finding where lines extended from those points intersect with other paths (as above). The missing piece here is determining the appropriate direction to head: I need to know the normal vector. I don't use it here: I use the normal vector at a particular point plus or minus a random amount to determine the overall angle of the quadrangles for a given region.

The normal vector is 90 degrees from the tangent vector and the tangent vector is just the value of the derivative parametric function at a given t value, normalized to produce a unit vector, so we already know how to do that too!

Circular Arcs

In theory computing intersections and splits on circular arcs is easier than for cubic and quadratic Bézier edges. The problem is that SVG circular arcs are somewhat messy: defined by end points and a radius (with the center of the circle implicit) plus a couple of flags. One flag determines the sweep direction and the other whether this is a large arc (more that 180°). And in fact, SVG circular arcs are really elliptical arcs that happen to have the the major and minor radii equal. The flags determine both which part of the circle (ellipse) is drawn and which of the possible circles intersecting through those two end points is the actual one.

All of this is messy and difficult to work with as is, so to do most things I normalize my arc representations to have the circle they are part of plus the starting and ending angles and normalizing the angle values so that the start angle is less than the ending angle. The normal vector is just the angle from the center. Splitting the arc is just doing basic arithmetic on the arc angles. Getting the point at a particular angle from the center is center + (r cos(θ), r sin(θ)). When this gets rendered out as SVG, the sweep flag is always off and the large arc flag can be determined by looking at the angle values.

References

[1] A Primer on Bézier Curves, filled with truly excellent information about how to do various things with Bézier curves. Whatever you need to do: start here.
[2] Math StackExchange, has a lot of well-answered questions about Bézier curves and geometric problems in general.