.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/plot_3_structural_properties.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_plot_3_structural_properties.py: Structural poperties of the MEG-informed parcellations ====================================================== .. GENERATED FROM PYTHON SOURCE LINES 7-8 Import the required packages .. GENERATED FROM PYTHON SOURCE LINES 8-24 .. code-block:: Python import os.path as op import numpy as np import pickle import matplotlib.pyplot as plt import matplotlib.patches as patches from mne import (read_forward_solution, pick_types_forward, convert_forward_solution, read_labels_from_annot, spatial_src_adjacency, read_source_spaces) from mne.viz import get_brain_class from mne.datasets import sample import megicparc .. GENERATED FROM PYTHON SOURCE LINES 25-26 Define input parameters for the flame algorithm running in megicperc .. GENERATED FROM PYTHON SOURCE LINES 26-37 .. code-block:: Python gamma_tot = np.arange(0, 1.01, 0.2) knn_tot = [10, 20, 30, 40] theta = 0.05 parc = 'aparc' sensors_meg = 'grad' folder_fl = op.join('..', 'data', 'data_mne_sample') string_target_file = op.join(folder_fl, '{:s}_flame_grad_k{:d}_gamma{:1.2f}_theta{:1.2f}.pkl') .. GENERATED FROM PYTHON SOURCE LINES 38-39 Load lead-field matrix, source-space and anatomy-based Atlas .. GENERATED FROM PYTHON SOURCE LINES 39-61 .. code-block:: Python data_path = sample.data_path() subjects_dir = op.join(data_path, 'subjects') subject = 'sample' fwd_file = op.join(data_path, 'MEG', subject, 'sample_audvis-meg-eeg-oct-6-fwd.fif') src_file = op.join(folder_fl, 'source_space_distance-src.fif') fwd = read_forward_solution(fwd_file) fwd = pick_types_forward(fwd, meg=sensors_meg, eeg=False, ref_meg=False, exclude='bads') fwd = convert_forward_solution(fwd, surf_ori=True, force_fixed=True, use_cps=True) src = read_source_spaces(src_file) fwd['src'] = src n_vert = fwd['src'][0]['nuse'] + fwd['src'][1]['nuse'] label_lh = read_labels_from_annot(subject=subject, parc=parc, hemi='lh', subjects_dir=subjects_dir) label_rh = read_labels_from_annot(subject=subject, parc=parc, hemi='rh', subjects_dir=subjects_dir) label = label_lh + label_rh .. rst-class:: sphx-glr-script-out .. code-block:: none Reading forward solution from /u/29/sommars1/unix/mne_data/MNE-sample-data/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif... Reading a source space... Computing patch statistics... Patch information added... Distance information added... [done] Reading a source space... Computing patch statistics... Patch information added... Distance information added... [done] 2 source spaces read Desired named matrix (kind = 3523) not available Read MEG forward solution (7498 sources, 306 channels, free orientations) Desired named matrix (kind = 3523) not available Read EEG forward solution (7498 sources, 60 channels, free orientations) Forward solutions combined: MEG, EEG Source spaces transformed to the forward solution coordinate frame 203 out of 366 channels remain after picking Average patch normals will be employed in the rotation to the local surface coordinates.... Converting to surface-based source orientations... [done] Reading a source space... Computing patch statistics... Patch information added... Distance information added... [done] Reading a source space... Computing patch statistics... Patch information added... Distance information added... [done] 2 source spaces read Reading labels from parcellation... read 34 labels from /u/29/sommars1/unix/mne_data/MNE-sample-data/subjects/sample/label/lh.aparc.annot Reading labels from parcellation... read 34 labels from /u/29/sommars1/unix/mne_data/MNE-sample-data/subjects/sample/label/rh.aparc.annot .. GENERATED FROM PYTHON SOURCE LINES 62-63 Initialization .. GENERATED FROM PYTHON SOURCE LINES 63-69 .. code-block:: Python num_roi = np.zeros((np.size(knn_tot), np.size(gamma_tot))) n_cc = {x : [] for x in \ ['k%d_gamma%1.2f'%(knn, gamma) for knn in knn_tot for gamma in gamma_tot]} coverage = np.zeros((np.size(knn_tot), np.size(gamma_tot))) .. GENERATED FROM PYTHON SOURCE LINES 70-71 Analysis with anatomy-based parcels. .. GENERATED FROM PYTHON SOURCE LINES 71-83 .. code-block:: Python proj_anrois = megicparc.labels_to_array(label, fwd['src']) # Number of regions if proj_anrois['outliers'] > 0: aux_num_roi = len(proj_anrois['parcel']) - 1 else: aux_num_roi = len(proj_anrois['parcel']) num_roi_an = aux_num_roi # Coverage coverage_an = 100 * (n_vert - proj_anrois['outliers']) / n_vert .. GENERATED FROM PYTHON SOURCE LINES 84-85 Analysis with meg-informed parcels .. GENERATED FROM PYTHON SOURCE LINES 85-108 .. code-block:: Python bmesh_adj_mat = np.array(spatial_src_adjacency(fwd['src']).todense()) for idx_g, gamma in enumerate(gamma_tot): for idx_k, knn in enumerate(knn_tot): target_file = string_target_file.format( subject, knn, gamma, theta) print('Loading %s'%target_file) with open(target_file, 'rb') as aux_lf: flame_data = pickle.load(aux_lf) # Number of regions num_roi[idx_k, idx_g] = flame_data['centroids'] # Coverage coverage[idx_k, idx_g] = \ 100 * (n_vert - flame_data['outliers']) / n_vert # Number of connected components n_cc['k%d_gamma%1.2f'%(knn, gamma)] += \ [len(megicparc.compute_connected_components( bmesh_adj_mat, flame_data['parcel'][ir])) \ for ir in range(flame_data['centroids'])] .. rst-class:: sphx-glr-script-out .. code-block:: none -- number of adjacent vertices : 8196 /m/home/home2/29/sommars1/data/Documents/megicparc/examples/plot_3_structural_properties.py:86: RuntimeWarning: 8.5% of original source space vertices have been omitted, tri-based adjacency will have holes. Consider using distance-based adjacency or morphing data to all source space vertices. bmesh_adj_mat = np.array(spatial_src_adjacency(fwd['src']).todense()) Loading ../data/data_mne_sample/sample_flame_grad_k10_gamma0.00_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k20_gamma0.00_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k30_gamma0.00_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k40_gamma0.00_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k10_gamma0.20_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k20_gamma0.20_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k30_gamma0.20_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k40_gamma0.20_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k10_gamma0.40_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k20_gamma0.40_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k30_gamma0.40_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k40_gamma0.40_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k10_gamma0.60_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k20_gamma0.60_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k30_gamma0.60_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k40_gamma0.60_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k10_gamma0.80_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k20_gamma0.80_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k30_gamma0.80_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k40_gamma0.80_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k10_gamma1.00_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k20_gamma1.00_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k30_gamma1.00_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k40_gamma1.00_theta0.05.pkl .. GENERATED FROM PYTHON SOURCE LINES 109-111 Compute mean and SEM of the number of connected components over the regions of a MEG-informed parcellations .. GENERATED FROM PYTHON SOURCE LINES 111-121 .. code-block:: Python n_cc_ave = np.zeros((np.size(knn_tot), np.size(gamma_tot))) n_cc_sem = np.zeros((np.size(knn_tot), np.size(gamma_tot))) for idx_g, gamma in enumerate(gamma_tot): for idx_k, knn in enumerate(knn_tot): aux_cc = np.asanyarray(n_cc['k%d_gamma%1.2f'%(knn, gamma)]) n_cc_ave[idx_k, idx_g] = np.mean(aux_cc) n_cc_sem[idx_k, idx_g] = np.std(aux_cc, ddof=1) / np.sqrt(np.size(aux_cc)) .. GENERATED FROM PYTHON SOURCE LINES 122-123 Additional analysis: impact of the anatomical constraints .. GENERATED FROM PYTHON SOURCE LINES 123-161 .. code-block:: Python k_nn_ac = [30] theta_ac = [0, 0.05] gamma_ac = np.arange(1, 0.4, -0.2) tolls = [0.05] spread = {} num_roi_constr = {} for idx_t, theta in enumerate(theta_ac): for idx_k, knn in enumerate(k_nn_ac): for idx_g, gamma in enumerate(gamma_ac): # 6.1. Reload parcellation target_file = string_target_file.format(subject, knn, gamma, theta) print('Loading %s'%target_file) with open(target_file, 'rb') as aux_lf: flame_data = pickle.load(aux_lf) num_roi_constr['k%d_gamma%1.2f_theta%1.2f'%(knn, gamma, theta)] = \ flame_data['centroids'] # 6.2. Compute cotingency matrix conting_mat = np.array([ np.array([ np.intersect1d(p_fl, p_an).shape[0] for p_an in proj_anrois['parcel']]) for p_fl in flame_data['parcel'][0:flame_data['centroids']]]) conting_mat = conting_mat / \ np.sum(conting_mat, axis=1)[:, np.newaxis] # 6.3. Compute spread of flame regions over anatomical spread['k%d_gamma%1.2f_theta%1.2f'%(knn, gamma, theta)] = \ np.zeros((flame_data['centroids'], np.size(tolls))) for idx_t, toll in enumerate(tolls): spread['k%d_gamma%1.2f_theta%1.2f'%(knn, gamma, theta)] [:, idx_t] = \ np.sum(conting_mat > toll, axis=1) .. rst-class:: sphx-glr-script-out .. code-block:: none Loading ../data/data_mne_sample/sample_flame_grad_k30_gamma1.00_theta0.00.pkl Loading ../data/data_mne_sample/sample_flame_grad_k30_gamma0.80_theta0.00.pkl Loading ../data/data_mne_sample/sample_flame_grad_k30_gamma0.60_theta0.00.pkl Loading ../data/data_mne_sample/sample_flame_grad_k30_gamma1.00_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k30_gamma0.80_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k30_gamma0.60_theta0.05.pkl .. GENERATED FROM PYTHON SOURCE LINES 162-163 General parameters for plotting .. GENERATED FROM PYTHON SOURCE LINES 163-177 .. code-block:: Python plt.rc('text', usetex=True) plt.rc('font', family='serif', size=22) plt.rcParams['axes.grid'] = True plt.rcParams['lines.linewidth'] = 3 plt.rcParams['errorbar.capsize'] = 3.5 plt.rcParams["figure.figsize"] = [12, 8] plt.rcParams['patch.force_edgecolor'] = True colors = np.array([np.array([0, 1, 0]), np.array([0, 0, 1]), np.array([1, 0, 0]), np.array([0, 0.9655, 1])]) .. GENERATED FROM PYTHON SOURCE LINES 178-179 Plot 1. Number of regions .. GENERATED FROM PYTHON SOURCE LINES 179-195 .. code-block:: Python f_numreg = plt.figure() for idx_k, knn in enumerate(knn_tot): plt.plot(gamma_tot, num_roi[idx_k,], color=colors[idx_k], label='$k$ = %d'%knn) plt.plot(gamma_tot, num_roi_an*np.ones(np.size(gamma_tot)), 'k--', label='DK') plt.legend(loc=[0.68, 0.47], fontsize=18) plt.xlabel(r'Weight of the spatial distances', fontsize=25) plt.ylabel(r'Number of parcels', fontsize=25) plt.ylim(0, 350) plt.xlim(-0.05, 1.05) f_numreg.set_size_inches(8.5, 6) plt.tight_layout(pad=1.5) .. image-sg:: /auto_examples/images/sphx_glr_plot_3_structural_properties_001.png :alt: plot 3 structural properties :srcset: /auto_examples/images/sphx_glr_plot_3_structural_properties_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 196-197 Plot 2. Number of connected components .. GENERATED FROM PYTHON SOURCE LINES 197-212 .. code-block:: Python f_cc = plt.figure() for idx_k, knn in enumerate(knn_tot): plt.errorbar(gamma_tot, n_cc_ave[idx_k, :], yerr=n_cc_sem[idx_k, :], label='$k$ = %d'%knn, color=colors[idx_k]) plt.legend(fontsize=18) plt.xlabel(r'Weight of the spatial distances', fontsize=25) plt.ylabel(r'Number connected components', fontsize=25) plt.ylim(0, 30) plt.xlim(-0.05, 1.05) #plt.show() f_cc.set_size_inches(8.5, 6) plt.tight_layout(pad=1.5) .. image-sg:: /auto_examples/images/sphx_glr_plot_3_structural_properties_002.png :alt: plot 3 structural properties :srcset: /auto_examples/images/sphx_glr_plot_3_structural_properties_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 213-214 . Plot 3. (Additional) Impact of anatomical constraints .. GENERATED FROM PYTHON SOURCE LINES 214-289 .. code-block:: Python bins = np.arange(-0.5, 69, 1) alpha = 0.7 knn = k_nn_ac[0] toll = tolls[0] fac, axac = plt.subplots(1, np.size(gamma_ac), figsize=[13, 5]) #fac.suptitle('$k =$ %d'%knn) idx_t = 0 for idx_g, gamma in enumerate(gamma_ac): aux_th1 = spread['k%d_gamma%1.2f_theta%1.2f'%(knn, gamma, theta_ac[1])][:, idx_t] axac[idx_g].hist(aux_th1, bins=bins, alpha=alpha, weights=np.ones_like(aux_th1)/np.shape(aux_th1)[0]*100, label=r'$\theta = $ %1.2f'%theta_ac[1]) aux_th2 = spread['k%d_gamma%1.2f_theta%1.2f'%(knn, gamma, theta_ac[0])] [:, idx_t] axac[idx_g].hist(aux_th2, bins=bins, alpha=alpha, weights=np.ones_like(aux_th2)/np.shape(aux_th2)[0]*100, label=r'$\theta = $ %1.2f'%theta_ac[0]) axac[idx_g].set_ylim(0, 100) axac[idx_g].set_xlim(-0.5, 10.5) #axac[idx_t, idx_g].legend() if idx_t == 0: axac[idx_g].set_title(r'$\gamma = %1.1f$ '%(gamma)) axac[idx_g].set_xticks(np.arange(1, 11, 2)) axac[idx_g].text(6.5, 85, '%d ;'%(num_roi_constr['k%d_gamma%1.2f_theta0.05'%(knn, gamma)]), ha='center', color='dodgerblue', fontweight='bold') axac[idx_g].text(8.5, 85, '%d'%(num_roi_constr['k%d_gamma%1.2f_theta0.00'%(knn, gamma)]), ha='center', color='darkorange', fontweight='bold') axac[idx_g].text(7.2, 75, 'parcels', ha='center') rect = patches.FancyBboxPatch((5, 72), 4.5, 22, boxstyle='round', facecolor='white', edgecolor='lightgray', zorder=2) axac[idx_g].add_patch(rect) else: axac[idx_g].set_xticks(np.arange(1, 11, 2)) if idx_g == 0: axac[idx_g].set_ylabel( r'Number of parcels ($\%$)', fontsize=30) else: axac[idx_g].set_yticklabels([]) axac[-1].legend(bbox_to_anchor=(1, 0.5), loc='center left', borderaxespad=0.2) fac.text(0.45, 0.02, 'Number of intersecting atlas parcels', ha='center', fontsize=30) plt.tight_layout(pad=1.5) gamma_tplot = np.arange(0, 1, 0.1) theta_tplot = np.array([0.01, 0.025, 0.05, 0.1, 0.2]) f_tg = np.array([theta * (gamma_tplot / (1-gamma_tplot)) \ for theta in theta_tplot]) colors_tplot = np.array([np.array([0, 1, 0]), np.array([0, 0, 1]), np.array([1, 0, 0]), np.array([0, 0.9655, 1]), np.array([1, 0.4828, 0.8621])]) f_tplot = plt.figure() for idx_t, theta in enumerate(theta_tplot): plt.plot(gamma_tplot, f_tg[idx_t], label=r'$\theta$ = %1.3f'%theta, color=colors_tplot[idx_t]) plt.legend() plt.xlim([0, 0.9]) plt.ylim([0, 2]) plt.xlabel(r'$\gamma$', fontsize=30) plt.ylabel(r'$\theta$ $\frac{\gamma}{1-\gamma}$', fontsize=30) f_tplot.set_size_inches(8, 6.5) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/images/sphx_glr_plot_3_structural_properties_003.png :alt: $\gamma = 1.0$ , $\gamma = 0.8$ , $\gamma = 0.6$ :srcset: /auto_examples/images/sphx_glr_plot_3_structural_properties_003.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_plot_3_structural_properties_004.png :alt: plot 3 structural properties :srcset: /auto_examples/images/sphx_glr_plot_3_structural_properties_004.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 290-291 Plot 5. Spatial cohesion .. GENERATED FROM PYTHON SOURCE LINES 291-326 .. code-block:: Python Brain = get_brain_class() knn_plot = 30 gamma_plot = [0, 0.4, 0.6, 1] theta_plot = 0.05 cart_coord = [-0.03, 0.025, 0.08] src = fwd['src'] V = np.concatenate((src[0]['rr'][src[0]['vertno']], src[1]['rr'][src[1]['vertno']]), axis=0) nvert = src[0]['nuse'] + src[1]['nuse'] idx_vert = np.argmin(np.linalg.norm(V - cart_coord, axis=1)) for gamma in gamma_plot: target_file = string_target_file.format( subject, knn_plot, gamma, theta_plot) print('Loading %s'%target_file) with open(target_file, 'rb') as aux_lf: flame_data = pickle.load(aux_lf) flame_labels = megicparc.store_flame_labels(flame_data, src, subject) parcels_vector = np.zeros(nvert, dtype='int') - 1 for ir in range(flame_data['centroids']): parcels_vector[flame_data['parcel'][ir]] = ir sel_roi = parcels_vector[idx_vert] index = [sel_roi + 1] brain = megicparc.plot_flame_labels(index, flame_labels, src, subject, subjects_dir, surf='inflated', color = [1, 0.64, 0.], plot_region=True, plot_points=False, plot_borders=False) brain.show_view(azimuth=173, elevation= 60) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/images/sphx_glr_plot_3_structural_properties_005.png :alt: plot 3 structural properties :srcset: /auto_examples/images/sphx_glr_plot_3_structural_properties_005.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_plot_3_structural_properties_006.png :alt: plot 3 structural properties :srcset: /auto_examples/images/sphx_glr_plot_3_structural_properties_006.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_plot_3_structural_properties_007.png :alt: plot 3 structural properties :srcset: /auto_examples/images/sphx_glr_plot_3_structural_properties_007.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_plot_3_structural_properties_008.png :alt: plot 3 structural properties :srcset: /auto_examples/images/sphx_glr_plot_3_structural_properties_008.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none Loading ../data/data_mne_sample/sample_flame_grad_k30_gamma0.00_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k30_gamma0.40_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k30_gamma0.60_theta0.05.pkl Loading ../data/data_mne_sample/sample_flame_grad_k30_gamma1.00_theta0.05.pkl .. GENERATED FROM PYTHON SOURCE LINES 327-328 Plot 7. Reduced source space on top of DK atlas .. GENERATED FROM PYTHON SOURCE LINES 328-342 .. code-block:: Python brain_redV = Brain(subject, hemi='both', surf='inflated', background='white', subjects_dir=subjects_dir, alpha=1) # Superimpose anatomical regions for ir in range(len(label)): brain_redV.add_label(label[ir], hemi=label[ir].hemi, alpha=0.9) # Superimpose centroids megicparc.plot_flame_centroids(flame_data, src, subject, subjects_dir, brain_redV) brain_redV.show_view(azimuth=170, elevation=90) "" .. image-sg:: /auto_examples/images/sphx_glr_plot_3_structural_properties_009.png :alt: plot 3 structural properties :srcset: /auto_examples/images/sphx_glr_plot_3_structural_properties_009.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none '' .. rst-class:: sphx-glr-timing **Total running time of the script:** (1 minutes 38.736 seconds) .. _sphx_glr_download_auto_examples_plot_3_structural_properties.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_3_structural_properties.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_3_structural_properties.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_