Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-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. 

260 

261 See Also 

262 -------- 

263 label, center_of_mass 

264 

265 Notes 

266 ----- 

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

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

269 

270 Examples 

271 -------- 

272 >>> from scipy import ndimage 

273 >>> import numpy as np 

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

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

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

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

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

279 >>> a 

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

281 [2, 2, 2, 0, 0, 0], 

282 [0, 0, 1, 1, 0, 0], 

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

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

285 [0, 0, 0, 0, 0, 0]]) 

286 >>> ndimage.find_objects(a) 

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

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

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

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

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

292 

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

294 >>> a[loc] 

295 array([[1, 1, 0], 

296 [1, 1, 0], 

297 [0, 0, 1]]) 

298 

299 """ 

300 input = numpy.asarray(input) 

301 if numpy.iscomplexobj(input): 

302 raise TypeError('Complex type not supported') 

303 

304 if max_label < 1: 

305 max_label = input.max() 

306 

307 return _nd_image.find_objects(input, max_label) 

308 

309 

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

311 """ 

312 Find indices of each distinct value in given array. 

313 

314 Parameters 

315 ---------- 

316 arr : ndarray of ints 

317 Array containing integer values. 

318 ignore_value : int, optional 

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

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

321 is None. 

322 

323 Returns 

324 ------- 

325 indices : dictionary 

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

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

328 index tuples covering all occurrences of the value within the 

329 array. 

330 

331 This dictionary can occupy significant memory, usually several times 

332 the size of the input array. 

333 

334 Notes 

335 ----- 

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

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

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

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

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

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

342 the indices saved for all distinct values. 

343 

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

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

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

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

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

349 

350 Some other closely related functionality, with different strengths and 

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

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

353 ``skimage.measure.regionprops()``. 

354 

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

356 REVERSE_INDICES option (as per the IDL documentation for the 

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

358 function). 

359 

360 .. versionadded:: 1.10.0 

361 

362 See Also 

363 -------- 

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

365 standard_deviation, numpy.where, numpy.unique 

366 

367 Examples 

368 -------- 

369 >>> import numpy as np 

370 >>> from scipy import ndimage 

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

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

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

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

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

376 >>> a 

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

378 [2, 2, 2, 0, 0, 0], 

379 [0, 0, 1, 1, 0, 0], 

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

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

382 [0, 0, 0, 0, 0, 0]]) 

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

384 

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

386 value in the input array. 

387 

388 >>> val_indices.keys() 

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

390 

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

392 with that value. 

393 

394 >>> ndx1 = val_indices[1] 

395 >>> ndx1 

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

397 

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

399 array with the same shape. 

400 

401 >>> a[ndx1] 

402 array([1, 1, 1, 1, 1]) 

403 

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

405 would no longer have an entry for zero. 

406 

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

408 >>> val_indices.keys() 

409 dict_keys([1, 2, 3]) 

410 

411 """ 

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

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

414 # with the same dtype as arr. 

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

416 ignoreIsNone = (ignore_value is None) 

417 if not ignoreIsNone: 

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

419 

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

421 return val_indices 

422 

423 

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

425 """ 

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

427 

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

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

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

431 second argument. 

432 

433 Parameters 

434 ---------- 

435 input : array_like 

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

437 labels : array_like or None 

438 Labels to objects in `input`. 

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

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

441 index : int, sequence of ints or None 

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

443 If a scalar, a single value is returned. 

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

445 func : callable 

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

447 out_dtype : dtype 

448 Dtype to use for `result`. 

449 default : int, float or None 

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

451 in `labels`. 

452 pass_positions : bool, optional 

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

454 Default is False. 

455 

456 Returns 

457 ------- 

458 result : ndarray 

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

460 

461 Examples 

462 -------- 

463 >>> import numpy as np 

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

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

466 ... [0, 0, 0, 7], 

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

468 >>> from scipy import ndimage 

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

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

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

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

473 

474 Falling back to `default`: 

475 

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

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

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

479 

480 Passing positions: 

481 

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

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

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

485 ... 

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

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

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

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

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

491 

492 """ 

493 

494 as_scalar = numpy.isscalar(index) 

495 input = numpy.asarray(input) 

496 

497 if pass_positions: 

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

499 

500 if labels is None: 

501 if index is not None: 

502 raise ValueError("index without defined labels") 

503 if not pass_positions: 

504 return func(input.ravel()) 

505 else: 

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

507 

508 try: 

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

510 except ValueError as e: 

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

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

513 

514 if index is None: 

515 if not pass_positions: 

516 return func(input[labels > 0]) 

517 else: 

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

519 

520 index = numpy.atleast_1d(index) 

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

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

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

524 (index.dtype, labels.dtype)) 

525 

526 index = index.astype(labels.dtype) 

527 

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

529 lo = index.min() 

530 hi = index.max() 

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

532 

533 # this also ravels the arrays 

534 labels = labels[mask] 

535 input = input[mask] 

536 if pass_positions: 

537 positions = positions[mask] 

538 

539 # sort everything by labels 

540 label_order = labels.argsort() 

541 labels = labels[label_order] 

542 input = input[label_order] 

543 if pass_positions: 

544 positions = positions[label_order] 

545 

546 index_order = index.argsort() 

547 sorted_index = index[index_order] 

548 

549 def do_map(inputs, output): 

550 """labels must be sorted""" 

551 nidx = sorted_index.size 

552 

553 # Find boundaries for each stretch of constant labels 

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

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

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

557 

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

559 if l == h: 

560 continue 

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

562 

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

564 temp[:] = default 

565 if not pass_positions: 

566 do_map([input], temp) 

567 else: 

568 do_map([input, positions], temp) 

569 

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

571 output[index_order] = temp 

572 if as_scalar: 

573 output = output[0] 

574 

575 return output 

576 

577 

578def _safely_castable_to_int(dt): 

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

580 int_size = np.dtype(int).itemsize 

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

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

583 return safe 

584 

585 

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

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

588 

589 Parameters 

590 ---------- 

591 input : array_like, N-D 

592 The input data to be analyzed. 

593 labels : array_like (N-D), optional 

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

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

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

597 the single labeled group. 

598 index : label or sequence of labels, optional 

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

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

601 `labels` is greater than 0. 

602 centered : bool, optional 

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

604 also returned. Default is False. 

605 

606 Returns 

607 ------- 

608 counts : int or ndarray of ints 

609 The number of elements in each labeled group. 

610 sums : scalar or ndarray of scalars 

611 The sums of the values in each labeled group. 

612 sums_c : scalar or ndarray of scalars, optional 

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

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

615 

616 """ 

617 def single_group(vals): 

618 if centered: 

619 vals_c = vals - vals.mean() 

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

621 else: 

622 return vals.size, vals.sum() 

623 

624 if labels is None: 

625 return single_group(input) 

626 

627 # ensure input and labels match sizes 

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

629 

630 if index is None: 

631 return single_group(input[labels > 0]) 

632 

633 if numpy.isscalar(index): 

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

635 

636 def _sum_centered(labels): 

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

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

639 # themselves). 

640 means = sums / counts 

641 centered_input = input - means[labels] 

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

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

644 weights=(centered_input * 

645 centered_input.conjugate()).ravel()) 

646 return bc 

647 

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

649 # label is larger than the number of values. 

650 

651 if (not _safely_castable_to_int(labels.dtype) or 

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

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

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

655 # label indices. 

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

657 counts = numpy.bincount(new_labels) 

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

659 if centered: 

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

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

662 # passing it _sum_centered. 

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

664 idxs = numpy.searchsorted(unique_labels, index) 

665 # make all of idxs valid 

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

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

668 else: 

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

670 # many, so call bincount directly. 

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

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

673 if centered: 

674 sums_c = _sum_centered(labels) 

675 # make sure all index values are valid 

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

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

678 idxs[~found] = 0 

679 

680 counts = counts[idxs] 

681 counts[~found] = 0 

682 sums = sums[idxs] 

683 sums[~found] = 0 

684 

685 if not centered: 

686 return (counts, sums) 

687 else: 

688 sums_c = sums_c[idxs] 

689 sums_c[~found] = 0 

690 return (counts, sums, sums_c) 

691 

692 

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

694 """ 

695 Calculate the sum of the values of the array. 

696 

697 Notes 

698 ----- 

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

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

701 docstring for more details. 

702 

703 """ 

704 return sum_labels(input, labels, index) 

705 

706 

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

708 """ 

709 Calculate the sum of the values of the array. 

710 

711 Parameters 

712 ---------- 

713 input : array_like 

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

715 are summed together. 

716 labels : array_like of ints, optional 

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

718 `input`. 

719 index : array_like, optional 

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

721 the objects to be measured. 

722 

723 Returns 

724 ------- 

725 sum : ndarray or scalar 

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

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

728 a scalar is returned. 

729 

730 See Also 

731 -------- 

732 mean, median 

733 

734 Examples 

735 -------- 

736 >>> from scipy import ndimage 

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

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

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

740 [1.0, 5.0] 

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

742 1 

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

744 6 

745 

746 

747 """ 

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

749 return sum 

750 

751 

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

753 """ 

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

755 

756 Parameters 

757 ---------- 

758 input : array_like 

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

760 regions. 

761 labels : array_like, optional 

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

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

764 which the mean of the elements is computed. 

765 index : int or sequence of ints, optional 

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

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

768 greater than 0 is calculated. 

769 

770 Returns 

771 ------- 

772 out : list 

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

774 regions labeled by the labels in `index`. 

775 

776 See Also 

777 -------- 

778 variance, standard_deviation, minimum, maximum, sum, label 

779 

780 Examples 

781 -------- 

782 >>> from scipy import ndimage 

783 >>> import numpy as np 

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

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

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

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

788 >>> labels 

789 array([[0, 0, 0, 0, 0], 

790 [0, 0, 0, 0, 0], 

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

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

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

794 >>> index 

795 array([0, 1]) 

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

797 [10.285714285714286, 21.0] 

798 

799 """ 

800 

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

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

803 

804 

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

806 """ 

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

808 specified sub-regions. 

809 

810 Parameters 

811 ---------- 

812 input : array_like 

813 Nd-image data to process. 

814 labels : array_like, optional 

815 Labels defining sub-regions in `input`. 

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

817 index : int or sequence of ints, optional 

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

819 `labels` is non-zero are used. 

820 

821 Returns 

822 ------- 

823 variance : float or ndarray 

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

825 specified. 

826 

827 See Also 

828 -------- 

829 label, standard_deviation, maximum, minimum, extrema 

830 

831 Examples 

832 -------- 

833 >>> import numpy as np 

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

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

836 ... [0, 0, 0, 7], 

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

838 >>> from scipy import ndimage 

839 >>> ndimage.variance(a) 

840 7.609375 

841 

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

843 

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

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

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

847 

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

849 

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

851 6.1875 

852 

853 """ 

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

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

856 

857 

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

859 """ 

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

861 optionally at specified sub-regions. 

862 

863 Parameters 

864 ---------- 

865 input : array_like 

866 N-D image data to process. 

867 labels : array_like, optional 

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

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

870 index : int or sequence of ints, optional 

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

872 `labels` is non-zero are used. 

873 

874 Returns 

875 ------- 

876 standard_deviation : float or ndarray 

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

878 `index` are specified. 

879 

880 See Also 

881 -------- 

882 label, variance, maximum, minimum, extrema 

883 

884 Examples 

885 -------- 

886 >>> import numpy as np 

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

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

889 ... [0, 0, 0, 7], 

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

891 >>> from scipy import ndimage 

892 >>> ndimage.standard_deviation(a) 

893 2.7585095613392387 

894 

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

896 

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

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

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

900 

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

902 

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

904 2.4874685927665499 

905 

906 """ 

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

908 

909 

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

911 find_min_positions=False, find_max_positions=False, 

912 find_median=False): 

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

914 median.""" 

915 

916 input = numpy.asanyarray(input) 

917 

918 find_positions = find_min_positions or find_max_positions 

919 positions = None 

920 if find_positions: 

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

922 

923 def single_group(vals, positions): 

924 result = [] 

925 if find_min: 

926 result += [vals.min()] 

927 if find_min_positions: 

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

929 if find_max: 

930 result += [vals.max()] 

931 if find_max_positions: 

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

933 if find_median: 

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

935 return result 

936 

937 if labels is None: 

938 return single_group(input, positions) 

939 

940 # ensure input and labels match sizes 

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

942 

943 if index is None: 

944 mask = (labels > 0) 

945 masked_positions = None 

946 if find_positions: 

947 masked_positions = positions[mask] 

948 return single_group(input[mask], masked_positions) 

949 

950 if numpy.isscalar(index): 

951 mask = (labels == index) 

952 masked_positions = None 

953 if find_positions: 

954 masked_positions = positions[mask] 

955 return single_group(input[mask], masked_positions) 

956 

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

958 # label is larger than the number of values. 

959 if (not _safely_castable_to_int(labels.dtype) or 

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

961 # remap labels, and indexes 

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

963 idxs = numpy.searchsorted(unique_labels, index) 

964 

965 # make all of idxs valid 

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

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

968 else: 

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

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

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

972 

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

974 

975 if find_median: 

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

977 else: 

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

979 input = input.ravel()[order] 

980 labels = labels.ravel()[order] 

981 if find_positions: 

982 positions = positions.ravel()[order] 

983 

984 result = [] 

985 if find_min: 

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

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

988 result += [mins[idxs]] 

989 if find_min_positions: 

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

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

992 result += [minpos[idxs]] 

993 if find_max: 

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

995 maxs[labels] = input 

996 result += [maxs[idxs]] 

997 if find_max_positions: 

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

999 maxpos[labels] = positions 

1000 result += [maxpos[idxs]] 

1001 if find_median: 

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

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

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

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

1006 hi[labels] = locs 

1007 lo = lo[idxs] 

1008 hi = hi[idxs] 

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

1010 # hi is an index to the largest value. 

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

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

1013 step = (hi - lo) // 2 

1014 lo += step 

1015 hi -= step 

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