贝赛尔曲线 和 椭圆 [一]

import sympy as sp

# 椭圆曲线长/短半轴
a, b = sp.symbols("a b")

# 假定 第一象限的 1/4 的椭圆 , 可以用 3阶 贝赛尔曲线模拟, 我们假定4个控制点为 p1(0,b),p2(x1,b),p3(a,y1),p4(a,0)
p1, p2, p3, p4 = sp.symbols("p1 p2 p3 p4")

# 贝赛曲线参数 p1*(1-t)**3 + 3*p2*(1-t)*(1-t)*t + 3*p3*(1-t)*t*t +p4*t*t*t
t = sp.symbols("t")

# 我们假定 (y1/x1)=(b/a) = k ,  x1=a*k,y1=b*k
x1, y1, k = sp.symbols("x1 y1 k")

# bezier 曲线主函数
f1 = (
    (1 - t) ** 3 * p1
    + 3 * (1 - t) ** 2 * t * p2
    + 3 * (1 - t) * t * t * p3
    + t**3 * p4
)

# x 的 方程
fx = f1.subs(p1, 0).subs(p2, x1).subs(p3, a).subs(p4, a)
# y 的 方程
fy = f1.subs(p1, b).subs(p2, b).subs(p3, y1).subs(p4, 0)


# 用 x1=a*k, y1*b*k 来 转化 fx,fy 为 k 的方程
fx = fx.subs(x1, a * k)
fy = fy.subs(y1, b * k)

# 椭圆的方程为 (x*x)/(a*a)+(y*y)/(b*b)=1
# 那么, 两边同时 -1 ,把左边 设为一个函数 f(x,y) ,当 f(x,y)为 0 时, 点(x,y) 刚好在圆上


fxy = fx * fx / a / a + fy * fy / b / b - 1


# 把fxy 按k做整理
fxy = sp.apart(fxy, k)

print(fxy)

print("-" * 80)
# fxy,是关于 k 的一元二次方程,那么,我们解方程,可以求得 k 的两个解
k1, k2 = sp.solve(fxy, k)
print(k1)
print(k2)
print("-" * 80)

# 观察k1,k2的表达式, 当 0<=t<=1 时, x1>0,y1>0,a>0,b>0 , k>0 可以求得

k_b = 4 * t**2 - 4 * t - 1
k_d = sp.sqrt(-4 * t**2 + 4 * t + 7)
k_a = 3 * (2 * t**2 - 2 * t + 1)

kk = (k_b + k_d) / k_a


# 此时,我们假定, t=1/3, 可以求到 k
# (sqrt(71)*3-17)/15
# 0.551896621301938

kk = kk.subs(t, 1 / 3)


# 我们来反向验证一下.
# 假定 k= (sqrt(71)*3-17)/15
kv = (sp.sqrt(71) * 3 - 17) / 15

# 分别用 t 来表示 x,y,
px = sp.apart(fx.subs(k, kv), t)
py = sp.apart(fy.subs(k, kv), t)
# 把 px, py 代入 椭圆 的 函数
pxy = px * px / a / a + py * py / b / b - 1
# 函数整理
pxy = sp.apart(pxy, t)
# 因式分解
pxy = sp.factor(pxy)
print(pxy)

# f(t) = -4*t**2*(-76 + 9*sqrt(71))*(t - 1)**2*(3*t - 2)*(3*t - 1)/25
# 有4个因子, t, t-1,3t-1,3t-2 ,

# 那么 当 t in [0,1/3,2/3,1]时, k = (sqrt(71)*3-17)/15

# 这个贝赛尔曲线, 在第一象限 有 4 个交点.
# p1( 0,            b)
# pa( fx(k,t=1/3),  fy(k,t=1/3) )
# pb( fx(k,t=2/3),  fy(k,t=2/3) )
# p4( a,            0)

# 控制点为
# p1(0,   b)
# p2(a*k, b)
# p3(a,   b*k )
# p4(a,   0)


# 由于贝赛尔曲线可以反向, 所以 把  p1 与 p4,  p2与p3, 同时交换位置,所得的曲线 轨迹相同, 方向 相反. 同样成立.

标签：p2,椭圆,subs,p4,赛尔,曲线,p1,y1,sp
来源： https://www.cnblogs.com/unm001/p/16542364.html