Can anyone suggest a fix or an alternate route to find the solutions to this system? In particular I only care about solutions (s,t) in [0,1]x[0,1]
Note: I'm looking for the intersection of two cubic Bezier curves here. I need the method to be guaranteed to find all solutions and hopefully within a reasonable amount of time (for my use this means a few seconds per pair of curves).
I tried using sympy but both solve() and solve_poly_system() returned empty lists.
Here's my code:
from sympy.solvers import solve_poly_system, solve
from sympy.abc import s,t
#here are two cubics. I'm looking for their intersection in [0,1]x[0,1]:
cub1 = 600*s**3 - 1037*s**2 + 274*s + 1237*t**3 - 2177*t**2 + 642*t + 77
cub2 = -534*s**3 + 582*s**2 + 437*s + 740*t**3 - 1817*t**2 + 1414*t - 548
#I know such a solution exists (from plotting these curves) and fsolve finds an approximation of it no problem:
from scipy.optimize import fsolve
fcub1 = lambda (s,t): 600*s**3 - 1037*s**2 + 274*s + 1237*t**3 - 2177*t**2 + 642*t + 77
fcub2 = lambda (s,t):-534*s**3 + 582*s**2 + 437*s + 740*t**3 - 1817*t**2 + 1414*t - 548
F = lambda x: [fcub1(x),fcub2(x)]
print 'fsolve gives (s,t) = ' + str(fsolve(F,(0.5,0.5)))
print 'F(s,t) = ' + str(F(fsolve(F,(0.5,0.5))))
#solve returns an empty list
print solve([cub1,cub2])
#solve_poly_system returns a DomainError: can't compute a Groebner basis over RR
print solve_poly_system([cub1,cub2])
This outputs:
fsolve gives (s,t) = [ 0.35114023 0.50444115]
F(s,t) = [4.5474735088646412e-13, 0.0]
[]
[]
Thanks for reading!
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.