scikit-image

1
import numpy as np
2

3
from ._find_contours_cy import _get_contour_segments
4

5
from collections import deque
6

7
_param_options = ('high', 'low')
8

9

10
def find_contours(
11
    image, level=None, fully_connected='low', positive_orientation='low', *, mask=None
12
):
13
    """Find iso-valued contours in a 2D array for a given level value.
14

15
    Uses the "marching squares" method to compute the iso-valued contours of
16
    the input 2D array for a particular level value. Array values are linearly
17
    interpolated to provide better precision for the output contours.
18

19
    Parameters
20
    ----------
21
    image : (M, N) ndarray of double
22
        Input image in which to find contours.
23
    level : float, optional
24
        Value along which to find contours in the array. By default, the level
25
        is set to (max(image) + min(image)) / 2
26

27
        .. versionchanged:: 0.18
28
            This parameter is now optional.
29
    fully_connected : str, {'low', 'high'}
30
         Indicates whether array elements below the given level value are to be
31
         considered fully-connected (and hence elements above the value will
32
         only be face connected), or vice-versa. (See notes below for details.)
33
    positive_orientation : str, {'low', 'high'}
34
         Indicates whether the output contours will produce positively-oriented
35
         polygons around islands of low- or high-valued elements. If 'low' then
36
         contours will wind counter-clockwise around elements below the
37
         iso-value. Alternately, this means that low-valued elements are always
38
         on the left of the contour. (See below for details.)
39
    mask : (M, N) ndarray of bool or None
40
        A boolean mask, True where we want to draw contours.
41
        Note that NaN values are always excluded from the considered region
42
        (``mask`` is set to ``False`` wherever ``array`` is ``NaN``).
43

44
    Returns
45
    -------
46
    contours : list of (K, 2) ndarrays
47
        Each contour is a ndarray of ``(row, column)`` coordinates along the contour.
48

49
    See Also
50
    --------
51
    skimage.measure.marching_cubes
52

53
    Notes
54
    -----
55
    The marching squares algorithm is a special case of the marching cubes
56
    algorithm [1]_.  A simple explanation is available here:
57

58
    https://users.polytech.unice.fr/~lingrand/MarchingCubes/algo.html
59

60
    There is a single ambiguous case in the marching squares algorithm: when
61
    a given ``2 x 2``-element square has two high-valued and two low-valued
62
    elements, each pair diagonally adjacent. (Where high- and low-valued is
63
    with respect to the contour value sought.) In this case, either the
64
    high-valued elements can be 'connected together' via a thin isthmus that
65
    separates the low-valued elements, or vice-versa. When elements are
66
    connected together across a diagonal, they are considered 'fully
67
    connected' (also known as 'face+vertex-connected' or '8-connected'). Only
68
    high-valued or low-valued elements can be fully-connected, the other set
69
    will be considered as 'face-connected' or '4-connected'. By default,
70
    low-valued elements are considered fully-connected; this can be altered
71
    with the 'fully_connected' parameter.
72

73
    Output contours are not guaranteed to be closed: contours which intersect
74
    the array edge or a masked-off region (either where mask is False or where
75
    array is NaN) will be left open. All other contours will be closed. (The
76
    closed-ness of a contours can be tested by checking whether the beginning
77
    point is the same as the end point.)
78

79
    Contours are oriented. By default, array values lower than the contour
80
    value are to the left of the contour and values greater than the contour
81
    value are to the right. This means that contours will wind
82
    counter-clockwise (i.e. in 'positive orientation') around islands of
83
    low-valued pixels. This behavior can be altered with the
84
    'positive_orientation' parameter.
85

86
    The order of the contours in the output list is determined by the position
87
    of the smallest ``x,y`` (in lexicographical order) coordinate in the
88
    contour.  This is a side effect of how the input array is traversed, but
89
    can be relied upon.
90

91
    .. warning::
92

93
       Array coordinates/values are assumed to refer to the *center* of the
94
       array element. Take a simple example input: ``[0, 1]``. The interpolated
95
       position of 0.5 in this array is midway between the 0-element (at
96
       ``x=0``) and the 1-element (at ``x=1``), and thus would fall at
97
       ``x=0.5``.
98

99
    This means that to find reasonable contours, it is best to find contours
100
    midway between the expected "light" and "dark" values. In particular,
101
    given a binarized array, *do not* choose to find contours at the low or
102
    high value of the array. This will often yield degenerate contours,
103
    especially around structures that are a single array element wide. Instead,
104
    choose a middle value, as above.
105

106
    References
107
    ----------
108
    .. [1] Lorensen, William and Harvey E. Cline. Marching Cubes: A High
109
           Resolution 3D Surface Construction Algorithm. Computer Graphics
110
           (SIGGRAPH 87 Proceedings) 21(4) July 1987, p. 163-170).
111
           :DOI:`10.1145/37401.37422`
112

113
    Examples
114
    --------
115
    >>> a = np.zeros((3, 3))
116
    >>> a[0, 0] = 1
117
    >>> a
118
    array([[1., 0., 0.],
119
           [0., 0., 0.],
120
           [0., 0., 0.]])
121
    >>> find_contours(a, 0.5)
122
    [array([[0. , 0.5],
123
           [0.5, 0. ]])]
124
    """
125
    if fully_connected not in _param_options:
126
        raise ValueError(
127
            'Parameters "fully_connected" must be either ' '"high" or "low".'
128
        )
129
    if positive_orientation not in _param_options:
130
        raise ValueError(
131
            'Parameters "positive_orientation" must be either ' '"high" or "low".'
132
        )
133
    if image.shape[0] < 2 or image.shape[1] < 2:
134
        raise ValueError("Input array must be at least 2x2.")
135
    if image.ndim != 2:
136
        raise ValueError('Only 2D arrays are supported.')
137
    if mask is not None:
138
        if mask.shape != image.shape:
139
            raise ValueError('Parameters "array" and "mask"' ' must have same shape.')
140
        if not np.can_cast(mask.dtype, bool, casting='safe'):
141
            raise TypeError('Parameter "mask" must be a binary array.')
142
        mask = mask.astype(np.uint8, copy=False)
143
    if level is None:
144
        level = (np.nanmin(image) + np.nanmax(image)) / 2.0
145

146
    segments = _get_contour_segments(
147
        image.astype(np.float64), float(level), fully_connected == 'high', mask=mask
148
    )
149
    contours = _assemble_contours(segments)
150
    if positive_orientation == 'high':
151
        contours = [c[::-1] for c in contours]
152
    return contours
153

154

155
def _assemble_contours(segments):
156
    current_index = 0
157
    contours = {}
158
    starts = {}
159
    ends = {}
160
    for from_point, to_point in segments:
161
        # Ignore degenerate segments.
162
        # This happens when (and only when) one vertex of the square is
163
        # exactly the contour level, and the rest are above or below.
164
        # This degenerate vertex will be picked up later by neighboring
165
        # squares.
166
        if from_point == to_point:
167
            continue
168

169
        tail, tail_num = starts.pop(to_point, (None, None))
170
        head, head_num = ends.pop(from_point, (None, None))
171

172
        if tail is not None and head is not None:
173
            # We need to connect these two contours.
174
            if tail is head:
175
                # We need to closed a contour: add the end point
176
                head.append(to_point)
177
            else:  # tail is not head
178
                # We need to join two distinct contours.
179
                # We want to keep the first contour segment created, so that
180
                # the final contours are ordered left->right, top->bottom.
181
                if tail_num > head_num:
182
                    # tail was created second. Append tail to head.
183
                    head.extend(tail)
184
                    # Remove tail from the detected contours
185
                    contours.pop(tail_num, None)
186
                    # Update starts and ends
187
                    starts[head[0]] = (head, head_num)
188
                    ends[head[-1]] = (head, head_num)
189
                else:  # tail_num <= head_num
190
                    # head was created second. Prepend head to tail.
191
                    tail.extendleft(reversed(head))
192
                    # Remove head from the detected contours
193
                    starts.pop(head[0], None)  # head[0] can be == to_point!
194
                    contours.pop(head_num, None)
195
                    # Update starts and ends
196
                    starts[tail[0]] = (tail, tail_num)
197
                    ends[tail[-1]] = (tail, tail_num)
198
        elif tail is None and head is None:
199
            # We need to add a new contour
200
            new_contour = deque((from_point, to_point))
201
            contours[current_index] = new_contour
202
            starts[from_point] = (new_contour, current_index)
203
            ends[to_point] = (new_contour, current_index)
204
            current_index += 1
205
        elif head is None:  # tail is not None
206
            # tail first element is to_point: the new segment should be
207
            # prepended.
208
            tail.appendleft(from_point)
209
            # Update starts
210
            starts[from_point] = (tail, tail_num)
211
        else:  # tail is None and head is not None:
212
            # head last element is from_point: the new segment should be
213
            # appended
214
            head.append(to_point)
215
            # Update ends
216
            ends[to_point] = (head, head_num)
217

218
    return [np.array(contour) for _, contour in sorted(contours.items())]
219