def spikes_mask(in_file, in_mask=None, out_file=None):
    """
    Utility function to calculate a mask in which check
    for :abbr:`EM (electromagnetic)` spikes.
    """

    import os.path as op
    import nibabel as nb
    import numpy as np
    from nilearn.image import mean_img
    from nilearn.plotting import plot_roi
    from scipy import ndimage as nd

    if out_file is None:
        fname, ext = op.splitext(op.basename(in_file))
        if ext == '.gz':
            fname, ext2 = op.splitext(fname)
            ext = ext2 + ext
        out_file = op.abspath('{}_spmask{}'.format(fname, ext))
        out_plot = op.abspath('{}_spmask.pdf'.format(fname))

    in_4d_nii = nb.load(in_file)
    orientation = nb.aff2axcodes(in_4d_nii.affine)

    if in_mask:
        mask_data = nb.load(in_mask).get_data()
        a = np.where(mask_data != 0)
        bbox = np.max(a[0]) - np.min(a[0]), np.max(a[1]) - \
            np.min(a[1]), np.max(a[2]) - np.min(a[2])
        longest_axis = np.argmax(bbox)

        # Input here is a binarized and intersected mask data from previous section
        dil_mask = nd.binary_dilation(
            mask_data, iterations=int(mask_data.shape[longest_axis] / 9))

        rep = list(mask_data.shape)
        rep[longest_axis] = -1
        new_mask_2d = dil_mask.max(axis=longest_axis).reshape(rep)

        rep = [1, 1, 1]
        rep[longest_axis] = mask_data.shape[longest_axis]
        new_mask_3d = np.logical_not(np.tile(new_mask_2d, rep))
    else:
        new_mask_3d = np.zeros(in_4d_nii.shape[:3]) == 1

    if orientation[0] in ['L', 'R']:
        new_mask_3d[0:2, :, :] = True
        new_mask_3d[-3:-1, :, :] = True
    else:
        new_mask_3d[:, 0:2, :] = True
        new_mask_3d[:, -3:-1, :] = True

    mask_nii = nb.Nifti1Image(new_mask_3d.astype(np.uint8), in_4d_nii.get_affine(),
                              in_4d_nii.get_header())
    mask_nii.to_filename(out_file)

    plot_roi(mask_nii, mean_img(in_4d_nii), output_file=out_plot)
    return out_file, out_plot

 def get_static_svg(self):
     '''Generate static svg of atlas (cannot manipulate in d3)'''
     svg_data = []
     with make_tmp_folder() as temp_dir:
         output_file='%s/atlas.svg' %(temp_dir)
         plotting.plot_roi(self.mr,annotate=False,draw_cross=False,cmap="nipy_spectral",black_bg=False, output_file=output_file)
         svg_file = open(output_file,'r')
         svg_data = svg_file.readlines()
         svg_file.close()
     return svg_data[4:]

def compute_all_subjects_mask():
    """ Computes the mask of all the subjects and the sesssions
    """
    masker = MultiNiftiMasker(mask_strategy='epi', memory=CACHE_DIR,
                              memory_level=2, n_jobs=10, verbose=5)
               
    imgs = dataset.func1 + dataset.func2
    masker.fit(imgs)
    masker.mask_img_.to_filename('all_subjects.nii.gz')
    plot_roi(masker.mask_img_)

def make_roi_image(nifti_file,png_img_file=None):
    """Make roi (mask) image"""
    nifti_file = str(nifti_file)
    mask_brain = plot_roi(nifti_file)
    if png_img_file:    
        mask_brain.savefig(png_img_file)
    plt.close('all')
    return mask_brain

def stripped_brain_overlay(in_file, overlay_file, out_file):
    import os.path
    import nibabel as nb
    import matplotlib as mpl
    mpl.use('Agg')
    from nilearn.plotting import plot_roi
    vmax = nb.load(in_file).get_data().reshape(-1).max()
    mask_display = plot_roi(
        in_file, overlay_file, output_file=out_file, title=out_file,
        display_mode="ortho", dim=-1, alpha=.3, vmax=vmax + 1)
    #  mask_display.bg_img(overlay_file)
    #  mask_display.title(out_file, x=0.01, y=0.99, size=15, color=None,
    #                     bgcolor=None, alpha=1)
    #  mask_display.display_mode = "yx"
    mask_display
    return os.path.abspath(out_file)

### Build a mask ##############################################################
# Thresholding
log_p_values[log_p_values < 5] = 0
plot_stat_map(new_img_like(fmri_img, log_p_values),
              mean_img, title='Thresholded p-values', annotate=False,
              colorbar=False, cut_coords=cut_coords)

# Binarization and intersection with VT mask
# (intersection corresponds to an "AND conjunction")
bin_p_values = (log_p_values != 0)
mask_vt_filename = haxby_dataset.mask_vt[0]
vt = nibabel.load(mask_vt_filename).get_data().astype(bool)
bin_p_values_and_vt = np.logical_and(bin_p_values, vt)

plot_roi(new_img_like(fmri_img, bin_p_values_and_vt.astype(np.int)),
         mean_img, title='Intersection with ventral temporal mask',
         cut_coords=cut_coords)

# Dilation
from scipy import ndimage
dil_bin_p_values_and_vt = ndimage.binary_dilation(bin_p_values_and_vt)
plot_roi(new_img_like(fmri_img, dil_bin_p_values_and_vt.astype(np.int)),
         mean_img, title='Dilated mask', cut_coords=cut_coords,
         annotate=False)

# Identification of connected components
plt.figure()
labels, n_labels = ndimage.label(dil_bin_p_values_and_vt)
first_roi_data = (labels == 1).astype(np.int)
second_roi_data = (labels == 2).astype(np.int)
fig_id = plt.subplot(2, 1, 1)

    get_tmaps=True)
localizer_anat_filename = localizer_dataset.anats[1]
localizer_cmap_filename = localizer_dataset.cmaps[1]
localizer_tmap_filename = localizer_dataset.tmaps[1]

###############################################################################
# demo the different plotting functions

# Plotting statistical maps
plotting.plot_stat_map(localizer_cmap_filename, bg_img=localizer_anat_filename,
                       threshold=3, title="plot_stat_map",
                       cut_coords=(36, -27, 66))

# Plotting glass brain
plotting.plot_glass_brain(localizer_tmap_filename, title='plot_glass_brain',
                          threshold=3)

# Plotting anatomical maps
plotting.plot_anat(haxby_anat_filename, title="plot_anat")

# Plotting ROIs (here the mask)
plotting.plot_roi(haxby_mask_filename, bg_img=haxby_anat_filename,
                  title="plot_roi")

# Plotting EPI haxby
mean_haxby_img = image.mean_img(haxby_func_filename)
plotting.plot_epi(mean_haxby_img, title="plot_epi")

import matplotlib.pyplot as plt
plt.show()

# To visualize results, we need to transform the clustering's labels back
# to a neuroimaging volume. For this, we use the NiftiMasker's
# inverse_transform method.
from nilearn.plotting import plot_roi, plot_epi, show

# Unmask the labels

# Avoid 0 label
labels = ward.labels_ + 1
labels_img = nifti_masker.inverse_transform(labels)

from nilearn.image import mean_img
mean_func_img = mean_img(func_filename)


first_plot = plot_roi(labels_img, mean_func_img, title="Ward parcellation",
                      display_mode='xz')

# common cut coordinates for all plots
cut_coords = first_plot.cut_coords

##################################################################
# labels_img is a Nifti1Image object, it can be saved to file with the
# following code:
labels_img.to_filename('parcellation.nii')


##################################################################
# Second, we illustrate the effect that the clustering has on the
# signal. We show the original data, and the approximation provided by
# the clustering by averaging the signal on each parcel.
#

for train, test in cv:
    svc.fit(fmri_masked[train], target[train])
    prediction = svc.predict(fmri_masked[test])
    cv_scores.append(np.sum(prediction == target[test])
                     / float(np.size(target[test])))

print(cv_scores)

### Unmasking #################################################################

# Retrieve the SVC discriminating weights
coef_ = svc.coef_

# Reverse masking thanks to the Nifti Masker
coef_img = nifti_masker.inverse_transform(coef_)

# Save the coefficients as a Nifti image
coef_img.to_filename('haxby_svc_weights.nii')

### Visualization #############################################################
import matplotlib.pyplot as plt
from nilearn.image.image import mean_img
from nilearn.plotting import plot_roi, plot_stat_map

mean_epi = mean_img(func_filename)
plot_stat_map(coef_img, mean_epi, title="SVM weights", display_mode="yx")

plot_roi(nifti_masker.mask_img_, mean_epi, title="Mask", display_mode="yx")

plt.show()

                           target_shape=shape, interpolation='nearest')
roi_masker = input_data.NiftiLabelsMasker(labels_img=atlas,
                                          mask_img=mask_filename)
roi_masker.fit(mask_filename)  # just to have it fitted

cut_coords = (30, -45, -12)

roi_score_img = roi_masker.inverse_transform(rsa[np.newaxis])
plotting.plot_stat_map(roi_score_img, title='RSA', cut_coords=cut_coords)
plt.savefig('rsa.png')
roi_score_img = roi_masker.inverse_transform(lm[np.newaxis])
plotting.plot_stat_map(roi_score_img, title='Linear model',
                       cut_coords=cut_coords)
plt.savefig('lm.png')

plotting.plot_roi(atlas, title="Harvard Oxford atlas", cut_coords=cut_coords)
# print labels

from scipy.stats import fligner
X = roi_masker.transform(func_filename)
y, session = np.loadtxt(haxby_dataset.session_target).astype('int').T
conditions = np.recfromtxt(haxby_dataset.conditions_target)['f0']
non_rest = conditions != b'rest'
conditions = conditions[non_rest]
y, session = y[non_rest], session[non_rest]
y = y[session < 4]

var_stat = np.zeros(X.shape[1])
for j, x in enumerate(X.T):
    _, var_stat[j] = fligner(
        x[y == 8], x[y == 1], x[y == 2], x[y == 3],

###########################################################################
# Visualize: Brain parcellations (Ward)
# -------------------------------------
#
# First, we display the parcellations of the brain image stored in attribute
# `labels_img_`
ward_labels_img = ward.labels_img_

# Now, ward_labels_img are Nifti1Image object, it can be saved to file
# with the following code:
ward_labels_img.to_filename('ward_parcellation.nii.gz')

from nilearn import plotting
from nilearn.image import mean_img, index_img

first_plot = plotting.plot_roi(ward_labels_img, title="Ward parcellation",
                               display_mode='xz')

# Grab cut coordinates from this plot to use as a common for all plots
cut_coords = first_plot.cut_coords
###########################################################################
# Compressed representation of Ward clustering
# --------------------------------------------
#
# Second, we illustrate the effect that the clustering has on the signal.
# We show the original data, and the approximation provided by the
# clustering by averaging the signal on each parcel.

# Grab number of voxels from attribute mask image (mask_img_).
import numpy as np
original_voxels = np.sum(ward.mask_img_.get_data())

###########################################################################
# Convert the fMRI volume's to a data matrix
# ..........................................
#
# We will use the :class:`nilearn.input_data.NiftiMasker` to extract the
# fMRI data on a mask and convert it to data series.
#
# The mask is a mask of the Ventral Temporal streaming coming from the
# Haxby study:
mask_filename = haxby_dataset.mask_vt[0]

# Let's visualize it, using the subject's anatomical image as a
# background
from nilearn import plotting
plotting.plot_roi(mask_filename, bg_img=haxby_dataset.anat[0],
                 cmap='Paired')

###########################################################################
# Now we use the NiftiMasker.
#
# We first create a masker, giving it the options that we care
# about. Here we use standardizing of the data, as it is often important
# for decoding
from nilearn.input_data import NiftiMasker
masker = NiftiMasker(mask_img=mask_filename, standardize=True)

# We give the masker a filename and retrieve a 2D array ready
# for machine learning with scikit-learn
fmri_masked = masker.fit_transform(fmri_filename)

###########################################################################

              title="p-values", cut_coords=(coronal, sagittal, axial))

### Build a mask ##############################################################

# Thresholding
pvalues[pvalues < 5] = 0
plot_stat_map(nibabel.Nifti1Image(pvalues, fmri_img.get_affine()), mean_img,
              title='Thresholded p-values',
              cut_coords=(coronal, sagittal, axial))

# Binarization and intersection with VT mask
bin_pvalues = (pvalues != 0)
vt = nibabel.load(haxby_files.mask_vt[0]).get_data().astype(bool)
bin_pvalues_and_vt = np.logical_and(bin_pvalues, vt)
plot_roi(nibabel.Nifti1Image(bin_pvalues_and_vt.astype(np.int), 
                             fmri_img.get_affine()), 
         mean_img, title='Intersection with ventral temporal mask',
         cut_coords=(coronal, sagittal, axial))

# Dilation
from scipy import ndimage
dil_bin_pvalues_and_vt = ndimage.binary_dilation(bin_pvalues_and_vt)
plot_roi(nibabel.Nifti1Image(dil_bin_pvalues_and_vt.astype(np.int), 
                             fmri_img.get_affine()), 
         mean_img, title='Dilated mask', cut_coords=(coronal, sagittal, axial))

# Identification of connected components
labels, n_labels = ndimage.label(dil_bin_pvalues_and_vt)
plot_roi(nibabel.Nifti1Image(labels, fmri_img.get_affine()), 
         mean_img, title='Connected components',
         cut_coords=(coronal, sagittal, axial))
plt.show()

"""

##############################################################################
# The original Yeo atlas
# -----------------------

# First we fetch the Yeo atlas
from nilearn import datasets

atlas_yeo_2011 = datasets.fetch_atlas_yeo_2011()
atlas_yeo = atlas_yeo_2011.thick_7

# Let's now plot it
from nilearn import plotting

plotting.plot_roi(atlas_yeo, title='Original Yeo atlas',
                  cut_coords=(8, -4, 9), colorbar=True, cmap='Paired')

##############################################################################
# The original Yeo atlas has 7 labels, that is indicated in the colorbar.
# The colorbar also shows the correspondence between the color and the label
#
# Note that these 7 labels correspond actually to networks that comprise
# several regions. We are going to split them up.

##############################################################################
# Relabeling the atlas into separated regions
# ---------------------------------------------
#
# Now we use the connected_label_regions to break appart the networks
# of the Yeo atlas into separated regions
from nilearn.regions import connected_label_regions

# -*- coding: utf-8 -*-
"""
Created on Wed Apr  1 09:16:38 2015

@author: mr243268
"""
import os
from loader import load_dynacomp, set_figure_base_dir
from nilearn.plotting import plot_roi

dataset = load_dynacomp()

FIG_DIR = set_figure_base_dir('rois')

for i in range(len(dataset.subjects)):
    for k in sorted(dataset.rois[i].keys()):
        output_file = os.path.join(FIG_DIR, k)
        plot_roi(dataset.rois[i][k], title=k, output_file=output_file)
    break

extraction.fit()
regions_img = extraction.regions_img_

################################################################################
# Visualization
# Show region extraction results by importing image & plotting utilities
from nilearn import plotting
from nilearn.image import index_img
from nilearn.plotting import find_xyz_cut_coords

# Showing region extraction results using 4D maps visualization tool
plotting.plot_prob_atlas(regions_img, display_mode='z', cut_coords=2,
                         view_type='contours', title="Regions extracted.")

# To reduce the complexity, we choose to display all the regions
# extracted from network 3
import numpy as np

DMN_network = index_img(atlas_networks, 3)
plotting.plot_roi(DMN_network, display_mode='z', cut_coords=1,
                  title='Network 3')

regions_indices_network3 = np.where(np.array(extraction.index_) == 3)
for index in regions_indices_network3[0]:
    cur_img = index_img(extraction.regions_img_, index)
    coords = find_xyz_cut_coords(cur_img)
    plotting.plot_roi(cur_img, display_mode='z', cut_coords=coords[2:3],
                      title="Blob of network3")
            
plotting.show()

def plot_mask(pet_files, pet_imgs):
    for pi, pf in zip(pet_imgs, pet_files):
        mask_path = os.path.join('figures', 'mask',
                                 pf.split('/')[-1].split('.')[0]) 
        plot_roi(masker.mask_img_, pi, output_file=mask_path,
                 title=pf.split('/')[-1].split('.')[0])

pet_files = []
pet_img = []
for idx, row in data.iterrows():
    pet_file = glob.glob(os.path.join(BASE_DIR,
                                      'I' + str(row.Image_ID_y), 'wI*.nii'))
    if len(pet_file)>0:
        pet_files.append(pet_file[0])
        img = nib.load(pet_file[0])
        pet_img.append(img)

masker = NiftiMasker(mask_strategy='epi',
                     mask_args=dict(opening=1))
masker.fit(pet_files)

plot_roi(masker.mask_img_, pet_file[0])

pet_masked = masker.transform_niimgs(pet_files, n_jobs=6)
pet_masked = np.vstack(pet_masked)
nb_vox = pet_masked.shape[1]

groups = [['AD', 'Normal'], ['AD', 'MCI'], ['MCI', 'LMCI'], ['MCI', 'Normal']]

for gr in groups:
    gr1_idx = data[data.DX_Group == gr[0]].index.values
    gr2_idx = data[data.DX_Group == gr[1]].index.values
    
    gr1_f = pet_masked[gr1_idx, :]
    gr2_f = pet_masked[gr2_idx, :]
    
    t_masked, p_masked = stats.ttest_ind(gr1_f, gr2_f)

                                               get_anats=True)


###############################################################################
# demo the different plotting function

# Plotting statistical maps
plotting.plot_stat_map(localizer.cmaps[3], bg_img=localizer.anats[3],
                       threshold=3, title="plot_stat_map",
                       cut_coords=(36, -27, 66))

# Plotting anatomical maps
plotting.plot_anat(haxby.anat[0], title="plot_anat")

# Plotting ROIs (here the mask)
plotting.plot_roi(haxby.mask_vt[0], bg_img=haxby.anat[0], title="plot_roi")

# Plotting EPI haxby
mean_haxby_img = image.mean_img(haxby.func[0])
plotting.plot_epi(mean_haxby_img, title="plot_epi")


###############################################################################
# demo the different display_mode

plotting.plot_stat_map(localizer.cmaps[3], display_mode='ortho',
                       cut_coords=(36, -27, 60),
                       title="display_mode='ortho', cut_coords=(36, -27, 60)")

plotting.plot_stat_map(localizer.cmaps[3], display_mode='z', cut_coords=5,
                       title="display_mode='z', cut_coords=5")

l, n = find_region_names_using_cut_coords(dmn_coords, atlas_img,
                                          labels=labels)

# where 'l' indicates new labels generated according to given atlas labels and
# coordinates
new_labels = l

# where 'n' indicates brain regions names labelled, generated according to given
# labels
region_names_involved = n

######################################################################
# Let's visualize
from nilearn.image import load_img
from nilearn import plotting
from nilearn.image import new_img_like

atlas_img = load_img(atlas_img)
affine = atlas_img.get_affine()
atlas_data = atlas_img.get_data()

for i, this_label in enumerate(new_labels):
    this_data = (atlas_data == this_label)
    this_data = this_data.astype(int)
    this_img = new_img_like(atlas_img, this_data, affine)
    plotting.plot_roi(this_img, cut_coords=dmn_coords[i],
                      title=region_names_involved[i])

plotting.show()