均匀三次B样条曲线.续 - 技术交流 - TCAX 字幕特效制作工具官方论坛 | ASS | TCAS | Python | Aegisub | Lua - Powered by Discuz!

TCAX 字幕特效制作工具官方论坛 | ASS | TCAS | Python | Aegisub | Lua

 找回密码
 加入社区
查看: 3102|回复: 12

均匀三次B样条曲线.续 [复制链接]

Rank: 4

发表于 2012-8-1 13:06:37 |显示全部楼层
懒得码字了,直接上代码+测试。
注:附上了完整的TCAX下的测试工程,要执行,请把tcCurve.py放到TCAX的util目录下。

实现
tcCurve.py (10.02 KB, 下载次数: 381)
  1. ###########################################################################################################################
  2. ###########################################################################################################################
  3. ###################### Uniform Cubic B-Spline Python Implementation, Object-Oriented Version With Const Speed Utility

  4. class UCBSpline:
  5.     # pre-calculate coefficients
  6.     def __init__(self, p):
  7.         self.p = p          # tuple or list of control points (x, y)
  8.         n = len(p) - 1      # control points p0 ~ pn, i.e., p[0] ~ p[n]
  9.         m = 3 + n + 1       # knots vector u0 ~um, i.e., u[0] ~ u[m]
  10.         self.m = m
  11.         step = 1 / (m - 2 * 3)
  12.         self.u = 3 * [0] + [_u * step for _u in range(m - 2 * 3 + 1)] + [1] * 3     # knots vector
  13.         num = n + 1 - 3     # number of curve segments
  14.         self.num = num
  15.         self.a1 = []
  16.         self.a2 = []
  17.         self.a3 = []
  18.         self.a4 = []
  19.         self.b1 = []
  20.         self.b2 = []
  21.         self.b3 = []
  22.         self.b4 = []
  23.         for i in range(num):
  24.             self.a1.append((-p[i][0] + 3 * p[i + 1][0] - 3 * p[i + 2][0] + p[i + 3][0]) / 6.0)
  25.             self.a2.append((3 * p[i][0] - 6 * p[i + 1][0] + 3 * p[i + 2][0]) / 6.0)
  26.             self.a3.append((-p[i][0] + p[i + 2][0]) / 2.0)
  27.             self.a4.append((p[i][0] + 4 * p[i + 1][0] + p[i + 2][0]) / 6.0)
  28.             self.b1.append((-p[i][1] + 3 * p[i + 1][1] - 3 * p[i + 2][1] + p[i + 3][1]) / 6.0)
  29.             self.b2.append((3 * p[i][1] - 6 * p[i + 1][1] + 3 * p[i + 2][1]) / 6.0)
  30.             self.b3.append((-p[i][1] + p[i + 2][1]) / 2.0)
  31.             self.b4.append((p[i][1] + 4 * p[i + 1][1] + p[i + 2][1]) / 6.0)
  32.         self.SIMPSON = 100  # number of sub-intervals of composite simpson integration, must be even, usually 100 is good enough
  33.         self.NEWTON = 0.0001  # epsilon, effective digits of the root calculated by Newton's method
  34.     # calculate the value at t
  35.     def __call__(self, t):
  36.         if t < 0 or t >= 1:
  37.             return None
  38.         else:       # find the curve segment
  39.             for i in range(3, self.m - 3):
  40.                 if self.u[i] <= t and t < self.u[i + 1]:
  41.                     break
  42.         idx = i - 3 # get coefficients
  43.         a1 = self.a1[idx]
  44.         a2 = self.a2[idx]
  45.         a3 = self.a3[idx]
  46.         a4 = self.a4[idx]
  47.         b1 = self.b1[idx]
  48.         b2 = self.b2[idx]
  49.         b3 = self.b3[idx]
  50.         b4 = self.b4[idx]
  51.         t = (t - self.u[i]) / (self.u[i + 1] - self.u[i])   # map to the current curve segment
  52.         return (a1 * t * t * t + a2 * t * t + a3 * t + a4,
  53.                 b1 * t * t * t + b2 * t * t + b3 * t + b4, 255)    # 255 is the point's alpha value, to meet the TCAX's point structure requirement ((x, y, a), (x, y, a), ...)
  54.     # total arc-length
  55.     def length(self):
  56.         from math import sqrt
  57.         sqr = lambda x: x * x
  58.         l = 0
  59.         for i in range(self.num):
  60.             a1 = self.a1[i]
  61.             a2 = self.a2[i]
  62.             a3 = self.a3[i]
  63.             a4 = self.a4[i]
  64.             b1 = self.b1[i]
  65.             b2 = self.b2[i]
  66.             b3 = self.b3[i]
  67.             b4 = self.b4[i]
  68.             f_grand = lambda t: sqrt(sqr(3 * a1 * t * t + 2 * a2 * t + a3) + sqr(3 * b1 * t * t + 2 * b2 * t + b3))       # integrand of UCB
  69.             l += CompositeSimpson(f_grand, 0, 1, self.SIMPSON)
  70.         return l
  71.     # find t2 according to a specific arc-length s and a starting point t1
  72.     def t2OfLength(self, t1, s):
  73.         if t1 < 0 or t1 >= 1:
  74.             return None
  75.         else:       # find the curve segment
  76.             for i in range(3, self.m - 3):
  77.                 if self.u[i] <= t1 and t1 < self.u[i + 1]:
  78.                     break
  79.         idx = i - 3 # get coefficients
  80.         a1 = self.a1[idx]
  81.         a2 = self.a2[idx]
  82.         a3 = self.a3[idx]
  83.         a4 = self.a4[idx]
  84.         b1 = self.b1[idx]
  85.         b2 = self.b2[idx]
  86.         b3 = self.b3[idx]
  87.         b4 = self.b4[idx]
  88.         t1 = (t1 - self.u[i]) / (self.u[i + 1] - self.u[i])   # map to the current curve segment
  89.         # Assumption: t1 and t are in the same piece of curve
  90.         from math import sqrt
  91.         sqr = lambda x: x * x
  92.         f_grand = lambda t: sqrt(sqr(3 * a1 * t * t + 2 * a2 * t + a3) + sqr(3 * b1 * t * t + 2 * b2 * t + b3))       # integrand of UCB
  93.         s1 = CompositeSimpson(f_grand, t1, 1, self.SIMPSON)     # length from the starting point t1 to the end of the current curve segment
  94.         if s1 >= s:    # t1 and t are in the same segment of curve
  95.             f_grand_d = lambda t: ((3 * a1 * t * t + 2 * a2 * t + a3) * (6 * a1 * t + 2 * a2) + (3 * b1 * t * t + 2 * b2 * t + b3) * (6 * b1 * t + 2 * b2 * t)) / f_grand(t)    # derivative of f_grand
  96.             # using Newton's method to find the root of function f
  97.             f = lambda t: CompositeSimpson(f_grand, t1, t, self.SIMPSON) - s
  98.             fd = lambda t: CompositeSimpsonDerivative(f_grand, f_grand_d, t1, t, self.SIMPSON)   # derivative of f
  99.             return NewtonMethod(f, fd, t1, self.NEWTON) * (self.u[i + 1] - self.u[i]) + self.u[i]   # find the t and map to the whole curve
  100.         # otherwise, t is in next segment(s)
  101.         else:
  102.             if idx + 1 == self.num:     # no more curve segment
  103.                 return 1
  104.             for i in range(idx + 1, self.num):
  105.                 a1 = self.a1[i]
  106.                 a2 = self.a2[i]
  107.                 a3 = self.a3[i]
  108.                 a4 = self.a4[i]
  109.                 b1 = self.b1[i]
  110.                 b2 = self.b2[i]
  111.                 b3 = self.b3[i]
  112.                 b4 = self.b4[i]
  113.                 f_grand = lambda t: sqrt(sqr(3 * a1 * t * t + 2 * a2 * t + a3) + sqr(3 * b1 * t * t + 2 * b2 * t + b3))       # integrand of UCB
  114.                 s0 = s1
  115.                 s1 += CompositeSimpson(f_grand, 0, 1, self.SIMPSON)
  116.                 if s1 >= s:     # we find at which segment the t lies
  117.                     break
  118.             if i == self.num:   # t is beyond the whole curve
  119.                 return 1
  120.             f_grand_d = lambda t: ((3 * a1 * t * t + 2 * a2 * t + a3) * (6 * a1 * t + 2 * a2) + (3 * b1 * t * t + 2 * b2 * t + b3) * (6 * b1 * t + 2 * b2 * t)) / f_grand(t)    # derivative of f_grand
  121.             # using Newton's method to find the root of function f
  122.             f = lambda t: CompositeSimpson(f_grand, 0, t, self.SIMPSON) - (s - s0)
  123.             fd = lambda t: CompositeSimpsonDerivative(f_grand, f_grand_d, 0, t, self.SIMPSON)   # derivative of f
  124.             return NewtonMethod(f, fd, 0, self.NEWTON) * (self.u[i + 3 + 1] - self.u[i + 3]) + self.u[i + 3]   # find the t and map to the whole curve, note, here i is actually idx, should + 3 to be real i

  125. # f=function, a=initial value, b=end value, n=number of intervals of size h, n must be even
  126. def CompositeSimpson(f, a, b, n):
  127.     h = (b - a) / n
  128.     S = f(a)
  129.     for i in range(1, n, 2):
  130.         x = a + h * i
  131.         S += 4 * f(x)
  132.     for i in range(2, n - 1, 2):
  133.         x = a + h * i
  134.         S += 2 * f(x)
  135.     S += f(b)
  136.     return h * S / 3

  137. # f=function, fd=derivative of f, t1=initial value (constant), t=end value (variable), n=number of intervals of size h, n must be even
  138. def CompositeSimpsonDerivative(f, fd, t1, t, n):
  139.     h = (t - t1) / n
  140.     hd = 1 / n
  141.     S = f(t1)
  142.     for i in range(1, n, 2):
  143.         x = t1 + h * i
  144.         S += 4 * f(x)
  145.     for i in range(2, n - 1, 2):
  146.         x = t1 + h * i
  147.         S += 2 * f(x)
  148.     S += f(t)
  149.     part1 = hd * S / 3
  150.     # part2
  151.     S = 0
  152.     for i in range(1, n, 2):
  153.         x = t1 + h * i
  154.         xd = hd * i
  155.         S += 4 * fd(x) * xd
  156.     for i in range(2, n - 1, 2):
  157.         x = t1 + h * i
  158.         xd = hd * i
  159.         S += 2 * fd(x) * xd
  160.     S += fd(t)
  161.     part2 = h * S / 3
  162.     return part1 + part2

  163. # f is the function whose root needs to be found, fd is the derivative of function f, x0 is the starting point, e is an epsilon, i.e., effective digits
  164. def NewtonMethod(f, fd, x0, e):
  165.     x = x0
  166.     while True:
  167.         x_new = x - f(x) / fd(x)
  168.         if abs(x_new - x) < e:
  169.             return x_new
  170.         x = x_new


  171. ###########################################################################################################################
  172. ###########################################################################################################################
  173. ###################### Uniform Cubic B-Spline Python Implementation

  174. def UniformCubicBSpline(p, t):
  175.     n = len(p) - 1     # control points p0 ~ pn, i.e., p[0] ~ p[n]
  176.     m = 3 + n + 1      # knots vector u0 ~um, i.e., u[0] ~ u[m]
  177.     step = 1 / (m - 2 * 3)
  178.     u = 3 * [0] + [_u * step for _u in range(m - 2 * 3 + 1)] + [1] * 3
  179.     if t < 0 or t >= 1:
  180.         return None
  181.     else:
  182.         for i in range(3, m - 3):
  183.             if u[i] <= t and t < u[i + 1]:
  184.                 break
  185.     return UniformCubicBSpline_Calc(p[i - 3], p[i - 2], p[i - 1], p[i], (t - u[i]) / (u[i + 1] - u[i]))    # can use either calc or calc2

  186. def UniformCubicBSpline_Calc(p1, p2, p3, p4, t):
  187.     a1 = (-p1[0] + 3 * p2[0] - 3 * p3[0] + p4[0]) / 6.0
  188.     a2 = (3 * p1[0] - 6 * p2[0] + 3 * p3[0]) / 6.0
  189.     a3 = (-p1[0] + p3[0]) / 2.0
  190.     a4 = (p1[0] + 4 * p2[0] + p3[0]) / 6.0
  191.     b1 = (-p1[1] + 3 * p2[1] - 3 * p3[1] + p4[1]) / 6.0
  192.     b2 = (3 * p1[1] - 6 * p2[1] + 3 * p3[1]) / 6.0
  193.     b3 = (-p1[1] + p3[1]) / 2.0
  194.     b4 = (p1[1] + 4 * p2[1] + p3[1]) / 6.0
  195.     return (a1 * t * t * t + a2 * t * t + a3 * t + a4,
  196.             b1 * t * t * t + b2 * t * t + b3 * t + b4, 255)    # 255 is the point's alpha value, to meet the TCAX's point structure requirement ((x, y, a), (x, y, a), ...)

  197. def UniformCubicBSpline_Calc2(p1, p2, p3, p4, t):
  198.     A1 = (1 - t) * (1 - t) * (1 - t) / 6.0
  199.     A2 = (3 * t * t * t - 6 * t * t + 4) / 6.0
  200.     A3 = (-3 * t * t * t + 3 * t * t + 3 * t + 1) / 6.0
  201.     A4 = t * t * t / 6.0
  202.     return (A1 * p1[0] + A2 * p2[0] + A3 * p3[0] + A4 * p4[0],
  203.             A1 * p1[1] + A2 * p2[1] + A3 * p3[1] + A4 * p4[1], 255)    # 255 is the point's alpha value, to meet the TCAX's point structure requirement ((x, y, a), (x, y, a), ...)
复制代码
测试用例
test_curves.rar (1.48 KB, 下载次数: 5612)
  1. from tcaxPy import *
  2. from util.cairo import *
  3. from util.tcCurve import *

  4. def tcaxPy_User():
  5.     surface = ImageSurface(FORMAT_ARGB32, GetVal(val_ResolutionX), GetVal(val_ResolutionY))
  6.     ctx = Context(surface)
  7.     # 控制点列表1
  8.     P = [(200, 500), (300, 100), (700, 100), (400, 500), (600, 700), (800, 400), (700, 200), (900, 300), (660, 400)]
  9.     # 控制点列表2
  10. #    P = [(640, 0)]
  11. #    for i in range(10):
  12. #        P.extend([(240 + randint(-40, 40), 200 + randint(-50, 50)), (1040 + randint(-40, 40), 200 + randint(-50, 50)), (640, 0)])
  13. #        P.extend([(340 + randint(-60, 60), 600 + randint(-50, 50)), (940 + randint(-60, 60), 600 + randint(-50, 50)), (640, 0)])
  14.     # 控制点列表3
  15. #    P = [(0, 100), (100, 0), (200, 0), (300, 100), (400, 200), (500, 200), (600, 100), (400, 400), (700, 50), (800, 200)]
  16.     # 控制点列表4
  17. #    P = [(300, 100), (400, 200), (500, 200), (600, 100)]
  18.     # 控制点列表5
  19. #    P = [(400, 100), (400, 500), (500, 500), (500, 100)]

  20.     # 绘制控制点
  21.     ctx.set_source_rgb(0, 0, 0)
  22.     num = len(P)
  23.     for i in range(num):
  24.         ctx.arc(P[i][0], P[i][1], 2, 0, 2 * pi)
  25.         ctx.fill()

  26.     # 计算曲线
  27.     ucb = UCBSpline(P)
  28.     L = ucb.length()        # 曲线总长度
  29.     STEP_N = int(L / 10 + 0.5)  # 取样点数
  30.     step_size = 1 / STEP_N
  31.     print('控制点: {0}    曲线总长度: {1:.02f}    所取点数: {2}'.format(len(P), L, STEP_N))
  32.     t1 = 0
  33.     points = []
  34.     for i in range(STEP_N):
  35.         # 普通1, 使用UniformCubicBSpline函数
  36.         #points.append(UniformCubicBSpline(P, i * step_size))
  37.         # 普通2
  38.         #points.append(ucb(i * step_size))
  39.         # 匀速版本
  40.         points.append(ucb(t1))
  41.         t1 = ucb.t2OfLength(t1, step_size * L)

  42.     # 绘制独立点
  43.     ctx.set_source_rgb(0, 1, 0)
  44.     num = len(points)
  45.     for i in range(num):
  46.         ctx.arc(points[i][0], points[i][1], 2, 0, 2 * pi)
  47.         ctx.fill()
  48.     # 绘制连续曲线
  49. #    ctx.move_to(points[0][0], points[0][1])
  50. #    for i in range(1, num):
  51. #        ctx.line_to(points[i][0], points[i][1])
  52. #    ctx.stroke()

  53.     # save the result
  54.     surface.write_to_png(abspath('1.png'))
  55.     PIX = surface.get_pix()
  56.     PIX = PixBlur(PIX, 10)
  57.     SavePix(abspath('2.png'), PIX)
  58.     #sys.exit()
复制代码
若干效果截图

1.png


1 (2).png


计算效率初步测试

Unnamed QQ Screenshot20120801025852.png


1 (3).png



回帖推荐

ミルク 发表于12楼  查看完整内容

根据这篇论文,重新理解了下UCBS,并简化了代码(更好理解)。
2

查看全部评分

Super Moderator

{\clip\t(\clip)}∀.I.R.nesSub

Rank: 6Rank: 6

发表于 2012-8-1 17:33:50 |显示全部楼层
本帖最后由 BurySakura 于 2012-8-1 17:43 编辑

效率相对我那个匀速贝塞尔真心好。
但是这个怎么控制曲线开始,结束的位置?
还有点的导数在哪?

curve.rar (18.44 KB, 下载次数: 470)
test.png

test1.png


psp.希望牛奶提供的大概思路,我看我那曲线是否能改下提高效率。

Rank: 4

发表于 2012-8-1 17:44:27 |显示全部楼层
BurySakura 发表于 2012-8-1 17:33
效率相对我那个匀速贝塞尔真心好。
但是这个怎么控制开始,结束点的位置?
还有点的导数在哪?


开始点和结束点好像只起一个引导作用。可以重叠两个点为开始点,效果就稍微接近一点。
  1. P = [(200, 500), (200, 500), (300, 100), (700, 100), (400, 500), (600, 700), (800, 400), (700, 200), (900, 300), (660, 400), (660, 400)]
复制代码
1.png


导数应该是 f_grand_d, fd, CompositeSimpsonDerivative

Rank: 4

发表于 2012-8-1 17:51:51 |显示全部楼层
补充:对于有n(n比较大)个控制点的情况,直接使用贝塞尔曲线,会造成曲线阶数很高,从而导致效率下降(阶数过高的贝塞尔曲线形状也不容易控制)。
可以用Composite Bezier Spline(复合贝塞尔样条曲线),降低阶数以提高效率。

p.s. 如果你的代码对于4, 5个控制点的情况效率不错,就说明优化的价值(提升的空间)不高了

Super Moderator

{\clip\t(\clip)}∀.I.R.nesSub

Rank: 6Rank: 6

发表于 2012-8-1 18:04:19 |显示全部楼层
本帖最后由 BurySakura 于 2012-8-1 18:17 编辑
ミルク 发表于 2012-8-1 17:51
补充:对于有n(n比较大)个控制点的情况,直接使用贝塞尔曲线,会造成曲线阶数很高,从而导致效率下降(阶 ...


看代码太累,说说你的实现方法。
现在在考虑用Cairo的curve去替代,反正用你编译那个能获得像素信息。

Rank: 4

发表于 2012-8-1 18:30:04 |显示全部楼层
BurySakura 发表于 2012-8-1 18:04
看代码太累,说说你的实现方法。
现在在考虑用Cairo的curve去替代,反正用你编译那个能获得像素信息。 ...

我没用复合贝塞尔样条曲线。关于B样条曲线,这个帖子有些说明了 http://www.tcax.org/forum.php?mod=viewthread&tid=459

匀速的做法和你一样的吧。用Simpson求积分,用Newton求根。

Super Moderator

{\clip\t(\clip)}∀.I.R.nesSub

Rank: 6Rank: 6

发表于 2012-8-1 22:53:55 |显示全部楼层
本帖最后由 BurySakura 于 2012-8-1 23:01 编辑

效率真心美,就是开始跟结束是永远的痛。

Uniform Cubic B-Spline.jpg

Uniform Cubic B-Spline01.jpg

Rank: 4

发表于 2012-8-1 23:17:55 |显示全部楼层
BurySakura 发表于 2012-8-1 22:53
效率真心美,就是开始跟结束是永远的痛。

两个都是用匀速UCBS的么。。

确实我也觉得开始和结束需要再研究研究。。(我也希望能有办法让曲线通过至少开始点和结束点(或者至少是第二点及倒数第二点。。

再看看有没有啥办法。。。

Rank: 4

发表于 2012-8-2 01:14:33 |显示全部楼层
http://www.tcax.org/forum.php?mo ... 4NDExODd8Mnw0NTE%3D

这篇论文貌似有提到。

Unnamed QQ Screenshot20120802011703.png


Unnamed QQ Screenshot20120802011905.png


Unnamed QQ Screenshot20120802011931.png


不过由Uniform变成Non-Uniform的了。。。
明天再看看

Super Moderator

{\clip\t(\clip)}∀.I.R.nesSub

Rank: 6Rank: 6

发表于 2012-8-2 10:01:25 |显示全部楼层
于是最后还是要走NURBS。

Super Moderator

{\clip\t(\clip)}∀.I.R.nesSub

Rank: 6Rank: 6

发表于 2012-8-2 14:54:27 |显示全部楼层
Uniform Cubic B-Spline curve test

未命名11.jpg

Rank: 4

发表于 2012-8-2 19:00:48 |显示全部楼层
A Practical Review of Uniform B-Splines.pdf (53.94 KB, 下载次数: 517)

根据这篇论文,重新理解了下UCBS,并简化了代码(更好理解)。

Cubic B-Spline.png


tcCurve.py (9.27 KB, 下载次数: 343)

Rank: 4

发表于 2012-8-3 18:18:04 |显示全部楼层
后续东西,我都发日记楼里了。

http://www.tcax.org/forum.php?mod=viewthread&tid=469

您需要登录后才可以回帖 登录 | 加入社区

GitHub|TCAX 主页

GMT+8, 2018-8-15 08:00

Powered by Discuz! X2

© 2001-2011 Comsenz Inc.

回顶部
RealH