Add Python examples: embedding + parcellations

Author: OualidBenkarim

LICENSE

@@ -0,0 +1 @@
+from brainspace._version import __version__

brainspace/__init__.py

@@ -0,0 +1,3 @@
+"""BrainSpace version"""
+
+__version__ = '0.1.0'

brainspace/_version.py

@@ -6,9 +6,13 @@
 # License: BSD 3 clause
 
 
+import warnings
+
 import numpy as np
 from scipy.stats import rankdata
+
 from sklearn.metrics.pairwise import cosine_similarity, rbf_kernel
+
 from .utils import dominant_set
 
 
@@ -44,10 +48,15 @@ def compute_affinity(x, kernel=None, sparsity=.9, gamma=None):
         Affinity matrix.
     """
 
-    if sparsity:
+    if sparsity is not None:
         x = dominant_set(x, k=1-sparsity, is_thresh=False, as_sparse=False)
 
     if kernel is None:
+        mask_neg = x < 0
+        if mask_neg.any():
+            x[mask_neg] = 0
+            warnings.warn('The matrix contains negative values and will '
+                          'be zeroed-out.')
         return x
 
     if kernel in {'pearson', 'spearman'}:
@@ -69,4 +78,10 @@ def compute_affinity(x, kernel=None, sparsity=.9, gamma=None):
     else:
         raise ValueError("Unknown kernel '{0}'.".format(kernel))
 
+    mask_neg = a < 0
+    if mask_neg.any():
+        a[mask_neg] = 0
+        warnings.warn('The matrix contains negative values and will '
+                      'be zeroed-out.')
+
     return a

brainspace/gradient/kernels.py

@@ -533,43 +533,89 @@ def propagate_labeling(surf, labeling, no_label=np.nan, mask=None, alpha=0.99,
     return new_labeling
 
 
-def smooth_array(surf, point_data, n_iter=10, mask=None, include_self=True,
-                 kernel='gaussian', sigma=1):
+@append_vtk(to='point')
+def smooth_array(surf, point_data, n_iter=5, mask=None, kernel='gaussian',
+                 relax=0.2, sigma=None, append=False, array_name=None):
+    """Propagate labeling on surface points.
+
+    Parameters
+    ----------
+    surf : vtkPolyData or BSPolyData
+        Input surface.
+    point_data : str, 1D ndarray
+        Input array to smooth. If str, it must be in the point data
+        attributes of `surf`. If ndarray, use this array.
+    n_iter : int, optional
+        Number of smoothing iterations. Default is 5.
+    mask : 1D ndarray, optional
+        Binary mask. If specified, smoothing is only performed on points
+        within the mask. Default is None.
+    kernel : {'uniform', 'gaussian', 'inverse_distance'}, optional
+        Smoothing kernel. Default is 'gaussian'.
+    relax : float, optional
+        Relaxation factor, contribution of neighboring points such that
+        ``0 < relax < 1``. Default is 0.2.
+    sigma :  float, optional
+        Gaussian kernel width. If None, use standard deviation of egde lengths.
+        Default is None.
+    append : bool, optional
+        If True, append array to point data attributes of input surface and
+        return surface. Otherwise, only return array. Default is False.
+    array_name : str, optional
+        Array name to append to surface's point data attributes. Only used if
+        ``append == True``. Default is None.
+
+    Returns
+    -------
+    output : vtkPolyData, BSPolyData or ndarray
+        A 1D array with the smoothed data. Return array if
+        ``append == False``. Otherwise, return input surface with the
+        new array.
+
+    """
+
+    if relax <= 0 or relax >= 1:
+        raise ValueError('Relaxation factor must be between 0 and 1.')
 
     if isinstance(point_data, str):
-        pd = surf.get_array(name=point_data, at='p')
+        point_data = surf.get_array(name=point_data, at='p')
+
+    if mask is not None:
+        pd = point_data[mask]
+    else:
+        pd = point_data
 
-    # adj = me.get_immediate_adjacency(surf)
-    # d = me.get_immediate_distance(surf)
-    # d.setdiag(0)
     if kernel == 'uniform':
-        w = me.get_immediate_adjacency(surf, include_self=include_self)
+        w = me.get_immediate_adjacency(surf, include_self=False, mask=mask,
+                                       dtype=np.float)
     elif kernel == 'gaussian':
-        w = me.get_immediate_distance(surf, metric='sqeuclidean')
-        # w.setdiag(0)
+        w = me.get_immediate_distance(surf, metric='sqeuclidean', mask=mask)
+        if sigma is None:
+            # sigma = w.data.mean() + 3 * w.data.std()
+            sigma = w.data.std()
         w.data *= -.5 / (sigma*sigma)
-        w.data[:] = np.exp(w.data)
-        if include_self:
-            w.setdiag(1)
+        np.exp(w.data, w.data)
     elif kernel == 'inverse_distance':
-        w = me.get_immediate_distance(surf, metric='euclidean')
+        w = me.get_immediate_distance(surf, metric='euclidean', mask=mask)
         w.data **= -1
-        if include_self:
-            w.setdiag(0)
+    else:
+        raise ValueError("Unknown kernel: {0}".format(kernel))
 
-    lam = 1
-    mu = lam
-    norm = w.sum(axis=1)
-    alpha = np.ones(pd.shape)
-    alpha[norm > 0] = 1 - lam
+    w = w.tocoo(copy=False)
+    ws = w.sum(axis=1).A1
+    w.data *= relax/ws[w.row]
 
-    beta = np.zeros(pd.shape)
-    beta[norm > 0] = lam / norm[norm > 0]
+    retain = np.ones(pd.shape)
+    retain[ws > 0] -= relax
 
     spd = pd.copy()
     for i in range(n_iter):
-        spd = alpha * spd + beta * w.dot(spd)
-
-    return spd
+        wp = w.dot(spd)
+        spd *= retain
+        spd += wp
 
+    if mask is not None:
+        spd = map_to_mask(spd, mask=mask)
+        spd[mask] = point_data[mask]
 
+    return spd

brainspace/mesh/array_operations.py

@@ -343,6 +343,7 @@ def get_edges(surf, mask=None):
 
     See Also
     --------
+    :func:`get_edge_length`
     :func:`get_points`
     :func:`get_cells`
 
@@ -355,6 +356,91 @@ def get_edges(surf, mask=None):
     return edges
 
 
+def get_edge_length(surf, metric='euclidean', mask=None):
+    """Get surface edge lengths.
+
+    Parameters
+    ----------
+    surf : vtkDataSet or BSDataSet
+        Input surface.
+    metric : {'euclidean', 'sqeuclidean'}, optional
+        Distance metric. Default is 'euclidean'.
+    mask : 1D ndarray, optional
+        Binary mask. If specified, only use points within the mask.
+        Default is None.
+
+    Returns
+    -------
+    edges : ndarray, shape (n_edges, 2)
+        Array of edges. Each element is a point id.
+
+    See Also
+    --------
+    :func:`get_edges`
+    :func:`get_immediate_distance`
+
+    """
+
+    points = get_points(surf, mask=mask)
+    edges = get_edges(surf, mask=mask)
+
+    dif = points[edges[:, 0]] - points[edges[:, 1]]
+    d = np.einsum('ij,ij->i', dif, dif)
+    if metric == 'euclidean':
+        d **= .5
+    return d
+
+
+def get_border_cells(surf):
+    """Get cells in boundary.
+
+    Cells in boundary have one boundary edge.
+
+    Parameters
+    ----------
+    surf : vtkDataSet or BSDataSet
+        Input surface.
+
+    Returns
+    -------
+    edges : 1D ndarray
+        Array of cells in border.
+
+    See Also
+    --------
+    :func:`get_edges`
+    :func:`get_immediate_distance`
+
+    """
+
+    ce = get_cell_edge_neighbors(surf, include_self=False)
+    return np.where(ce.getnnz(axis=1) < 3)[0]
+
+
+def get_border_edges(surf):
+    """Get edges in border.
+
+    Parameters
+    ----------
+    surf : vtkDataSet or BSDataSet
+        Input surface.
+
+    Returns
+    -------
+    edges : 2D ndarray, shape = (n_edges, 2)
+        Array of edges in border. Each element is a point id.
+
+    See Also
+    --------
+    :func:`get_edges`
+    :func:`get_immediate_distance`
+
+    """
+
+    ce = get_cell_edge_neighbors(surf, include_self=False)
+    return np.where(ce.getnnz(axis=1) < 3)[0]
+
+
 def get_immediate_distance(surf, metric='euclidean', mask=None,
                            dtype=np.float32):
     """Get immediate distance matrix.
@@ -389,8 +475,8 @@ def get_immediate_distance(surf, metric='euclidean', mask=None,
 
     """
 
-    n_pts = surf.GetNumberOfPoints() if mask is None else np.count_nonzero(mask)
     points = get_points(surf, mask=mask)
+    n_pts = points.shape[0]
     edges = get_edges(surf, mask=mask)
 
     dif = points[edges[:, 0]] - points[edges[:, 1]]

brainspace/mesh/mesh_elements.py

@@ -141,7 +141,7 @@ def AddRenderer(self, row=None, col=None, renderer=None, **kwargs):
         if row is None or isinstance(row, tuple):
             row = slice(None) if row is None else slice(*row)
         else:
-            row = slice(row, row+1)
+            row = slice(row, row + 1)
         if col is None or isinstance(col, tuple):
             col = slice(None) if col is None else slice(*col)
         else:
@@ -214,11 +214,11 @@ def show(self, interactive=True, embed_nb=False, scale=None,
         return None
 
     def close(self, *args):
-        try:
-            if hasattr(self, 'panel'):
-                del self.panel
-        except:
-            pass
+        # try:
+        #     if hasattr(self, 'panel'):
+        #         del self.panel
+        # except:
+        #     pass
         self.ren_win.Finalize()
         # self.iren.RemoveAllObservers()
         self.iren.TerminateApp()
@@ -256,8 +256,8 @@ def screenshot(self, filename=None, scale=None, transparent_bg=True):
 
     def Render(self):
         self.ren_win.Render()
-        if hasattr(self, 'panel'):
-            self.panel.param.trigger('object')
+        # if hasattr(self, 'panel'):
+        #     self.panel.param.trigger('object')
 
 
 def _get_qt_app():

brainspace/plotting/base.py

@@ -35,13 +35,12 @@ Authors
 License
 -------
 
-The BrainSpace source code is available under the Revised BSD (3-Clause) license
+The BrainSpace source code is available under the BSD (3-Clause) license.
 
 
 Support
 -------
 
 If you have problems installing the software or questions about usage,
-documentation or something else related to BrainSpace, you can post to our
-mailing list Mail.
-
+documentation or something else related to BrainSpace, you can post to the
+Issues section of our repository.

docs/index.rst

@@ -8,6 +8,7 @@
 # Always prefer setuptools over distutils
 from setuptools import setup, find_packages
 from os import path
+from io import open as io_open
 
 here = path.abspath(path.dirname(__file__))
 
@@ -16,6 +17,11 @@
 with open(path.join(here, 'README.rst'), encoding='utf-8') as f:
     long_description = f.read()
 
+__version__ = None
+version_file = path.join(here, 'brainspace/_version.py')
+with io_open(version_file, mode='r') as fd:
+    exec(fd.read())
+
 # Arguments marked as "Required" below must be included for upload to PyPI.
 # Fields marked as "Optional" may be commented out.
 
@@ -39,7 +45,7 @@
     # For a discussion on single-sourcing the version across setup.py and the
     # project code, see
     # https://packaging.python.org/en/latest/single_source_version.html
-    version='0.1',  # Required
+    version=__version__,  # Required
 
     # This is a one-line description or tagline of what your project does. This
     # corresponds to the "Summary" metadata field:
@@ -76,7 +82,7 @@
 
     # This should be your name or the name of the organization which owns the
     # project.
-    author='The Python Packaging Authority',  # Optional
+    author='BrainSpace developers',  # Optional
 
     # This should be a valid email address corresponding to the author listed
     # above.
@@ -101,7 +107,7 @@
         'Topic :: Scientific/Engineering',
 
         # Pick your license as you wish
-        'License :: OSI Approved :: MIT License',
+        'License :: OSI Approved :: BSD License',
 
         # Specify the Python versions you support here. In particular, ensure
         # that you indicate whether you support Python 2, Python 3 or both.