diff --git a/Lib/fontTools/misc/arrayTools.py b/Lib/fontTools/misc/arrayTools.py
index e76ced7f8..138ad8f36 100644
--- a/Lib/fontTools/misc/arrayTools.py
+++ b/Lib/fontTools/misc/arrayTools.py
@@ -228,6 +228,19 @@ def rectCenter(rect):
     (xMin, yMin, xMax, yMax) = rect
     return (xMin+xMax)/2, (yMin+yMax)/2
 
+def rectArea(rect):
+    """Determine rectangle area.
+
+    Args:
+        rect: Bounding rectangle, expressed as tuples
+            ``(xMin, yMin, xMax, yMax)``.
+
+    Returns:
+        The area of the rectangle.
+    """
+    (xMin, yMin, xMax, yMax) = rect
+    return (yMax - yMin) * (xMax - xMin)
+
 def intRect(rect):
     """Round a rectangle to integer values.
 
diff --git a/Lib/fontTools/misc/bezierTools.py b/Lib/fontTools/misc/bezierTools.py
index 659de34e2..63bfb0903 100644
--- a/Lib/fontTools/misc/bezierTools.py
+++ b/Lib/fontTools/misc/bezierTools.py
@@ -2,9 +2,13 @@
 """fontTools.misc.bezierTools.py -- tools for working with Bezier path segments.
 """
 
-from fontTools.misc.arrayTools import calcBounds
+from fontTools.misc.arrayTools import calcBounds, sectRect, rectArea
+from fontTools.misc.transform import Offset, Identity
 from fontTools.misc.py23 import *
 import math
+from collections import namedtuple
+
+Intersection = namedtuple("Intersection", ["pt", "t1", "t2"])
 
 
 __all__ = [
@@ -25,6 +29,14 @@ __all__ = [
     "splitCubicAtT",
     "solveQuadratic",
     "solveCubic",
+    "quadraticPointAtT",
+    "cubicPointAtT",
+    "linePointAtT",
+    "segmentPointAtT",
+    "lineLineIntersections",
+    "curveLineIntersections",
+    "curveCurveIntersections",
+    "segmentSegmentIntersections",
 ]
 
 
@@ -42,23 +54,31 @@ def calcCubicArcLength(pt1, pt2, pt3, pt4, tolerance=0.005):
     Returns:
         Arc length value.
     """
-    return calcCubicArcLengthC(complex(*pt1), complex(*pt2), complex(*pt3), complex(*pt4), tolerance)
+    return calcCubicArcLengthC(
+        complex(*pt1), complex(*pt2), complex(*pt3), complex(*pt4), tolerance
+    )
 
 
 def _split_cubic_into_two(p0, p1, p2, p3):
-    mid = (p0 + 3 * (p1 + p2) + p3) * .125
-    deriv3 = (p3 + p2 - p1 - p0) * .125
-    return ((p0, (p0 + p1) * .5, mid - deriv3, mid),
-            (mid, mid + deriv3, (p2 + p3) * .5, p3))
+    mid = (p0 + 3 * (p1 + p2) + p3) * 0.125
+    deriv3 = (p3 + p2 - p1 - p0) * 0.125
+    return (
+        (p0, (p0 + p1) * 0.5, mid - deriv3, mid),
+        (mid, mid + deriv3, (p2 + p3) * 0.5, p3),
+    )
+
 
 def _calcCubicArcLengthCRecurse(mult, p0, p1, p2, p3):
-	arch = abs(p0-p3)
-	box = abs(p0-p1) + abs(p1-p2) + abs(p2-p3)
-	if arch * mult >= box:
-		return (arch + box) * .5
-	else:
-		one,two = _split_cubic_into_two(p0,p1,p2,p3)
-		return _calcCubicArcLengthCRecurse(mult, *one) + _calcCubicArcLengthCRecurse(mult, *two)
+    arch = abs(p0 - p3)
+    box = abs(p0 - p1) + abs(p1 - p2) + abs(p2 - p3)
+    if arch * mult >= box:
+        return (arch + box) * 0.5
+    else:
+        one, two = _split_cubic_into_two(p0, p1, p2, p3)
+        return _calcCubicArcLengthCRecurse(mult, *one) + _calcCubicArcLengthCRecurse(
+            mult, *two
+        )
+
 
 def calcCubicArcLengthC(pt1, pt2, pt3, pt4, tolerance=0.005):
     """Calculates the arc length for a cubic Bezier segment.
@@ -70,7 +90,7 @@ def calcCubicArcLengthC(pt1, pt2, pt3, pt4, tolerance=0.005):
     Returns:
         Arc length value.
     """
-    mult = 1. + 1.5 * tolerance # The 1.5 is a empirical hack; no math
+    mult = 1.0 + 1.5 * tolerance  # The 1.5 is a empirical hack; no math
     return _calcCubicArcLengthCRecurse(mult, pt1, pt2, pt3, pt4)
 
 
@@ -85,7 +105,7 @@ def _dot(v1, v2):
 def _intSecAtan(x):
     # In : sympy.integrate(sp.sec(sp.atan(x)))
     # Out: x*sqrt(x**2 + 1)/2 + asinh(x)/2
-    return x * math.sqrt(x**2 + 1)/2 + math.asinh(x)/2
+    return x * math.sqrt(x ** 2 + 1) / 2 + math.asinh(x) / 2
 
 
 def calcQuadraticArcLength(pt1, pt2, pt3):
@@ -141,16 +161,16 @@ def calcQuadraticArcLengthC(pt1, pt2, pt3):
     d = d1 - d0
     n = d * 1j
     scale = abs(n)
-    if scale == 0.:
-        return abs(pt3-pt1)
-    origDist = _dot(n,d0)
+    if scale == 0.0:
+        return abs(pt3 - pt1)
+    origDist = _dot(n, d0)
     if abs(origDist) < epsilon:
-        if _dot(d0,d1) >= 0:
-            return abs(pt3-pt1)
+        if _dot(d0, d1) >= 0:
+            return abs(pt3 - pt1)
         a, b = abs(d0), abs(d1)
-        return (a*a + b*b) / (a+b)
-    x0 = _dot(d,d0) / origDist
-    x1 = _dot(d,d1) / origDist
+        return (a * a + b * b) / (a + b)
+    x0 = _dot(d, d0) / origDist
+    x1 = _dot(d, d1) / origDist
     Len = abs(2 * (_intSecAtan(x1) - _intSecAtan(x0)) * origDist / (scale * (x1 - x0)))
     return Len
 
@@ -190,13 +210,17 @@ def approximateQuadraticArcLengthC(pt1, pt2, pt3):
     # to be integrated with the best-matching fifth-degree polynomial
     # approximation of it.
     #
-    #https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Legendre_quadrature
+    # https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Legendre_quadrature
 
     # abs(BezierCurveC[2].diff(t).subs({t:T})) for T in sorted(.5, .5±sqrt(3/5)/2),
     # weighted 5/18, 8/18, 5/18 respectively.
-    v0 = abs(-0.492943519233745*pt1 + 0.430331482911935*pt2 + 0.0626120363218102*pt3)
-    v1 = abs(pt3-pt1)*0.4444444444444444
-    v2 = abs(-0.0626120363218102*pt1 - 0.430331482911935*pt2 + 0.492943519233745*pt3)
+    v0 = abs(
+        -0.492943519233745 * pt1 + 0.430331482911935 * pt2 + 0.0626120363218102 * pt3
+    )
+    v1 = abs(pt3 - pt1) * 0.4444444444444444
+    v2 = abs(
+        -0.0626120363218102 * pt1 - 0.430331482911935 * pt2 + 0.492943519233745 * pt3
+    )
 
     return v0 + v1 + v2
 
@@ -220,14 +244,18 @@ def calcQuadraticBounds(pt1, pt2, pt3):
         (0.0, 0.0, 100, 100)
     """
     (ax, ay), (bx, by), (cx, cy) = calcQuadraticParameters(pt1, pt2, pt3)
-    ax2 = ax*2.0
-    ay2 = ay*2.0
+    ax2 = ax * 2.0
+    ay2 = ay * 2.0
     roots = []
     if ax2 != 0:
-        roots.append(-bx/ax2)
+        roots.append(-bx / ax2)
     if ay2 != 0:
-        roots.append(-by/ay2)
-    points = [(ax*t*t + bx*t + cx, ay*t*t + by*t + cy) for t in roots if 0 <= t < 1] + [pt1, pt3]
+        roots.append(-by / ay2)
+    points = [
+        (ax * t * t + bx * t + cx, ay * t * t + by * t + cy)
+        for t in roots
+        if 0 <= t < 1
+    ] + [pt1, pt3]
     return calcBounds(points)
 
 
@@ -256,7 +284,9 @@ def approximateCubicArcLength(pt1, pt2, pt3, pt4):
         >>> approximateCubicArcLength((0, 0), (50, 0), (100, -50), (-50, 0)) # cusp
         154.80848416537057
     """
-    return approximateCubicArcLengthC(complex(*pt1), complex(*pt2), complex(*pt3), complex(*pt4))
+    return approximateCubicArcLengthC(
+        complex(*pt1), complex(*pt2), complex(*pt3), complex(*pt4)
+    )
 
 
 def approximateCubicArcLengthC(pt1, pt2, pt3, pt4):
@@ -276,11 +306,21 @@ def approximateCubicArcLengthC(pt1, pt2, pt3, pt4):
 
     # abs(BezierCurveC[3].diff(t).subs({t:T})) for T in sorted(0, .5±(3/7)**.5/2, .5, 1),
     # weighted 1/20, 49/180, 32/90, 49/180, 1/20 respectively.
-    v0 = abs(pt2-pt1)*.15
-    v1 = abs(-0.558983582205757*pt1 + 0.325650248872424*pt2 + 0.208983582205757*pt3 + 0.024349751127576*pt4)
-    v2 = abs(pt4-pt1+pt3-pt2)*0.26666666666666666
-    v3 = abs(-0.024349751127576*pt1 - 0.208983582205757*pt2 - 0.325650248872424*pt3 + 0.558983582205757*pt4)
-    v4 = abs(pt4-pt3)*.15
+    v0 = abs(pt2 - pt1) * 0.15
+    v1 = abs(
+        -0.558983582205757 * pt1
+        + 0.325650248872424 * pt2
+        + 0.208983582205757 * pt3
+        + 0.024349751127576 * pt4
+    )
+    v2 = abs(pt4 - pt1 + pt3 - pt2) * 0.26666666666666666
+    v3 = abs(
+        -0.024349751127576 * pt1
+        - 0.208983582205757 * pt2
+        - 0.325650248872424 * pt3
+        + 0.558983582205757 * pt4
+    )
+    v4 = abs(pt4 - pt3) * 0.15
 
     return v0 + v1 + v2 + v3 + v4
 
@@ -313,7 +353,13 @@ def calcCubicBounds(pt1, pt2, pt3, pt4):
     yRoots = [t for t in solveQuadratic(ay3, by2, cy) if 0 <= t < 1]
     roots = xRoots + yRoots
 
-    points = [(ax*t*t*t + bx*t*t + cx * t + dx, ay*t*t*t + by*t*t + cy * t + dy) for t in roots] + [pt1, pt4]
+    points = [
+        (
+            ax * t * t * t + bx * t * t + cx * t + dx,
+            ay * t * t * t + by * t * t + cy * t + dy,
+        )
+        for t in roots
+    ] + [pt1, pt4]
     return calcBounds(points)
 
 
@@ -356,8 +402,8 @@ def splitLine(pt1, pt2, where, isHorizontal):
     pt1x, pt1y = pt1
     pt2x, pt2y = pt2
 
-    ax = (pt2x - pt1x)
-    ay = (pt2y - pt1y)
+    ax = pt2x - pt1x
+    ay = pt2y - pt1y
 
     bx = pt1x
     by = pt1y
@@ -410,8 +456,9 @@ def splitQuadratic(pt1, pt2, pt3, where, isHorizontal):
         ((50, 50), (75, 50), (100, 0))
     """
     a, b, c = calcQuadraticParameters(pt1, pt2, pt3)
-    solutions = solveQuadratic(a[isHorizontal], b[isHorizontal],
-        c[isHorizontal] - where)
+    solutions = solveQuadratic(
+        a[isHorizontal], b[isHorizontal], c[isHorizontal] - where
+    )
     solutions = sorted([t for t in solutions if 0 <= t < 1])
     if not solutions:
         return [(pt1, pt2, pt3)]
@@ -446,8 +493,9 @@ def splitCubic(pt1, pt2, pt3, pt4, where, isHorizontal):
         ((92.5259, 25), (95.202, 17.5085), (97.7062, 9.17517), (100, 1.77636e-15))
     """
     a, b, c, d = calcCubicParameters(pt1, pt2, pt3, pt4)
-    solutions = solveCubic(a[isHorizontal], b[isHorizontal], c[isHorizontal],
-        d[isHorizontal] - where)
+    solutions = solveCubic(
+        a[isHorizontal], b[isHorizontal], c[isHorizontal], d[isHorizontal] - where
+    )
     solutions = sorted([t for t in solutions if 0 <= t < 1])
     if not solutions:
         return [(pt1, pt2, pt3, pt4)]
@@ -512,17 +560,17 @@ def _splitQuadraticAtT(a, b, c, *ts):
     cx, cy = c
     for i in range(len(ts) - 1):
         t1 = ts[i]
-        t2 = ts[i+1]
-        delta = (t2 - t1)
+        t2 = ts[i + 1]
+        delta = t2 - t1
         # calc new a, b and c
-        delta_2 = delta*delta
+        delta_2 = delta * delta
         a1x = ax * delta_2
         a1y = ay * delta_2
-        b1x = (2*ax*t1 + bx) * delta
-        b1y = (2*ay*t1 + by) * delta
-        t1_2 = t1*t1
-        c1x = ax*t1_2 + bx*t1 + cx
-        c1y = ay*t1_2 + by*t1 + cy
+        b1x = (2 * ax * t1 + bx) * delta
+        b1y = (2 * ay * t1 + by) * delta
+        t1_2 = t1 * t1
+        c1x = ax * t1_2 + bx * t1 + cx
+        c1y = ay * t1_2 + by * t1 + cy
 
         pt1, pt2, pt3 = calcQuadraticPoints((a1x, a1y), (b1x, b1y), (c1x, c1y))
         segments.append((pt1, pt2, pt3))
@@ -540,24 +588,26 @@ def _splitCubicAtT(a, b, c, d, *ts):
     dx, dy = d
     for i in range(len(ts) - 1):
         t1 = ts[i]
-        t2 = ts[i+1]
-        delta = (t2 - t1)
+        t2 = ts[i + 1]
+        delta = t2 - t1
 
-        delta_2 = delta*delta
-        delta_3 = delta*delta_2
-        t1_2 = t1*t1
-        t1_3 = t1*t1_2
+        delta_2 = delta * delta
+        delta_3 = delta * delta_2
+        t1_2 = t1 * t1
+        t1_3 = t1 * t1_2
 
         # calc new a, b, c and d
         a1x = ax * delta_3
         a1y = ay * delta_3
-        b1x = (3*ax*t1 + bx) * delta_2
-        b1y = (3*ay*t1 + by) * delta_2
-        c1x = (2*bx*t1 + cx + 3*ax*t1_2) * delta
-        c1y = (2*by*t1 + cy + 3*ay*t1_2) * delta
-        d1x = ax*t1_3 + bx*t1_2 + cx*t1 + dx
-        d1y = ay*t1_3 + by*t1_2 + cy*t1 + dy
-        pt1, pt2, pt3, pt4 = calcCubicPoints((a1x, a1y), (b1x, b1y), (c1x, c1y), (d1x, d1y))
+        b1x = (3 * ax * t1 + bx) * delta_2
+        b1y = (3 * ay * t1 + by) * delta_2
+        c1x = (2 * bx * t1 + cx + 3 * ax * t1_2) * delta
+        c1y = (2 * by * t1 + cy + 3 * ay * t1_2) * delta
+        d1x = ax * t1_3 + bx * t1_2 + cx * t1 + dx
+        d1y = ay * t1_3 + by * t1_2 + cy * t1 + dy
+        pt1, pt2, pt3, pt4 = calcCubicPoints(
+            (a1x, a1y), (b1x, b1y), (c1x, c1y), (d1x, d1y)
+        )
         segments.append((pt1, pt2, pt3, pt4))
     return segments
 
@@ -569,8 +619,7 @@ def _splitCubicAtT(a, b, c, d, *ts):
 from math import sqrt, acos, cos, pi
 
 
-def solveQuadratic(a, b, c,
-        sqrt=sqrt):
+def solveQuadratic(a, b, c, sqrt=sqrt):
     """Solve a quadratic equation.
 
     Solves *a*x*x + b*x + c = 0* where a, b and c are real.
@@ -590,13 +639,13 @@ def solveQuadratic(a, b, c,
             roots = []
         else:
             # We have a linear equation with 1 root.
-            roots = [-c/b]
+            roots = [-c / b]
     else:
         # We have a true quadratic equation.  Apply the quadratic formula to find two roots.
-        DD = b*b - 4.0*a*c
+        DD = b * b - 4.0 * a * c
         if DD >= 0.0:
             rDD = sqrt(DD)
-            roots = [(-b+rDD)/2.0/a, (-b-rDD)/2.0/a]
+            roots = [(-b + rDD) / 2.0 / a, (-b - rDD) / 2.0 / a]
         else:
             # complex roots, ignore
             roots = []
@@ -646,52 +695,52 @@ def solveCubic(a, b, c, d):
         # returns unreliable results, so we fall back to quad.
         return solveQuadratic(b, c, d)
     a = float(a)
-    a1 = b/a
-    a2 = c/a
-    a3 = d/a
+    a1 = b / a
+    a2 = c / a
+    a3 = d / a
 
-    Q = (a1*a1 - 3.0*a2)/9.0
-    R = (2.0*a1*a1*a1 - 9.0*a1*a2 + 27.0*a3)/54.0
+    Q = (a1 * a1 - 3.0 * a2) / 9.0
+    R = (2.0 * a1 * a1 * a1 - 9.0 * a1 * a2 + 27.0 * a3) / 54.0
 
-    R2 = R*R
-    Q3 = Q*Q*Q
+    R2 = R * R
+    Q3 = Q * Q * Q
     R2 = 0 if R2 < epsilon else R2
     Q3 = 0 if abs(Q3) < epsilon else Q3
 
     R2_Q3 = R2 - Q3
 
-    if R2 == 0. and Q3 == 0.:
-        x = round(-a1/3.0, epsilonDigits)
+    if R2 == 0.0 and Q3 == 0.0:
+        x = round(-a1 / 3.0, epsilonDigits)
         return [x, x, x]
-    elif R2_Q3 <= epsilon * .5:
+    elif R2_Q3 <= epsilon * 0.5:
         # The epsilon * .5 above ensures that Q3 is not zero.
-        theta = acos(max(min(R/sqrt(Q3), 1.0), -1.0))
-        rQ2 = -2.0*sqrt(Q)
-        a1_3 = a1/3.0
-        x0 = rQ2*cos(theta/3.0) - a1_3
-        x1 = rQ2*cos((theta+2.0*pi)/3.0) - a1_3
-        x2 = rQ2*cos((theta+4.0*pi)/3.0) - a1_3
+        theta = acos(max(min(R / sqrt(Q3), 1.0), -1.0))
+        rQ2 = -2.0 * sqrt(Q)
+        a1_3 = a1 / 3.0
+        x0 = rQ2 * cos(theta / 3.0) - a1_3
+        x1 = rQ2 * cos((theta + 2.0 * pi) / 3.0) - a1_3
+        x2 = rQ2 * cos((theta + 4.0 * pi) / 3.0) - a1_3
         x0, x1, x2 = sorted([x0, x1, x2])
         # Merge roots that are close-enough
         if x1 - x0 < epsilon and x2 - x1 < epsilon:
-            x0 = x1 = x2 = round((x0 + x1 + x2) / 3., epsilonDigits)
+            x0 = x1 = x2 = round((x0 + x1 + x2) / 3.0, epsilonDigits)
         elif x1 - x0 < epsilon:
-            x0 = x1 = round((x0 + x1) / 2., epsilonDigits)
+            x0 = x1 = round((x0 + x1) / 2.0, epsilonDigits)
             x2 = round(x2, epsilonDigits)
         elif x2 - x1 < epsilon:
             x0 = round(x0, epsilonDigits)
-            x1 = x2 = round((x1 + x2) / 2., epsilonDigits)
+            x1 = x2 = round((x1 + x2) / 2.0, epsilonDigits)
         else:
             x0 = round(x0, epsilonDigits)
             x1 = round(x1, epsilonDigits)
             x2 = round(x2, epsilonDigits)
         return [x0, x1, x2]
     else:
-        x = pow(sqrt(R2_Q3)+abs(R), 1/3.0)
-        x = x + Q/x
+        x = pow(sqrt(R2_Q3) + abs(R), 1 / 3.0)
+        x = x + Q / x
         if R >= 0.0:
             x = -x
-        x = round(x - a1/3.0, epsilonDigits)
+        x = round(x - a1 / 3.0, epsilonDigits)
         return [x]
 
 
@@ -699,6 +748,7 @@ def solveCubic(a, b, c, d):
 # Conversion routines for points to parameters and vice versa
 #
 
+
 def calcQuadraticParameters(pt1, pt2, pt3):
     x2, y2 = pt2
     x3, y3 = pt3
@@ -753,10 +803,399 @@ def calcCubicPoints(a, b, c, d):
     return (x1, y1), (x2, y2), (x3, y3), (x4, y4)
 
 
+#
+# Point at time
+#
+
+
+def linePointAtT(pt1, pt2, t):
+    """Finds the point at time `t` on a line.
+
+    Args:
+        pt1, pt2: Coordinates of the line as 2D tuples.
+        t: The time along the line.
+
+    Returns:
+        A 2D tuple with the coordinates of the point.
+    """
+    return ((pt1[0] * (1 - t) + pt2[0] * t), (pt1[1] * (1 - t) + pt2[1] * t))
+
+
+def quadraticPointAtT(pt1, pt2, pt3, t):
+    """Finds the point at time `t` on a quadratic curve.
+
+    Args:
+        pt1, pt2, pt3: Coordinates of the curve as 2D tuples.
+        t: The time along the curve.
+
+    Returns:
+        A 2D tuple with the coordinates of the point.
+    """
+    x = (1 - t) * (1 - t) * pt1[0] + 2 * (1 - t) * t * pt2[0] + t * t * pt3[0]
+    y = (1 - t) * (1 - t) * pt1[1] + 2 * (1 - t) * t * pt2[1] + t * t * pt3[1]
+    return (x, y)
+
+
+def cubicPointAtT(pt1, pt2, pt3, pt4, t):
+    """Finds the point at time `t` on a cubic curve.
+
+    Args:
+        pt1, pt2, pt3, pt4: Coordinates of the curve as 2D tuples.
+        t: The time along the curve.
+
+    Returns:
+        A 2D tuple with the coordinates of the point.
+    """
+    x = (
+        (1 - t) * (1 - t) * (1 - t) * pt1[0]
+        + 3 * (1 - t) * (1 - t) * t * pt2[0]
+        + 3 * (1 - t) * t * t * pt3[0]
+        + t * t * t * pt4[0]
+    )
+    y = (
+        (1 - t) * (1 - t) * (1 - t) * pt1[1]
+        + 3 * (1 - t) * (1 - t) * t * pt2[1]
+        + 3 * (1 - t) * t * t * pt3[1]
+        + t * t * t * pt4[1]
+    )
+    return (x, y)
+
+
+def segmentPointAtT(seg, t):
+    if len(seg) == 2:
+        return linePointAtT(*seg, t)
+    elif len(seg) == 3:
+        return quadraticPointAtT(*seg, t)
+    elif len(seg) == 4:
+        return cubicPointAtT(*seg, t)
+    raise ValueError("Unknown curve degree")
+
+
+#
+# Intersection finders
+#
+
+
+def _line_t_of_pt(s, e, pt):
+    sx, sy = s
+    ex, ey = e
+    px, py = pt
+    if not math.isclose(sx, ex):
+        return (px - sx) / (ex - sx)
+    if not math.isclose(sy, ey):
+        return (py - sy) / (ey - sy)
+    # Line is a point!
+    return -1
+
+
+def _both_points_are_on_same_side_of_origin(a, b, origin):
+    xDiff = (a[0] - origin[0]) * (b[0] - origin[0])
+    yDiff = (a[1] - origin[1]) * (b[1] - origin[1])
+    return not (xDiff <= 0.0 and yDiff <= 0.0)
+
+
+def lineLineIntersections(s1, e1, s2, e2):
+    """Finds intersections between two line segments.
+
+    Args:
+        s1, e1: Coordinates of the first line as 2D tuples.
+        s2, e2: Coordinates of the second line as 2D tuples.
+
+    Returns:
+        A list of ``Intersection`` objects, each object having ``pt``, ``t1``
+        and ``t2`` attributes containing the intersection point, time on first
+        segment and time on second segment respectively.
+
+    Examples::
+
+        >>> a = lineLineIntersections( (310,389), (453, 222), (289, 251), (447, 367))
+        >>> len(a)
+        1
+        >>> intersection = a[0]
+        >>> intersection.pt
+        (374.44882952482897, 313.73458370177315)
+        >>> (intersection.t1, intersection.t2)
+        (0.45069111555824454, 0.5408153767394238)
+    """
+    s1x, s1y = s1
+    e1x, e1y = e1
+    s2x, s2y = s2
+    e2x, e2y = e2
+    if (
+        math.isclose(s2x, e2x) and math.isclose(s1x, e1x) and not math.isclose(s1x, s2x)
+    ):  # Parallel vertical
+        return []
+    if (
+        math.isclose(s2y, e2y) and math.isclose(s1y, e1y) and not math.isclose(s1y, s2y)
+    ):  # Parallel horizontal
+        return []
+    if math.isclose(s2x, e2x) and math.isclose(s2y, e2y):  # Line segment is tiny
+        return []
+    if math.isclose(s1x, e1x) and math.isclose(s1y, e1y):  # Line segment is tiny
+        return []
+    if math.isclose(e1x, s1x):
+        x = s1x
+        slope34 = (e2y - s2y) / (e2x - s2x)
+        y = slope34 * (x - s2x) + s2y
+        pt = (x, y)
+        return [
+            Intersection(
+                pt=pt, t1=_line_t_of_pt(s1, e1, pt), t2=_line_t_of_pt(s2, e2, pt)
+            )
+        ]
+    if math.isclose(s2x, e2x):
+        x = s2x
+        slope12 = (e1y - s1y) / (e1x - s1x)
+        y = slope12 * (x - s1x) + s1y
+        pt = (x, y)
+        return [
+            Intersection(
+                pt=pt, t1=_line_t_of_pt(s1, e1, pt), t2=_line_t_of_pt(s2, e2, pt)
+            )
+        ]
+
+    slope12 = (e1y - s1y) / (e1x - s1x)
+    slope34 = (e2y - s2y) / (e2x - s2x)
+    if math.isclose(slope12, slope34):
+        return []
+    x = (slope12 * s1x - s1y - slope34 * s2x + s2y) / (slope12 - slope34)
+    y = slope12 * (x - s1x) + s1y
+    pt = (x, y)
+    if _both_points_are_on_same_side_of_origin(
+        pt, e1, s1
+    ) and _both_points_are_on_same_side_of_origin(pt, s2, e2):
+        return [
+            Intersection(
+                pt=pt, t1=_line_t_of_pt(s1, e1, pt), t2=_line_t_of_pt(s2, e2, pt)
+            )
+        ]
+    return []
+
+
+def _alignment_transformation(segment):
+    # Returns a transformation which aligns a segment horizontally at the
+    # origin. Apply this transformation to curves and root-find to find
+    # intersections with the segment.
+    start = segment[0]
+    end = segment[-1]
+    angle = math.atan2(end[1] - start[1], end[0] - start[0])
+    return Identity.rotate(-angle).translate(-start[0], -start[1])
+
+
+def _curve_line_intersections_t(curve, line):
+    aligned_curve = _alignment_transformation(line).transformPoints(curve)
+    if len(curve) == 3:
+        a, b, c = calcQuadraticParameters(*aligned_curve)
+        intersections = solveQuadratic(a[1], b[1], c[1])
+    elif len(curve) == 4:
+        a, b, c, d = calcCubicParameters(*aligned_curve)
+        intersections = solveCubic(a[1], b[1], c[1], d[1])
+    else:
+        raise ValueError("Unknown curve degree")
+    return sorted([i for i in intersections if 0.0 <= i <= 1])
+
+
+def curveLineIntersections(curve, line):
+    """Finds intersections between a curve and a line.
+
+    Args:
+        curve: List of coordinates of the curve segment as 2D tuples.
+        line: List of coordinates of the line segment as 2D tuples.
+
+    Returns:
+        A list of ``Intersection`` objects, each object having ``pt``, ``t1``
+        and ``t2`` attributes containing the intersection point, time on first
+        segment and time on second segment respectively.
+
+    Examples::
+        >>> curve = [ (100, 240), (30, 60), (210, 230), (160, 30) ]
+        >>> line  = [ (25, 260), (230, 20) ]
+        >>> intersections = curveLineIntersections(curve, line)
+        >>> len(intersections)
+        3
+        >>> intersections[0].pt
+        (84.90010344084885, 189.87306176459828)
+    """
+    if len(curve) == 3:
+        pointFinder = quadraticPointAtT
+    elif len(curve) == 4:
+        pointFinder = cubicPointAtT
+    else:
+        raise ValueError("Unknown curve degree")
+    intersections = []
+    for t in _curve_line_intersections_t(curve, line):
+        pt = pointFinder(*curve, t)
+        intersections.append(Intersection(pt=pt, t1=t, t2=_line_t_of_pt(*line, pt)))
+    return intersections
+
+
+def _curve_bounds(c):
+    if len(c) == 3:
+        return calcQuadraticBounds(*c)
+    elif len(c) == 4:
+        return calcCubicBounds(*c)
+    raise ValueError("Unknown curve degree")
+
+
+def _split_segment_at_t(c, t):
+    if len(c) == 2:
+        s, e = c
+        midpoint = linePointAtT(s, e, t)
+        return [(s, midpoint), (midpoint, e)]
+    if len(c) == 3:
+        return splitQuadraticAtT(*c, t)
+    elif len(c) == 4:
+        return splitCubicAtT(*c, t)
+    raise ValueError("Unknown curve degree")
+
+
+def _curve_curve_intersections_t(
+    curve1, curve2, precision=1e-3, range1=None, range2=None
+):
+    bounds1 = _curve_bounds(curve1)
+    bounds2 = _curve_bounds(curve2)
+
+    if not range1:
+        range1 = (0.0, 1.0)
+    if not range2:
+        range2 = (0.0, 1.0)
+
+    # If bounds don't intersect, go home
+    intersects, _ = sectRect(bounds1, bounds2)
+    if not intersects:
+        return []
+
+    def midpoint(r):
+        return 0.5 * (r[0] + r[1])
+
+    # If they do overlap but they're tiny, approximate
+    if rectArea(bounds1) < precision and rectArea(bounds2) < precision:
+        return [(midpoint(range1), midpoint(range2))]
+
+    c11, c12 = _split_segment_at_t(curve1, 0.5)
+    c11_range = (range1[0], midpoint(range1))
+    c12_range = (midpoint(range1), range1[1])
+
+    c21, c22 = _split_segment_at_t(curve2, 0.5)
+    c21_range = (range2[0], midpoint(range2))
+    c22_range = (midpoint(range2), range2[1])
+
+    found = []
+    found.extend(
+        _curve_curve_intersections_t(
+            c11, c21, precision, range1=c11_range, range2=c21_range
+        )
+    )
+    found.extend(
+        _curve_curve_intersections_t(
+            c12, c21, precision, range1=c12_range, range2=c21_range
+        )
+    )
+    found.extend(
+        _curve_curve_intersections_t(
+            c11, c22, precision, range1=c11_range, range2=c22_range
+        )
+    )
+    found.extend(
+        _curve_curve_intersections_t(
+            c12, c22, precision, range1=c12_range, range2=c22_range
+        )
+    )
+
+    unique_key = lambda ts: (int(ts[0] / precision), int(ts[1] / precision))
+    seen = set()
+    unique_values = []
+
+    for ts in found:
+        key = unique_key(ts)
+        if key in seen:
+            continue
+        seen.add(key)
+        unique_values.append(ts)
+
+    return unique_values
+
+
+def curveCurveIntersections(curve1, curve2):
+    """Finds intersections between a curve and a curve.
+
+    Args:
+        curve1: List of coordinates of the first curve segment as 2D tuples.
+        curve2: List of coordinates of the second curve segment as 2D tuples.
+
+    Returns:
+        A list of ``Intersection`` objects, each object having ``pt``, ``t1``
+        and ``t2`` attributes containing the intersection point, time on first
+        segment and time on second segment respectively.
+
+    Examples::
+        >>> curve1 = [ (10,100), (90,30), (40,140), (220,220) ]
+        >>> curve2 = [ (5,150), (180,20), (80,250), (210,190) ]
+        >>> intersections = curveCurveIntersections(curve1, curve2)
+        >>> len(intersections)
+        3
+        >>> intersections[0].pt
+        (81.7831487395506, 109.88904552375288)
+    """
+    intersection_ts = _curve_curve_intersections_t(curve1, curve2)
+    return [
+        Intersection(pt=segmentPointAtT(curve1, ts[0]), t1=ts[0], t2=ts[1])
+        for ts in intersection_ts
+    ]
+
+
+def segmentSegmentIntersections(seg1, seg2):
+    """Finds intersections between two segments.
+
+    Args:
+        seg1: List of coordinates of the first segment as 2D tuples.
+        seg2: List of coordinates of the second segment as 2D tuples.
+
+    Returns:
+        A list of ``Intersection`` objects, each object having ``pt``, ``t1``
+        and ``t2`` attributes containing the intersection point, time on first
+        segment and time on second segment respectively.
+
+    Examples::
+        >>> curve1 = [ (10,100), (90,30), (40,140), (220,220) ]
+        >>> curve2 = [ (5,150), (180,20), (80,250), (210,190) ]
+        >>> intersections = segmentSegmentIntersections(curve1, curve2)
+        >>> len(intersections)
+        3
+        >>> intersections[0].pt
+        (81.7831487395506, 109.88904552375288)
+        >>> curve3 = [ (100, 240), (30, 60), (210, 230), (160, 30) ]
+        >>> line  = [ (25, 260), (230, 20) ]
+        >>> intersections = segmentSegmentIntersections(curve3, line)
+        >>> len(intersections)
+        3
+        >>> intersections[0].pt
+        (84.90010344084885, 189.87306176459828)
+
+    """
+    # Arrange by degree
+    swapped = False
+    if len(seg2) > len(seg1):
+        seg2, seg1 = seg1, seg2
+        swapped = True
+    if len(seg1) > 2:
+        if len(seg2) > 2:
+            intersections = curveCurveIntersections(seg1, seg2)
+        else:
+            intersections = curveLineIntersections(seg1, seg2)
+    elif len(seg1) == 2 and len(seg2) == 2:
+        intersections = lineLineIntersections(*seg1, *seg2)
+    else:
+        raise ValueError("Couldn't work out which intersection function to use")
+    if not swapped:
+        return intersections
+    return [Intersection(pt=i.pt, t1=i.t2, t2=i.t1) for i in intersections]
+
+
 def _segmentrepr(obj):
     """
-        >>> _segmentrepr([1, [2, 3], [], [[2, [3, 4], [0.1, 2.2]]]])
-        '(1, (2, 3), (), ((2, (3, 4), (0.1, 2.2))))'
+    >>> _segmentrepr([1, [2, 3], [], [[2, [3, 4], [0.1, 2.2]]]])
+    '(1, (2, 3), (), ((2, (3, 4), (0.1, 2.2))))'
     """
     try:
         it = iter(obj)
@@ -773,7 +1212,9 @@ def printSegments(segments):
     for segment in segments:
         print(_segmentrepr(segment))
 
+
 if __name__ == "__main__":
     import sys
     import doctest
+
     sys.exit(doctest.testmod().failed)