Yves的解决方案很好。这是我的代码，如果有帮助的话：
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):
return max(0, min(b, d) - max(a, c))
def cubicBoundingBoxesIntersect(cubs):
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):
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
self.t2 = t2
def cubicXcubicIntersections(cubs):
Tol_deC = 1
maxIts = 100
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))
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)
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
此外，以防有人需要，我编写了 Casteljau 算法的代码，用于拆分任意次数的贝塞尔曲线。在上述代码中，我只是用 halveCubic 替换了它，这只是将 Casteljau 方法显式化并限制为三次情况（t=0.5）。
def splitBezier(points,t):
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)
我希望你能帮助某些人，感谢您的帮助，Yves！
solve_poly_system的源代码中，有一行检查顺序：if a <= 2 and b <= 2 and c <= 2 and d <= 2:- HYRY