doc magical 85 factor (=radius of head in mm), improve bvef handling (m…

…ne-tools#6600)

* DOC: magical 85, improve bvef handling

* fix brainvision radius scaling

* update tests in montage for BVEF

* incorporate Alex' feedback

* add whatsnew

* FIX: scale BV coords from mm to m before making the montage

* add test

* fix typo

Co-Authored-By: Alexandre Gramfort <[email protected]>

doc/whats_new.rst

@@ -55,6 +55,8 @@ Changelog
 Bug
 ~~~
 
+- Fix wrong assumptions about units in BrainVision montages and add test asserting units in "mm" or "auto", by `Stefan Appelhoff`_
+
 - Fix scaling issue with signals in mV in EDF files read with :func:`mne.io.read_raw_edf` by `Alex Gramfort`_
 
 - Fix bug in :func:`mne.io.read_raw_brainvision` so that recording date timestamps are also recognized if channel reference data is negative, by `Stefan Appelhoff`_

mne/channels/montage.py

@@ -6,6 +6,7 @@
 #          Jona Sassenhagen <[email protected]>
 #          Teon Brooks <[email protected]>
 #          Christian Brodbeck <[email protected]>
+#          Stefan Appelhoff <[email protected]>
 #
 # License: Simplified BSD
 
@@ -134,9 +135,10 @@ def read_montage(kind, ch_names=None, path=None, unit='m', transform=False):
         mne/channels/data/montages folder in your mne-python installation.
     unit : 'm' | 'cm' | 'mm' | 'auto'
         Unit of the input file. When 'auto' the montage is normalized to
-        a sphere of radius equal to the average brain size. Defaults to 'auto'.
+        a sphere with a radius capturing the average head size (8.5cm).
+        Defaults to 'auto'.
     transform : bool
-        If True, points will be transformed to Neuromag space. The fidicuals,
+        If True, points will be transformed to Neuromag space. The fiducials,
         'nasion', 'lpa', 'rpa' must be specified in the montage file. Useful
         for points captured using Polhemus FastSCAN. Default is False.
 
@@ -278,7 +280,9 @@ def read_montage(kind, ch_names=None, path=None, unit='m', transform=False):
         ch_names_ = data[:, 0].tolist()
         az = np.deg2rad(data[:, 2].astype(float))
         pol = np.deg2rad(data[:, 1].astype(float))
-        pos = _sph_to_cart(np.array([np.ones(len(az)) * 85., az, pol]).T)
+        rad = np.ones(len(az))  # spherical head model
+        rad *= 85.  # scale up to realistic head radius (8.5cm == 85mm)
+        pos = _sph_to_cart(np.array([rad, az, pol]).T)
     elif ext == '.csd':
         # CSD toolbox
         try:  # newer version
@@ -305,7 +309,9 @@ def read_montage(kind, ch_names=None, path=None, unit='m', transform=False):
         az = np.deg2rad(np.array([h if a >= 0. else 180 + h
                                   for h, a in zip(horiz, az)]))
         pol = radius * np.pi
-        pos = _sph_to_cart(np.array([np.ones(len(az)) * 85., az, pol]).T)
+        rad = np.ones(len(az))  # spherical head model
+        rad *= 85.  # scale up to realistic head radius (8.5cm == 85mm)
+        pos = _sph_to_cart(np.array([rad, az, pol]).T)
     elif ext == '.hpts':
         # MNE-C specified format for generic digitizer data
         fid_names = ['1', '2', '3']
@@ -322,19 +328,31 @@ def read_montage(kind, ch_names=None, path=None, unit='m', transform=False):
         pos[:, [0, 1]] = pos[:, [1, 0]] * [-1, 1]
     elif ext == '.bvef':
         # 'BrainVision Electrodes File' format
+        # Based on BrainVision Analyzer coordinate system: Defined between
+        # standard electrode positions: X-axis from T7 to T8, Y-axis from Oz to
+        # Fpz, Z-axis orthogonal from XY-plane through Cz, fit to a sphere if
+        # idealized (when radius=1), specified in millimeters
+        if unit not in ['auto', 'mm']:
+            raise ValueError('`unit` must be "auto" or "mm" for .bvef files.')
         root = ElementTree.parse(fname).getroot()
         ch_names_ = [s.text for s in root.findall("./Electrode/Name")]
         theta = [float(s.text) for s in root.findall("./Electrode/Theta")]
         pol = np.deg2rad(np.array(theta))
         phi = [float(s.text) for s in root.findall("./Electrode/Phi")]
         az = np.deg2rad(np.array(phi))
-        pos = _sph_to_cart(np.array([np.ones(len(az)) * 85., az, pol]).T)
+        rad = [float(s.text) for s in root.findall("./Electrode/Radius")]
+        rad = np.array(rad)  # specified in mm
+        if set(rad) == set([1]):
+            # idealized montage (spherical head model), scale up to realistic
+            # head radius (8.5cm == 85mm)
+            rad = np.array(rad) * 85.
+        pos = _sph_to_cart(np.array([rad, az, pol]).T)
     else:
         raise ValueError('Currently the "%s" template is not supported.' %
                          kind)
     selection = np.arange(len(pos))
 
-    if unit == 'auto':  # rescale to 0.085
+    if unit == 'auto':  # rescale to realistic head radius in meters: 0.085
         pos -= np.mean(pos, axis=0)
         pos = 0.085 * (pos / np.linalg.norm(pos, axis=1).mean())
     elif unit == 'mm':

mne/channels/tests/test_montage.py

@@ -1,4 +1,5 @@
 # Author: Teon Brooks <[email protected]>
+#         Stefan Appelhoff <[email protected]>
 #
 # License: BSD (3-clause)
 
@@ -196,17 +197,18 @@ def test_montage():
         elp=[[-48.20043, 57.55106, 39.86971], [0.0, 60.73848, 59.4629],
              [48.1426, 57.58403, 39.89198], [41.64599, 66.91489, 31.8278]],
         hpts=[[-95, -3, -3], [-1, -1., -3.], [-2, -2, 2.], [0, 0, 0]],
-        bvef=[[-26.266444, 80.839803, 5.204748e-15],
-              [3.680313e-15, 60.104076, 60.104076],
-              [-46.325632, 57.207392, 42.500000],
-              [-68.766444, 49.961746, 5.204748e-15]],
+        bvef=[[-2.62664445e-02,  8.08398039e-02,  5.20474890e-18],
+              [3.68031324e-18,  6.01040764e-02,  6.01040764e-02],
+              [-4.63256329e-02,  5.72073923e-02,  4.25000000e-02],
+              [-6.87664445e-02,  4.99617464e-02,  5.20474890e-18]],
     )
     for key, text in inputs.items():
         kind = key.split('_')[-1]
         fname = op.join(tempdir, 'test.' + kind)
         with open(fname, 'w') as fid:
             fid.write(text)
-        montage = read_montage(fname)
+        unit = 'mm' if kind == 'bvef' else 'm'
+        montage = read_montage(fname, unit=unit)
         if kind in ('sfp', 'txt'):
             assert ('very_very_very_long_name' in montage.ch_names)
         assert_equal(len(montage.ch_names), 4)
@@ -232,6 +234,11 @@ def test_montage():
             assert_array_less(0.01, distance_from_centroid)
         assert_array_almost_equal(poss[key], montage.pos, 4, err_msg=key)
 
+    # Bvef is either auto or mm in terms of "units"
+    with pytest.raises(ValueError, match='be "auto" or "mm" for .bvef files.'):
+        bvef_file = op.join(tempdir, 'test.' + 'bvef')
+        read_montage(bvef_file, unit='m')
+
     # Test reading in different letter case.
     ch_names = ["F3", "FZ", "F4", "FC3", "FCz", "FC4", "C3", "CZ", "C4", "CP3",
                 "CPZ", "CP4", "P3", "PZ", "P4", "O1", "OZ", "O2"]

mne/io/brainvision/brainvision.py

@@ -1,6 +1,5 @@
 # -*- coding: utf-8 -*-
 """Conversion tool from Brain Vision EEG to FIF."""
-
 # Authors: Teon Brooks <[email protected]>
 #          Christian Brodbeck <[email protected]>
 #          Eric Larson <[email protected]>
@@ -525,24 +524,36 @@ def _get_vhdr_info(vhdr_fname, eog, misc, scale, montage):
                               else FIFF.FIFF_UNIT_NONE)
     misc = list(misc_chs.keys()) if misc == 'auto' else misc
 
-    # create montage
+    # create montage: 'Coordinates' section in VHDR file corresponds to "BVEF"
+    # BrainVision Electrode File. The data are based on BrainVision Analyzer
+    # coordinate system: Defined between standard electrode positions: X-axis
+    # from T7 to T8, Y-axis from Oz to Fpz, Z-axis orthogonal from XY-plane
+    # through Cz, fit to a sphere if idealized (when radius=1), specified in mm
     if cfg.has_section('Coordinates') and montage in (None, 'deprecated'):
         from ...transforms import _sph_to_cart
         from ...channels.montage import Montage
         montage_pos = list()
         montage_names = list()
         to_misc = list()
+        # Go through channels
         for ch in cfg.items('Coordinates'):
             ch_name = ch_dict[ch[0]]
             montage_names.append(ch_name)
-            radius, theta, phi = [float(c) for c in ch[1].split(',')]
             # 1: radius, 2: theta, 3: phi
+            rad, theta, phi = [float(c) for c in ch[1].split(',')]
             pol = np.deg2rad(theta)
             az = np.deg2rad(phi)
-            pos = _sph_to_cart(np.array([[radius * 85., az, pol]]))[0]
+            # Coordinates could be "idealized" (spherical head model)
+            if rad == 1:
+                # scale up to realistic head radius (8.5cm == 85mm)
+                rad *= 85.
+            pos = _sph_to_cart(np.array([[rad, az, pol]]))[0]
             if (pos == 0).all() and ch_name not in list(eog) + misc:
                 to_misc.append(ch_name)
             montage_pos.append(pos)
+        # Make a montage, normalizing from BrainVision units "mm" to "m", the
+        # unit used for montages in MNE
+        montage_pos = np.array(montage_pos) / 1e3
         montage_sel = np.arange(len(montage_pos))
         montage = Montage(montage_pos, montage_names, 'Brainvision',
                           montage_sel)

mne/io/brainvision/tests/test_brainvision.py

@@ -1,4 +1,5 @@
 # -*- coding: utf-8 -*-
+"""Test reading of BrainVision format."""
 # Author: Teon Brooks <[email protected]>
 #         Stefan Appelhoff <[email protected]>
 #
@@ -422,6 +423,28 @@ def test_brainvision_vectorized_data():
     assert_array_almost_equal(raw._data[:, :2], first_two_samples_all_chs)
 
 
+def test_coodinates_extraction():
+    """Test reading of [Coordinates] section if present."""
+    # vhdr 2 has a Coordinates section
+    with pytest.warns(RuntimeWarning, match='coordinate information'):
+        raw = read_raw_brainvision(vhdr_v2_path)
+
+    # Basic check of extracted coordinates
+    assert raw.info['dig'] is not None
+    diglist = raw.info['dig']
+    coords = np.array([dig['r'] for dig in diglist])
+    assert coords.shape[0] == len(raw.ch_names)
+    assert coords.shape[1] == 3
+
+    # Make sure the scaling seems right
+    # a coordinate more than 20cm away from origin is implausible
+    assert coords.max() < 0.2
+
+    # vhdr 1 does not have a Coordinates section
+    raw2 = read_raw_brainvision(vhdr_path)
+    assert raw2.info['dig'] is None
+
+
 @testing.requires_testing_data
 def test_brainvision_neuroone_export():
     """Test Brainvision file exported with neuroone system."""
@@ -448,7 +471,7 @@ def test_read_vmrk_annotations():
     try:
         temp.close()
         unlink(temp.name)
-    except FileNotFoundError:
+    except IOError:
         pass
 
 

mne/transforms.py

@@ -689,16 +689,28 @@ def _cart_to_sph(cart):
     return out
 
 
-def _sph_to_cart(sph):
-    """Convert spherical coordinates to Cartesion coordinates."""
-    assert sph.ndim == 2 and sph.shape[1] == 3
-    sph = np.atleast_2d(sph)
-    out = np.empty((len(sph), 3))
-    out[:, 2] = sph[:, 0] * np.cos(sph[:, 2])
-    xy = sph[:, 0] * np.sin(sph[:, 2])
-    out[:, 0] = xy * np.cos(sph[:, 1])
-    out[:, 1] = xy * np.sin(sph[:, 1])
-    return out
+def _sph_to_cart(sph_pts):
+    """Convert spherical coordinates to Cartesion coordinates.
+
+    Parameters
+    ----------
+    sph_pts : ndarray, shape (n_points, 3)
+        Array containing points in spherical coordinates (rad, azimuth, polar)
+
+    Returns
+    -------
+    cart_pts : ndarray, shape (n_points, 3)
+        Array containing points in Cartesian coordinates (x, y, z)
+
+    """
+    assert sph_pts.ndim == 2 and sph_pts.shape[1] == 3
+    sph_pts = np.atleast_2d(sph_pts)
+    cart_pts = np.empty((len(sph_pts), 3))
+    cart_pts[:, 2] = sph_pts[:, 0] * np.cos(sph_pts[:, 2])
+    xy = sph_pts[:, 0] * np.sin(sph_pts[:, 2])
+    cart_pts[:, 0] = xy * np.cos(sph_pts[:, 1])
+    cart_pts[:, 1] = xy * np.sin(sph_pts[:, 1])
+    return cart_pts
 
 
 def _get_n_moments(order):