伊夫的解决方案运行良好。这里是我的情况下，代码可以帮助任何人：
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
而且万一有人想它，我写了编码德卡斯特里奥算法为分裂任意程度的贝塞尔曲线。在上面的代码中，我只是用halveCubic来代替它，这只是de Casteljau的方法，但是明确地将其限制为立方案（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)
我希望能帮助那里的人！ 再次感谢你帮助伊夫！
我认为这不能由sympy解决。例外是因为您在评估中使用了浮点值。在'solve_poly_system'的源代码中有一行检查顺序：'如果a <= 2且b <= 2且c <= 2且d <= 2：' – HYRY
@HYRY：哦，你说得对DomainError。我的意思是删除所有的小数来检查，但我错过了一个！ 另外，你是否说你相信solve_poly_system只能用于每个多项式的阶数<3的系统？ – AndyP