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.get_data()[50:55,40:55] = 3.
        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(ir.get_data())

   def __init__(self, data, affine, axis_names, metadata={}, 
                lps=True):
      """ Creates a new nipy image with an affine mapping.

      Parameters
      ----------

      data : ndarray
         ndarray representing the data.

      affine : 4x4 ndarray
         affine transformation to the reference coordinate system

      axis_names : [string]
         names of the axes in the coordinate system.
      """

      if len(axis_names) < 3:
         raise ValueError('XYZImage must have a minimum of 3 axes')

      # The first three axes are assumed to be the
      # spatial ones
      xyz_transform = XYZTransform(affine, axis_names[:3], lps)
      nonspatial_names = axis_names[3:]
        
      if nonspatial_names:
         nonspatial_affine_transform = AffineTransform.from_start_step(nonspatial_names, nonspatial_names, [0]*(data.ndim-3), [1]*(data.ndim-3))
         full_dimensional_affine_transform = cmap_product(xyz_transform, nonspatial_affine_transform)
      else:
         full_dimensional_affine_transform = xyz_transform 

      self._xyz_transform = xyz_transform

      Image.__init__(self, data, full_dimensional_affine_transform,
                     metadata=metadata)

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 __init__(self, data, affine, coord_sys, metadata=None):
        """ Creates a new nipy image with an affine mapping.

            Parameters
            ----------

            data : ndarray
                ndarray representing the data.
            affine : 4x4 ndarray
                affine transformation to the reference coordinate system
            coord_system : string
                name of the reference coordinate system.
        """

        function_domain = CoordinateSystem(['axis%d' % i for i in range(3)], 
                                        name=coord_sys)
        function_range = CoordinateSystem(['x','y','z'], name='world')
        spatial_coordmap = AffineTransform(function_domain, function_range,
                                           affine)

        nonspatial_names = ['axis%d' % i for i in range(3, data.ndim)]
        if nonspatial_names:
            nonspatial_coordmap = AffineTransform.from_start_step(nonspatial_names, nonspatial_names, [0]*(data.ndim-3), [1]*(data.ndim-3))
            full_coordmap = cmap_product(coordmap, nonspatial_coordmap)
        else:
            full_coordmap = spatial_coordmap 

        self._spatial_coordmap = spatial_coordmap

        self.coord_sys = coord_sys
        Image.__init__(self, data, full_coordmap) 
        if metadata is not None:
            self.metadata = metadata

    def get_nifti(self, topo_view, base_nifti=None, **kwargs):
        """
        Process the nifti

        Parameters
        ----------
        topo_view: array-like
            Topological view to create nifti. 3D.

        Returns
        -------
        image: nipy image
            Nifti image from topological view.
        """
        if base_nifti is None:
            assert self.base_nifti is not None, ("`base.nii` not in dataset "
                                                 "directory. You may need to "
                                                 "reprocess.")
            base_nifti = self.base_nifti
            image = Image.from_image(base_nifti, data=topo_view)
        else:
            if isinstance(base_nifti, str):
                base_nifti = load_image(base_nifti)
            base2new_affine = np.linalg.inv(
                base_nifti.affine).dot(self.base_nifti.affine)
            cmap = AffineTransform("kji", "zxy", base2new_affine)
            image = Image.from_image(base_nifti, data=topo_view, coordmap=cmap)

        return image

def load(filename):
    """Load an image from the given filename.

    Parameters
    ----------
    filename : string
        Should resolve to a complete filename path.

    Returns
    -------
    image : An `Image` object
        If successful, a new `Image` object is returned.

    See Also
    --------
    save_image : function for saving images
    fromarray : function for creating images from numpy arrays

    Examples
    --------

    >>> from nipy.io.api import load_image
    >>> from nipy.testing import anatfile
    >>> img = load_image(anatfile)
    >>> img.shape
    (33, 41, 25)
    """
    img = formats.load(filename)
    aff = img.get_affine()
    shape = img.get_shape()
    hdr = img.get_header()

    # Get info from NIFTI header, if present, to tell which axes are
    # which.  This is a NIFTI-specific kludge, that might be abstracted
    # out into the image backend in a general way.  Similarly for
    # getting zooms

    # axis_renames is a dictionary: dict([(int, str)])
    # that has keys in range(3)
    # the axes of the Image are renamed from 'ijk'
    # using these names

    try:
        axis_renames = hdr.get_axis_renames()
    except (TypeError, AttributeError):
        axis_renames = {}

    try:
        zooms = hdr.get_zooms()
    except AttributeError:
        zooms = np.ones(len(shape))

    # affine_transform is a 3-d transform

    affine_transform3d, affine_transform = \
        affine_transform_from_array(aff, 'ijk', pixdim=zooms[3:])
    img = Image(img.get_data(), affine_transform.renamed_domain(axis_renames))
    img.header = hdr
    return img

def test_labels1():
    img = load_image(funcfile)
    data = img.get_data()
    parcelmap = Image(img[0].get_data(), AfT("kji", "zyx", np.eye(4)))
    parcelmap = (parcelmap.get_data() * 100).astype(np.int32)
    v = 0
    for i, d in axis0_generator(data, parcels(parcelmap)):
        v += d.shape[1]
    assert_equal(v, parcelmap.size)

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.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 load(filename):
    """Load an image from the given filename.

    Parameters
    ----------
    filename : string
        Should resolve to a complete filename path.

    Returns
    -------
    image : An `Image` object
        If successful, a new `Image` object is returned.

    See Also
    --------
    save_image : function for saving images
    fromarray : function for creating images from numpy arrays

    Examples
    --------

    >>> from nipy.io.api import load_image
    >>> from nipy.testing import anatfile
    >>> img = load_image(anatfile)
    >>> img.shape
    (33, 41, 25)
    """
    img = nib.load(filename)
    aff = img.get_affine()
    shape = img.get_shape()
    hdr = img.get_header()
    # If the header implements it, get a list of names, one per axis,
    # and put this into the coordinate map.  In fact, no image format
    # implements this at the moment, so in practice, the following code
    # is not currently called. 
    axis_renames = {}
    try:
        axis_names = hdr.axis_names
    except AttributeError:
        pass
    else:
        # axis_renames is a dictionary: dict([(int, str)]) that has keys
        # in range(3). The axes of the Image are renamed from 'ijk' using
        # these names
        for i in range(min([len(axis_names), 3])):
            name = axis_names[i]
            if not (name is None or name == ''):
                axis_renames[i] = name
    zooms = hdr.get_zooms()
    # affine_transform is a 3-d transform
    affine_transform3d, affine_transform = \
        affine_transform_from_array(aff, 'ijk', pixdim=zooms[3:])
    img = Image(img.get_data(), affine_transform.renamed_domain(axis_renames))
    img.header = hdr
    return img

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_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.get_data()[50:55,40:55,30:33] = 3.
    a = np.identity(4)
    g2 = ArrayCoordMap.from_shape(g, shape)[10]
    ir = resample(i, g2.coordmap, a, g2.shape)
    assert_array_almost_equal(ir.get_data(), i[10].get_data())

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_rollaxis():
    data = np.random.standard_normal((3,4,7,5))
    im = Image(data, AffineTransform.from_params('ijkl', 'xyzt', np.diag([1,2,3,4,1])))

    # for the inverse we must specify an integer
    yield assert_raises, ValueError, image.rollaxis, im, 'i', True

    # Check that rollaxis preserves diagonal affines, as claimed

    yield assert_almost_equal, image.rollaxis(im, 1).affine, np.diag([2,1,3,4,1])
    yield assert_almost_equal, image.rollaxis(im, 2).affine, np.diag([3,1,2,4,1])
    yield assert_almost_equal, image.rollaxis(im, 3).affine, np.diag([4,1,2,3,1])

    # Check that ambiguous axes raise an exception
    # 'l' appears both as an axis and a reference coord name
    # and in different places

    im_amb = Image(data, AffineTransform.from_params('ijkl', 'xylt', np.diag([1,2,3,4,1])))
    yield assert_raises, ValueError, image.rollaxis, im_amb, 'l'

    # But if it's unambiguous, then
    # 'l' can appear both as an axis and a reference coord name

    im_unamb = Image(data, AffineTransform.from_params('ijkl', 'xyzl', np.diag([1,2,3,4,1])))
    im_rolled = image.rollaxis(im_unamb, 'l')
    yield assert_almost_equal, im_rolled.get_data(), \
        im_unamb.get_data().transpose([3,0,1,2])

    for i, o, n in zip('ijkl', 'xyzt', range(4)):
        im_i = image.rollaxis(im, i)
        im_o = image.rollaxis(im, o)
        im_n = image.rollaxis(im, n)

        yield assert_almost_equal, im_i.get_data(), \
                                  im_o.get_data()

        yield assert_almost_equal, im_i.affine, \
            im_o.affine

        yield assert_almost_equal, im_n.get_data(), \
            im_o.get_data()

        for _im in [im_n, im_o, im_i]:
            im_n_inv = image.rollaxis(_im, n, inverse=True)

            yield assert_almost_equal, im_n_inv.affine, \
                im.affine

            yield assert_almost_equal, im_n_inv.get_data(), \
                im.get_data()

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 load(filename):
    """Load an image from the given filename.

    Parameters
    ----------
    filename : string
        Should resolve to a complete filename path.

    Returns
    -------
    image : An `Image` object
        If successful, a new `Image` object is returned.

    See Also
    --------
    save_image : function for saving images
    fromarray : function for creating images from numpy arrays

    Examples
    --------

    >>> from nipy.io.api import load_image
    >>> from nipy.testing import anatfile
    >>> img = load_image(anatfile)
    >>> img.shape
    (33, 41, 25)
    """
    img = formats.load(filename)
    aff = img.get_affine()
    shape = img.get_shape()
    hdr = img.get_header()
    # Get info from NIFTI header, if present, to tell which axes are
    # which.  This is a NIFTI-specific kludge, that might be abstracted
    # out into the image backend in a general way.  Similarly for
    # getting zooms
    try:
        fps = hdr.get_dim_info()
    except (TypeError, AttributeError):
        fps = (None, None, None)
    ijk = ijk_from_fps(fps)
    try:
        zooms = hdr.get_zooms()
    except AttributeError:
        zooms = np.ones(len(shape))
    aff = _match_affine(aff, len(shape), zooms)
    coordmap = coordmap_from_affine(aff, ijk)
    img = Image(img.get_data(), coordmap)
    img.header = hdr
    return img

def sources_to_nifti(CHECKPOINT, MASKMAT, BASENIFTI, ONAME, savepath, voxels, win):
    bnifti = load_image(BASENIFTI)
    mask = loadmat(MASKMAT)['mask']
    model = np.load(CHECKPOINT) # Numpy array of sources from Infomax ICA

    for i in range(len(model)): # Goes component by component

        W = model[i,:].reshape([voxels,win])

        f = zeros(len(mask))
        idx = where(mask==1)
        data = zeros((bnifti.shape[0],bnifti.shape[1],bnifti.shape[2],W.shape[1]))

        f[idx[0].tolist()] = detrend(W)/std(W)

        for j in range(0,W.shape[1]):
            data[:,:,:,j] = reshape(f,(bnifti.shape[0],bnifti.shape[1],bnifti.shape[2] ), order='F')

        img = Image.from_image(bnifti,data=data)

        os.chdir(savepath)

        fn = ONAME + "%s.nii" % (str(i)) # Where result should be saved and under what name

        save_image(img,fn)

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.get_data()[50:55,40:55,30:33] = 3.
    # 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))
    assert_array_almost_equal(ir.get_data()[44:49,32:47,20:23], 3.)

def test_synchronized_order():

    data = np.random.standard_normal((3,4,7,5))
    im = Image(data, AffineTransform.from_params('ijkl', 'xyzt', np.diag([1,2,3,4,1])))

    im_scrambled = im.reordered_axes('iljk').reordered_reference('xtyz')
    im_unscrambled = image.synchronized_order(im_scrambled, im)
    
    yield assert_equal, im_unscrambled.coordmap, im.coordmap
    yield assert_almost_equal, im_unscrambled.get_data(), im.get_data()
    yield assert_equal, im_unscrambled, im
    yield assert_true, im_unscrambled == im
    yield assert_false, im_unscrambled != im

    # the images don't have to be the same shape

    data2 = np.random.standard_normal((3,11,9,4))
    im2 = Image(data, AffineTransform.from_params('ijkl', 'xyzt', np.diag([1,2,3,4,1])))

    im_scrambled2 = im2.reordered_axes('iljk').reordered_reference('xtyz')
    im_unscrambled2 = image.synchronized_order(im_scrambled2, im)

    yield assert_equal, im_unscrambled2.coordmap, im.coordmap

    # or the same coordmap

    data3 = np.random.standard_normal((3,11,9,4))
    im3 = Image(data, AffineTransform.from_params('ijkl', 'xyzt', np.diag([1,9,3,-2,1])))

    im_scrambled3 = im3.reordered_axes('iljk').reordered_reference('xtyz')
    im_unscrambled3 = image.synchronized_order(im_scrambled3, im)

    yield assert_equal, im_unscrambled3.axes, im.axes
    yield assert_equal, im_unscrambled3.reference, im.reference