How to use the brainspace.mesh.mesh_elements function in brainspace

# Embedding
    evs, _ = diffusion_mapping(a, n_components=30, alpha=0, diffusion_time=1,
                               random_state=random_state)
    evs = normalize(evs)  # To find spherical clusters

    if is_size:
        n_clusters = evs.shape[0] // n_clusters

    # Find clusters
    if approach == 'kmeans':
        _, cluster_labs, _ = k_means(evs, n_clusters=n_clusters,
                                     random_state=random_state, n_jobs=n_jobs,
                                     n_init=n_init)
    else:
        conn = me.get_immediate_adjacency(surf, include_self=False)
        if mask is not None:
            conn = conn[mask, :][:, mask]
        hc = AgglomerativeClustering(n_clusters=n_clusters, connectivity=conn,
                                     linkage='ward')
        cluster_labs = hc.fit(evs).labels_

    # Valid clusters start from 1
    cluster_labs += 1

    # Find centers
    if with_centers:
        points = surf.Points if mask is None else surf.Points[mask]
        centroids = reduce_by_labels(points, cluster_labs, red_op='mean',
                                     axis=1)

        centroid_labs = np.zeros_like(cluster_labs)

References
    ----------
    * Alexander-Bloch A, Shou H, Liu S, Satterthwaite TD, Glahn DC,
      Shinohara RT, Vandekar SN and Raznahan A (2018). On testing for spatial
      correspondence between maps of human brain structure and function.
      NeuroImage, 178:540-51.
    * Blaser R and Fryzlewicz P (2016). Random Rotation Ensembles.
      Journal of Machine Learning Research, 17(4): 1–26.
    * https://netneurotools.readthedocs.io

    """

    # Handle if user provides spheres
    if not isinstance(points_lh, np.ndarray):
        points_lh = me.get_points(points_lh)

    if points_rh is not None:
        if not isinstance(points_rh, np.ndarray):
            points_rh = me.get_points(points_rh)

    pts = {'lh': points_lh}
    if points_rh is not None:
        pts['rh'] = points_rh

        # for reflecting across Y-Z plane
        reflect = np.array([[-1, 0, 0], [0, 1, 0], [0, 0, 1]])

    idx = {k: np.arange(p.shape[0]) for k, p in pts.items()}
    spin = {k: np.empty((n_rep, p.shape[0]), dtype=int)
            for k, p in pts.items()}
    if not unique:

def _get_pids_naive(source, target, k=1, source_mask=None, target_mask=None,
                    return_weights=True, n_jobs=1):
    """Resampling based on k nearest points."""

    sp = me.get_points(source, mask=source_mask)
    tp = me.get_points(target, mask=target_mask)

    tree = cKDTree(sp, leafsize=20, compact_nodes=False, copy_data=False,
                   balanced_tree=False)
    dist, pids = tree.query(tp, k=k, eps=0, n_jobs=n_jobs)

    if return_weights:
        return pids, 1 / dist
    return pids

def _get_pids_naive(source, target, k=1, source_mask=None, target_mask=None,
                    return_weights=True, n_jobs=1):
    """Resampling based on k nearest points."""

    sp = me.get_points(source, mask=source_mask)
    tp = me.get_points(target, mask=target_mask)

    tree = cKDTree(sp, leafsize=20, compact_nodes=False, copy_data=False,
                   balanced_tree=False)
    dist, pids = tree.query(tp, k=k, eps=0, n_jobs=n_jobs)

    if return_weights:
        return pids, 1 / dist
    return pids

Shinohara RT, Vandekar SN and Raznahan A (2018). On testing for spatial
      correspondence between maps of human brain structure and function.
      NeuroImage, 178:540-51.
    * Blaser R and Fryzlewicz P (2016). Random Rotation Ensembles.
      Journal of Machine Learning Research, 17(4): 1–26.
    * https://netneurotools.readthedocs.io

    """

    # Handle if user provides spheres
    if not isinstance(points_lh, np.ndarray):
        points_lh = me.get_points(points_lh)

    if points_rh is not None:
        if not isinstance(points_rh, np.ndarray):
            points_rh = me.get_points(points_rh)

    pts = {'lh': points_lh}
    if points_rh is not None:
        pts['rh'] = points_rh

        # for reflecting across Y-Z plane
        reflect = np.array([[-1, 0, 0], [0, 1, 0], [0, 0, 1]])

    idx = {k: np.arange(p.shape[0]) for k, p in pts.items()}
    spin = {k: np.empty((n_rep, p.shape[0]), dtype=int)
            for k, p in pts.items()}
    if not unique:
        tree = {k: cKDTree(p, leafsize=20) for k, p in pts.items()}

    rs = check_random_state(random_state)

Returns
    -------
    output : vtkPolyData, BSPolyData or ndarray
        Return ndarray if ``append == False``. Otherwise, return input surface
        with the new array.

    """

    if red_func not in ['sum', 'mean', 'mode', 'one_third', 'min', 'max'] and \
            not callable(red_func):
        ValueError('Unknown reduction function \'{0}\'.'.format(red_func))

    if isinstance(cell_data, str):
        cell_data = surf.get_array(name=cell_data, at='c')

    pc = me.get_point2cell_connectivity(surf)
    if isinstance(red_func, str) and red_func != 'mode':

        if red_func in ['sum', 'mean', 'one_third']:
            pd = pc * cell_data
            if red_func == 'mean':
                nnz_row = pc.getnnz(axis=1)
                nnz_row[nnz_row == 0] = 1  # Avoid NaN
                pd = pd / nnz_row
            elif red_func == 'one_third':
                pd = pd / 3
        else:
            pd1 = pc.multiply(cell_data)
            if red_func == 'max':
                pd = np.maximum.reduceat(pd1.data, pc.indptr[:-1])
            else:  # min
                pd = np.minimum.reduceat(pd1.data, pc.indptr[:-1])

append : bool, optional
        If True, append array to cell data attributes of input surface and
        return surface. Otherwise, only return array. Default is False.
    key : str, optional
        Array name to append to surface's point data attributes. Only used if
        ``append == True``. Default is 'point_ncells'.

    Returns
    -------
    output : vtkPolyData, BSPolyData or ndarray
        Return ndarray if ``append == False``. Otherwise, return input surface
        with the new array.

    """

    return me.get_point2cell_connectivity(surf).getnnz(axis=1)