scikit-image

1
import numpy as np
2
from scipy import signal
3

4

5
def approximate_polygon(coords, tolerance):
6
    """Approximate a polygonal chain with the specified tolerance.
7

8
    It is based on the Douglas-Peucker algorithm.
9

10
    Note that the approximated polygon is always within the convex hull of the
11
    original polygon.
12

13
    Parameters
14
    ----------
15
    coords : (K, 2) array
16
        Coordinate array.
17
    tolerance : float
18
        Maximum distance from original points of polygon to approximated
19
        polygonal chain. If tolerance is 0, the original coordinate array
20
        is returned.
21

22
    Returns
23
    -------
24
    coords : (L, 2) array
25
        Approximated polygonal chain where L <= K.
26

27
    References
28
    ----------
29
    .. [1] https://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm
30
    """
31
    if tolerance <= 0:
32
        return coords
33

34
    chain = np.zeros(coords.shape[0], 'bool')
35
    # pre-allocate distance array for all points
36
    dists = np.zeros(coords.shape[0])
37
    chain[0] = True
38
    chain[-1] = True
39
    pos_stack = [(0, chain.shape[0] - 1)]
40
    end_of_chain = False
41

42
    while not end_of_chain:
43
        start, end = pos_stack.pop()
44
        # determine properties of current line segment
45
        r0, c0 = coords[start, :]
46
        r1, c1 = coords[end, :]
47
        dr = r1 - r0
48
        dc = c1 - c0
49
        segment_angle = -np.arctan2(dr, dc)
50
        segment_dist = c0 * np.sin(segment_angle) + r0 * np.cos(segment_angle)
51

52
        # select points in-between line segment
53
        segment_coords = coords[start + 1 : end, :]
54
        segment_dists = dists[start + 1 : end]
55

56
        # check whether to take perpendicular or euclidean distance with
57
        # inner product of vectors
58

59
        # vectors from points -> start and end
60
        dr0 = segment_coords[:, 0] - r0
61
        dc0 = segment_coords[:, 1] - c0
62
        dr1 = segment_coords[:, 0] - r1
63
        dc1 = segment_coords[:, 1] - c1
64
        # vectors points -> start and end projected on start -> end vector
65
        projected_lengths0 = dr0 * dr + dc0 * dc
66
        projected_lengths1 = -dr1 * dr - dc1 * dc
67
        perp = np.logical_and(projected_lengths0 > 0, projected_lengths1 > 0)
68
        eucl = np.logical_not(perp)
69
        segment_dists[perp] = np.abs(
70
            segment_coords[perp, 0] * np.cos(segment_angle)
71
            + segment_coords[perp, 1] * np.sin(segment_angle)
72
            - segment_dist
73
        )
74
        segment_dists[eucl] = np.minimum(
75
            # distance to start point
76
            np.sqrt(dc0[eucl] ** 2 + dr0[eucl] ** 2),
77
            # distance to end point
78
            np.sqrt(dc1[eucl] ** 2 + dr1[eucl] ** 2),
79
        )
80

81
        if np.any(segment_dists > tolerance):
82
            # select point with maximum distance to line
83
            new_end = start + np.argmax(segment_dists) + 1
84
            pos_stack.append((new_end, end))
85
            pos_stack.append((start, new_end))
86
            chain[new_end] = True
87

88
        if len(pos_stack) == 0:
89
            end_of_chain = True
90

91
    return coords[chain, :]
92

93

94
# B-Spline subdivision
95
_SUBDIVISION_MASKS = {
96
    # degree: (mask_even, mask_odd)
97
    #         extracted from (degree + 2)th row of Pascal's triangle
98
    1: ([1, 1], [1, 1]),
99
    2: ([3, 1], [1, 3]),
100
    3: ([1, 6, 1], [0, 4, 4]),
101
    4: ([5, 10, 1], [1, 10, 5]),
102
    5: ([1, 15, 15, 1], [0, 6, 20, 6]),
103
    6: ([7, 35, 21, 1], [1, 21, 35, 7]),
104
    7: ([1, 28, 70, 28, 1], [0, 8, 56, 56, 8]),
105
}
106

107

108
def subdivide_polygon(coords, degree=2, preserve_ends=False):
109
    """Subdivision of polygonal curves using B-Splines.
110

111
    Note that the resulting curve is always within the convex hull of the
112
    original polygon. Circular polygons stay closed after subdivision.
113

114
    Parameters
115
    ----------
116
    coords : (K, 2) array
117
        Coordinate array.
118
    degree : {1, 2, 3, 4, 5, 6, 7}, optional
119
        Degree of B-Spline. Default is 2.
120
    preserve_ends : bool, optional
121
        Preserve first and last coordinate of non-circular polygon. Default is
122
        False.
123

124
    Returns
125
    -------
126
    coords : (L, 2) array
127
        Subdivided coordinate array.
128

129
    References
130
    ----------
131
    .. [1] http://mrl.nyu.edu/publications/subdiv-course2000/coursenotes00.pdf
132
    """
133
    if degree not in _SUBDIVISION_MASKS:
134
        raise ValueError("Invalid B-Spline degree. Only degree 1 - 7 is " "supported.")
135

136
    circular = np.all(coords[0, :] == coords[-1, :])
137

138
    method = 'valid'
139
    if circular:
140
        # remove last coordinate because of wrapping
141
        coords = coords[:-1, :]
142
        # circular convolution by wrapping boundaries
143
        method = 'same'
144

145
    mask_even, mask_odd = _SUBDIVISION_MASKS[degree]
146
    # divide by total weight
147
    mask_even = np.array(mask_even, float) / (2**degree)
148
    mask_odd = np.array(mask_odd, float) / (2**degree)
149

150
    even = signal.convolve2d(
151
        coords.T, np.atleast_2d(mask_even), mode=method, boundary='wrap'
152
    )
153
    odd = signal.convolve2d(
154
        coords.T, np.atleast_2d(mask_odd), mode=method, boundary='wrap'
155
    )
156

157
    out = np.zeros((even.shape[1] + odd.shape[1], 2))
158
    out[1::2] = even.T
159
    out[::2] = odd.T
160

161
    if circular:
162
        # close polygon
163
        out = np.vstack([out, out[0, :]])
164

165
    if preserve_ends and not circular:
166
        out = np.vstack([coords[0, :], out, coords[-1, :]])
167

168
    return out
169