pcorn.py 5.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160
  1. # Utilities for piecewise cornu representation of curves
  2. from math import *
  3. import clothoid
  4. import cornu
  5. class Segment:
  6. def __init__(self, z0, z1, th0, th1):
  7. self.z0 = z0
  8. self.z1 = z1
  9. self.th0 = th0
  10. self.th1 = th1
  11. self.compute()
  12. def __repr__(self):
  13. return '[' + `self.z0` + `self.z1` + ' ' + `self.th0` + ' ' + `self.th1` + ']'
  14. def compute(self):
  15. dx = self.z1[0] - self.z0[0]
  16. dy = self.z1[1] - self.z0[1]
  17. chord = hypot(dy, dx)
  18. chth = atan2(dy, dx)
  19. k0, k1 = clothoid.solve_clothoid(self.th0, self.th1)
  20. charc = clothoid.compute_chord(k0, k1)
  21. self.chord = chord
  22. self.chth = chth
  23. self.k0, self.k1 = k0, k1
  24. self.charc = charc
  25. self.arclen = chord / charc
  26. self.thmid = self.chth - self.th0 + 0.5 * self.k0 - 0.125 * self.k1
  27. self.setup_xy_fresnel()
  28. def setup_xy_fresnel(self):
  29. k0, k1 = self.k0, self.k1
  30. if k1 == 0: k1 = 1e-6 # hack
  31. if k1 != 0:
  32. sqrk1 = sqrt(2 * abs(k1))
  33. t0 = (k0 - .5 * k1) / sqrk1
  34. t1 = (k0 + .5 * k1) / sqrk1
  35. (y0, x0) = cornu.eval_cornu(t0)
  36. (y1, x1) = cornu.eval_cornu(t1)
  37. chord_th = atan2(y1 - y0, x1 - x0)
  38. chord = hypot(y1 - y0, x1 - x0)
  39. scale = self.chord / chord
  40. if k1 >= 0:
  41. th = self.chth - chord_th
  42. self.mxx = scale * cos(th)
  43. self.myx = scale * sin(th)
  44. self.mxy = -self.myx
  45. self.myy = self.mxx
  46. else:
  47. th = self.chth + chord_th
  48. self.mxx = scale * cos(th)
  49. self.myx = scale * sin(th)
  50. self.mxy = self.myx
  51. self.myy = -self.mxx
  52. # rotate -chord_th, flip top/bottom, rotate self.chth
  53. self.x0 = self.z0[0] - (self.mxx * x0 + self.mxy * y0)
  54. self.y0 = self.z0[1] - (self.myx * x0 + self.myy * y0)
  55. def th(self, s):
  56. u = s / self.arclen - 0.5
  57. return self.thmid + (0.5 * self.k1 * u + self.k0) * u
  58. def xy(self, s):
  59. # using fresnel integrals; polynomial approx might be better
  60. u = s / self.arclen - 0.5
  61. k0, k1 = self.k0, self.k1
  62. if k1 == 0: k1 = 1e-6 # hack
  63. if k1 != 0:
  64. sqrk1 = sqrt(2 * abs(k1))
  65. t = (k0 + u * k1) / sqrk1
  66. (y, x) = cornu.eval_cornu(t)
  67. return [self.x0 + self.mxx * x + self.mxy * y,
  68. self.y0 + self.myx * x + self.myy * y]
  69. def find_extrema(self):
  70. # find solutions of th(s) = 0 mod pi/2
  71. # todo: find extra solutions when there's an inflection
  72. th0 = self.thmid + 0.125 * self.k1 - 0.5 * self.k0
  73. th1 = self.thmid + 0.125 * self.k1 + 0.5 * self.k0
  74. twooverpi = 2 / pi
  75. n0 = int(floor(th0 * twooverpi))
  76. n1 = int(floor(th1 * twooverpi))
  77. if th1 > th0: signum = 1
  78. else: signum = -1
  79. result = []
  80. for i in range(n0, n1, signum):
  81. th = pi/2 * (i + 0.5 * (signum + 1))
  82. a = .5 * self.k1
  83. b = self.k0
  84. c = self.thmid - th
  85. if a == 0:
  86. u1 = -c/b
  87. u2 = 1000
  88. else:
  89. sqrtdiscrim = sqrt(b * b - 4 * a * c)
  90. u1 = (-b - sqrtdiscrim) / (2 * a)
  91. u2 = (-b + sqrtdiscrim) / (2 * a)
  92. if u1 >= -0.5 and u1 < 0.5:
  93. result.append(self.arclen * (u1 + 0.5))
  94. if u2 >= -0.5 and u2 < 0.5:
  95. result.append(self.arclen * (u2 + 0.5))
  96. return result
  97. class Curve:
  98. def __init__(self, segs):
  99. self.segs = segs
  100. self.compute()
  101. def compute(self):
  102. arclen = 0
  103. sstarts = []
  104. for seg in self.segs:
  105. sstarts.append(arclen)
  106. arclen += seg.arclen
  107. self.arclen = arclen
  108. self.sstarts = sstarts
  109. def th(self, s, deltas = False):
  110. u = s / self.arclen
  111. s = self.arclen * (u - floor(u))
  112. if s == 0 and not deltas: s = self.arclen
  113. i = 0
  114. while i < len(self.segs) - 1:
  115. # binary search would make a lot of sense here
  116. snext = self.sstarts[i + 1]
  117. if s < snext or (not deltas and s == snext):
  118. break
  119. i += 1
  120. return self.segs[i].th(s - self.sstarts[i])
  121. def xy(self, s):
  122. u = s / self.arclen
  123. s = self.arclen * (u - floor(u))
  124. i = 0
  125. while i < len(self.segs) - 1:
  126. # binary search would make a lot of sense here
  127. if s <= self.sstarts[i + 1]:
  128. break
  129. i += 1
  130. return self.segs[i].xy(s - self.sstarts[i])
  131. def find_extrema(self):
  132. result = []
  133. for i in range(len(self.segs)):
  134. seg = self.segs[i]
  135. for s in seg.find_extrema():
  136. result.append(s + self.sstarts[i])
  137. return result
  138. def find_breaks(self):
  139. result = []
  140. for i in range(len(self.segs)):
  141. pseg = self.segs[(i + len(self.segs) - 1) % len(self.segs)]
  142. seg = self.segs[i]
  143. th = clothoid.mod_2pi(pseg.chth + pseg.th1 - (seg.chth - seg.th0))
  144. print '% pseg', pseg.chth + pseg.th1, 'seg', seg.chth - seg.th0
  145. pisline = pseg.k0 == 0 and pseg.k1 == 0
  146. sisline = seg.k0 == 0 and seg.k1 == 0
  147. if fabs(th) > 1e-3 or (pisline and not sisline) or (sisline and not pisline):
  148. result.append(self.sstarts[i])
  149. return result