Repository: MICA-MNI/BrainSpace (language: python)

Issue:
0.1.4 release for VTK pre-release compatibility
A release (or pre-release) on PyPI that includes #75 would be appreciated. We run tests where we `pip install --pre` our whole dependency tree, and currently VTK having a pre-release is breaking with 0.1.3.

This would be easiest for me, but I know it's a chore for you, so if you don't have time to do it, we'll find a different solution.


Discussion / hints:
@OualidBenkarim do you have time for this soon? Otherwise I can have a look mid to late next week.
Completely forgot about this :( It will be done by Monday
Thank you!

Fix patch:
diff --git a/.github/workflows/tagged_release.yml b/.github/workflows/tagged_release.yml
new file mode 100644
index 00000000..ad6d010a
--- /dev/null
+++ b/.github/workflows/tagged_release.yml
@@ -0,0 +1,45 @@
+name: "tagged_release"
+
+on:
+  push:
+    branches:
+      - master
+    tags:
+      - "v*"
+
+jobs:
+  tagged-release:
+    name: "Tagged Release"
+    runs-on: "ubuntu-latest"
+
+    steps:
+    - uses: actions/checkout@v2
+
+    - name: Set up Python 3.8.
+      uses: actions/setup-python@v2
+      with:
+        python-version: 3.8
+
+    - name: Install Python BrainSpace.
+      run: |
+        python -m pip install --upgrade pip
+        python -m pip install --upgrade setuptools wheel
+
+    - name: Build binary wheel and tarball.
+      run: |
+        python setup.py sdist bdist_wheel
+
+    - name: Create Github release.
+      uses: "marvinpinto/action-automatic-releases@latest"
+      with:
+        repo_token: "${{ secrets.GITHUB_TOKEN }}"
+        prerelease: false
+        files: |
+          LICENSE
+          dist/*
+
+    - name: Publish to PyPi.
+      uses: pypa/gh-action-pypi-publish@release/v1
+      with:
+        password: ${{ secrets.PYPI_TOKEN }}
+        verbose: true
\ No newline at end of file
diff --git a/CHANGELOG.txt b/CHANGELOG.txt
index c6aea702..550cd236 100644
--- a/CHANGELOG.txt
+++ b/CHANGELOG.txt
@@ -1,11 +1,16 @@
-v0.1.4, April 2022 -- Bug Fixes
+v0.1.4, August 2022 -- Support for vtk 9.2
 MATLAB:
   - Fixed the uploading of tutorial data to FileExchange.
-  - Fixed a bug that caused the detection of connected affinity matrices to error.
+  - Fixed a bug that caused the detection of connected affinity matrices to
+  error.
   
 Python:
-  - GradientMaps.fit() now returns an error when the affinity matrix contains NaN or Inf.
-  - Added a warning when plot_hemispheres() is run without a display. 
+  - GradientMaps.fit() now returns an error when the affinity matrix contains
+  NaN or Inf.
+  - Added a warning when plot_hemispheres() is run without a display.
+  - Added support for vtk 9.2
+  - For plotting, array_name=None now plots surfaces without any array data.
+
 
 v0.1.3, January 2022 -- Fix dependency error
 Python:
diff --git a/brainspace/_version.py b/brainspace/_version.py
index 5647b8ec..bbb6467a 100644
--- a/brainspace/_version.py
+++ b/brainspace/_version.py
@@ -1,3 +1,3 @@
 """BrainSpace version"""
 
-__version__ = '0.1.3'
+__version__ = '0.1.4'
diff --git a/brainspace/gradient/gradient.py b/brainspace/gradient/gradient.py
index 4e8e68b8..d4765e77 100644
--- a/brainspace/gradient/gradient.py
+++ b/brainspace/gradient/gradient.py
@@ -29,7 +29,7 @@ def _fit_one(x, app, kernel, n_components, random_state, gamma=None,
         Inverse kernel width. Only used if ``kernel`` == 'gaussian'.
         If None, ``gamma=1/n_feat``. Default is None.
     sparsity : float, optional
-        Proportion of smallest elements to zero-out for each row.
+        Proportion of the smallest elements to zero-out for each row.
         Default is 0.9.
     kwargs : kwds, optional
         Additional keyword parameters passed to the embedding approach.
@@ -45,9 +45,9 @@ def _fit_one(x, app, kernel, n_components, random_state, gamma=None,
     a = compute_affinity(x, kernel=kernel, sparsity=sparsity, gamma=gamma)
 
     if np.isnan(a).any() or np.isinf(a).any():
-        raise ValueError(
-            "Affinity matrix contains NaN or Inf values. Common causes of this include NaNs/Infs or rows of zeros in the input matrix."
-        )
+        raise ValueError('Affinity matrix contains NaN or Inf values. Common '
+                         'causes of this include NaNs/Infs or rows of zeros '
+                         'in the input matrix.')
 
     kwds_emb = {'n_components': n_components, 'random_state': random_state}
     kwds_emb.update(kwargs)
@@ -87,7 +87,7 @@ class GradientMaps(BaseEstimator):
         If None, use input matrix. Default is None.
     alignment : {'procrustes', 'joint'}, object or None
         Alignment approach. Only used when two or more datasets are provided.
-        If None, no alignment is peformed. If `object`, it accepts an instance
+        If None, no alignment is performed. If `object`, it accepts an instance
         of :class:`.ProcrustesAlignment`. Default is None.
 
         - If 'procrustes', datasets are aligned using generalized procrustes
@@ -134,7 +134,7 @@ def fit(self, x, gamma=None, sparsity=0.9, n_iter=10, reference=None,
             Inverse kernel width. Only used if ``kernel == 'gaussian'``.
             If None, ``gamma=1/n_feat``. Default is None.
         sparsity : float, optional
-            Proportion of smallest elements to zero-out for each row.
+            Proportion of the smallest elements to zero-out for each row.
             Default is 0.9.
         n_iter : int, optional
             Number of iterations for procrustes alignment. Default is 10.
diff --git a/brainspace/gradient/utils.py b/brainspace/gradient/utils.py
index 4717ce64..199f7b68 100644
--- a/brainspace/gradient/utils.py
+++ b/brainspace/gradient/utils.py
@@ -38,8 +38,8 @@ def is_symmetric(x, tol=1E-10):
     if ssp.issparse(x):
         if x.format not in ['csr', 'csc', 'coo']:
             x = x.tocoo(copy=False)
-        dif1 = x - x.T
-        return np.all(np.abs(dif1.data) < tol)
+        dif = x - x.T
+        return np.all(np.abs(dif.data) < tol)
 
     return np.allclose(x, x.T, atol=tol)
 
@@ -141,7 +141,7 @@ def _dominant_set_dense(s, k, is_thresh=False, norm=False, copy=True):
 
 
 def dominant_set(s, k, is_thresh=False, norm=False, copy=True, as_sparse=True):
-    """Keep largest elements for each row. Zero-out the rest.
+    """Keep the largest elements for each row. Zero-out the rest.
 
     Parameters
     ----------
@@ -183,3 +183,94 @@ def dominant_set(s, k, is_thresh=False, norm=False, copy=True, as_sparse=True):
         return _dominant_set_sparse(s, k, is_thresh=is_thresh, norm=norm)
 
     return _dominant_set_dense(s, k, is_thresh=is_thresh, norm=norm, copy=copy)
+
+
+def ravel_symmetric(x, with_diagonal=False):
+    """Return the flattened upper triangular part of a symmetric matrix.
+
+    Parameters
+    ----------
+    x : 2D ndarray or sparse matrix, shape=(n, n)
+        Input array.
+    with_diagonal : bool, optional
+        If True, also return diagonal elements. Default is False.
+
+    Returns
+    -------
+    output : 1D ndarray, shape (n_feat,)
+        The flattened upper triangular part of `x`. If with_diagonal
+        is True, ``n_feat = n * (n + 1) / 2`` and
+        ``n_feat = n * (n - 1) / 2`` otherwise.
+
+    """
+
+    n = x.shape[0]
+    k = 0 if with_diagonal else -1
+    mask_lt = np.tri(n, k=k, dtype=np.bool)
+
+    if ssp.issparse(x) and not ssp.isspmatrix_csc(x):
+        x = x.tocsc(copy=False)
+
+    return x[mask_lt.T]
+
+
+def unravel_symmetric(x, size, as_sparse=False, part='both', fmt='csr'):
+    """Build symmetric matrix from array with upper triangular elements.
+
+    Parameters
+    ----------
+    x : 1D ndarray
+        Input data with elements to go in the upper triangular part.
+    size : int
+        Number of rows/columns of matrix.
+    as_sparse : bool, optional
+        Return a sparse matrix. Default is False.
+    part: {'both', 'upper', 'lower'}, optional
+        Build matrix with elements if both or just on triangular part.
+        Default is both.
+    fmt: str, optional
+        Format of sparse matrix. Only used if ``as_sparse=True``.
+        Default is 'csr'.
+
+    Returns
+    -------
+    sym : 2D ndarray or sparse matrix, shape = (size, size)
+        Array with the lower/upper or both (symmetric) triangular parts
+        built from `x`.
+
+    """
+
+    k = 1
+    if (size * (size + 1) // 2) == x.size:
+        k = 0
+    elif (size * (size - 1) // 2) != x.size:
+        raise ValueError('Cannot unravel data. Wrong size.')
+
+    shape = (size, size)
+    if as_sparse:
+        mask = x != 0
+        x = x[mask]
+
+        idx = np.triu_indices(size, k=k)
+        idx = [idx1[mask] for idx1 in idx]
+        if part == 'lower':
+            idx = idx[::-1]
+        elif part == 'both':
+            idx = np.concatenate(idx), np.concatenate(idx[::-1])
+            x = np.tile(x, 2)
+
+        xs = ssp.coo_matrix((x, idx), shape=shape)
+        if fmt != 'coo':
+            xs = xs.asformat(fmt, copy=False)
+
+    else:
+        mask_lt = np.tri(size, k=-k, dtype=np.bool)
+        xs = np.zeros(shape, dtype=x.dtype)
+
+        xs[mask_lt.T] = x
+        if part == 'both':
+            xs.T[mask_lt.T] = x
+        elif part == 'lower':
+            xs = xs.T
+
+    return xs
diff --git a/brainspace/mesh/array_operations.py b/brainspace/mesh/array_operations.py
index 7b24e236..35c56e8f 100755
--- a/brainspace/mesh/array_operations.py
+++ b/brainspace/mesh/array_operations.py
@@ -10,7 +10,7 @@
 
 import numpy as np
 from scipy.stats import mode
-from scipy.spatial import cKDTree
+from scipy.spatial import KDTree
 from scipy.sparse.csgraph import laplacian, connected_components
 
 from sklearn.utils.extmath import weighted_mode
@@ -624,31 +624,41 @@ def propagate_labeling(surf, labeling, no_label=np.nan, mask=None, alpha=0.99,
     # Graph matrix
     if mode == 'connectivity':
         adj = me.get_ring_adjacency(surf, n_ring=n_ring, include_self=False,
-                                    dtype=np.float)
+                                    dtype=np.float, mask=mask)
     else:
-        adj = me.get_ring_distance(surf, n_ring=n_ring, dtype=np.float)
+        adj = me.get_ring_distance(surf, n_ring=n_ring, dtype=np.float,
+                                   mask=mask)
         adj.data[:] = np.exp(-adj.data/n_ring**2)
 
     if mask is not None:
         adj = adj[mask][:, mask]
 
-    graph_matrix = -alpha * laplacian(adj, normed=True)
-    diag_mask = (graph_matrix.row == graph_matrix.col)
-    graph_matrix.data[diag_mask] = 0.0
+    # graph_matrix = -alpha * laplacian(adj, normed=True)
+    graph = laplacian(adj, normed=True)
+    graph.data *= -alpha
+    graph.data[graph.row == graph.col] = 0
+    graph.eliminate_zeros()
+    graph = graph.tocsr(copy=False)
+    # diag_mask = (graph_matrix.row == graph_matrix.col)
+    # graph_matrix.data[diag_mask] = 0.0
 
     # Label distributions and label static
     lab_dist = np.zeros((n_pts, n_labs))
     lab_dist[np.argwhere(labeled)[:, 0], idx_lab] = 1
 
-    lab_static = lab_dist.copy()
-    lab_static *= 1 - alpha
+    # lab_static = lab_dist.copy()
+    # lab_static *= 1 - alpha
+    lab_static = (1 - alpha) * lab_dist
 
     # propagation
     lab_dist_perv = lab_dist
     for i in range(n_iter):
-        lab_dist = graph_matrix.dot(lab_dist) + lab_static
+        lab_dist = graph.dot(lab_dist)
+        lab_dist += lab_static
+        # lab_dist = graph_matrix.dot(lab_dist) + lab_static
 
-        if np.linalg.norm(lab_dist - lab_dist_perv, 'fro') < tol:
+        lab_dist_perv -= lab_dist
+        if np.linalg.norm(lab_dist_perv, 'fro') < tol:
             break
 
         lab_dist_perv = lab_dist
@@ -823,15 +833,14 @@ def _get_pids_sphere(source, target, source_mask=None, target_mask=None):
 
 
 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."""
+                    return_weights=True):
+    """Resampling based on the 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)
+    tree = KDTree(sp, leafsize=20)
+    dist, pids = tree.query(tp, k=k, eps=0)
 
     if return_weights:
         return pids, 1 / dist
@@ -899,6 +908,9 @@ def resample_pointdata(source, target, data, is_sphere=False, source_mask=None,
     """
     opt = ['mean', 'mode', 'weighted_mean', 'weighted_mode']
 
+    if n_jobs != 1:
+        warnings.warn('The n_jobs parameter is deprecated and will be removed '
+                      'in a future version', DeprecationWarning)
     is_list = True
     if not isinstance(data, list):
         data = [data]
@@ -918,7 +930,7 @@ def resample_pointdata(source, target, data, is_sphere=False, source_mask=None,
             use_weights = True
 
         pids = _get_pids_naive(source, target, k=k, source_mask=source_mask,
-                               target_mask=target_mask, n_jobs=n_jobs,
+                               target_mask=target_mask,
                                return_weights=use_weights)
         if use_weights:
             pids, w = pids
@@ -942,7 +954,7 @@ def resample_pointdata(source, target, data, is_sphere=False, source_mask=None,
         if k == 1:
             feat = d[pids]
         elif red_func[i] == 'mean':
-            feat = np.mean(d[pids], axis=1)
+            feat = np.nanmean(d[pids], axis=1)
         elif red_func[i] == 'weighted_mean':
             feat = np.average(d[pids], weights=w, axis=1)
         elif red_func[i] == 'mode':
@@ -954,7 +966,7 @@ def resample_pointdata(source, target, data, is_sphere=False, source_mask=None,
             raise ValueError('Unknown red_func: {0}'.format(red_func[i]))
 
         if target_mask is not None:
-            feat = map_to_mask(feat, mask=target_mask, fill=fill)
+            feat = map_to_mask(feat, mask=target_mask, axis=1, fill=fill)
         resampled[i] = feat
 
     if append and key is not None:
diff --git a/brainspace/mesh/mesh_cluster.py b/brainspace/mesh/mesh_cluster.py
index d8051cbc..ab2cc925 100755
--- a/brainspace/mesh/mesh_cluster.py
+++ b/brainspace/mesh/mesh_cluster.py
@@ -90,7 +90,8 @@ def cluster_points(surf, n_clusters=100, is_size=False, mask=None,
     # Find clusters
     if approach == 'kmeans':
         if n_jobs != 1:
-            warn("The n_jobs parameter is deprecated and will be removed in a future version", DeprecationWarning)
+            warn('The n_jobs parameter is deprecated and will be removed '
+                 'in a future version', DeprecationWarning)
         _, cluster_labs, _ = k_means(evs, n_clusters=n_clusters,
                                      random_state=random_state, n_init=n_init)
     else:
diff --git a/brainspace/mesh/mesh_io.py b/brainspace/mesh/mesh_io.py
index 5c0e2d67..5c53995e 100755
--- a/brainspace/mesh/mesh_io.py
+++ b/brainspace/mesh/mesh_io.py
@@ -8,6 +8,7 @@
 import tempfile
 import gzip
 import os
+
 from vtk import (vtkPLYReader, vtkPLYWriter, vtkXMLPolyDataReader,
                  vtkXMLPolyDataWriter, vtkPolyDataReader, vtkPolyDataWriter)
 
@@ -62,6 +63,31 @@ def _select_writer(otype):
     return writer
 
 
+def _uncompress(ipth, opth, block_size=65536):
+    """Uncompresses files. Currently only supports gzip.
+
+    Parameters
+    ----------
+    ipth : str
+        Input filename.
+    opth : str
+        Output filename.
+    block_size : int, optional
+        Size of blocks of the input that are read at a time, by default 65536
+
+    """
+    if ipth.split('.')[-1] == 'gz':
+        with gzip.open(ipth, 'rb') as i_file, open(opth, 'wb') as o_file:
+            while True:
+                block = i_file.read(block_size)
+                if not block:
+                    break
+                else:
+                    o_file.write(block)
+    else:
+        ValueError('Unknown file format.')
+
+
 def read_surface(ipth, itype=None, return_data=True, update=True):
     """Read surface data.
 
@@ -101,11 +127,14 @@ def read_surface(ipth, itype=None, return_data=True, update=True):
 
     if itype == 'gz':
         extension = ipth.split('.')[-2]
-        tmp = tempfile.NamedTemporaryFile(suffix='.' + extension, delete = False)
+        tmp = tempfile.NamedTemporaryFile(suffix='.' + extension, delete=False)
         tmp_name = tmp.name
-        tmp.close() # Close and reopen because windows throws permission errors when both reading and writing.
+        # Close and reopen because windows throws permission errors when both
+        # reading and writing.
+        tmp.close()
         _uncompress(ipth, tmp_name)
-        result = read_surface(tmp_name, extension, return_data=return_data, update=update)
+        result = read_surface(tmp_name, extension, return_data=return_data,
+                              update=update)
         os.unlink(tmp_name)
         return result
 
@@ -178,27 +207,3 @@ def convert_surface(ipth, opth, itype=None, otype=None, oformat=None):
     """
     reader = read_surface(ipth, itype=itype, return_data=False, update=False)
     write_surface(reader, opth, oformat=oformat, otype=otype)
-
-def _uncompress(ipth, opth, block_size=65536):
-    """Uncompresses files. Currently only supports gzip.
-
-    Parameters
-    ----------
-    ipth : str
-        Input filename.
-    opth : str
-        Output filename.
-    block_size : int, optional
-        Size of blocks of the input that are read at a time, by default 65536
-    
-    """
-    if ipth.split('.')[-1] == 'gz':
-        with gzip.open(ipth, 'rb') as i_file, open(opth, 'wb') as o_file:
-            while True:
-                block = i_file.read(block_size)
-                if not block:
-                    break
-                else:
-                    o_file.write(block)
-    else:
-        ValueError('Unknown file format.')
\ No newline at end of file
diff --git a/brainspace/mesh/mesh_operations.py b/brainspace/mesh/mesh_operations.py
index 90102259..40e250c3 100755
--- a/brainspace/mesh/mesh_operations.py
+++ b/brainspace/mesh/mesh_operations.py
@@ -21,7 +21,8 @@
 from ..vtk_interface import wrap_vtk, serial_connect, get_output
 from ..vtk_interface.pipeline import connect
 from ..vtk_interface.decorators import wrap_input, append_vtk
-from ..utils.parcellation import relabel_consecutive, map_to_mask
+from ..utils.parcellation import (relabel_consecutive, map_to_mask,
+                                  reduce_by_labels)
 
 ASSOC_CELLS = vtkDataObject.FIELD_ASSOCIATION_CELLS
 ASSOC_POINTS = vtkDataObject.FIELD_ASSOCIATION_POINTS
@@ -67,6 +68,9 @@ def _surface_selection(surf, array, low=-np.inf, upp=np.inf, use_cell=False):
     if array.ndim > 1:
         raise ValueError('Arrays has more than one dimension.')
 
+    if not use_cell:
+        order_name = surf.append_array(np.arange(surf.n_points), at='p')
+
     if low == -np.inf:
         low = array.min()
     if upp == np.inf:
@@ -95,6 +99,12 @@ def _surface_selection(surf, array, low=-np.inf, upp=np.inf, use_cell=False):
         surf.remove_array(name=array_name, at=at)
         surf_sel.remove_array(name=array_name, at=at)
 
+    if not use_cell:
+        surf_sel = sort_polydata_points(surf_sel, order_name)
+        surf_sel.remove_array(name=order_name, at='p')
+
+        surf.remove_array(name=order_name, at='p')
+
     return surf_sel
 
 
@@ -431,6 +441,36 @@ def get_connected_components(surf, labeling=None, mask=None, fill=0,
     return cc
 
 
+@wrap_input(0)
+def sort_polydata_points(surf, labeling, append_data=True):
+
+    if isinstance(labeling, str):
+        labeling = surf.get_array(labeling, at='p')
+
+    lab_con = relabel_consecutive(labeling)
+
+    idx_sorted = np.argsort(lab_con)
+    new_pts = surf.Points[idx_sorted]
+    # new_cells = relabel(surf.GetCells2D().ravel(), lab_con).reshape(-1, 3)
+    new_cells = lab_con[surf.GetCells2D()]
+    s = build_polydata(new_pts, cells=new_cells)
+
+    if append_data is None or append_data is False:
+        return s
+
+    if append_data is True:
+        append_data = {'p', 'c', 'f'}
+    elif isinstance(append_data, str):
+        append_data = {append_data}
+
+    for at in append_data:
+        for v, k in zip(*surf.get_array(at=at, return_name=True)):
+            if at in {'p', 'point'}:
+                v = v[idx_sorted]
+            s.append_array(v, name=k, at=at)
+    return s
+
+
 @wrap_input(0)
 def split_surface(surf, labeling=None):
     """ Split surface according to the labeling.
@@ -462,15 +502,17 @@ def split_surface(surf, labeling=None):
         labeling = surf.get_array(labeling, at='p')
 
     ulab = np.unique(labeling)
-    return {l: mask_points(surf, labeling == l) for l in ulab}
+    if ulab.size == 1:
+        return {1: surf}
+
+    return {k: mask_points(surf, labeling == k) for k in ulab}
 
 
 @wrap_input(0)
-def downsample_with_parcellation(surf, labeling, name='parcel',
-                                 check_connected=True):
-    """ Downsample surface according to the labeling.
+def downsample_with_parcellation(surf, labeling, name='parcel'):
+    """ Downsample surface according to labeling.
 
-    Such that, each parcel centroid is used as a point in the new donwsampled
+    Such that, each parcel centroid is used as a point in the new downsampled
     surface. Connectivity is based on neighboring parcels.
 
     Parameters
@@ -483,10 +525,6 @@ def downsample_with_parcellation(surf, labeling, name='parcel',
     name : str, optional
         Name of the downsampled parcellation appended to the PointData of the
         new surface. Default is 'parcel'.
-    check_connected : bool, optional
-        Whether to check if the points in each parcel are connected.
-        Downsampling may produce inconsistent results if some parcels have more
-        than one connected component. Default is True.
 
     Returns
     -------
@@ -498,28 +536,22 @@ def downsample_with_parcellation(surf, labeling, name='parcel',
     if isinstance(labeling, str):
         labeling = surf.get_array(labeling, at='p')
 
-    labeling_small = np.unique(labeling)
-    nlabs = labeling_small.size
-
-    labeling_con = relabel_consecutive(labeling)
+    cc = get_connected_components(surf, labeling=labeling) - 1
+    lab_small = reduce_by_labels(labeling, cc, red_op='min')
+    nlabs = lab_small.size
 
     adj = get_immediate_adjacency(surf)
-    adj_neigh = adj.multiply(labeling_con).tocsr()
-
+    adj_neigh = adj.multiply(cc).tocsr()
     adj_small = np.zeros((nlabs, nlabs), dtype=np.bool)
     for i in range(nlabs):
-        arow = adj_neigh[labeling_con == i]
+        arow = adj_neigh[cc == i]
         for j in range(i + 1, nlabs):
             adj_small[j, i] = adj_small[i, j] = np.any(arow.data == j)
 
     points = np.empty((nlabs, 3))
     cells = []
     for i in range(nlabs):
-        m = labeling_con == i
-
-        if check_connected and csg.connected_components(adj[m][:, m])[0] > 1:
-            warnings.warn("Parcel %d is not fully connected. Downsampling may "
-                          "produce inconsistent results." % labeling_small[i])
+        m = cc == i
 
         neigh = np.unique(adj_neigh[m].data)
         neigh = neigh[neigh != i]
@@ -537,5 +569,15 @@ def downsample_with_parcellation(surf, labeling, name='parcel',
 
     cells = np.unique(np.sort(np.vstack(cells), axis=1), axis=0)
     surf_small = build_polydata(points, cells=cells)
-    surf_small.append_array(labeling_small, name=name, at='p')
+    surf_small.append_array(lab_small, name=name, at='p')
+
+    if nlabs == np.unique(labeling).size:
+        return surf_small
+
+    d = split_surface(surf_small)
+    if len(d) > 1:
+        for k, v in d.items():
+            d[k] = sort_polydata_points(v, 'parcel')
+        surf_small = combine_surfaces(*list(d.values()))
+
     return surf_small
diff --git a/brainspace/null_models/spin.py b/brainspace/null_models/spin.py
index 07a474b4..fee64d84 100644
--- a/brainspace/null_models/spin.py
+++ b/brainspace/null_models/spin.py
@@ -9,7 +9,7 @@
 import numpy as np
 from scipy.optimize import linear_sum_assignment
 from scipy.spatial.distance import cdist
-from scipy.spatial import cKDTree
+from scipy.spatial import KDTree
 
 from sklearn.utils import check_random_state
 from sklearn.base import BaseEstimator
@@ -80,7 +80,8 @@ def _generate_spins(points_lh, points_rh=None, unique=False, n_rep=100,
     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()}
+        # tree = {k: cKDTree(p, leafsize=20) for k, p in pts.items()}
+        tree = {k: KDTree(p, leafsize=20) for k, p in pts.items()}
 
     rs = check_random_state(random_state)
 
@@ -105,7 +106,7 @@ def _generate_spins(points_lh, points_rh=None, unique=False, n_rep=100,
                 row, col = linear_sum_assignment(dist)
                 spin[k][i, idx[k]] = idx[k][col]
             else:
-                _, spin[k][i] = tree[k].query(p @ rot[k], k=1, n_jobs=1)
+                _, spin[k][i] = tree[k].query(p @ rot[k], k=1)
 
     return spin
 
@@ -289,11 +290,15 @@ def randomize(self, x_lh, x_rh=None):
         Returns
         -------
         rand_lh : ndarray, shape = (n_rep, n_lh, n_feat)
-            Permutations of `x_rh`. If ``n_feat == 1``, shape = (n_rep, n_lh).
+            Permutations of `x_lh`. If ``n_feat == 1``, shape = (n_rep, n_lh).
         rand_lh : ndarray, shape = (n_rep, n_rh, n_feat)
             Permutations of `x_rh`. If ``n_feat == 1``, shape = (n_rep, n_rh).
             None if `x_rh` is None. Only if `spin_rh_` is not None.
 
+        See Also
+        --------
+        :func:`.SpinPermutations.randomize_gen`
+
         """
 
         rand_lh = x_lh[self.spin_lh_]
@@ -305,3 +310,38 @@ def randomize(self, x_lh, x_rh=None):
             rand_rh = x_rh[self.spin_rh_]
 
         return rand_lh, rand_rh
+
+    def randomize_gen(self, x_lh, x_rh=None):
+        """ Generate random samples from `x_lh` and `x_rh`.
+
+        This is the generator version of :func:`.SpinPermutations.randomize`.
+
+        Parameters
+        ----------
+        x_lh : ndarray, shape = (n_lh,) or (n_lh, n_feat)
+            Array of variables arranged in columns, where `n_feat` is the
+            number of variables.
+        x_rh : ndarray, shape = (n_rh,) or (n_rh, n_feat), optional
+            Array of variables arranged in columns for the right hemisphere.
+            Default is None.
+
+        Yields
+        -------
+        rand_lh : ndarray, shape = (n_lh, n_feat)
+            Permutation of `x_lh`. If ``n_feat == 1``, shape = (n_lh,).
+        rand_lh : ndarray, shape = (n_rh, n_feat)
+            Permutation of `x_rh`. If ``n_feat == 1``, shape = (n_rh,).
+            None if `x_rh` is None. Only if `spin_rh_` is not None.
+
+        See Also
+        --------
+        :func:`.SpinPermutations.randomize`
+        """
+
+        for i in range(self.n_rep):
+            rand_lh = x_lh[self.spin_lh_[i]]
+            if self.spin_rh_ is not None and x_rh is not None:
+                rand_rh = x_rh[self.spin_rh_[i]]
+                yield rand_lh, rand_rh
+            else:
+                yield rand_lh
diff --git a/brainspace/plotting/base.py b/brainspace/plotting/base.py
index d2b93ac8..5463bad5 100644
--- a/brainspace/plotting/base.py
+++ b/brainspace/plotting/base.py
@@ -345,6 +345,10 @@ def close(self):
         del self.ren_win
         self.ren_win = None
         if self.iren:
+            try:
+                self.iren.SetDone(True)
+            except:
+                pass
             self.iren.TerminateApp()
             del self.iren
             self.iren = None
@@ -357,6 +361,10 @@ def quit(self, *args):
         else:
             self.ren_win.Finalize()
             if self.iren:
+                try:
+                    self.iren.SetDone(True)
+                except:
+                    pass
                 self.iren.TerminateApp()
             if self.use_qt:
                 self.app_window.close()
diff --git a/brainspace/plotting/surface_plotting.py b/brainspace/plotting/surface_plotting.py
index af8f8f9b..0ff838fd 100644
--- a/brainspace/plotting/surface_plotting.py
+++ b/brainspace/plotting/surface_plotting.py
@@ -102,9 +102,10 @@ def build_plotter(surfs, layout, array_name=None, view=None, color_bar=None,
     array_name : array-like, optional
         Names of point data array to plot for each layout entry.
         Use a tuple with multiple array names to plot multiple arrays
-        (overlays) per layout entry. Default is None.
+        (overlays) per layout entry. If None, plot surfaces without any array
+        data. Default is None.
     view : array-like, optional
-        View for each each layout entry. Possible views are {'lateral',
+        View for each layout entry. Possible views are {'lateral',
         'medial', 'ventral', 'dorsal'}. If None, use default view.
         Default is None.
     color_bar : {'left', 'right', 'top', 'bottom'} or None, optional
@@ -252,7 +253,11 @@ def build_plotter(surfs, layout, array_name=None, view=None, color_bar=None,
 
         s = surfs[layout[i, j]]
         for ia, name in enumerate(array_name[i, j]):
-            if name is False or name is None:
+            # change to: if None or True, plot surface without array data
+            # if false: do not plot anything
+
+            # if name is False or name is None:
+            if name is False:
                 continue
 
             sp = specs[ia, i, j]
@@ -339,9 +344,10 @@ def plot_surf(surfs, layout, array_name=None, view=None, color_bar=None,
     array_name : array-like, optional
         Names of point data array to plot for each layout entry.
         Use a tuple with multiple array names to plot multiple arrays
-        (overlays) per layout entry. Default is None.
+        (overlays) per layout entry. If None, plot surfaces without any array
+        data. Default is None.
     view : array-like, optional
-        View for each each layout entry. Possible views are {'lateral',
+        View for each layout entry. Possible views are {'lateral',
         'medial', 'ventral', 'dorsal'}. If None, use default view.
         Default is None.
     color_bar : {'left', 'right', 'top', 'bottom'} or None, optional
@@ -384,6 +390,8 @@ def plot_surf(surfs, layout, array_name=None, view=None, color_bar=None,
     scale : tuple, optional
         Scale (magnification). Only used if ``screenshot==True``.
         Default is None.
+    return_plotter: bool, optional
+        If True, return plotter instead of returning image or rendering.
     suppress_warnings : bool, optional
         Whether to suppress warnings. Default is False.
     kwargs : keyword-valued args
@@ -421,17 +429,17 @@ def plot_surf(surfs, layout, array_name=None, view=None, color_bar=None,
     """
 
     if os.name != "nt":
-        # Windows doesn't have a DISPLAY variable... Not sure how to deal with that.
+        # Windows doesn't have a DISPLAY variable...
+        # Not sure how to deal with that.
         # So just skip this check for Windows for now.
-        if ("DISPLAY" not in os.environ or not os.environ["DISPLAY"]) and not suppress_warnings:
+        if ("DISPLAY" not in os.environ or not os.environ["DISPLAY"]) \
+                and not suppress_warnings:
             warnings.warn(
-                (
-                    "Running plot_hemispheres without a display may result in a crash. " +
-                    "For a workaround please consult https://github.com/MICA-MNI/BrainSpace/issues/66. " +
-                    "To suppress this warning set suppress_warnings=True."
-                ), 
-                RuntimeWarning
-            )
+                'Running plot_hemispheres without a display may result in a '
+                'crash. For a workaround please consult '
+                'https://github.com/MICA-MNI/BrainSpace/issues/66. '
+                'To suppress this warning set suppress_warnings=True.',
+                RuntimeWarning)
 
     if screenshot and filename is None:
         raise ValueError('Filename is required.')
@@ -472,7 +480,7 @@ def plot_hemispheres(surf_lh, surf_rh, array_name=None, color_bar=False,
     array_name : str, list of str, ndarray or list of ndarray, optional
         Name of point data array to plot. If ndarray, the array is split for
         the left and right hemispheres. If list, plot one row per array.
-        Default is None.
+        If None, plot surfaces without any array data. Default is None.
     color_bar : bool, optional
         Plot color bar for each array (row). Default is False.
     color_range : {'sym'}, tuple or sequence.
@@ -482,9 +490,9 @@ def plot_hemispheres(surf_lh, surf_rh, array_name=None, color_bar=False,
         Label text for column/row. Possible keys are {'left', 'right',
         'top', 'bottom'}, which indicate the location. Default is None.
     layout_style : str
-        Layout style for hemispheres. If 'row', layout is a single row 
-        alternating lateral and medial views, from left to right. If 'grid', 
-        layout is a 2x2 grid, with lateral views in the top row, medial 
+        Layout style for hemispheres. If 'row', layout is a single row
+        alternating lateral and medial views, from left to right. If 'grid',
+        layout is a 2x2 grid, with lateral views in the top row, medial
         views in the bottom row, and left and right columns. Default is 'row'.
     nan_color : tuple
         Color for nan values. Default is (0, 0, 0, 1).
@@ -514,7 +522,6 @@ def plot_hemispheres(surf_lh, surf_rh, array_name=None, color_bar=False,
     kwargs : keyword-valued args
         Additional arguments passed to the plotter.
 
-
     Returns
     -------
     figure : Ipython Image or None
@@ -539,6 +546,7 @@ def plot_hemispheres(surf_lh, surf_rh, array_name=None, color_bar=False,
         elif array_name.ndim == 1:
             array_name = [array_name]
 
+    to_remove = []
     if isinstance(array_name, list):
         layout = [layout] * len(array_name)
         array_name2 = []
@@ -548,12 +556,13 @@ def plot_hemispheres(surf_lh, surf_rh, array_name=None, color_bar=False,
                 name = surf_lh.append_array(an[:n_pts_lh], at='p')
                 surf_rh.append_array(an[n_pts_lh:], name=name, at='p')
                 array_name2.append(name)
+                to_remove.append(name)
             else:
                 array_name2.append(an)
         array_name = np.asarray(array_name2)[:, None]
 
     if layout_style == 'grid':
-        
+
         # create 2x2 grid for each array_name and stack altogether
         n_arrays = len(array_name)
         array_names, layouts = [], []
@@ -562,7 +571,7 @@ def plot_hemispheres(surf_lh, surf_rh, array_name=None, color_bar=False,
             layouts.append(np.array(l).reshape(2, 2).T.tolist())
         array_name = np.vstack(array_names)
         layout = np.vstack(layouts)
-        
+
         view = [['lateral', 'medial'], ['medial', 'lateral']] * n_arrays
         share = 'both'
     else:
@@ -572,11 +581,24 @@ def plot_hemispheres(surf_lh, surf_rh, array_name=None, color_bar=False,
     if isinstance(cmap, list):
         cmap = np.asarray(cmap)[:, None]
 
+    # when embed_nb=True, interactive (with panel) only supports one renderer,
+    # here we have at least 4
+    if embed_nb:
+        interactive = False
+
     kwds = {'view': view, 'share': share}
     kwds.update(kwargs)
-    return plot_surf(surfs, layout, array_name=array_name, color_bar=color_bar,
-                     color_range=color_range, label_text=label_text, cmap=cmap,
-                     nan_color=nan_color, zoom=zoom, background=background,
-                     size=size, interactive=interactive, embed_nb=embed_nb,
-                     screenshot=screenshot, filename=filename, scale=scale,
-                     transparent_bg=transparent_bg, **kwds)
+    res = plot_surf(surfs, layout, array_name=array_name, color_bar=color_bar,
+                    color_range=color_range, label_text=label_text, cmap=cmap,
+                    nan_color=nan_color, zoom=zoom, background=background,
+                    size=size, interactive=interactive, embed_nb=embed_nb,
+                    screenshot=screenshot, filename=filename, scale=scale,
+                    transparent_bg=transparent_bg, **kwds)
+
+    # remove arrays added to surfaces if any
+    # cannot do it if return_plotter=True
+    if not kwargs.get('return_plotter', False):
+        surf_lh.remove_array(name=to_remove, at='p')
+        surf_rh.remove_array(name=to_remove, at='p')
+
+    return res
diff --git a/brainspace/plotting/utils.py b/brainspace/plotting/utils.py
index 6c232cde..ff260670 100644
--- a/brainspace/plotting/utils.py
+++ b/brainspace/plotting/utils.py
@@ -364,8 +364,29 @@ def _get_ranges(layout, surfs, array_name, share, cbar_range, nvals=256):
         elif share == 'c':
             ax, idx = 1, iter_prod(idx, [slice(None)], range(ncol))
 
-        specs['min'][:] = np.nanmin(specs['min'], axis=ax, keepdims=True)
-        specs['max'][:] = np.nanmax(specs['max'], axis=ax, keepdims=True)
+        # avoid warning when all nan
+        def _nanmin(a):
+            if np.isnan(a).all():
+                return np.nan
+            return np.nanmin(a)
+
+        def _nanmax(a):
+            if np.isnan(a).all():
+                return np.nan
+            return np.nanmax(a)
+
+        if ax == (1, 2):
+            specs['min'][:] = np.atleast_3d(_nanmin(specs['min']))
+            specs['max'][:] = np.atleast_3d(_nanmax(specs['max']))
+        else:
+            min_val = np.apply_along_axis(_nanmin, ax, specs['min'])
+            max_val = np.apply_along_axis(_nanmax, ax, specs['max'])
+            if ax == 1:
+                specs['min'][:] = min_val[None]
+                specs['max'][:] = max_val[None]
+            else:
+                specs['min'][:] = np.atleast_3d(min_val)
+                specs['max'][:] = np.atleast_3d(max_val)
 
         for i in idx:
             if specs['disc'][i].all():
diff --git a/brainspace/utils/parcellation.py b/brainspace/utils/parcellation.py
index 904eb137..a93a27f5 100755
--- a/brainspace/utils/parcellation.py
+++ b/brainspace/utils/parcellation.py
@@ -74,9 +74,9 @@ def find_label_correspondence(lab1, lab2):
 
     Parameters
     ----------
-    lab1 : ndarray, shape = (n_lab,)
+    lab1 : ndarray, shape = (n_labels,)
         First array of labels.
-    lab2 : ndarray, shape = (n_lab,)
+    lab2 : ndarray, shape = (n_labels,)
         Second array of labels.
 
     Returns
@@ -108,14 +108,14 @@ def relabel_by_overlap(lab, ref_lab):
 
     Parameters
     ----------
-    lab : ndarray, shape = (n_lab,)
+    lab : ndarray, shape = (n_labels,)
         Array of labels.
-    ref_lab : ndarray, shape = (n_lab,)
+    ref_lab : ndarray, shape = (n_labels,)
         Reference array of labels.
 
     Returns
     -------
-    new_lab : ndarray, shape = (n_lab,)
+    new_lab : ndarray, shape = (n_labels,)
         Array relabeled using the reference array.
 
     Notes
@@ -147,7 +147,7 @@ def map_to_mask(values, mask, fill=0, axis=0):
 
     Parameters
     ----------
-    values : ndarray, shape = (n_rows, n_cols) or (n_cols,)
+    values : ndarray, shape = (n_values,) or (n_samples, n_values)
         Source array of values.
     mask : ndarray, shape = (n_mask,)
         Mask of boolean values. Data is mapped to mask.
@@ -161,22 +161,23 @@ def map_to_mask(values, mask, fill=0, axis=0):
     -------
     output : ndarray
         Values mapped to mask. If `values` is 1D, shape (n_mask,).
-        When `values` is 2D, shape (n_rows, n_mask) if ``axis == 0`` and
-        (n_mask, n_cols) otherwise.
+        When `values` is 2D, shape (n_samples, n_mask) if ``axis == 0`` and
+        (n_mask, n_samples) otherwise.
 
     """
 
     if np.issubdtype(values.dtype, np.integer) and not np.isfinite(fill):
         raise ValueError("Cannot use non-finite 'fill' with integer arrays.")
 
-    if values.ndim == 1:
-        axis = 0
+    if axis == 1 and values.ndim > 1:
+        values = values.T
 
     values2d = np.atleast_2d(values)
-    n = values2d.shape[axis]
-    mapped = np.full((n, mask.size), fill, dtype=values.dtype)
-    mapped[:, mask] = values2d if axis == 0 else values2d.T
+    if values2d.shape[1] > np.count_nonzero(mask):
+        raise ValueError("Mask cannot allocate values.")
 
+    mapped = np.full((values2d.shape[0], mask.size), fill, dtype=values.dtype)
+    mapped[:, mask] = values2d
     if values.ndim == 1:
         return mapped[0]
     if axis == 1:
@@ -184,7 +185,8 @@ def map_to_mask(values, mask, fill=0, axis=0):
     return mapped
 
 
-def map_to_labels(source_val, target_lab, mask=None, fill=0, source_lab=None):
+def map_to_labels(source_val, target_lab, mask=None, fill=0, axis=0,
+                  source_lab=None):
     """Map data in source to target according to their labels.
 
     Target labels are sorted in ascending order, such that the smallest label
@@ -193,45 +195,67 @@ def map_to_labels(source_val, target_lab, mask=None, fill=0, source_lab=None):
 
     Parameters
     ----------
-    source_val : ndarray, shape = (n_val,)
+    source_val : ndarray, shape = (n_values,) or (n_samples, n_values)
         Source array of values.
-    target_lab : ndarray, shape = (n_lab,)
+    target_lab : ndarray, shape = (n_labels,)
         Target labels.
-    mask : ndarray, shape = (n_lab,), optional
+    mask : int or ndarray, shape = (n_labels,), optional
         If mask is not None, only consider target labels in mask.
+        If int, it indicates the label denoting label to discard.
         Default is None.
     fill : float, optional
         Value used to fill elements outside the mask. Only used if mask is not
         None. Default is 0.
-    source_lab : ndarray, shape = (n_val,), optional
+    axis : {0, 1}, optional
+        If ``axis == 0``, map rows. Otherwise, map columns. Default is 0.
+    source_lab : ndarray, shape = (n_values,), optional
         Source labels for source values. If None, use unique labels in
         `target_lab` in ascending order. Default is None.
 
     Returns
     -------
-    target_val : ndarray, shape = (n_lab,)
+    target_val : ndarray, shape = (n_labels,) or (n_samples, n_labels)
         Target array with corresponding source values.
 
     """
 
     if mask is not None:
-        target_lab2 = target_lab[mask]
-        labs2 = map_to_labels(source_val, target_lab2, source_lab=source_lab)
-        return map_to_mask(labs2, mask, fill=fill)
+        if not isinstance(mask, np.ndarray):
+            mask = target_lab != mask
+
+        target_lab = target_lab[mask]
+        mapped = map_to_labels(source_val, target_lab, axis=axis,
+                               source_lab=source_lab)
+        return map_to_mask(mapped, mask, fill=fill, axis=axis)
 
-    if source_lab is None:
-        uq_tl, idx_tl = np.unique(target_lab, return_inverse=True)
-        return source_val[idx_tl]
+    if axis == 1 and source_val.ndim > 1:
+        source_val = source_val.T
 
-    if source_lab.size != source_val.size:
-        raise ValueError('Source values and labels must have same size.')
+    values2d = np.atleast_2d(source_val)
 
-    uq_sl, idx_sl = np.unique(source_lab, return_inverse=True)
-    if source_lab.size != uq_sl.size:
-        raise ValueError('Source labels must have distinct labels.')
+    ulab, idx = np.unique(target_lab, return_inverse=True)
+    if ulab.size > values2d.shape[1]:
+        raise ValueError('There are more target labels than source values.')
 
-    source_val = source_val[idx_sl]
-    return source_val[target_lab]
+    if source_lab is not None:
+        if source_lab.size != values2d.shape[1]:
+            raise ValueError('Source values and labels must have same size.')
+
+        if not np.isin(ulab, source_lab).all():
+            raise ValueError('Cannot find target labels in source labels.')
+
+        uq_sl, idx_sl = np.unique(source_lab, return_inverse=True)
+        if source_lab.size != uq_sl.size:
+            raise ValueError('Source labels must have distinct labels.')
+
+        values2d = values2d[:, idx_sl]
+
+    mapped = values2d[:, idx]
+    if source_val.ndim == 1:
+        return mapped[0]
+    if axis == 1:
+        return mapped.T
+    return mapped
 
 
 def _get_redop(red_op, weights=None, axis=None):
@@ -259,7 +283,7 @@ def fred(x, w): return np.min(x, axis=axis)
 
 
 def reduce_by_labels(values, labels, weights=None, target_labels=None,
-                     red_op='mean', axis=0, dtype=np.float):
+                     red_op='mean', axis=0, dtype=None):
     """Summarize data in `values` according to `labels`.
 
     Parameters
@@ -282,7 +306,9 @@ def reduce_by_labels(values, labels, weights=None, target_labels=None,
         a 1D array of values, array of weights (or None) and return a scalar
         value. Default is 'mean'.
     dtype : dtype, optional
-        Data type of output array. Default is float.
+        Data type of output array. When None, if `red_op` in
+        {'min', 'max', 'sum', 'mode'}, output is same type as `values`,
+        otherwise output is float. Default is None.
     axis : {0, 1}, optional
         If ``axis == 0``, apply to each row (reduce number of columns per row).
         Otherwise, apply to each column (reduce number of rows per column).
@@ -315,17 +341,22 @@ def reduce_by_labels(values, labels, weights=None, target_labels=None,
     else:
         fred = red_op
 
+    if dtype is None:
+        dtype = np.float
+        if red_op in {'min', 'max', 'sum', 'mode'}:
+            dtype = values.dtype
+
     mapped = np.empty((v2d.shape[0], uq_tl.size), dtype=dtype)
-    for ilab, lab in enumerate(uq_tl):
+    for i, lab in enumerate(uq_tl):
         mask = labels == lab
         wm = None if weights is None else weights[:, mask]
 
         if isinstance(red_op, str):
-            mapped[:, ilab] = fred(v2d[:, mask], wm)
+            mapped[:, i] = fred(v2d[:, mask], wm)
 
         else:
             for idx in range(v2d.shape[0]):
-                mapped[idx, ilab] = fred(v2d[idx, mask], wm)
+                mapped[idx, i] = fred(v2d[idx, mask], wm)
 
     if idx_back is not None:
         mapped = mapped[:, idx_back]
diff --git a/brainspace/vtk_interface/decorators.py b/brainspace/vtk_interface/decorators.py
index dafa893a..63341369 100644
--- a/brainspace/vtk_interface/decorators.py
+++ b/brainspace/vtk_interface/decorators.py
@@ -23,7 +23,7 @@ def wrap_input(*xargs, skip=False):
     ----------
     xargs : sequence of int and str
         Positional indices (integers) and keys as strings (for keyword
-        args) to wrap. If no specified, try to wrap all args.
+        args) to wrap. If none specified, try to wrap all args.
     skip : bool, optional
         Wrap all arguments except those in `xargs`. Default is False.
 
@@ -252,19 +252,19 @@ def append_vtk(to='point'):
 
     Parameters
     ----------
-    to : {'point', 'cell', 'field'}, optional
+    to: {'point', 'cell', 'field'}, optional
         Append data to PointData, CellData or FieldData. Default is 'point'.
 
     Returns
     -------
-    wrapped_func : callable
+    wrapped_func: callable
         Wrapped function.
 
     Notes
     -----
     All functions using this decorator must:
 
-    - Return an ndarray. The size of the array must be consistent with
+    - Return a ndarray. The size of the array must be consistent with
       the data it will be appended to (e.g., number of points if
       ``to == 'point'``), except for FieldData.
 
diff --git a/brainspace/vtk_interface/wrappers/__init__.py b/brainspace/vtk_interface/wrappers/__init__.py
index 9b3cd360..2d6fab62 100644
--- a/brainspace/vtk_interface/wrappers/__init__.py
+++ b/brainspace/vtk_interface/wrappers/__init__.py
@@ -37,7 +37,7 @@
            'BSPolyDataMapper',
            'BSLabeledContourMapper',
            'BSLabeledDataMapper',
-           'BSLabeledPlacementMapper',
+           'BSLabelPlacementMapper',
            'BSPolyDataMapper2D',
            'BSTextMapper2D',
            'BSWindowToImageFilter',
diff --git a/brainspace/vtk_interface/wrappers/algorithm.py b/brainspace/vtk_interface/wrappers/algorithm.py
index 23e1366a..e63e00c3 100644
--- a/brainspace/vtk_interface/wrappers/algorithm.py
+++ b/brainspace/vtk_interface/wrappers/algorithm.py
@@ -69,7 +69,7 @@ def SetLookupTable(self, obj=None, **kwargs):
             Lookup table. If None, a LookupTable is created.
             Default is None.
         kwargs : optional keyword arguments
-            Arguments are use to set the lookup table.
+            Arguments are used to set the lookup table.
         """
 
         if obj is None:
@@ -88,7 +88,7 @@ def SetLookupTableWithEnabling(self, obj=None, **kwargs):
             Lookup table. If None, the lut is created.
             Default is None.
         kwargs : optional keyword arguments
-            Arguments are use to set the lookup table.
+            Arguments are used to set the lookup table.
         """
 
         obj = BSLookupTableWithEnabling(vtkobject=obj, **kwargs)
@@ -104,7 +104,7 @@ def SetWindowLevelLookupTable(self, obj=None, **kwargs):
             Lookup table. If None, the lut is created.
             Default is None.
         kwargs : optional keyword arguments
-            Arguments are use to set the lookup table.
+            Arguments are used to set the lookup table.
         """
 
         obj = BSWindowLevelLookupTable(vtkobject=obj, **kwargs)
@@ -120,7 +120,7 @@ def SetColorTransferFunction(self, obj=None, **kwargs):
             Lookup table. If None, the color transfer function is created.
             Default is None.
         kwargs : optional keyword arguments
-            Arguments are use to set the lookup table.
+            Arguments are used to set the lookup table.
         """
 
         obj = BSColorTransferFunction(vtkobject=obj, **kwargs)
@@ -137,7 +137,7 @@ def SetDiscretizableColorTransferFunction(self, obj=None, **kwargs):
             Lookup table. If None, the discretizable color transfer function
             is created. Default is None.
         kwargs : optional keyword arguments
-            Arguments are use to set the lookup table.
+            Arguments are used to set the lookup table.
         """
 
         obj = BSDiscretizableColorTransferFunction(vtkobject=obj, **kwargs)
@@ -219,7 +219,7 @@ def SetTextProperty(self, obj=None, **kwargs):
             Label text property. If None, the property is created.
             Default is None.
         kwargs : optional keyword arguments
-            Arguments are use to set the property.
+            Arguments are used to set the property.
         """
 
         obj = BSTextProperty(vtkobject=obj, **kwargs)
@@ -272,7 +272,7 @@ def SetTextProperty(self, obj=None, **kwargs):
             Label text property. If None, the property is created.
             Default is None.
         kwargs : optional keyword arguments
-            Arguments are use to set the property.
+            Arguments are used to set the property.
         """
 
         if obj is None:
diff --git a/brainspace/vtk_interface/wrappers/base.py b/brainspace/vtk_interface/wrappers/base.py
index 7fa9fe17..ccea1c2c 100644
--- a/brainspace/vtk_interface/wrappers/base.py
+++ b/brainspace/vtk_interface/wrappers/base.py
@@ -9,12 +9,20 @@
 import numpy as np
 
 import vtk
+from vtk import (vtkAbstractArray, vtkStringArray, vtkIdList, vtkVariantArray,
+                 vtkDataArray)
 from vtk.numpy_interface import dataset_adapter as dsa
 from vtk.util.vtkConstants import VTK_STRING
 
 from .utils import call_vtk, get_vtk_methods, is_numpy_string, is_vtk_string
 
 
+try:
+    from vtk import vtkUnicodeStringArray
+except ImportError:
+    vtkUnicodeStringArray = vtkStringArray
+
+
 class VTKMethodWrapper:
     def __init__(self, name):
         self.name = name
@@ -386,7 +394,8 @@ def BSWrapVTKObject(obj):
 
 def _string_to_numpy(a):
     dtype = np.string_
-    if isinstance(a, vtk.vtkUnicodeStringArray):
+    if isinstance(a, vtkUnicodeStringArray) \
+            and vtkUnicodeStringArray != vtkStringArray:
         dtype = np.unicode_
     shape = a.GetNumberOfTuples(), a.GetNumberOfComponents()
     an = [a.GetValue(i) for i in range(a.GetNumberOfValues())]
@@ -395,9 +404,9 @@ def _string_to_numpy(a):
 
 def _numpy_to_string(a, array_type=None):
     if np.issubdtype(a.dtype, np.string_) or array_type == VTK_STRING:
-        av = vtk.vtkStringArray()
+        av = vtkStringArray()
     else:
-        av = vtk.vtkUnicodeStringArray()
+        av = vtkUnicodeStringArray()
     av.SetNumberOfComponents(1 if a.ndim == 1 else a.shape[1])
     av.SetNumberOfValues(a.size)
     for i, s in enumerate(a.ravel()):
@@ -412,7 +421,7 @@ def _variant_to_numpy(a):
 
 
 def _numpy_to_variant(a):
-    av = vtk.vtkVariantArray()
+    av = vtkVariantArray()
     av.SetNumberOfComponents(1 if a.ndim == 1 else a.shape[1])
     av.SetNumberOfValues(a.size)
     for i, s in enumerate(a.ravel()):
@@ -426,13 +435,13 @@ def _idlist_to_numpy(a):
 
 
 def wrap_vtk_array(a):
-    if isinstance(a, vtk.vtkIdList):
+    if isinstance(a, vtkIdList):
         return _idlist_to_numpy(a)
-    if isinstance(a, (vtk.vtkStringArray, vtk.vtkUnicodeStringArray)):
+    if isinstance(a, (vtkStringArray, vtkUnicodeStringArray)):
         return _string_to_numpy(a)
-    if isinstance(a, vtk.vtkVariantArray):
+    if isinstance(a, vtkVariantArray):
         return _variant_to_numpy(a)
-    if isinstance(a, vtk.vtkDataArray):
+    if isinstance(a, vtkDataArray):
         return dsa.vtkDataArrayToVTKArray(a)
     raise ValueError('Unsupported array type: {0}'.format(type(a)))
 
@@ -459,7 +468,7 @@ def wrap_vtk(obj, **kwargs):
 
     Returns
     -------
-    wrapper : BSVTKObjectWrapper
+    wrapper: BSVTKObjectWrapper
         The wrapped object.
     """
 
@@ -467,7 +476,7 @@ def wrap_vtk(obj, **kwargs):
     if len(kwargs) > 0:
         wobj.setVTK(**kwargs)
 
-    if isinstance(obj, (vtk.vtkAbstractArray, vtk.vtkIdList)):
+    if isinstance(obj, (vtkAbstractArray, vtkIdList)):
         try:
             return wrap_vtk_array(obj)
         except:


Test patch:
diff --git a/brainspace/tests/test_plotting.py b/brainspace/tests/test_plotting.py
index 77d465a2..370f1370 100644
--- a/brainspace/tests/test_plotting.py
+++ b/brainspace/tests/test_plotting.py
@@ -47,14 +47,16 @@ def plotter_multiple_renderers():
     return p
 
 
-@pytest.mark.skipif(ipy is None and pn is None, reason="Requires panel")
+@pytest.mark.skipif(ipy is None or pn is None, reason="Requires panel")
 def test_plotter_panel():
+
+    cls = (pn.pane.VTK, pn.pane.vtk.vtk.VTKRenderWindowSynchronized)
     p = plotter_single_renderer()
     img = p.to_panel()
-    assert isinstance(img, pn.pane.VTK)
+    assert isinstance(img, cls)
 
     img = p.show(embed_nb=True, interactive=True)
-    assert isinstance(img, pn.pane.VTK)
+    assert isinstance(img, cls)
     p.close()
 
     p = plotter_multiple_renderers()