Solving a system of cubic polynomials (to find intersection of Bezier curves)

Question 1

For the intersection of Béziers, there are better approaches. (http://pomax.github.io/bezierinfo/#curveintersection, http://www.tsplines.com/technology/edu/CurveIntersection.pdf).

My recommendation for a simple solution: implement the Bezier subdivision algorithm (http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/Bezier/bezier-sub.html). For both curves, compute the bounding box of the control vertices. If they overlap, an intersection is possible, subdivide and repeat the process with the halves (there will be four comparisons this time). Continue recursively.

You mustn't fear exponential explosion (1, 4, 16, 256... comparisons), as rather soon many boxes will stop overlapping.

Note that in theory you can work with the convex hulls of the control points, but in practice a simple bounding box is quite enough and much easier to work with.

Question 2

Yves's solution worked well. Here's my code in case it helps anyone:

from math import sqrt
def cubicCurve(P,t):
    return P[0]*(1-t)**3 + 3*P[1]*t*(1-t)**2 + 3*P[2]*(1-t)*t**2 + P[3]*t**3
def cubicMinMax_x(points):
    local_extremizers = [0,1]
    a = [p.real for p in points]
    delta = a[1]**2 - (a[0] + a[1]*a[2] + a[2]**2 + (a[0] - a[1])*a[3])
    if delta>=0:
        sqdelta = sqrt(delta)/(a[0] - 3*a[1] + 3*a[2] - a[3])
        tau = a[0] - 2*a[1] + a[2]
        r1 = tau+sqdelta
        r2 = tau-sqdelta
        if 0<r1<1:
            local_extremizers.append(r1)
        if 0<r2<1:
            local_extremizers.append(r2)
    localExtrema = [cubicCurve(a,t) for t in local_extremizers]
    return min(localExtrema),max(localExtrema)
def cubicMinMax_y(points):
    return cubicMinMax_x([-1j*p for p in points])
def intervalIntersectionWidth(a,b,c,d): #returns width of the intersection of intervals [a,b] and [c,d]  (thinking of these as intervals on the real number line)
    return max(0, min(b, d) - max(a, c))
def cubicBoundingBoxesIntersect(cubs):#INPUT: 2-tuple of cubics (given bu control points) #OUTPUT: boolean
    x1min,x1max = cubicMinMax_x(cubs[0])
    y1min,y1max = cubicMinMax_y(cubs[0])
    x2min,x2max = cubicMinMax_x(cubs[1])
    y2min,y2max = cubicMinMax_y(cubs[1])
    if intervalIntersectionWidth(x1min,x1max,x2min,x2max) and intervalIntersectionWidth(y1min,y1max,y2min,y2max):
        return True
    else:
        return False
def cubicBoundingBoxArea(cub_points):#INPUT: 2-tuple of cubics (given bu control points) #OUTPUT: boolean
    xmin,xmax = cubicMinMax_x(cub_points)
    ymin,ymax = cubicMinMax_y(cub_points)
    return (xmax-xmin)*(ymax-ymin)
def halveCubic(P):
    return ([P[0], (P[0]+P[1])/2, (P[0]+2*P[1]+P[2])/4, (P[0]+3*P[1]+3*P[2]+P[3])/8],[(P[0]+3*P[1]+3*P[2]+P[3])/8,(P[1]+2*P[2]+P[3])/4,(P[2]+P[3])/2,P[3]])
class Pair(object):
    def __init__(self,cub1,cub2,t1,t2):
        self.cub1 = cub1
        self.cub2 = cub2
        self.t1 = t1 #the t value to get the mid point of this curve from cub1
        self.t2 = t2 #the t value to get the mid point of this curve from cub2
def cubicXcubicIntersections(cubs):
#INPUT: a tuple cubs=([P0,P1,P2,P3], [Q0,Q1,Q2,Q3]) defining the two cubic to check for intersections between.  See cubicCurve fcn for definition of P0,...,P3
#OUTPUT: a list of tuples (t,s) in [0,1]x[0,1] such that cubicCurve(cubs[0],t) - cubicCurve(cubs[1],s) < Tol_deC
#Note: This will return exactly one such tuple for each intersection (assuming Tol_deC is small enough)
    Tol_deC = 1 ##### This should be set based on your accuracy needs.  Making it smaller will have relatively little effect on performance.  Mine is set to 1 because this is the area of a pixel in my setup and so the curve (drawn by hand/mouse) is only accurate up to a pixel at most. 
    maxIts =  100 ##### This should be something like maxIts = 1-log(Tol_deC/length)/log(2), where length is the length of the longer of the two cubics, but I'm not actually sure how close to being parameterized by arclength these curves are... so I guess I'll leave that as an exercise for the interested reader :)
    pair_list = [Pair(cubs[0],cubs[1],0.5,0.5)]
    intersection_list = []
    k=0
    while pair_list != []:
        newPairs = []
        delta = 0.5**(k+2)
        for pair in pair_list:
            if cubicBoundingBoxesIntersect((pair.cub1,pair.cub2)):
                if cubicBoundingBoxArea(pair.cub1)<Tol_deC and cubicBoundingBoxArea(pair.cub2)<Tol_deC:
                    intersection_list.append((pair.t1,pair.t2)) #this is the point in the middle of the pair
                    for otherPair in pair_list:
                        if pair.cub1==otherPair.cub1 or pair.cub2==otherPair.cub2 or pair.cub1==otherPair.cub2 or pair.cub2==otherPair.cub1:
                            pair_list.remove(otherPair) #this is just an ad-hoc fix to keep it from repeating intersection points
                else:
                    (c11,c12) = halveCubic(pair.cub1)
                    (t11,t12) = (pair.t1-delta,pair.t1+delta)
                    (c21,c22) = halveCubic(pair.cub2)
                    (t21,t22) = (pair.t2-delta,pair.t2+delta)
                    newPairs += [Pair(c11,c21,t11,t21), Pair(c11,c22,t11,t22), Pair(c12,c21,t12,t21), Pair(c12,c22,t12,t22)]
        pair_list = newPairs
        k += 1
        if k > maxIts:
            raise Exception ("cubicXcubicIntersections has reached maximum iterations without terminating... either there's a problem/bug or you can fix by raising the max iterations or lowering Tol_deC")
    return intersection_list

Also just in case anyone wants it, I wrote coded de Casteljau's algorithm for splitting Bezier curves of arbitrary degree. In the above code I just replace it with halveCubic, which is just de Casteljau's method but made explicit and restricted to the cubic case (with t=0.5).

def splitBezier(points,t):
#returns 2 tuples of control points for the two resulting Bezier curves
    points_left=[]
    points_right=[]
    (points_left,points_right) =  splitBezier_deCasteljau_recursion((points_left,points_right),points,t)
    points_right.reverse()
    return (points_left,points_right)
def splitBezier_deCasteljau_recursion(cub_lr,points,t):
    (cub_left,cub_right)=cub_lr
    if len(points)==1:
        cub_left.append(points[0])
        cub_right.append(points[0])
    else:
        n = len(points)-1
        newPoints=[None]*n
        cub_left.append(points[0])
        cub_right.append(points[n])
        for i in range(n):
            newPoints[i] = (1-t)*points[i] + t*points[i+1]
        (cub_left, cub_right) = splitBezier_deCasteljau_recursion((cub_left,cub_right),newPoints,t)
    return (cub_left, cub_right)

I hope that helps someone out there! Thanks again for you help Yves!

Question 3

What about simplifying the system ?

By a suitable linear combination, you can eliminate one of s^3 or t^3, and solve the remaining second degree equation. By plugging the result in the other equation, you get a single equation in a single unknown.

Or solve the algebraic equation obtained by resultants:

36011610661302281 - 177140805507270756*s - 201454039857766711*s^2 + 1540826307929388607*s^3 + 257712262726095899*s^4 - 4599101672917940010*s^5 + 1114665205197856508*s^6 + 6093758014794453276*s^7 - 5443785088068396888*s^8 + 1347614193395309112*s^9 = 0

(http://www.dr-mikes-maths.com/polynomial-reduction.html)