Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

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

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: 

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("insufficient bit-depth in requested output type") 

229 

230 if caller_provided_output: 

231 # result was written in-place 

232 return max_label 

233 else: 

234 return output, max_label 

235 

236 

237def find_objects(input, max_label=0): 

238 """ 

239 Find objects in a labeled array. 

240 

241 Parameters 

242 ---------- 

243 input : ndarray of ints 

244 Array containing objects defined by different labels. Labels with 

245 value 0 are ignored. 

246 max_label : int, optional 

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

248 given, the positions of all objects are returned. 

249 

250 Returns 

251 ------- 

252 object_slices : list of tuples 

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

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

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

256 None is returned instead of a slice. 

257 

258 See Also 

259 -------- 

260 label, center_of_mass 

261 

262 Notes 

263 ----- 

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

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

266 

267 Examples 

268 -------- 

269 >>> from scipy import ndimage 

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

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

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

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

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

275 >>> a 

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

277 [2, 2, 2, 0, 0, 0], 

278 [0, 0, 1, 1, 0, 0], 

279 [0, 0, 1, 1, 0, 0], 

280 [0, 0, 0, 0, 1, 0], 

281 [0, 0, 0, 0, 0, 0]]) 

282 >>> ndimage.find_objects(a) 

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

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

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

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

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

288 

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

290 >>> a[loc] 

291 array([[1, 1, 0], 

292 [1, 1, 0], 

293 [0, 0, 1]]) 

294 

295 """ 

296 input = numpy.asarray(input) 

297 if numpy.iscomplexobj(input): 

298 raise TypeError('Complex type not supported') 

299 

300 if max_label < 1: 

301 max_label = input.max() 

302 

303 return _nd_image.find_objects(input, max_label) 

304 

305 

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

307 """ 

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

309 

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

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

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

313 second argument. 

314 

315 Parameters 

316 ---------- 

317 input : array_like 

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

319 labels : array_like or None 

320 Labels to objects in `input`. 

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

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

323 index : int, sequence of ints or None 

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

325 If a scalar, a single value is returned. 

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

327 func : callable 

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

329 out_dtype : dtype 

330 Dtype to use for `result`. 

331 default : int, float or None 

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

333 in `labels`. 

334 pass_positions : bool, optional 

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

336 Default is False. 

337 

338 Returns 

339 ------- 

340 result : ndarray 

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

342 

343 Examples 

344 -------- 

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

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

347 ... [0, 0, 0, 7], 

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

349 >>> from scipy import ndimage 

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

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

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

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

354 

355 Falling back to `default`: 

356 

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

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

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

360 

361 Passing positions: 

362 

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

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

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

366 ... 

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

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

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

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

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

372 

373 """ 

374 

375 as_scalar = numpy.isscalar(index) 

376 input = numpy.asarray(input) 

377 

378 if pass_positions: 

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

380 

381 if labels is None: 

382 if index is not None: 

383 raise ValueError("index without defined labels") 

384 if not pass_positions: 

385 return func(input.ravel()) 

386 else: 

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

388 

389 try: 

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

391 except ValueError: 

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

393 "(excepting dimensions with width 1)") 

394 

395 if index is None: 

396 if not pass_positions: 

397 return func(input[labels > 0]) 

398 else: 

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

400 

401 index = numpy.atleast_1d(index) 

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

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

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

405 (index.dtype, labels.dtype)) 

406 

407 index = index.astype(labels.dtype) 

408 

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

410 lo = index.min() 

411 hi = index.max() 

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

413 

414 # this also ravels the arrays 

415 labels = labels[mask] 

416 input = input[mask] 

417 if pass_positions: 

418 positions = positions[mask] 

419 

420 # sort everything by labels 

421 label_order = labels.argsort() 

422 labels = labels[label_order] 

423 input = input[label_order] 

424 if pass_positions: 

425 positions = positions[label_order] 

426 

427 index_order = index.argsort() 

428 sorted_index = index[index_order] 

429 

430 def do_map(inputs, output): 

431 """labels must be sorted""" 

432 nidx = sorted_index.size 

433 

434 # Find boundaries for each stretch of constant labels 

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

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

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

438 

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

440 if l == h: 

441 continue 

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

443 

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

445 temp[:] = default 

446 if not pass_positions: 

447 do_map([input], temp) 

448 else: 

449 do_map([input, positions], temp) 

450 

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

452 output[index_order] = temp 

453 if as_scalar: 

454 output = output[0] 

455 

456 return output 

457 

458 

459def _safely_castable_to_int(dt): 

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

461 int_size = np.dtype(int).itemsize 

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

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

464 return safe 

465 

466 

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

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

469 

470 Parameters 

471 ---------- 

472 input : array_like, N-D 

473 The input data to be analyzed. 

474 labels : array_like (N-D), optional 

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

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

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

478 the single labeled group. 

479 index : label or sequence of labels, optional 

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

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

482 `labels` is greater than 0. 

483 centered : bool, optional 

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

485 also returned. Default is False. 

486 

487 Returns 

488 ------- 

489 counts : int or ndarray of ints 

490 The number of elements in each labeled group. 

491 sums : scalar or ndarray of scalars 

492 The sums of the values in each labeled group. 

493 sums_c : scalar or ndarray of scalars, optional 

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

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

496 

497 """ 

498 def single_group(vals): 

499 if centered: 

500 vals_c = vals - vals.mean() 

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

502 else: 

503 return vals.size, vals.sum() 

504 

505 if labels is None: 

506 return single_group(input) 

507 

508 # ensure input and labels match sizes 

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

510 

511 if index is None: 

512 return single_group(input[labels > 0]) 

513 

514 if numpy.isscalar(index): 

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

516 

517 def _sum_centered(labels): 

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

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

520 # themselves). 

521 means = sums / counts 

522 centered_input = input - means[labels] 

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

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

525 weights=(centered_input * 

526 centered_input.conjugate()).ravel()) 

527 return bc 

528 

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

530 # label is larger than the number of values. 

531 

532 if (not _safely_castable_to_int(labels.dtype) or 

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

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

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

536 # label indices. 

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

538 counts = numpy.bincount(new_labels) 

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

540 if centered: 

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

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

543 # passing it _sum_centered. 

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

545 idxs = numpy.searchsorted(unique_labels, index) 

546 # make all of idxs valid 

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

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

549 else: 

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

551 # many, so call bincount directly. 

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

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

554 if centered: 

555 sums_c = _sum_centered(labels) 

556 # make sure all index values are valid 

557 idxs = numpy.asanyarray(index, numpy.int).copy() 

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

559 idxs[~found] = 0 

560 

561 counts = counts[idxs] 

562 counts[~found] = 0 

563 sums = sums[idxs] 

564 sums[~found] = 0 

565 

566 if not centered: 

567 return (counts, sums) 

568 else: 

569 sums_c = sums_c[idxs] 

570 sums_c[~found] = 0 

571 return (counts, sums, sums_c) 

572 

573 

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

575 """ 

576 Calculate the sum of the values of the array. 

577 

578 Parameters 

579 ---------- 

580 input : array_like 

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

582 are summed together. 

583 labels : array_like of ints, optional 

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

585 `input`. 

586 index : array_like, optional 

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

588 the objects to be measured. 

589 

590 Returns 

591 ------- 

592 sum : ndarray or scalar 

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

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

595 a scalar is returned. 

596 

597 See also 

598 -------- 

599 mean, median 

600 

601 Examples 

602 -------- 

603 >>> from scipy import ndimage 

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

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

606 >>> ndimage.sum(input, labels, index=[1,2]) 

607 [1.0, 5.0] 

608 >>> ndimage.sum(input, labels, index=1) 

609 1 

610 >>> ndimage.sum(input, labels) 

611 6 

612 

613 

614 """ 

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

616 return sum 

617 

618 

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

620 """ 

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

622 

623 Parameters 

624 ---------- 

625 input : array_like 

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

627 regions. 

628 labels : array_like, optional 

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

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

631 which the mean of the elements is computed. 

632 index : int or sequence of ints, optional 

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

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

635 greater than 0 is calculated. 

636 

637 Returns 

638 ------- 

639 out : list 

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

641 regions labeled by the labels in `index`. 

642 

643 See also 

644 -------- 

645 variance, standard_deviation, minimum, maximum, sum, label 

646 

647 Examples 

648 -------- 

649 >>> from scipy import ndimage 

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

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

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

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

654 >>> labels 

655 array([[0, 0, 0, 0, 0], 

656 [0, 0, 0, 0, 0], 

657 [0, 0, 0, 0, 0], 

658 [0, 0, 0, 1, 1], 

659 [0, 0, 0, 1, 1]]) 

660 >>> index 

661 array([0, 1]) 

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

663 [10.285714285714286, 21.0] 

664 

665 """ 

666 

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

668 return sum / numpy.asanyarray(count).astype(numpy.float) 

669 

670 

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

672 """ 

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

674 specified sub-regions. 

675 

676 Parameters 

677 ---------- 

678 input : array_like 

679 Nd-image data to process. 

680 labels : array_like, optional 

681 Labels defining sub-regions in `input`. 

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

683 index : int or sequence of ints, optional 

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

685 `labels` is non-zero are used. 

686 

687 Returns 

688 ------- 

689 variance : float or ndarray 

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

691 specified. 

692 

693 See Also 

694 -------- 

695 label, standard_deviation, maximum, minimum, extrema 

696 

697 Examples 

698 -------- 

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

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

701 ... [0, 0, 0, 7], 

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

703 >>> from scipy import ndimage 

704 >>> ndimage.variance(a) 

705 7.609375 

706 

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

708 

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

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

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

712 

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

714 

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

716 6.1875 

717 

718 """ 

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

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

721 

722 

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

724 """ 

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

726 optionally at specified sub-regions. 

727 

728 Parameters 

729 ---------- 

730 input : array_like 

731 N-D image data to process. 

732 labels : array_like, optional 

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

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

735 index : int or sequence of ints, optional 

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

737 `labels` is non-zero are used. 

738 

739 Returns 

740 ------- 

741 standard_deviation : float or ndarray 

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

743 `index` are specified. 

744 

745 See Also 

746 -------- 

747 label, variance, maximum, minimum, extrema 

748 

749 Examples 

750 -------- 

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

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

753 ... [0, 0, 0, 7], 

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

755 >>> from scipy import ndimage 

756 >>> ndimage.standard_deviation(a) 

757 2.7585095613392387 

758 

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

760 

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

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

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

764 

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

766 

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

768 2.4874685927665499 

769 

770 """ 

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

772 

773 

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

775 find_min_positions=False, find_max_positions=False, 

776 find_median=False): 

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

778 median.""" 

779 

780 input = numpy.asanyarray(input) 

781 

782 find_positions = find_min_positions or find_max_positions 

783 positions = None 

784 if find_positions: 

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

786 

787 def single_group(vals, positions): 

788 result = [] 

789 if find_min: 

790 result += [vals.min()] 

791 if find_min_positions: 

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

793 if find_max: 

794 result += [vals.max()] 

795 if find_max_positions: 

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

797 if find_median: 

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

799 return result 

800 

801 if labels is None: 

802 return single_group(input, positions) 

803 

804 # ensure input and labels match sizes 

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

806 

807 if index is None: 

808 mask = (labels > 0) 

809 masked_positions = None 

810 if find_positions: 

811 masked_positions = positions[mask] 

812 return single_group(input[mask], masked_positions) 

813 

814 if numpy.isscalar(index): 

815 mask = (labels == index) 

816 masked_positions = None 

817 if find_positions: 

818 masked_positions = positions[mask] 

819 return single_group(input[mask], masked_positions) 

820 

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

822 # label is larger than the number of values. 

823 if (not _safely_castable_to_int(labels.dtype) or 

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

825 # remap labels, and indexes 

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

827 idxs = numpy.searchsorted(unique_labels, index) 

828 

829 # make all of idxs valid 

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

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

832 else: 

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

834 idxs = numpy.asanyarray(index, numpy.int).copy() 

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

836 

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

838 

839 if find_median: 

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

841 else: 

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

843 input = input.ravel()[order] 

844 labels = labels.ravel()[order] 

845 if find_positions: 

846 positions = positions.ravel()[order] 

847 

848 result = [] 

849 if find_min: 

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

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

852 result += [mins[idxs]] 

853 if find_min_positions: 

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

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

856 result += [minpos[idxs]] 

857 if find_max: 

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

859 maxs[labels] = input 

860 result += [maxs[idxs]] 

861 if find_max_positions: 

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

863 maxpos[labels] = positions 

864 result += [maxpos[idxs]] 

865 if find_median: 

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

867 lo = numpy.zeros(labels.max() + 2, numpy.int) 

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

869 hi = numpy.zeros(labels.max() + 2, numpy.int) 

870 hi[labels] = locs 

871 lo = lo[idxs] 

872 hi = hi[idxs] 

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

874 # hi is an index to the largest value. 

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

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

877 step = (hi - lo) // 2 

878 lo += step 

879 hi -= step 

880 result += [(input[lo] + input[hi]) / 2.0] 

881 

882 return result 

883 

884 

885def minimum(input, labels=None, index=None): 

886 """ 

887 Calculate the minimum of the values of an array over labeled regions. 

888 

889 Parameters 

890 ---------- 

891 input : array_like 

892 Array_like of values. For each region specified by `labels`, the 

893 minimal values of `input` over the region is computed. 

894 labels : array_like, optional 

895 An array_like of integers marking different regions over which the 

896 minimum value of `input` is to be computed. `labels` must have the 

897 same shape as `input`. If `labels` is not specified, the minimum 

898 over the whole array is returned. 

899 index : array_like, optional 

900 A list of region labels that are taken into account for computing the 

901 minima. If index is None, the minimum over all elements where `labels` 

902 is non-zero is returned. 

903 

904 Returns 

905 ------- 

906 minimum : float or list of floats 

907 List of minima of `input` over the regions determined by `labels` and 

908 whose index is in `index`. If `index` or `labels` are not specified, a 

909 float is returned: the minimal value of `input` if `labels` is None, 

910 and the minimal value of elements where `labels` is greater than zero 

911 if `index` is None. 

912 

913 See also 

914 -------- 

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

916 standard_deviation 

917 

918 Notes 

919 ----- 

920 The function returns a Python list and not a NumPy array, use 

921 `np.array` to convert the list to an array. 

922 

923 Examples 

924 -------- 

925 >>> from scipy import ndimage 

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

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

928 ... [0, 0, 0, 7], 

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

930 >>> labels, labels_nb = ndimage.label(a) 

931 >>> labels 

932 array([[1, 1, 0, 0], 

933 [1, 1, 0, 2], 

934 [0, 0, 0, 2], 

935 [3, 3, 0, 0]]) 

936 >>> ndimage.minimum(a, labels=labels, index=np.arange(1, labels_nb + 1)) 

937 [1.0, 4.0, 3.0] 

938 >>> ndimage.minimum(a) 

939 0.0 

940 >>> ndimage.minimum(a, labels=labels) 

941 1.0 

942 

943 """ 

944 return _select(input, labels, index, find_min=True)[0] 

945 

946 

947def maximum(input, labels=None, index=None): 

948 """ 

949 Calculate the maximum of the values of an array over labeled regions. 

950 

951 Parameters 

952 ---------- 

953 input : array_like 

954 Array_like of values. For each region specified by `labels`, the 

955 maximal values of `input` over the region is computed. 

956 labels : array_like, optional 

957 An array of integers marking different regions over which the 

958 maximum value of `input` is to be computed. `labels` must have the 

959 same shape as `input`. If `labels` is not specified, the maximum 

960 over the whole array is returned. 

961 index : array_like, optional 

962 A list of region labels that are taken into account for computing the 

963 maxima. If index is None, the maximum over all elements where `labels` 

964 is non-zero is returned. 

965 

966 Returns 

967 ------- 

968 output : float or list of floats 

969 List of maxima of `input` over the regions determined by `labels` and 

970 whose index is in `index`. If `index` or `labels` are not specified, a 

971 float is returned: the maximal value of `input` if `labels` is None, 

972 and the maximal value of elements where `labels` is greater than zero 

973 if `index` is None. 

974 

975 See also 

976 -------- 

977 label, minimum, median, maximum_position, extrema, sum, mean, variance, 

978 standard_deviation 

979 

980 Notes 

981 ----- 

982 The function returns a Python list and not a NumPy array, use 

983 `np.array` to convert the list to an array. 

984 

985 Examples 

986 -------- 

987 >>> a = np.arange(16).reshape((4,4)) 

988 >>> a 

989 array([[ 0, 1, 2, 3], 

990 [ 4, 5, 6, 7], 

991 [ 8, 9, 10, 11], 

992 [12, 13, 14, 15]]) 

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

994 >>> labels[:2,:2] = 1 

995 >>> labels[2:, 1:3] = 2 

996 >>> labels 

997 array([[1, 1, 0, 0], 

998 [1, 1, 0, 0], 

999 [0, 2, 2, 0], 

1000 [0, 2, 2, 0]]) 

1001 >>> from scipy import ndimage 

1002 >>> ndimage.maximum(a) 

1003 15.0 

1004 >>> ndimage.maximum(a, labels=labels, index=[1,2]) 

1005 [5.0, 14.0] 

1006 >>> ndimage.maximum(a, labels=labels) 

1007 14.0 

1008 

1009 >>> b = np.array([[1, 2, 0, 0], 

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

1011 ... [0, 0, 0, 7], 

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

1013 >>> labels, labels_nb = ndimage.label(b) 

1014 >>> labels 

1015 array([[1, 1, 0, 0], 

1016 [1, 1, 0, 2], 

1017 [0, 0, 0, 2], 

1018 [3, 3, 0, 0]]) 

1019 >>> ndimage.maximum(b, labels=labels, index=np.arange(1, labels_nb + 1)) 

1020 [5.0, 7.0, 9.0] 

1021 

1022 """ 

1023 return _select(input, labels, index, find_max=True)[0] 

1024 

1025 

1026def median(input, labels=None, index=None): 

1027 """ 

1028 Calculate the median of the values of an array over labeled regions. 

1029 

1030 Parameters 

1031 ---------- 

1032 input : array_like 

1033 Array_like of values. For each region specified by `labels`, the 

1034 median value of `input` over the region is computed. 

1035 labels : array_like, optional 

1036 An array_like of integers marking different regions over which the 

1037 median value of `input` is to be computed. `labels` must have the 

1038 same shape as `input`. If `labels` is not specified, the median 

1039 over the whole array is returned. 

1040 index : array_like, optional 

1041 A list of region labels that are taken into account for computing the 

1042 medians. If index is None, the median over all elements where `labels` 

1043 is non-zero is returned. 

1044 

1045 Returns 

1046 ------- 

1047 median : float or list of floats 

1048 List of medians of `input` over the regions determined by `labels` and 

1049 whose index is in `index`. If `index` or `labels` are not specified, a 

1050 float is returned: the median value of `input` if `labels` is None, 

1051 and the median value of elements where `labels` is greater than zero 

1052 if `index` is None. 

1053 

1054 See also 

1055 -------- 

1056 label, minimum, maximum, extrema, sum, mean, variance, standard_deviation 

1057 

1058 Notes 

1059 ----- 

1060 The function returns a Python list and not a NumPy array, use 

1061 `np.array` to convert the list to an array. 

1062 

1063 Examples 

1064 -------- 

1065 >>> from scipy import ndimage 

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

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

1068 ... [0, 0, 0, 7], 

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

1070 >>> labels, labels_nb = ndimage.label(a) 

1071 >>> labels 

1072 array([[1, 1, 0, 2], 

1073 [1, 1, 0, 2], 

1074 [0, 0, 0, 2], 

1075 [3, 3, 0, 0]]) 

1076 >>> ndimage.median(a, labels=labels, index=np.arange(1, labels_nb + 1)) 

1077 [2.5, 4.0, 6.0] 

1078 >>> ndimage.median(a) 

1079 1.0 

1080 >>> ndimage.median(a, labels=labels) 

1081 3.0 

1082 

1083 """ 

1084 return _select(input, labels, index, find_median=True)[0] 

1085 

1086 

1087def minimum_position(input, labels=None, index=None): 

1088 """ 

1089 Find the positions of the minimums of the values of an array at labels. 

1090 

1091 Parameters 

1092 ---------- 

1093 input : array_like 

1094 Array_like of values. 

1095 labels : array_like, optional 

1096 An array of integers marking different regions over which the 

1097 position of the minimum value of `input` is to be computed. 

1098 `labels` must have the same shape as `input`. If `labels` is not 

1099 specified, the location of the first minimum over the whole 

1100 array is returned. 

1101 

1102 The `labels` argument only works when `index` is specified. 

1103 index : array_like, optional 

1104 A list of region labels that are taken into account for finding the 

1105 location of the minima. If `index` is None, the ``first`` minimum 

1106 over all elements where `labels` is non-zero is returned. 

1107 

1108 The `index` argument only works when `labels` is specified. 

1109 

1110 Returns 

1111 ------- 

1112 output : list of tuples of ints 

1113 Tuple of ints or list of tuples of ints that specify the location 

1114 of minima of `input` over the regions determined by `labels` and 

1115 whose index is in `index`. 

1116 

1117 If `index` or `labels` are not specified, a tuple of ints is 

1118 returned specifying the location of the first minimal value of `input`. 

1119 

1120 See also 

1121 -------- 

1122 label, minimum, median, maximum_position, extrema, sum, mean, variance, 

1123 standard_deviation 

1124 

1125 Examples 

1126 -------- 

1127 >>> a = np.array([[10, 20, 30], 

1128 ... [40, 80, 100], 

1129 ... [1, 100, 200]]) 

1130 >>> b = np.array([[1, 2, 0, 1], 

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

1132 ... [0, 0, 0, 7], 

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

1134 

1135 >>> from scipy import ndimage 

1136 

1137 >>> ndimage.minimum_position(a) 

1138 (2, 0) 

1139 >>> ndimage.minimum_position(b) 

1140 (0, 2) 

1141 

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

1143 

1144 >>> label, pos = ndimage.label(a) 

1145 >>> ndimage.minimum_position(a, label, index=np.arange(1, pos+1)) 

1146 [(2, 0)] 

1147 

1148 >>> label, pos = ndimage.label(b) 

1149 >>> ndimage.minimum_position(b, label, index=np.arange(1, pos+1)) 

1150 [(0, 0), (0, 3), (3, 1)] 

1151 

1152 """ 

1153 dims = numpy.array(numpy.asarray(input).shape) 

1154 # see numpy.unravel_index to understand this line. 

1155 dim_prod = numpy.cumprod([1] + list(dims[:0:-1]))[::-1] 

1156 

1157 result = _select(input, labels, index, find_min_positions=True)[0] 

1158 

1159 if numpy.isscalar(result): 

1160 return tuple((result // dim_prod) % dims) 

1161 

1162 return [tuple(v) for v in (result.reshape(-1, 1) // dim_prod) % dims] 

1163 

1164 

1165def maximum_position(input, labels=None, index=None): 

1166 """ 

1167 Find the positions of the maximums of the values of an array at labels. 

1168 

1169 For each region specified by `labels`, the position of the maximum 

1170 value of `input` within the region is returned. 

1171 

1172 Parameters 

1173 ---------- 

1174 input : array_like 

1175 Array_like of values. 

1176 labels : array_like, optional 

1177 An array of integers marking different regions over which the 

1178 position of the maximum value of `input` is to be computed. 

1179 `labels` must have the same shape as `input`. If `labels` is not 

1180 specified, the location of the first maximum over the whole 

1181 array is returned. 

1182 

1183 The `labels` argument only works when `index` is specified. 

1184 index : array_like, optional 

1185 A list of region labels that are taken into account for finding the 

1186 location of the maxima. If `index` is None, the first maximum 

1187 over all elements where `labels` is non-zero is returned. 

1188 

1189 The `index` argument only works when `labels` is specified. 

1190 

1191 Returns 

1192 ------- 

1193 output : list of tuples of ints 

1194 List of tuples of ints that specify the location of maxima of 

1195 `input` over the regions determined by `labels` and whose index 

1196 is in `index`. 

1197 

1198 If `index` or `labels` are not specified, a tuple of ints is 

1199 returned specifying the location of the ``first`` maximal value 

1200 of `input`. 

1201 

1202 See also 

1203 -------- 

1204 label, minimum, median, maximum_position, extrema, sum, mean, variance, 

1205 standard_deviation 

1206  

1207 Examples 

1208 -------- 

1209 >>> from scipy import ndimage 

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

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

1212 ... [0, 0, 0, 7], 

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

1214 >>> ndimage.maximum_position(a) 

1215 (3, 0) 

1216 

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

1218 

1219 >>> lbl = np.array([[0, 1, 2, 3], 

1220 ... [0, 1, 2, 3], 

1221 ... [0, 1, 2, 3], 

1222 ... [0, 1, 2, 3]]) 

1223 >>> ndimage.maximum_position(a, lbl, 1) 

1224 (1, 1) 

1225  

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

1227 

1228 >>> ndimage.maximum_position(a, lbl) 

1229 (2, 3) 

1230  

1231 If there are no maxima, the position of the first element is returned: 

1232  

1233 >>> ndimage.maximum_position(a, lbl, 2) 

1234 (0, 2) 

1235 

1236 """ 

1237 dims = numpy.array(numpy.asarray(input).shape) 

1238 # see numpy.unravel_index to understand this line. 

1239 dim_prod = numpy.cumprod([1] + list(dims[:0:-1]))[::-1] 

1240 

1241 result = _select(input, labels, index, find_max_positions=True)[0] 

1242 

1243 if numpy.isscalar(result): 

1244 return tuple((result // dim_prod) % dims) 

1245 

1246 return [tuple(v) for v in (result.reshape(-1, 1) // dim_prod) % dims] 

1247 

1248 

1249def extrema(input, labels=None, index=None): 

1250 """ 

1251 Calculate the minimums and maximums of the values of an array 

1252 at labels, along with their positions. 

1253 

1254 Parameters 

1255 ---------- 

1256 input : ndarray 

1257 N-D image data to process. 

1258 labels : ndarray, optional 

1259 Labels of features in input. 

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

1261 index : int or sequence of ints, optional 

1262 Labels to include in output. If None (default), all values where 

1263 non-zero `labels` are used. 

1264 

1265 Returns 

1266 ------- 

1267 minimums, maximums : int or ndarray 

1268 Values of minimums and maximums in each feature. 

1269 min_positions, max_positions : tuple or list of tuples 

1270 Each tuple gives the N-D coordinates of the corresponding minimum 

1271 or maximum. 

1272 

1273 See Also 

1274 -------- 

1275 maximum, minimum, maximum_position, minimum_position, center_of_mass 

1276 

1277 Examples 

1278 -------- 

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

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

1281 ... [0, 0, 0, 7], 

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

1283 >>> from scipy import ndimage 

1284 >>> ndimage.extrema(a) 

1285 (0, 9, (0, 2), (3, 0)) 

1286 

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

1288 

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

1290 >>> ndimage.extrema(a, lbl, index=np.arange(1, nlbl+1)) 

1291 (array([1, 4, 3]), 

1292 array([5, 7, 9]), 

1293 [(0, 0), (1, 3), (3, 1)], 

1294 [(1, 0), (2, 3), (3, 0)]) 

1295 

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

1297 

1298 >>> ndimage.extrema(a, lbl) 

1299 (1, 9, (0, 0), (3, 0)) 

1300 

1301 """ 

1302 dims = numpy.array(numpy.asarray(input).shape) 

1303 # see numpy.unravel_index to understand this line. 

1304 dim_prod = numpy.cumprod([1] + list(dims[:0:-1]))[::-1] 

1305 

1306 minimums, min_positions, maximums, max_positions = _select(input, labels, 

1307 index, 

1308 find_min=True, 

1309 find_max=True, 

1310 find_min_positions=True, 

1311 find_max_positions=True) 

1312 

1313 if numpy.isscalar(minimums): 

1314 return (minimums, maximums, tuple((min_positions // dim_prod) % dims), 

1315 tuple((max_positions // dim_prod) % dims)) 

1316 

1317 min_positions = [tuple(v) for v in (min_positions.reshape(-1, 1) // dim_prod) % dims] 

1318 max_positions = [tuple(v) for v in (max_positions.reshape(-1, 1) // dim_prod) % dims] 

1319 

1320 return minimums, maximums, min_positions, max_positions 

1321 

1322 

1323def center_of_mass(input, labels=None, index=None): 

1324 """ 

1325 Calculate the center of mass of the values of an array at labels. 

1326 

1327 Parameters 

1328 ---------- 

1329 input : ndarray 

1330 Data from which to calculate center-of-mass. The masses can either 

1331 be positive or negative. 

1332 labels : ndarray, optional 

1333 Labels for objects in `input`, as generated by `ndimage.label`. 

1334 Only used with `index`. Dimensions must be the same as `input`. 

1335 index : int or sequence of ints, optional 

1336 Labels for which to calculate centers-of-mass. If not specified, 

1337 all labels greater than zero are used. Only used with `labels`. 

1338 

1339 Returns 

1340 ------- 

1341 center_of_mass : tuple, or list of tuples 

1342 Coordinates of centers-of-mass. 

1343 

1344 Examples 

1345 -------- 

1346 >>> a = np.array(([0,0,0,0], 

1347 ... [0,1,1,0], 

1348 ... [0,1,1,0], 

1349 ... [0,1,1,0])) 

1350 >>> from scipy import ndimage 

1351 >>> ndimage.measurements.center_of_mass(a) 

1352 (2.0, 1.5) 

1353 

1354 Calculation of multiple objects in an image 

1355 

1356 >>> b = np.array(([0,1,1,0], 

1357 ... [0,1,0,0], 

1358 ... [0,0,0,0], 

1359 ... [0,0,1,1], 

1360 ... [0,0,1,1])) 

1361 >>> lbl = ndimage.label(b)[0] 

1362 >>> ndimage.measurements.center_of_mass(b, lbl, [1,2]) 

1363 [(0.33333333333333331, 1.3333333333333333), (3.5, 2.5)] 

1364 

1365 Negative masses are also accepted, which can occur for example when 

1366 bias is removed from measured data due to random noise. 

1367 

1368 >>> c = np.array(([-1,0,0,0], 

1369 ... [0,-1,-1,0], 

1370 ... [0,1,-1,0], 

1371 ... [0,1,1,0])) 

1372 >>> ndimage.measurements.center_of_mass(c) 

1373 (-4.0, 1.0) 

1374 

1375 If there are division by zero issues, the function does not raise an 

1376 error but rather issues a RuntimeWarning before returning inf and/or NaN. 

1377 

1378 >>> d = np.array([-1, 1]) 

1379 >>> ndimage.measurements.center_of_mass(d) 

1380 (inf,) 

1381 """ 

1382 normalizer = sum(input, labels, index) 

1383 grids = numpy.ogrid[[slice(0, i) for i in input.shape]] 

1384 

1385 results = [sum(input * grids[dir].astype(float), labels, index) / normalizer 

1386 for dir in range(input.ndim)] 

1387 

1388 if numpy.isscalar(results[0]): 

1389 return tuple(results) 

1390 

1391 return [tuple(v) for v in numpy.array(results).T] 

1392 

1393 

1394def histogram(input, min, max, bins, labels=None, index=None): 

1395 """ 

1396 Calculate the histogram of the values of an array, optionally at labels. 

1397 

1398 Histogram calculates the frequency of values in an array within bins 

1399 determined by `min`, `max`, and `bins`. The `labels` and `index` 

1400 keywords can limit the scope of the histogram to specified sub-regions 

1401 within the array. 

1402 

1403 Parameters 

1404 ---------- 

1405 input : array_like 

1406 Data for which to calculate histogram. 

1407 min, max : int 

1408 Minimum and maximum values of range of histogram bins. 

1409 bins : int 

1410 Number of bins. 

1411 labels : array_like, optional 

1412 Labels for objects in `input`. 

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

1414 index : int or sequence of ints, optional 

1415 Label or labels for which to calculate histogram. If None, all values 

1416 where label is greater than zero are used 

1417 

1418 Returns 

1419 ------- 

1420 hist : ndarray 

1421 Histogram counts. 

1422 

1423 Examples 

1424 -------- 

1425 >>> a = np.array([[ 0. , 0.2146, 0.5962, 0. ], 

1426 ... [ 0. , 0.7778, 0. , 0. ], 

1427 ... [ 0. , 0. , 0. , 0. ], 

1428 ... [ 0. , 0. , 0.7181, 0.2787], 

1429 ... [ 0. , 0. , 0.6573, 0.3094]]) 

1430 >>> from scipy import ndimage 

1431 >>> ndimage.measurements.histogram(a, 0, 1, 10) 

1432 array([13, 0, 2, 1, 0, 1, 1, 2, 0, 0]) 

1433 

1434 With labels and no indices, non-zero elements are counted: 

1435 

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

1437 >>> ndimage.measurements.histogram(a, 0, 1, 10, lbl) 

1438 array([0, 0, 2, 1, 0, 1, 1, 2, 0, 0]) 

1439 

1440 Indices can be used to count only certain objects: 

1441 

1442 >>> ndimage.measurements.histogram(a, 0, 1, 10, lbl, 2) 

1443 array([0, 0, 1, 1, 0, 0, 1, 1, 0, 0]) 

1444 

1445 """ 

1446 _bins = numpy.linspace(min, max, bins + 1) 

1447 

1448 def _hist(vals): 

1449 return numpy.histogram(vals, _bins)[0] 

1450 

1451 return labeled_comprehension(input, labels, index, _hist, object, None, 

1452 pass_positions=False) 

1453 

1454 

1455def watershed_ift(input, markers, structure=None, output=None): 

1456 """ 

1457 Apply watershed from markers using image foresting transform algorithm. 

1458 

1459 Parameters 

1460 ---------- 

1461 input : array_like 

1462 Input. 

1463 markers : array_like 

1464 Markers are points within each watershed that form the beginning 

1465 of the process. Negative markers are considered background markers 

1466 which are processed after the other markers. 

1467 structure : structure element, optional 

1468 A structuring element defining the connectivity of the object can be 

1469 provided. If None, an element is generated with a squared 

1470 connectivity equal to one. 

1471 output : ndarray, optional 

1472 An output array can optionally be provided. The same shape as input. 

1473 

1474 Returns 

1475 ------- 

1476 watershed_ift : ndarray 

1477 Output. Same shape as `input`. 

1478 

1479 References 

1480 ---------- 

1481 .. [1] A.X. Falcao, J. Stolfi and R. de Alencar Lotufo, "The image 

1482 foresting transform: theory, algorithms, and applications", 

1483 Pattern Analysis and Machine Intelligence, vol. 26, pp. 19-29, 2004. 

1484 

1485 """ 

1486 input = numpy.asarray(input) 

1487 if input.dtype.type not in [numpy.uint8, numpy.uint16]: 

1488 raise TypeError('only 8 and 16 unsigned inputs are supported') 

1489 

1490 if structure is None: 

1491 structure = morphology.generate_binary_structure(input.ndim, 1) 

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

1493 if structure.ndim != input.ndim: 

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

1495 for ii in structure.shape: 

1496 if ii != 3: 

1497 raise RuntimeError('structure dimensions must be equal to 3') 

1498 

1499 if not structure.flags.contiguous: 

1500 structure = structure.copy() 

1501 markers = numpy.asarray(markers) 

1502 if input.shape != markers.shape: 

1503 raise RuntimeError('input and markers must have equal shape') 

1504 

1505 integral_types = [numpy.int0, 

1506 numpy.int8, 

1507 numpy.int16, 

1508 numpy.int32, 

1509 numpy.int_, 

1510 numpy.int64, 

1511 numpy.intc, 

1512 numpy.intp] 

1513 

1514 if markers.dtype.type not in integral_types: 

1515 raise RuntimeError('marker should be of integer type') 

1516 

1517 if isinstance(output, numpy.ndarray): 

1518 if output.dtype.type not in integral_types: 

1519 raise RuntimeError('output should be of integer type') 

1520 else: 

1521 output = markers.dtype 

1522 

1523 output = _ni_support._get_output(output, input) 

1524 _nd_image.watershed_ift(input, markers, structure, output) 

1525 return output