Python scipy.ndimage.affine_transform() Examples

def map_image(img, out_affine, out_shape, ras2ras=np.array([[1.0, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]),
              order=1):
    """
    Function to map image to new voxel space (RAS orientation)

    :param nibabel.MGHImage img: the src 3D image with data and affine set
    :param np.ndarray out_affine: trg image affine
    :param np.ndarray out_shape: the trg shape information
    :param np.ndarray ras2ras: ras2ras an additional maping that should be applied (default=id to just reslice)
    :param int order: order of interpolation (0=nearest,1=linear(default),2=quadratic,3=cubic)
    :return: mapped Image data array
    """
    from scipy.ndimage import affine_transform
    from numpy.linalg import inv

    # compute vox2vox from src to trg
    vox2vox = inv(out_affine) @ ras2ras @ img.affine

    # here we apply the inverse vox2vox (to pull back the src info to the target image)
    new_data = affine_transform(img.get_data(), inv(vox2vox), output_shape=out_shape, order=order)
    return new_data 

def transform_image(image, angle=0.0, scale=1.0, aniso=1.0, translation=(0, 0), order=1):
    dx, dy = translation
    scale = 1.0/scale
    c = cos(angle)
    s = sin(angle)
    sm = np.array([[scale / aniso, 0], [0, scale * aniso]], 'f')
    m = np.array([[c, -s], [s, c]], 'f')
    m = np.dot(sm, m)
    w, h = image.shape
    c = np.array([w, h]) / 2.0
    d = c - np.dot(m, c) + np.array([dx * w, dy * h])
    return ndi.affine_transform(image, m, offset=d, order=order, mode="nearest", output=dtype("f")) 

def transform_coordinate_space(modality_1, modality_2):
    """
    Accepts nifty objects
    Transfers coordinate space from modality_2 to modality_1
    """
    aff_t1 = modality_1.affine
    aff_t2 = modality_2.affine
    inv_af_2 = np.linalg.inv(aff_t2)

    out_shape = modality_1.get_fdata().shape

    # desired transformation
    T = inv_af_2.dot(aff_t1)
    transformed = ndimage.affine_transform(modality_2.get_fdata(), T, output_shape=out_shape)

    return transformed 

def test_affine_transform26(self):
        # test homogeneous coordinates
        data = numpy.array([[4, 1, 3, 2],
                            [7, 6, 8, 5],
                            [3, 5, 3, 6]])
        for order in range(0, 6):
            if (order > 1):
                filtered = ndimage.spline_filter(data, order=order)
            else:
                filtered = data
            tform_original = numpy.eye(2)
            offset_original = -numpy.ones((2, 1))
            tform_h1 = numpy.hstack((tform_original, offset_original))
            tform_h2 = numpy.vstack((tform_h1, [[0, 0, 1]]))
            out1 = ndimage.affine_transform(filtered, tform_original,
                                            offset_original.ravel(),
                                            order=order, prefilter=False)
            out2 = ndimage.affine_transform(filtered, tform_h1, order=order,
                                            prefilter=False)
            out3 = ndimage.affine_transform(filtered, tform_h2, order=order,
                                            prefilter=False)
            for out in [out1, out2, out3]:
                assert_array_almost_equal(out, [[0, 0, 0, 0],
                                                [0, 4, 1, 3],
                                                [0, 7, 6, 8]]) 

def test_uint64_max():
    # Test interpolation respects uint64 max.  Reported to fail at least on
    # win32 (due to the 32 bit visual C compiler using signed int64 when
    # converting between uint64 to double) and Debian on s390x.
    # Interpolation is always done in double precision floating point, so we
    # use the largest uint64 value for which int(float(big)) still fits in
    # a uint64.
    big = 2**64-1025
    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_equal(x[1], int(float(big)))
    assert_equal(x[2], int(float(big)))
    # Tests zoom / shift
    x = ndimage.shift(arr, 0.1)
    assert_equal(x[1], int(float(big)))
    assert_equal(x[2], int(float(big))) 

def test_affine_transform09(self):
        data = numpy.array([[4, 1, 3, 2],
                               [7, 6, 8, 5],
                               [3, 5, 3, 6]])
        for order in range(0, 6):
            if (order > 1):
                filtered = ndimage.spline_filter(data,
                                                           order=order)
            else:
                filtered = data
            out = ndimage.affine_transform(filtered,[[1, 0],
                                                               [0, 1]],
                                  [-1, -1], order=order, prefilter=False)
            assert_array_almost_equal(out, [[0, 0, 0, 0],
                                       [0, 4, 1, 3],
                                       [0, 7, 6, 8]]) 

def affine_transform(image, affine_value):
    shape = image.shape
    alpha_affine = min(shape[0], shape[1]) * affine_value
    random_state = np.random.RandomState(None)
    # Random affine
    shape_size = shape[:2]
    center_square = np.float32(shape_size) // 2
    square_size = min(shape_size) // 3
    pts1 = np.float32(
        [center_square + square_size, [center_square[0] + square_size, center_square[1] - square_size],
         center_square - square_size])
    pts2 = pts1 + random_state.uniform(-alpha_affine, alpha_affine, size=pts1.shape).astype(np.float32)
    M = calcAffineMatrix(pts1, pts2)
    R = M[0:2, 0:2]
    Off = M[:, 2]
    for aD in range(shape[2]):
        image[:, :, aD] = scipy_affine_transform(image[:, :, aD], R, offset=Off)
    return image 

def test_affine_transform09(self):
        data = numpy.array([[4, 1, 3, 2],
                            [7, 6, 8, 5],
                            [3, 5, 3, 6]])
        for order in range(0, 6):
            if (order > 1):
                filtered = ndimage.spline_filter(data, order=order)
            else:
                filtered = data
            out = ndimage.affine_transform(filtered, [[1, 0], [0, 1]],
                                           [-1, -1], order=order,
                                           prefilter=False)
            assert_array_almost_equal(out, [[0, 0, 0, 0],
                                            [0, 4, 1, 3],
                                            [0, 7, 6, 8]]) 

def test_affine_transform08(self):
        data = numpy.array([[4, 1, 3, 2],
                            [7, 6, 8, 5],
                            [3, 5, 3, 6]])
        for order in range(0, 6):
            out = ndimage.affine_transform(data, [[1, 0], [0, 1]],
                                           [-1, -1], order=order)
            assert_array_almost_equal(out, [[0, 0, 0, 0],
                                            [0, 4, 1, 3],
                                            [0, 7, 6, 8]]) 

def test_affine_transform11(self):
        data = [1, 5, 2, 6, 3, 7, 4, 4]
        for order in range(0, 6):
            out = ndimage.affine_transform(data, [[2]], 0, (4,), order=order)
            assert_array_almost_equal(out, [1, 2, 3, 4]) 

def test_affine_transform12(self):
        data = [1, 2, 3, 4]
        for order in range(0, 6):
            out = ndimage.affine_transform(data, [[0.5]], 0, (8,), order=order)
            assert_array_almost_equal(out[::2], [1, 2, 3, 4]) 

def test_affine_transform13(self):
        data = [[1, 2, 3, 4],
                [5, 6, 7, 8],
                [9.0, 10, 11, 12]]
        for order in range(0, 6):
            out = ndimage.affine_transform(data, [[1, 0], [0, 2]], 0, (3, 2),
                                           order=order)
            assert_array_almost_equal(out, [[1, 3], [5, 7], [9, 11]]) 

def test_affine_transform14(self):
        data = [[1, 2, 3, 4],
                [5, 6, 7, 8],
                [9, 10, 11, 12]]
        for order in range(0, 6):
            out = ndimage.affine_transform(data, [[2, 0], [0, 1]], 0, (1, 4),
                                           order=order)
            assert_array_almost_equal(out, [[1, 2, 3, 4]]) 

def test_affine_transform15(self):
        data = [[1, 2, 3, 4],
                [5, 6, 7, 8],
                [9, 10, 11, 12]]
        for order in range(0, 6):
            out = ndimage.affine_transform(data, [[2, 0], [0, 2]], 0, (1, 2),
                                           order=order)
            assert_array_almost_equal(out, [[1, 3]]) 

def test_affine_transform17(self):
        data = [[1, 2, 3, 4],
                [5, 6, 7, 8],
                [9, 10, 11, 12]]
        for order in range(0, 6):
            out = ndimage.affine_transform(data, [[0.5, 0], [0, 1]], 0,
                                           (6, 4), order=order)
            assert_array_almost_equal(out[::2, ...], data) 

def test_affine_transform18(self):
        data = [[1, 2, 3, 4],
                [5, 6, 7, 8],
                [9, 10, 11, 12]]
        for order in range(0, 6):
            out = ndimage.affine_transform(data, [[0.5, 0], [0, 0.5]], 0,
                                           (6, 8), order=order)
            assert_array_almost_equal(out[::2, ::2], data) 

def test_affine_transform19(self):
        data = numpy.array([[1, 2, 3, 4],
                            [5, 6, 7, 8],
                            [9, 10, 11, 12]], numpy.float64)
        for order in range(0, 6):
            out = ndimage.affine_transform(data, [[0.5, 0], [0, 0.5]], 0,
                                           (6, 8), order=order)
            out = ndimage.affine_transform(out, [[2.0, 0], [0, 2.0]], 0,
                                           (3, 4), order=order)
            assert_array_almost_equal(out, data) 

def test_affine_transform20(self):
        data = [[1, 2, 3, 4],
                [5, 6, 7, 8],
                [9, 10, 11, 12]]
        for order in range(0, 6):
            out = ndimage.affine_transform(data, [[0], [2]], 0, (2,),
                                           order=order)
            assert_array_almost_equal(out, [1, 3]) 

def test_affine_transform21(self):
        data = [[1, 2, 3, 4],
                [5, 6, 7, 8],
                [9, 10, 11, 12]]
        for order in range(0, 6):
            out = ndimage.affine_transform(data, [[2], [0]], 0, (2,),
                                           order=order)
            assert_array_almost_equal(out, [1, 9]) 

def test_affine_transform23(self):
        # shift and offset interaction; see issue #1547
        data = numpy.array([4, 1, 3, 2])
        for order in range(0, 6):
            out = ndimage.affine_transform(data, [[0.5]], [-1], (8,),
                                           order=order)
            assert_array_almost_equal(out[::2], [0, 4, 1, 3]) 

def test_affine_transform24(self):
        # consistency between diagonal and non-diagonal case; see issue #1547
        data = numpy.array([4, 1, 3, 2])
        for order in range(0, 6):
            with suppress_warnings() as sup:
                sup.filter(UserWarning,
                           "The behaviour of affine_transform with a one-dimensional array .* has changed")
                out1 = ndimage.affine_transform(data, [2], -1, order=order)
            out2 = ndimage.affine_transform(data, [[2]], -1, order=order)
            assert_array_almost_equal(out1, out2) 

def test_affine_transform25(self):
        # consistency between diagonal and non-diagonal case; see issue #1547
        data = numpy.array([4, 1, 3, 2])
        for order in range(0, 6):
            with suppress_warnings() as sup:
                sup.filter(UserWarning,
                           "The behaviour of affine_transform with a one-dimensional array .* has changed")
                out1 = ndimage.affine_transform(data, [0.5], -1, order=order)
            out2 = ndimage.affine_transform(data, [[0.5]], -1, order=order)
            assert_array_almost_equal(out1, out2) 

def test_uint64_max():
    # Test interpolation respects uint64 max.  Reported to fail at least on
    # win32 (due to the 32 bit visual C compiler using signed int64 when
    # converting between uint64 to double) and Debian on s390x.
    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)) 

def test_affine_transform27(self):
        # test valid homogeneous transformation matrix
        data = numpy.array([[4, 1, 3, 2],
                            [7, 6, 8, 5],
                            [3, 5, 3, 6]])
        tform_h1 = numpy.hstack((numpy.eye(2), -numpy.ones((2, 1))))
        tform_h2 = numpy.vstack((tform_h1, [[5, 2, 1]]))
        assert_raises(ValueError, ndimage.affine_transform, data, tform_h2) 

def test_affine_transform_multi_d_endianness_with_output_parameter(self):
        # affine transform given output ndarray or dtype with either endianness
        # see issue #4127
        data = numpy.array([1])
        for out in [data.dtype, data.dtype.newbyteorder(),
                    numpy.empty_like(data),
                    numpy.empty_like(data).astype(data.dtype.newbyteorder())]:
            returned = ndimage.affine_transform(data, [[1]], output=out)
            result = out if returned is None else returned
            assert_array_almost_equal(result, [1]) 

def test_zoom_affine01(self):
        data = [[1, 2, 3, 4],
                [5, 6, 7, 8],
                [9, 10, 11, 12]]
        for order in range(0, 6):
            with suppress_warnings() as sup:
                sup.filter(UserWarning,
                           "The behaviour of affine_transform with a one-dimensional array .* has changed")
                out = ndimage.affine_transform(data, [0.5, 0.5], 0,
                                               (6, 8), order=order)
            assert_array_almost_equal(out[::2, ::2], data) 

def test_affine(x):
    M = np.eye(4)
    np.random.seed(42)
    M[:3] += .1 * np.random.uniform(-1, 1, (3, 4))    
    out1 = affine(x, M, interpolation = "nearest")
    out2 = ndimage.affine_transform(x, M, order=0, prefilter=False)
    return out1,out2 

def test_affine_reshape(x):
    M = np.eye(4)
    np.random.seed(42)
    M[:3] += .1 * np.random.uniform(-1, 1, (3, 4))
    output_shape = (33,45,97)
    out1 = affine(x, M, interpolation = "linear", output_shape = output_shape)
    out2 = ndimage.affine_transform(x, M, order=1, prefilter=False, output_shape = output_shape)
    return out1,out2 

def pitch_time_deformation_spec(spec, stretch_rate=None, shift_rate=None):
    """ Based on https://github.com/f0k/ismir2018/blob/master/experiments/augment.py """
    shape = spec.shape
    random = (np.random.rand(2) - .5) * 2

    warnings.filterwarnings("ignore", "The behaviour of affine_transform with a "
                            "one-dimensional array supplied for the matrix parameter "
                            "has changed in scipy 0.18.0.", module="scipy.ndimage")

    if stretch_rate is None:
        stretch_rate = 1 + random[0] * config.MAX_STRETCHING

    if shift_rate is None:
        shift_rate = 1 + random[1] * config.MAX_SHIFTING_PITCH

    new_length = int(shape[1] * stretch_rate)

    # We can do shifting/stretching and cropping in a single affine
    # transform (including setting the upper bands to zero if we shift
    # down the signal so far that we exceed the input's nyquist rate)

    spec = spline_filter(spec.T, 2).astype(spec.dtype)
    spec = affine_transform(spec, (1 / stretch_rate, 1 / shift_rate),
                            output_shape=(new_length, shape[0]),
                            offset=(0, 0), mode='constant', order=2,
                            prefilter=False)

    # clip possible negative values introduced by the interpolation
    # delete the last frame if it is empty

    if (spec[-1] == 0).all():
        spec = spec[:-1]
    return np.maximum(spec.T, 0), stretch_rate 

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.int64,
           np.intp, np.uintp, 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)