1# Copyright (C) 2003-2005 Peter J. Verveer 2# 3# Redistribution and use in source and binary forms, with or without 4# modification, are permitted provided that the following conditions 5# are met: 6# 7# 1. Redistributions of source code must retain the above copyright 8# notice, this list of conditions and the following disclaimer. 9# 10# 2. Redistributions in binary form must reproduce the above 11# copyright notice, this list of conditions and the following 12# disclaimer in the documentation and/or other materials provided 13# with the distribution. 14# 15# 3. The name of the author may not be used to endorse or promote 16# products derived from this software without specific prior 17# written permission. 18# 19# THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS 20# OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 21# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 22# ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY 23# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 24# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE 25# GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 26# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 27# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 28# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 29# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 30 31import numpy 32import numpy as np 33from . import _ni_support 34from . import _ni_label 35from . import _nd_image 36from . import morphology 37 38__all__ = ['label', 'find_objects', 'labeled_comprehension', 'sum', 'mean', 39 'variance', 'standard_deviation', 'minimum', 'maximum', 'median', 40 'minimum_position', 'maximum_position', 'extrema', 'center_of_mass', 41 'histogram', 'watershed_ift', 'sum_labels'] 42 43 44def label(input, structure=None, output=None): 45 """ 46 Label features in an array. 47 48 Parameters 49 ---------- 50 input : array_like 51 An array-like object to be labeled. Any non-zero values in `input` are 52 counted as features and zero values are considered the background. 53 structure : array_like, optional 54 A structuring element that defines feature connections. 55 `structure` must be centrosymmetric 56 (see Notes). 57 If no structuring element is provided, 58 one is automatically generated with a squared connectivity equal to 59 one. That is, for a 2-D `input` array, the default structuring element 60 is:: 61 62 [[0,1,0], 63 [1,1,1], 64 [0,1,0]] 65 66 output : (None, data-type, array_like), optional 67 If `output` is a data type, it specifies the type of the resulting 68 labeled feature array. 69 If `output` is an array-like object, then `output` will be updated 70 with the labeled features from this function. This function can 71 operate in-place, by passing output=input. 72 Note that the output must be able to store the largest label, or this 73 function will raise an Exception. 74 75 Returns 76 ------- 77 label : ndarray or int 78 An integer ndarray where each unique feature in `input` has a unique 79 label in the returned array. 80 num_features : int 81 How many objects were found. 82 83 If `output` is None, this function returns a tuple of 84 (`labeled_array`, `num_features`). 85 86 If `output` is a ndarray, then it will be updated with values in 87 `labeled_array` and only `num_features` will be returned by this 88 function. 89 90 See Also 91 -------- 92 find_objects : generate a list of slices for the labeled features (or 93 objects); useful for finding features' position or 94 dimensions 95 96 Notes 97 ----- 98 A centrosymmetric matrix is a matrix that is symmetric about the center. 99 See [1]_ for more information. 100 101 The `structure` matrix must be centrosymmetric to ensure 102 two-way connections. 103 For instance, if the `structure` matrix is not centrosymmetric 104 and is defined as:: 105 106 [[0,1,0], 107 [1,1,0], 108 [0,0,0]] 109 110 and the `input` is:: 111 112 [[1,2], 113 [0,3]] 114 115 then the structure matrix would indicate the 116 entry 2 in the input is connected to 1, 117 but 1 is not connected to 2. 118 119 Examples 120 -------- 121 Create an image with some features, then label it using the default 122 (cross-shaped) structuring element: 123 124 >>> from scipy.ndimage import label, generate_binary_structure 125 >>> a = np.array([[0,0,1,1,0,0], 126 ... [0,0,0,1,0,0], 127 ... [1,1,0,0,1,0], 128 ... [0,0,0,1,0,0]]) 129 >>> labeled_array, num_features = label(a) 130 131 Each of the 4 features are labeled with a different integer: 132 133 >>> num_features 134 4 135 >>> labeled_array 136 array([[0, 0, 1, 1, 0, 0], 137 [0, 0, 0, 1, 0, 0], 138 [2, 2, 0, 0, 3, 0], 139 [0, 0, 0, 4, 0, 0]]) 140 141 Generate a structuring element that will consider features connected even 142 if they touch diagonally: 143 144 >>> s = generate_binary_structure(2,2) 145 146 or, 147 148 >>> s = [[1,1,1], 149 ... [1,1,1], 150 ... [1,1,1]] 151 152 Label the image using the new structuring element: 153 154 >>> labeled_array, num_features = label(a, structure=s) 155 156 Show the 2 labeled features (note that features 1, 3, and 4 from above are 157 now considered a single feature): 158 159 >>> num_features 160 2 161 >>> labeled_array 162 array([[0, 0, 1, 1, 0, 0], 163 [0, 0, 0, 1, 0, 0], 164 [2, 2, 0, 0, 1, 0], 165 [0, 0, 0, 1, 0, 0]]) 166 167 References 168 ---------- 169 170 .. [1] James R. Weaver, "Centrosymmetric (cross-symmetric) 171 matrices, their basic properties, eigenvalues, and 172 eigenvectors." The American Mathematical Monthly 92.10 173 (1985): 711-717. 174 175 """ 176 input = numpy.asarray(input) 177 if numpy.iscomplexobj(input): 178 raise TypeError('Complex type not supported') 179 if structure is None: 180 structure = morphology.generate_binary_structure(input.ndim, 1) 181 structure = numpy.asarray(structure, dtype=bool) 182 if structure.ndim != input.ndim: 183 raise RuntimeError('structure and input must have equal rank') 184 for ii in structure.shape: 185 if ii != 3: 186 raise ValueError('structure dimensions must be equal to 3') 187 188 # Use 32 bits if it's large enough for this image. 189 # _ni_label.label() needs two entries for background and 190 # foreground tracking 191 need_64bits = input.size >= (2**31 - 2) 192 193 if isinstance(output, numpy.ndarray): 194 if output.shape != input.shape: 195 raise ValueError("output shape not correct") 196 caller_provided_output = True 197 else: 198 caller_provided_output = False 199 if output is None: 200 output = np.empty(input.shape, np.intp if need_64bits else np.int32) 201 else: 202 output = np.empty(input.shape, output) 203 204 # handle scalars, 0-D arrays 205 if input.ndim == 0 or input.size == 0: 206 if input.ndim == 0: 207 # scalar 208 maxlabel = 1 if (input != 0) else 0 209 output[...] = maxlabel 210 else: 211 # 0-D 212 maxlabel = 0 213 if caller_provided_output: 214 return maxlabel 215 else: 216 return output, maxlabel 217 218 try: 219 max_label = _ni_label._label(input, structure, output) 220 except _ni_label.NeedMoreBits as e: 221 # Make another attempt with enough bits, then try to cast to the 222 # new type. 223 tmp_output = np.empty(input.shape, np.intp if need_64bits else np.int32) 224 max_label = _ni_label._label(input, structure, tmp_output) 225 output[...] = tmp_output[...] 226 if not np.all(output == tmp_output): 227 # refuse to return bad results 228 raise RuntimeError( 229 "insufficient bit-depth in requested output type" 230 ) from e 231 232 if caller_provided_output: 233 # result was written in-place 234 return max_label 235 else: 236 return output, max_label 237 238 239def find_objects(input, max_label=0): 240 """ 241 Find objects in a labeled array. 242 243 Parameters 244 ---------- 245 input : ndarray of ints 246 Array containing objects defined by different labels. Labels with 247 value 0 are ignored. 248 max_label : int, optional 249 Maximum label to be searched for in `input`. If max_label is not 250 given, the positions of all objects are returned. 251 252 Returns 253 ------- 254 object_slices : list of tuples 255 A list of tuples, with each tuple containing N slices (with N the 256 dimension of the input array). Slices correspond to the minimal 257 parallelepiped that contains the object. If a number is missing, 258 None is returned instead of a slice. 259 260 See Also 261 -------- 262 label, center_of_mass 263 264 Notes 265 ----- 266 This function is very useful for isolating a volume of interest inside 267 a 3-D array, that cannot be "seen through". 268 269 Examples 270 -------- 271 >>> from scipy import ndimage 272 >>> a = np.zeros((6,6), dtype=int) 273 >>> a[2:4, 2:4] = 1 274 >>> a[4, 4] = 1 275 >>> a[:2, :3] = 2 276 >>> a[0, 5] = 3 277 >>> a 278 array([[2, 2, 2, 0, 0, 3], 279 [2, 2, 2, 0, 0, 0], 280 [0, 0, 1, 1, 0, 0], 281 [0, 0, 1, 1, 0, 0], 282 [0, 0, 0, 0, 1, 0], 283 [0, 0, 0, 0, 0, 0]]) 284 >>> ndimage.find_objects(a) 285 [(slice(2, 5, None), slice(2, 5, None)), (slice(0, 2, None), slice(0, 3, None)), (slice(0, 1, None), slice(5, 6, None))] 286 >>> ndimage.find_objects(a, max_label=2) 287 [(slice(2, 5, None), slice(2, 5, None)), (slice(0, 2, None), slice(0, 3, None))] 288 >>> ndimage.find_objects(a == 1, max_label=2) 289 [(slice(2, 5, None), slice(2, 5, None)), None] 290 291 >>> loc = ndimage.find_objects(a)[0] 292 >>> a[loc] 293 array([[1, 1, 0], 294 [1, 1, 0], 295 [0, 0, 1]]) 296 297 """ 298 input = numpy.asarray(input) 299 if numpy.iscomplexobj(input): 300 raise TypeError('Complex type not supported') 301 302 if max_label < 1: 303 max_label = input.max() 304 305 return _nd_image.find_objects(input, max_label) 306 307 308def labeled_comprehension(input, labels, index, func, out_dtype, default, pass_positions=False): 309 """ 310 Roughly equivalent to [func(input[labels == i]) for i in index]. 311 312 Sequentially applies an arbitrary function (that works on array_like input) 313 to subsets of an N-D image array specified by `labels` and `index`. 314 The option exists to provide the function with positional parameters as the 315 second argument. 316 317 Parameters 318 ---------- 319 input : array_like 320 Data from which to select `labels` to process. 321 labels : array_like or None 322 Labels to objects in `input`. 323 If not None, array must be same shape as `input`. 324 If None, `func` is applied to raveled `input`. 325 index : int, sequence of ints or None 326 Subset of `labels` to which to apply `func`. 327 If a scalar, a single value is returned. 328 If None, `func` is applied to all non-zero values of `labels`. 329 func : callable 330 Python function to apply to `labels` from `input`. 331 out_dtype : dtype 332 Dtype to use for `result`. 333 default : int, float or None 334 Default return value when a element of `index` does not exist 335 in `labels`. 336 pass_positions : bool, optional 337 If True, pass linear indices to `func` as a second argument. 338 Default is False. 339 340 Returns 341 ------- 342 result : ndarray 343 Result of applying `func` to each of `labels` to `input` in `index`. 344 345 Examples 346 -------- 347 >>> a = np.array([[1, 2, 0, 0], 348 ... [5, 3, 0, 4], 349 ... [0, 0, 0, 7], 350 ... [9, 3, 0, 0]]) 351 >>> from scipy import ndimage 352 >>> lbl, nlbl = ndimage.label(a) 353 >>> lbls = np.arange(1, nlbl+1) 354 >>> ndimage.labeled_comprehension(a, lbl, lbls, np.mean, float, 0) 355 array([ 2.75, 5.5 , 6. ]) 356 357 Falling back to `default`: 358 359 >>> lbls = np.arange(1, nlbl+2) 360 >>> ndimage.labeled_comprehension(a, lbl, lbls, np.mean, float, -1) 361 array([ 2.75, 5.5 , 6. , -1. ]) 362 363 Passing positions: 364 365 >>> def fn(val, pos): 366 ... print("fn says: %s : %s" % (val, pos)) 367 ... return (val.sum()) if (pos.sum() % 2 == 0) else (-val.sum()) 368 ... 369 >>> ndimage.labeled_comprehension(a, lbl, lbls, fn, float, 0, True) 370 fn says: [1 2 5 3] : [0 1 4 5] 371 fn says: [4 7] : [ 7 11] 372 fn says: [9 3] : [12 13] 373 array([ 11., 11., -12., 0.]) 374 375 """ 376 377 as_scalar = numpy.isscalar(index) 378 input = numpy.asarray(input) 379 380 if pass_positions: 381 positions = numpy.arange(input.size).reshape(input.shape) 382 383 if labels is None: 384 if index is not None: 385 raise ValueError("index without defined labels") 386 if not pass_positions: 387 return func(input.ravel()) 388 else: 389 return func(input.ravel(), positions.ravel()) 390 391 try: 392 input, labels = numpy.broadcast_arrays(input, labels) 393 except ValueError as e: 394 raise ValueError("input and labels must have the same shape " 395 "(excepting dimensions with width 1)") from e 396 397 if index is None: 398 if not pass_positions: 399 return func(input[labels > 0]) 400 else: 401 return func(input[labels > 0], positions[labels > 0]) 402 403 index = numpy.atleast_1d(index) 404 if np.any(index.astype(labels.dtype).astype(index.dtype) != index): 405 raise ValueError("Cannot convert index values from <%s> to <%s> " 406 "(labels' type) without loss of precision" % 407 (index.dtype, labels.dtype)) 408 409 index = index.astype(labels.dtype) 410 411 # optimization: find min/max in index, and select those parts of labels, input, and positions 412 lo = index.min() 413 hi = index.max() 414 mask = (labels >= lo) & (labels <= hi) 415 416 # this also ravels the arrays 417 labels = labels[mask] 418 input = input[mask] 419 if pass_positions: 420 positions = positions[mask] 421 422 # sort everything by labels 423 label_order = labels.argsort() 424 labels = labels[label_order] 425 input = input[label_order] 426 if pass_positions: 427 positions = positions[label_order] 428 429 index_order = index.argsort() 430 sorted_index = index[index_order] 431 432 def do_map(inputs, output): 433 """labels must be sorted""" 434 nidx = sorted_index.size 435 436 # Find boundaries for each stretch of constant labels 437 # This could be faster, but we already paid N log N to sort labels. 438 lo = numpy.searchsorted(labels, sorted_index, side='left') 439 hi = numpy.searchsorted(labels, sorted_index, side='right') 440 441 for i, l, h in zip(range(nidx), lo, hi): 442 if l == h: 443 continue 444 output[i] = func(*[inp[l:h] for inp in inputs]) 445 446 temp = numpy.empty(index.shape, out_dtype) 447 temp[:] = default 448 if not pass_positions: 449 do_map([input], temp) 450 else: 451 do_map([input, positions], temp) 452 453 output = numpy.zeros(index.shape, out_dtype) 454 output[index_order] = temp 455 if as_scalar: 456 output = output[0] 457 458 return output 459 460 461def _safely_castable_to_int(dt): 462 """Test whether the NumPy data type `dt` can be safely cast to an int.""" 463 int_size = np.dtype(int).itemsize 464 safe = ((np.issubdtype(dt, np.signedinteger) and dt.itemsize <= int_size) or 465 (np.issubdtype(dt, np.unsignedinteger) and dt.itemsize < int_size)) 466 return safe 467 468 469def _stats(input, labels=None, index=None, centered=False): 470 """Count, sum, and optionally compute (sum - centre)^2 of input by label 471 472 Parameters 473 ---------- 474 input : array_like, N-D 475 The input data to be analyzed. 476 labels : array_like (N-D), optional 477 The labels of the data in `input`. This array must be broadcast 478 compatible with `input`; typically, it is the same shape as `input`. 479 If `labels` is None, all nonzero values in `input` are treated as 480 the single labeled group. 481 index : label or sequence of labels, optional 482 These are the labels of the groups for which the stats are computed. 483 If `index` is None, the stats are computed for the single group where 484 `labels` is greater than 0. 485 centered : bool, optional 486 If True, the centered sum of squares for each labeled group is 487 also returned. Default is False. 488 489 Returns 490 ------- 491 counts : int or ndarray of ints 492 The number of elements in each labeled group. 493 sums : scalar or ndarray of scalars 494 The sums of the values in each labeled group. 495 sums_c : scalar or ndarray of scalars, optional 496 The sums of mean-centered squares of the values in each labeled group. 497 This is only returned if `centered` is True. 498 499 """ 500 def single_group(vals): 501 if centered: 502 vals_c = vals - vals.mean() 503 return vals.size, vals.sum(), (vals_c * vals_c.conjugate()).sum() 504 else: 505 return vals.size, vals.sum() 506 507 if labels is None: 508 return single_group(input) 509 510 # ensure input and labels match sizes 511 input, labels = numpy.broadcast_arrays(input, labels) 512 513 if index is None: 514 return single_group(input[labels > 0]) 515 516 if numpy.isscalar(index): 517 return single_group(input[labels == index]) 518 519 def _sum_centered(labels): 520 # `labels` is expected to be an ndarray with the same shape as `input`. 521 # It must contain the label indices (which are not necessarily the labels 522 # themselves). 523 means = sums / counts 524 centered_input = input - means[labels] 525 # bincount expects 1-D inputs, so we ravel the arguments. 526 bc = numpy.bincount(labels.ravel(), 527 weights=(centered_input * 528 centered_input.conjugate()).ravel()) 529 return bc 530 531 # Remap labels to unique integers if necessary, or if the largest 532 # label is larger than the number of values. 533 534 if (not _safely_castable_to_int(labels.dtype) or 535 labels.min() < 0 or labels.max() > labels.size): 536 # Use numpy.unique to generate the label indices. `new_labels` will 537 # be 1-D, but it should be interpreted as the flattened N-D array of 538 # label indices. 539 unique_labels, new_labels = numpy.unique(labels, return_inverse=True) 540 counts = numpy.bincount(new_labels) 541 sums = numpy.bincount(new_labels, weights=input.ravel()) 542 if centered: 543 # Compute the sum of the mean-centered squares. 544 # We must reshape new_labels to the N-D shape of `input` before 545 # passing it _sum_centered. 546 sums_c = _sum_centered(new_labels.reshape(labels.shape)) 547 idxs = numpy.searchsorted(unique_labels, index) 548 # make all of idxs valid 549 idxs[idxs >= unique_labels.size] = 0 550 found = (unique_labels[idxs] == index) 551 else: 552 # labels are an integer type allowed by bincount, and there aren't too 553 # many, so call bincount directly. 554 counts = numpy.bincount(labels.ravel()) 555 sums = numpy.bincount(labels.ravel(), weights=input.ravel()) 556 if centered: 557 sums_c = _sum_centered(labels) 558 # make sure all index values are valid 559 idxs = numpy.asanyarray(index, numpy.int_).copy() 560 found = (idxs >= 0) & (idxs < counts.size) 561 idxs[~found] = 0 562 563 counts = counts[idxs] 564 counts[~found] = 0 565 sums = sums[idxs] 566 sums[~found] = 0 567 568 if not centered: 569 return (counts, sums) 570 else: 571 sums_c = sums_c[idxs] 572 sums_c[~found] = 0 573 return (counts, sums, sums_c) 574 575 576def sum(input, labels=None, index=None): 577 """ 578 Calculate the sum of the values of the array. 579 580 Notes 581 ----- 582 This is an alias for `ndimage.sum_labels` kept for backwards compatibility 583 reasons, for new code please prefer `sum_labels`. See the `sum_labels` 584 docstring for more details. 585 586 """ 587 return sum_labels(input, labels, index) 588 589 590def sum_labels(input, labels=None, index=None): 591 """ 592 Calculate the sum of the values of the array. 593 594 Parameters 595 ---------- 596 input : array_like 597 Values of `input` inside the regions defined by `labels` 598 are summed together. 599 labels : array_like of ints, optional 600 Assign labels to the values of the array. Has to have the same shape as 601 `input`. 602 index : array_like, optional 603 A single label number or a sequence of label numbers of 604 the objects to be measured. 605 606 Returns 607 ------- 608 sum : ndarray or scalar 609 An array of the sums of values of `input` inside the regions defined 610 by `labels` with the same shape as `index`. If 'index' is None or scalar, 611 a scalar is returned. 612 613 See also 614 -------- 615 mean, median 616 617 Examples 618 -------- 619 >>> from scipy import ndimage 620 >>> input = [0,1,2,3] 621 >>> labels = [1,1,2,2] 622 >>> ndimage.sum(input, labels, index=[1,2]) 623 [1.0, 5.0] 624 >>> ndimage.sum(input, labels, index=1) 625 1 626 >>> ndimage.sum(input, labels) 627 6 628 629 630 """ 631 count, sum = _stats(input, labels, index) 632 return sum 633 634 635def mean(input, labels=None, index=None): 636 """ 637 Calculate the mean of the values of an array at labels. 638 639 Parameters 640 ---------- 641 input : array_like 642 Array on which to compute the mean of elements over distinct 643 regions. 644 labels : array_like, optional 645 Array of labels of same shape, or broadcastable to the same shape as 646 `input`. All elements sharing the same label form one region over 647 which the mean of the elements is computed. 648 index : int or sequence of ints, optional 649 Labels of the objects over which the mean is to be computed. 650 Default is None, in which case the mean for all values where label is 651 greater than 0 is calculated. 652 653 Returns 654 ------- 655 out : list 656 Sequence of same length as `index`, with the mean of the different 657 regions labeled by the labels in `index`. 658 659 See also 660 -------- 661 variance, standard_deviation, minimum, maximum, sum, label 662 663 Examples 664 -------- 665 >>> from scipy import ndimage 666 >>> a = np.arange(25).reshape((5,5)) 667 >>> labels = np.zeros_like(a) 668 >>> labels[3:5,3:5] = 1 669 >>> index = np.unique(labels) 670 >>> labels 671 array([[0, 0, 0, 0, 0], 672 [0, 0, 0, 0, 0], 673 [0, 0, 0, 0, 0], 674 [0, 0, 0, 1, 1], 675 [0, 0, 0, 1, 1]]) 676 >>> index 677 array([0, 1]) 678 >>> ndimage.mean(a, labels=labels, index=index) 679 [10.285714285714286, 21.0] 680 681 """ 682 683 count, sum = _stats(input, labels, index) 684 return sum / numpy.asanyarray(count).astype(numpy.float64) 685 686 687def variance(input, labels=None, index=None): 688 """ 689 Calculate the variance of the values of an N-D image array, optionally at 690 specified sub-regions. 691 692 Parameters 693 ---------- 694 input : array_like 695 Nd-image data to process. 696 labels : array_like, optional 697 Labels defining sub-regions in `input`. 698 If not None, must be same shape as `input`. 699 index : int or sequence of ints, optional 700 `labels` to include in output. If None (default), all values where 701 `labels` is non-zero are used. 702 703 Returns 704 ------- 705 variance : float or ndarray 706 Values of variance, for each sub-region if `labels` and `index` are 707 specified. 708 709 See Also 710 -------- 711 label, standard_deviation, maximum, minimum, extrema 712 713 Examples 714 -------- 715 >>> a = np.array([[1, 2, 0, 0], 716 ... [5, 3, 0, 4], 717 ... [0, 0, 0, 7], 718 ... [9, 3, 0, 0]]) 719 >>> from scipy import ndimage 720 >>> ndimage.variance(a) 721 7.609375 722 723 Features to process can be specified using `labels` and `index`: 724 725 >>> lbl, nlbl = ndimage.label(a) 726 >>> ndimage.variance(a, lbl, index=np.arange(1, nlbl+1)) 727 array([ 2.1875, 2.25 , 9. ]) 728 729 If no index is given, all non-zero `labels` are processed: 730 731 >>> ndimage.variance(a, lbl) 732 6.1875 733 734 """ 735 count, sum, sum_c_sq = _stats(input, labels, index, centered=True) 736 return sum_c_sq / np.asanyarray(count).astype(float) 737 738 739def standard_deviation(input, labels=None, index=None): 740 """ 741 Calculate the standard deviation of the values of an N-D image array, 742 optionally at specified sub-regions. 743 744 Parameters 745 ---------- 746 input : array_like 747 N-D image data to process. 748 labels : array_like, optional 749 Labels to identify sub-regions in `input`. 750 If not None, must be same shape as `input`. 751 index : int or sequence of ints, optional 752 `labels` to include in output. If None (default), all values where 753 `labels` is non-zero are used. 754 755 Returns 756 ------- 757 standard_deviation : float or ndarray 758 Values of standard deviation, for each sub-region if `labels` and 759 `index` are specified. 760 761 See Also 762 -------- 763 label, variance, maximum, minimum, extrema 764 765 Examples 766 -------- 767 >>> a = np.array([[1, 2, 0, 0], 768 ... [5, 3, 0, 4], 769 ... [0, 0, 0, 7], 770 ... [9, 3, 0, 0]]) 771 >>> from scipy import ndimage 772 >>> ndimage.standard_deviation(a) 773 2.7585095613392387 774 775 Features to process can be specified using `labels` and `index`: 776 777 >>> lbl, nlbl = ndimage.label(a) 778 >>> ndimage.standard_deviation(a, lbl, index=np.arange(1, nlbl+1)) 779 array([ 1.479, 1.5 , 3. ]) 780 781 If no index is given, non-zero `labels` are processed: 782 783 >>> ndimage.standard_deviation(a, lbl) 784 2.4874685927665499 785 786 """ 787 return numpy.sqrt(variance(input, labels, index)) 788 789 790def _select(input, labels=None, index=None, find_min=False, find_max=False, 791 find_min_positions=False, find_max_positions=False, 792 find_median=False): 793 """Returns min, max, or both, plus their positions (if requested), and 794 median.""" 795 796 input = numpy.asanyarray(input) 797 798 find_positions = find_min_positions or find_max_positions 799 positions = None 800 if find_positions: 801 positions = numpy.arange(input.size).reshape(input.shape) 802 803 def single_group(vals, positions): 804 result = [] 805 if find_min: 806 result += [vals.min()] 807 if find_min_positions: 808 result += [positions[vals == vals.min()][0]] 809 if find_max: 810 result += [vals.max()] 811 if find_max_positions: 812 result += [positions[vals == vals.max()][0]] 813 if find_median: 814 result += [numpy.median(vals)] 815 return result 816 817 if labels is None: 818 return single_group(input, positions) 819 820 # ensure input and labels match sizes 821 input, labels = numpy.broadcast_arrays(input, labels) 822 823 if index is None: 824 mask = (labels > 0) 825 masked_positions = None 826 if find_positions: 827 masked_positions = positions[mask] 828 return single_group(input[mask], masked_positions) 829 830 if numpy.isscalar(index): 831 mask = (labels == index) 832 masked_positions = None 833 if find_positions: 834 masked_positions = positions[mask] 835 return single_group(input[mask], masked_positions) 836 837 # remap labels to unique integers if necessary, or if the largest 838 # label is larger than the number of values. 839 if (not _safely_castable_to_int(labels.dtype) or 840 labels.min() < 0 or labels.max() > labels.size): 841 # remap labels, and indexes 842 unique_labels, labels = numpy.unique(labels, return_inverse=True) 843 idxs = numpy.searchsorted(unique_labels, index) 844 845 # make all of idxs valid 846 idxs[idxs >= unique_labels.size] = 0 847 found = (unique_labels[idxs] == index) 848 else: 849 # labels are an integer type, and there aren't too many 850 idxs = numpy.asanyarray(index, numpy.int_).copy() 851 found = (idxs >= 0) & (idxs <= labels.max()) 852 853 idxs[~ found] = labels.max() + 1 854 855 if find_median: 856 order = numpy.lexsort((input.ravel(), labels.ravel())) 857 else: 858 order = input.ravel().argsort() 859 input = input.ravel()[order] 860 labels = labels.ravel()[order] 861 if find_positions: 862 positions = positions.ravel()[order] 863 864 result = [] 865 if find_min: 866 mins = numpy.zeros(labels.max() + 2, input.dtype) 867 mins[labels[::-1]] = input[::-1] 868 result += [mins[idxs]] 869 if find_min_positions: 870 minpos = numpy.zeros(labels.max() + 2, int) 871 minpos[labels[::-1]] = positions[::-1] 872 result += [minpos[idxs]] 873 if find_max: 874 maxs = numpy.zeros(labels.max() + 2, input.dtype) 875 maxs[labels] = input 876 result += [maxs[idxs]] 877 if find_max_positions: 878 maxpos = numpy.zeros(labels.max() + 2, int) 879 maxpos[labels] = positions 880 result += [maxpos[idxs]] 881 if find_median: 882 locs = numpy.arange(len(labels)) 883 lo = numpy.zeros(labels.max() + 2, numpy.int_) 884 lo[labels[::-1]] = locs[::-1] 885 hi = numpy.zeros(labels.max() + 2, numpy.int_) 886 hi[labels] = locs 887 lo = lo[idxs] 888 hi = hi[idxs] 889 # lo is an index to the lowest value in input for each label, 890 # hi is an index to the largest value. 891 # move them to be either the same ((hi - lo) % 2 == 0) or next 892 # to each other ((hi - lo) % 2 == 1), then average. 893 step = (hi - lo) // 2 894 lo += step 895 hi -= step 896 if (np.issubdtype(input.dtype, np.integer) 897 or np.issubdtype(input.dtype, np.bool_)): 898 # avoid integer overflow or boolean addition (gh-12836) 899 result += [(input[lo].astype('d') + input[hi].astype('d')) / 2.0] 900 else: 901 result += [(input[lo] + input[hi]) / 2.0] 902 903 return result 904 905 906def minimum(input, labels=None, index=None): 907 """ 908 Calculate the minimum of the values of an array over labeled regions. 909 910 Parameters 911 ---------- 912 input : array_like 913 Array_like of values. For each region specified by `labels`, the 914 minimal values of `input` over the region is computed. 915 labels : array_like, optional 916 An array_like of integers marking different regions over which the 917 minimum value of `input` is to be computed. `labels` must have the 918 same shape as `input`. If `labels` is not specified, the minimum 919 over the whole array is returned. 920 index : array_like, optional 921 A list of region labels that are taken into account for computing the 922 minima. If index is None, the minimum over all elements where `labels` 923 is non-zero is returned. 924 925 Returns 926 ------- 927 minimum : float or list of floats 928 List of minima of `input` over the regions determined by `labels` and 929 whose index is in `index`. If `index` or `labels` are not specified, a 930 float is returned: the minimal value of `input` if `labels` is None, 931 and the minimal value of elements where `labels` is greater than zero 932 if `index` is None. 933 934 See also 935 -------- 936 label, maximum, median, minimum_position, extrema, sum, mean, variance, 937 standard_deviation 938 939 Notes 940 ----- 941 The function returns a Python list and not a NumPy array, use 942 `np.array` to convert the list to an array. 943 944 Examples 945 -------- 946 >>> from scipy import ndimage 947 >>> a = np.array([[1, 2, 0, 0], 948 ... [5, 3, 0, 4], 949 ... [0, 0, 0, 7], 950 ... [9, 3, 0, 0]]) 951 >>> labels, labels_nb = ndimage.label(a) 952 >>> labels 953 array([[1, 1, 0, 0], 954 [1, 1, 0, 2], 955 [0, 0, 0, 2], 956 [3, 3, 0, 0]]) 957 >>> ndimage.minimum(a, labels=labels, index=np.arange(1, labels_nb + 1)) 958 [1.0, 4.0, 3.0] 959 >>> ndimage.minimum(a) 960 0.0 961 >>> ndimage.minimum(a, labels=labels) 962 1.0 963 964 """ 965 return _select(input, labels, index, find_min=True)[0] 966 967 968def maximum(input, labels=None, index=None): 969 """ 970 Calculate the maximum of the values of an array over labeled regions. 971 972 Parameters 973 ---------- 974 input : array_like 975 Array_like of values. For each region specified by `labels`, the 976 maximal values of `input` over the region is computed. 977 labels : array_like, optional 978 An array of integers marking different regions over which the 979 maximum value of `input` is to be computed. `labels` must have the 980 same shape as `input`. If `labels` is not specified, the maximum 981 over the whole array is returned. 982 index : array_like, optional 983 A list of region labels that are taken into account for computing the 984 maxima. If index is None, the maximum over all elements where `labels` 985 is non-zero is returned. 986 987 Returns 988 ------- 989 output : float or list of floats 990 List of maxima of `input` over the regions determined by `labels` and 991 whose index is in `index`. If `index` or `labels` are not specified, a 992 float is returned: the maximal value of `input` if `labels` is None, 993 and the maximal value of elements where `labels` is greater than zero 994 if `index` is None. 995 996 See also 997 -------- 998 label, minimum, median, maximum_position, extrema, sum, mean, variance, 999 standard_deviation 1000 1001 Notes 1002 ----- 1003 The function returns a Python list and not a NumPy array, use 1004 `np.array` to convert the list to an array. 1005 1006 Examples 1007 -------- 1008 >>> a = np.arange(16).reshape((4,4)) 1009 >>> a 1010 array([[ 0, 1, 2, 3], 1011 [ 4, 5, 6, 7], 1012 [ 8, 9, 10, 11], 1013 [12, 13, 14, 15]]) 1014 >>> labels = np.zeros_like(a) 1015 >>> labels[:2,:2] = 1 1016 >>> labels[2:, 1:3] = 2 1017 >>> labels 1018 array([[1, 1, 0, 0], 1019 [1, 1, 0, 0], 1020 [0, 2, 2, 0], 1021 [0, 2, 2, 0]]) 1022 >>> from scipy import ndimage 1023 >>> ndimage.maximum(a) 1024 15.0 1025 >>> ndimage.maximum(a, labels=labels, index=[1,2]) 1026 [5.0, 14.0] 1027 >>> ndimage.maximum(a, labels=labels) 1028 14.0 1029 1030 >>> b = np.array([[1, 2, 0, 0], 1031 ... [5, 3, 0, 4], 1032 ... [0, 0, 0, 7], 1033 ... [9, 3, 0, 0]]) 1034 >>> labels, labels_nb = ndimage.label(b) 1035 >>> labels 1036 array([[1, 1, 0, 0], 1037 [1, 1, 0, 2], 1038 [0, 0, 0, 2], 1039 [3, 3, 0, 0]]) 1040 >>> ndimage.maximum(b, labels=labels, index=np.arange(1, labels_nb + 1)) 1041 [5.0, 7.0, 9.0] 1042 1043 """ 1044 return _select(input, labels, index, find_max=True)[0] 1045 1046 1047def median(input, labels=None, index=None): 1048 """ 1049 Calculate the median of the values of an array over labeled regions. 1050 1051 Parameters 1052 ---------- 1053 input : array_like 1054 Array_like of values. For each region specified by `labels`, the 1055 median value of `input` over the region is computed. 1056 labels : array_like, optional 1057 An array_like of integers marking different regions over which the 1058 median value of `input` is to be computed. `labels` must have the 1059 same shape as `input`. If `labels` is not specified, the median 1060 over the whole array is returned. 1061 index : array_like, optional 1062 A list of region labels that are taken into account for computing the 1063 medians. If index is None, the median over all elements where `labels` 1064 is non-zero is returned. 1065 1066 Returns 1067 ------- 1068 median : float or list of floats 1069 List of medians of `input` over the regions determined by `labels` and 1070 whose index is in `index`. If `index` or `labels` are not specified, a 1071 float is returned: the median value of `input` if `labels` is None, 1072 and the median value of elements where `labels` is greater than zero 1073 if `index` is None. 1074 1075 See also 1076 -------- 1077 label, minimum, maximum, extrema, sum, mean, variance, standard_deviation 1078 1079 Notes 1080 ----- 1081 The function returns a Python list and not a NumPy array, use 1082 `np.array` to convert the list to an array. 1083 1084 Examples 1085 -------- 1086 >>> from scipy import ndimage 1087 >>> a = np.array([[1, 2, 0, 1], 1088 ... [5, 3, 0, 4], 1089 ... [0, 0, 0, 7], 1090 ... [9, 3, 0, 0]]) 1091 >>> labels, labels_nb = ndimage.label(a) 1092 >>> labels 1093 array([[1, 1, 0, 2], 1094 [1, 1, 0, 2], 1095 [0, 0, 0, 2], 1096 [3, 3, 0, 0]]) 1097 >>> ndimage.median(a, labels=labels, index=np.arange(1, labels_nb + 1)) 1098 [2.5, 4.0, 6.0] 1099 >>> ndimage.median(a) 1100 1.0 1101 >>> ndimage.median(a, labels=labels) 1102 3.0 1103 1104 """ 1105 return _select(input, labels, index, find_median=True)[0] 1106 1107 1108def minimum_position(input, labels=None, index=None): 1109 """ 1110 Find the positions of the minimums of the values of an array at labels. 1111 1112 Parameters 1113 ---------- 1114 input : array_like 1115 Array_like of values. 1116 labels : array_like, optional 1117 An array of integers marking different regions over which the 1118 position of the minimum value of `input` is to be computed. 1119 `labels` must have the same shape as `input`. If `labels` is not 1120 specified, the location of the first minimum over the whole 1121 array is returned. 1122 1123 The `labels` argument only works when `index` is specified. 1124 index : array_like, optional 1125 A list of region labels that are taken into account for finding the 1126 location of the minima. If `index` is None, the ``first`` minimum 1127 over all elements where `labels` is non-zero is returned. 1128 1129 The `index` argument only works when `labels` is specified. 1130 1131 Returns 1132 ------- 1133 output : list of tuples of ints 1134 Tuple of ints or list of tuples of ints that specify the location 1135 of minima of `input` over the regions determined by `labels` and 1136 whose index is in `index`. 1137 1138 If `index` or `labels` are not specified, a tuple of ints is 1139 returned specifying the location of the first minimal value of `input`. 1140 1141 See also 1142 -------- 1143 label, minimum, median, maximum_position, extrema, sum, mean, variance, 1144 standard_deviation 1145 1146 Examples 1147 -------- 1148 >>> a = np.array([[10, 20, 30], 1149 ... [40, 80, 100], 1150 ... [1, 100, 200]]) 1151 >>> b = np.array([[1, 2, 0, 1], 1152 ... [5, 3, 0, 4], 1153 ... [0, 0, 0, 7], 1154 ... [9, 3, 0, 0]]) 1155 1156 >>> from scipy import ndimage 1157 1158 >>> ndimage.minimum_position(a) 1159 (2, 0) 1160 >>> ndimage.minimum_position(b) 1161 (0, 2) 1162 1163 Features to process can be specified using `labels` and `index`: 1164 1165 >>> label, pos = ndimage.label(a) 1166 >>> ndimage.minimum_position(a, label, index=np.arange(1, pos+1)) 1167 [(2, 0)] 1168 1169 >>> label, pos = ndimage.label(b) 1170 >>> ndimage.minimum_position(b, label, index=np.arange(1, pos+1)) 1171 [(0, 0), (0, 3), (3, 1)] 1172 1173 """ 1174 dims = numpy.array(numpy.asarray(input).shape) 1175 # see numpy.unravel_index to understand this line. 1176 dim_prod = numpy.cumprod([1] + list(dims[:0:-1]))[::-1] 1177 1178 result = _select(input, labels, index, find_min_positions=True)[0] 1179 1180 if numpy.isscalar(result): 1181 return tuple((result // dim_prod) % dims) 1182 1183 return [tuple(v) for v in (result.reshape(-1, 1) // dim_prod) % dims] 1184 1185 1186def maximum_position(input, labels=None, index=None): 1187 """ 1188 Find the positions of the maximums of the values of an array at labels. 1189 1190 For each region specified by `labels`, the position of the maximum 1191 value of `input` within the region is returned. 1192 1193 Parameters 1194 ---------- 1195 input : array_like 1196 Array_like of values. 1197 labels : array_like, optional 1198 An array of integers marking different regions over which the 1199 position of the maximum value of `input` is to be computed. 1200 `labels` must have the same shape as `input`. If `labels` is not 1201 specified, the location of the first maximum over the whole 1202 array is returned. 1203 1204 The `labels` argument only works when `index` is specified. 1205 index : array_like, optional 1206 A list of region labels that are taken into account for finding the 1207 location of the maxima. If `index` is None, the first maximum 1208 over all elements where `labels` is non-zero is returned. 1209 1210 The `index` argument only works when `labels` is specified. 1211 1212 Returns 1213 ------- 1214 output : list of tuples of ints 1215 List of tuples of ints that specify the location of maxima of 1216 `input` over the regions determined by `labels` and whose index 1217 is in `index`. 1218 1219 If `index` or `labels` are not specified, a tuple of ints is 1220 returned specifying the location of the ``first`` maximal value 1221 of `input`. 1222 1223 See also 1224 -------- 1225 label, minimum, median, maximum_position, extrema, sum, mean, variance, 1226 standard_deviation 1227 1228 Examples 1229 -------- 1230 >>> from scipy import ndimage 1231 >>> a = np.array([[1, 2, 0, 0], 1232 ... [5, 3, 0, 4], 1233 ... [0, 0, 0, 7], 1234 ... [9, 3, 0, 0]]) 1235 >>> ndimage.maximum_position(a) 1236 (3, 0) 1237 1238 Features to process can be specified using `labels` and `index`: 1239 1240 >>> lbl = np.array([[0, 1, 2, 3], 1241 ... [0, 1, 2, 3], 1242 ... [0, 1, 2, 3], 1243 ... [0, 1, 2, 3]]) 1244 >>> ndimage.maximum_position(a, lbl, 1) 1245 (1, 1) 1246 1247 If no index is given, non-zero `labels` are processed: 1248 1249 >>> ndimage.maximum_position(a, lbl) 1250 (2, 3) 1251 1252 If there are no maxima, the position of the first element is returned: 1253 1254 >>> ndimage.maximum_position(a, lbl, 2) 1255 (0, 2) 1256 1257 """ 1258 dims = numpy.array(numpy.asarray(input).shape) 1259 # see numpy.unravel_index to understand this line. 1260 dim_prod = numpy.cumprod([1] + list(dims[:0:-1]))[::-1] 1261 1262 result = _select(input, labels, index, find_max_positions=True)[0] 1263 1264 if numpy.isscalar(result): 1265 return tuple((result // dim_prod) % dims) 1266 1267 return [tuple(v) for v in (result.reshape(-1, 1) // dim_prod) % dims] 1268 1269 1270def extrema(input, labels=None, index=None): 1271 """ 1272 Calculate the minimums and maximums of the values of an array 1273 at labels, along with their positions. 1274 1275 Parameters 1276 ---------- 1277 input : ndarray 1278 N-D image data to process. 1279 labels : ndarray, optional 1280 Labels of features in input. 1281 If not None, must be same shape as `input`. 1282 index : int or sequence of ints, optional 1283 Labels to include in output. If None (default), all values where 1284 non-zero `labels` are used. 1285 1286 Returns 1287 ------- 1288 minimums, maximums : int or ndarray 1289 Values of minimums and maximums in each feature. 1290 min_positions, max_positions : tuple or list of tuples 1291 Each tuple gives the N-D coordinates of the corresponding minimum 1292 or maximum. 1293 1294 See Also 1295 -------- 1296 maximum, minimum, maximum_position, minimum_position, center_of_mass 1297 1298 Examples 1299 -------- 1300 >>> a = np.array([[1, 2, 0, 0], 1301 ... [5, 3, 0, 4], 1302 ... [0, 0, 0, 7], 1303 ... [9, 3, 0, 0]]) 1304 >>> from scipy import ndimage 1305 >>> ndimage.extrema(a) 1306 (0, 9, (0, 2), (3, 0)) 1307 1308 Features to process can be specified using `labels` and `index`: 1309 1310 >>> lbl, nlbl = ndimage.label(a) 1311 >>> ndimage.extrema(a, lbl, index=np.arange(1, nlbl+1)) 1312 (array([1, 4, 3]), 1313 array([5, 7, 9]), 1314 [(0, 0), (1, 3), (3, 1)], 1315 [(1, 0), (2, 3), (3, 0)]) 1316 1317 If no index is given, non-zero `labels` are processed: 1318 1319 >>> ndimage.extrema(a, lbl) 1320 (1, 9, (0, 0), (3, 0)) 1321 1322 """ 1323 dims = numpy.array(numpy.asarray(input).shape) 1324 # see numpy.unravel_index to understand this line. 1325 dim_prod = numpy.cumprod([1] + list(dims[:0:-1]))[::-1] 1326 1327 minimums, min_positions, maximums, max_positions = _select(input, labels, 1328 index, 1329 find_min=True, 1330 find_max=True, 1331 find_min_positions=True, 1332 find_max_positions=True) 1333 1334 if numpy.isscalar(minimums): 1335 return (minimums, maximums, tuple((min_positions // dim_prod) % dims), 1336 tuple((max_positions // dim_prod) % dims)) 1337 1338 min_positions = [tuple(v) for v in (min_positions.reshape(-1, 1) // dim_prod) % dims] 1339 max_positions = [tuple(v) for v in (max_positions.reshape(-1, 1) // dim_prod) % dims] 1340 1341 return minimums, maximums, min_positions, max_positions 1342 1343 1344def center_of_mass(input, labels=None, index=None): 1345 """ 1346 Calculate the center of mass of the values of an array at labels. 1347 1348 Parameters 1349 ---------- 1350 input : ndarray 1351 Data from which to calculate center-of-mass. The masses can either 1352 be positive or negative. 1353 labels : ndarray, optional 1354 Labels for objects in `input`, as generated by `ndimage.label`. 1355 Only used with `index`. Dimensions must be the same as `input`. 1356 index : int or sequence of ints, optional 1357 Labels for which to calculate centers-of-mass. If not specified, 1358 all labels greater than zero are used. Only used with `labels`. 1359 1360 Returns 1361 ------- 1362 center_of_mass : tuple, or list of tuples 1363 Coordinates of centers-of-mass. 1364 1365 Examples 1366 -------- 1367 >>> a = np.array(([0,0,0,0], 1368 ... [0,1,1,0], 1369 ... [0,1,1,0], 1370 ... [0,1,1,0])) 1371 >>> from scipy import ndimage 1372 >>> ndimage.measurements.center_of_mass(a) 1373 (2.0, 1.5) 1374 1375 Calculation of multiple objects in an image 1376 1377 >>> b = np.array(([0,1,1,0], 1378 ... [0,1,0,0], 1379 ... [0,0,0,0], 1380 ... [0,0,1,1], 1381 ... [0,0,1,1])) 1382 >>> lbl = ndimage.label(b)[0] 1383 >>> ndimage.measurements.center_of_mass(b, lbl, [1,2]) 1384 [(0.33333333333333331, 1.3333333333333333), (3.5, 2.5)] 1385 1386 Negative masses are also accepted, which can occur for example when 1387 bias is removed from measured data due to random noise. 1388 1389 >>> c = np.array(([-1,0,0,0], 1390 ... [0,-1,-1,0], 1391 ... [0,1,-1,0], 1392 ... [0,1,1,0])) 1393 >>> ndimage.measurements.center_of_mass(c) 1394 (-4.0, 1.0) 1395 1396 If there are division by zero issues, the function does not raise an 1397 error but rather issues a RuntimeWarning before returning inf and/or NaN. 1398 1399 >>> d = np.array([-1, 1]) 1400 >>> ndimage.measurements.center_of_mass(d) 1401 (inf,) 1402 """ 1403 normalizer = sum(input, labels, index) 1404 grids = numpy.ogrid[[slice(0, i) for i in input.shape]] 1405 1406 results = [sum(input * grids[dir].astype(float), labels, index) / normalizer 1407 for dir in range(input.ndim)] 1408 1409 if numpy.isscalar(results[0]): 1410 return tuple(results) 1411 1412 return [tuple(v) for v in numpy.array(results).T] 1413 1414 1415def histogram(input, min, max, bins, labels=None, index=None): 1416 """ 1417 Calculate the histogram of the values of an array, optionally at labels. 1418 1419 Histogram calculates the frequency of values in an array within bins 1420 determined by `min`, `max`, and `bins`. The `labels` and `index` 1421 keywords can limit the scope of the histogram to specified sub-regions 1422 within the array. 1423 1424 Parameters 1425 ---------- 1426 input : array_like 1427 Data for which to calculate histogram. 1428 min, max : int 1429 Minimum and maximum values of range of histogram bins. 1430 bins : int 1431 Number of bins. 1432 labels : array_like, optional 1433 Labels for objects in `input`. 1434 If not None, must be same shape as `input`. 1435 index : int or sequence of ints, optional 1436 Label or labels for which to calculate histogram. If None, all values 1437 where label is greater than zero are used 1438 1439 Returns 1440 ------- 1441 hist : ndarray 1442 Histogram counts. 1443 1444 Examples 1445 -------- 1446 >>> a = np.array([[ 0. , 0.2146, 0.5962, 0. ], 1447 ... [ 0. , 0.7778, 0. , 0. ], 1448 ... [ 0. , 0. , 0. , 0. ], 1449 ... [ 0. , 0. , 0.7181, 0.2787], 1450 ... [ 0. , 0. , 0.6573, 0.3094]]) 1451 >>> from scipy import ndimage 1452 >>> ndimage.measurements.histogram(a, 0, 1, 10) 1453 array([13, 0, 2, 1, 0, 1, 1, 2, 0, 0]) 1454 1455 With labels and no indices, non-zero elements are counted: 1456 1457 >>> lbl, nlbl = ndimage.label(a) 1458 >>> ndimage.measurements.histogram(a, 0, 1, 10, lbl) 1459 array([0, 0, 2, 1, 0, 1, 1, 2, 0, 0]) 1460 1461 Indices can be used to count only certain objects: 1462 1463 >>> ndimage.measurements.histogram(a, 0, 1, 10, lbl, 2) 1464 array([0, 0, 1, 1, 0, 0, 1, 1, 0, 0]) 1465 1466 """ 1467 _bins = numpy.linspace(min, max, bins + 1) 1468 1469 def _hist(vals): 1470 return numpy.histogram(vals, _bins)[0] 1471 1472 return labeled_comprehension(input, labels, index, _hist, object, None, 1473 pass_positions=False) 1474 1475 1476def watershed_ift(input, markers, structure=None, output=None): 1477 """ 1478 Apply watershed from markers using image foresting transform algorithm. 1479 1480 Parameters 1481 ---------- 1482 input : array_like 1483 Input. 1484 markers : array_like 1485 Markers are points within each watershed that form the beginning 1486 of the process. Negative markers are considered background markers 1487 which are processed after the other markers. 1488 structure : structure element, optional 1489 A structuring element defining the connectivity of the object can be 1490 provided. If None, an element is generated with a squared 1491 connectivity equal to one. 1492 output : ndarray, optional 1493 An output array can optionally be provided. The same shape as input. 1494 1495 Returns 1496 ------- 1497 watershed_ift : ndarray 1498 Output. Same shape as `input`. 1499 1500 References 1501 ---------- 1502 .. [1] A.X. Falcao, J. Stolfi and R. de Alencar Lotufo, "The image 1503 foresting transform: theory, algorithms, and applications", 1504 Pattern Analysis and Machine Intelligence, vol. 26, pp. 19-29, 2004. 1505 1506 """ 1507 input = numpy.asarray(input) 1508 if input.dtype.type not in [numpy.uint8, numpy.uint16]: 1509 raise TypeError('only 8 and 16 unsigned inputs are supported') 1510 1511 if structure is None: 1512 structure = morphology.generate_binary_structure(input.ndim, 1) 1513 structure = numpy.asarray(structure, dtype=bool) 1514 if structure.ndim != input.ndim: 1515 raise RuntimeError('structure and input must have equal rank') 1516 for ii in structure.shape: 1517 if ii != 3: 1518 raise RuntimeError('structure dimensions must be equal to 3') 1519 1520 if not structure.flags.contiguous: 1521 structure = structure.copy() 1522 markers = numpy.asarray(markers) 1523 if input.shape != markers.shape: 1524 raise RuntimeError('input and markers must have equal shape') 1525 1526 integral_types = [numpy.int0, 1527 numpy.int8, 1528 numpy.int16, 1529 numpy.int32, 1530 numpy.int_, 1531 numpy.int64, 1532 numpy.intc, 1533 numpy.intp] 1534 1535 if markers.dtype.type not in integral_types: 1536 raise RuntimeError('marker should be of integer type') 1537 1538 if isinstance(output, numpy.ndarray): 1539 if output.dtype.type not in integral_types: 1540 raise RuntimeError('output should be of integer type') 1541 else: 1542 output = markers.dtype 1543 1544 output = _ni_support._get_output(output, input) 1545 _nd_image.watershed_ift(input, markers, structure, output) 1546 return output 1547