ENH: Add option to smooth to nearest (#275)

* ENH: Add option to smooth to nearest

* DOC: Update docstring

* FIX: Flake

* FIX: Undordered vertices are the worst

* FIX: Orders are hard

* DOC: Changes

Author: larsoner

.gitignore

@@ -5,6 +5,7 @@
 .#*
 *.swp
 *.orig
+*.mov
 build
 
 dist/

CHANGES

@@ -1,5 +1,24 @@
-PySurfer Changes
-================
+Changelog
+=========
+
+.. currentmodule:: surfer
+
+Development version (0.10.dev0)
+-------------------------------
+
+- Added an option to smooth to nearest vertex in :meth:`Brain.add_data` using
+  ``smoothing_steps='nearest'``
+- Added options for using offscreen mode
+- Improved integration with Jupyter notebook
+- Avoided view changes when using :meth:`Brain.add_foci`
+
+Version 0.9
+-----------
+
+- Fixed transparency issues with colormaps with
+  :meth:`Brain.scale_data_colormap`
+- Added an example of using custom colors
+- Added options for choosing units for :class:`Brain` (``m`` or ``mm``)
 
 Version 0.8
 -----------

doc/changes.rst

@@ -0,0 +1 @@
+.. include:: ../CHANGES

doc/index.rst

@@ -19,6 +19,7 @@ More Information
     auto_examples/index.rst
     documentation/index.rst
     python_reference.rst
+    changes.rst
 
 Authors
 -------

examples/plot_meg_inverse_solution.py

@@ -50,9 +50,10 @@ def time_label(t):
     # colormap to use
     colormap = 'hot'
 
-    # add data and set the initial time displayed to 100 ms
+    # add data and set the initial time displayed to 100 ms,
+    # plotted using the nearest relevant colors
     brain.add_data(data, colormap=colormap, vertices=vertices,
-                   smoothing_steps=5, time=time, time_label=time_label,
+                   smoothing_steps='nearest', time=time, time_label=time_label,
                    hemi=hemi, initial_time=0.1, verbose=False)
 
 # scale colormap

surfer/tests/test_utils.py

@@ -1,6 +1,10 @@
+from distutils.version import LooseVersion
 import numpy as np
+import scipy
+from scipy import sparse
+import pytest
 import matplotlib as mpl
-from numpy.testing import assert_array_almost_equal, assert_array_equal
+from numpy.testing import assert_allclose, assert_array_equal
 
 from surfer import utils
 
@@ -44,12 +48,12 @@ def test_surface():
         x = surface.x
         surface.apply_xfm(xfm)
         x_ = surface.x
-        assert_array_almost_equal(x + 2, x_)
+        assert_allclose(x + 2, x_)
 
         # normals
         nn = _slow_compute_normals(surface.coords, surface.faces[:10000])
         nn_fast = utils._compute_normals(surface.coords, surface.faces[:10000])
-        assert_array_almost_equal(nn, nn_fast)
+        assert_allclose(nn, nn_fast)
         assert 50 < np.linalg.norm(surface.coords, axis=-1).mean() < 100  # mm
     surface = utils.Surface('fsaverage', 'lh', 'inflated',
                             subjects_dir=subj_dir, units='m')
@@ -99,3 +103,18 @@ def test_create_color_lut():
     # Test that we can ask for a specific number of colors
     cmap_out = utils.create_color_lut("Reds", 12)
     assert cmap_out.shape == (12, 4)
+
+
+def test_smooth():
+    """Test smoothing support."""
+    adj_mat = sparse.csc_matrix(np.repeat(np.repeat(np.eye(2), 2, 0), 2, 1))
+    vertices = np.array([0, 2])
+    want = np.repeat(np.eye(2), 2, axis=0)
+    smooth = utils.smoothing_matrix(vertices, adj_mat).toarray()
+    assert_allclose(smooth, want)
+    if LooseVersion(scipy.__version__) < LooseVersion('1.3'):
+        with pytest.raises(RuntimeError, match='nearest.*requires'):
+            utils.smoothing_matrix(vertices, adj_mat, 'nearest')
+    else:
+        smooth = utils.smoothing_matrix(vertices, adj_mat, 'nearest').toarray()
+        assert_allclose(smooth, want)

surfer/utils.py

@@ -579,10 +579,36 @@ def smoothing_matrix(vertices, adj_mat, smoothing_steps=20, verbose=None):
     smooth_mat : sparse matrix
         smoothing matrix with size N x len(vertices)
     """
+    if smoothing_steps == 'nearest':
+        mat = _nearest(vertices, adj_mat)
+    else:
+        mat = _smooth(vertices, adj_mat, smoothing_steps)
+    return mat
+
+
+def _nearest(vertices, adj_mat):
+    import scipy
+    from scipy.sparse.csgraph import dijkstra
+    if LooseVersion(scipy.__version__) < LooseVersion('1.3'):
+        raise RuntimeError('smoothing_steps="nearest" requires SciPy >= 1.3')
+    # Vertices can be out of order, so sort them to start ...
+    order = np.argsort(vertices)
+    vertices = vertices[order]
+    _, _, sources = dijkstra(adj_mat, False, indices=vertices, min_only=True,
+                             return_predecessors=True)
+    col = np.searchsorted(vertices, sources)
+    # ... then get things back to the correct configuration.
+    col = order[col]
+    row = np.arange(len(col))
+    data = np.ones(len(col))
+    mat = sparse.coo_matrix((data, (row, col)))
+    assert mat.shape == (adj_mat.shape[0], len(vertices)), mat.shape
+    return mat
+
+
+def _smooth(vertices, adj_mat, smoothing_steps):
     from scipy import sparse
-
     logger.info("Updating smoothing matrix, be patient..")
-
     e = adj_mat.copy()
     e.data[e.data == 2] = 1
     n_vertices = e.shape[0]

surfer/viz.py

@@ -1005,9 +1005,12 @@ def add_data(self, array, min=None, max=None, thresh=None,
             alpha level to control opacity of the overlay.
         vertices : numpy array
             vertices for which the data is defined (needed if len(data) < nvtx)
-        smoothing_steps : int or None
-            number of smoothing steps (smoothing is used if len(data) < nvtx)
-            Default : 20
+        smoothing_steps : int | str | None
+            Number of smoothing steps (if data come from surface subsampling).
+            Can be None to use the fewest steps that result in all vertices
+            taking on data values, or "nearest" such that each high resolution
+            vertex takes the value of the its nearest (on the sphere)
+            low-resolution vertex. Default is 20.
         time : numpy array
             time points in the data array (if data is 2D or 3D)
         time_label : str | callable | None
@@ -2114,13 +2117,17 @@ def data_time_index(self):
         raise RuntimeError("Brain instance has no data overlay")
 
     @verbose
-    def set_data_smoothing_steps(self, smoothing_steps, verbose=None):
+    def set_data_smoothing_steps(self, smoothing_steps=20, verbose=None):
         """Set the number of smoothing steps
 
         Parameters
         ----------
-        smoothing_steps : int
-            Number of smoothing steps
+        smoothing_steps : int | str | None
+            Number of smoothing steps (if data come from surface subsampling).
+            Can be None to use the fewest steps that result in all vertices
+            taking on data values, or "nearest" such that each high resolution
+            vertex takes the value of the its nearest (on the sphere)
+            low-resolution vertex. Default is 20.
         verbose : bool, str, int, or None
             If not None, override default verbose level (see surfer.verbose).
         """