bezfigs.py 5.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183
  1. import sys
  2. from math import *
  3. import fromcubic
  4. import tocubic
  5. import cornu
  6. def eps_prologue(x0, y0, x1, y1, draw_box = False):
  7. print '%!PS-Adobe-3.0 EPSF'
  8. print '%%BoundingBox:', x0, y0, x1, y1
  9. print '%%EndComments'
  10. print '%%EndProlog'
  11. print '%%Page: 1 1'
  12. if draw_box:
  13. print x0, y0, 'moveto', x0, y1, 'lineto', x1, y1, 'lineto', x1, y0, 'lineto closepath stroke'
  14. def eps_trailer():
  15. print '%%EOF'
  16. def fit_cubic_superfast(z0, z1, arclen, th0, th1, aab):
  17. chord = hypot(z1[0] - z0[0], z1[1] - z0[1])
  18. cth0, sth0 = cos(th0), sin(th0)
  19. cth1, sth1 = -cos(th1), -sin(th1)
  20. armlen = .66667 * arclen
  21. a = armlen * aab
  22. b = armlen - a
  23. bz = [z0, (z0[0] + cth0 * a, z0[1] + sth0 * a),
  24. (z1[0] + cth1 * b, z1[1] + sth1 * b), z1]
  25. return bz
  26. def fit_cubic(z0, z1, arclen, th_fn, fast, aabmin = 0, aabmax = 1.):
  27. chord = hypot(z1[0] - z0[0], z1[1] - z0[1])
  28. if (arclen < 1.000001 * chord):
  29. return [z0, z1], 0
  30. th0 = th_fn(0)
  31. th1 = th_fn(arclen)
  32. imax = 4
  33. jmax = 10
  34. if fast:
  35. imax = 1
  36. jmax = 0
  37. for i in range(imax):
  38. for j in range(jmax + 1):
  39. if jmax == 0:
  40. aab = 0.5 * (aabmin + aabmax)
  41. else:
  42. aab = aabmin + (aabmax - aabmin) * j / jmax
  43. if fast == 2:
  44. bz = fit_cubic_superfast(z0, z1, arclen, th0, th1, aab)
  45. else:
  46. bz = tocubic.fit_cubic_arclen(z0, z1, arclen, th0, th1, aab)
  47. score = tocubic.measure_bz_rk4(bz, arclen, th_fn)
  48. print '% aab =', aab, 'score =', score
  49. sys.stdout.flush()
  50. if j == 0 or score < best_score:
  51. best_score = score
  52. best_aab = aab
  53. best_bz = bz
  54. daab = .06 * (aabmax - aabmin)
  55. aabmin = max(0, best_aab - daab)
  56. aabmax = min(1, best_aab + daab)
  57. print '%--- best_aab =', best_aab
  58. return best_bz, best_score
  59. def cornu_to_cubic(t0, t1, figno):
  60. if figno == 1:
  61. aabmin = 0
  62. aabmax = 0.4
  63. elif figno == 2:
  64. aabmin = 0.5
  65. aabmax = 1.
  66. else:
  67. aabmin = 0
  68. aabmax = 1.
  69. fast = 0
  70. if figno == 3:
  71. fast = 1
  72. elif figno == 4:
  73. fast = 2
  74. def th_fn(s):
  75. return (s + t0) ** 2
  76. y0, x0 = cornu.eval_cornu(t0)
  77. y1, x1 = cornu.eval_cornu(t1)
  78. bz, score = fit_cubic((x0, y0), (x1, y1), t1 - t0, th_fn, fast, aabmin, aabmax)
  79. return bz, score
  80. def plot_k_of_bz(bz):
  81. dbz = tocubic.bz_deriv(bz)
  82. ddbz = tocubic.bz_deriv(dbz)
  83. cmd = 'moveto'
  84. ss = [0]
  85. def arclength_deriv(x, ss):
  86. dx, dy = tocubic.bz_eval(dbz, x)
  87. return [hypot(dx, dy)]
  88. dt = 0.01
  89. t = 0
  90. for i in range(101):
  91. dx, dy = tocubic.bz_eval(dbz, t)
  92. ddx, ddy = tocubic.bz_eval(ddbz, t)
  93. k = (ddy * dx - dy * ddx) / (dx * dx + dy * dy) ** 1.5
  94. print 100 + 500 * ss[0], 100 + 200 * k, cmd
  95. cmd = 'lineto'
  96. dsdx = arclength_deriv(t, ss)
  97. tocubic.rk4(ss, dsdx, t, .01, arclength_deriv)
  98. t += dt
  99. print 'stroke'
  100. def plot_k_nominal(s0, s1):
  101. k0 = 2 * s0
  102. k1 = 2 * s1
  103. print 'gsave 0.5 setlinewidth'
  104. print 100, 100 + 200 * k0, 'moveto'
  105. print 100 + 500 * (s1 - s0), 100 + 200 * k1, 'lineto'
  106. print 'stroke grestore'
  107. def simple_bez():
  108. eps_prologue(95, 126, 552, 508, 0)
  109. tocubic.plot_prolog()
  110. print '/ss 1.5 def'
  111. print '/circle { ss 0 moveto currentpoint exch ss sub exch ss 0 360 arc } bind def'
  112. bz, score = cornu_to_cubic(.5, 1.1, 2)
  113. fromcubic.plot_bzs([[bz]], (-400, 100), 1000, True)
  114. print 'stroke'
  115. print '/Times-Roman 12 selectfont'
  116. print '95 130 moveto ((x0, y0)) show'
  117. print '360 200 moveto ((x1, y1)) show'
  118. print '480 340 moveto ((x2, y2)) show'
  119. print '505 495 moveto ((x3, y3)) show'
  120. print 'showpage'
  121. eps_trailer()
  122. def fast_bez(figno):
  123. if figno == 3:
  124. y1 = 520
  125. else:
  126. y1 = 550
  127. eps_prologue(95, 140, 552, y1, 0)
  128. tocubic.plot_prolog()
  129. print '/ss 1.5 def'
  130. print '/circle { ss 0 moveto currentpoint exch ss sub exch ss 0 360 arc } bind def'
  131. bz, score = cornu_to_cubic(.5, 1.1, figno)
  132. fromcubic.plot_bzs([[bz]], (-400, 100), 1000, True)
  133. print 'stroke'
  134. plot_k_nominal(.5, 1.1)
  135. plot_k_of_bz(bz)
  136. print 'showpage'
  137. eps_trailer()
  138. def bezfig(s1):
  139. eps_prologue(95, 38, 510, 550, 0)
  140. #print '0.5 0.5 scale 500 100 translate'
  141. tocubic.plot_prolog()
  142. print '/ss 1.5 def'
  143. print '/circle { ss 0 moveto currentpoint exch ss sub exch ss 0 360 arc } bind def'
  144. bz, score = cornu_to_cubic(.5, 0.85, 1)
  145. fromcubic.plot_bzs([[bz]], (-400, 0), 1000, True)
  146. print 'stroke'
  147. plot_k_nominal(.5, 0.85)
  148. plot_k_of_bz(bz)
  149. bz, score = cornu_to_cubic(.5, 0.85, 2)
  150. fromcubic.plot_bzs([[bz]], (-400, 100), 1000, True)
  151. print 'stroke'
  152. print 'gsave 0 50 translate'
  153. plot_k_nominal(.5, .85)
  154. plot_k_of_bz(bz)
  155. print 'grestore'
  156. print 'showpage'
  157. import sys
  158. if __name__ == '__main__':
  159. figno = int(sys.argv[1])
  160. if figno == 0:
  161. simple_bez()
  162. elif figno == 1:
  163. bezfig(1.0)
  164. elif figno == 2:
  165. bezfig(0.85)
  166. else:
  167. fast_bez(figno)
  168. #fast_bez(4)