~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/ndimage/morphology.py

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

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