1"""extrema.py - local minima and maxima 2 3This module provides functions to find local maxima and minima of an image. 4Here, local maxima (minima) are defined as connected sets of pixels with equal 5gray level which is strictly greater (smaller) than the gray level of all 6pixels in direct neighborhood of the connected set. In addition, the module 7provides the related functions h-maxima and h-minima. 8 9Soille, P. (2003). Morphological Image Analysis: Principles and Applications 10(2nd ed.), Chapter 6. Springer-Verlag New York, Inc. 11""" 12import numpy as np 13 14from .._shared.utils import deprecate_kwarg, warn 15from ..util import dtype_limits, invert, crop 16from . import grayreconstruct, _util 17from ._extrema_cy import _local_maxima 18 19 20def _add_constant_clip(image, const_value): 21 """Add constant to the image while handling overflow issues gracefully. 22 """ 23 min_dtype, max_dtype = dtype_limits(image, clip_negative=False) 24 25 if const_value > (max_dtype - min_dtype): 26 raise ValueError("The added constant is not compatible" 27 "with the image data type.") 28 29 result = image + const_value 30 result[image > max_dtype-const_value] = max_dtype 31 return(result) 32 33 34def _subtract_constant_clip(image, const_value): 35 """Subtract constant from image while handling underflow issues. 36 """ 37 min_dtype, max_dtype = dtype_limits(image, clip_negative=False) 38 39 if const_value > (max_dtype-min_dtype): 40 raise ValueError("The subtracted constant is not compatible" 41 "with the image data type.") 42 43 result = image - const_value 44 result[image < (const_value + min_dtype)] = min_dtype 45 return(result) 46 47 48@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'}, removed_version="1.0", 49 deprecated_version="0.19") 50def h_maxima(image, h, footprint=None): 51 """Determine all maxima of the image with height >= h. 52 53 The local maxima are defined as connected sets of pixels with equal 54 gray level strictly greater than the gray level of all pixels in direct 55 neighborhood of the set. 56 57 A local maximum M of height h is a local maximum for which 58 there is at least one path joining M with an equal or higher local maximum 59 on which the minimal value is f(M) - h (i.e. the values along the path 60 are not decreasing by more than h with respect to the maximum's value) 61 and no path to an equal or higher local maximum for which the minimal 62 value is greater. 63 64 The global maxima of the image are also found by this function. 65 66 Parameters 67 ---------- 68 image : ndarray 69 The input image for which the maxima are to be calculated. 70 h : unsigned integer 71 The minimal height of all extracted maxima. 72 footprint : ndarray, optional 73 The neighborhood expressed as an n-D array of 1's and 0's. 74 Default is the ball of radius 1 according to the maximum norm 75 (i.e. a 3x3 square for 2D images, a 3x3x3 cube for 3D images, etc.) 76 77 Returns 78 ------- 79 h_max : ndarray 80 The local maxima of height >= h and the global maxima. 81 The resulting image is a binary image, where pixels belonging to 82 the determined maxima take value 1, the others take value 0. 83 84 See Also 85 -------- 86 skimage.morphology.extrema.h_minima 87 skimage.morphology.extrema.local_maxima 88 skimage.morphology.extrema.local_minima 89 90 References 91 ---------- 92 .. [1] Soille, P., "Morphological Image Analysis: Principles and 93 Applications" (Chapter 6), 2nd edition (2003), ISBN 3540429883. 94 95 Examples 96 -------- 97 >>> import numpy as np 98 >>> from skimage.morphology import extrema 99 100 We create an image (quadratic function with a maximum in the center and 101 4 additional constant maxima. 102 The heights of the maxima are: 1, 21, 41, 61, 81 103 104 >>> w = 10 105 >>> x, y = np.mgrid[0:w,0:w] 106 >>> f = 20 - 0.2*((x - w/2)**2 + (y-w/2)**2) 107 >>> f[2:4,2:4] = 40; f[2:4,7:9] = 60; f[7:9,2:4] = 80; f[7:9,7:9] = 100 108 >>> f = f.astype(int) 109 110 We can calculate all maxima with a height of at least 40: 111 112 >>> maxima = extrema.h_maxima(f, 40) 113 114 The resulting image will contain 3 local maxima. 115 """ 116 117 # Check for h value that is larger then range of the image. If this 118 # is True then there are no h-maxima in the image. 119 if h > np.ptp(image): 120 return np.zeros(image.shape, dtype=np.uint8) 121 122 # Check for floating point h value. For this to work properly 123 # we need to explicitly convert image to float64. 124 # 125 # FIXME: This could give incorrect results if image is int64 and 126 # has a very high dynamic range. The dtype of image is 127 # changed to float64, and different integer values could 128 # become the same float due to rounding. 129 # 130 # >>> ii64 = np.iinfo(np.int64) 131 # >>> a = np.array([ii64.max, ii64.max - 2]) 132 # >>> a[0] == a[1] 133 # False 134 # >>> b = a.astype(np.float64) 135 # >>> b[0] == b[1] 136 # True 137 # 138 if np.issubdtype(type(h), np.floating) and \ 139 np.issubdtype(image.dtype, np.integer): 140 if ((h % 1) != 0): 141 warn('possible precision loss converting image to ' 142 'floating point. To silence this warning, ' 143 'ensure image and h have same data type.', 144 stacklevel=2) 145 image = image.astype(float) 146 else: 147 h = image.dtype.type(h) 148 149 if (h == 0): 150 raise ValueError("h = 0 is ambiguous, use local_maxima() " 151 "instead?") 152 153 if np.issubdtype(image.dtype, np.floating): 154 # The purpose of the resolution variable is to allow for the 155 # small rounding errors that inevitably occur when doing 156 # floating point arithmetic. We want shifted_img to be 157 # guaranteed to be h less than image. If we only subtract h 158 # there may be pixels were shifted_img ends up being 159 # slightly greater than image - h. 160 # 161 # The resolution is scaled based on the pixel values in the 162 # image because floating point precision is relative. A 163 # very large value of 1.0e10 will have a large precision, 164 # say +-1.0e4, and a very small value of 1.0e-10 will have 165 # a very small precision, say +-1.0e-16. 166 # 167 resolution = 2 * np.finfo(image.dtype).resolution * np.abs(image) 168 shifted_img = image - h - resolution 169 else: 170 shifted_img = _subtract_constant_clip(image, h) 171 172 rec_img = grayreconstruct.reconstruction(shifted_img, image, 173 method='dilation', 174 footprint=footprint) 175 residue_img = image - rec_img 176 return (residue_img >= h).astype(np.uint8) 177 178 179@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'}, removed_version="1.0", 180 deprecated_version="0.19") 181def h_minima(image, h, footprint=None): 182 """Determine all minima of the image with depth >= h. 183 184 The local minima are defined as connected sets of pixels with equal 185 gray level strictly smaller than the gray levels of all pixels in direct 186 neighborhood of the set. 187 188 A local minimum M of depth h is a local minimum for which 189 there is at least one path joining M with an equal or lower local minimum 190 on which the maximal value is f(M) + h (i.e. the values along the path 191 are not increasing by more than h with respect to the minimum's value) 192 and no path to an equal or lower local minimum for which the maximal 193 value is smaller. 194 195 The global minima of the image are also found by this function. 196 197 Parameters 198 ---------- 199 image : ndarray 200 The input image for which the minima are to be calculated. 201 h : unsigned integer 202 The minimal depth of all extracted minima. 203 footprint : ndarray, optional 204 The neighborhood expressed as an n-D array of 1's and 0's. 205 Default is the ball of radius 1 according to the maximum norm 206 (i.e. a 3x3 square for 2D images, a 3x3x3 cube for 3D images, etc.) 207 208 Returns 209 ------- 210 h_min : ndarray 211 The local minima of depth >= h and the global minima. 212 The resulting image is a binary image, where pixels belonging to 213 the determined minima take value 1, the others take value 0. 214 215 See Also 216 -------- 217 skimage.morphology.extrema.h_maxima 218 skimage.morphology.extrema.local_maxima 219 skimage.morphology.extrema.local_minima 220 221 References 222 ---------- 223 .. [1] Soille, P., "Morphological Image Analysis: Principles and 224 Applications" (Chapter 6), 2nd edition (2003), ISBN 3540429883. 225 226 Examples 227 -------- 228 >>> import numpy as np 229 >>> from skimage.morphology import extrema 230 231 We create an image (quadratic function with a minimum in the center and 232 4 additional constant maxima. 233 The depth of the minima are: 1, 21, 41, 61, 81 234 235 >>> w = 10 236 >>> x, y = np.mgrid[0:w,0:w] 237 >>> f = 180 + 0.2*((x - w/2)**2 + (y-w/2)**2) 238 >>> f[2:4,2:4] = 160; f[2:4,7:9] = 140; f[7:9,2:4] = 120; f[7:9,7:9] = 100 239 >>> f = f.astype(int) 240 241 We can calculate all minima with a depth of at least 40: 242 243 >>> minima = extrema.h_minima(f, 40) 244 245 The resulting image will contain 3 local minima. 246 """ 247 if h > np.ptp(image): 248 return np.zeros(image.shape, dtype=np.uint8) 249 250 if np.issubdtype(type(h), np.floating) and \ 251 np.issubdtype(image.dtype, np.integer): 252 if ((h % 1) != 0): 253 warn('possible precision loss converting image to ' 254 'floating point. To silence this warning, ' 255 'ensure image and h have same data type.', 256 stacklevel=2) 257 image = image.astype(float) 258 else: 259 h = image.dtype.type(h) 260 261 if (h == 0): 262 raise ValueError("h = 0 is ambiguous, use local_minima() " 263 "instead?") 264 265 if np.issubdtype(image.dtype, np.floating): 266 resolution = 2 * np.finfo(image.dtype).resolution * np.abs(image) 267 shifted_img = image + h + resolution 268 else: 269 shifted_img = _add_constant_clip(image, h) 270 271 rec_img = grayreconstruct.reconstruction(shifted_img, image, 272 method='erosion', 273 footprint=footprint) 274 residue_img = rec_img - image 275 return (residue_img >= h).astype(np.uint8) 276 277 278@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'}, removed_version="1.0", 279 deprecated_version="0.19") 280def local_maxima(image, footprint=None, connectivity=None, indices=False, 281 allow_borders=True): 282 """Find local maxima of n-dimensional array. 283 284 The local maxima are defined as connected sets of pixels with equal gray 285 level (plateaus) strictly greater than the gray levels of all pixels in the 286 neighborhood. 287 288 Parameters 289 ---------- 290 image : ndarray 291 An n-dimensional array. 292 footprint : ndarray, optional 293 The footprint (structuring element) used to determine the neighborhood 294 of each evaluated pixel (``True`` denotes a connected pixel). It must 295 be a boolean array and have the same number of dimensions as `image`. 296 If neither `footprint` nor `connectivity` are given, all adjacent 297 pixels are considered as part of the neighborhood. 298 connectivity : int, optional 299 A number used to determine the neighborhood of each evaluated pixel. 300 Adjacent pixels whose squared distance from the center is less than or 301 equal to `connectivity` are considered neighbors. Ignored if 302 `footprint` is not None. 303 indices : bool, optional 304 If True, the output will be a tuple of one-dimensional arrays 305 representing the indices of local maxima in each dimension. If False, 306 the output will be a boolean array with the same shape as `image`. 307 allow_borders : bool, optional 308 If true, plateaus that touch the image border are valid maxima. 309 310 Returns 311 ------- 312 maxima : ndarray or tuple[ndarray] 313 If `indices` is false, a boolean array with the same shape as `image` 314 is returned with ``True`` indicating the position of local maxima 315 (``False`` otherwise). If `indices` is true, a tuple of one-dimensional 316 arrays containing the coordinates (indices) of all found maxima. 317 318 Warns 319 ----- 320 UserWarning 321 If `allow_borders` is false and any dimension of the given `image` is 322 shorter than 3 samples, maxima can't exist and a warning is shown. 323 324 See Also 325 -------- 326 skimage.morphology.local_minima 327 skimage.morphology.h_maxima 328 skimage.morphology.h_minima 329 330 Notes 331 ----- 332 This function operates on the following ideas: 333 334 1. Make a first pass over the image's last dimension and flag candidates 335 for local maxima by comparing pixels in only one direction. 336 If the pixels aren't connected in the last dimension all pixels are 337 flagged as candidates instead. 338 339 For each candidate: 340 341 2. Perform a flood-fill to find all connected pixels that have the same 342 gray value and are part of the plateau. 343 3. Consider the connected neighborhood of a plateau: if no bordering sample 344 has a higher gray level, mark the plateau as a definite local maximum. 345 346 Examples 347 -------- 348 >>> from skimage.morphology import local_maxima 349 >>> image = np.zeros((4, 7), dtype=int) 350 >>> image[1:3, 1:3] = 1 351 >>> image[3, 0] = 1 352 >>> image[1:3, 4:6] = 2 353 >>> image[3, 6] = 3 354 >>> image 355 array([[0, 0, 0, 0, 0, 0, 0], 356 [0, 1, 1, 0, 2, 2, 0], 357 [0, 1, 1, 0, 2, 2, 0], 358 [1, 0, 0, 0, 0, 0, 3]]) 359 360 Find local maxima by comparing to all neighboring pixels (maximal 361 connectivity): 362 363 >>> local_maxima(image) 364 array([[False, False, False, False, False, False, False], 365 [False, True, True, False, False, False, False], 366 [False, True, True, False, False, False, False], 367 [ True, False, False, False, False, False, True]]) 368 >>> local_maxima(image, indices=True) 369 (array([1, 1, 2, 2, 3, 3]), array([1, 2, 1, 2, 0, 6])) 370 371 Find local maxima without comparing to diagonal pixels (connectivity 1): 372 373 >>> local_maxima(image, connectivity=1) 374 array([[False, False, False, False, False, False, False], 375 [False, True, True, False, True, True, False], 376 [False, True, True, False, True, True, False], 377 [ True, False, False, False, False, False, True]]) 378 379 and exclude maxima that border the image edge: 380 381 >>> local_maxima(image, connectivity=1, allow_borders=False) 382 array([[False, False, False, False, False, False, False], 383 [False, True, True, False, True, True, False], 384 [False, True, True, False, True, True, False], 385 [False, False, False, False, False, False, False]]) 386 """ 387 image = np.asarray(image, order="C") 388 if image.size == 0: 389 # Return early for empty input 390 if indices: 391 # Make sure that output is a tuple of 1 empty array per dimension 392 return np.nonzero(image) 393 else: 394 return np.zeros(image.shape, dtype=bool) 395 396 if allow_borders: 397 # Ensure that local maxima are always at least one smaller sample away 398 # from the image border 399 image = np.pad(image, 1, mode='constant', constant_values=image.min()) 400 401 # Array of flags used to store the state of each pixel during evaluation. 402 # See _extrema_cy.pyx for their meaning 403 flags = np.zeros(image.shape, dtype=np.uint8) 404 _util._set_border_values(flags, value=3) 405 406 if any(s < 3 for s in image.shape): 407 # Warn and skip if any dimension is smaller than 3 408 # -> no maxima can exist & footprint can't be applied 409 warn( 410 "maxima can't exist for an image with any dimension smaller 3 " 411 "if borders aren't allowed", 412 stacklevel=3 413 ) 414 else: 415 footprint = _util._resolve_neighborhood(footprint, connectivity, 416 image.ndim) 417 neighbor_offsets = _util._offsets_to_raveled_neighbors( 418 image.shape, footprint, center=((1,) * image.ndim) 419 ) 420 421 try: 422 _local_maxima(image.ravel(), flags.ravel(), neighbor_offsets) 423 except TypeError: 424 if image.dtype == np.float16: 425 # Provide the user with clearer error message 426 raise TypeError("dtype of `image` is float16 which is not " 427 "supported, try upcasting to float32") 428 else: 429 raise # Otherwise raise original message 430 431 if allow_borders: 432 # Revert padding performed at the beginning of the function 433 flags = crop(flags, 1) 434 else: 435 # No padding was performed but set edge values back to 0 436 _util._set_border_values(flags, value=0) 437 438 if indices: 439 return np.nonzero(flags) 440 else: 441 return flags.view(bool) 442 443 444@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'}, removed_version="1.0", 445 deprecated_version="0.19") 446def local_minima(image, footprint=None, connectivity=None, indices=False, 447 allow_borders=True): 448 """Find local minima of n-dimensional array. 449 450 The local minima are defined as connected sets of pixels with equal gray 451 level (plateaus) strictly smaller than the gray levels of all pixels in the 452 neighborhood. 453 454 Parameters 455 ---------- 456 image : ndarray 457 An n-dimensional array. 458 footprint : ndarray, optional 459 The footprint (structuring element) used to determine the neighborhood 460 of each evaluated pixel (``True`` denotes a connected pixel). It must 461 be a boolean array and have the same number of dimensions as `image`. 462 If neither `footprint` nor `connectivity` are given, all adjacent 463 pixels are considered as part of the neighborhood. 464 connectivity : int, optional 465 A number used to determine the neighborhood of each evaluated pixel. 466 Adjacent pixels whose squared distance from the center is less than or 467 equal to `connectivity` are considered neighbors. Ignored if 468 `footprint` is not None. 469 indices : bool, optional 470 If True, the output will be a tuple of one-dimensional arrays 471 representing the indices of local minima in each dimension. If False, 472 the output will be a boolean array with the same shape as `image`. 473 allow_borders : bool, optional 474 If true, plateaus that touch the image border are valid minima. 475 476 Returns 477 ------- 478 minima : ndarray or tuple[ndarray] 479 If `indices` is false, a boolean array with the same shape as `image` 480 is returned with ``True`` indicating the position of local minima 481 (``False`` otherwise). If `indices` is true, a tuple of one-dimensional 482 arrays containing the coordinates (indices) of all found minima. 483 484 See Also 485 -------- 486 skimage.morphology.local_maxima 487 skimage.morphology.h_maxima 488 skimage.morphology.h_minima 489 490 Notes 491 ----- 492 This function operates on the following ideas: 493 494 1. Make a first pass over the image's last dimension and flag candidates 495 for local minima by comparing pixels in only one direction. 496 If the pixels aren't connected in the last dimension all pixels are 497 flagged as candidates instead. 498 499 For each candidate: 500 501 2. Perform a flood-fill to find all connected pixels that have the same 502 gray value and are part of the plateau. 503 3. Consider the connected neighborhood of a plateau: if no bordering sample 504 has a smaller gray level, mark the plateau as a definite local minimum. 505 506 Examples 507 -------- 508 >>> from skimage.morphology import local_minima 509 >>> image = np.zeros((4, 7), dtype=int) 510 >>> image[1:3, 1:3] = -1 511 >>> image[3, 0] = -1 512 >>> image[1:3, 4:6] = -2 513 >>> image[3, 6] = -3 514 >>> image 515 array([[ 0, 0, 0, 0, 0, 0, 0], 516 [ 0, -1, -1, 0, -2, -2, 0], 517 [ 0, -1, -1, 0, -2, -2, 0], 518 [-1, 0, 0, 0, 0, 0, -3]]) 519 520 Find local minima by comparing to all neighboring pixels (maximal 521 connectivity): 522 523 >>> local_minima(image) 524 array([[False, False, False, False, False, False, False], 525 [False, True, True, False, False, False, False], 526 [False, True, True, False, False, False, False], 527 [ True, False, False, False, False, False, True]]) 528 >>> local_minima(image, indices=True) 529 (array([1, 1, 2, 2, 3, 3]), array([1, 2, 1, 2, 0, 6])) 530 531 Find local minima without comparing to diagonal pixels (connectivity 1): 532 533 >>> local_minima(image, connectivity=1) 534 array([[False, False, False, False, False, False, False], 535 [False, True, True, False, True, True, False], 536 [False, True, True, False, True, True, False], 537 [ True, False, False, False, False, False, True]]) 538 539 and exclude minima that border the image edge: 540 541 >>> local_minima(image, connectivity=1, allow_borders=False) 542 array([[False, False, False, False, False, False, False], 543 [False, True, True, False, True, True, False], 544 [False, True, True, False, True, True, False], 545 [False, False, False, False, False, False, False]]) 546 """ 547 return local_maxima( 548 image=invert(image), 549 footprint=footprint, 550 connectivity=connectivity, 551 indices=indices, 552 allow_borders=allow_borders 553 ) 554