MRG: Add units kwarg (#238)

* ENH: Add units kwarg

* FIX: Tweak for CircleCI

* FIX: Wrong check

* FIX: One more

* FIX: Another

* FIX: Try running another example

Author: larsoner

.circleci/config.yml

@@ -61,7 +61,7 @@ jobs:
               if [[ $(cat build.txt) == "html_dev" ]] || [[ $(cat build.txt) == "html_stable" ]]; then
                 SUBJECTS_DIR=~/mne_data/MNE-sample-data/subjects python -c "import mne; mne.datasets._download_all_example_data()";
               fi;
-        - run: cd doc && sphinx-build -D plot_gallery=1 -D sphinx_gallery_conf.filename_pattern=^\(\(?\!plot_fmri_activation_volume\|plot_morphometry\|plot_label\|plot_probabilistic_label\|plot_resting_correlations\|plot_transparent_brain\|rotate_animation\|save_movie\|save_views\).\)*\$ -b html -d _build/doctrees . _build/html
+        - run: cd doc && sphinx-build -D plot_gallery=1 -D sphinx_gallery_conf.filename_pattern=^\(\(?\!plot_fmri_activation_volume\|plot_morphometry\|plot_label\.py\|plot_probabilistic_label\|plot_resting_correlations\|plot_transparent_brain\|rotate_animation\|save_movie\|save_views\).\)*\$ -b html -d _build/doctrees . _build/html
 
         - store_artifacts:
             path: doc/_build/html/

.gitignore

@@ -27,3 +27,4 @@ pysurfer.egg-info
 examples/example_data/coord-lh.label
 .ipynb_checkpoints/
 .cache/
+.pytest_cache/

examples/plot_label_foci.py

@@ -16,7 +16,8 @@
 """
 Bring up the visualization.
 """
-brain = Brain(subject_id, "lh", "inflated", cortex=("gray", -2, 7, True))
+brain = Brain(subject_id, "lh", "inflated", cortex=("gray", -2, 7, True),
+              units='m')
 
 """
 First we'll identify a stereotaxic focus in the MNI coordinate system. This
@@ -62,4 +63,4 @@
 """
 Set the camera position to show the extent of the labels.
 """
-brain.show_view(dict(elevation=40, distance=430))
+brain.show_view(dict(elevation=40, distance=0.430))

examples/plot_meg_inverse_solution.py

@@ -13,72 +13,59 @@
 
 print(__doc__)
 
-"""
-define subject, surface and hemisphere(s) to plot
-"""
+###############################################################################
+# Set up some useful variables and make the plot.
+
+# define subject, surface and hemisphere(s) to plot:
+
 subject_id, surf = 'fsaverage', 'inflated'
 hemi = 'lh'
 
-"""
-create Brain object for visualization
-"""
+# create Brain object for visualization
 brain = Brain(subject_id, hemi, surf, size=(400, 400), background='w',
-              interaction='terrain', cortex='bone')
+              interaction='terrain', cortex='bone', units='m')
 
-"""
-label for time annotation in milliseconds
-"""
+# label for time annotation in milliseconds
 
 
 def time_label(t):
     return 'time=%0.2f ms' % (t * 1e3)
 
 
-"""
-read MNE dSPM inverse solution
-"""
+# Read MNE dSPM inverse solution and plot
+
 for hemi in ['lh']:  # , 'rh']:
     stc_fname = os.path.join('example_data', 'meg_source_estimate-' +
                              hemi + '.stc')
     stc = read_stc(stc_fname)
 
-    """
-    data and vertices for which the data is defined
-    """
+    # data and vertices for which the data is defined
     data = stc['data']
     vertices = stc['vertices']
 
-    """
-    time points (in seconds)
-    """
+    # time points (in seconds)
     time = np.linspace(stc['tmin'], stc['tmin'] + data.shape[1] * stc['tstep'],
                        data.shape[1], endpoint=False)
 
-    """
-    colormap to use
-    """
+    # 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
     brain.add_data(data, colormap=colormap, vertices=vertices,
                    smoothing_steps=5, time=time, time_label=time_label,
                    hemi=hemi, initial_time=0.1, verbose=False)
 
-"""
-scale colormap
-"""
+# scale colormap
 brain.scale_data_colormap(fmin=13, fmid=18, fmax=22, transparent=True,
                           verbose=False)
 
-"""
-To change the time displayed to 80 ms uncomment this line
-"""
+###############################################################################
+# To change the time displayed to 80 ms uncomment this line:
+
 # brain.set_time(0.08)
 
-"""
-uncomment these lines to use the interactive TimeViewer GUI
-"""
+###############################################################################
+# uncomment these lines to use the interactive TimeViewer GUI
+
 # from surfer import TimeViewer
 # viewer = TimeViewer(brain)

examples/plot_vector_meg_inverse_solution.py

@@ -19,7 +19,7 @@
 subject_id, surf = 'fsaverage', 'white'
 hemi = 'lh'
 brain = Brain(subject_id, hemi, surf, size=(800, 800), interaction='terrain',
-              cortex='0.5', alpha=0.5, show_toolbar=True)
+              cortex='0.5', alpha=0.5, show_toolbar=True, units='m')
 
 # Read the MNE dSPM inverse solution
 

setup.cfg

@@ -1,2 +1,2 @@
 [tool:pytest]
-addopts = --showlocals --durations=10 --doctest-modules -rs --cov --cov-report= --doctest-ignore-import-errors
+addopts = --showlocals --durations=10 --doctest-modules -rs --cov-report= --doctest-ignore-import-errors

surfer/tests/test_utils.py

@@ -51,6 +51,11 @@ def test_surface():
         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 50 < np.linalg.norm(surface.coords, axis=-1).mean() < 100  # mm
+    surface = utils.Surface('fsaverage', 'lh', 'inflated',
+                            subjects_dir=subj_dir, units='m')
+    surface.load_geometry()
+    assert 0.05 < np.linalg.norm(surface.coords, axis=-1).mean() < 0.1  # m
 
 
 def test_huge_cross():

surfer/tests/test_viz.py

@@ -6,6 +6,7 @@
 from tempfile import mkdtemp, mktemp
 import warnings
 
+import pytest
 from nose.tools import assert_equal, assert_in, assert_not_in
 from nose.plugins.skip import SkipTest
 from mayavi import mlab
@@ -240,6 +241,9 @@ def test_label():
     brain.add_label("V2", color="#FF6347", alpha=.6)
     brain.add_label("entorhinal", color=(.2, 1, .5), alpha=.6)
     brain.set_surf('white')
+    brain.show_view(dict(elevation=40, distance=430), distance=430)
+    with pytest.raises(ValueError, match='!='):
+        brain.show_view(dict(elevation=40, distance=430), distance=431)
 
     # remove labels
     brain.remove_labels('V1')

surfer/utils.py

@@ -71,27 +71,18 @@ class Surface(object):
         Which hemisphere to load
     surf : string
         Name of the surface to load (eg. inflated, orig ...)
-    data_path : string
-        Path where to look for data
-    x: 1d array
-        x coordinates of vertices
-    y: 1d array
-        y coordinates of vertices
-    z: 1d array
-        z coordinates of vertices
-    coords : 2d array of shape [n_vertices, 3]
-        The vertices coordinates
-    faces : 2d array
-        The faces ie. the triangles
-    nn : 2d array
-        Normalized surface normals for vertices.
     subjects_dir : str | None
         If not None, this directory will be used as the subjects directory
         instead of the value set using the SUBJECTS_DIR environment variable.
+    offset : float | None
+        If float, align inside edge of each hemisphere to center + offset.
+        If None, do not change coordinates (default).
+    units : str
+        Can be 'm' or 'mm' (default).
     """
 
     def __init__(self, subject_id, hemi, surf, subjects_dir=None,
-                 offset=None):
+                 offset=None, units='mm'):
         """Surface
 
         Parameters
@@ -116,6 +107,7 @@ def __init__(self, subject_id, hemi, surf, subjects_dir=None,
         self.coords = None
         self.faces = None
         self.nn = None
+        self.units = _check_units(units)
 
         subjects_dir = _get_subjects_dir(subjects_dir)
         self.data_path = op.join(subjects_dir, subject_id)
@@ -124,6 +116,8 @@ def load_geometry(self):
         surf_path = op.join(self.data_path, "surf",
                             "%s.%s" % (self.hemi, self.surf))
         coords, faces = nib.freesurfer.read_geometry(surf_path)
+        if self.units == 'm':
+            coords /= 1000.
         if self.offset is not None:
             if self.hemi == 'lh':
                 coords[:, 0] -= (np.max(coords[:, 0]) + self.offset)
@@ -140,11 +134,6 @@ def load_geometry(self):
             self.faces[:] = faces
             self.nn[:] = nn
 
-    def save_geometry(self):
-        surf_path = op.join(self.data_path, "surf",
-                            "%s.%s" % (self.hemi, self.surf))
-        nib.freesurfer.write_geometry(surf_path, self.coords, self.faces)
-
     @property
     def x(self):
         return self.coords[:, 0]
@@ -392,6 +381,12 @@ def dec(*args, **kwargs):
 ###############################################################################
 # USEFUL FUNCTIONS
 
+def _check_units(units):
+    if units not in ('m', 'mm'):
+        raise ValueError('Units must be "m" or "mm", got %r' % (units,))
+    return units
+
+
 def find_closest_vertices(surface_coords, point_coords):
     """Return the vertices on a surface mesh closest to some given coordinates.
 
@@ -414,25 +409,29 @@ def find_closest_vertices(surface_coords, point_coords):
     return np.argmin(cdist(surface_coords, point_coords), axis=0)
 
 
-def tal_to_mni(coords):
+def tal_to_mni(coords, units='mm'):
     """Convert Talairach coords to MNI using the Lancaster transform.
 
     Parameters
     ----------
     coords : n x 3 numpy array
         Array of Talairach coordinates
+    units : str
+        Can be 'm' or 'mm' (default).
 
     Returns
     -------
     mni_coords : n x 3 numpy array
-        Array of coordinates converted to MNI space
-
+        Array of coordinates converted to MNI space.
     """
     coords = np.atleast_2d(coords)
     xfm = np.array([[1.06860, -0.00396, 0.00826,  1.07816],
                     [0.00640,  1.05741, 0.08566,  1.16824],
                     [-0.01281, -0.08863, 1.10792, -4.17805],
                     [0.00000,  0.00000, 0.00000,  1.00000]])
+    units = _check_units(units)
+    if units == 'm':
+        xfm[:3, 3] /= 1000.
     mni_coords = np.dot(np.c_[coords, np.ones(coords.shape[0])], xfm.T)[:, :3]
     return mni_coords
 
@@ -597,7 +596,8 @@ def smoothing_matrix(vertices, adj_mat, smoothing_steps=20, verbose=None):
 
 @verbose
 def coord_to_label(subject_id, coord, label, hemi='lh', n_steps=30,
-                   map_surface='white', coord_as_vert=False, verbose=None):
+                   map_surface='white', coord_as_vert=False, units='mm',
+                   verbose=None):
     """Create label from MNI coordinate
 
     Parameters
@@ -616,18 +616,24 @@ def coord_to_label(subject_id, coord, label, hemi='lh', n_steps=30,
         The surface name used to find the closest point
     coord_as_vert : bool
         whether the coords parameter should be interpreted as vertex ids
+    units : str
+        Can be 'm' or 'mm' (default).
     verbose : bool, str, int, or None
         If not None, override default verbose level (see surfer.verbose).
     """
-    geo = Surface(subject_id, hemi, map_surface)
+    geo = Surface(subject_id, hemi, map_surface, units=units)
     geo.load_geometry()
 
+    coords = geo.coords
+    # work in mm from here on
+    if geo.units == 'm':
+        coords = coords * 1000
     if coord_as_vert:
-        coord = geo.coords[coord]
+        coord = coords[coord]
 
-    n_vertices = len(geo.coords)
+    n_vertices = len(coords)
     adj_mat = mesh_edges(geo.faces)
-    foci_vtxs = find_closest_vertices(geo.coords, [coord])
+    foci_vtxs = find_closest_vertices(coords, [coord])
     data = np.zeros(n_vertices)
     data[foci_vtxs] = 1.
     smooth_mat = smoothing_matrix(np.arange(n_vertices), adj_mat, 1)
@@ -641,7 +647,7 @@ def coord_to_label(subject_id, coord, label, hemi='lh', n_steps=30,
     f.write('#label at %s from subject %s\n' % (coord, subject_id))
     f.write('%d\n' % len(idx))
     for i in idx:
-        x, y, z = geo.coords[i]
+        x, y, z = coords[i]
         f.write('%d  %f  %f  %f 0.000000\n' % (i, x, y, z))
 
 
@@ -696,7 +702,7 @@ def has_fsaverage(subjects_dir=None, raise_error=True, return_why=False):
 
 def has_imageio():
     try:
-        import imageio  # NOQA
+        import imageio  # noqa, analysis:ignore
     except ImportError:
         return False
     else: