bezier.py 6.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195
  1. # Purpose: general bezier curve
  2. # Created: 26.03.2010
  3. # Copyright (c) 2010-2018 Manfred Moitzi
  4. # License: MIT License
  5. from typing import TYPE_CHECKING, List, Iterable, Tuple, Dict
  6. from ezdxf.math.vector import Vector
  7. if TYPE_CHECKING:
  8. from ezdxf.eztypes import Vertex
  9. """
  10. Bezier curves
  11. =============
  12. https://www.cl.cam.ac.uk/teaching/2000/AGraphHCI/SMEG/node3.html
  13. A Bezier curve is a weighted sum of n+1 control points, P0, P1, ..., Pn, where the weights are the Bernstein
  14. polynomials.
  15. The Bezier curve of order n+1 (degree n) has n+1 control points. These are the first three orders of Bezier curve
  16. definitions.
  17. (75) linear P(t) = (1-t)*P0 + t*P1
  18. (76) quadratic P(t) = (1-t)^2*P0 + 2*(t-1)*t*P1 + t^2*P2
  19. (77) cubic P(t) = (1-t)^3*P0 + 3*(t-1)^2*t*P1 + 3*(t-1)*t^2*P2 + t^3*P3
  20. Ways of thinking about Bezier curves
  21. ------------------------------------
  22. There are several useful ways in which you can think about Bezier curves. Here are the ones that I use.
  23. Linear interpolation
  24. ~~~~~~~~~~~~~~~~~~~~
  25. Equation (75) is obviously a linear interpolation between two points. Equation (76) can be rewritten as a linear
  26. interpolation between linear interpolations between points.
  27. Weighted average
  28. ~~~~~~~~~~~~~~~~
  29. A Bezier curve can be seen as a weighted average of all of its control points. Because all of the weights are
  30. positive, and because the weights sum to one, the Bezier curve is guaranteed to lie within the convex hull of its
  31. control points.
  32. Refinement of the control polygon
  33. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  34. A Bezier curve can be seen as some sort of refinement of the polygon made by connecting its control points in order.
  35. The Bezier curve starts and ends at the two end points and its shape is determined by the relative positions of the
  36. n-1 other control points, although it will generally not pass through these other control points. The tangent vectors
  37. at the start and end of the curve pass through the end point and the immediately adjacent point.
  38. Continuity
  39. ----------
  40. You should note that each Bezier curve is independent of any other Bezier curve. If we wish two Bezier curves to join
  41. with any type of continuity, then we must explicitly position the control points of the second curve so that they bear
  42. the appropriate relationship with the control points in the first curve.
  43. Any Bezier curve is infinitely differentiable within itself, and is therefore continuous to any degree.
  44. """
  45. class Bezier:
  46. """
  47. A general Bezier curve.
  48. Accepts 2d points as definition points, but output ist always 3d (z-axis is 0).
  49. """
  50. def __init__(self, defpoints: Iterable['Vertex']):
  51. self._defpoints = [Vector(p) for p in defpoints] # type: List[Vector]
  52. @property
  53. def control_points(self) -> List[Vector]:
  54. return self._defpoints
  55. def approximate(self, segments: int = 20) -> Iterable[Vector]:
  56. step = 1.0 / float(segments)
  57. for point_index in range(segments + 1):
  58. yield self.point(point_index * step)
  59. def point(self, t: float) -> Vector:
  60. """
  61. Returns point at BezierCurve(t) as tuple (x, y, z)
  62. Args:
  63. t: parameter in range [0, 1]
  64. Returns: Vector(x, y, z)
  65. """
  66. if t < 0. or t > 1.:
  67. raise ValueError('parameter t in range [0, 1]')
  68. if (1.0 - t) < 5e-6:
  69. t = 1.0
  70. point = [0., 0., 0.]
  71. defpoints = self._defpoints
  72. len_defpoints = len(defpoints)
  73. for axis in (0, 1, 2):
  74. for i in range(len_defpoints):
  75. bsf = bernstein_basis(len_defpoints - 1, i, t)
  76. point[axis] += bsf * defpoints[i][axis]
  77. return Vector(point)
  78. class DBezier(Bezier):
  79. """
  80. Calculate the Points and Derivative of a Bezier curve.
  81. """
  82. def point(self, t: float) -> Tuple[Vector, Vector, Vector]:
  83. """
  84. Returns (point, derivative1, derivative2) at BezierCurve(t)
  85. Args:
  86. t: parameter in range [0, 1]
  87. Returns: (point, derivative1, derivative2)
  88. point -- Vector(x, y, z)
  89. derivative1 -- Vector(x, y, z)
  90. derivative2 -- Vector(x, y, z)
  91. """
  92. if t < 0. or t > 1.:
  93. raise ValueError('parameter t in range [0, 1]')
  94. if (1.0 - t) < 5e-6:
  95. t = 1.0
  96. defpoints = self._defpoints
  97. npts = len(defpoints)
  98. npts0 = npts - 1
  99. point = [0., 0., 0.]
  100. d1 = [0., 0., 0.]
  101. d2 = [0., 0., 0.]
  102. for axis in (0, 1, 2):
  103. if t == 0.0:
  104. d1[axis] = npts0 * (defpoints[1][axis] - defpoints[0][axis])
  105. d2[axis] = npts0 * (npts0 - 1) * (defpoints[0][axis] - 2. * defpoints[1][axis] + defpoints[2][axis])
  106. for i in range(len(defpoints)):
  107. tempbasis = bernstein_basis(npts0, i, t)
  108. point[axis] += tempbasis * defpoints[i][axis]
  109. if 0.0 < t < 1.0:
  110. d1[axis] += ((i - npts0 * t) / (t * (1. - t))) * tempbasis * defpoints[i][axis]
  111. temp1 = (i - npts0 * t) ** 2
  112. temp2 = temp1 - npts0 * t ** 2 - i * (1. - 2. * t)
  113. d2[axis] += (temp2 / (t ** 2 * (1. - t) ** 2)) * tempbasis * defpoints[i][axis]
  114. if t == 1.0:
  115. d1[axis] = npts0 * (defpoints[npts0][axis] - defpoints[npts0 - 1][axis])
  116. d2[axis] = npts0 * (npts0 - 1) * (defpoints[npts0][axis] - 2 * defpoints[npts0 - 1][axis] +
  117. defpoints[npts0 - 2][axis])
  118. return Vector(point), Vector(d1), Vector(d2)
  119. def bernstein_basis(n: int, i: int, t: float) -> float:
  120. # handle the special cases to avoid domain problem with pow
  121. if t == 0.0 and i == 0:
  122. ti = 1.0
  123. else:
  124. ti = pow(t, i)
  125. if n == i and t == 1.0:
  126. tni = 1.0
  127. else:
  128. tni = pow((1.0 - t), (n - i))
  129. Ni = factorial(n) / (factorial(i) * factorial(n - i))
  130. return Ni * ti * tni
  131. class _Factorial:
  132. _values = {0: 1.0} # type: Dict[int, float]
  133. def __init__(self, maxvalue: int = 33):
  134. value = 1.
  135. for i in range(1, maxvalue):
  136. value *= i
  137. self._values[i] = value
  138. def get(self, value: int) -> float:
  139. value = int(value)
  140. try:
  141. return self._values[value]
  142. except KeyError:
  143. result = self.get(value - 1) * value
  144. self._values[value] = result
  145. return result
  146. _factorial = _Factorial()
  147. factorial = _factorial.get