Skip to content

Commit

Permalink
BUG: ndimage: Correctly handle the intp type on 32-bit platforms.
Browse files Browse the repository at this point in the history
  • Loading branch information
@stefanv
stefanv committed Aug 29, 2011
2 parents 936e795 + 2bbcb3f commit 7e7ac99
Show file tree
Hide file tree
Showing 3 changed files with 83 additions and 5 deletions.
2 changes: 2 additions & 0 deletions scipy/ndimage/src/nd_image.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,8 @@ typedef enum
tInt32=PyArray_INT32,
tUInt32=PyArray_UINT32,
tInt64=PyArray_INT64,
tuintp=PyArray_UINTP,
tintp=PyArray_INTP,
tUInt64=PyArray_UINT64,
tFloat32=PyArray_FLOAT32,
tFloat64=PyArray_FLOAT64,
Expand Down
27 changes: 22 additions & 5 deletions scipy/ndimage/src/ni_interpolation.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -307,10 +307,11 @@ int NI_SplineFilter1D(PyArrayObject *input, int order, int axis,
return PyErr_Occurred() ? 0 : 1;
}

/* copy row of coordinate array from location at _p to _coor */
#define CASE_MAP_COORDINATES(_p, _coor, _rank, _stride, _type) \
case t ## _type: \
{ \
int _hh; \
npy_intp _hh; \
for(_hh = 0; _hh < _rank; _hh++) { \
_coor[_hh] = *(_type*)_p; \
_p += _stride; \
Expand Down Expand Up @@ -487,13 +488,15 @@ NI_GeometricTransform(PyArrayObject *input, int (*map)(npy_intp*, double*,
CASE_MAP_COORDINATES(p, icoor, irank, cstride, UInt8);
CASE_MAP_COORDINATES(p, icoor, irank, cstride, UInt16);
CASE_MAP_COORDINATES(p, icoor, irank, cstride, UInt32);
CASE_MAP_COORDINATES(p, icoor, irank, cstride, uintp);
#if HAS_UINT64
CASE_MAP_COORDINATES(p, icoor, irank, cstride, UInt64);
#endif
CASE_MAP_COORDINATES(p, icoor, irank, cstride, Int8);
CASE_MAP_COORDINATES(p, icoor, irank, cstride, Int16);
CASE_MAP_COORDINATES(p, icoor, irank, cstride, Int32);
CASE_MAP_COORDINATES(p, icoor, irank, cstride, Int64);
CASE_MAP_COORDINATES(p, icoor, irank, cstride, intp);
CASE_MAP_COORDINATES(p, icoor, irank, cstride, Float32);
CASE_MAP_COORDINATES(p, icoor, irank, cstride, Float64);
default:
Expand Down Expand Up @@ -580,13 +583,15 @@ NI_GeometricTransform(PyArrayObject *input, int (*map)(npy_intp*, double*,
CASE_INTERP_COEFF(coeff, pi, idxs[hh], UInt8);
CASE_INTERP_COEFF(coeff, pi, idxs[hh], UInt16);
CASE_INTERP_COEFF(coeff, pi, idxs[hh], UInt32);
CASE_INTERP_COEFF(coeff, pi, idxs[hh], uintp);
#if HAS_UINT64
CASE_INTERP_COEFF(coeff, pi, idxs[hh], UInt64);
#endif
CASE_INTERP_COEFF(coeff, pi, idxs[hh], Int8);
CASE_INTERP_COEFF(coeff, pi, idxs[hh], Int16);
CASE_INTERP_COEFF(coeff, pi, idxs[hh], Int32);
CASE_INTERP_COEFF(coeff, pi, idxs[hh], Int64);
CASE_INTERP_COEFF(coeff, pi, idxs[hh], intp);
CASE_INTERP_COEFF(coeff, pi, idxs[hh], Float32);
CASE_INTERP_COEFF(coeff, pi, idxs[hh], Float64);
default:
Expand All @@ -610,14 +615,19 @@ NI_GeometricTransform(PyArrayObject *input, int (*map)(npy_intp*, double*,
CASE_INTERP_OUT_UINT(po, t, UInt8, 0, MAX_UINT8);
CASE_INTERP_OUT_UINT(po, t, UInt16, 0, MAX_UINT16);
CASE_INTERP_OUT_UINT(po, t, UInt32, 0, MAX_UINT32);
CASE_INTERP_OUT_UINT(po, t, uintp, 0, MAX_UINTP);
#if HAS_UINT64
/* FIXME */
CASE_INTERP_OUT_UINT(po, t, UInt64, 0, MAX_UINT32);
/* There was a bug in numpy as of (at least) <= 1.6.1 such that
* MAX_UINT64 was incorrectly defined, leading to a compiler error.
* NPY_MAX_UINT64 is correctly defined
*/
CASE_INTERP_OUT_UINT(po, t, UInt64, 0, NPY_MAX_UINT64);
#endif
CASE_INTERP_OUT_INT(po, t, Int8, MIN_INT8, MAX_INT8);
CASE_INTERP_OUT_INT(po, t, Int16, MIN_INT16, MAX_INT16);
CASE_INTERP_OUT_INT(po, t, Int32, MIN_INT32, MAX_INT32);
CASE_INTERP_OUT_INT(po, t, Int64, MIN_INT64, MAX_INT64);
CASE_INTERP_OUT_INT(po, t, intp, MIN_INTP, MAX_INTP);
CASE_INTERP_OUT(po, t, Float32);
CASE_INTERP_OUT(po, t, Float64);
default:
Expand Down Expand Up @@ -876,13 +886,15 @@ int NI_ZoomShift(PyArrayObject *input, PyArrayObject* zoom_ar,
CASE_INTERP_COEFF(coeff, pi, idxs[hh], UInt8);
CASE_INTERP_COEFF(coeff, pi, idxs[hh], UInt16);
CASE_INTERP_COEFF(coeff, pi, idxs[hh], UInt32);
CASE_INTERP_COEFF(coeff, pi, idxs[hh], uintp);
#if HAS_UINT64
CASE_INTERP_COEFF(coeff, pi, idxs[hh], UInt64);
#endif
CASE_INTERP_COEFF(coeff, pi, idxs[hh], Int8);
CASE_INTERP_COEFF(coeff, pi, idxs[hh], Int16);
CASE_INTERP_COEFF(coeff, pi, idxs[hh], Int32);
CASE_INTERP_COEFF(coeff, pi, idxs[hh], Int64);
CASE_INTERP_COEFF(coeff, pi, idxs[hh], intp);
CASE_INTERP_COEFF(coeff, pi, idxs[hh], Float32);
CASE_INTERP_COEFF(coeff, pi, idxs[hh], Float64);
default:
Expand All @@ -906,14 +918,19 @@ int NI_ZoomShift(PyArrayObject *input, PyArrayObject* zoom_ar,
CASE_INTERP_OUT_UINT(po, t, UInt8, 0, MAX_UINT8);
CASE_INTERP_OUT_UINT(po, t, UInt16, 0, MAX_UINT16);
CASE_INTERP_OUT_UINT(po, t, UInt32, 0, MAX_UINT32);
CASE_INTERP_OUT_UINT(po, t, uintp, 0, MAX_UINTP);
#if HAS_UINT64
/* FIXME */
CASE_INTERP_OUT_UINT(po, t, UInt64, 0, MAX_UINT32);
/* There was a bug in numpy as of (at least) <= 1.6.1 such that
* MAX_UINT64 was incorrectly defined, leading to a compiler error.
* NPY_MAX_UINT64 is correctly defined
*/
CASE_INTERP_OUT_UINT(po, t, UInt64, 0, NPY_MAX_UINT64);
#endif
CASE_INTERP_OUT_INT(po, t, Int8, MIN_INT8, MAX_INT8);
CASE_INTERP_OUT_INT(po, t, Int16, MIN_INT16, MAX_INT16);
CASE_INTERP_OUT_INT(po, t, Int32, MIN_INT32, MAX_INT32);
CASE_INTERP_OUT_INT(po, t, Int64, MIN_INT64, MAX_INT64);
CASE_INTERP_OUT_INT(po, t, intp, MIN_INTP, MAX_INTP);
CASE_INTERP_OUT(po, t, Float32);
CASE_INTERP_OUT(po, t, Float64);
default:
Expand Down
59 changes: 59 additions & 0 deletions scipy/ndimage/tests/test_datatypes.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,59 @@
""" Testing data types for ndimage calls
"""

import numpy as np

from scipy import ndimage

from numpy.testing import (assert_array_almost_equal,
assert_array_equal)

from nose.tools import assert_true, assert_equal, assert_raises


def test_map_coordinates_dts():
# check that ndimage accepts different data types for interpolation
data = np.array([[4, 1, 3, 2],
[7, 6, 8, 5],
[3, 5, 3, 6]])
shifted_data = np.array([[0, 0, 0, 0],
[0, 4, 1, 3],
[0, 7, 6, 8]])
idx = np.indices(data.shape)
dts = (np.uint8, np.uint16, np.uint32, np.uint64, np.int8, np.int16,
np.int32, np.intp, np.int64, np.float32, np.float64)
for order in range(0, 6):
for data_dt in dts:
these_data = data.astype(data_dt)
for coord_dt in dts:
# affine mapping
mat = np.eye(2, dtype=coord_dt)
off = np.zeros((2,), dtype=coord_dt)
out = ndimage.affine_transform(these_data, mat, off)
assert_array_almost_equal(these_data, out)
# map coordinates
coords_m1 = idx.astype(coord_dt) - 1
coords_p10 = idx.astype(coord_dt) + 10
out = ndimage.map_coordinates(these_data, coords_m1, order=order)
assert_array_almost_equal(out, shifted_data)
# check constant fill works
out = ndimage.map_coordinates(these_data, coords_p10, order=order)
assert_array_almost_equal(out, np.zeros((3,4)))
# check shift and zoom
out = ndimage.shift(these_data, 1)
assert_array_almost_equal(out, shifted_data)
out = ndimage.zoom(these_data, 1)
assert_array_almost_equal(these_data, out)


def test_uint64_max():
# Test interpolation respects uint64 max
big = 2**64-1
arr = np.array([big, big, big], dtype=np.uint64)
# Tests geometric transform (map_coordinates, affine_transform)
inds = np.indices(arr.shape) - 0.1
x = ndimage.map_coordinates(arr, inds)
assert_true(x[1] > (2**63))
# Tests zoom / shift
x = ndimage.shift(arr, 0.1)
assert_true(x[1] > (2**63))

0 comments on commit 7e7ac99

Please sign in to comment.