fromcubic.py 6.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208
  1. # Convert piecewise cubic into piecewise clothoid representation.
  2. from math import *
  3. import clothoid
  4. import pcorn
  5. import tocubic
  6. import offset
  7. def read_bz(f):
  8. result = []
  9. for l in f.xreadlines():
  10. s = l.split()
  11. if len(s) > 0:
  12. cmd = s[-1]
  13. #print s[:-1], cmd
  14. if cmd == 'm':
  15. sp = []
  16. result.append(sp)
  17. curpt = [float(x) for x in s[0:2]]
  18. startpt = curpt
  19. elif cmd == 'l':
  20. newpt = [float(x) for x in s[0:2]]
  21. sp.append((curpt, newpt))
  22. curpt = newpt
  23. elif cmd == 'c':
  24. c1 = [float(x) for x in s[0:2]]
  25. c2 = [float(x) for x in s[2:4]]
  26. newpt = [float(x) for x in s[4:6]]
  27. sp.append((curpt, c1, c2, newpt))
  28. curpt = newpt
  29. return result
  30. def plot_bzs(bzs, z0, scale, fancy = False):
  31. for sp in bzs:
  32. for i in range(len(sp)):
  33. bz = sp[i]
  34. tocubic.plot_bz(bz, z0, scale, i == 0)
  35. print 'stroke'
  36. if fancy:
  37. for i in range(len(sp)):
  38. bz = sp[i]
  39. x0, y0 = z0[0] + scale * bz[0][0], z0[1] + scale * bz[0][1]
  40. print 'gsave', x0, y0, 'translate circle fill grestore'
  41. if len(bz) == 4:
  42. x1, y1 = z0[0] + scale * bz[1][0], z0[1] + scale * bz[1][1]
  43. x2, y2 = z0[0] + scale * bz[2][0], z0[1] + scale * bz[2][1]
  44. x3, y3 = z0[0] + scale * bz[3][0], z0[1] + scale * bz[3][1]
  45. print 'gsave 0.5 setlinewidth', x0, y0, 'moveto'
  46. print x1, y1, 'lineto stroke'
  47. print x2, y2, 'moveto'
  48. print x3, y3, 'lineto stroke grestore'
  49. print 'gsave', x1, y1, 'translate 0.75 dup scale circle fill grestore'
  50. print 'gsave', x2, y2, 'translate 0.75 dup scale circle fill grestore'
  51. print 'gsave', x3, y3, 'translate 0.75 dup scale circle fill grestore'
  52. def measure_bz_cloth(seg, bz, n = 100):
  53. bz_arclen = tocubic.bz_arclength_rk4(bz)
  54. arclen_ratio = seg.arclen / bz_arclen
  55. dbz = tocubic.bz_deriv(bz)
  56. def measure_derivs(x, ys):
  57. dx, dy = tocubic.bz_eval(dbz, x)
  58. ds = hypot(dx, dy)
  59. s = ys[0] * arclen_ratio
  60. dscore = ds * (tocubic.mod_2pi(atan2(dy, dx) - seg.th(s)) ** 2)
  61. #print s, atan2(dy, dx), seg.th(s)
  62. return [ds, dscore]
  63. dt = 1./n
  64. t = 0
  65. ys = [0, 0]
  66. for i in range(n):
  67. dydx = measure_derivs(t, ys)
  68. tocubic.rk4(ys, dydx, t, dt, measure_derivs)
  69. t += dt
  70. return ys[1]
  71. def cubic_bz_to_pcorn(bz, thresh):
  72. dx = bz[3][0] - bz[0][0]
  73. dy = bz[3][1] - bz[0][1]
  74. dx1 = bz[1][0] - bz[0][0]
  75. dy1 = bz[1][1] - bz[0][1]
  76. dx2 = bz[3][0] - bz[2][0]
  77. dy2 = bz[3][1] - bz[2][1]
  78. chth = atan2(dy, dx)
  79. th0 = tocubic.mod_2pi(chth - atan2(dy1, dx1))
  80. th1 = tocubic.mod_2pi(atan2(dy2, dx2) - chth)
  81. seg = pcorn.Segment(bz[0], bz[3], th0, th1)
  82. err = measure_bz_cloth(seg, bz)
  83. if err < thresh:
  84. return [seg]
  85. else:
  86. # de Casteljau
  87. x01, y01 = 0.5 * (bz[0][0] + bz[1][0]), 0.5 * (bz[0][1] + bz[1][1])
  88. x12, y12 = 0.5 * (bz[1][0] + bz[2][0]), 0.5 * (bz[1][1] + bz[2][1])
  89. x23, y23 = 0.5 * (bz[2][0] + bz[3][0]), 0.5 * (bz[2][1] + bz[3][1])
  90. xl2, yl2 = 0.5 * (x01 + x12), 0.5 * (y01 + y12)
  91. xr1, yr1 = 0.5 * (x12 + x23), 0.5 * (y12 + y23)
  92. xm, ym = 0.5 * (xl2 + xr1), 0.5 * (yl2 + yr1)
  93. bzl = [bz[0], (x01, y01), (xl2, yl2), (xm, ym)]
  94. bzr = [(xm, ym), (xr1, yr1), (x23, y23), bz[3]]
  95. segs = cubic_bz_to_pcorn(bzl, 0.5 * thresh)
  96. segs.extend(cubic_bz_to_pcorn(bzr, 0.5 * thresh))
  97. return segs
  98. def bzs_to_pcorn(bzs, thresh = 1e-9):
  99. result = []
  100. for sp in bzs:
  101. rsp = []
  102. for bz in sp:
  103. if len(bz) == 2:
  104. dx = bz[1][0] - bz[0][0]
  105. dy = bz[1][1] - bz[0][1]
  106. th = atan2(dy, dx)
  107. rsp.append(pcorn.Segment(bz[0], bz[1], 0, 0))
  108. else:
  109. rsp.extend(cubic_bz_to_pcorn(bz, thresh))
  110. result.append(rsp)
  111. return result
  112. def plot_segs(segs):
  113. for i in range(len(segs)):
  114. seg = segs[i]
  115. if i == 0:
  116. print seg.z0[0], seg.z0[1], 'moveto'
  117. print seg.z1[0], seg.z1[1], 'lineto'
  118. print 'stroke'
  119. for i in range(len(segs)):
  120. seg = segs[i]
  121. if i == 0:
  122. print 'gsave', seg.z0[0], seg.z0[1], 'translate circle fill grestore'
  123. print 'gsave', seg.z1[0], seg.z1[1], 'translate circle fill grestore'
  124. import sys
  125. def test_to_pcorn():
  126. C1 = 0.55228
  127. bz = [(100, 100), (100 + 400 * C1, 100), (500, 500 - 400 * C1), (500, 500)]
  128. for i in range(0, 13):
  129. thresh = .1 ** i
  130. segs = cubic_bz_to_pcorn(bz, thresh)
  131. plot_segs(segs)
  132. print >> sys.stderr, thresh, len(segs)
  133. print '0 20 translate'
  134. if __name__ == '__main__':
  135. f = file(sys.argv[1])
  136. bzs = read_bz(f)
  137. rsps = bzs_to_pcorn(bzs, 1)
  138. #print rsps
  139. tocubic.plot_prolog()
  140. print 'grestore'
  141. print '1 -1 scale 0 -720 translate'
  142. print '/ss 1.5 def'
  143. print '/circle { ss 0 moveto currentpoint exch ss sub exch ss 0 360 arc } bind def'
  144. tot = 0
  145. for segs in rsps:
  146. curve = pcorn.Curve(segs)
  147. #curve = offset.offset(curve, 10)
  148. print '%', curve.arclen
  149. print '%', curve.sstarts
  150. if 0:
  151. print 'gsave 1 0 0 setrgbcolor'
  152. cmd = 'moveto'
  153. for i in range(100):
  154. s = i * .01 * curve.arclen
  155. x, y = curve.xy(s)
  156. th = curve.th(s)
  157. sth = 5 * sin(th)
  158. cth = 5 * cos(th)
  159. print x, y, cmd
  160. cmd = 'lineto'
  161. print 'closepath stroke grestore'
  162. for i in range(100):
  163. s = i * .01 * curve.arclen
  164. x, y = curve.xy(s)
  165. th = curve.th(s)
  166. sth = 5 * sin(th)
  167. cth = 5 * cos(th)
  168. if 0:
  169. print x - cth, y - sth, 'moveto'
  170. print x + cth, y + sth, 'lineto stroke'
  171. if 1:
  172. for s in curve.find_breaks():
  173. print 'gsave 0 1 0 setrgbcolor'
  174. x, y = curve.xy(s)
  175. print x, y, 'translate 2 dup scale circle fill'
  176. print 'grestore'
  177. #plot_segs(segs)
  178. print 'gsave 0 0 0 setrgbcolor'
  179. optim = 3
  180. thresh = 1e-2
  181. new_bzs = tocubic.pcorn_curve_to_bzs(curve, optim, thresh)
  182. tot += len(new_bzs)
  183. plot_bzs([new_bzs], (0, 0), 1, True)
  184. print 'grestore'
  185. print 'grestore'
  186. print '/Helvetica 12 selectfont'
  187. print '36 720 moveto (thresh=%g optim=%d) show' % (thresh, optim)
  188. print '( tot segs=%d) show' % tot
  189. print 'showpage'
  190. #plot_bzs(bzs, (100, 100), 1)