def test_load_spaces():
    # Test spaces get read correctly
    shape = np.array((6, 5, 4, 3, 2))
    zooms = np.array((2, 3, 4, 5, 6))
    data = np.random.normal(size=shape)
    # Default with no affine in header, or in image
    ni_img = nib.Nifti1Image(data, None)
    hdr = get_header(ni_img)
    hdr.set_zooms(zooms)
    # Expected affine is from the pixdims and the center of the image.  Default
    # is also flipped X.
    offsets = (1 - shape[:3]) / 2. * zooms[:3] * (-1, 1, 1)
    exp_aff = from_matvec(np.diag([-2, 3, 4, 5, 6]),
                          list(offsets) + [0, 0])
    in_cs = CS('ijktu', name='voxels')
    exp_cmap = AT(in_cs, unknown_csm(5), exp_aff)
    assert_equal(nifti2nipy(ni_img).coordmap, exp_cmap)
    an_aff = from_matvec(np.diag([1.1, 2.2, 3.3]), [10, 11, 12])
    exp_aff = from_matvec(np.diag([1.1, 2.2, 3.3, 5, 6]), [10, 11, 12, 0, 0])
    for label, csm in (('scanner', scanner_csm),
                       ('aligned', aligned_csm),
                       ('talairach', talairach_csm),
                       ('mni', mni_csm)):
        hdr.set_sform(an_aff, label)
        assert_equal(nifti2nipy(ni_img).coordmap, AT(in_cs, csm(5), exp_aff))

def test_mm_scaling():
    # Test the micron and meter scale the affine right
    data = np.random.normal(size=list(range(4)))
    xyz_aff = from_matvec(np.diag([2, 3, 4]), [11, 12, 13])
    exp_aff = from_matvec(np.diag([2, 3, 4, 1]), [11, 12, 13, 0])
    in_cs = CS('ijkt', name='voxels')
    out_cs = aligned_csm(4)
    # No space scaling
    ni_img = nib.Nifti1Image(data, xyz_aff)
    hdr = get_header(ni_img)
    assert_equal(hdr.get_xyzt_units(), ('unknown', 'unknown'))
    assert_equal(nifti2nipy(ni_img).coordmap, AT(in_cs, out_cs, exp_aff))
    # mm is assumed
    hdr.set_xyzt_units('mm')
    assert_equal(nifti2nipy(ni_img).coordmap, AT(in_cs, out_cs, exp_aff))
    # microns !
    hdr.set_xyzt_units('micron')
    scaler = np.diag([1 / 1000., 1 / 1000., 1 / 1000., 1, 1])
    assert_equal(nifti2nipy(ni_img).coordmap,
                 AT(in_cs, out_cs, np.dot(scaler, exp_aff)))
    # mm again !  This test implicitly asserts that the nifti image affine is
    # not being changed by the conversion routine, otherwise we'd pick up the
    # microns scaling above.
    hdr.set_xyzt_units('mm')
    assert_equal(nifti2nipy(ni_img).coordmap, AT(in_cs, out_cs, exp_aff))
    # meters !
    hdr.set_xyzt_units('meter')
    scaler = np.diag([1000., 1000., 1000., 1, 1])
    assert_equal(nifti2nipy(ni_img).coordmap,
                 AT(in_cs, out_cs, np.dot(scaler, exp_aff)))

def test_adapt_affine():
    # Adapt affine to missing or extra input dimensions
    aff_3d = from_matvec(np.arange(9).reshape((3, 3)), [11, 12, 13])
    # For 4x4 affine, 3D image, no-op
    assert_array_equal(adapt_affine(aff_3d, 3), aff_3d)
    # For 4x4 affine, 4D image, add extra identity dimension
    assert_array_equal(adapt_affine(aff_3d, 4),
                       [[ 0,  1,  2,  0,  11],
                        [ 3,  4,  5,  0,  12],
                        [ 6,  7,  8,  0,  13],
                        [ 0,  0,  0,  1,   0],
                        [ 0,  0,  0,  0,   1]])
    # For 5x5 affine, 4D image, identity
    aff_4d = from_matvec(np.arange(16).reshape((4, 4)), [11, 12, 13, 14])
    assert_array_equal(adapt_affine(aff_4d, 4), aff_4d)
    # For 4x4 affine, 2D image, dropped column
    assert_array_equal(adapt_affine(aff_3d, 2),
                       [[ 0,  1,  11],
                        [ 3,  4,  12],
                        [ 6,  7,  13],
                        [ 0,  0,   1]])
    # For 4x4 affine, 1D image, 2 dropped columnn
    assert_array_equal(adapt_affine(aff_3d, 1),
                       [[ 0,  11],
                        [ 3,  12],
                        [ 6,  13],
                        [ 0,   1]])
    # For 3x3 affine, 2D image, identity
    aff_2d = from_matvec(np.arange(4).reshape((2, 2)), [11, 12])
    assert_array_equal(adapt_affine(aff_2d, 2), aff_2d)

def test_save_toffset():
    # Check toffset only gets set for time
    shape = (2, 3, 4, 5, 6, 7)
    data = np.random.normal(size = shape)
    aff = from_matvec(np.diag([2., 3, 4, 5, 6, 7]),
                              [11, 12, 13, 14, 15, 16])
    xyz_names = talairach_csm(3).coord_names
    in_cs = CS('ijklmn')
    for t_name in 't', 'time':
        cmap = AT(in_cs, CS(xyz_names + (t_name, 'q', 'r')), aff)
        ni_img = nipy2nifti(Image(data, cmap))
        assert_equal(get_header(ni_img)['toffset'], 14)
    for time_like in ('hz', 'ppm', 'rads'):
        cmap = AT(in_cs, CS(xyz_names + (time_like, 'q', 'r')), aff)
        ni_img = nipy2nifti(Image(data, cmap))
        assert_equal(get_header(ni_img)['toffset'], 0)
    # Check that non-matching time causes a nifti error when toffset !=0
    shape_shifted = (2, 3, 4, 6, 5, 7)
    for t_name in 't', 'time':
        # No toffset, this is OK
        cmap = AT(CS(('i', 'j', 'k', 'u', t_name, 'v')),
                  CS(xyz_names + ('u', t_name, 'v')),
                  np.diag([3, 4, 5, 6, 0, 7, 1]))
        assert_equal(nipy2nifti(Image(data, cmap)).shape, shape_shifted)
        # toffset with 0 on TR (time) diagonal
        aff_z1 = from_matvec(np.diag([2., 3, 4, 5, 0, 7]),
                             [11, 12, 13, 14, 15, 16])
        cmap = AT(CS(('i', 'j', 'k', 'u', t_name, 'v')),
                  CS(xyz_names + ('u', t_name, 'v')),
                  aff_z1)
        # Default is to fix the zero
        assert_equal(nipy2nifti(Image(data, cmap)).shape,
                     shape_shifted)
        assert_equal(nipy2nifti(Image(data, cmap), fix0=True).shape,
                     shape_shifted)
        # Unless fix0 is False
        assert_raises(NiftiError, nipy2nifti, Image(data, cmap), fix0=False)
        # Fix doesn't work if there is more than one zero row and column
        aff_z2 = from_matvec(np.diag([2., 3, 4, 0, 0, 7]),
                             [11, 12, 13, 14, 15, 16])
        cmap = AT(CS(('i', 'j', 'k', 'u', t_name, 'v')),
                  CS(xyz_names + ('u', t_name, 'v')),
                  aff_z2)
        assert_raises(NiftiError, nipy2nifti, Image(data, cmap), fix0=True)
    # zeros on the diagonal are not a problem for non-time, with toffset,
    # because we don't need to set the 'time' part of the translation vector,
    # and therefore we don't need to know which *output axis* is time-like
    for t_name in 'hz', 'ppm', 'rads':
        cmap = AT(CS(('i', 'j', 'k', 'u', t_name, 'v')),
                  CS(xyz_names + ('u', t_name, 'v')),
                  aff_z1)
        assert_equal(nipy2nifti(Image(data, cmap), fix0=False).shape,
                     shape_shifted)
        cmap = AT(CS(('i', 'j', 'k', 'u', t_name, 'v')),
                  CS(xyz_names + ('u', t_name, 'v')),
                  aff_z2)
        assert_equal(nipy2nifti(Image(data, cmap), fix0=False).shape,
                     shape_shifted)

def test_save_toffset():
    # Check toffset only gets set for time
    shape = (2, 3, 4, 5, 6, 7)
    data = np.random.normal(size = shape)
    aff = from_matvec(np.diag([2., 3, 4, 5, 6, 7]),
                              [11, 12, 13, 14, 15, 16])
    xyz_names = talairach_csm(3).coord_names
    in_cs = CS('ijklmn')
    for t_name in 't', 'time':
        cmap = AT(in_cs, CS(xyz_names + (t_name, 'q', 'r')), aff)
        ni_img = nipy2nifti(Image(data, cmap))
        assert_equal(ni_img.get_header()['toffset'], 14)
    for time_like in ('hz', 'ppm', 'rads'):
        cmap = AT(in_cs, CS(xyz_names + (time_like, 'q', 'r')), aff)
        ni_img = nipy2nifti(Image(data, cmap))
        assert_equal(ni_img.get_header()['toffset'], 0)
    # Check that non-matching time causes a nifti error when toffset !=0
    shape_shifted = (2, 3, 4, 6, 5, 7)
    for t_name in 't', 'time':
        # No toffset, this is OK
        cmap = AT(CS(('i', 'j', 'k', 'u', t_name, 'v')),
                  CS(xyz_names + ('u', t_name, 'v')),
                  np.diag([3, 4, 5, 6, 0, 7, 1]))
        assert_equal(nipy2nifti(Image(data, cmap)).shape, shape_shifted)
        # toffset, non-matching error
        aff_z1 = from_matvec(np.diag([2., 3, 4, 5, 0, 7]),
                             [11, 12, 13, 14, 15, 16])
        cmap = AT(CS(('i', 'j', 'k', 'u', t_name, 'v')),
                  CS(xyz_names + ('u', t_name, 'v')),
                  aff_z1)
        assert_raises(NiftiError, nipy2nifti, Image(data, cmap))
        # Unless fix0 set
        assert_equal(nipy2nifti(Image(data, cmap), fix0=True).shape,
                     shape_shifted)
        # Even this doesn't work if there is more than one zero row and column
        aff_z2 = from_matvec(np.diag([2., 3, 4, 0, 0, 7]),
                             [11, 12, 13, 14, 15, 16])
        cmap = AT(CS(('i', 'j', 'k', 'u', t_name, 'v')),
                  CS(xyz_names + ('u', t_name, 'v')),
                  aff_z2)
        assert_raises(NiftiError, nipy2nifti, Image(data, cmap), fix0=True)
    # No problem for non-time
    for t_name in 'hz', 'ppm', 'rads':
        cmap = AT(CS(('i', 'j', 'k', 'u', t_name, 'v')),
                  CS(xyz_names + ('u', t_name, 'v')),
                  aff)
        assert_equal(nipy2nifti(Image(data, cmap), fix0=True).shape,
                     shape_shifted)

def test_load_dim_info():
    # Test freq, phase, slice get set correctly on load
    data = np.random.normal(size=list(range(3)))
    xyz_aff = from_matvec(np.diag([2, 3, 4]), [11, 12, 13])
    in_cs = CS('ijk', name='voxels')
    out_cs = aligned_csm(3)
    # Just confirm that the default leads to no axis renaming
    ni_img = nib.Nifti1Image(data, xyz_aff)
    hdr = get_header(ni_img)
    assert_equal(hdr.get_dim_info(), (None, None, None))
    assert_equal(nifti2nipy(ni_img).coordmap, AT(in_cs, out_cs, xyz_aff))
    # But now...
    hdr.set_dim_info(freq=1)
    assert_equal(nifti2nipy(ni_img).coordmap,
                 AT(CS(('i', 'freq', 'k'), "voxels"), out_cs, xyz_aff))
    hdr.set_dim_info(freq=2)
    assert_equal(nifti2nipy(ni_img).coordmap,
                 AT(CS(('i', 'j', 'freq'), "voxels"), out_cs, xyz_aff))
    hdr.set_dim_info(phase=1)
    assert_equal(nifti2nipy(ni_img).coordmap,
                 AT(CS(('i', 'phase', 'k'), "voxels"), out_cs, xyz_aff))
    hdr.set_dim_info(slice=0)
    assert_equal(nifti2nipy(ni_img).coordmap,
                 AT(CS(('slice', 'j', 'k'), "voxels"), out_cs, xyz_aff))
    hdr.set_dim_info(freq=1, phase=0, slice=2)
    assert_equal(nifti2nipy(ni_img).coordmap,
                 AT(CS(('phase', 'freq', 'slice'), "voxels"), out_cs, xyz_aff))

def test_orthogonal_dims():
    # Test whether conversion to nifti raises an error for non-orthogonal
    # non-spatial dimensions
    # This affine is all nicely diagonal
    aff = from_matvec(np.diag([2., 3, 4, 5, 6]), [10, 11, 12, 13, 14])
    data = np.random.normal(size=(3, 4, 5, 6, 7))
    img = Image(data, vox2mni(aff))
    def as3d(aff):
        return from_matvec(aff[:3, :3], aff[:3, -1])
    assert_array_equal(get_affine(nipy2nifti(img)), as3d(aff))
    # Non-orthogonal spatial dimensions OK
    aff[:3, :3] = np.random.normal(size=(3, 3))
    img = Image(data, vox2mni(aff))
    assert_array_equal(get_affine(nipy2nifti(img)), as3d(aff))
    # Space must be orthogonal to time etc
    aff[0, 3] = 0.1
    assert_raises(NiftiError, nipy2nifti, img)
    aff[0, 3] = 0
    assert_array_equal(get_affine(nipy2nifti(img)), as3d(aff))
    aff[3, 0] = 0.1
    assert_raises(NiftiError, nipy2nifti, img)
    aff[3, 0] = 0
    assert_array_equal(get_affine(nipy2nifti(img)), as3d(aff))
    aff[4, 0] = 0.1
    assert_raises(NiftiError, nipy2nifti, img)

def test_against_spm_resample():
    # Test resampling against images resampled with SPM12
    # anatomical.nii has a diagonal -2, 2 2 affine;
    # functional.nii has a diagonal -4, 4 4 affine;
    # These are a bit boring, so first add some rotations and translations to
    # the anatomical image affine, and then resample to the first volume in the
    # functional, and compare to the same thing in SPM.
    # See ``make_moved_anat.py`` script in this directory for input to SPM.
    anat = nib.load(pjoin(DATA_DIR, 'anatomical.nii'))
    func = nib.load(pjoin(DATA_DIR, 'functional.nii'))
    some_rotations = euler2mat(0.1, 0.2, 0.3)
    extra_affine = from_matvec(some_rotations, [3, 4, 5])
    moved_anat = nib.Nifti1Image(anat.get_data().astype(float),
                                 extra_affine.dot(anat.affine),
                                 anat.header)
    one_func = nib.Nifti1Image(func.dataobj[..., 0],
                               func.affine,
                               func.header)
    moved2func = resample_from_to(moved_anat, one_func, order=1, cval=np.nan)
    spm_moved = nib.load(pjoin(DATA_DIR, 'resampled_anat_moved.nii'))
    assert_spm_resampling_close(moved_anat, moved2func, spm_moved)
    # Next we resample the rotated anatomical image to output space, and compare
    # to the same operation done with SPM (our own version of 'reorient.m' by
    # John Ashburner).
    moved2output = resample_to_output(moved_anat, 4, order=1, cval=np.nan)
    spm2output = nib.load(pjoin(DATA_DIR, 'reoriented_anat_moved.nii'))
    assert_spm_resampling_close(moved_anat, moved2output, spm2output);

def test_save3():
    # A test to ensure that when a file is saved, the affine
    # and the data agree. In this case, things don't agree:
    # i) the pixdim is off
    # ii) makes the affine off
    step = np.array([3.45,2.3,4.5,6.9])
    shape = (13,5,7,3)
    mni_xyz = mni_csm(3).coord_names
    cmap = AT(CS('jkli'),
              CS(('t',) + mni_xyz[::-1]),
              from_matvec(np.diag([0,3,5,1]), step))
    data = np.random.standard_normal(shape)
    img = api.Image(data, cmap)
    # with InTemporaryDirectory():
    with InTemporaryDirectory():
        save_image(img, TMP_FNAME)
        tmp = load_image(TMP_FNAME)
        # Detach image from file so we can delete it
        data = tmp.get_data().copy()
        img2 = api.Image(data, tmp.coordmap, tmp.metadata)
        del tmp
    assert_equal(tuple([img.shape[l] for l in [3,2,1,0]]), img2.shape)
    a = np.transpose(np.asarray(img), [3,2,1,0])
    assert_false(np.allclose(img.affine, img2.affine))
    assert_true(np.allclose(a, img2.get_data()))

def xslice(x, y_spec, z_spec, world):
    """
    Return an LPS slice through a 3d box with x fixed.

    Parameters
    ----------
    x : float
       The value at which x is fixed.
    y_spec : sequence
       A sequence with 2 values of form ((float, float), int). The
       (float, float) components are the min and max y values; the int
       is the number of points.
    z_spec : sequence
       As for `y_spec` but for z
    world : str or CoordinateSystem CoordSysMaker or XYZSpace
        World 3D space to which resulting coordmap refers

    Returns
    -------
    affine_transform : AffineTransform
       An affine transform that describes an plane in
       LPS coordinates with x fixed.

    Examples
    --------
    >>> y_spec = ([-114,114], 115) # voxels of size 2 in y, starting at -114, ending at 114
    >>> z_spec = ([-70,100], 86) # voxels of size 2 in z, starting at -70, ending at 100
    >>> x30 = xslice(30, y_spec, z_spec, 'scanner')
    >>> x30([0,0])
    array([  30., -114.,  -70.])
    >>> x30([114,85])
    array([  30.,  114.,  100.])
    >>> x30
    AffineTransform(
       function_domain=CoordinateSystem(coord_names=('i_y', 'i_z'), name='slice', coord_dtype=float64),
       function_range=CoordinateSystem(coord_names=('scanner-x=L->R', 'scanner-y=P->A', 'scanner-z=I->S'), name='scanner', coord_dtype=float64),
       affine=array([[   0.,    0.,   30.],
                     [   2.,    0., -114.],
                     [   0.,    2.,  -70.],
                     [   0.,    0.,    1.]])
    )
    >>> bounding_box(x30, (y_spec[1], z_spec[1]))
    ((30.0, 30.0), (-114.0, 114.0), (-70.0, 100.0))
    """
    affine_range = get_world_cs(world)
    (ymin, ymax), yno = y_spec
    y_tick = (ymax-ymin) / (yno - 1.0)
    (zmin, zmax), zno = z_spec
    z_tick = (zmax-zmin) / (zno - 1.0)
    origin = [x, ymin, zmin]
    colvectors = np.asarray([[0, 0],
                             [y_tick, 0],
                             [0, z_tick]])
    T = from_matvec(colvectors, origin)
    affine_domain = CoordinateSystem(['i_y', 'i_z'], 'slice')
    return AffineTransform(affine_domain,
                           affine_range,
                           T)

def yslice(y, x_spec, z_spec, world):
    """ Return a slice through a 3d box with y fixed.

    Parameters
    ----------
    y : float
       The value at which y is fixed.
    x_spec : sequence
       A sequence with 2 values of form ((float, float), int). The
       (float, float) components are the min and max x values; the int
       is the number of points.
    z_spec : sequence
       As for `x_spec` but for z
    world : str or CoordinateSystem CoordSysMaker or XYZSpace
        World 3D space to which resulting coordmap refers

    Returns
    -------
    affine_transform : AffineTransform
       An affine transform that describes an plane in
       LPS coordinates with y fixed.

    Examples
    --------
    >>> x_spec = ([-92,92], 93) # voxels of size 2 in x, starting at -92, ending at 92
    >>> z_spec = ([-70,100], 86) # voxels of size 2 in z, starting at -70, ending at 100
    >>> y70 = yslice(70, x_spec, z_spec, 'mni')
    >>> y70
    AffineTransform(
       function_domain=CoordinateSystem(coord_names=('i_x', 'i_z'), name='slice', coord_dtype=float64),
       function_range=CoordinateSystem(coord_names=('mni-x=L->R', 'mni-y=P->A', 'mni-z=I->S'), name='mni', coord_dtype=float64),
       affine=array([[  2.,   0., -92.],
                     [  0.,   0.,  70.],
                     [  0.,   2., -70.],
                     [  0.,   0.,   1.]])
    )
    >>> y70([0,0])
    array([-92.,  70., -70.])
    >>> y70([92,85])
    array([  92.,   70.,  100.])
    >>> bounding_box(y70, (x_spec[1], z_spec[1]))
    ((-92.0, 92.0), (70.0, 70.0), (-70.0, 100.0))
    """
    affine_range = get_world_cs(world)
    (xmin, xmax), xno = x_spec
    x_tick = (xmax-xmin) / (xno - 1.0)
    (zmin, zmax), zno = z_spec
    z_tick = (zmax-zmin) / (zno - 1.0)
    origin = [xmin, y, zmin]
    colvectors = np.asarray([[x_tick, 0],
                             [0, 0],
                             [0, z_tick]])
    T = from_matvec(colvectors, origin)
    affine_domain = CoordinateSystem(['i_x', 'i_z'], 'slice')
    return AffineTransform(affine_domain,
                           affine_range,
                           T)

def zslice(z, x_spec, y_spec, world):
    """ Return a slice through a 3d box with z fixed.

    Parameters
    ----------
    z : float
       The value at which z is fixed.
    x_spec : sequence
       A sequence with 2 values of form ((float, float), int). The
       (float, float) components are the min and max x values; the int
       is the number of points.
    y_spec : sequence
       As for `x_spec` but for y
    world : str or CoordinateSystem CoordSysMaker or XYZSpace
        World 3D space to which resulting coordmap refers

    Returns
    -------
    affine_transform : AffineTransform
       An affine transform that describes a plane in LPS coordinates with z
       fixed.

    Examples
    --------
    >>> x_spec = ([-92,92], 93) # voxels of size 2 in x, starting at -92, ending at 92
    >>> y_spec = ([-114,114], 115) # voxels of size 2 in y, starting at -114, ending at 114
    >>> z40 = zslice(40, x_spec, y_spec, 'unknown')
    >>> z40
    AffineTransform(
       function_domain=CoordinateSystem(coord_names=('i_x', 'i_y'), name='slice', coord_dtype=float64),
       function_range=CoordinateSystem(coord_names=('unknown-x=L->R', 'unknown-y=P->A', 'unknown-z=I->S'), name='unknown', coord_dtype=float64),
       affine=array([[   2.,    0.,  -92.],
                     [   0.,    2., -114.],
                     [   0.,    0.,   40.],
                     [   0.,    0.,    1.]])
    )
    >>> z40([0,0])
    array([ -92., -114.,   40.])
    >>> z40([92,114])
    array([  92.,  114.,   40.])
    >>> bounding_box(z40, (x_spec[1], y_spec[1]))
    ((-92.0, 92.0), (-114.0, 114.0), (40.0, 40.0))
    """
    affine_range = get_world_cs(world)
    (xmin, xmax), xno = x_spec
    x_tick = (xmax-xmin) / (xno - 1.0)
    (ymin, ymax), yno = y_spec
    y_tick = (ymax-ymin) / (yno - 1.0)
    origin = [xmin, ymin, z]
    colvectors = np.asarray([[x_tick, 0],
                             [0, y_tick],
                             [0, 0]])
    T = from_matvec(colvectors, origin)
    affine_domain = CoordinateSystem(['i_x', 'i_y'], 'slice')
    return AffineTransform(affine_domain,
                           affine_range,
                           T)

def test_is_xyz_affable():
    # Whether there exists an xyz affine for this coordmap
    affine = np.diag([2,4,5,6,1])
    cmap = AffineTransform(VARS['d_cs_r4'], VARS['r_cs_r4'], affine)
    assert_true(is_xyz_affable(cmap))
    assert_false(is_xyz_affable(cmap.reordered_range([3,0,1,2])))
    assert_false(is_xyz_affable(cmap.reordered_domain([3,0,1,2])))
    # Can pass in own validator
    my_valtor = dict(blind='x', leading='y', ditch='z')
    r_cs = CS(('blind', 'leading', 'ditch'), 'fall')
    affine = from_matvec(np.arange(9).reshape((3, 3)), [11, 12, 13])
    cmap = AffineTransform(VARS['d_cs_r3'], r_cs, affine)
    # No xyz affine if we don't use our custom dictionary
    assert_false(is_xyz_affable(cmap))
    # Is if we do
    assert_true(is_xyz_affable(cmap, my_valtor))

def test_save4():
    # Same as test_save3 except we have reordered the 'ijk' input axes.
    shape = (13,5,7,3)
    step = np.array([3.45,2.3,4.5,6.9])
    # When the input coords are in the 'ljki' order, the affines get
    # rearranged.  Note that the 'start' below, must be 0 for
    # non-spatial dimensions, because we have no way to store them in
    # most cases.  For example, a 'start' of [1,5,3,1] would be lost on
    # reload
    mni_xyz = mni_csm(3).coord_names
    cmap = AT(CS('tkji'),
              CS((('t',) + mni_xyz[::-1])),
              from_matvec(np.diag([2., 3, 5, 1]), step))
    data = np.random.standard_normal(shape)
    img = api.Image(data, cmap)
    with InTemporaryDirectory():
        save_image(img, TMP_FNAME)
        tmp = load_image(TMP_FNAME)
        data = tmp.get_data().copy()
        # Detach image from file so we can delete it
        img2 = api.Image(data, tmp.coordmap, tmp.metadata)
        del tmp
    P = np.array([[0,0,0,1,0],
                  [0,0,1,0,0],
                  [0,1,0,0,0],
                  [1,0,0,0,0],
                  [0,0,0,0,1]])
    res = np.dot(P, np.dot(img.affine, P.T))
    # the step part of the affine should be set correctly
    assert_array_almost_equal(res[:4,:4], img2.affine[:4,:4])
    # start in the spatial dimensions should be set correctly
    assert_array_almost_equal(res[:3,-1], img2.affine[:3,-1])
    # start in the time dimension should be 3.45 as in img, because NIFTI stores
    # the time offset in hdr[``toffset``]
    assert_not_equal(res[3,-1], img2.affine[3,-1])
    assert_equal(res[3,-1], 3.45)
    # shapes should be reversed because img has coordinates reversed
    assert_equal(img.shape[::-1], img2.shape)
    # data should be transposed because coordinates are reversed
    assert_array_almost_equal(
           np.transpose(np.asarray(img2),[3,2,1,0]),
           np.asarray(img))
    # coordinate names should be reversed as well
    assert_equal(img2.coordmap.function_domain.coord_names,
                 img.coordmap.function_domain.coord_names[::-1])
    assert_equal(img2.coordmap.function_domain.coord_names,
                 ('i', 'j', 'k', 't'))

def test_other_axes():
    # With a diagonal affine, we can do PCA on any axis
    ncomp = 5
    img = data_dict['fmridata']
    in_coords = list(img.axes.coord_names)
    img_data = img.get_data()
    for axis_no, axis_name in enumerate('ijkt'):
        p = pca_image(img, axis_name, ncomp=ncomp)
        n = img.shape[axis_no]
        bv_key = 'basis_vectors over ' + axis_name
        assert_equal(_rank(p), n - 1)
        assert_equal(p[bv_key].shape, (n, n - 1))
        # We get the expected data back
        dp = pca_array(img_data, axis_no, ncomp=ncomp)
        # We have to make sure the signs are the same; on Windows it seems the
        # signs can flip even between two runs on the same data
        pos_p = img_res2pos1(p, bv_key)
        pos_dp = res2pos1(dp)
        img_bps = pos_p['basis_projections']
        assert_almost_equal(pos_dp['basis_vectors'], pos_p[bv_key])
        assert_almost_equal(pos_dp['basis_projections'], img_bps.get_data())
        # And we've replaced the expected axis
        exp_coords = in_coords[:]
        exp_coords[exp_coords.index(axis_name)] = 'PCA components'
        assert_equal(img_bps.axes.coord_names, exp_coords)
    # If the affine is not diagonal, we'll get an error
    aff = from_matvec(np.arange(16).reshape(4,4))
    nd_cmap = AffineTransform(img.axes, img.reference, aff)
    nd_img = Image(img_data, nd_cmap)
    for axis_name in 'ijkt':
        assert_raises(AxisError, pca_image, nd_img, axis_name)
    # Only for the non-diagonal parts
    aff = np.array([[1, 2, 0, 0, 10],
                    [2, 1, 0, 0, 11],
                    [0, 0, 3, 0, 12],
                    [0, 0, 0, 4, 13],
                    [0, 0, 0, 0, 1]])
    nd_cmap = AffineTransform(img.axes, img.reference, aff)
    nd_img = Image(img_data, nd_cmap)
    for axis_name in 'ij':
        assert_raises(AxisError, pca_image, nd_img, axis_name)
    for axis_name in 'kt':
        p = pca_image(img, axis_name, ncomp=ncomp)
        exp_coords = in_coords[:]
        exp_coords[exp_coords.index(axis_name)] = 'PCA components'
        assert_equal(p['basis_projections'].axes.coord_names, exp_coords)

""" Make anatomical image with altered affine

* Add some rotations and translations to affine;
* Save as ``.nii`` file so SPM can read it.

See ``resample_using_spm.m`` for processing of this generated image by SPM.
"""

import numpy as np

import nibabel as nib
from nibabel.eulerangles import euler2mat
from nibabel.affines import from_matvec

img = nib.load('anatomical.nii')
some_rotations = euler2mat(0.1, 0.2, 0.3)
extra_affine = from_matvec(some_rotations, [3, 4, 5])
moved_anat = nib.Nifti1Image(img.dataobj,
                             extra_affine.dot(img.affine),
                             img.header)
moved_anat.set_data_dtype(np.float32)
nib.save(moved_anat, 'anat_moved.nii')

def test_xyz_affine():
    # Getting an xyz affine from coordmaps
    aff3d = from_matvec(np.arange(9).reshape((3,3)), [15,16,17])
    cmap3d = AffineTransform(VARS['d_cs_r3'], VARS['r_cs_r3'], aff3d)
    rzs = np.c_[np.arange(12).reshape((4,3)), [0,0,0,12]]
    aff4d = from_matvec(rzs, [15,16,17,18])
    cmap4d = AffineTransform(VARS['d_cs_r4'], VARS['r_cs_r4'], aff4d)
    # Simplest case of 3D affine -> affine unchanged
    assert_array_equal(xyz_affine(cmap3d), aff3d)
    # 4D (5, 5) affine -> 3D equivalent
    assert_array_equal(xyz_affine(cmap4d), aff3d)
    # Any dimensions not spatial, AxesError
    r_cs = CS(('mni-x', 'mni-y', 'mni-q'), 'mni')
    funny_cmap = AffineTransform(VARS['d_cs_r3'],r_cs, aff3d)
    assert_raises(AxesError, xyz_affine, funny_cmap)
    r_cs = CS(('mni-x', 'mni-q', 'mni-z'), 'mni')
    funny_cmap = AffineTransform(VARS['d_cs_r3'],r_cs, aff3d)
    assert_raises(AxesError, xyz_affine, funny_cmap)
    # We insist that the coordmap is in output xyz order
    permutations = (0, 2, 1), (1, 0, 2), (1, 2, 0), (2, 0, 1), (2, 1, 0)
    for perm in permutations:
        assert_raises(AxesError, xyz_affine, cmap3d.reordered_range(perm))
    # The input order doesn't matter, as long as the xyz axes map to the first
    # three input axes
    for perm in permutations:
        assert_array_equal(xyz_affine(
            cmap3d.reordered_domain(perm)), aff3d[:, perm + (-1,)])
    # But if the corresponding input axes not in the first three, an axis error
    wrong_inputs = cmap4d.reordered_domain([0, 1, 3, 2])
    assert_raises(AxesError, xyz_affine, wrong_inputs)
    # xyzs must be orthogonal to dropped axis
    for i in range(3):
        aff = aff4d.copy()
        aff[i,3] = 1
        cmap = AffineTransform(VARS['d_cs_r4'], VARS['r_cs_r4'], aff)
        assert_raises(AffineError, xyz_affine, cmap)
        # And if reordered
        assert_raises(AxesError, xyz_affine, cmap.reordered_range([2,0,1,3]))
    # Non-square goes to square
    aff54 = np.array([[0, 1, 2, 15],
                      [3, 4, 5, 16],
                      [6, 7, 8, 17],
                      [0, 0, 0, 18],
                      [0, 0, 0, 1]])
    cmap = AffineTransform(VARS['d_cs_r3'], VARS['r_cs_r4'], aff54)
    assert_array_equal(xyz_affine(cmap), aff3d)
    aff57 = np.array([[0, 1, 2, 0, 0, 0, 15],
                      [3, 4, 5, 0, 0, 0, 16],
                      [6, 7, 8, 0, 0, 0, 17],
                      [0, 0, 0, 0, 0, 0, 18],
                      [0, 0, 0, 0, 0, 0, 1]])
    d_cs_r6 = CS('ijklmn', 'voxels')
    cmap = AffineTransform(d_cs_r6, VARS['r_cs_r4'], aff57)
    assert_array_equal(xyz_affine(cmap), aff3d)
    # Non-affine raises SpaceTypeError
    cmap_cmap = CoordinateMap(VARS['d_cs_r4'], VARS['r_cs_r4'], lambda x:x*3)
    assert_raises(SpaceTypeError, xyz_affine, cmap_cmap)
    # Not enough dimensions - SpaceTypeError
    d_cs_r2 = CS('ij', 'voxels')
    r_cs_r2 = CS(VARS['r_names'][:2], 'mni')
    cmap = AffineTransform(d_cs_r2, r_cs_r2,
                           np.array([[2,0,10],[0,3,11],[0,0,1]]))
    assert_raises(AxesError, xyz_affine, cmap)
    # Can pass in own validator
    my_valtor = dict(blind='x', leading='y', ditch='z')
    r_cs = CS(('blind', 'leading', 'ditch'), 'fall')
    cmap = AffineTransform(VARS['d_cs_r3'],r_cs, aff3d)
    assert_raises(AxesError, xyz_affine, cmap)
    assert_array_equal(xyz_affine(cmap, my_valtor), aff3d)
    # Slices in x, y, z coordmaps raise error because of missing spatial
    # dimensions
    arr = np.arange(120).reshape((2, 3, 4, 5))
    aff = np.diag([2, 3, 4, 5, 1])
    img = Image(arr, vox2mni(aff))
    assert_raises(AxesError, xyz_affine, img[1].coordmap)
    assert_raises(AxesError, xyz_affine, img[:,1].coordmap)
    assert_raises(AxesError, xyz_affine, img[:,:,1].coordmap)

#.........这里部分代码省略.........

    If ``dim_info`` is set coherently, set input axis names to 'slice', 'freq',
    'phase' from ``dim_info``.

    Get the output spatial coordinate names from the 'scanner', 'aligned',
    'talairach', 'mni' XYZ spaces (see :mod:`nipy.core.reference.spaces`).

    We construct the N-D affine by taking the XYZ affine and adding scaling
    diagonal elements from ``pixdim``.

    If the space units in NIFTI ``xyzt_units`` are 'microns' or 'meters' we
    adjust the affine to mm units, but warn because this might be a mistake.

    If the time units in NIFTI `xyzt_units` are 'msec' or 'usec', scale the time
    axis ``pixdim`` values accordingly.

    Ignore the intent-related fields for now, but warn that we are doing so if
    there appears to be specific information in there.
    """
    hdr = ni_img.get_header()
    affine = ni_img.get_affine()
    # Affine will not be None from a loaded image, but just in case
    if affine is None:
        affine = hdr.get_best_affine()
    else:
        affine = affine.copy()
    data = ni_img.get_data()
    shape = list(ni_img.shape)
    ndim = len(shape)
    if ndim < 3:
        raise NiftiError("With less than 3 dimensions we cannot be sure "
                         "which input and output dimensions you intend for "
                         "the coordinate map.  Please fix this image with "
                         "nibabel or some other tool")
    # For now we only warn if intent is set to an unexpected value
    intent, _, _ = hdr.get_intent()
    if intent != 'none':
        warnings.warn('Ignoring intent field meaning "%s"' % intent,
                        UserWarning)
    # Which space?
    world_label = hdr.get_value_label('sform_code')
    if world_label == 'unknown':
        world_label = hdr.get_value_label('qform_code')
    world_space = XFORM2SPACE.get(world_label, ncrs.unknown_space)
    # Get information from dim_info
    input_names3 = list('ijk')
    freq, phase, slice = hdr.get_dim_info()
    if not freq is None:
        input_names3[freq] = 'freq'
    if not phase is None:
        input_names3[phase] = 'phase'
    if not slice is None:
        input_names3[slice] = 'slice'
    # Add to mm scaling, with warning
    space_units, time_like_units = hdr.get_xyzt_units()
    if space_units in ('micron', 'meter'):
        warnings.warn('"%s" space scaling in NIFTI ``xyt_units field; '
                      'applying scaling to affine, but this may not be what '
                      'you want' % space_units, UserWarning)
        if space_units == 'micron':
            affine[:3] /= 1000.
        elif space_units == 'meter':
            affine[:3] *= 1000.
    input_cs3 = CS(input_names3, name='voxels')
    output_cs3 = world_space.to_coordsys_maker()(3)
    cmap3 = AT(input_cs3, output_cs3, affine)
    if ndim == 3:
        return Image(data, cmap3, {'header': hdr})
    space_units, time_like_units = hdr.get_xyzt_units()
    units_info = TIME_LIKE_UNITS.get(time_like_units, None)
    n_ns = ndim - 3
    ns_zooms = list(hdr.get_zooms()[3:])
    ns_trans = [0] * n_ns
    ns_names = tuple('uvw')
    # Have we got a time axis?
    if (shape[3] == 1 and ndim > 4 and units_info is None):
        # Squeeze length 1 no-time axis
        shape.pop(3)
        ns_zooms.pop(0)
        ns_trans.pop(0)
        data = data.reshape(shape)
        n_ns -= 1
    else: # have time-like
        if units_info is None:
            units_info = TIME_LIKE_UNITS['sec']
        time_name = units_info['name']
        if units_info['scaling'] != 1:
            ns_zooms[0] *= units_info['scaling']
        if time_name == 't':
            # Get time offset
            ns_trans[0] = hdr['toffset']
        ns_names = (time_name,) + ns_names
    output_cs = CS(ns_names[:n_ns])
    input_cs = CS(ns_names[:n_ns])
    aff = from_matvec(np.diag(ns_zooms), ns_trans)
    ns_cmap = AT(input_cs, output_cs, aff)
    cmap = cm_product(cmap3, ns_cmap,
                      input_name=cmap3.function_domain.name,
                      output_name=cmap3.function_range.name)
    return Image(data, cmap, {'header': hdr})

def resample(image, target, mapping, shape, order=3, mode='constant',
             cval=0.0):
    """ Resample `image` to `target` CoordinateMap

    Use a "world-to-world" mapping `mapping` and spline interpolation of a 
    `order`.

    Here, "world-to-world" refers to the fact that mapping should be a
    callable that takes a physical coordinate in "target" and gives a
    physical coordinate in "image".

    Parameters
    ----------
    image : Image instance
       image that is to be resampled.
    target : CoordinateMap
       coordinate map for output image.
    mapping : callable or tuple or array
       transformation from target.function_range to
       image.coordmap.function_range, i.e. 'world-to-world mapping'. Can
       be specified in three ways: a callable, a tuple (A, b)
       representing the mapping y=dot(A,x)+b or a representation of this
       mapping as an affine array, in homogeneous coordinates.
    shape : sequence of int
       shape of output array, in target.function_domain.
    order : int, optional
       what order of interpolation to use in ``scipy.ndimage``.
    mode : str, optional
        Points outside the boundaries of the input are filled according to the
        given mode ('constant', 'nearest', 'reflect' or 'wrap'). Default is
        'constant'.
    cval : scalar, optional
        Value used for points outside the boundaries of the input if
        mode='constant'. Default is 0.0.

    Returns
    -------
    output : Image instance
       Image has interpolated data and output.coordmap == target.
    """
    if not callable(mapping):
        if type(mapping) is type(()):
            mapping = from_matvec(*mapping)
        # image world to target world mapping
        TW2IW = AffineTransform(target.function_range,
                                image.coordmap.function_range,
                                mapping)
    else:
        if isinstance(mapping, AffineTransform):
            TW2IW = mapping
        else:
            TW2IW = CoordinateMap(target.function_range,
                                  image.coordmap.function_range,
                                  mapping)
    # target voxel to image world mapping
    TV2IW = compose(TW2IW, target)
    # CoordinateMap describing mapping from target voxel to
    # image world coordinates
    if not isinstance(TV2IW, AffineTransform):
        # interpolator evaluates image at values image.coordmap.function_range,
        # i.e. physical coordinates rather than voxel coordinates
        grid = ArrayCoordMap.from_shape(TV2IW, shape)
        interp = ImageInterpolator(image, order=order, mode=mode, cval=cval)
        idata = interp.evaluate(grid.transposed_values)
        del(interp)
    else: # it is an affine transform, but, what if we compose?
        TV2IV = compose(image.coordmap.inverse(), TV2IW)
        if isinstance(TV2IV, AffineTransform): # still affine
            A, b = to_matvec(TV2IV.affine)
            idata = affine_transform(image.get_data(), A,
                                     offset=b,