Coverage for /usr/lib/python3/dist-packages/scipy/ndimage/_measurements.py: 9%

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', 'value_indices'] 

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 >>> import numpy as np 

126 >>> a = np.array([[0,0,1,1,0,0], 

127 ... [0,0,0,1,0,0], 

128 ... [1,1,0,0,1,0], 

129 ... [0,0,0,1,0,0]]) 

130 >>> labeled_array, num_features = label(a) 

131 

132 Each of the 4 features are labeled with a different integer: 

133 

134 >>> num_features 

135 4 

136 >>> labeled_array 

137 array([[0, 0, 1, 1, 0, 0], 

138 [0, 0, 0, 1, 0, 0], 

139 [2, 2, 0, 0, 3, 0], 

140 [0, 0, 0, 4, 0, 0]]) 

141 

142 Generate a structuring element that will consider features connected even 

143 if they touch diagonally: 

144 

145 >>> s = generate_binary_structure(2,2) 

146 

147 or, 

148 

149 >>> s = [[1,1,1], 

150 ... [1,1,1], 

151 ... [1,1,1]] 

152 

153 Label the image using the new structuring element: 

154 

155 >>> labeled_array, num_features = label(a, structure=s) 

156 

157 Show the 2 labeled features (note that features 1, 3, and 4 from above are 

158 now considered a single feature): 

159 

160 >>> num_features 

161 2 

162 >>> labeled_array 

163 array([[0, 0, 1, 1, 0, 0], 

164 [0, 0, 0, 1, 0, 0], 

165 [2, 2, 0, 0, 1, 0], 

166 [0, 0, 0, 1, 0, 0]]) 

167 

168 References 

169 ---------- 

170 

171 .. [1] James R. Weaver, "Centrosymmetric (cross-symmetric) 

172 matrices, their basic properties, eigenvalues, and 

173 eigenvectors." The American Mathematical Monthly 92.10 

174 (1985): 711-717. 

175 

176 """ 

177 input = numpy.asarray(input) 

178 if numpy.iscomplexobj(input): 

179 raise TypeError('Complex type not supported') 

180 if structure is None: 

181 structure = _morphology.generate_binary_structure(input.ndim, 1) 

182 structure = numpy.asarray(structure, dtype=bool) 

183 if structure.ndim != input.ndim: 

184 raise RuntimeError('structure and input must have equal rank') 

185 for ii in structure.shape: 

186 if ii != 3: 

187 raise ValueError('structure dimensions must be equal to 3') 

188 

189 # Use 32 bits if it's large enough for this image. 

190 # _ni_label.label() needs two entries for background and 

191 # foreground tracking 

192 need_64bits = input.size >= (2**31 - 2) 

193 

194 if isinstance(output, numpy.ndarray): 

195 if output.shape != input.shape: 

196 raise ValueError("output shape not correct") 

197 caller_provided_output = True 

198 else: 

199 caller_provided_output = False 

200 if output is None: 

201 output = np.empty(input.shape, np.intp if need_64bits else np.int32) 

202 else: 

203 output = np.empty(input.shape, output) 

204 

205 # handle scalars, 0-D arrays 

206 if input.ndim == 0 or input.size == 0: 

207 if input.ndim == 0: 

208 # scalar 

209 maxlabel = 1 if (input != 0) else 0 

210 output[...] = maxlabel 

211 else: 

212 # 0-D 

213 maxlabel = 0 

214 if caller_provided_output: 

215 return maxlabel 

216 else: 

217 return output, maxlabel 

218 

219 try: 

220 max_label = _ni_label._label(input, structure, output) 

221 except _ni_label.NeedMoreBits as e: 

222 # Make another attempt with enough bits, then try to cast to the 

223 # new type. 

224 tmp_output = np.empty(input.shape, np.intp if need_64bits else np.int32) 

225 max_label = _ni_label._label(input, structure, tmp_output) 

226 output[...] = tmp_output[...] 

227 if not np.all(output == tmp_output): 

228 # refuse to return bad results 

229 raise RuntimeError( 

230 "insufficient bit-depth in requested output type" 

231 ) from e 

232 

233 if caller_provided_output: 

234 # result was written in-place 

235 return max_label 

236 else: 

237 return output, max_label 

238 

239 

240def find_objects(input, max_label=0): 

241 """ 

242 Find objects in a labeled array. 

243 

244 Parameters 

245 ---------- 

246 input : ndarray of ints 

247 Array containing objects defined by different labels. Labels with 

248 value 0 are ignored. 

249 max_label : int, optional 

250 Maximum label to be searched for in `input`. If max_label is not 

251 given, the positions of all objects are returned. 

252 

253 Returns 

254 ------- 

255 object_slices : list of tuples 

256 A list of tuples, with each tuple containing N slices (with N the 

257 dimension of the input array). Slices correspond to the minimal 

258 parallelepiped that contains the object. If a number is missing, 

259 None is returned instead of a slice. The label ``l`` corresponds to 

260 the index ``l-1`` in the returned list. 

261 

262 See Also 

263 -------- 

264 label, center_of_mass 

265 

266 Notes 

267 ----- 

268 This function is very useful for isolating a volume of interest inside 

269 a 3-D array, that cannot be "seen through". 

270 

271 Examples 

272 -------- 

273 >>> from scipy import ndimage 

274 >>> import numpy as np 

275 >>> a = np.zeros((6,6), dtype=int) 

276 >>> a[2:4, 2:4] = 1 

277 >>> a[4, 4] = 1 

278 >>> a[:2, :3] = 2 

279 >>> a[0, 5] = 3 

280 >>> a 

281 array([[2, 2, 2, 0, 0, 3], 

282 [2, 2, 2, 0, 0, 0], 

283 [0, 0, 1, 1, 0, 0], 

284 [0, 0, 1, 1, 0, 0], 

285 [0, 0, 0, 0, 1, 0], 

286 [0, 0, 0, 0, 0, 0]]) 

287 >>> ndimage.find_objects(a) 

288 [(slice(2, 5, None), slice(2, 5, None)), (slice(0, 2, None), slice(0, 3, None)), (slice(0, 1, None), slice(5, 6, None))] 

289 >>> ndimage.find_objects(a, max_label=2) 

290 [(slice(2, 5, None), slice(2, 5, None)), (slice(0, 2, None), slice(0, 3, None))] 

291 >>> ndimage.find_objects(a == 1, max_label=2) 

292 [(slice(2, 5, None), slice(2, 5, None)), None] 

293 

294 >>> loc = ndimage.find_objects(a)[0] 

295 >>> a[loc] 

296 array([[1, 1, 0], 

297 [1, 1, 0], 

298 [0, 0, 1]]) 

299 

300 """ 

301 input = numpy.asarray(input) 

302 if numpy.iscomplexobj(input): 

303 raise TypeError('Complex type not supported') 

304 

305 if max_label < 1: 

306 max_label = input.max() 

307 

308 return _nd_image.find_objects(input, max_label) 

309 

310 

311def value_indices(arr, *, ignore_value=None): 

312 """ 

313 Find indices of each distinct value in given array. 

314 

315 Parameters 

316 ---------- 

317 arr : ndarray of ints 

318 Array containing integer values. 

319 ignore_value : int, optional 

320 This value will be ignored in searching the `arr` array. If not 

321 given, all values found will be included in output. Default 

322 is None. 

323 

324 Returns 

325 ------- 

326 indices : dictionary 

327 A Python dictionary of array indices for each distinct value. The 

328 dictionary is keyed by the distinct values, the entries are array 

329 index tuples covering all occurrences of the value within the 

330 array. 

331 

332 This dictionary can occupy significant memory, usually several times 

333 the size of the input array. 

334 

335 Notes 

336 ----- 

337 For a small array with few distinct values, one might use 

338 `numpy.unique()` to find all possible values, and ``(arr == val)`` to 

339 locate each value within that array. However, for large arrays, 

340 with many distinct values, this can become extremely inefficient, 

341 as locating each value would require a new search through the entire 

342 array. Using this function, there is essentially one search, with 

343 the indices saved for all distinct values. 

344 

345 This is useful when matching a categorical image (e.g. a segmentation 

346 or classification) to an associated image of other data, allowing 

347 any per-class statistic(s) to then be calculated. Provides a 

348 more flexible alternative to functions like ``scipy.ndimage.mean()`` 

349 and ``scipy.ndimage.variance()``. 

350 

351 Some other closely related functionality, with different strengths and 

352 weaknesses, can also be found in ``scipy.stats.binned_statistic()`` and 

353 the `scikit-image <https://scikit-image.org/>`_ function 

354 ``skimage.measure.regionprops()``. 

355 

356 Note for IDL users: this provides functionality equivalent to IDL's 

357 REVERSE_INDICES option (as per the IDL documentation for the 

358 `HISTOGRAM <https://www.l3harrisgeospatial.com/docs/histogram.html>`_ 

359 function). 

360 

361 .. versionadded:: 1.10.0 

362 

363 See Also 

364 -------- 

365 label, maximum, median, minimum_position, extrema, sum, mean, variance, 

366 standard_deviation, numpy.where, numpy.unique 

367 

368 Examples 

369 -------- 

370 >>> import numpy as np 

371 >>> from scipy import ndimage 

372 >>> a = np.zeros((6, 6), dtype=int) 

373 >>> a[2:4, 2:4] = 1 

374 >>> a[4, 4] = 1 

375 >>> a[:2, :3] = 2 

376 >>> a[0, 5] = 3 

377 >>> a 

378 array([[2, 2, 2, 0, 0, 3], 

379 [2, 2, 2, 0, 0, 0], 

380 [0, 0, 1, 1, 0, 0], 

381 [0, 0, 1, 1, 0, 0], 

382 [0, 0, 0, 0, 1, 0], 

383 [0, 0, 0, 0, 0, 0]]) 

384 >>> val_indices = ndimage.value_indices(a) 

385 

386 The dictionary `val_indices` will have an entry for each distinct 

387 value in the input array. 

388 

389 >>> val_indices.keys() 

390 dict_keys([0, 1, 2, 3]) 

391 

392 The entry for each value is an index tuple, locating the elements 

393 with that value. 

394 

395 >>> ndx1 = val_indices[1] 

396 >>> ndx1 

397 (array([2, 2, 3, 3, 4]), array([2, 3, 2, 3, 4])) 

398 

399 This can be used to index into the original array, or any other 

400 array with the same shape. 

401 

402 >>> a[ndx1] 

403 array([1, 1, 1, 1, 1]) 

404 

405 If the zeros were to be ignored, then the resulting dictionary 

406 would no longer have an entry for zero. 

407 

408 >>> val_indices = ndimage.value_indices(a, ignore_value=0) 

409 >>> val_indices.keys() 

410 dict_keys([1, 2, 3]) 

411 

412 """ 

413 # Cope with ignore_value being None, without too much extra complexity 

414 # in the C code. If not None, the value is passed in as a numpy array 

415 # with the same dtype as arr. 

416 ignore_value_arr = numpy.zeros((1,), dtype=arr.dtype) 

417 ignoreIsNone = (ignore_value is None) 

418 if not ignoreIsNone: 

419 ignore_value_arr[0] = ignore_value_arr.dtype.type(ignore_value) 

420 

421 val_indices = _nd_image.value_indices(arr, ignoreIsNone, ignore_value_arr) 

422 return val_indices 

423 

424 

425def labeled_comprehension(input, labels, index, func, out_dtype, default, pass_positions=False): 

426 """ 

427 Roughly equivalent to [func(input[labels == i]) for i in index]. 

428 

429 Sequentially applies an arbitrary function (that works on array_like input) 

430 to subsets of an N-D image array specified by `labels` and `index`. 

431 The option exists to provide the function with positional parameters as the 

432 second argument. 

433 

434 Parameters 

435 ---------- 

436 input : array_like 

437 Data from which to select `labels` to process. 

438 labels : array_like or None 

439 Labels to objects in `input`. 

440 If not None, array must be same shape as `input`. 

441 If None, `func` is applied to raveled `input`. 

442 index : int, sequence of ints or None 

443 Subset of `labels` to which to apply `func`. 

444 If a scalar, a single value is returned. 

445 If None, `func` is applied to all non-zero values of `labels`. 

446 func : callable 

447 Python function to apply to `labels` from `input`. 

448 out_dtype : dtype 

449 Dtype to use for `result`. 

450 default : int, float or None 

451 Default return value when a element of `index` does not exist 

452 in `labels`. 

453 pass_positions : bool, optional 

454 If True, pass linear indices to `func` as a second argument. 

455 Default is False. 

456 

457 Returns 

458 ------- 

459 result : ndarray 

460 Result of applying `func` to each of `labels` to `input` in `index`. 

461 

462 Examples 

463 -------- 

464 >>> import numpy as np 

465 >>> a = np.array([[1, 2, 0, 0], 

466 ... [5, 3, 0, 4], 

467 ... [0, 0, 0, 7], 

468 ... [9, 3, 0, 0]]) 

469 >>> from scipy import ndimage 

470 >>> lbl, nlbl = ndimage.label(a) 

471 >>> lbls = np.arange(1, nlbl+1) 

472 >>> ndimage.labeled_comprehension(a, lbl, lbls, np.mean, float, 0) 

473 array([ 2.75, 5.5 , 6. ]) 

474 

475 Falling back to `default`: 

476 

477 >>> lbls = np.arange(1, nlbl+2) 

478 >>> ndimage.labeled_comprehension(a, lbl, lbls, np.mean, float, -1) 

479 array([ 2.75, 5.5 , 6. , -1. ]) 

480 

481 Passing positions: 

482 

483 >>> def fn(val, pos): 

484 ... print("fn says: %s : %s" % (val, pos)) 

485 ... return (val.sum()) if (pos.sum() % 2 == 0) else (-val.sum()) 

486 ... 

487 >>> ndimage.labeled_comprehension(a, lbl, lbls, fn, float, 0, True) 

488 fn says: [1 2 5 3] : [0 1 4 5] 

489 fn says: [4 7] : [ 7 11] 

490 fn says: [9 3] : [12 13] 

491 array([ 11., 11., -12., 0.]) 

492 

493 """ 

494 

495 as_scalar = numpy.isscalar(index) 

496 input = numpy.asarray(input) 

497 

498 if pass_positions: 

499 positions = numpy.arange(input.size).reshape(input.shape) 

500 

501 if labels is None: 

502 if index is not None: 

503 raise ValueError("index without defined labels") 

504 if not pass_positions: 

505 return func(input.ravel()) 

506 else: 

507 return func(input.ravel(), positions.ravel()) 

508 

509 try: 

510 input, labels = numpy.broadcast_arrays(input, labels) 

511 except ValueError as e: 

512 raise ValueError("input and labels must have the same shape " 

513 "(excepting dimensions with width 1)") from e 

514 

515 if index is None: 

516 if not pass_positions: 

517 return func(input[labels > 0]) 

518 else: 

519 return func(input[labels > 0], positions[labels > 0]) 

520 

521 index = numpy.atleast_1d(index) 

522 if np.any(index.astype(labels.dtype).astype(index.dtype) != index): 

523 raise ValueError("Cannot convert index values from <%s> to <%s> " 

524 "(labels' type) without loss of precision" % 

525 (index.dtype, labels.dtype)) 

526 

527 index = index.astype(labels.dtype) 

528 

529 # optimization: find min/max in index, and select those parts of labels, input, and positions 

530 lo = index.min() 

531 hi = index.max() 

532 mask = (labels >= lo) & (labels <= hi) 

533 

534 # this also ravels the arrays 

535 labels = labels[mask] 

536 input = input[mask] 

537 if pass_positions: 

538 positions = positions[mask] 

539 

540 # sort everything by labels 

541 label_order = labels.argsort() 

542 labels = labels[label_order] 

543 input = input[label_order] 

544 if pass_positions: 

545 positions = positions[label_order] 

546 

547 index_order = index.argsort() 

548 sorted_index = index[index_order] 

549 

550 def do_map(inputs, output): 

551 """labels must be sorted""" 

552 nidx = sorted_index.size 

553 

554 # Find boundaries for each stretch of constant labels 

555 # This could be faster, but we already paid N log N to sort labels. 

556 lo = numpy.searchsorted(labels, sorted_index, side='left') 

557 hi = numpy.searchsorted(labels, sorted_index, side='right') 

558 

559 for i, l, h in zip(range(nidx), lo, hi): 

560 if l == h: 

561 continue 

562 output[i] = func(*[inp[l:h] for inp in inputs]) 

563 

564 temp = numpy.empty(index.shape, out_dtype) 

565 temp[:] = default 

566 if not pass_positions: 

567 do_map([input], temp) 

568 else: 

569 do_map([input, positions], temp) 

570 

571 output = numpy.zeros(index.shape, out_dtype) 

572 output[index_order] = temp 

573 if as_scalar: 

574 output = output[0] 

575 

576 return output 

577 

578 

579def _safely_castable_to_int(dt): 

580 """Test whether the NumPy data type `dt` can be safely cast to an int.""" 

581 int_size = np.dtype(int).itemsize 

582 safe = ((np.issubdtype(dt, np.signedinteger) and dt.itemsize <= int_size) or 

583 (np.issubdtype(dt, np.unsignedinteger) and dt.itemsize < int_size)) 

584 return safe 

585 

586 

587def _stats(input, labels=None, index=None, centered=False): 

588 """Count, sum, and optionally compute (sum - centre)^2 of input by label 

589 

590 Parameters 

591 ---------- 

592 input : array_like, N-D 

593 The input data to be analyzed. 

594 labels : array_like (N-D), optional 

595 The labels of the data in `input`. This array must be broadcast 

596 compatible with `input`; typically, it is the same shape as `input`. 

597 If `labels` is None, all nonzero values in `input` are treated as 

598 the single labeled group. 

599 index : label or sequence of labels, optional 

600 These are the labels of the groups for which the stats are computed. 

601 If `index` is None, the stats are computed for the single group where 

602 `labels` is greater than 0. 

603 centered : bool, optional 

604 If True, the centered sum of squares for each labeled group is 

605 also returned. Default is False. 

606 

607 Returns 

608 ------- 

609 counts : int or ndarray of ints 

610 The number of elements in each labeled group. 

611 sums : scalar or ndarray of scalars 

612 The sums of the values in each labeled group. 

613 sums_c : scalar or ndarray of scalars, optional 

614 The sums of mean-centered squares of the values in each labeled group. 

615 This is only returned if `centered` is True. 

616 

617 """ 

618 def single_group(vals): 

619 if centered: 

620 vals_c = vals - vals.mean() 

621 return vals.size, vals.sum(), (vals_c * vals_c.conjugate()).sum() 

622 else: 

623 return vals.size, vals.sum() 

624 

625 if labels is None: 

626 return single_group(input) 

627 

628 # ensure input and labels match sizes 

629 input, labels = numpy.broadcast_arrays(input, labels) 

630 

631 if index is None: 

632 return single_group(input[labels > 0]) 

633 

634 if numpy.isscalar(index): 

635 return single_group(input[labels == index]) 

636 

637 def _sum_centered(labels): 

638 # `labels` is expected to be an ndarray with the same shape as `input`. 

639 # It must contain the label indices (which are not necessarily the labels 

640 # themselves). 

641 means = sums / counts 

642 centered_input = input - means[labels] 

643 # bincount expects 1-D inputs, so we ravel the arguments. 

644 bc = numpy.bincount(labels.ravel(), 

645 weights=(centered_input * 

646 centered_input.conjugate()).ravel()) 

647 return bc 

648 

649 # Remap labels to unique integers if necessary, or if the largest 

650 # label is larger than the number of values. 

651 

652 if (not _safely_castable_to_int(labels.dtype) or 

653 labels.min() < 0 or labels.max() > labels.size): 

654 # Use numpy.unique to generate the label indices. `new_labels` will 

655 # be 1-D, but it should be interpreted as the flattened N-D array of 

656 # label indices. 

657 unique_labels, new_labels = numpy.unique(labels, return_inverse=True) 

658 counts = numpy.bincount(new_labels) 

659 sums = numpy.bincount(new_labels, weights=input.ravel()) 

660 if centered: 

661 # Compute the sum of the mean-centered squares. 

662 # We must reshape new_labels to the N-D shape of `input` before 

663 # passing it _sum_centered. 

664 sums_c = _sum_centered(new_labels.reshape(labels.shape)) 

665 idxs = numpy.searchsorted(unique_labels, index) 

666 # make all of idxs valid 

667 idxs[idxs >= unique_labels.size] = 0 

668 found = (unique_labels[idxs] == index) 

669 else: 

670 # labels are an integer type allowed by bincount, and there aren't too 

671 # many, so call bincount directly. 

672 counts = numpy.bincount(labels.ravel()) 

673 sums = numpy.bincount(labels.ravel(), weights=input.ravel()) 

674 if centered: 

675 sums_c = _sum_centered(labels) 

676 # make sure all index values are valid 

677 idxs = numpy.asanyarray(index, numpy.int_).copy() 

678 found = (idxs >= 0) & (idxs < counts.size) 

679 idxs[~found] = 0 

680 

681 counts = counts[idxs] 

682 counts[~found] = 0 

683 sums = sums[idxs] 

684 sums[~found] = 0 

685 

686 if not centered: 

687 return (counts, sums) 

688 else: 

689 sums_c = sums_c[idxs] 

690 sums_c[~found] = 0 

691 return (counts, sums, sums_c) 

692 

693 

694def sum(input, labels=None, index=None): 

695 """ 

696 Calculate the sum of the values of the array. 

697 

698 Notes 

699 ----- 

700 This is an alias for `ndimage.sum_labels` kept for backwards compatibility 

701 reasons, for new code please prefer `sum_labels`. See the `sum_labels` 

702 docstring for more details. 

703 

704 """ 

705 return sum_labels(input, labels, index) 

706 

707 

708def sum_labels(input, labels=None, index=None): 

709 """ 

710 Calculate the sum of the values of the array. 

711 

712 Parameters 

713 ---------- 

714 input : array_like 

715 Values of `input` inside the regions defined by `labels` 

716 are summed together. 

717 labels : array_like of ints, optional 

718 Assign labels to the values of the array. Has to have the same shape as 

719 `input`. 

720 index : array_like, optional 

721 A single label number or a sequence of label numbers of 

722 the objects to be measured. 

723 

724 Returns 

725 ------- 

726 sum : ndarray or scalar 

727 An array of the sums of values of `input` inside the regions defined 

728 by `labels` with the same shape as `index`. If 'index' is None or scalar, 

729 a scalar is returned. 

730 

731 See Also 

732 -------- 

733 mean, median 

734 

735 Examples 

736 -------- 

737 >>> from scipy import ndimage 

738 >>> input = [0,1,2,3] 

739 >>> labels = [1,1,2,2] 

740 >>> ndimage.sum_labels(input, labels, index=[1,2]) 

741 [1.0, 5.0] 

742 >>> ndimage.sum_labels(input, labels, index=1) 

743 1 

744 >>> ndimage.sum_labels(input, labels) 

745 6 

746 

747 

748 """ 

749 count, sum = _stats(input, labels, index) 

750 return sum 

751 

752 

753def mean(input, labels=None, index=None): 

754 """ 

755 Calculate the mean of the values of an array at labels. 

756 

757 Parameters 

758 ---------- 

759 input : array_like 

760 Array on which to compute the mean of elements over distinct 

761 regions. 

762 labels : array_like, optional 

763 Array of labels of same shape, or broadcastable to the same shape as 

764 `input`. All elements sharing the same label form one region over 

765 which the mean of the elements is computed. 

766 index : int or sequence of ints, optional 

767 Labels of the objects over which the mean is to be computed. 

768 Default is None, in which case the mean for all values where label is 

769 greater than 0 is calculated. 

770 

771 Returns 

772 ------- 

773 out : list 

774 Sequence of same length as `index`, with the mean of the different 

775 regions labeled by the labels in `index`. 

776 

777 See Also 

778 -------- 

779 variance, standard_deviation, minimum, maximum, sum, label 

780 

781 Examples 

782 -------- 

783 >>> from scipy import ndimage 

784 >>> import numpy as np 

785 >>> a = np.arange(25).reshape((5,5)) 

786 >>> labels = np.zeros_like(a) 

787 >>> labels[3:5,3:5] = 1 

788 >>> index = np.unique(labels) 

789 >>> labels 

790 array([[0, 0, 0, 0, 0], 

791 [0, 0, 0, 0, 0], 

792 [0, 0, 0, 0, 0], 

793 [0, 0, 0, 1, 1], 

794 [0, 0, 0, 1, 1]]) 

795 >>> index 

796 array([0, 1]) 

797 >>> ndimage.mean(a, labels=labels, index=index) 

798 [10.285714285714286, 21.0] 

799 

800 """ 

801 

802 count, sum = _stats(input, labels, index) 

803 return sum / numpy.asanyarray(count).astype(numpy.float64) 

804 

805 

806def variance(input, labels=None, index=None): 

807 """ 

808 Calculate the variance of the values of an N-D image array, optionally at 

809 specified sub-regions. 

810 

811 Parameters 

812 ---------- 

813 input : array_like 

814 Nd-image data to process. 

815 labels : array_like, optional 

816 Labels defining sub-regions in `input`. 

817 If not None, must be same shape as `input`. 

818 index : int or sequence of ints, optional 

819 `labels` to include in output. If None (default), all values where 

820 `labels` is non-zero are used. 

821 

822 Returns 

823 ------- 

824 variance : float or ndarray 

825 Values of variance, for each sub-region if `labels` and `index` are 

826 specified. 

827 

828 See Also 

829 -------- 

830 label, standard_deviation, maximum, minimum, extrema 

831 

832 Examples 

833 -------- 

834 >>> import numpy as np 

835 >>> a = np.array([[1, 2, 0, 0], 

836 ... [5, 3, 0, 4], 

837 ... [0, 0, 0, 7], 

838 ... [9, 3, 0, 0]]) 

839 >>> from scipy import ndimage 

840 >>> ndimage.variance(a) 

841 7.609375 

842 

843 Features to process can be specified using `labels` and `index`: 

844 

845 >>> lbl, nlbl = ndimage.label(a) 

846 >>> ndimage.variance(a, lbl, index=np.arange(1, nlbl+1)) 

847 array([ 2.1875, 2.25 , 9. ]) 

848 

849 If no index is given, all non-zero `labels` are processed: 

850 

851 >>> ndimage.variance(a, lbl) 

852 6.1875 

853 

854 """ 

855 count, sum, sum_c_sq = _stats(input, labels, index, centered=True) 

856 return sum_c_sq / np.asanyarray(count).astype(float) 

857 

858 

859def standard_deviation(input, labels=None, index=None): 

860 """ 

861 Calculate the standard deviation of the values of an N-D image array, 

862 optionally at specified sub-regions. 

863 

864 Parameters 

865 ---------- 

866 input : array_like 

867 N-D image data to process. 

868 labels : array_like, optional 

869 Labels to identify sub-regions in `input`. 

870 If not None, must be same shape as `input`. 

871 index : int or sequence of ints, optional 

872 `labels` to include in output. If None (default), all values where 

873 `labels` is non-zero are used. 

874 

875 Returns 

876 ------- 

877 standard_deviation : float or ndarray 

878 Values of standard deviation, for each sub-region if `labels` and 

879 `index` are specified. 

880 

881 See Also 

882 -------- 

883 label, variance, maximum, minimum, extrema 

884 

885 Examples 

886 -------- 

887 >>> import numpy as np 

888 >>> a = np.array([[1, 2, 0, 0], 

889 ... [5, 3, 0, 4], 

890 ... [0, 0, 0, 7], 

891 ... [9, 3, 0, 0]]) 

892 >>> from scipy import ndimage 

893 >>> ndimage.standard_deviation(a) 

894 2.7585095613392387 

895 

896 Features to process can be specified using `labels` and `index`: 

897 

898 >>> lbl, nlbl = ndimage.label(a) 

899 >>> ndimage.standard_deviation(a, lbl, index=np.arange(1, nlbl+1)) 

900 array([ 1.479, 1.5 , 3. ]) 

901 

902 If no index is given, non-zero `labels` are processed: 

903 

904 >>> ndimage.standard_deviation(a, lbl) 

905 2.4874685927665499 

906 

907 """ 

908 return numpy.sqrt(variance(input, labels, index)) 

909 

910 

911def _select(input, labels=None, index=None, find_min=False, find_max=False, 

912 find_min_positions=False, find_max_positions=False, 

913 find_median=False): 

914 """Returns min, max, or both, plus their positions (if requested), and 

915 median.""" 

916 

917 input = numpy.asanyarray(input) 

918 

919 find_positions = find_min_positions or find_max_positions 

920 positions = None 

921 if find_positions: 

922 positions = numpy.arange(input.size).reshape(input.shape) 

923 

924 def single_group(vals, positions): 

925 result = [] 

926 if find_min: 

927 result += [vals.min()] 

928 if find_min_positions: 

929 result += [positions[vals == vals.min()][0]] 

930 if find_max: 

931 result += [vals.max()] 

932 if find_max_positions: 

933 result += [positions[vals == vals.max()][0]] 

934 if find_median: 

935 result += [numpy.median(vals)] 

936 return result 

937 

938 if labels is None: 

939 return single_group(input, positions) 

940 

941 # ensure input and labels match sizes 

942 input, labels = numpy.broadcast_arrays(input, labels) 

943 

944 if index is None: 

945 mask = (labels > 0) 

946 masked_positions = None 

947 if find_positions: 

948 masked_positions = positions[mask] 

949 return single_group(input[mask], masked_positions) 

950 

951 if numpy.isscalar(index): 

952 mask = (labels == index) 

953 masked_positions = None 

954 if find_positions: 

955 masked_positions = positions[mask] 

956 return single_group(input[mask], masked_positions) 

957 

958 # remap labels to unique integers if necessary, or if the largest 

959 # label is larger than the number of values. 

960 if (not _safely_castable_to_int(labels.dtype) or 

961 labels.min() < 0 or labels.max() > labels.size): 

962 # remap labels, and indexes 

963 unique_labels, labels = numpy.unique(labels, return_inverse=True) 

964 idxs = numpy.searchsorted(unique_labels, index) 

965 

966 # make all of idxs valid 

967 idxs[idxs >= unique_labels.size] = 0 

968 found = (unique_labels[idxs] == index) 

969 else: 

970 # labels are an integer type, and there aren't too many 

971 idxs = numpy.asanyarray(index, numpy.int_).copy() 

972 found = (idxs >= 0) & (idxs <= labels.max()) 

973 

974 idxs[~ found] = labels.max() + 1 

975 

976 if find_median: 

977 order = numpy.lexsort((input.ravel(), labels.ravel())) 

978 else: 

979 order = input.ravel().argsort() 

980 input = input.ravel()[order] 

981 labels = labels.ravel()[order] 

982 if find_positions: 

983 positions = positions.ravel()[order] 

984 

985 result = [] 

986 if find_min: 

987 mins = numpy.zeros(labels.max() + 2, input.dtype) 

988 mins[labels[::-1]] = input[::-1] 

989 result += [mins[idxs]] 

990 if find_min_positions: 

991 minpos = numpy.zeros(labels.max() + 2, int) 

992 minpos[labels[::-1]] = positions[::-1] 

993 result += [minpos[idxs]] 

994 if find_max: 

995 maxs = numpy.zeros(labels.max() + 2, input.dtype) 

996 maxs[labels] = input 

997 result += [maxs[idxs]] 

998 if find_max_positions: 

999 maxpos = numpy.zeros(labels.max() + 2, int) 

1000 maxpos[labels] = positions 

1001 result += [maxpos[idxs]] 

1002 if find_median: 

1003 locs = numpy.arange(len(labels)) 

1004 lo = numpy.zeros(labels.max() + 2, numpy.int_) 

1005 lo[labels[::-1]] = locs[::-1] 

1006 hi = numpy.zeros(labels.max() + 2, numpy.int_) 

1007 hi[labels] = locs 

1008 lo = lo[idxs] 

1009 hi = hi[idxs] 

1010 # lo is an index to the lowest value in input for each label, 

1011 # hi is an index to the largest value. 

1012 # move them to be either the same ((hi - lo) % 2 == 0) or next 

1013 # to each other ((hi - lo) % 2 == 1), then average. 

1014 step = (hi - lo) // 2 

1015 lo += step 

1016 hi -= step 

1017 if (np.issubdtype(input.dtype, np.integer)