def test_slice_from_3d():
    # Resample a 3d image, returning a zslice, yslice and xslice
    #
    # This example creates a coordmap that coincides with
    # a given z, y, or x slice of an image, and checks that
    # resampling agrees with the data in the given slice.
    shape = (100,90,80)
    g = AffineTransform.from_params('ijk',
                                    'xyz',
                                    np.diag([0.5,0.5,0.5,1]))
    img = Image(np.ones(shape), g)
    img.get_data()[50:55,40:55,30:33] = 3
    I = np.identity(4)
    zsl = slices.zslice(26,
                        ((0,49.5), 100),
                        ((0,44.5), 90),
                        img.reference)
    ir = resample(img, zsl, I, (100, 90))
    assert_array_almost_equal(ir.get_data(), img[:,:,53].get_data())
    ysl = slices.yslice(22,
                        ((0,49.5), 100),
                        ((0,39.5), 80),
                        img.reference)
    ir = resample(img, ysl, I, (100, 80))
    assert_array_almost_equal(ir.get_data(), img[:,45,:].get_data())
    xsl = slices.xslice(15.5,
                        ((0,44.5), 90),
                        ((0,39.5), 80),
                        img.reference)
    ir = resample(img, xsl, I, (90, 80))
    assert_array_almost_equal(ir.get_data(), img[32,:,:].get_data())

def resamp(nipy_img, onto_aff, onto_shape, order=3, cval=-1, w2wmap=None, debug=0):
    """
    utility function to resample an image
    input: 
	nipy image
    	onto affine
	onto shape
	order : scipy resample interpolation order
	cval: value outside when resampling to bigger box
    """
    arraycoo = 'ijklmnopq'[:len(onto_shape)]
    spacecoo = 'xyztrsuvw'[:len(onto_shape)] 
    if debug: print('\narraycoo: ', arraycoo,  
                     '\nspacecoo: ', spacecoo, '\nonto_aff\n', onto_aff)
    
    dmaker = CoordSysMaker(arraycoo, 'generic-array')
    rmaker = CoordSysMaker(spacecoo, 'generic-scanner')
    cm_maker = cmap.CoordMapMaker(dmaker, rmaker)
    cmap_out = cm_maker.make_affine(onto_aff)
    if debug: print('cmap_out:\n',cmap_out)
    
    if w2wmap == None: w2wmap = np.eye(onto_aff.shape[0])
    if debug: print('w2wmap:\n',w2wmap)
        
    return resample(nipy_img, cmap_out, w2wmap, onto_shape, order=order, cval=cval)

def test_nonaffine():
    # resamples an image along a curve through the image.
    #
    # FIXME: use the reference.evaluate.Grid to perform this nicer
    # FIXME: Remove pylab references
    def curve(x):  # function accept N by 1, returns N by 2
        return np.vstack([5 * np.sin(x.T), 5 * np.cos(x.T)]).T + [52, 47]

    for names in (("xy", "ij", "t", "u"), ("ij", "xy", "t", "s")):
        in_names, out_names, tin_names, tout_names = names
        g = AffineTransform.from_params(in_names, out_names, np.identity(3))
        img = Image(np.ones((100, 90)), g)
        img[50:55, 40:55] = 3.0
        tcoordmap = AffineTransform.from_start_step(tin_names, tout_names, [0], [np.pi * 1.8 / 100])
        ir = resample(img, tcoordmap, curve, (100,))
    if gui_review:
        import pylab

        pylab.figure(num=3)
        pylab.imshow(img, interpolation="nearest")
        d = curve(np.linspace(0, 1.8 * np.pi, 100))
        pylab.plot(d[0], d[1])
        pylab.gca().set_ylim([0, 99])
        pylab.gca().set_xlim([0, 89])
        pylab.figure(num=4)
        pylab.plot(np.asarray(ir))

def resample_image(source_file, target_file, outdir, w2wmap=None, order=3,
                   cval=0, verbose=0):
    """ Resample the source image to match the target image using Nipy.

    Parameters
    ----------
    source_file: str (mandatory)
        the image to resample.
    target_file: str (mandatory)
        the reference image.
    outdir: str (mandatory)
        the folder where the resampled image will be saved.
    w2wmap: array (4, 4) or callable
        physical to physical transformation.
    verbose: int (optional, default 0)
        the verbosity level.

    Returns
    -------
    resampled_file: str
        the resampled image.
    """
    # Get target image information
    target_image = nipy.load_image(target_file)
    onto_shape = target_image.shape[:3]
    onto_aff = xyz_affine(target_image.affine, xyz=[0, 1, 2], verbose=verbose)

    # Define index and physical coordinate systems
    arraycoo = "ijklmnopq"[:len(onto_shape)]
    spacecoo = "xyztrsuvw"[:len(onto_shape)]
    if verbose > 0:
        print("\narraycoo: ", arraycoo, "\nspacecoo: ", spacecoo,
              "\nonto_aff\n", onto_aff)
    dmaker = CoordSysMaker(arraycoo, 'generic-array')
    rmaker = CoordSysMaker(spacecoo, 'generic-scanner')
    cm_maker = cmap.CoordMapMaker(dmaker, rmaker)
    cmap_out = cm_maker.make_affine(onto_aff)
    if verbose > 0:
        print("cmap_out:\n", cmap_out)

    # Define the default physical to physical transformation
    if w2wmap is None:
        w2wmap = np.eye(onto_aff.shape[0])
    if verbose > 0:
        print("w2wmap:\n", w2wmap)

    # Resample
    source_image = nipy.load_image(source_file)
    resampled_image = resample(
        source_image, cmap_out, w2wmap, onto_shape, order=order, cval=cval)

    # Save the resampled image
    resampled_file = os.path.join(
        outdir, "resampled_{0}".format(os.path.basename(source_file)))
    nipy.save_image(resampled_image, resampled_file)

    return resampled_file

def test_resample2d2():
    g = AffineTransform.from_params("ij", "xy", np.diag([0.5, 0.5, 1]))
    i = Image(np.ones((100, 90)), g)
    i[50:55, 40:55] = 3.0
    a = np.identity(3)
    a[:2, -1] = 4.0
    A = np.identity(2)
    b = np.ones(2) * 4
    ir = resample(i, i.coordmap, (A, b), (100, 90))
    yield assert_array_almost_equal, ir[42:47, 32:47], 3.0

def test_resample2d3():
    # Same as test_resample2d, only a different way of specifying
    # the transform: here it is an (A,b) pair
    g = AffineTransform.from_params('ij', 'xy', np.diag([0.5,0.5,1]))
    i = Image(np.ones((100,90)), g)
    i.get_data()[50:55,40:55] = 3.
    a = np.identity(3)
    a[:2,-1] = 4.
    ir = resample(i, i.coordmap, a, (100,90))
    assert_array_almost_equal(ir.get_data()[42:47,32:47], 3.)

def test_resample2d2():
    g = AffineTransform.from_params('ij', 'xy', np.diag([0.5,0.5,1]))
    i = Image(np.ones((100,90)), g)
    i.get_data()[50:55,40:55] = 3.
    a = np.identity(3)
    a[:2,-1] = 4.
    A = np.identity(2)
    b = np.ones(2)*4
    ir = resample(i, i.coordmap, (A, b), (100,90))
    assert_array_almost_equal(ir.get_data()[42:47,32:47], 3.)

def test_resample2d3():
    # Same as test_resample2d, only a different way of specifying
    # the transform: here it is an (A,b) pair
    g = AffineTransform.from_params("ij", "xy", np.diag([0.5, 0.5, 1]))
    i = Image(np.ones((100, 90)), g)
    i[50:55, 40:55] = 3.0
    a = np.identity(3)
    a[:2, -1] = 4.0
    ir = resample(i, i.coordmap, a, (100, 90))
    yield assert_array_almost_equal, ir[42:47, 32:47], 3.0

def test_slice_from_3d():
    # Resample a 3d image, returning a zslice, yslice and xslice
    #
    # This example creates a coordmap that coincides with
    # the 10th slice of an image, and checks that
    # resampling agrees with the data in the 10th slice.
    shape = (100, 90, 80)
    g = AffineTransform.from_params("ijk", "xyz", np.diag([0.5, 0.5, 0.5, 1]))
    i = Image(np.ones(shape), g)
    i[50:55, 40:55, 30:33] = 3
    a = np.identity(4)
    zsl = slices.zslice(26, ((0, 44.5), 90), ((0, 39.5), 80), i.reference)
    ir = resample(i, zsl, a, (90, 80))
    yield assert_true(np.allclose(np.asarray(ir), np.asarray(i[53])))
    ysl = slices.yslice(22, (0, 49.5), (0, 39.5), i.reference, (100, 80))
    ir = resample(i, ysl.coordmap, a, ysl.shape)
    yield assert_true(np.allclose(np.asarray(ir), np.asarray(i[:, 45])))
    xsl = slices.xslice(15.5, (0, 49.5), (0, 44.5), i.reference, (100, 90))
    ir = resample(i, xsl.coordmap, a, xsl.shape)
    yield assert_true(np.allclose(np.asarray(ir), np.asarray(i[:, :, 32])))

def test_rotate2d():
    # Rotate an image in 2d on a square grid, should result in transposed image
    g = AffineTransform.from_params('ij', 'xy', np.diag([0.7,0.5,1]))
    g2 = AffineTransform.from_params('ij', 'xy', np.diag([0.5,0.7,1]))
    i = Image(np.ones((100,100)), g)
    # This sets the image data by writing into the array
    i.get_data()[50:55,40:55] = 3.
    a = np.array([[0,1,0],
                  [1,0,0],
                  [0,0,1]], np.float)
    ir = resample(i, g2, a, (100, 100))
    assert_array_almost_equal(ir.get_data().T, i.get_data())

def test_resample_outvalue():
    # Test resampling with different modes, constant values, datatypes, orders

    def func(xyz):
        return xyz + np.asarray([1,0,0])

    coordmap =  vox2mni(np.eye(4))
    arr = np.arange(3 * 3 * 3).reshape(3, 3, 3)
    aff = np.eye(4)
    aff[0, 3] = 1.  # x translation
    for mapping, dt, order in product(
        [aff, func],
        [np.int8, np.intp, np.int32, np.int64, np.float32, np.float64],
        [0, 1, 3]):
        img = Image(arr.astype(dt), coordmap)
        # Test constant value of 0
        img2 = resample(img, coordmap, mapping, img.shape,
                        order=order, mode='constant', cval=0.)
        exp_arr = np.zeros(arr.shape)
        exp_arr[:-1, :, :] = arr[1:, :, :]
        assert_array_almost_equal(img2.get_data(), exp_arr)
        # Test constant value of 1
        img2 = resample(img, coordmap, mapping, img.shape,
                        order=order, mode='constant', cval=1.)
        exp_arr[-1, :, :] = 1
        assert_array_almost_equal(img2.get_data(), exp_arr)
        # Test nearest neighbor
        img2 = resample(img, coordmap, mapping, img.shape,
                        order=order, mode='nearest')
        exp_arr[-1, :, :] = arr[-1, :, :]
        assert_array_almost_equal(img2.get_data(), exp_arr)
    # Test img2img
    target_coordmap = vox2mni(aff)
    target = Image(arr, target_coordmap)
    img2 = resample_img2img(img, target, 3, 'nearest')
    assert_array_almost_equal(img2.get_data(), exp_arr)
    img2 = resample_img2img(img, target, 3, 'constant', cval=1.)
    exp_arr[-1, :, :] = 1
    assert_array_almost_equal(img2.get_data(), exp_arr)

def test_2d_from_3d():
    # Resample a 3d image on a 2d affine grid
    # This example creates a coordmap that coincides with
    # the 10th slice of an image, and checks that
    # resampling agrees with the data in the 10th slice.
    shape = (100, 90, 80)
    g = AffineTransform.from_params("ijk", "xyz", np.diag([0.5, 0.5, 0.5, 1]))
    i = Image(np.ones(shape), g)
    i[50:55, 40:55, 30:33] = 3.0
    a = np.identity(4)
    g2 = ArrayCoordMap.from_shape(g, shape)[10]
    ir = resample(i, g2.coordmap, a, g2.shape)
    yield assert_array_almost_equal, np.asarray(ir), np.asarray(i[10])

def test_rotate3d():
    # Rotate / transpose a 3d image on a non-square grid

    g = AffineTransform.from_params("ijk", "xyz", np.diag([0.5, 0.6, 0.7, 1]))
    g2 = AffineTransform.from_params("ijk", "xyz", np.diag([0.5, 0.7, 0.6, 1]))

    shape = (100, 90, 80)
    i = Image(np.ones(shape), g)
    i[50:55, 40:55, 30:33] = 3.0

    a = np.array([[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1.0]])

    ir = resample(i, g2, a, (100, 80, 90))
    yield assert_array_almost_equal, np.transpose(np.asarray(ir), (0, 2, 1)), i

def test_resample2d1():
    # Tests the same as test_resample2d, only using a callable instead of
    # an AffineTransform instance
    g = AffineTransform.from_params('ij', 'xy', np.diag([0.5,0.5,1]))
    i = Image(np.ones((100,90)), g)
    i.get_data()[50:55,40:55] = 3.
    a = np.identity(3)
    a[:2,-1] = 4.
    A = np.identity(2)
    b = np.ones(2)*4
    def mapper(x):
        return np.dot(x, A.T) + b
    ir = resample(i, i.coordmap, mapper, (100,90))
    assert_array_almost_equal(ir.get_data()[42:47,32:47], 3.)

def test_rotate3d():
    # Rotate / transpose a 3d image on a non-square grid
    g = AffineTransform.from_params('ijk', 'xyz', np.diag([0.5,0.6,0.7,1]))
    g2 = AffineTransform.from_params('ijk', 'xyz', np.diag([0.5,0.7,0.6,1]))
    shape = (100,90,80)
    i = Image(np.ones(shape), g)
    i.get_data()[50:55,40:55,30:33] = 3.
    a = np.array([[1,0,0,0],
                  [0,0,1,0],
                  [0,1,0,0],
                  [0,0,0,1.]])
    ir = resample(i, g2, a, (100,80,90))
    assert_array_almost_equal(np.transpose(ir.get_data(), (0,2,1)),
                              i.get_data())

def test_rotate2d2():
    # Rotate an image in 2d on a non-square grid,
    # should result in transposed image

    g = AffineTransform.from_params("ij", "xy", np.diag([0.7, 0.5, 1]))
    g2 = AffineTransform.from_params("ij", "xy", np.diag([0.5, 0.7, 1]))

    i = Image(np.ones((100, 80)), g)
    i[50:55, 40:55] = 3.0

    a = np.array([[0, 1, 0], [1, 0, 0], [0, 0, 1]], np.float)

    ir = resample(i, g2, a, (80, 100))
    yield assert_array_almost_equal, np.asarray(ir).T, i

def test_resample2d():
    g = AffineTransform.from_params('ij', 'xy', np.diag([0.5,0.5,1]))
    i = Image(np.ones((100,90)), g)
    i.get_data()[50:55,40:55] = 3.
    # This mapping describes a mapping from the "target" physical
    # coordinates to the "image" physical coordinates.  The 3x3 matrix
    # below indicates that the "target" physical coordinates are related
    # to the "image" physical coordinates by a shift of -4 in each
    # coordinate.  Or, to find the "image" physical coordinates, given
    # the "target" physical coordinates, we add 4 to each "target
    # coordinate".  The resulting resampled image should show the
    # overall image shifted -8,-8 voxels towards the origin
    a = np.identity(3)
    a[:2,-1] = 4.
    ir = resample(i, i.coordmap, a, (100,90))
    assert_array_almost_equal(ir.get_data()[42:47,32:47], 3.)

def test_resample3d():
    g = AffineTransform.from_params("ijk", "xyz", np.diag([0.5, 0.5, 0.5, 1]))
    shape = (100, 90, 80)
    i = Image(np.ones(shape), g)
    i[50:55, 40:55, 30:33] = 3.0
    # This mapping describes a mapping from the "target" physical
    # coordinates to the "image" physical coordinates.  The 4x4 matrix
    # below indicates that the "target" physical coordinates are related
    # to the "image" physical coordinates by a shift of -4 in each
    # coordinate.  Or, to find the "image" physical coordinates, given
    # the "target" physical coordinates, we add 4 to each "target
    # coordinate".  The resulting resampled image should show the
    # overall image shifted [-6,-8,-10] voxels towards the origin
    a = np.identity(4)
    a[:3, -1] = [3, 4, 5]
    ir = resample(i, i.coordmap, a, (100, 90, 80))
    yield assert_array_almost_equal, ir[44:49, 32:47, 20:23], 3.0

def test_rotate2d3():
    # Another way to rotate/transpose the image, similar to
    # test_rotate2d2 and test_rotate2d except the world of the
    # output coordmap are the same as the world of the
    # original image. That is, the data is transposed on disk, but the
    # output coordinates are still 'x,'y' order, not 'y', 'x' order as
    # above

    # this functionality may or may not be used a lot. if data is to
    # be transposed but one wanted to keep the NIFTI order of output
    # coords this would do the trick

    g = AffineTransform.from_params("xy", "ij", np.diag([0.5, 0.7, 1]))
    i = Image(np.ones((100, 80)), g)
    i[50:55, 40:55] = 3.0

    a = np.identity(3)
    g2 = AffineTransform.from_params("xy", "ij", np.array([[0, 0.5, 0], [0.7, 0, 0], [0, 0, 1]]))
    ir = resample(i, g2, a, (80, 100))
    v2v = compose(g.inverse(), g2)
    yield assert_array_almost_equal, np.asarray(ir).T, i

def test_rotate2d3():
    # Another way to rotate/transpose the image, similar to
    # test_rotate2d2 and test_rotate2d, except the world of the
    # output coordmap is the same as the world of the
    # original image. That is, the data is transposed on disk, but the
    # output coordinates are still 'x,'y' order, not 'y', 'x' order as
    # above

    # this functionality may or may not be used a lot. if data is to
    # be transposed but one wanted to keep the NIFTI order of output
    # coords this would do the trick
    g = AffineTransform.from_params('xy', 'ij', np.diag([0.5,0.7,1]))
    i = Image(np.ones((100,80)), g)
    # This sets the image data by writing into the array
    i.get_data()[50:55,40:55] = 3.
    a = np.identity(3)
    g2 = AffineTransform.from_params('xy', 'ij', np.array([[0,0.5,0],
                                                  [0.7,0,0],
                                                  [0,0,1]]))
    ir = resample(i, g2, a, (80,100))
    assert_array_almost_equal(ir.get_data().T, i.get_data())