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 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 setUp(self):
     names = ['zspace', 'yspace', 'xspace']
     shape = (10,20,30)
     self.img = Image(np.zeros(shape), 
                      AffineTransform.from_start_step(names, names, (0,)*3, (1,)*3))
     self.img2 = Image(np.ones(shape), 
                       AffineTransform.from_start_step(names, names, (0,)*3, (1,)*3))
                    
     shape = (3,5,4)
     self.img3 = Image(np.zeros(shape), 
                       AffineTransform.from_start_step(names, names, (0,)*3, (1,)*3))
     self.img4 = Image(np.zeros(shape), 
                       AffineTransform.from_start_step(names, names, (0,)*3, (1,)*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_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_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 __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_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_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 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_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())

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_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_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 output_resid(outfile, fmri_image, clobber=False):
    """
    Create an output file of the residuals parameter from the OLS pass of
    fmristat.

    Uses affine part of the first image to output resids unless
    fmri_image is an Image.

    Parameters
    ----------
    outfile :
    fmri_image : ``FmriImageList`` or 4D image
       If ``FmriImageList``, needs attributes ``volume_start_times``,
       supports len(), and object[0] has attributes ``affine``,
       ``coordmap`` and ``shape``, from which we create a new 4D
       coordmap and shape
       If 4D image, use the images coordmap and shape
    clobber : bool
       if True, overwrite previous output

    Returns
    -------
    regression_output : 
    """

    if isinstance(fmri_image, FmriImageList):
        n = len(fmri_image.list)
        T = np.zeros((5,5))
        g = fmri_image[0].coordmap
        T[1:,1:] = fmri_image[0].affine
        T[0,0] = (fmri_image.volume_start_times[1:] - 
                  fmri_image.volume_start_times[:-1]).mean()
        # FIXME: NIFTI specific naming here
        innames = ["l"] + list(g.input_coords.coord_names)
        outnames = ["t"] + list(g.output_coords.coord_names)
        cmap = AffineTransform.from_params(innames,
                                  outnames, T)
        shape = (n,) + fmri_image[0].shape
    elif isinstance(fmri_image, Image):
        cmap = fmri_image.coordmap
        shape = fmri_image.shape
    else:
        raise ValueError, "expecting FmriImageList or 4d Image"

    outim = ModelOutputImage(outfile, cmap, shape, clobber=clobber)
    return regression.RegressionOutput(outim, regression.output_resid)

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