Path operations

[raphlinus]

I'm doing deep research into robust path operations. That may well end up being three blog posts:

• Robust path operations on polylines (Draft of path ops blog #78 is draft)
• Cubic/cubic Bézier intersection
• Putting it all together

These blog posts should be rich in figures. In addition, I am tempted to make interactive diagrams where the control points of the Béziers can be moved, and code updates the various bounds and so on.

Perhaps the order of the last two should be reversed. A key concept is: what should the output of curve/curve intersection look like, to be usable in a path ops context? I propose:

• The two input curves are y-monotonic and have endpoints match in y coordinate (y0..y1)
  • thus, curves can be represented as $x = f(y)$ and $x = g(y)$, $y_0 \leq y \leq y1$
• Output is a sequence of segments that cover the y0..y1 range
• Segment is left-order, ambiguous, or right-order
  • Ambiguous also contains an optional crossing point, only valid for intersections of odd multiplicity
• Order means g(y) - f(y) > epsilon (or g(y) - f(y) < -epsilon)
• Ambiguous means Fréchet distance < epsilon
  • May require some refinement when slope is really high (Fréchet distance at endpoints may be overestimate)
  • Also note these two epsilons are not the same, this one should be higher

How ambiguous sections are handled depends on the path op. For intersection or union, when two curves are near the path op selects one of the curves. Cut over at the crossing point (or don't cut if there's no crossing point). For difference, both curve segments are in the output; merge them (summing winding numbers), possibly adding horizontal segments of length ~ epsilon at boundaries of region.

Remainder of this outline sketches the algorithm for cubic/cubic intersection. It should be much better than the standard "fat line" Bézier clipping algorithm.

Outer loop is subdivision, but works in 1 dimension (y). (Note: not obvious for pedagogical reasons whether to do $y = f(x)$ or $x = f(y)$. Note Bentley-Ottmann paper is the former)

For both cubics, compute an estimated parabola of the form $x = ay^2 + by + c$. One strategy: fix error = 0 at endpoints and make integral of error from y0 to y1 = 0 (signed area). Another strategy: least squares using moments. Probably only a modest constant factor between the two.

Given cubic and parabola estimate, compute bounds so $d_{min} \leq f(y) - (ay^2 + by + c) \leq d_{max}$. Rigorous technique: hybrid curve technique. Given cubic (p0, p1, p2, p3), make two quadratics (p0, lerp(p0, p1, 1.5), p3) and (p0, lerp(p3, p2, 1.5), p3), cubic is enclosed between them. Bounds of those quadratics wrt parabola can be determined by cubic equation. Maybe there's a better constant factor but this solution is good enough because scaling is right.

Claim: cubic scaling. As $\Delta y$ scales, $|f(y) - (ay^2 + by + c)|$ scales by $\Delta y^3$.

Given two estimated parabolas (a, b, c, d_min, d_max), easy to represent upper and lower bounds of difference: $d_{min0} - d_{max1} + (a_1 - a_0)y^2 + (b_1 - b_0)y + c_1 - c_0 \leq g(x) - f(x) \leq d_{max1} - d_{min0} + (a_1 - a_0)y^2 + (b_1 - b_0)y + c_1 - c_0$. Possible ranges where this could be 0 can be determined by two quadratic solves.

Happy case: these ranges are significantly smaller than original y0..y1 range. Keep iterating. Convergence should be cubic.

Also happy case: bounds may prove regions where $|g(x) - f(x) &lt; \epsilon|$. Output as ambiguous.

Less happy case: subdivide, for example at y = 0.5(y0+y1). Given cubic scaling, expect errors to be 1/8.

Big question: when does subdivision stop? Should be $\epsilon^{-1/3}$. Roughly 100 segments for $\epsilon = 10^{-6}$. Not terrible, not great. Can we do better? Adversarial case is when curves are very close to each other. Note: fat line algorithm also does very badly in this case.

Fréchet bound techniques

This section is a bit more speculative. It should be applied when bounds on parabolas do not significantly reduce range.

First, skew f(y) and g(y) linearly $f'(y) = f(y) + a_0y + b_0$ and $g'(y) = g(y) + a_1y + b_1$ so $f'(y_0) = f'(y_1) = g'(y_0) = g'(y_1) = 0$.

Compute upper bound on Fréchet distance between (skewed) f and g = d. Also compute max slope dx/dy of f and g (straightforward to get upper bound using interval techniques) = m. Claim: max $|g'(y) - f'(y)| \leq d \sqrt{1+m^2} = \Delta x$. Now we know $(a_0 - a_1)y + b_0 - b_1 - \Delta x \leq g(y) - f(x) \leq (a_0 - a_1)y + b_0 - b_1 + \Delta x$. Again these bounds can be used to refine results. Hypothesis is that this handles exactly the case that's otherwise tough: curves very close to each other, near second-order intersection.

How to compute upper bound on Fréchet distance? Hypothesis: sample both cubics along 6 points equally spaced by arc length. For some $c$ (hypothesis is that $c = 2$ is sufficient$, Fréchet distance is less than c times maximum distance between pair of sampled points. Holds up in random testing, will try to prove.

References

• Curve intersection using Bézier clipping, Sederberg and Nishita 1990
• Intersection algorithms based on geometric intervals, North 2007, Masters thesis
• Comparison of Distance Measures for Planar Curves, Alt et al, 2004
• Bézier clipping is quadratically convergent, Schultz 2007