Python affines.from_matvec函数代码示例

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)