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 save_niigz(file_path, vol, affine=None, header=None):
    """Saves a volume into a Nifti (.nii.gz) file.

    Parameters
    ----------
    vol: Numpy 3D or 4D array
        Volume with the data to be saved.
    file_path: string
        Output file name path
    affine: 4x4 Numpy array
        Array with the affine transform of the file.
    header: nibabel.nifti1.Nifti1Header, optional
        Header for the file, optional but recommended.

    Note
    ----
        affine and header only work for numpy volumes.

    """
    if isinstance(vol, np.ndarray):
        log.debug('Saving numpy nifti file: ' + file_path)
        ni = nib.Nifti1Image(vol, affine, header)
        nib.save(ni, file_path)

    elif isinstance(vol, nib.Nifti1Image):
        log.debug('Saving nibabel nifti file: ' + file_path)
        nib.save(vol, file_path)

    elif isinstance(vol, niim.Image):
        log.debug('Saving nibabel nifti file: ' + file_path)
        save_image(vol, file_path)

def space_time_realign(Images,TR=2,numslices=None,SliceTime='asc_alt_2',RefScan=None):
    '''
    4D simultaneous slice timing and spatial realignment. Adapted from
    Alexis Roche's example script, and extend to be used for multiplex
    imaging sequences
    
    Inputs:
    
        Images: list of images, input as a list of strings
        
        numslices: for non-multiplex sequence, default to be the number of
            slices in the image. For multiplex sequence, enter as a tuple,
            such that the first element is the number of planes acquired in
            parallel between each other, and the second element is the number
            of slices of each parallel plane/slab
        
        SliceTime:enter as a string to specify how the slices are ordered.
            Choices are the following
            1).'ascending': sequential ascending acquisition
            2).'descending': sequential descending acquisition
            3).'asc_alt_2': ascending interleaved, starting at first slice
            4).'asc_alt_2_1': ascending interleaved, starting at the second
                slice
            5).'desc_alt_2': descending interleaved, starting at last slice
            6).'asc_alt_siemens': ascending interleaved, starting at the first
                slice if odd number of slices, or second slice if even number
                of slices
            7).'asc_alt_half': ascending interleaved by half the volume
            8).'desc_alt_half': descending interleaved by half the volume
        
        RefScan: reference volume for spatial realignment movement estimation
    '''
    
    # load images    
    runs = [load_image(run) for run in Images]
    # parse data info
    if numslices is None:
        numslices = runs[0].shape[2]
        numplanes = 1
    elif isinstance(numslices,tuple):
        numslices = numslices[0]
        numplanes = numplanes[1]
    # parse slice timing according to the input
    slice_timing = getattr(timefuncs,SliceTime)(TR,numslices)
    #repeat the slice timing for multiplex seqquence
    slice_timing = np.tile(slice_timing,numplanes)
    # Spatio-temporal realigner assuming interleaved ascending slice order
    R = SpaceTimeRealign(runs, tr=TR, slice_times=slice_timing, slice_info=2,
                         affine_class='Rigid')
    
    print('Slice times: %s' % slice_timing)
    # Estimate motion within- and between-sessions
    R.estimate(refscan=RefScan)
    # Resample data on a regular space+time lattice using 4d interpolation
    print('Saving results ...')
    for i in range(len(runs)):
        corr_run = R.resample(i)
        fname = os.path.join(os.path.split(Images[i])[0],'ra' + os.path.split(Images[i])[1])
        save_image(corr_run, fname)
        print(fname)

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

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

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

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

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

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

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

    return resampled_file

def sample_map_allen_space(image_path, annot_csv_path, save_path, type='well_id'):
    #assign values to samples in MNI space

    I=nipy.load_image(image_path)
    image_name=os.path.basename(image_path)
    df=pd.DataFrame.from_csv(annot_csv_path)
    coordinate, well_id=(np.array( df['mri_voxel_x']) , np.array(df['mri_voxel_y']), np.array(df['mri_voxel_z'] )), df[type]
    I._data[np.where(I._data!=0)]=0
    I._data[coordinate]=well_id
    nipy.save_image(I, os.path.join(save_path, image_name))

def tsdiffana(args):
    """ Generate tsdiffana plots from command line params `args`

    Parameters
    ----------
    args : object
        object with attributes

        * filename : str - 4D image filename
        * out_file : str - graphics file to write to instead of leaving
          graphics on screen
        * time_axis : str - name or number of time axis in `filename`
        * slice_axis : str - name or number of slice axis in `filename`
        * write_results : bool - if True, write images and plots to files
        * out_path : None or str - path to which to write results
        * out_fname_label : None or filename - suffix of output results files

    Returns
    -------
    axes : Matplotlib axes
       Axes on which we have done the plots.
    """
    if args.out_file is not None and args.write_results:
        raise ValueError("Cannot have OUT_FILE and WRITE_RESULTS options "
                         "together")
    img, time_axis, slice_axis = parse_fname_axes(args.filename,
                                                  args.time_axis,
                                                  args.slice_axis)
    results = time_slice_diffs_image(img, time_axis, slice_axis)
    axes = plot_tsdiffs(results)
    if args.out_file is None and not args.write_results:
        # interactive mode
        return axes
    if args.out_file is not None:
        # plot only mode
        axes[0].figure.savefig(args.out_file)
        return axes
    # plot and images mode
    froot, ext, addext = splitext_addext(args.filename)
    fpath, fbase = psplit(froot)
    fpath = fpath if args.out_path is None else args.out_path
    fbase = fbase if args.out_fname_label is None else args.out_fname_label
    axes[0].figure.savefig(pjoin(fpath, 'tsdiff_' + fbase + '.png'))
    # Save image volumes
    for key, prefix in (('slice_diff2_max_vol', 'dv2_max_'),
                        ('diff2_mean_vol', 'dv2_mean_')):
        fname = pjoin(fpath, prefix + fbase + ext + addext)
        nipy.save_image(results[key], fname)
    # Save time courses into npz
    np.savez(pjoin(fpath, 'tsdiff_' + fbase + '.npz'),
             volume_means=results['volume_means'],
             slice_mean_diff2=results['slice_mean_diff2'],
            )
    return axes

def save_niigz(filepath, vol, header=None, affine=None):
    """Saves a volume into a Nifti (.nii.gz) file.

    Parameters
    ----------
    vol: Numpy 3D or 4D array
        Volume with the data to be saved.

    file_path: string
        Output file name path

    affine: (optional) 4x4 Numpy array
        Array with the affine transform of the file.
        This is needed if vol is a np.ndarray.

    header: (optional) nibabel.nifti1.Nifti1Header, optional
        Header for the file, optional but recommended.
        This is needed if vol is a np.ndarray.

    Note
    ----
    affine and header only work for numpy volumes.
    """
    # delayed import because could not install nipy on Python 3 on OSX
    we_have_nipy = False
    try:
        import nipy.core.image as     niim
        from   nipy            import save_image
    except:
        pass
    else:
        we_have_nipy = True

    if isinstance(vol, np.ndarray):
        log.debug('Saving numpy nifti file: {}.'.format(filepath))
        ni = nib.Nifti1Image(vol, affine, header)
        nib.save(ni, filepath)

    elif isinstance(vol, nib.Nifti1Image):
        log.debug('Saving nibabel nifti file: {}.'.format(filepath))
        nib.save(vol, filepath)

    elif we_have_nipy and isinstance(vol, niim.Image):
        log.debug('Saving nipy nifti file: {}.'.format(filepath))
        save_image(vol, filepath)

    #elif isinstance(vol, NeuroImage):
    #    log.debug('Saving boyle.NeuroImage nifti file: {}.'.format(filepath))
    #    nib.save(vol.img, filepath)

    else:
        raise ValueError('Could not recognise input vol filetype. Got: {}.'.format(repr_imgs(vol)))

def peelTemplateBrain():
    ns=181
    nr=217
    nc=181
    gt_template=np.fromfile('data/phantom_1.0mm_normal_crisp.rawb', dtype=np.ubyte).reshape((ns,nr,nc))
    t1_template=np.fromfile('data/t1/t1_icbm_normal_1mm_pn0_rf0.rawb', dtype=np.ubyte).reshape((ns,nr,nc))
    t2_template=np.fromfile('data/t2/t2_icbm_normal_1mm_pn0_rf0.rawb', dtype=np.ubyte).reshape((ns,nr,nc))
    #t1_template*=((1<=gt_template)*(gt_template<=3)+(gt_template==8))
    t1_template*=((1<=gt_template)*(gt_template<=3))
    t2_template*=((1<=gt_template)*(gt_template<=3))
    affine_transform=AffineTransform('ijk', ['aligned-z=I->S','aligned-y=P->A', 'aligned-x=L->R'], np.eye(4))
    t1_template=Image(t1_template, affine_transform)
    t2_template=Image(t2_template, affine_transform)
    nipy.save_image(t1_template,'data/t1/t1_icbm_normal_1mm_pn0_rf0_peeled.nii.gz')
    nipy.save_image(t2_template,'data/t2/t2_icbm_normal_1mm_pn0_rf0_peeled.nii.gz')

def time_space_realign(run_fnames, TR, time_to_space, slice_axis):
    run_imgs = [load_image(run) for run in run_fnames]
    # Spatio-temporal realigner
    R = FmriRealign4d(run_imgs,
                      tr=TR,
                      slice_order=time_to_space,
                      slice_info=(slice_axis, 1))
    # Estimate motion within- and between-sessions
    R.estimate(refscan=None)
    # Save back out
    for i, fname in enumerate(run_fnames):
        corr_run = R.resample(i)
        pth, name = os.path.split(fname)
        processed_fname = os.path.join(pth, 'ra' + name)
        save_image(corr_run, processed_fname)

def generateTestingPair(betaGT):
    betaGTRads=np.array(betaGT, dtype=np.float64)
    betaGTRads[0:3]=np.copy(np.pi*betaGTRads[0:3]/180.0)
    ns=181
    nr=217
    nc=181
    left=np.fromfile('data/t2/t2_icbm_normal_1mm_pn0_rf0.rawb', dtype=np.ubyte).reshape(ns,nr,nc)
    left=left.astype(np.float64)
    right=np.fromfile('data/t1/t1_icbm_normal_1mm_pn0_rf0.rawb', dtype=np.ubyte).reshape(ns,nr,nc)
    right=right.astype(np.float64)
    right=rcommon.applyRigidTransformation3D(right, betaGTRads)
    affine_transform=AffineTransform('ijk', ['aligned-z=I->S','aligned-y=P->A', 'aligned-x=L->R'], np.eye(4))
    left=Image(left, affine_transform)
    right=Image(right, affine_transform)
    nipy.save_image(left,'moving.nii')
    nipy.save_image(right,'fixed.nii')

    def save_niftis(self, X):
        base_nifti = nipy.load_image(self.base_nifti_file)

        if self.pca is not None and self.pca_components:
            X = self.pca.inverse_transform(X)

        images = []
        out_files = []
        for i, x in enumerate(X):
            image = self.make_image(x, base_nifti, do_pca=False)
            out_file = path.join(self.tmp_path, 'tmp_image_%d.nii.gz' % i)
            nipy.save_image(image, out_file)
            images.append(image)
            out_files.append(out_file)

        return images, out_files

def main():
    try:
        DATA_PATH = sys.argv[1]
    except IndexError:
        raise RuntimeError("Pass data path on command line")
    subjects = get_subjects(DATA_PATH)
    for name in sorted(subjects):
        subject = subjects[name]
        print("Smoothing subject " + name)
        for run in subject['functionals']:
            fname = run['filename']
            pth, fpart = os.path.split(fname)
            ra_fname = os.path.join(pth, 'ra' + fpart)
            sra_fname = os.path.join(pth, 'sra' + fpart)
            img = load_image(ra_fname)
            save_image(smooth_image(img, 8.), sra_fname)

    def _run_interface(self, runtime):

        all_ims = []

        for image in self.inputs.in_file:
            im = nb.load(image)
            im.affine = im.get_affine()
            all_ims.append(im)

        if not isdefined(self.inputs.tr_slices):
            TR_slices = None
        else:
            TR_slices = self.inputs.tr_slices

        R = FR4d(all_ims, tr=self.inputs.tr,
            slice_order=self.inputs.slice_order,
            interleaved=self.inputs.interleaved,
            tr_slices=TR_slices,
            time_interp=self.inputs.time_interp,
            start=self.inputs.start)

        R.estimate(loops=self.inputs.loops,
                   between_loops=self.inputs.between_loops,
                   speedup=self.inputs.speedup)

        corr_run = R.resample()
        self._out_file_path = []
        self._par_file_path = []

        for j, corr in enumerate(corr_run):
            self._out_file_path.append(os.path.abspath('corr_%s.nii.gz' %
            (split_filename(self.inputs.in_file[j])[1])))
            save_image(corr, self._out_file_path[j])

            self._par_file_path.append(os.path.abspath('%s.par' %
            (os.path.split(self.inputs.in_file[j])[1])))
            mfile = open(self._par_file_path[j], 'w')
            motion = R._transforms[j]
            #output a .par file that looks like fsl.mcflirt's .par file
            for i, mo in enumerate(motion):
                params = ['%.10f' % item for item in np.hstack((mo.rotation,
                                                             mo.translation))]
                string = ' '.join(params) + '\n'
                mfile.write(string)
            mfile.close()

        return runtime

    def _run_interface(self, runtime):
        from nipy import save_image, load_image
        all_ims = [load_image(fname) for fname in self.inputs.in_file]

        if not isdefined(self.inputs.slice_times):
            from nipy.algorithms.registration.groupwise_registration import \
                SpaceRealign
            R = SpaceRealign(all_ims)
        else:
            from nipy.algorithms.registration import SpaceTimeRealign
            R = SpaceTimeRealign(
                all_ims,
                tr=self.inputs.tr,
                slice_times=self.inputs.slice_times,
                slice_info=self.inputs.slice_info,
            )

        R.estimate(refscan=None)

        corr_run = R.resample()
        self._out_file_path = []
        self._par_file_path = []

        for j, corr in enumerate(corr_run):
            self._out_file_path.append(
                os.path.abspath('corr_%s.nii.gz' %
                                (split_filename(self.inputs.in_file[j])[1])))
            save_image(corr, self._out_file_path[j])

            self._par_file_path.append(
                os.path.abspath('%s.par' %
                                (os.path.split(self.inputs.in_file[j])[1])))
            mfile = open(self._par_file_path[j], 'w')
            motion = R._transforms[j]
            # nipy does not encode euler angles. return in original form of
            # translation followed by rotation vector see:
            # http://en.wikipedia.org/wiki/Rodrigues'_rotation_formula
            for i, mo in enumerate(motion):
                params = [
                    '%.10f' % item
                    for item in np.hstack((mo.translation, mo.rotation))
                ]
                string = ' '.join(params) + '\n'
                mfile.write(string)
            mfile.close()

        return runtime

    def _run_interface(self, runtime):
        from nipy.algorithms.registration import FmriRealign4d as FR4d

        all_ims = [load_image(fname) for fname in self.inputs.in_file]

        if not isdefined(self.inputs.tr_slices):
            TR_slices = None
        else:
            TR_slices = self.inputs.tr_slices

        R = FR4d(
            all_ims,
            tr=self.inputs.tr,
            slice_order=self.inputs.slice_order,
            tr_slices=TR_slices,
            time_interp=self.inputs.time_interp,
            start=self.inputs.start,
        )

        R.estimate(
            loops=list(self.inputs.loops),
            between_loops=list(self.inputs.between_loops),
            speedup=list(self.inputs.speedup),
        )

        corr_run = R.resample()
        self._out_file_path = []
        self._par_file_path = []

        for j, corr in enumerate(corr_run):
            self._out_file_path.append(os.path.abspath("corr_%s.nii.gz" % (split_filename(self.inputs.in_file[j])[1])))
            save_image(corr, self._out_file_path[j])

            self._par_file_path.append(os.path.abspath("%s.par" % (os.path.split(self.inputs.in_file[j])[1])))
            mfile = open(self._par_file_path[j], "w")
            motion = R._transforms[j]
            # nipy does not encode euler angles. return in original form of
            # translation followed by rotation vector see:
            # http://en.wikipedia.org/wiki/Rodrigues'_rotation_formula
            for i, mo in enumerate(motion):
                params = ["%.10f" % item for item in np.hstack((mo.translation, mo.rotation))]
                string = " ".join(params) + "\n"
                mfile.write(string)
            mfile.close()

        return runtime

def tissue_classification(img,mask=None,niters=25,beta=0.5,ngb_size=6,probc=None,probg=None,probw=None):
    import numpy as np

    from nipy import load_image, save_image
    from nipy.core.image.image_spaces import (make_xyz_image,
                                              xyz_affine)
    from nipy.algorithms.segmentation import BrainT1Segmentation
    import os
    # Input image
    img = load_image(img)

    # Input mask image
    mask_img = mask
    if mask_img == None:
        mask_img = img
    else:
        mask_img = load_image(mask_img)

    # Other optional arguments
    #niters = int(get_argument('niters', 25))
    #beta = float(get_argument('beta', 0.5))
    #ngb_size = int(get_argument('ngb_size', 6))

    # Perform tissue classification
    mask = mask_img.get_data() > 0
    S = BrainT1Segmentation(img.get_data(), mask=mask, model='5k',
        niters=niters, beta=beta, ngb_size=ngb_size)

    # Save label image
    outfile = os.path.abspath('hard_classif.nii')
    save_image(make_xyz_image(S.label, xyz_affine(img), 'scanner'),
        outfile)
    print('Label image saved in: %s' % outfile)

    # Compute fuzzy Dice indices if a 3-class fuzzy model is provided
    if not probc == None and\
       not probg == None and\
       not probw == None:
        print('Computing Dice index')
        gold_ppm = np.zeros(S.ppm.shape)
        gold_ppm_img = (probc, probg, probw)
        for k in range(3):
            img = load_image(gold_ppm_img[k])
            gold_ppm[..., k] = img.get_data()
        d = fuzzy_dice(gold_ppm, S.ppm, np.where(mask_img.get_data() > 0))
        print('Fuzzy Dice indices: %s' % d)

def segment_file(input_filename, output_filename,
                 model_name, threshold, verbosity,
                 use_tta):
    """Segment a volume file.

    :param input_filename: the input filename.
    :param output_filename: the output filename.
    :param model_name: the name of model to use.
    :param threshold: threshold to apply in predictions (if None,
                      no threshold will be applied)
    :param verbosity: the verbosity level.
    :param use_tta: whether it should use TTA (test-time augmentation)
                    or not.
    :return: the output filename.
    """
    nii_original = nipy.load_image(input_filename)
    pixdim = nii_original.header["pixdim"][3]
    target_resample = "0.25x0.25x{:.5f}".format(pixdim)

    nii_resampled = resampling.resample_nipy(nii_original, new_size=target_resample, new_size_type='mm',
                                             interpolation='linear', verbose=verbosity)
    nii_resampled = nipy2nifti(nii_resampled)
    pred_slices = segment_volume(nii_resampled, model_name, threshold,
                                 use_tta)

    original_res = "{:.5f}x{:.5f}x{:.5f}".format(
        nii_original.header["pixdim"][1],
        nii_original.header["pixdim"][2],
        nii_original.header["pixdim"][3])

    volume_affine = nii_resampled.affine
    volume_header = nii_resampled.header
    nii_segmentation = nib.Nifti1Image(pred_slices, volume_affine,
                                       volume_header)
    nii_segmentation = nifti2nipy(nii_segmentation)

    nii_resampled_original = resampling.resample_nipy(nii_segmentation, new_size=original_res, new_size_type='mm',
                                                      interpolation='linear', verbose=verbosity)
    res_data = nii_resampled_original.get_data()

    # Threshold after resampling, only if specified
    if threshold is not None:
        res_data = threshold_predictions(res_data, 0.5)

    nipy.save_image(nii_resampled_original, output_filename)
    return output_filename

def tissue_segmentation(save_path,segmented_image, type, info):

    if info=='freesurfer' and type=="GM":
        # FreeSurfer gray matter regions
        # SupraGMV=LeftCerebralCortex+RightCerebralCortex+SubcorticalGMV
        region=np.array([42,3,10,11,12,13,17,18,26,49,50,51,52,53,54,58])


        S=nipy.load_image(segmented_image)
        name="GM_"+ os.path.basename(segmented_image)
        index=np.in1d(S._data, region).reshape(S._data.shape) # coordinates of voxels with values = values from region array
        S._data[index]=1 # set voxels values to 1
        index=np.logical_not(index) #get coordinates of voxels with values not = values from region array
        S._data[index]=0 # set voxels values to 0
        nipy.save_image(S, os.path.join(save_path,name))
    else:
        raise ValueError('Not implemented!')

def smooth_mask_nipy(infile, outfile, fwhm=14):
    """uses nipy to smooth an image using gaussian filter of fwhm"""
    img = nipy.load_image(infile)
    lf = nipy.algorithms.kernel_smooth.LinearFilter(img.coordmap,
                                                    img.shape,
                                                    fwhm)
    simg = lf.smooth(img)
    outimg = nipy.save_image(simg, outfile)
    return outimg

def test_nipy_diagnose():
    # Test nipy diagnose script
    fimg = load_image(funcfile)
    ncomps = 12
    with InTemporaryDirectory() as tmpdir:
        cmd = ['nipy_diagnose', funcfile,
               '--ncomponents={0}'.format(ncomps),
               '--out-path=' + tmpdir]
        run_command(cmd)
        for out_fname in ('components_functional.png',
                          'pcnt_var_functional.png',
                          'tsdiff_functional.png',
                          'vectors_components_functional.npz'):
            assert_true(isfile(out_fname))
        for out_img in ('max_functional.nii.gz',
                        'mean_functional.nii.gz',
                        'min_functional.nii.gz',
                        'std_functional.nii.gz'):
            img = load_image(out_img)
            assert_equal(img.shape, fimg.shape[:-1])
            del img
        pca_img = load_image('pca_functional.nii.gz')
        assert_equal(pca_img.shape, fimg.shape[:-1] + (ncomps,))
        vecs_comps = np.load('vectors_components_functional.npz')
        vec_diff = vecs_comps['slice_mean_diff2'].copy()# just in case
        assert_equal(vec_diff.shape, (fimg.shape[-1]-1, fimg.shape[2]))
        del pca_img, vecs_comps
    with InTemporaryDirectory() as tmpdir:
        # Check we can pass in slice and time flags
        s0_img = rollimg(fimg, 'k')
        save_image(s0_img, 'slice0.nii')
        cmd = ['nipy_diagnose', 'slice0.nii',
               '--ncomponents={0}'.format(ncomps),
               '--out-path=' + tmpdir,
               '--time-axis=t',
               '--slice-axis=0']
        run_command(cmd)
        pca_img = load_image('pca_slice0.nii')
        assert_equal(pca_img.shape, s0_img.shape[:-1] + (ncomps,))
        vecs_comps = np.load('vectors_components_slice0.npz')
        assert_almost_equal(vecs_comps['slice_mean_diff2'], vec_diff)
        del pca_img, vecs_comps