The Bezier Clipping Algorithm


Bezier clipping or fatline clipping is a fast and efficient algorithm to find intersections between two bezier curves. Paper.js contains a javascript implementation of the bezier clipping algorithm based on the paper published by Sederberg T., Nishita T. ( The implementation currently assumes cubic bezier curves.

This post attempts to explain the algorithm and its run time characteristics and behaviour. Especially what our implementation of the algorithm does and what it does not.


If you want to run this notebook locally; I would highly recommend it.

Along with all the usual stuff —IPython, Numpy, matplotlib, you would need the following.

Copy into the same folder as this Ipython notebook. (STATUTORY WARNING: code contains nonidiomatic, javascript-like python).

The file is a direct translation of paper.js implementation of bezier clipping into python from javascript.This javascript $\leftrightarrow$ python change should not affect the behaviour of the algorithm discussed below in any significant way. The running time, memory used etc. might change however. The file also contains some not-so-interesting matplotlib code used to draw stuff.

Import the things we need.

In [1]:
%matplotlib inline
import numpy as np
from numpy.random import random

import matplotlib
import matplotlib.pyplot as plt
from matplotlib.path import Path
import matplotlib.patches as patches

from IPython.html.widgets import interact, interactive, fixed
from IPython.html import widgets
from IPython.display import display
In [2]:
%load_ext autoreload
%autoreload 1
%aimport bezclip_manual

# copy into the same folder as this Ipython notebook
import bezclip_manual as bezutil

Some configuration options.

In [3]:
# NOTE: Comment the following line if you are viewing this on a low DPI screen
%config InlineBackend.figure_format = 'retina'
matplotlib.rcParams['']='Source Code Pro'

1 Overview of the bezier clipping algorithm [TOP]

For the purpose of explanation we have to make some basic assumptions. Also I am assuming you are familiar with how bezier curves are defined and their basic properties.

  1. When we talk about curves, it means specifically a $3^{rd}$ degree bezier curve (cubic bezier curve) expressed in parametric form.
  2. A curve is defined by 4 control points or control vertices $(P_0, P_1, P_2, P_3)$
  3. Then the explicit mathematical equation for curve $B$ $$ \begin{equation} \tag{1} {\bf B(t)} = \sum_{i=0}^{3} \underbrace{{3 \choose i}}_{\text{binomial c.}} (1-t)^{(3-i)}.t^i.\underbrace{P_i}_{\text{CVs}} \end{equation}\ ~~~~~~|\ t = [0, 1] $$
  4. $t$ is called the curve time parameter. For the curve above, $t$ varies between 0 and 1 inclusive.
  5. The baseline of a curve, denoted by $\bf L$ is just the line connecting the start and end points of the curve.

1.1 Intersections of bezier curves [TOP]

We can largely group all possible intersections between two bezier curves A and B, into the following.

  1. Crossings – The intersections between A and B are well spaced w.r.t. the error limit ($\epsilon$) chosen in the intersection finding algorithm.
  2. Tangential intersections – Either A and B come very close to ($<\epsilon$) each other and then diverge, or A and B has a pair of intersections which are very close to each other ($\simeq\epsilon$)

The original bezier clipping algorithm, as detailed by Sederberg et. al. handles only crossings, though they extend the algorithm by introducing various concepts such as Focus of a curve. In this notebook, an intersection means a crossing, we are not going to deal with tangential intersections now.

1.2 Fatline of a bezier curve [TOP]

Fatline is just a wide line that completely contains the curve, it is represented by two lines parallel to the curve's baseline $\bf L$ of the curve, at distances $\bf d_{min}$ and $\bf d_{max}$ on either side of $L$.

In [4]:
# v is the array containing <x and y> of 4 coefficients of a cubic bezier curve.
# v = [x0, y0, x1, y1, x2, y2, x3, y3]
def getFatline(v):
    "This method returns a tuple containing maximum and minimum offset width for the 'fatline'."
    # Starting point of the curve
    q0x = v[0] 
    q0y = v[1]
    # End point of the curve
    q3x = v[6] 
    q3y = v[7]
    # Calculate the fat-line L, for Q is the baseline l and two
    # offsets which completely encloses the curve P.
    d1 = bezutil.getSignedDistance(q0x, q0y, q3x, q3y, v[2], v[3]) or 0
    d2 = bezutil.getSignedDistance(q0x, q0y, q3x, q3y, v[4], v[5]) or 0
    factor =  3. / 4. if d1 * d2 > 0 else 4. / 9. # Get a tighter fit
    dMin = factor * min(0, d1, d2)
    dMax = factor * max(0, d1, d2)
    # The width of the 'fatline' is |dMin| + |dMax|
    return (dMin, dMax)

Using the above method we can draw a fatline around a bezier curve.

In [5]:
curve = np.array([-50., 0., -25., 20., 25., -20., 50., 0.])
def drawFatlineI(x2=curve[2], y2=curve[3], x3=curve[4], y3=curve[5]):
    fig = plt.figure(figsize=(6, 6))
    ax = fig.add_axes([0., 0., 1., 1.])
    ax.set_title('The fatline of a curve');

    crv = np.array([curve[0], curve[1], x2, y2, x3, y3, curve[6], curve[7]])
    # Get the fatline
    fatDist = getFatline(crv)
    dmin, dmax = bezutil.drawCurve(ax, crv, '#e74c3c', fatcolor='#9b59b6', fatdist=fatDist, setsize=True)[1]
    # Annotate dMin, dMax etc.
    bp = curve[6:]
    dminE = dmin[2:]
    dmaxE = dmax[2:]
    ax.annotate("", xy=bp, xytext=dminE, arrowprops=dict(arrowstyle="<->", shrinkA=0))
    mid = (bp + dminE)/2.
    ax.text(mid[0]+1, mid[1]-1, '$d_{min}$')
    ax.annotate("", xy=bp, xytext=dmaxE, arrowprops=dict(arrowstyle="<->", shrinkA=0))
    mid = (bp + dmaxE)/2.
    ax.text(mid[0]+1, mid[1]-1, '$d_{max}$')
step = 5
         x2 = widgets.FloatSliderWidget(min=-100.0, max=100.0, step=step, value=curve[2], description="x2"),
         y2 = widgets.FloatSliderWidget(min=-100.0, max=100.0, step=step, value=curve[3], description="y2"),
         x3 = widgets.FloatSliderWidget(min=-100.0, max=100.0, step=step, value=curve[4], description="x3"),
         y3 = widgets.FloatSliderWidget(min=-100.0, max=100.0, step=step, value=curve[5], description="y3"));

1.3 Finding intersections using fatline [TOP]

Let's try and overlap two fatlines. Obviously, the intersection point will be contained within the intersection of both fatlines.

In [6]:
curve1 = np.array([-50., 0., -25., 20., 25., -20., 50., 0.])
curve2 = np.array([10., -40.,  -5., -25.,  -30., 25.,  0., 50.])

curve2r = bezutil.rotate(curve2, rotation)
fig = plt.figure(figsize=(8, 6))
ax = fig.add_axes([0.0, 0.0, 0.75, 1.])
fatDist1 = getFatline(curve1)
fatDist2 = getFatline(curve2r)
bezutil.drawCurve(ax, curve1, '#e74c3c', fatdist=fatDist1, setsize=True)
bezutil.drawCurve(ax, curve2r, '#34495e', fatdist=fatDist2)
locs = bezutil.getCrossings(curve1, curve2r)
for t1, p1, t2, p2 in locs:
    ax.plot([p1.x], [p1.y], 'x', markeredgecolor='w', markeredgewidth=2, markersize=5.0)
# ax.legend(["Curve A", "Curve B"], borderaxespad=0)
count = len(locs)
print str.format('{} crossing{} found.', count, '' if count == 1 else 's')
1 crossing found.

Let's say we pick Curve A and discard the parts of curve that is not within Curve B's fatline. Now we are left with a much smaller Curve A'. We use the fatline of this new curve to clip Curve B. We can narrow down on the intersection point by alternatively clipping the remnants of Curve A with the remnants of Curve B. We can stop once we have clipped both curves to nearly a small point.

NOTE: It might take a bit more time to render the following piece of code, since it renders all steps of the clipping process

In [19]:
curve1 = np.array([-50., 0., -25., 20., 25., -20., 50., 0.])
curve2 = np.array([10., -40.,  -5., -25.,  -30., 25.,  0., 50.])
fatDist1 = getFatline(curve1)
fatDist2 = getFatline(curve2)
nrow, ncol = (6, 3)
gs = plt.GridSpec(nrow, ncol)

@interact(rotation=(0, 2*np.pi), steps=(0, 32, 1))
def plotClipI(rotation=0, steps=7):
    nrow, ncol = (6, 3)
    curve2r = bezutil.rotate(curve2, rotation)
    fig = plt.figure(figsize=(10, 18))
    axes = [plt.subplot(gs[i, j]) for i in range(nrow) for j in range(ncol)]
    for ax in axes:
    locs = bezutil.getCrossings(curve1, curve2r, steps=steps, fig=fig, subplots=axes)
    print "* 'Length' mentioned in the plots is the bezier time parameter t (in the range [0, 1])"

print "Note: This might take some time."
* 'Length' mentioned in the plots is the bezier time parameter t (in the range [0, 1])
Note: This might take some time.

The plots above shows how the algorithm recursively converges to the intersections.

The algorithm always converges to the intersections between two curves, if there are any; also note that the algorithm manages to narrow down on the actual intersection point quite rapidly, and the speed at which it narrows it down accelerates as we are closer to the actual intersection point.

In other words, the algorithm is globally convergent; and on average bezier clipping algorithm converges quadratically. We will examine these properties of the algorithm later in detail.

1.4 The clipping process in detail [TOP]

If we look closely at one of the clipping steps, the clipping is not really exact. In the plot below, Curve B is being clipped by the fatline of Curve A, and parts of curve B extends quite a bit from both sides of the fatline. If we clip enough curves, we can determine that the algorithm always underestimates the possible clip. This has one advantage that we always leave some room for error, and we will not miss the intersection point because of loss of precision doing arithmetic on small IEEE floating point numbers.

In [8]:
curve1 = np.array([-50., 0., -25., 20., 25., -20., 50., 0.])
curve2 = np.array([10., -40.,  -5., -25.,  -30., 25.,  0., 50.])
fatDist1 = getFatline(curve1)
fig = plt.figure(figsize=(5, 5))
ax = fig.add_axes([0., 0., 1., 1.])
locs = bezutil.getCrossings(curve1, curve2, steps=0, size=30, fig=fig, subplots=[ax])
ax.set_title('Detail of fatline clipping');

Let's consider the plot above, Curve $A$'s fatline is clipping Curve $B$. Fatline is a much simpler representation than the actual bezier curve. It is just two parallel lines at $(d_{min},\ d_{max})$ from the baseline ${\bf L}~of~curve~{\bf A}$. And we need a simpler representation for curve $\bf B$ as well.

For $B$, we would apply a simple transformation into something called a non-parametric bezier curve.

In [9]:
# Function to transform a curve into a non-parametric curve
# The non-parametric version of the curve has:
# 1. it's control vertices evenly spaced in [0, 1] along the x axis
# 2. The distance of the control vertices from a given fatline along the yaxis
# this method also return the associated details of the non-parametric curve such as,
# the convex-hull of the control points, and the corresponding tMin and tMax
def nonParamCurveDetails(curve, L, fatdist):
    # Distances of control points from the fatline base L
    dp0 = bezutil.getSignedDistance(L[0], L[1], L[2], L[3], curve[0], curve[1])
    dp1 = bezutil.getSignedDistance(L[0], L[1], L[2], L[3], curve[2], curve[3])
    dp2 = bezutil.getSignedDistance(L[0], L[1], L[2], L[3], curve[4], curve[5])
    dp3 = bezutil.getSignedDistance(L[0], L[1], L[2], L[3], curve[6], curve[7])
    # note that the parameter values 't' are equally spaced in [0, 1]
    # Calculate the convexhull for the non-parametric bezier curve.
    # The algorithm used is an optimized (since we don't have to take care of
    # general cases)version of Andrew's Monotone chain algorithm
    top, bottom = bezutil.getConvexHull(dp0, dp1, dp2, dp3)
    # Now we convex-hull at dMin and dMax
    dMin, dMax = fatdist
    tMinClip = bezutil.clipConvexHull(top, bottom, dMin, dMax);
    tMaxClip = bezutil.clipConvexHull(top, bottom, dMin, dMax);
    # done
    # merge both convexhulls into a polygon for plotting
#     top.append(top[0])
    return (np.array([0., dp0, 1./3., dp1, 2./3., dp2, 1., dp3]),
            np.array(top).reshape(-1), (tMinClip, tMaxClip,),)

# setup the plot
curve1 = np.array([-50., 0., -25., 30., 25., -20., 50., 0.])
curve2 = np.array([10., -40.,  -5., -30.,  -30., 35.,  0., 50.])
gs = plt.GridSpec(1, 2)

@interact(rotation=(0.0, 2*np.pi))
def drawNonParamCurveI(rotation = 0):
    fig = plt.figure(figsize=(12, 6))
    # Calculate the non-parametric curve for A
    curve1r = bezutil.rotate(curve1, rotation)
    fatDist = getFatline(curve1r)
    L = curve1r[np.array([0, 1, 6, 7])]
    nonParamB, chull, tclips = nonParamCurveDetails(curve2, L, fatdist=fatDist)
    # call the visualisation proxy. Not included here, 
    # just because it is quite long and has more matplotlib stuff than bezier clipping!
    bezutil.drawNonParamCurveI(fig, gs, curve1r, curve2, fatDist, nonParamB, chull, tclips)
    tMin, tMax = tclips
    print str.format('Part of curve B, where ({:3.2f} < t < {:3.2f}) are within the fatline.', tMin, tMax);
    print str.format('tDiff after clipping = {:3.2f}', tMax - tMin);
Part of curve B, where (0.28 < t < 0.69) are within the fatline.
tDiff after clipping = 0.42

The plot on the right shows the non-parametric version of Curve $B$ (the dark blue curve on the left plot). It shows the distance of curve $\bf B$ from the fatline base $\bf L$ vs. the bezier time parameter $\bf t_i$. Note that the control points of the non-parametric curve are placed at equally spaced intervals on t, where t = [0, 1]. The non-parametric curve is just a mapping between the distance from a given line ($L$) and the curve time parameter of a bezier curve – for derivations see the paper mentioned in the introduction.

Now we can just clip this distance curve or non-parametric curve at $d_{min}$ and $d_{max}$* and look at the corresponding parameter $t_i$ on the x-axis of the second plot and keep parts of $B$ at $(t_{min} < t < t_{max})$.

But there is a small issue with this idea. The non-parametric curve is still a bezier curve, and it is expensive and difficult to calculate the intersection between a line and the non-parametric curve. If we recursively execute such an expensive operation, the algorithm might turn out to be the slowest in it's class.

We can use the convex-hull property of a bezier curve to make clipping much simpler and faster. We would just clip away parts of convex-hull of the non-parametric curve that lies above the horizontal line $y=d_{min}$ and the parts below $y=d_{max})$. Then we can choose $t_{min}$ and $t_{max}$ as the minimum and maximum values along the x-axis ($t_i$) corresponding to the clipped convex-hull. Convex-hull only contains straight lines and this is much easier to deal with.

A suitably low value $\epsilon$ is chosen, so that if $t_{diff}$ on both curves is less than that value, the algorithm stops clipping and returns the average of $t_{min}$ and $t_{max}$ as the intersection between A and B. In paperjs $\epsilon$ is defined as $1 \times 10 ^{-5}$.

This method can be optimised further to speed it up, such as generating convex-hulls faster since we already have some information about the nature of the non-parametric curve (control points equally spaced about the x-axis, only 4 points in case of cubic curves etc.).

In the plot above, you can rotate the fatline and notice how the algorithm behaves when the angle between the fatline and the curve changes.

*Also you may have noticed the relationship between the clipped bezier overhangs from the fatline edges (on the left plot) and the distance from the clipping points on the convex-hull (red markers on the right plot) and the nearest point on the non-parametric curve along $d_{min}$ or $d_{max}$.

1.5 Multiple intersections [TOP]

Like many subdivision algorithms, bezier clipping converges to one intersection at a time. If there are multiple intersections, the fatline would fail to clip away the curve effectively.

The way we handle this is to check if there is no significant reduction in $(t_{diff} = t_{max} - t_{min})$ between successive clipping operation, subdivide one of the curves in half and proceed with the algorithm. Sederberg et. al. recommends a reduction by at least 20% ($\Delta_t < 0.8$).

The value 0.8 was selected based on testing and various performance tradeoffs. Any value reasonably within the range (0.1, 0.9) may work with the algorithm, however if the value is too small, say <= 0.5, the algorithm will eagerly subdivide cases where it doesn't need to; or if this value is too high (> 0.85), the algorithm will try to recurse more without much success until, $\Delta t$ becomes more than this value and then subdivide.

1.6 The complete algorithm [TOP]

method clipBezier(A :curve, B: curve):

  1. Compute the fatline for ${\bf A}$ $(L, d_{min}, d_{max})$
  2. Compute the distance $d_i$ of control points of ${\bf B}$ from the fatline base
  3. Create a non-parametric curve from the distances $d_i$ obtained above
  4. Clip the convex-hull of the non-parametric curve, calculate $t_{min}$ and $t_{max}$
    1. If either the entire convex-hull is above $d_{max}$ or below $d_{min}$, the curves do not intersect
  5. Let ${\bf B^\prime} \leftarrow B(t_{min}, t_{max})$ ; part of $B$ roughly contained by the fatline
  6. If $(t_{max}-t_{min}) < errorEpsilon$:
    1. add $(t_{max}-t_{min})/2$ to the list of intersections
    2. Return
  7. Else if $\Delta t_{current} > 0.8\ and\ \Delta t_{previous} > 0.8$:
    1. $parts \leftarrow subdivide(longest\ of\ (A, B^\prime))$; let $X$ be the shortest of $A, B^\prime$
    2. clipBezier($parts[0]$, $X$)
    3. clipBezier($parts[1]$, $X$)
  8. Else:
    1. clipBezier($B^\prime$, $A$)

2 The run-time behaviour of the bezier clipping algorithm [TOP]

2.1 Generating test data [TOP]

A set of around 600 test cases are included in file, since generating this test data may take some time (especially for number of intersections > 4).

The test data has the following distribution.

In [10]:
fig, ax = plt.subplots(1)
x = np.arange(7, 0, -1) - 0.1, bezutil.NTestCases, 0.2, color='#f67088', edgecolor='#f67088');
ax.set_xlabel('Number of intersections')
ax.set_ylabel('Number of test cases')
ax.set_title('Test data distribution');

Please note that this distribution of test data is not representative of what the algorithm might encounter in real world, where I believe, a significant majority of the cases have 1 or 2 intersections between them.

One caveat of using randomly generated test cases is that, the curvature variations within the curve is often not wild enough to produce curve pairs intersecting at more than 6 or 7 locations. In order to asses the behaviour of the algorithm for different inputs, included below are some hand-rolled test cases, with number of intersections ranging from $1~to~9$. As you may remember, the maximum possible number of intersections between two cubic bezier curves is 9.

In [11]:
cases1_9 = [[[ 40.1, -23.1, -5.399999999999999, -26, 7.200000000000003, 52.1, -32.7, 28.6 ],
            [ -40.4, -34.8, -6.899999999999999, -31.3, -1, 38, 29.5, 41.5 ]],
            [[ -3.6, -42.6, -35.6, 2.200000000000003, -43.99999999999999, 37, 44.6, 11.1 ],
            [ -42.2, -13.1, 54.10000000000001, -49.2, 48, -15.899999999999999, 6, 40.5 ]],
            [[ -41.6, -22.2, 22, -50.2, -34.6, 92.7, 32.9, 19.8 ],
            [ -23.1, -39.9, -35.2, 27.9, 42.4, 33.00000000000001, 42.4, 33.00000000000001 ]],
            [[ -44.9, -41.8, 13.100000000000001, 36.10000000000001, 99.29999999999998, 75.3, 1.3000000000000043, -41.8 ],
            [ -44.9, -10.6, -3.1000000000000014, -21.3, 69.9, -57.8, 23.7, 42.7 ]],
            [[ -40.9, -15.6, 77, -85.8, -75.7, 84, 41.7, 18.5 ],
            [ 9.4, 0.4, -77, -100.1, -16.6, 11.200000000000003, 41.5, 43 ]],
            [[ 42.4, -29.6, -202.6, 25.4, 189.3, -24.700000000000003, -37.300000000000004, 35.49999999999999 ],
            [ -30.8, -26.3, -16.2, 40.900000000000006, -2.1000000000000014, 90.19999999999999, 25.7, -40.5 ]],
            [[ 24.2, 42.4, -162.60000000000002, -162.5, 125.4, 89.5, -31.000000000000004, -30.6 ],
            [ -39.8, -18.7, 75.5, -37.4, -58.1, 69.4, 42.4, 29.9 ]],
            [[ -41.9, 7, 160.9, -155.1, -172.4, 155, 38.00000000000001, -11.2 ],
            [ -31.5, 8.2, 184.1, 102.8, -160, -133.5, 16.4, 20 ]],
            [[ 25.3, 21.4, -93.4, -180.5, 90.9, 177.2, -31, -15.8 ],
            [ 26.9, -22.6, -196.3, 48.300000000000004, 193.4, -52, -21.8, 24 ]]];

# plot the test cases
gs = plt.GridSpec(3, 3)
fig = plt.figure(figsize=(6, 6))
axes = [plt.subplot(gs[i, j]) for i in range(3) for j in range(3)]
colors = ['#e74c3c', '#34495e']
for i, ax in enumerate(axes):
    bezutil.drawCurve(ax, cases1_9[i][0], color=colors[0], points=False)
    bezutil.drawCurve(ax, cases1_9[i][1], color=colors[1], points=False)
    ax.set_xlim(left=-50, right=50)
    ax.set_ylim(bottom=-50, top=50)
    locs = bezutil.getCrossingsPlain(cases1_9[i][0], cases1_9[i][1])
    n = len(locs)
    ax.set_title(str.format("{} crossing{}", n, 's' if n>1 else ''))
    for loc in locs:
        ax.plot(loc[1].x, loc[1].y, marker='o', color='w', markeredgecolor='k')

2.2 Rate of convergence [TOP]

To compare the performance, let's look at the amortised rate of convergence of the algorithm.

The bezier clipping algorithm converges to $L$ with a rate of convergence $\bf q$ if, there exists a number $\mu ~|~ \mu \in (0, 1)$ such that,

$$ \begin{equation} \tag{2} \lim_{x\rightarrow\infty} {\frac {|x_{k+1}-L|}{|x_k-L|^q}} = \mu \end{equation} $$

Here we use $t_{diff}$, which should ideally converge to 0 at an intersection, so $L = 0$.

In [12]:
results = map(lambda c: bezutil.getCrossingsInstr(c[0], c[1]), bezutil.TestCases)

Reduce the results and plot them.

In [13]:
def reduceTestResults(tDiffs):
    """Returns (nIntersections, maxDepth, nSpans, tReduction) for a test result"""
    maxDepth = max(tDiffs)[0]
    tDiffSum = [0. for x in range(maxDepth+1)]
    nIx = nDegen = 0
    for t in tDiffs:
        nIx += 1 if t[4] else 0;
        nDegen += 1 if t[3] else 0;
        tDiffSum[t[0]] += t[5]
    tRedux = []
    for i in range(2, len(tDiffSum)):
        A = np.ceil(np.log10(1./ (tDiffSum[i-2] if tDiffSum[i-2] > 0 else 1.)))
        B = np.ceil(np.log10(1./ (tDiffSum[i] if tDiffSum[i] > 0 else 1.)))
        redux = B / A
        if redux == np.inf or redux == -np.inf:
            redux = 0
    return np.concatenate(([nIx, maxDepth, (nIx+nDegen)], tRedux))

# NOTE: Results are arranged num(IX) from 7 to 1
V = map(lambda o: reduceTestResults(o[1]), results)
splitLen = np.cumsum(bezutil.NTestCases)[:-1]
V = np.split(V, splitLen)
# Group result arrays by number of intersections and normalize the shape
dataForIx = []
for v in V:
    padz = np.array(map(len, v))
    maxD = max(padz)
    padz = maxD - padz
    for i, pz in enumerate(padz):
        v[i] = np.pad(v[i], (0, pz), 'constant', constant_values=[np.nan])
    dataForIx.append(np.concatenate(v).reshape((-1, maxD)))

resultsSum = []
for tIx in dataForIx:
    algDetail ,tOrder = np.split(tIx, [3], axis = 1)
    # The summary
    algDetail = np.median(algDetail, axis=0)
    tOrder = np.nanmean(tOrder, axis=0)
    resultsSum.append((algDetail, tOrder))

palette = ['#f67088', '#c59331', '#82a831', '#34af89', '#37aaba', '#8096f4', '#f45deb'];

def plotResult(ax, res, degree=3, num=50, color=0, plotOriginal=True, showTable=False, label='', mlabel=''):
    xv = np.arange(0, len(res))
    val = np.linspace(0, len(res), num=num)
    p = np.poly1d(np.polyfit(xv, res, degree))
    if plotOriginal:
        ax.plot(res, '.', color=palette[color]);
    ax.plot(val, p(val), lw=2, color=palette[color], label=label);
    mean = np.mean(res)
    ax.axhline(mean, color=palette[color], label=mlabel)
    # raw data as a table
    if showTable:
        table = bezutil.ListTable()
        tD = [[d, '{:2.1}'.format(t)] for d, t in enumerate(res)]
        tD = ['$Depth$'];
        tD.extend(range(1, len(res)+1))
        tD = ['$Rate of convergence$']
        tD.extend(['{:2.1f}'.format(t) for d, t in enumerate(res)])
def iPlot(n_Intersections = 1):
    n_Intersections -= 1
    fig, ax = plt.subplots(1)
    res = resultsSum[n_Intersections][1]
    mean = np.mean(res)
    plotResult(ax, res, color=n_Intersections, mlabel='Mean = {:2.2f}'.format(mean), showTable=True)
    ax.set_xlabel('$Recursion\ depth$', fontdict={'size': 15});
    ax.set_ylabel('$Rate\ of\ convergence\ (q)$', fontdict={'size': 15});
    ax.set_title('Rate of convergence vs recursion depth');
$Rate of convergence$

Try changing the plot for number of intersections between 1 and 7. The average rate of convergence per depth is around 2 is all cases from depth > ~10. Below is a plot of all test cases overlaid.

In [14]:
fig, ax = plt.subplots(1)
ax.set_xlabel('$Recursion\ depth$', fontdict={'size': 15});
ax.set_ylabel('$Rate\ of\ convergence\ (q)$', fontdict={'size': 15});
ax.set_title('Rate of convergence vs recursion depth');
means = []
for i, res in enumerate(resultsSum):
    plotResult(ax, res[1], color=i, plotOriginal=False, label='${}\ intersections$'.format(i+1))
ax.legend(loc=2, borderaxespad=0., bbox_to_anchor=(1.01, 1));
print 'The average rate of convergence is {:2.1f}'.format(np.mean(means))
The average rate of convergence is 1.6