.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/example_spangy.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_example_spangy.py: .. _example_spangy: =================================== example of SPANGY (spectral decomposition) tools in slam =================================== .. GENERATED FROM PYTHON SOURCE LINES 8-14 .. code-block:: Python # Authors: # License: MIT # sphinx_gallery_thumbnail_number = 2 .. GENERATED FROM PYTHON SOURCE LINES 15-19 NOTE: there is no visualization tool in slam, but we provide at the end of this script exemplare code to do the visualization with an external solution ############################################################################## .. GENERATED FROM PYTHON SOURCE LINES 21-22 importation of slam modules .. GENERATED FROM PYTHON SOURCE LINES 22-27 .. code-block:: Python import numpy as np import slam.io as sio import slam.curvature as scurv import slam.spangy as spgy .. GENERATED FROM PYTHON SOURCE LINES 28-29 LOAD MESH .. GENERATED FROM PYTHON SOURCE LINES 29-37 .. code-block:: Python mesh = sio.load_mesh( '../examples/data/example_mesh.gii') vertices = mesh.vertices num_vertices = len(vertices) print('{} vertices'.format(num_vertices)) N = 1500 # N should be < to the number of vertices. .. rst-class:: sphx-glr-script-out .. code-block:: none 2328 vertices .. GENERATED FROM PYTHON SOURCE LINES 38-39 Compute eigenpairs and mass matrix .. GENERATED FROM PYTHON SOURCE LINES 39-41 .. code-block:: Python eigVal, eigVects, lap_b = spgy.eigenpairs(mesh, N) .. rst-class:: sphx-glr-script-out .. code-block:: none Computing Laplacian Computing mesh weights of type fem -edge length threshold needed for 0 values = 0.0 % -number of Nan in weights: 0 = 0.0 % -number of Negative values in weights: 936 = 6.706792777300086 % -nb Nan in Laplacian : 0 -nb Inf in Laplacian : 0 .. GENERATED FROM PYTHON SOURCE LINES 42-43 CURVATURE .. GENERATED FROM PYTHON SOURCE LINES 43-47 .. code-block:: Python PrincipalCurvatures, PrincipalDir1, PrincipalDir2 = \ scurv.curvatures_and_derivatives(mesh) mean_curv = 0.5 * (PrincipalCurvatures[0, :] + PrincipalCurvatures[1, :]) .. rst-class:: sphx-glr-script-out .. code-block:: none Calculating vertex normals .... Please wait Finished calculating vertex normals Calculating curvature tensors ... Please wait Finished Calculating curvature tensors Calculating Principal Components ... Please wait Finished Calculating principal components .. GENERATED FROM PYTHON SOURCE LINES 48-49 WHOLE BRAIN MEAN-CURVATURE SPECTRUM .. GENERATED FROM PYTHON SOURCE LINES 49-79 .. code-block:: Python grouped_spectrum, group_indices, coefficients, nlevels \ = spgy.spectrum(mean_curv, lap_b, eigVects, eigVal) levels = len(group_indices) # a. Whole brain parameters mL_in_MM3 = 1000 CM2_in_MM2 = 100 volume = mesh.volume surface_area = mesh.area afp = np.sum(grouped_spectrum[1:]) print('** a. Whole brain parameters **') print('Volume = %d mL, Area = %d cm², Analyze Folding Power = %f,' % (np.floor(volume / mL_in_MM3), np.floor(surface_area / CM2_in_MM2), afp)) # b. Band number of parcels print('** b. Band number of parcels **') print('B4 = %f, B5 = %f, B6 = %f' % (0, 0, 0)) # c. Band power print('** c. Band power **') print('B4 = %f, B5 = %f, B6 = %f' % (grouped_spectrum[4], grouped_spectrum[5], grouped_spectrum[6])) # d. Band relative power print('** d. Band relative power **') print('B4 = %0.5f, B5 = %0.5f , B6 = %0.5f' % (grouped_spectrum[4] / afp, grouped_spectrum[5] / afp, grouped_spectrum[6] / afp)) .. rst-class:: sphx-glr-script-out .. code-block:: none ** a. Whole brain parameters ** Volume = 30 mL, Area = 64 cm², Analyze Folding Power = 37.227117, ** b. Band number of parcels ** B4 = 0.000000, B5 = 0.000000, B6 = 0.000000 ** c. Band power ** B4 = 13.947520, B5 = 5.143363, B6 = 0.897768 ** d. Band relative power ** B4 = 0.37466, B5 = 0.13816 , B6 = 0.02412 .. GENERATED FROM PYTHON SOURCE LINES 80-81 LOCAL SPECTRAL BANDS .. GENERATED FROM PYTHON SOURCE LINES 81-100 .. code-block:: Python loc_dom_band, frecomposed = spgy.local_dominance_map(coefficients, mean_curv, levels, group_indices, eigVects) # WHOLE BRAIN MEAN-CURVATURE<=0 & MEAN-CURVATURE>0 SPECTRI # -------------------------------------------------------- # Define negative mean curvature subsignal mean_curv_sulci = np.zeros((mean_curv.shape)) mean_curv_sulci[mean_curv <= 0] = mean_curv[mean_curv <= 0] grouped_spectrum_sulci, group_indices_sulci, coefficients_sulci, _ \ = spgy.spectrum(mean_curv_sulci, lap_b, eigVects, eigVal) # Define positive mean curvature subsignal mean_curv_gyri = np.zeros((mean_curv.shape)) mean_curv_gyri[mean_curv > 0] = mean_curv[mean_curv > 0] grouped_spectrum_gyri, group_indices_gyri, coefficients_gyri, _ \ = spgy.spectrum(mean_curv_gyri, lap_b, eigVects, eigVal) .. GENERATED FROM PYTHON SOURCE LINES 101-301 VISUALIZATION USING EXTERNAL TOOLS ############################################################################ import slam.plot as splt import matplotlib.pyplot as plt from matplotlib.colors import ListedColormap import pyvista as pv ############################################################################### # Plot of mean curvature on the mesh visb_sc = splt.visbrain_plot( mesh=mesh, tex=mean_curv, caption='Mean Curvature', cmap='jet') visb_sc.preview() # Plot coefficients and bands for all mean curvature signal fig, (ax1, ax2) = plt.subplots(1, 2) ax1.scatter(np.sqrt(eigVal/2*np.pi), coefficients, marker='+', s=10, linewidths=0.5) #ax1.plot(np.sqrt(eigVal[1:]) / (2 * np.pi), # coefficients[1:]) # remove B0 coefficients #ax1.scatter(np.sqrt(eigVal[1:]/2*np.pi), # coefficients[1:], marker='+', s=10, linewidths=0.5) # remove B0 coefficients ax1.set_xlabel('Frequency (m⁻¹)') ax1.set_ylabel('Coefficients') # print(grouped_spectrum) ax2.bar(np.arange(0, levels), grouped_spectrum) #ax2.bar(np.arange(1, nlevels), grouped_spectrum[1:]) # remove B0 ax2.set_xlabel('Spangy Frequency Bands') ax2.set_ylabel('Power Spectrum') plt.show() # Plot of spectral dominant bands on the mesh visb_sc = splt.visbrain_plot(mesh=mesh, tex=loc_dom_band, caption='Local Dominant Band', cmap='jet') visb_sc.preview() # Plot mean curvature coefficients & compacted power spectrum characterizing # either Sulci either Gyri folding pattern # --------------------------------------------------------------------------- coefficients_colors_sulci \ = plot_global_coefficients_and_bands_sulci_or_gyri( group_indices_sulci, coefficients_sulci, colormap_sulci ) coefficients_colors_gyri \ = plot_global_coefficients_and_bands_sulci_or_gyri( group_indices_gyri, coefficients_gyri, colormap_gyri ) fig, axs = plt.subplots(2, 2) # GLOBAL FOLDING PATTERN OF SULCI #axs[0,0].plot(np.sqrt(eigVal/2*np.pi), # coefficients_sulci, color=coefficients_colors_sulci) axs[0,0].scatter(np.sqrt(eigVal/2*np.pi), coefficients_sulci, marker='+', s=10, linewidths=0.5, color=coefficients_colors_sulci) #axs[0,0].scatter(np.sqrt(eigVal[1:]/2*np.pi), # coefficients_sulci[1:], marker='+', s=10, linewidths=0.5, # color=coefficients_colors_sulci[1:]) # remove B0 coefficient axs[0,0].set_xlabel('Frequency (m⁻¹)') axs[0,0].set_ylabel('Coefficients mean_curv<=0') axs[0,1].bar(np.arange(0, nlevels), grouped_spectrum_sulci.squeeze(), color=colormap_sulci) #axs[0,1].bar(np.arange(1, nlevels), # grouped_spectrum_sulci[1:].squeeze(), # color=colormap_sulci[1:]) # remove B0 axs[0,1].set_xlabel('Spangy frequency bands') axs[0,1].set_ylabel('Power spectrum mean_curv<=0') # GLOBAL FOLDING PATTERN OF GYRI #axs[1,0].plot(np.sqrt(eigVal/2*np.pi), # coefficients_gyri, color=coefficients_colors_gyri) axs[1,0].scatter(np.sqrt(eigVal/2*np.pi), coefficients_gyri, marker='+', s=10, linewidths=0.5, color=coefficients_colors_gyri) #axs[1,0].scatter(np.sqrt(eigVal[1:]/2*np.pi), # coefficients_gyri[1:], marker='+', s=10, linewidths=0.5, # color=coefficients_colors_gyri[1:]) # remove B0 coefficient axs[1,0].set_xlabel('Frequency (m⁻¹)') axs[1,0].set_ylabel('Coefficients mean_curv>0') axs[1,1].bar(np.arange(0, nlevels), grouped_spectrum_gyri.squeeze(), color=colormap_gyri) #axs[1,1].bar(np.arange(1, nlevels), # grouped_spectrum_gyri[1:].squeeze(), color=colormap_gyri[1:]) # remove B0 axs[1,1].set_xlabel('Spangy frequency bands') axs[1,1].set_ylabel('Power spectrum mean_curv>0') plt.show() # LOCAL SPECTRAL BANDS # -------------------- # Plot of spectral dominant bands on the mesh, with automatized colormap number_of_displayed_bands = len(np.unique(loc_dom_band)) band_values = np.linspace(np.min(loc_dom_band), np.max(loc_dom_band), number_of_displayed_bands+1) # 13 = 6 positive bands * 2 + Band 0 --> to generalize band_colors = np.empty((number_of_displayed_bands+1, 4)) limit = np.max(loc_dom_band) - number_of_high_folding_pattern_bands + 1 # limit band number under which band is displayed in black band_colors[band_values < limit] = colors[0, :] # low frequency bands # associate one color per high frequency sulcus band for i in range(number_of_high_folding_pattern_bands): band_colors[band_values == -(limit+i)] = colors[i+1, :] # associate one color per high frequency gyrus band for i in range(number_of_high_folding_pattern_bands): band_colors[band_values == (limit+i)] = colors[i+1+number_of_high_folding_pattern_bands, :] localbands_colormap = ListedColormap(band_colors) #loc_dom_band = loc_dom_band.astype(int) p = pv.Plotter() p.add_mesh(mesh, scalars=loc_dom_band, show_edges=False, cmap=localbands_colormap, show_scalar_bar=False) p.add_text("Local Dominant Bands", font_size=14) p.add_scalar_bar('Band n°', fmt="%.0f") p.show() ################################################# # DETAILED AND GENERALIZED DISPLAY OPTIONS # ################################################# # AUTOMATIZED COLORMAP DEPENDING ON NUMBER OF BANDS # ------------------------------------------------- # Define automatically the sulci/gyri folding pattern colormap, # in correlation with nlevels --> colors, colormap_sulci, colormap_gyri number_of_high_folding_pattern_bands = int(nlevels/2) colors = np.zeros((2*number_of_high_folding_pattern_bands+1, 4)) # 4 colors for 7 bands (black + B4 + B5 + B6) color_variations_coef = np.linspace(0, 1, number_of_high_folding_pattern_bands) # low frequency bands colors[0, :] = np.array([0, 0, 0, 1]) # to display in grey: np.array([0.75, 0.75, 0.75, 1]) # colors of high-folding-pattern sulci bands for i in range(number_of_high_folding_pattern_bands): colors[i+1, :] = np.array([0, color_variations_coef[i], 1, 1]) # colors of high-folding-pattern gyri bands for i in range(number_of_high_folding_pattern_bands): colors[i+1+number_of_high_folding_pattern_bands, :] \ = np.array([1, color_variations_coef[i], 0, 1]) # allocate dedicated color to each bar of the global power spectrum: # Colormap for Sulci folding patterns colormap_sulci = [] # len = 7 for i in np.arange(0, nlevels-number_of_high_folding_pattern_bands): colormap_sulci.append(colors[0, :]) # B0, B-1, B-2, B-3 for i in np.arange(0, number_of_high_folding_pattern_bands): colormap_sulci.append(colors[i+1, :]) # B-4, B-5, B-6 # Colormap for Gyri folding patterns colormap_gyri = [] # len = 7 for i in np.arange(0, nlevels-number_of_high_folding_pattern_bands): colormap_gyri.append(colors[0, :]) # B0, B1, B2, B3 for i in np.arange(0, number_of_high_folding_pattern_bands): colormap_gyri.append( colors[i+1+number_of_high_folding_pattern_bands, :]) # B4, B5, B6 # GLOBAL COEFFICIENTS AND BANDS # ----------------------------- def plot_global_coefficients_and_bands_sulci_or_gyri( group_indices, coefficients, colormap): coefficients_colors = [] # len = N (e.g. 1500 eigenpairs) # coefficient corresponding to band B0 coefficients_colors.append(colormap[0]) # coefficients corresponding to other bands i.e. B1..B6 & B-1..B-6 band_last_eigenvalue = [] for i in range(len(group_indices)): band_last_eigenvalue.append(group_indices[i][1]) # i, first eigen value of the band Bi for i in np.arange(1, len(coefficients)): # len(coefficients) = N j = 0 while i > band_last_eigenvalue[j]: j = j+1 if i < band_last_eigenvalue[j]: break coefficients_colors.append(colormap[j]) return coefficients_colors .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 17.094 seconds) .. _sphx_glr_download_auto_examples_example_spangy.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example_spangy.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_spangy.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: example_spangy.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_