← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/ndimage/interpolation.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
http://bugs.debian.org/489149
* 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:
Lib/ndimage/interpolation.py
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 types
32
 
import math
33
 
import numpy
34
 
import _ni_support
35
 
import _nd_image
36
 
 
37
 
 
38
 
def spline_filter1d(input, order = 3, axis = -1, output = numpy.float64,
39
 
                    output_type = None):
40
 
    """Calculates a one-dimensional spline filter along the given axis.
41
 
 
42
 
    The lines of the array along the given axis are filtered by a
43
 
    spline filter. The order of the spline must be >= 2 and <= 5.
44
 
    """
45
 
    if order < 0 or order > 5:
46
 
        raise RuntimeError, 'spline order not supported'
47
 
    input = numpy.asarray(input)
48
 
    if numpy.iscomplexobj(input):
49
 
        raise TypeError, 'Complex type not supported'
50
 
    output, return_value = _ni_support._get_output(output, input,
51
 
                                                    output_type)
52
 
    if order in [0, 1]:
53
 
        output[...] = numpy.array(input)
54
 
    else:
55
 
        axis = _ni_support._check_axis(axis, input.ndim)
56
 
        _nd_image.spline_filter1d(input, order, axis, output)
57
 
    return return_value
58
 
 
59
 
 
60
 
def spline_filter(input, order = 3, output = numpy.float64,
61
 
                  output_type = None):
62
 
    """Multi-dimensional spline filter.
63
 
 
64
 
    Note: The multi-dimensional filter is implemented as a sequence of
65
 
    one-dimensional spline filters. The intermediate arrays are stored
66
 
    in the same data type as the output. Therefore, for output types
67
 
    with a limited precision, the results may be imprecise because
68
 
    intermediate results may be stored with insufficient precision.
69
 
    """
70
 
    if order < 2 or order > 5:
71
 
        raise RuntimeError, 'spline order not supported'
72
 
    input = numpy.asarray(input)
73
 
    if numpy.iscomplexobj(input):
74
 
        raise TypeError, 'Complex type not supported'
75
 
    output, return_value = _ni_support._get_output(output, input,
76
 
                                                    output_type)
77
 
    if order not in [0, 1] and input.ndim > 0:
78
 
        for axis in range(input.ndim):
79
 
            spline_filter1d(input, order, axis, output = output)
80
 
            input = output
81
 
    else:
82
 
        output[...] = input[...]
83
 
    return return_value
84
 
 
85
 
def geometric_transform(input, mapping, output_shape = None,
86
 
                        output_type = None, output = None, order = 3,
87
 
                        mode = 'constant', cval = 0.0, prefilter = True,
88
 
                        extra_arguments = (), extra_keywords = {}):
89
 
    """Apply an arbritrary geometric transform.
90
 
 
91
 
    The given mapping function is used to find, for each point in the
92
 
    output, the corresponding coordinates in the input. The value of the
93
 
    input at those coordinates is determined by spline interpolation of
94
 
    the requested order.
95
 
 
96
 
    mapping must be a callable object that accepts a tuple of length
97
 
    equal to the output array rank and returns the corresponding input
98
 
    coordinates as a tuple of length equal to the input array
99
 
    rank. Points outside the boundaries of the input are filled
100
 
    according to the given mode ('constant', 'nearest', 'reflect' or
101
 
    'wrap'). The output shape can optionally be given. If not given,
102
 
    it is equal to the input shape. The parameter prefilter determines
103
 
    if the input is pre-filtered before interpolation (necessary for
104
 
    spline interpolation of order > 1).  If False it is assumed that
105
 
    the input is already filtered. The extra_arguments and
106
 
    extra_keywords arguments can be used to provide extra arguments
107
 
    and keywords that are passed to the mapping function at each call.
108
 
 
109
 
    Example usage:
110
 
      >>> a = arange(12.).reshape((4,3))
111
 
      >>> def shift_func(output_coordinates):
112
 
      ...     return (output_coordinates[0]-0.5, output_coordinates[1]-0.5)
113
 
      ...
114
 
      >>> print geometric_transform(a,shift_func)
115
 
      array([[ 0.    ,  0.    ,  0.    ],
116
 
             [ 0.    ,  1.3625,  2.7375],
117
 
             [ 0.    ,  4.8125,  6.1875],
118
 
             [ 0.    ,  8.2625,  9.6375]])
119
 
    """
120
 
    if order < 0 or order > 5:
121
 
        raise RuntimeError, 'spline order not supported'
122
 
    input = numpy.asarray(input)
123
 
    if numpy.iscomplexobj(input):
124
 
        raise TypeError, 'Complex type not supported'
125
 
    if output_shape == None:
126
 
        output_shape = input.shape
127
 
    if input.ndim < 1 or len(output_shape) < 1:
128
 
        raise RuntimeError, 'input and output rank must be > 0'
129
 
    mode = _ni_support._extend_mode_to_code(mode)
130
 
    if prefilter and order > 1:
131
 
        filtered = spline_filter(input, order, output = numpy.float64)
132
 
    else:
133
 
        filtered = input
134
 
    output, return_value = _ni_support._get_output(output, input,
135
 
                                        output_type, shape = output_shape)
136
 
    _nd_image.geometric_transform(filtered, mapping, None, None, None,
137
 
               output, order, mode, cval, extra_arguments, extra_keywords)
138
 
    return return_value
139
 
 
140
 
 
141
 
def map_coordinates(input, coordinates, output_type = None, output = None,
142
 
                order = 3, mode = 'constant', cval = 0.0, prefilter = True):
143
 
    """Apply an arbritrary coordinate transformation.
144
 
 
145
 
    The array of coordinates is used to find, for each point in the output,
146
 
    the corresponding coordinates in the input. The value of the input at
147
 
    that coordinates is determined by spline interpolation of the
148
 
    requested order.
149
 
 
150
 
    The shape of the output is derived from that of the coordinate
151
 
    array by dropping the first axis. The values of the array along
152
 
    the first axis are the coordinates in the input array at which the
153
 
    output value is found.  For example, if the input has dimensions
154
 
    (100,200,3), then the shape of coordinates will be (3,100,200,3),
155
 
    where coordinates[:,1,2,3] specify the input coordinate at which
156
 
    output[1,2,3] is found.
157
 
 
158
 
    Points outside the boundaries of the input are filled according to
159
 
    the given mode ('constant', 'nearest', 'reflect' or 'wrap'). The
160
 
    parameter prefilter determines if the input is pre-filtered before
161
 
    interpolation (necessary for spline interpolation of order >
162
 
    1). If False it is assumed that the input is already filtered.
163
 
 
164
 
    Example usage:
165
 
      >>> a = arange(12.).reshape((4,3))
166
 
      >>> print a
167
 
      [[  0.   1.   2.]
168
 
       [  3.   4.   5.]
169
 
       [  6.   7.   8.]
170
 
       [  9.  10.  11.]]
171
 
      >>> output = map_coordinates(a,[[0.5, 2], [0.5, 1]],order=1)
172
 
      >>> print output
173
 
      [ 2. 7.]
174
 
 
175
 
      Here, the interpolated value of a[0.5,0.5] gives output[0], while
176
 
      a[2,1] is output[1].
177
 
    """
178
 
    if order < 0 or order > 5:
179
 
        raise RuntimeError, 'spline order not supported'
180
 
    input = numpy.asarray(input)
181
 
    if numpy.iscomplexobj(input):
182
 
        raise TypeError, 'Complex type not supported'
183
 
    coordinates = numpy.asarray(coordinates)
184
 
    if numpy.iscomplexobj(coordinates):
185
 
        raise TypeError, 'Complex type not supported'
186
 
    output_shape = coordinates.shape[1:]
187
 
    if input.ndim < 1 or len(output_shape) < 1:
188
 
        raise RuntimeError, 'input and output rank must be > 0'
189
 
    if coordinates.shape[0] != input.ndim:
190
 
        raise RuntimeError, 'invalid shape for coordinate array'
191
 
    mode = _ni_support._extend_mode_to_code(mode)
192
 
    if prefilter and order > 1:
193
 
        filtered = spline_filter(input, order, output = numpy.float64)
194
 
    else:
195
 
        filtered = input
196
 
    output, return_value = _ni_support._get_output(output, input,
197
 
                                        output_type, shape = output_shape)
198
 
    _nd_image.geometric_transform(filtered, None, coordinates, None, None,
199
 
               output, order, mode, cval, None, None)
200
 
    return return_value
201
 
 
202
 
 
203
 
def affine_transform(input, matrix, offset = 0.0, output_shape = None,
204
 
                     output_type = None, output = None, order = 3,
205
 
                     mode = 'constant', cval = 0.0, prefilter = True):
206
 
    """Apply an affine transformation.
207
 
 
208
 
    The given matrix and offset are used to find for each point in the
209
 
    output the corresponding coordinates in the input by an affine
210
 
    transformation. The value of the input at those coordinates is
211
 
    determined by spline interpolation of the requested order. Points
212
 
    outside the boundaries of the input are filled according to the given
213
 
    mode. The output shape can optionally be given. If not given it is
214
 
    equal to the input shape. The parameter prefilter determines if the
215
 
    input is pre-filtered before interpolation, if False it is assumed
216
 
    that the input is already filtered.
217
 
 
218
 
    The matrix must be two-dimensional or can also be given as a
219
 
    one-dimensional sequence or array. In the latter case, it is
220
 
    assumed that the matrix is diagonal. A more efficient algorithms
221
 
    is then applied that exploits the separability of the problem.
222
 
    """
223
 
    if order < 0 or order > 5:
224
 
        raise RuntimeError, 'spline order not supported'
225
 
    input = numpy.asarray(input)
226
 
    if numpy.iscomplexobj(input):
227
 
        raise TypeError, 'Complex type not supported'
228
 
    if output_shape == None:
229
 
        output_shape = input.shape
230
 
    if input.ndim < 1 or len(output_shape) < 1:
231
 
        raise RuntimeError, 'input and output rank must be > 0'
232
 
    mode = _ni_support._extend_mode_to_code(mode)
233
 
    if prefilter and order > 1:
234
 
        filtered = spline_filter(input, order, output = numpy.float64)
235
 
    else:
236
 
        filtered = input
237
 
    output, return_value = _ni_support._get_output(output, input,
238
 
                                        output_type, shape = output_shape)
239
 
    matrix = numpy.asarray(matrix, dtype = numpy.float64)
240
 
    if matrix.ndim not in [1, 2] or matrix.shape[0] < 1:
241
 
        raise RuntimeError, 'no proper affine matrix provided'
242
 
    if matrix.shape[0] != input.ndim:
243
 
        raise RuntimeError, 'affine matrix has wrong number of rows'
244
 
    if matrix.ndim == 2 and matrix.shape[1] != output.ndim:
245
 
        raise RuntimeError, 'affine matrix has wrong number of columns'
246
 
    if not matrix.flags.contiguous:
247
 
        matrix = matrix.copy()
248
 
    offset = _ni_support._normalize_sequence(offset, input.ndim)
249
 
    offset = numpy.asarray(offset, dtype = numpy.float64)
250
 
    if offset.ndim != 1 or offset.shape[0] < 1:
251
 
        raise RuntimeError, 'no proper offset provided'
252
 
    if not offset.flags.contiguous:
253
 
        offset = offset.copy()
254
 
    if matrix.ndim == 1:
255
 
        _nd_image.zoom_shift(filtered, matrix, offset, output, order,
256
 
                             mode, cval)
257
 
    else:
258
 
        _nd_image.geometric_transform(filtered, None, None, matrix, offset,
259
 
                            output, order, mode, cval, None, None)
260
 
    return return_value
261
 
 
262
 
 
263
 
def shift(input, shift, output_type = None, output = None, order = 3,
264
 
          mode = 'constant', cval = 0.0, prefilter = True):
265
 
    """Shift an array.
266
 
 
267
 
    The array is shifted using spline interpolation of the requested
268
 
    order. Points outside the boundaries of the input are filled according
269
 
    to the given mode. The parameter prefilter determines if the input is
270
 
    pre-filtered before interpolation, if False it is assumed that the
271
 
    input is already filtered.
272
 
    """
273
 
    if order < 0 or order > 5:
274
 
        raise RuntimeError, 'spline order not supported'
275
 
    input = numpy.asarray(input)
276
 
    if numpy.iscomplexobj(input):
277
 
        raise TypeError, 'Complex type not supported'
278
 
    if input.ndim < 1:
279
 
        raise RuntimeError, 'input and output rank must be > 0'
280
 
    mode = _ni_support._extend_mode_to_code(mode)
281
 
    if prefilter and order > 1:
282
 
        filtered = spline_filter(input, order, output = numpy.float64)
283
 
    else:
284
 
        filtered = input
285
 
    output, return_value = _ni_support._get_output(output, input,
286
 
                                                    output_type)
287
 
    shift = _ni_support._normalize_sequence(shift, input.ndim)
288
 
    shift = [-ii for ii in shift]
289
 
    shift = numpy.asarray(shift, dtype = numpy.float64)
290
 
    if not shift.flags.contiguous:
291
 
        shift = shift.copy()
292
 
    _nd_image.zoom_shift(filtered, None, shift, output, order, mode, cval)
293
 
    return return_value
294
 
 
295
 
 
296
 
def zoom(input, zoom, output_type = None, output = None, order = 3,
297
 
         mode = 'constant', cval = 0.0, prefilter = True):
298
 
    """Zoom an array.
299
 
 
300
 
    The array is zoomed using spline interpolation of the requested order.
301
 
    Points outside the boundaries of the input are filled according to the
302
 
    given mode. The parameter prefilter determines if the input is pre-
303
 
    filtered before interpolation, if False it is assumed that the input
304
 
    is already filtered.
305
 
    """
306
 
    if order < 0 or order > 5:
307
 
        raise RuntimeError, 'spline order not supported'
308
 
    input = numpy.asarray(input)
309
 
    if numpy.iscomplexobj(input):
310
 
        raise TypeError, 'Complex type not supported'
311
 
    if input.ndim < 1:
312
 
        raise RuntimeError, 'input and output rank must be > 0'
313
 
    mode = _ni_support._extend_mode_to_code(mode)
314
 
    if prefilter and order > 1:
315
 
        filtered = spline_filter(input, order, output = numpy.float64)
316
 
    else:
317
 
        filtered = input
318
 
    zoom = _ni_support._normalize_sequence(zoom, input.ndim)
319
 
    output_shape = [int(ii * jj) for ii, jj in zip(input.shape, zoom)]
320
 
    zoom = [1.0 / ii for ii in zoom]
321
 
    output, return_value = _ni_support._get_output(output, input,
322
 
                                        output_type, shape = output_shape)
323
 
    zoom = numpy.asarray(zoom, dtype = numpy.float64)
324
 
    if not zoom.flags.contiguous:
325
 
        zoom = shift.copy()
326
 
    _nd_image.zoom_shift(filtered, zoom, None, output, order, mode, cval)
327
 
    return return_value
328
 
 
329
 
def _minmax(coor, minc, maxc):
330
 
    if coor[0] < minc[0]:
331
 
        minc[0] = coor[0]
332
 
    if coor[0] > maxc[0]:
333
 
        maxc[0] = coor[0]
334
 
    if coor[1] < minc[1]:
335
 
        minc[1] = coor[1]
336
 
    if coor[1] > maxc[1]:
337
 
        maxc[1] = coor[1]
338
 
    return minc, maxc
339
 
 
340
 
def rotate(input, angle, axes = (-1, -2), reshape = True,
341
 
           output_type = None, output = None, order = 3,
342
 
           mode = 'constant', cval = 0.0, prefilter = True):
343
 
    """Rotate an array.
344
 
 
345
 
    The array is rotated in the plane defined by the two axes given by the
346
 
    axes parameter using spline interpolation of the requested order. The
347
 
    angle is given in degrees. Points outside the boundaries of the input
348
 
    are filled according to the given mode. If reshape is true, the output
349
 
    shape is adapted so that the input array is contained completely in
350
 
    the output. The parameter prefilter determines if the input is pre-
351
 
    filtered before interpolation, if False it is assumed that the input
352
 
    is already filtered.
353
 
    """
354
 
    input = numpy.asarray(input)
355
 
    axes = list(axes)
356
 
    rank = input.ndim
357
 
    if axes[0] < 0:
358
 
        axes[0] += rank
359
 
    if axes[1] < 0:
360
 
        axes[1] += rank
361
 
    if axes[0] < 0 or axes[1] < 0 or axes[0] > rank or axes[1] > rank:
362
 
        raise RuntimeError, 'invalid rotation plane specified'
363
 
    if axes[0] > axes[1]:
364
 
        axes = axes[1], axes[0]
365
 
    angle = numpy.pi / 180 * angle
366
 
    m11 = math.cos(angle)
367
 
    m12 = math.sin(angle)
368
 
    m21 = -math.sin(angle)
369
 
    m22 = math.cos(angle)
370
 
    matrix = numpy.array([[m11, m12],
371
 
                             [m21, m22]], dtype = numpy.float64)
372
 
    iy = input.shape[axes[0]]
373
 
    ix = input.shape[axes[1]]
374
 
    if reshape:
375
 
        mtrx = numpy.array([[ m11, -m21],
376
 
                               [-m12,  m22]], dtype = numpy.float64)
377
 
        minc = [0, 0]
378
 
        maxc = [0, 0]
379
 
        coor = numpy.dot(mtrx, [0, ix])
380
 
        minc, maxc = _minmax(coor, minc, maxc)
381
 
        coor = numpy.dot(mtrx, [iy, 0])
382
 
        minc, maxc = _minmax(coor, minc, maxc)
383
 
        coor = numpy.dot(mtrx, [iy, ix])
384
 
        minc, maxc = _minmax(coor, minc, maxc)
385
 
        oy = int(maxc[0] - minc[0] + 0.5)
386
 
        ox = int(maxc[1] - minc[1] + 0.5)
387
 
    else:
388
 
        oy = input.shape[axes[0]]
389
 
        ox = input.shape[axes[1]]
390
 
    offset = numpy.zeros((2,), dtype = numpy.float64)
391
 
    offset[0] = float(oy) / 2.0 - 0.5
392
 
    offset[1] = float(ox) / 2.0 - 0.5
393
 
    offset = numpy.dot(matrix, offset)
394
 
    tmp = numpy.zeros((2,), dtype = numpy.float64)
395
 
    tmp[0] = float(iy) / 2.0 - 0.5
396
 
    tmp[1] = float(ix) / 2.0 - 0.5
397
 
    offset = tmp - offset
398
 
    output_shape = list(input.shape)
399
 
    output_shape[axes[0]] = oy
400
 
    output_shape[axes[1]] = ox
401
 
    output_shape = tuple(output_shape)
402
 
    output, return_value = _ni_support._get_output(output, input,
403
 
                                        output_type, shape = output_shape)
404
 
    if input.ndim <= 2:
405
 
        affine_transform(input, matrix, offset, output_shape, None, output,
406
 
                         order, mode, cval, prefilter)
407
 
    else:
408
 
        coordinates = []
409
 
        size = numpy.product(input.shape,axis=0)
410
 
        size /= input.shape[axes[0]]
411
 
        size /= input.shape[axes[1]]
412
 
        for ii in range(input.ndim):
413
 
            if ii not in axes:
414
 
                coordinates.append(0)
415
 
            else:
416
 
                coordinates.append(slice(None, None, None))
417
 
        iter_axes = range(input.ndim)
418
 
        iter_axes.reverse()
419
 
        iter_axes.remove(axes[0])
420
 
        iter_axes.remove(axes[1])
421
 
        os = (output_shape[axes[0]], output_shape[axes[1]])
422
 
        for ii in range(size):
423
 
            ia = input[tuple(coordinates)]
424
 
            oa = output[tuple(coordinates)]
425
 
            affine_transform(ia, matrix, offset, os, None, oa, order, mode,
426
 
                             cval, prefilter)
427
 
            for jj in iter_axes:
428
 
                if coordinates[jj] < input.shape[jj] - 1:
429
 
                    coordinates[jj] += 1
430
 
                    break
431
 
                else:
432
 
                    coordinates[jj] = 0
433
 
    return return_value