def mask_brain(FeatDir, fMask=""):
    '''
    This function creates a brain mask based on the subject's normalized
    T1_brain image. It can also incorporate an external mask.

    Input parameters:
          FeatDir:    The .feat directory where segmentation results reside.
          fMask:      The file name for the binary mask image. If not provided,
                      then it will be omitted.
    
    Returns:
          NONE:
    
    Output:
          It produces a mask image callsed mask_brain.nii.gz under the
          reg directory of the FeatDir.
    '''
    # directory and file names
    RegDir = os.path.join(FeatDir, 'reg')
    sBase = image_base(FeatDir)
    fT1 = sBase + '.nii.gz'
    fT1_bin = sBase + '_bin.nii.gz'
    ffmri = os.path.join(RegDir, 'func2standard_r.nii.gz')
    fIntersect = os.path.join(RegDir, 'mask_brain.nii.gz')

    # threshold the T1 image
    com_fslmaths = 'fslmaths ' + fT1
    com_fslmaths += ' -bin ' + fT1_bin
    res = os.system(com_fslmaths)

    # reslicing the binarized T1 image to the fMRI space
    fT1_bin_r = reslice_to_fMRI(fT1_bin, ffmri)

    # if an external mask is provided
    if fMask!="":
        # first, copy the mask to the reg directory
        MaskDir, fMaskImg = os.path.split(os.path.abspath(fMask))
        fMaskCopy = os.path.join(RegDir, fMaskImg)
        com_cp = 'cp ' + os.path.join(MaskDir, fMaskImg) + ' ' + RegDir
        res = os.system(com_cp)

        # then reslice the mask copy to the voxel size of fMRI data
        fMaskCopy_r = reslice_to_fMRI(fMaskCopy, ffmri)

        # multiplying the resliced external mask and the brain mask
        # first, load mask images
        img_mask = nib.load(fMaskCopy_r)
        X_mask = img_mask.get_data()
        img_brain = nib.load(fT1_bin_r)
        X_brain = img_brain.get_data()
        # then multiply masks
        X_prod = X_mask * X_brain
        # then saving the new mask image
        maskimg = nib.Nifti1Image(X_prod, img_mask.get_affine())
        nib.save(maskimg, fIntersect)

    # if an external mask is not provided
    else:
        com_cp = 'cp ' + fT1_bin_r + ' ' + fIntersect
        res = os.system(com_cp)

def intersect_masks(input_masks, output_filename=None, threshold=0.5, cc=True):
    """
    Given a list of input mask images, generate the output image which
    is the the threshold-level intersection of the inputs 

    Parameters
    ----------
    input_masks: list of strings or ndarrays
        paths of the input images nsubj set as len(input_mask_files), or
        individual masks.
    output_filename, string:
        Path of the output image, if None no file is saved.
    threshold: float within [0, 1[, optional
        gives the level of the intersection.
        threshold=1 corresponds to keeping the intersection of all
        masks, whereas threshold=0 is the union of all masks.
    cc: bool, optional
        If true, extract the main connected component
        
    Returns
    -------
    grp_mask, boolean array of shape the image shape
    """  
    grp_mask = None
    if threshold > 1:
        raise ValueError('The threshold should be < 1')
    if threshold <0:
        raise ValueError('The threshold should be > 0')
    threshold = min(threshold, 1 - 1.e-7)

    for this_mask in input_masks:
        if isinstance(this_mask, basestring):
            # We have a filename
            this_mask = load(this_mask).get_data()
        if grp_mask is None:
            grp_mask = this_mask.copy().astype(np.int)
        else:
            grp_mask += this_mask
    
    grp_mask = grp_mask > (threshold * len(list(input_masks)))
    
    if np.any(grp_mask > 0) and cc:
        grp_mask = largest_cc(grp_mask)
    
    if output_filename is not None:
        if isinstance(input_masks[0], basestring):
            nim = load(input_masks[0]) 
            header = nim.get_header()
            affine = nim.get_affine()
        else:
            header = dict()
            affine = np.eye(4)
        header['descrip'] = 'mask image'
        output_image = nifti1.Nifti1Image(grp_mask.astype(np.uint8),
                                            affine=affine,
                                            header=header,
                                         )
        save(output_image, output_filename)

    return grp_mask > 0

    def _run_interface(self, runtime):
        from dipy.reconst import dti

        # Load the 4D image files
        img = nb.load(self.inputs.in_file)
        data = img.get_data()
        affine = img.get_affine()

        # Load the gradient strengths and directions
        gtab = self._get_gradient_table()

        # Mask the data so that tensors are not fit for
        # unnecessary voxels
        mask = data[..., 0] > 50

        # Fit the tensors to the data
        tenmodel = dti.TensorModel(gtab)
        tenfit = tenmodel.fit(data, mask)

        # Calculate the mode of each voxel's tensor
        mode_data = tenfit.mode

        # Write as a 3D Nifti image with the original affine
        img = nb.Nifti1Image(mode_data, affine)
        out_file = self._gen_filename('mode')
        nb.save(img, out_file)
        IFLOGGER.info('Tensor mode image saved as {i}'.format(i=out_file))
        return runtime

    def to_image(self, path=None, data=None):
        """Write itself as a binary image, and returns it
        
        Parameters
        ==========
        path: string, path of the output image, if any
        data: array of shape self.size,
              data to put in the nonzer-region of the image
        """
        if data is None:
            wdata = np.zeros(self.shape, np.int8)
        else:
            wdata = np.zeros(self.shape, data.dtype)
        wdata[self.ijk[:, 0], self.ijk[:, 1], self.ijk[:, 2]] = 1
        if data is not None:
            if data.size != self.size:
                raise ValueError('incorrect data size')
            wdata[wdata > 0] = data

        nim = Nifti1Image(wdata, self.affine)
        nim.get_header()['descrip'] = 'mask image representing domain %s' \
            % self.id
        if path is not None:
            save(nim, path)
        return nim

def test_get_vox_dims():
    # setup
    affine = np.eye(4)
    np.fill_diagonal(affine, [-3, 3, 3])

    output_dir = os.path.join(OUTPUT_DIR, inspect.stack()[0][3])
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)

    # 3D vol
    vol = create_random_image(affine=affine)
    np.testing.assert_array_equal(get_vox_dims(vol), [3, 3, 3])

    # 3D image file
    saved_img_filename = os.path.join(output_dir, "vol.nii.gz")
    nibabel.save(vol, saved_img_filename)
    np.testing.assert_array_equal(get_vox_dims(vol), [3, 3, 3])

    # 4D niimg
    film = create_random_image(n_scans=10, affine=affine)
    np.testing.assert_array_equal(get_vox_dims(film), [3, 3, 3])

    # 4D image file
    film = create_random_image(n_scans=10, affine=affine)
    saved_img_filename = os.path.join(output_dir, "4D.nii.gz")
    nibabel.save(film, saved_img_filename)
    np.testing.assert_array_equal(get_vox_dims(film), [3, 3, 3])

def save_vol(vol, output_filename=None, output_dir=None, basename=None,
             concat=False, **kwargs):
    """
    Saves a single volume to disk.

    """

    if not output_filename is None:
        nibabel.save(vol, output_filename)

        return output_filename
    else:
        if output_dir is None:
            raise ValueError(
                'One of output_filename and ouput_dir must be provided')

    if not basename is None:
        if isinstance(basename, list):
            basename = basename[:1]
        else:
            basename = [basename]

    # delegate to legacy save_vols
    return save_vols([vol], output_dir, basenames=basename,
                      concat=False, **kwargs)[0]

def reconst_flow_core(flow, extra_args=[]):
    with TemporaryDirectory() as out_dir:
        data_path, bval_path, bvec_path = get_data('small_25')
        vol_img = nib.load(data_path)
        volume = vol_img.get_data()
        mask = np.ones_like(volume[:, :, :, 0])
        mask_img = nib.Nifti1Image(mask.astype(np.uint8), vol_img.affine)
        mask_path = join(out_dir, 'tmp_mask.nii.gz')
        nib.save(mask_img, mask_path)

        dti_flow = flow()

        args = [data_path, bval_path, bvec_path, mask_path]
        args.extend(extra_args)

        dti_flow.run(*args, out_dir=out_dir)

        fa_path = dti_flow.last_generated_outputs['out_fa']
        fa_data = nib.load(fa_path).get_data()
        assert_equal(fa_data.shape, volume.shape[:-1])

        tensor_path = dti_flow.last_generated_outputs['out_tensor']
        tensor_data = nib.load(tensor_path)
        assert_equal(tensor_data.shape[-1], 6)
        assert_equal(tensor_data.shape[:-1], volume.shape[:-1])

        ga_path = dti_flow.last_generated_outputs['out_ga']
        ga_data = nib.load(ga_path).get_data()
        assert_equal(ga_data.shape, volume.shape[:-1])

        rgb_path = dti_flow.last_generated_outputs['out_rgb']
        rgb_data = nib.load(rgb_path)
        assert_equal(rgb_data.shape[-1], 3)
        assert_equal(rgb_data.shape[:-1], volume.shape[:-1])

        md_path = dti_flow.last_generated_outputs['out_md']
        md_data = nib.load(md_path).get_data()
        assert_equal(md_data.shape, volume.shape[:-1])

        ad_path = dti_flow.last_generated_outputs['out_ad']
        ad_data = nib.load(ad_path).get_data()
        assert_equal(ad_data.shape, volume.shape[:-1])

        rd_path = dti_flow.last_generated_outputs['out_rd']
        rd_data = nib.load(rd_path).get_data()
        assert_equal(rd_data.shape, volume.shape[:-1])

        mode_path = dti_flow.last_generated_outputs['out_mode']
        mode_data = nib.load(mode_path).get_data()
        assert_equal(mode_data.shape, volume.shape[:-1])

        evecs_path = dti_flow.last_generated_outputs['out_evec']
        evecs_data = nib.load(evecs_path).get_data()
        assert_equal(evecs_data.shape[-2:], tuple((3, 3)))
        assert_equal(evecs_data.shape[:-2], volume.shape[:-1])

        evals_path = dti_flow.last_generated_outputs['out_eval']
        evals_data = nib.load(evals_path).get_data()
        assert_equal(evals_data.shape[-1], 3)
        assert_equal(evals_data.shape[:-1], volume.shape[:-1])

def nifti_image_files(outdir, filelist, shape):
    for f in filelist:
        hdr = nb.Nifti1Header()
        hdr.set_data_shape(shape)
        img = np.random.random(shape)
        nb.save(nb.Nifti1Image(img, np.eye(4), hdr),
                 os.path.join(outdir, f))

def labels_to_rd(labels_ct_native_nii, rd_nii, correct_cta, output_dir):
    """ Register the rd to the ct space.
    """

    # Output autocompletion
    labels_rescale_file = os.path.join(output_dir, "BrainMask_to_rd.nii.gz")

    # Load images
    labels_ct_native_im = nibabel.load(labels_ct_native_nii)
    labels_data = labels_ct_native_im.get_data()
    rd_im = nibabel.load(rd_nii)
    rd_data = rd_im.get_data()
    cta = labels_ct_native_im.get_affine()
    rda = rd_im.get_affine()

    # Correct the rda affine matrix
    print cta	
    #cta[1, 3] += 21
    cta[2, 2] = correct_cta
    print cta


    # Inverse affine transformation
    icta = inverse_affine(cta)
    t = numpy.dot(icta, rda)
    numpy.savetxt(os.path.join(output_dir,"t_icta_rda.txt"), t)
    # Matricial dot product
    labels_rescale = numpy.zeros(rd_data.shape)
    dot_image = numpy.zeros(rd_data.shape + (3, ))
    x = numpy.linspace(0, rd_data.shape[0] - 1, rd_data.shape[0])
    y = numpy.linspace(0, rd_data.shape[1] - 1, rd_data.shape[1])
    z = numpy.linspace(0, rd_data.shape[2] - 1, rd_data.shape[2])
    xg, yg, zg = numpy.meshgrid(y, x, z)
    print 'dot image shape: ', dot_image.shape
    print 'yg shape: ', yg.shape
    print 'xg shape: ', xg.shape
    print 'zg shape: ', zg.shape
    dot_image[..., 0] = yg
    dot_image[..., 1] = xg
    dot_image[..., 2] = zg
    dot_image = threed_dot(t, dot_image)

    cnt = 0
    print rd_data.size
    for x in range(rd_data.shape[0]):
        for y in range(rd_data.shape[1]):
            for z in range(rd_data.shape[2]):
            #for z in range(cut_brain_index, rd_data.shape[2]):
                if cnt % 100000 == 0:
                    print cnt
                cnt += 1
                voxel_labels = dot_image[x, y, z]
                if (voxel_labels > 0).all() and (voxel_labels < (numpy.asarray(labels_data.shape) - 1)).all():
                    voxel_labels = numpy.round(voxel_labels)
                    labels_rescale[x, y, z] = labels_data[voxel_labels[0], voxel_labels[1], voxel_labels[2]]

    labels_rescale_im = nibabel.Nifti1Image(labels_rescale, rda)
    nibabel.save(labels_rescale_im, labels_rescale_file)

    return labels_rescale_file

    def _run_interface(self, runtime):
        nii1 = nb.load(self.inputs.volume1)
        nii2 = nb.load(self.inputs.volume2)

        origdata1 = np.logical_not(np.logical_or(nii1.get_data() == 0, np.isnan(nii1.get_data())))
        origdata2 = np.logical_not(np.logical_or(nii2.get_data() == 0, np.isnan(nii2.get_data())))

        if isdefined(self.inputs.mask_volume):
            maskdata = nb.load(self.inputs.mask_volume).get_data()
            maskdata = np.logical_not(np.logical_or(maskdata == 0, np.isnan(maskdata)))
            origdata1 = np.logical_and(maskdata, origdata1)
            origdata2 = np.logical_and(maskdata, origdata2)

        for method in ("dice", "jaccard"):
            setattr(self, "_" + method, self._bool_vec_dissimilarity(origdata1, origdata2, method=method))

        self._volume = int(origdata1.sum() - origdata2.sum())

        both_data = np.zeros(origdata1.shape)
        both_data[origdata1] = 1
        both_data[origdata2] += 2

        nb.save(nb.Nifti1Image(both_data, nii1.get_affine(), nii1.get_header()), self.inputs.out_file)

        return runtime

def reorient_image(input_axes, in_file, output_dir):
    """ Rectify the orientation of an image.
    """
    # get the transformation to the RAS space
    rotation = swap_affine(input_axes)
    det = numpy.linalg.det(rotation)
    if det != 1:
        raise Exception("Determinant must be equal to "
                        "one got: {0}.".format(det))

    # load image
    image = nibabel.load(in_file)

    # get affine transform (qform or sform)
    affine = image.get_affine()

    # apply transformation
    transformation = numpy.dot(rotation, affine)
    image.set_qform(transformation)
    image.set_sform(transformation)

    # save result
    reoriented_file = os.path.join(output_dir, "im_reorient.nii.gz")
    nibabel.save(image, reoriented_file)

    return reoriented_file

def data_for_neuronal_network(seg_src_path, seg_dest_path):

    ih.createPathIfNotExists(seg_dest_path)

    seg_list = os.listdir(seg_src_path) # seg_list: list of all files in segmentation directory

    print ("List of input files:")
    print (seg_list)

    file_num = 0

    for f in seg_list:
        print("file number: " + str(file_num))
        print("read file...", f)

        path = seg_src_path + "/" + f
        print(path)

        image = nib.load(path)
        headers = image.header
        # print(headers)

        matrix = image.affine
       # matrix[np.where(matrix > 1)] = 1
       # matrix[np.where(matrix < 1)] = 0
        print(matrix)

        # manipulations.....
        #numbers

        dest_file = seg_dest_path + "/" + f
        nib.save(image, dest_file)

def save2nifti(data, affine, file_name):
    """
    Save data to a nifti file.

    """
    img = nib.Nifti1Image(data, affine)
    nib.save(img, file_name)

def saveImageToANewNiiWithHeaderFromOtherGivenExactFilePaths(labelImageCreatedByPredictions,
                                          fullFilenameToSaveWith,
                                          fullFilenameOfOriginalImageToCopyHeader,
                                          npDtype = np.dtype(np.float32),
                                          myLogger=None) :

    fullFilenameToSaveWith = os.path.abspath(fullFilenameToSaveWith) # Cleans the .././/...
    img_proxy_for_orig_image = nib.load(fullFilenameOfOriginalImageToCopyHeader)
    hdr_for_orig_image = img_proxy_for_orig_image.header
    
    affine_trans_to_ras = img_proxy_for_orig_image.affine
    
    newLabelImg = nib.Nifti1Image(labelImageCreatedByPredictions, affine_trans_to_ras) #Nifti Constructor. data is the image itself, dimensions x,y,z,time. The second argument is the affine RAS transf.
    newLabelImg.set_data_dtype(npDtype)

    dimensionsOfTheGivenArrayImageToSave = len(labelImageCreatedByPredictions.shape)
    newZooms = list(hdr_for_orig_image.get_zooms()[:dimensionsOfTheGivenArrayImageToSave])
    if len(newZooms) < dimensionsOfTheGivenArrayImageToSave : #Eg if original image was 3D, but I need to save a multichannel image.
	newZooms = newZooms + [1.0]*(dimensionsOfTheGivenArrayImageToSave - len(newZooms))
    newLabelImg.header.set_zooms(newZooms)

    if not fullFilenameToSaveWith.endswith(".nii.gz") :
            fullFilenameToSaveWith = fullFilenameToSaveWith + ".nii.gz"
    nib.save(newLabelImg, fullFilenameToSaveWith)

    if myLogger<>None :
        myLogger.print3("Image saved at: " + str(fullFilenameToSaveWith))
    else :
	print("Image saved at: " + str(fullFilenameToSaveWith)) 

def fill_nan(in_file, fill_value=0.):
    """Replace nan values with a given value

    Parameters
    ----------
    in_file : str
        Path to image file

    fill_value : float, optional
        Value replacing nan

    Returns
    -------
    out_file : str
        Path to the filled file
    """
    img = nibabel.load(in_file)
    data = img.get_data()
    if np.any(np.isnan(data)):
        data[np.isnan(data)] = fill_value
    img = nibabel.Nifti1Image(data, img.get_affine(), img.get_header())
    out_file, _ = os.path.splitext(in_file)
    out_file += '_no_nan.nii'
    nibabel.save(img, out_file)
    return out_file

    def export(self, nidm_version, export_dir):
        """
        Create prov graph.
        """
        if self.expl_mean_sq_file is None:
            # Create Contrast Explained Mean Square Map as fstat<num>.nii.gz
            # multiplied by sigmasquareds.nii.gz and save it in export_dir
            fstat_img = nib.load(self.stat_file)
            fstat = fstat_img.get_data()

            sigma_sq_img = nib.load(self.sigma_sq_file)
            sigma_sq = sigma_sq_img.get_data()

            expl_mean_sq = nib.Nifti1Image(
                fstat*sigma_sq, fstat_img.get_qform())

            self.filename = ("ContrastExplainedMeanSquareMap" +
                             self.num + ".nii.gz")
            self.expl_mean_sq_file = os.path.join(
                export_dir, self.filename)
            nib.save(expl_mean_sq, self.expl_mean_sq_file)

        self.file = NIDMFile(self.id, self.expl_mean_sq_file,
                             filename=self.filename,
                             sha=self.sha, fmt=self.fmt)

        # Contrast Explained Mean Square Map entity
        path, filename = os.path.split(self.expl_mean_sq_file)
        self.add_attributes((
            (PROV['type'], self.type),
            (NIDM_IN_COORDINATE_SPACE, self.coord_space.id),
            (PROV['label'], self.label)))

def sphere(file_roi, vol, mm_x=0, mm_y=0, mm_z=0, radius=4, dtype=np.uint8):
    img = nibabel.load(vol)
    header = img.get_header()
    voxels_size= header.get_zooms()
    coord_center_voxel = mm_to_voxel(vol, [[mm_x, mm_y, mm_z]])    
    shape_x, shape_y, shape_z = img.shape[0], img.shape[1], img.shape[2]
    roi = np.zeros((shape_x, shape_y, shape_z))
    
    radius = radius / voxels_size[0]
    if radius > 1 :
       radius -= 1 

    first_x = coord_center_voxel[0][0] - radius
    first_y = coord_center_voxel[0][1] - radius
    first_z = coord_center_voxel[0][2] - radius

    elem = morphology.ball(radius)
    
    #check voxel size is isotropic
    #check radius is at least one voxel
    for x in range(0,elem.shape[0]):
            for y in range(0, elem.shape[1]):
                    for z in range(0, elem.shape[2]):
                        x_ = int(first_x) + x
                        y_ = int(first_y) + y
                        z_ = int(first_z) + z
                        roi[x_, y_, z_] = elem[x,y,z]
  
    roi_img = nibabel.Nifti1Image(roi, img.get_affine(),img.get_header() )
    nibabel.save(roi_img, file_roi)

	def _run_interface(self, runtime):
		tracks, header = trk.read(self.inputs.in_file)
		if not isdefined(self.inputs.data_dims):
			data_dims = header['dim']
		else:
			data_dims = self.inputs.data_dims

		if not isdefined(self.inputs.voxel_dims):
			voxel_size = header['voxel_size']
		else:
			voxel_size = self.inputs.voxel_dims

		affine = header['vox_to_ras']

		streams = ((ii[0]) for ii in tracks)
		data = density_map(streams, data_dims, voxel_size)
		if data.max() < 2**15:
		   data = data.astype('int16')

		img = nb.Nifti1Image(data,affine)
		out_file = op.abspath(self.inputs.out_filename)
		nb.save(img, out_file)
		iflogger.info('Track density map saved as {i}'.format(i=out_file))
		iflogger.info('Data Dimensions {d}'.format(d=data_dims))
		iflogger.info('Voxel Dimensions {v}'.format(v=voxel_size))
		return runtime

	def _run_interface(self, runtime):
		import dipy.reconst.dti as dti
		import dipy.denoise.noise_estimate as ne
		from dipy.core.gradients import gradient_table
		from nipype.utils.filemanip import split_filename

		import nibabel as nib

		fname = self.inputs.in_file
		img = nib.load(fname)
		data = img.get_data()
		affine = img.get_affine()

		bvals = self.inputs.bval
		bvecs = self.inputs.bvec

		gtab = gradient_table(bvals, bvecs)
		sigma = ne.estimate_sigma(data)
		dti = dti.TensorModel(gtab,fit_method='RESTORE', sigma=sigma)
		dtifit = dti.fit(data)
		fa = dtifit.fa

		_, base, _ = split_filename(fname)
		nib.save(nib.Nifti1Image(fa, affine), base + '_FA.nii')

		return runtime

    def coefs_to_niftii(self, coefs, nii_name='coefs.nii.gz'):
        """

        """

        print(coefs)
        if self.n_features != len(coefs):
            print("Num features {} does not match coef size {}. Did you downsize?".format(
                self.n_features, len(coefs)))
            raise(ValueError)

        restored_coefs = self.restore_coefs(coefs)
        print("restoring... num non-zero",np.count_nonzero(restored_coefs) )

        try:
            restored_coefs = restored_coefs.reshape(np.array(self.nii_shape))
            print("Restored shape:", restored_coefs.shape)
        except IndexError:
            raise(IndexError, "Coefs don't match niftii shape")



        sample_img = nib.load('tt29.nii')
        affine = sample_img.affine.copy()
        header = sample_img.header.copy()

        header.set_data_dtype(np.float64)
        assert(header.get_data_dtype()==np.float64)

        print("Saving coefs as datatype", header.get_data_dtype())
        i = nib.Nifti1Image(restored_coefs.astype(np.float64), affine, header)
        nib.save(i, nii_name)

        return restored_coefs