Package ndimage :: Module morphology
[hide private]
[frames] | no frames]

Source Code for Module ndimage.morphology

  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   
 31  from __future__ import division 
 32  import numpy 
 33  import _ni_support 
 34  import _nd_image 
 35  import filters 
 36  import types 
 37   
 38   
39 -def _center_is_true(structure, origin):
40 structure = numpy.array(structure) 41 coor = tuple([oo + ss // 2 for ss, oo in zip(structure.shape, 42 origin)]) 43 return bool(structure[coor])
44
45 -def iterate_structure(structure, iterations, origin = None):
46 """Iterate a structure by dilating it with itself. 47 48 If origin is None, only the iterated structure is returned. If 49 not, a tuple of the iterated structure and the modified origin is 50 returned. 51 """ 52 structure = numpy.asarray(structure) 53 if iterations < 2: 54 return structure.copy() 55 ni = iterations - 1 56 shape = [ii + ni * (ii - 1) for ii in structure.shape] 57 pos = [ni * (structure.shape[ii] / 2) for ii in range(len(shape))] 58 slc = [slice(pos[ii], pos[ii] + structure.shape[ii], None) 59 for ii in range(len(shape))] 60 out = numpy.zeros(shape, bool) 61 out[slc] = structure != 0 62 out = binary_dilation(out, structure, iterations = ni) 63 if origin is None: 64 return out 65 else: 66 origin = _ni_support._normalize_sequence(origin, structure.ndim) 67 origin = [iterations * o for o in origin] 68 return out, origin
69
70 -def generate_binary_structure(rank, connectivity):
71 """Generate a binary structure for binary morphological operations. 72 73 The inputs are the rank of the array to which the structure will 74 be applied and the square of the connectivity of the structure. 75 """ 76 if connectivity < 1: 77 connectivity = 1 78 if rank < 1: 79 if connectivity < 1: 80 return numpy.array(0, dtype = bool) 81 else: 82 return numpy.array(1, dtype = bool) 83 output = numpy.zeros([3] * rank, bool) 84 output = numpy.fabs(numpy.indices([3] * rank) - 1) 85 output = numpy.add.reduce(output, 0) 86 return numpy.asarray(output <= connectivity, dtype = bool)
87 88
89 -def _binary_erosion(input, structure, iterations, mask, output, 90 border_value, origin, invert, brute_force):
91 input = numpy.asarray(input) 92 if numpy.iscomplexobj(input): 93 raise TypeError, 'Complex type not supported' 94 if structure is None: 95 structure = generate_binary_structure(input.ndim, 1) 96 else: 97 structure = numpy.asarray(structure) 98 structure = structure.astype(bool) 99 if structure.ndim != input.ndim: 100 raise RuntimeError, 'structure rank must equal input rank' 101 if not structure.flags.contiguous: 102 structure = structure.copy() 103 if numpy.product(structure.shape,axis=0) < 1: 104 raise RuntimeError, 'structure must not be empty' 105 if mask is not None: 106 mask = numpy.asarray(mask) 107 if mask.shape != input.shape: 108 raise RuntimeError, 'mask and input must have equal sizes' 109 origin = _ni_support._normalize_sequence(origin, input.ndim) 110 cit = _center_is_true(structure, origin) 111 if isinstance(output, numpy.ndarray): 112 if numpy.iscomplexobj(output): 113 raise TypeError, 'Complex output type not supported' 114 else: 115 output = bool 116 output, return_value = _ni_support._get_output(output, input) 117 118 119 if iterations == 1: 120 _nd_image.binary_erosion(input, structure, mask, output, 121 border_value, origin, invert, cit, 0) 122 return return_value 123 elif cit and not brute_force: 124 changed, coordinate_list = _nd_image.binary_erosion(input, 125 structure, mask, output, border_value, origin, invert, cit, 1) 126 structure = structure[tuple([slice(None, None, -1)] * 127 structure.ndim)] 128 for ii in range(len(origin)): 129 origin[ii] = -origin[ii] 130 if not structure.shape[ii] & 1: 131 origin[ii] -= 1 132 if mask is not None: 133 msk = numpy.asarray(mask) 134 msk = mask.astype(numpy.int8) 135 if msk is mask: 136 msk = mask.copy() 137 mask = msk 138 if not structure.flags.contiguous: 139 structure = structure.copy() 140 _nd_image.binary_erosion2(output, structure, mask, iterations - 1, 141 origin, invert, coordinate_list) 142 return return_value 143 else: 144 tmp_in = numpy.zeros(input.shape, bool) 145 if return_value is None: 146 tmp_out = output 147 else: 148 tmp_out = numpy.zeros(input.shape, bool) 149 if not iterations & 1: 150 tmp_in, tmp_out = tmp_out, tmp_in 151 changed = _nd_image.binary_erosion(input, structure, mask, 152 tmp_out, border_value, origin, invert, cit, 0) 153 ii = 1 154 while (ii < iterations) or (iterations < 1) and changed: 155 tmp_in, tmp_out = tmp_out, tmp_in 156 changed = _nd_image.binary_erosion(tmp_in, structure, mask, 157 tmp_out, border_value, origin, invert, cit, 0) 158 ii += 1 159 if return_value is not None: 160 return tmp_out
161 162
163 -def binary_erosion(input, structure = None, iterations = 1, mask = None, 164 output = None, border_value = 0, origin = 0, brute_force = False):
165 """Multi-dimensional binary erosion with the given structure. 166 167 An output array can optionally be provided. The origin parameter 168 controls the placement of the filter. If no structuring element is 169 provided an element is generated with a squared connectivity equal 170 to one. The border_value parameter gives the value of the array 171 outside the border. The erosion operation is repeated iterations 172 times. If iterations is less than 1, the erosion is repeated until 173 the result does not change anymore. If a mask is given, only those 174 elements with a true value at the corresponding mask element are 175 modified at each iteration. 176 """ 177 return _binary_erosion(input, structure, iterations, mask, 178 output, border_value, origin, 0, brute_force)
179
180 -def binary_dilation(input, structure = None, iterations = 1, mask = None, 181 output = None, border_value = 0, origin = 0, brute_force = False):
182 """Multi-dimensional binary dilation with the given structure. 183 184 An output array can optionally be provided. The origin parameter 185 controls the placement of the filter. If no structuring element is 186 provided an element is generated with a squared connectivity equal 187 to one. The dilation operation is repeated iterations times. If 188 iterations is less than 1, the dilation is repeated until the 189 result does not change anymore. If a mask is given, only those 190 elements with a true value at the corresponding mask element are 191 modified at each iteration. 192 """ 193 input = numpy.asarray(input) 194 if structure is None: 195 structure = generate_binary_structure(input.ndim, 1) 196 origin = _ni_support._normalize_sequence(origin, input.ndim) 197 structure = numpy.asarray(structure) 198 structure = structure[tuple([slice(None, None, -1)] * 199 structure.ndim)] 200 for ii in range(len(origin)): 201 origin[ii] = -origin[ii] 202 if not structure.shape[ii] & 1: 203 origin[ii] -= 1 204 return _binary_erosion(input, structure, iterations, mask, 205 output, border_value, origin, 1, brute_force)
206 207
208 -def binary_opening(input, structure = None, iterations = 1, output = None, 209 origin = 0):
210 """Multi-dimensional binary opening with the given structure. 211 212 An output array can optionally be provided. The origin parameter 213 controls the placement of the filter. If no structuring element is 214 provided an element is generated with a squared connectivity equal 215 to one. The iterations parameter gives the number of times the 216 erosions and then the dilations are done. 217 """ 218 input = numpy.asarray(input) 219 if structure is None: 220 rank = input.ndim 221 structure = generate_binary_structure(rank, 1) 222 tmp = binary_erosion(input, structure, iterations, None, None, 0, 223 origin) 224 return binary_dilation(tmp, structure, iterations, None, output, 0, 225 origin)
226 227
228 -def binary_closing(input, structure = None, iterations = 1, output = None, 229 origin = 0):
230 """Multi-dimensional binary closing with the given structure. 231 232 An output array can optionally be provided. The origin parameter 233 controls the placement of the filter. If no structuring element is 234 provided an element is generated with a squared connectivity equal 235 to one. The iterations parameter gives the number of times the 236 dilations and then the erosions are done. 237 """ 238 input = numpy.asarray(input) 239 if structure is None: 240 rank = input.ndim 241 structure = generate_binary_structure(rank, 1) 242 tmp = binary_dilation(input, structure, iterations, None, None, 0, 243 origin) 244 return binary_erosion(tmp, structure, iterations, None, output, 0, 245 origin)
246 247
248 -def binary_hit_or_miss(input, structure1 = None, structure2 = None, 249 output = None, origin1 = 0, origin2 = None):
250 """Multi-dimensional binary hit-or-miss transform. 251 252 An output array can optionally be provided. The origin parameters 253 controls the placement of the structuring elements. If the first 254 structuring element is not given one is generated with a squared 255 connectivity equal to one. If the second structuring element is 256 not provided, it set equal to the inverse of the first structuring 257 element. If the origin for the second structure is equal to None 258 it is set equal to the origin of the first. 259 """ 260 input = numpy.asarray(input) 261 if structure1 is None: 262 structure1 = generate_binary_structure(input.ndim, 1) 263 if structure2 is None: 264 structure2 = numpy.logical_not(structure1) 265 origin1 = _ni_support._normalize_sequence(origin1, input.ndim) 266 if origin2 is None: 267 origin2 = origin1 268 else: 269 origin2 = _ni_support._normalize_sequence(origin2, input.ndim) 270 271 tmp1 = _binary_erosion(input, structure1, 1, None, None, 0, origin1, 272 0, False) 273 inplace = isinstance(output, numpy.ndarray) 274 result = _binary_erosion(input, structure2, 1, None, output, 0, 275 origin2, 1, False) 276 if inplace: 277 numpy.logical_not(output, output) 278 numpy.logical_and(tmp1, output, output) 279 else: 280 numpy.logical_not(result, result) 281 return numpy.logical_and(tmp1, result)
282
283 -def binary_propagation(input, structure = None, mask = None, 284 output = None, border_value = 0, origin = 0):
285 """Multi-dimensional binary propagation with the given structure. 286 287 An output array can optionally be provided. The origin parameter 288 controls the placement of the filter. If no structuring element is 289 provided an element is generated with a squared connectivity equal 290 to one. If a mask is given, only those elements with a true value at 291 the corresponding mask element are. 292 293 This function is functionally equivalent to calling binary_dilation 294 with the number of iterations less then one: iterative dilation until 295 the result does not change anymore. 296 """ 297 return binary_dilation(input, structure, -1, mask, output, 298 border_value, origin)
299
300 -def binary_fill_holes(input, structure = None, output = None, origin = 0):
301 """Fill the holes in binary objects. 302 303 An output array can optionally be provided. The origin parameter 304 controls the placement of the filter. If no structuring element is 305 provided an element is generated with a squared connectivity equal 306 to one. 307 """ 308 mask = numpy.logical_not(input) 309 tmp = numpy.zeros(mask.shape, bool) 310 inplace = isinstance(output, numpy.ndarray) 311 if inplace: 312 binary_dilation(tmp, structure, -1, mask, output, 1, origin) 313 numpy.logical_not(output, output) 314 else: 315 output = binary_dilation(tmp, structure, -1, mask, None, 1, 316 origin) 317 numpy.logical_not(output, output) 318 return output
319
320 -def grey_erosion(input, size = None, footprint = None, structure = None, 321 output = None, mode = "reflect", cval = 0.0, origin = 0):
322 """Calculate a grey values erosion. 323 324 Either a size or a footprint, or the structure must be provided. An 325 output array can optionally be provided. The origin parameter 326 controls the placement of the filter. The mode parameter 327 determines how the array borders are handled, where cval is the 328 value when mode is equal to 'constant'. 329 """ 330 return filters._min_or_max_filter(input, size, footprint, structure, 331 output, mode, cval, origin, 1)
332 333
334 -def grey_dilation(input, size = None, footprint = None, structure = None, 335 output = None, mode = "reflect", cval = 0.0, origin = 0):
336 """Calculate a grey values dilation. 337 338 Either a size or a footprint, or the structure must be 339 provided. An output array can optionally be provided. The origin 340 parameter controls the placement of the filter. The mode parameter 341 determines how the array borders are handled, where cval is the 342 value when mode is equal to 'constant'. 343 """ 344 if structure is not None: 345 structure = numpy.asarray(structure) 346 structure = structure[tuple([slice(None, None, -1)] * 347 structure.ndim)] 348 if footprint is not None: 349 footprint = numpy.asarray(footprint) 350 footprint = footprint[tuple([slice(None, None, -1)] * 351 footprint.ndim)] 352 input = numpy.asarray(input) 353 origin = _ni_support._normalize_sequence(origin, input.ndim) 354 for ii in range(len(origin)): 355 origin[ii] = -origin[ii] 356 if footprint is not None: 357 sz = footprint.shape[ii] 358 else: 359 sz = size[ii] 360 if not sz & 1: 361 origin[ii] -= 1 362 return filters._min_or_max_filter(input, size, footprint, structure, 363 output, mode, cval, origin, 0)
364 365
366 -def grey_opening(input, size = None, footprint = None, structure = None, 367 output = None, mode = "reflect", cval = 0.0, origin = 0):
368 """Multi-dimensional grey valued opening. 369 370 Either a size or a footprint, or the structure must be provided. An 371 output array can optionally be provided. The origin parameter 372 controls the placement of the filter. The mode parameter 373 determines how the array borders are handled, where cval is the 374 value when mode is equal to 'constant'. 375 """ 376 tmp = grey_erosion(input, size, footprint, structure, None, mode, 377 cval, origin) 378 return grey_dilation(tmp, size, footprint, structure, output, mode, 379 cval, origin)
380 381
382 -def grey_closing(input, size = None, footprint = None, structure = None, 383 output = None, mode = "reflect", cval = 0.0, origin = 0):
384 """Multi-dimensional grey valued closing. 385 386 Either a size or a footprint, or the structure must be provided. An 387 output array can optionally be provided. The origin parameter 388 controls the placement of the filter. The mode parameter 389 determines how the array borders are handled, where cval is the 390 value when mode is equal to 'constant'. 391 """ 392 tmp = grey_dilation(input, size, footprint, structure, None, mode, 393 cval, origin) 394 return grey_erosion(tmp, size, footprint, structure, output, mode, 395 cval, origin)
396 397
398 -def morphological_gradient(input, size = None, footprint = None, 399 structure = None, output = None, mode = "reflect", 400 cval = 0.0, origin = 0):
401 """Multi-dimensional morphological gradient. 402 403 Either a size or a footprint, or the structure must be provided. An 404 output array can optionally be provided. The origin parameter 405 controls the placement of the filter. The mode parameter 406 determines how the array borders are handled, where cval is the 407 value when mode is equal to 'constant'. 408 """ 409 tmp = grey_dilation(input, size, footprint, structure, None, mode, 410 cval, origin) 411 if isinstance(output, numpy.ndarray): 412 grey_erosion(input, size, footprint, structure, output, mode, 413 cval, origin) 414 return numpy.subtract(tmp, output, output) 415 else: 416 return (tmp - grey_erosion(input, size, footprint, structure, 417 None, mode, cval, origin))
418 419
420 -def morphological_laplace(input, size = None, footprint = None, 421 structure = None, output = None, 422 mode = "reflect", cval = 0.0, origin = 0):
423 """Multi-dimensional morphological laplace. 424 425 Either a size or a footprint, or the structure must be provided. An 426 output array can optionally be provided. The origin parameter 427 controls the placement of the filter. The mode parameter 428 determines how the array borders are handled, where cval is the 429 value when mode is equal to 'constant'. 430 """ 431 tmp1 = grey_dilation(input, size, footprint, structure, None, mode, 432 cval, origin) 433 if isinstance(output, numpy.ndarray): 434 grey_erosion(input, size, footprint, structure, output, mode, 435 cval, origin) 436 numpy.add(tmp1, output, output) 437 del tmp1 438 numpy.subtract(output, input, output) 439 return numpy.subtract(output, input, output) 440 else: 441 tmp2 = grey_erosion(input, size, footprint, structure, None, mode, 442 cval, origin) 443 numpy.add(tmp1, tmp2, tmp2) 444 del tmp1 445 numpy.subtract(tmp2, input, tmp2) 446 numpy.subtract(tmp2, input, tmp2) 447 return tmp2
448 449
450 -def white_tophat(input, size = None, footprint = None, structure = None, 451 output = None, mode = "reflect", cval = 0.0, origin = 0):
452 """Multi-dimensional white tophat filter. 453 454 Either a size or a footprint, or the structure must be provided. An 455 output array can optionally be provided. The origin parameter 456 controls the placement of the filter. The mode parameter 457 determines how the array borders are handled, where cval is the 458 value when mode is equal to 'constant'. 459 """ 460 tmp = grey_erosion(input, size, footprint, structure, None, mode, 461 cval, origin) 462 if isinstance(output, numpy.ndarray): 463 grey_dilation(tmp, size, footprint, structure, output, mode, cval, 464 origin) 465 del tmp 466 return numpy.subtract(input, output, output) 467 else: 468 tmp = grey_dilation(tmp, size, footprint, structure, None, mode, 469 cval, origin) 470 return input - tmp
471 472
473 -def black_tophat(input, size = None, footprint = None, 474 structure = None, output = None, mode = "reflect", 475 cval = 0.0, origin = 0):
476 """Multi-dimensional black tophat filter. 477 478 Either a size or a footprint, or the structure must be provided. An 479 output array can optionally be provided. The origin parameter 480 controls the placement of the filter. The mode parameter 481 determines how the array borders are handled, where cval is the 482 value when mode is equal to 'constant'. 483 """ 484 tmp = grey_dilation(input, size, footprint, structure, None, mode, 485 cval, origin) 486 if isinstance(output, numpy.ndarray): 487 grey_erosion(tmp, size, footprint, structure, output, mode, cval, 488 origin) 489 del tmp 490 return numpy.subtract(output, input, output) 491 else: 492 tmp = grey_erosion(tmp, size, footprint, structure, None, mode, 493 cval, origin) 494 return tmp - input
495 496
497 -def distance_transform_bf(input, metric = "euclidean", sampling = None, 498 return_distances = True, return_indices = False, 499 distances = None, indices = None):
500 """Distance transform function by a brute force algorithm. 501 502 This function calculates the distance transform of the input, by 503 replacing each background element (zero values), with its 504 shortest distance to the foreground (any element non-zero). Three 505 types of distance metric are supported: 'euclidean', 'city_block' 506 and 'chessboard'. 507 508 In addition to the distance transform, the feature transform can 509 be calculated. In this case the index of the closest background 510 element is returned along the first axis of the result. 511 512 The return_distances, and return_indices flags can be used to 513 indicate if the distance transform, the feature transform, or both 514 must be returned. 515 516 Optionally the sampling along each axis can be given by the 517 sampling parameter which should be a sequence of length equal to 518 the input rank, or a single number in which the sampling is assumed 519 to be equal along all axes. This parameter is only used in the 520 case of the euclidean distance transform. 521 522 This function employs a slow brute force algorithm, see also the 523 function distance_transform_cdt for more efficient city_block and 524 chessboard algorithms. 525 526 the distances and indices arguments can be used to give optional 527 output arrays that must be of the correct size and type (float64 528 and int32). 529 """ 530 if (not return_distances) and (not return_indices): 531 msg = 'at least one of distances/indices must be specified' 532 raise RuntimeError, msg 533 tmp1 = numpy.asarray(input) != 0 534 struct = generate_binary_structure(tmp1.ndim, tmp1.ndim) 535 tmp2 = binary_dilation(tmp1, struct) 536 tmp2 = numpy.logical_xor(tmp1, tmp2) 537 tmp1 = tmp1.astype(numpy.int8) - tmp2.astype(numpy.int8) 538 del tmp2 539 metric = metric.lower() 540 if metric == 'euclidean': 541 metric = 1 542 elif metric == 'cityblock': 543 metric = 2 544 elif metric == 'chessboard': 545 metric = 3 546 else: 547 raise RuntimeError, 'distance metric not supported' 548 if sampling is not None: 549 sampling = _ni_support._normalize_sequence(sampling, tmp1.ndim) 550 sampling = numpy.asarray(sampling, dtype = numpy.float64) 551 if not sampling.flags.contiguous: 552 sampling = sampling.copy() 553 if return_indices: 554 ft = numpy.zeros(tmp1.shape, dtype = numpy.int32) 555 else: 556 ft = None 557 if return_distances: 558 if distances is None: 559 if metric == 1: 560 dt = numpy.zeros(tmp1.shape, dtype = numpy.float64) 561 else: 562 dt = numpy.zeros(tmp1.shape, dtype = numpy.uint32) 563 else: 564 if distances.shape != tmp1.shape: 565 raise RuntimeError, 'distances array has wrong shape' 566 if metric == 1: 567 if distances.dtype.type != numpy.float64: 568 raise RuntimeError, 'distances array must be float64' 569 else: 570 if distances.dtype.type != numpy.uint32: 571 raise RuntimeError, 'distances array must be uint32' 572 dt = distances 573 else: 574 dt = None 575 _nd_image.distance_transform_bf(tmp1, metric, sampling, dt, ft) 576 if return_indices: 577 if isinstance(indices, numpy.ndarray): 578 if indices.dtype.type != numpy.int32: 579 raise RuntimeError, 'indices must of int32 type' 580 if indices.shape != (tmp1.ndim,) + tmp1.shape: 581 raise RuntimeError, 'indices has wrong shape' 582 tmp2 = indices 583 else: 584 tmp2 = numpy.indices(tmp1.shape, dtype = numpy.int32) 585 ft = numpy.ravel(ft) 586 for ii in range(tmp2.shape[0]): 587 rtmp = numpy.ravel(tmp2[ii, ...])[ft] 588 rtmp.shape = tmp1.shape 589 tmp2[ii, ...] = rtmp 590 ft = tmp2 591 # construct and return the result 592 result = [] 593 if return_distances and not isinstance(distances, numpy.ndarray): 594 result.append(dt) 595 if return_indices and not isinstance(indices, numpy.ndarray): 596 result.append(ft) 597 if len(result) == 2: 598 return tuple(result) 599 elif len(result) == 1: 600 return result[0] 601 else: 602 return None
603
604 -def distance_transform_cdt(input, structure = 'chessboard', 605 return_distances = True, return_indices = False, 606 distances = None, indices = None):
607 """Distance transform for chamfer type of transforms. 608 609 The structure determines the type of chamfering that is done. If 610 the structure is equal to 'cityblock' a structure is generated 611 using generate_binary_structure with a squared distance equal to 612 1. If the structure is equal to 'chessboard', a structure is 613 generated using generate_binary_structure with a squared distance 614 equal to the rank of the array. These choices correspond to the 615 common interpretations of the cityblock and the chessboard 616 distance metrics in two dimensions. 617 618 In addition to the distance transform, the feature transform can 619 be calculated. In this case the index of the closest background 620 element is returned along the first axis of the result. 621 622 The return_distances, and return_indices flags can be used to 623 indicate if the distance transform, the feature transform, or both 624 must be returned. 625 626 The distances and indices arguments can be used to give optional 627 output arrays that must be of the correct size and type (both int32). 628 """ 629 if (not return_distances) and (not return_indices): 630 msg = 'at least one of distances/indices must be specified' 631 raise RuntimeError, msg 632 ft_inplace = isinstance(indices, numpy.ndarray) 633 dt_inplace = isinstance(distances, numpy.ndarray) 634 input = numpy.asarray(input) 635 if structure == 'cityblock': 636 rank = input.ndim 637 structure = generate_binary_structure(rank, 1) 638 elif structure == 'chessboard': 639 rank = input.ndim 640 structure = generate_binary_structure(rank, rank) 641 else: 642 try: 643 structure = numpy.asarray(structure) 644 except: 645 raise RuntimeError, 'invalid structure provided' 646 for s in structure.shape: 647 if s != 3: 648 raise RuntimeError, 'structure sizes must be equal to 3' 649 if not structure.flags.contiguous: 650 structure = structure.copy() 651 if dt_inplace: 652 if distances.dtype.type != numpy.int32: 653 raise RuntimeError, 'distances must be of int32 type' 654 if distances.shape != input.shape: 655 raise RuntimeError, 'distances has wrong shape' 656 dt = distances 657 dt[...] = numpy.where(input, -1, 0).astype(numpy.int32) 658 else: 659 dt = numpy.where(input, -1, 0).astype(numpy.int32) 660 rank = dt.ndim 661 if return_indices: 662 sz = numpy.product(dt.shape,axis=0) 663 ft = numpy.arange(sz, dtype = numpy.int32) 664 ft.shape = dt.shape 665 else: 666 ft = None 667 _nd_image.distance_transform_op(structure, dt, ft) 668 dt = dt[tuple([slice(None, None, -1)] * rank)] 669 if return_indices: 670 ft = ft[tuple([slice(None, None, -1)] * rank)] 671 _nd_image.distance_transform_op(structure, dt, ft) 672 dt = dt[tuple([slice(None, None, -1)] * rank)] 673 if return_indices: 674 ft = ft[tuple([slice(None, None, -1)] * rank)] 675 ft = numpy.ravel(ft) 676 if ft_inplace: 677 if indices.dtype.type != numpy.int32: 678 raise RuntimeError, 'indices must of int32 type' 679 if indices.shape != (dt.ndim,) + dt.shape: 680 raise RuntimeError, 'indices has wrong shape' 681 tmp = indices 682 else: 683 tmp = numpy.indices(dt.shape, dtype = numpy.int32) 684 for ii in range(tmp.shape[0]): 685 rtmp = numpy.ravel(tmp[ii, ...])[ft] 686 rtmp.shape = dt.shape 687 tmp[ii, ...] = rtmp 688 ft = tmp 689 690 # construct and return the result 691 result = [] 692 if return_distances and not dt_inplace: 693 result.append(dt) 694 if return_indices and not ft_inplace: 695 result.append(ft) 696 if len(result) == 2: 697 return tuple(result) 698 elif len(result) == 1: 699 return result[0] 700 else: 701 return None
702 703
704 -def distance_transform_edt(input, sampling = None, 705 return_distances = True, return_indices = False, 706 distances = None, indices = None):
707 """Exact euclidean distance transform. 708 709 In addition to the distance transform, the feature transform can 710 be calculated. In this case the index of the closest background 711 element is returned along the first axis of the result. 712 713 The return_distances, and return_indices flags can be used to 714 indicate if the distance transform, the feature transform, or both 715 must be returned. 716 717 Optionally the sampling along each axis can be given by the 718 sampling parameter which should be a sequence of length equal to 719 the input rank, or a single number in which the sampling is assumed 720 to be equal along all axes. 721 722 the distances and indices arguments can be used to give optional 723 output arrays that must be of the correct size and type (float64 724 and int32). 725 """ 726 if (not return_distances) and (not return_indices): 727 msg = 'at least one of distances/indices must be specified' 728 raise RuntimeError, msg 729 ft_inplace = isinstance(indices, numpy.ndarray) 730 dt_inplace = isinstance(distances, numpy.ndarray) 731 # calculate the feature transform 732 input = numpy.where(input, 1, 0).astype(numpy.int8) 733 if sampling is not None: 734 sampling = _ni_support._normalize_sequence(sampling, input.ndim) 735 sampling = numpy.asarray(sampling, dtype = numpy.float64) 736 if not sampling.flags.contiguous: 737 sampling = sampling.copy() 738 if ft_inplace: 739 ft = indices 740 if ft.shape != (input.ndim,) + input.shape: 741 raise RuntimeError, 'indices has wrong shape' 742 if ft.dtype.type != numpy.int32: 743 raise RuntimeError, 'indices must be of int32 type' 744 else: 745 ft = numpy.zeros((input.ndim,) + input.shape, 746 dtype = numpy.int32) 747 _nd_image.euclidean_feature_transform(input, sampling, ft) 748 # if requested, calculate the distance transform 749 if return_distances: 750 dt = ft - numpy.indices(input.shape, dtype = ft.dtype) 751 dt = dt.astype(numpy.float64) 752 if sampling is not None: 753 for ii in range(len(sampling)): 754 dt[ii, ...] *= sampling[ii] 755 numpy.multiply(dt, dt, dt) 756 if dt_inplace: 757 dt = numpy.add.reduce(dt, axis = 0) 758 if distances.shape != dt.shape: 759 raise RuntimeError, 'indices has wrong shape' 760 if distances.dtype.type != numpy.float64: 761 raise RuntimeError, 'indices must be of float64 type' 762 numpy.sqrt(dt, distances) 763 del dt 764 else: 765 dt = numpy.add.reduce(dt, axis = 0) 766 dt = numpy.sqrt(dt) 767 # construct and return the result 768 result = [] 769 if return_distances and not dt_inplace: 770 result.append(dt) 771 if return_indices and not ft_inplace: 772 result.append(ft) 773 if len(result) == 2: 774 return tuple(result) 775 elif len(result) == 1: 776 return result[0] 777 else: 778 return None
779