_polygon.py 5.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168
  1. import numpy as np
  2. from scipy import signal
  3. def approximate_polygon(coords, tolerance):
  4. """Approximate a polygonal chain with the specified tolerance.
  5. It is based on the Douglas-Peucker algorithm.
  6. Note that the approximated polygon is always within the convex hull of the
  7. original polygon.
  8. Parameters
  9. ----------
  10. coords : (K, 2) array
  11. Coordinate array.
  12. tolerance : float
  13. Maximum distance from original points of polygon to approximated
  14. polygonal chain. If tolerance is 0, the original coordinate array
  15. is returned.
  16. Returns
  17. -------
  18. coords : (L, 2) array
  19. Approximated polygonal chain where L <= K.
  20. References
  21. ----------
  22. .. [1] https://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm
  23. """
  24. if tolerance <= 0:
  25. return coords
  26. chain = np.zeros(coords.shape[0], 'bool')
  27. # pre-allocate distance array for all points
  28. dists = np.zeros(coords.shape[0])
  29. chain[0] = True
  30. chain[-1] = True
  31. pos_stack = [(0, chain.shape[0] - 1)]
  32. end_of_chain = False
  33. while not end_of_chain:
  34. start, end = pos_stack.pop()
  35. # determine properties of current line segment
  36. r0, c0 = coords[start, :]
  37. r1, c1 = coords[end, :]
  38. dr = r1 - r0
  39. dc = c1 - c0
  40. segment_angle = -np.arctan2(dr, dc)
  41. segment_dist = c0 * np.sin(segment_angle) + r0 * np.cos(segment_angle)
  42. # select points in-between line segment
  43. segment_coords = coords[start + 1 : end, :]
  44. segment_dists = dists[start + 1 : end]
  45. # check whether to take perpendicular or euclidean distance with
  46. # inner product of vectors
  47. # vectors from points -> start and end
  48. dr0 = segment_coords[:, 0] - r0
  49. dc0 = segment_coords[:, 1] - c0
  50. dr1 = segment_coords[:, 0] - r1
  51. dc1 = segment_coords[:, 1] - c1
  52. # vectors points -> start and end projected on start -> end vector
  53. projected_lengths0 = dr0 * dr + dc0 * dc
  54. projected_lengths1 = -dr1 * dr - dc1 * dc
  55. perp = np.logical_and(projected_lengths0 > 0, projected_lengths1 > 0)
  56. eucl = np.logical_not(perp)
  57. segment_dists[perp] = np.abs(
  58. segment_coords[perp, 0] * np.cos(segment_angle)
  59. + segment_coords[perp, 1] * np.sin(segment_angle)
  60. - segment_dist
  61. )
  62. segment_dists[eucl] = np.minimum(
  63. # distance to start point
  64. np.sqrt(dc0[eucl] ** 2 + dr0[eucl] ** 2),
  65. # distance to end point
  66. np.sqrt(dc1[eucl] ** 2 + dr1[eucl] ** 2),
  67. )
  68. if np.any(segment_dists > tolerance):
  69. # select point with maximum distance to line
  70. new_end = start + np.argmax(segment_dists) + 1
  71. pos_stack.append((new_end, end))
  72. pos_stack.append((start, new_end))
  73. chain[new_end] = True
  74. if len(pos_stack) == 0:
  75. end_of_chain = True
  76. return coords[chain, :]
  77. # B-Spline subdivision
  78. _SUBDIVISION_MASKS = {
  79. # degree: (mask_even, mask_odd)
  80. # extracted from (degree + 2)th row of Pascal's triangle
  81. 1: ([1, 1], [1, 1]),
  82. 2: ([3, 1], [1, 3]),
  83. 3: ([1, 6, 1], [0, 4, 4]),
  84. 4: ([5, 10, 1], [1, 10, 5]),
  85. 5: ([1, 15, 15, 1], [0, 6, 20, 6]),
  86. 6: ([7, 35, 21, 1], [1, 21, 35, 7]),
  87. 7: ([1, 28, 70, 28, 1], [0, 8, 56, 56, 8]),
  88. }
  89. def subdivide_polygon(coords, degree=2, preserve_ends=False):
  90. """Subdivision of polygonal curves using B-Splines.
  91. Note that the resulting curve is always within the convex hull of the
  92. original polygon. Circular polygons stay closed after subdivision.
  93. Parameters
  94. ----------
  95. coords : (K, 2) array
  96. Coordinate array.
  97. degree : {1, 2, 3, 4, 5, 6, 7}, optional
  98. Degree of B-Spline. Default is 2.
  99. preserve_ends : bool, optional
  100. Preserve first and last coordinate of non-circular polygon. Default is
  101. False.
  102. Returns
  103. -------
  104. coords : (L, 2) array
  105. Subdivided coordinate array.
  106. References
  107. ----------
  108. .. [1] http://mrl.nyu.edu/publications/subdiv-course2000/coursenotes00.pdf
  109. """
  110. if degree not in _SUBDIVISION_MASKS:
  111. raise ValueError("Invalid B-Spline degree. Only degree 1 - 7 is " "supported.")
  112. circular = np.all(coords[0, :] == coords[-1, :])
  113. method = 'valid'
  114. if circular:
  115. # remove last coordinate because of wrapping
  116. coords = coords[:-1, :]
  117. # circular convolution by wrapping boundaries
  118. method = 'same'
  119. mask_even, mask_odd = _SUBDIVISION_MASKS[degree]
  120. # divide by total weight
  121. mask_even = np.array(mask_even, float) / (2**degree)
  122. mask_odd = np.array(mask_odd, float) / (2**degree)
  123. even = signal.convolve2d(
  124. coords.T, np.atleast_2d(mask_even), mode=method, boundary='wrap'
  125. )
  126. odd = signal.convolve2d(
  127. coords.T, np.atleast_2d(mask_odd), mode=method, boundary='wrap'
  128. )
  129. out = np.zeros((even.shape[1] + odd.shape[1], 2))
  130. out[1::2] = even.T
  131. out[::2] = odd.T
  132. if circular:
  133. # close polygon
  134. out = np.vstack([out, out[0, :]])
  135. if preserve_ends and not circular:
  136. out = np.vstack([coords[0, :], out, coords[-1, :]])
  137. return out