forked from scikit-image/scikit-image
-
Notifications
You must be signed in to change notification settings - Fork 0
/
_canny.py
348 lines (296 loc) · 12.8 KB
/
_canny.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
"""
canny.py - Canny Edge detector
Reference: Canny, J., A Computational Approach To Edge Detection, IEEE Trans.
Pattern Analysis and Machine Intelligence, 8:679-714, 1986
Originally part of CellProfiler, code licensed under both GPL and BSD licenses.
Website: http://www.cellprofiler.org
Copyright (c) 2003-2009 Massachusetts Institute of Technology
Copyright (c) 2009-2011 Broad Institute
All rights reserved.
Original author: Lee Kamentsky
"""
import numpy as np
import scipy.ndimage as ndi
from ..util.dtype import dtype_limits
from .._shared.filters import gaussian
from .._shared.utils import check_nD
from ..filters._gaussian import gaussian
def _preprocess(image, mask, sigma, mode, cval):
"""Generate a smoothed image and an eroded mask.
The image is smoothed using a gaussian filter ignoring masked
pixels and the mask is eroded.
Parameters
----------
image : array
Image to be smoothed.
mask : array
Mask with 1's for significant pixels, 0's for masked pixels.
sigma : scalar or sequence of scalars
Standard deviation for Gaussian kernel. The standard
deviations of the Gaussian filter are given for each axis as a
sequence, or as a single number, in which case it is equal for
all axes.
mode : str, {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}
The ``mode`` parameter determines how the array borders are
handled, where ``cval`` is the value when mode is equal to
'constant'.
cval : float, optional
Value to fill past edges of input if `mode` is 'constant'.
Returns
-------
smoothed_image : ndarray
The smoothed array
eroded_mask : ndarray
The eroded mask.
Notes
-----
This function calculates the fractional contribution of masked pixels
by applying the function to the mask (which gets you the fraction of
the pixel data that's due to significant points). We then mask the image
and apply the function. The resulting values will be lower by the
bleed-over fraction, so you can recalibrate by dividing by the function
on the mask to recover the effect of smoothing from just the significant
pixels.
"""
gaussian_kwargs = dict(sigma=sigma, mode=mode, cval=cval,
preserve_range=False)
if mask is None:
# Smooth the masked image
smoothed_image = gaussian(image, **gaussian_kwargs)
eroded_mask = np.ones(image.shape, dtype=bool)
eroded_mask[:1, :] = 0
eroded_mask[-1:, :] = 0
eroded_mask[:, :1] = 0
eroded_mask[:, -1:] = 0
return smoothed_image, eroded_mask
masked_image = np.zeros_like(image)
masked_image[mask] = image[mask]
# Compute the fractional contribution of masked pixels by applying
# the function to the mask (which gets you the fraction of the
# pixel data that's due to significant points)
bleed_over = (
gaussian(mask.astype(float), **gaussian_kwargs) + np.finfo(float).eps
)
# Smooth the masked image
smoothed_image = gaussian(masked_image, **gaussian_kwargs)
# Lower the result by the bleed-over fraction, so you can
# recalibrate by dividing by the function on the mask to recover
# the effect of smoothing from just the significant pixels.
smoothed_image /= bleed_over
# Make the eroded mask. Setting the border value to zero will wipe
# out the image edges for us.
s = ndi.generate_binary_structure(2, 2)
eroded_mask = ndi.binary_erosion(mask, s, border_value=0)
return smoothed_image, eroded_mask
def _set_local_maxima(magnitude, pts, w_num, w_denum, row_slices,
col_slices, out):
"""Get the magnitudes shifted left to make a matrix of the points to
the right of pts. Similarly, shift left and down to get the points
to the top right of pts.
"""
r_0, r_1, r_2, r_3 = row_slices
c_0, c_1, c_2, c_3 = col_slices
c1 = magnitude[r_0, c_0][pts[r_1, c_1]]
c2 = magnitude[r_2, c_2][pts[r_3, c_3]]
m = magnitude[pts]
w = w_num[pts] / w_denum[pts]
c_plus = c2 * w + c1 * (1 - w) <= m
c1 = magnitude[r_1, c_1][pts[r_0, c_0]]
c2 = magnitude[r_3, c_3][pts[r_2, c_2]]
c_minus = c2 * w + c1 * (1 - w) <= m
out[pts] = c_plus & c_minus
return out
def _get_local_maxima(isobel, jsobel, magnitude, eroded_mask):
"""Edge thinning by non-maximum suppression.
Finds the normal to the edge at each point using the arctangent of the
ratio of the Y sobel over the X sobel - pragmatically, we can
look at the signs of X and Y and the relative magnitude of X vs Y
to sort the points into 4 categories: horizontal, vertical,
diagonal and antidiagonal.
Look in the normal and reverse directions to see if the values
in either of those directions are greater than the point in question.
Use interpolation (via _set_local_maxima) to get a mix of points
instead of picking the one that's the closest to the normal.
"""
abs_isobel = np.abs(isobel)
abs_jsobel = np.abs(jsobel)
eroded_mask = eroded_mask & (magnitude > 0)
# Normals' orientations
is_horizontal = eroded_mask & (abs_isobel >= abs_jsobel)
is_vertical = eroded_mask & (abs_isobel <= abs_jsobel)
is_up = (isobel >= 0)
is_down = (isobel <= 0)
is_right = (jsobel >= 0)
is_left = (jsobel <= 0)
#
# --------- Find local maxima --------------
#
# Assign each point to have a normal of 0-45 degrees, 45-90 degrees,
# 90-135 degrees and 135-180 degrees.
#
local_maxima = np.zeros(magnitude.shape, bool)
# ----- 0 to 45 degrees ------
# Mix diagonal and horizontal
pts_plus = is_up & is_right
pts_minus = is_down & is_left
pts = ((pts_plus | pts_minus) & is_horizontal)
# Get the magnitudes shifted left to make a matrix of the points to the
# right of pts. Similarly, shift left and down to get the points to the
# top right of pts.
local_maxima = _set_local_maxima(
magnitude, pts, abs_jsobel, abs_isobel,
[slice(1, None), slice(-1), slice(1, None), slice(-1)],
[slice(None), slice(None), slice(1, None), slice(-1)],
local_maxima)
# ----- 45 to 90 degrees ------
# Mix diagonal and vertical
#
pts = ((pts_plus | pts_minus) & is_vertical)
local_maxima = _set_local_maxima(
magnitude, pts, abs_isobel, abs_jsobel,
[slice(None), slice(None), slice(1, None), slice(-1)],
[slice(1, None), slice(-1), slice(1, None), slice(-1)],
local_maxima)
# ----- 90 to 135 degrees ------
# Mix anti-diagonal and vertical
#
pts_plus = is_down & is_right
pts_minus = is_up & is_left
pts = ((pts_plus | pts_minus) & is_vertical)
local_maxima = _set_local_maxima(
magnitude, pts, abs_isobel, abs_jsobel,
[slice(None), slice(None), slice(-1), slice(1, None)],
[slice(1, None), slice(-1), slice(1, None), slice(-1)],
local_maxima)
# ----- 135 to 180 degrees ------
# Mix anti-diagonal and anti-horizontal
#
pts = ((pts_plus | pts_minus) & is_horizontal)
local_maxima = _set_local_maxima(
magnitude, pts, abs_jsobel, abs_isobel,
[slice(-1), slice(1, None), slice(-1), slice(1, None)],
[slice(None), slice(None), slice(1, None), slice(-1)],
local_maxima)
return local_maxima
def canny(image, sigma=1., low_threshold=None, high_threshold=None,
mask=None, use_quantiles=False, *, mode='constant', cval=0.0):
"""Edge filter an image using the Canny algorithm.
Parameters
----------
image : 2D array
Grayscale input image to detect edges on; can be of any dtype.
sigma : float, optional
Standard deviation of the Gaussian filter.
low_threshold : float, optional
Lower bound for hysteresis thresholding (linking edges).
If None, low_threshold is set to 10% of dtype's max.
high_threshold : float, optional
Upper bound for hysteresis thresholding (linking edges).
If None, high_threshold is set to 20% of dtype's max.
mask : array, dtype=bool, optional
Mask to limit the application of Canny to a certain area.
use_quantiles : bool, optional
If ``True`` then treat low_threshold and high_threshold as
quantiles of the edge magnitude image, rather than absolute
edge magnitude values. If ``True`` then the thresholds must be
in the range [0, 1].
mode : str, {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}
The ``mode`` parameter determines how the array borders are
handled during Gaussian filtering, where ``cval`` is the value when
mode is equal to 'constant'.
cval : float, optional
Value to fill past edges of input if `mode` is 'constant'.
Returns
-------
output : 2D array (image)
The binary edge map.
See also
--------
skimage.sobel
Notes
-----
The steps of the algorithm are as follows:
* Smooth the image using a Gaussian with ``sigma`` width.
* Apply the horizontal and vertical Sobel operators to get the gradients
within the image. The edge strength is the norm of the gradient.
* Thin potential edges to 1-pixel wide curves. First, find the normal
to the edge at each point. This is done by looking at the
signs and the relative magnitude of the X-Sobel and Y-Sobel
to sort the points into 4 categories: horizontal, vertical,
diagonal and antidiagonal. Then look in the normal and reverse
directions to see if the values in either of those directions are
greater than the point in question. Use interpolation to get a mix of
points instead of picking the one that's the closest to the normal.
* Perform a hysteresis thresholding: first label all points above the
high threshold as edges. Then recursively label any point above the
low threshold that is 8-connected to a labeled point as an edge.
References
----------
.. [1] Canny, J., A Computational Approach To Edge Detection, IEEE Trans.
Pattern Analysis and Machine Intelligence, 8:679-714, 1986
:DOI:`10.1109/TPAMI.1986.4767851`
.. [2] William Green's Canny tutorial
https://en.wikipedia.org/wiki/Canny_edge_detector
Examples
--------
>>> from skimage import feature
>>> rng = np.random.default_rng()
>>> # Generate noisy image of a square
>>> im = np.zeros((256, 256))
>>> im[64:-64, 64:-64] = 1
>>> im += 0.2 * rng.random(im.shape)
>>> # First trial with the Canny filter, with the default smoothing
>>> edges1 = feature.canny(im)
>>> # Increase the smoothing for better results
>>> edges2 = feature.canny(im, sigma=3)
"""
# Regarding masks, any point touching a masked point will have a gradient
# that is "infected" by the masked point, so it's enough to erode the
# mask by one and then mask the output. We also mask out the border points
# because who knows what lies beyond the edge of the image?
check_nD(image, 2)
dtype_max = dtype_limits(image, clip_negative=False)[1]
if low_threshold is None:
low_threshold = 0.1
elif use_quantiles:
if not(0.0 <= low_threshold <= 1.0):
raise ValueError("Quantile thresholds must be between 0 and 1.")
else:
low_threshold /= dtype_max
if high_threshold is None:
high_threshold = 0.2
elif use_quantiles:
if not(0.0 <= high_threshold <= 1.0):
raise ValueError("Quantile thresholds must be between 0 and 1.")
else:
high_threshold /= dtype_max
if high_threshold < low_threshold:
raise ValueError("low_threshold should be lower then high_threshold")
# Image filtering
smoothed, eroded_mask = _preprocess(image, mask, sigma, mode, cval)
# Gradient magnitude estimation
jsobel = ndi.sobel(smoothed, axis=1)
isobel = ndi.sobel(smoothed, axis=0)
magnitude = np.hypot(isobel, jsobel)
if use_quantiles:
low_threshold, high_threshold = np.percentile(magnitude,
[100.0 * low_threshold,
100.0 * high_threshold])
# Non-maximum suppression
local_maxima = _get_local_maxima(isobel, jsobel, magnitude, eroded_mask)
# Double thresholding and edge traking
low_mask = local_maxima & (magnitude >= low_threshold)
#
# Segment the low-mask, then only keep low-segments that have
# some high_mask component in them
#
strel = np.ones((3, 3), bool)
labels, count = ndi.label(low_mask, strel)
if count == 0:
return low_mask
high_mask = local_maxima & (magnitude >= high_threshold)
nonzero_sums = np.unique(labels[high_mask])
good_label = np.zeros((count + 1,), bool)
good_label[nonzero_sums] = True
output_mask = good_label[labels]
return output_mask