[MRG] convert surface using python.

.travis.yml

@@ -104,6 +104,7 @@ install:
         python -c 'import mne; mne.datasets.testing.data_path(verbose=True)';
         if [ "${DEPS}" == "full" ]; then
           export FREESURFER_HOME=$(python -c 'import mne; print(mne.datasets.testing.data_path())');
+          export MNE_SKIP_FS_FLASH_CALL=1;
         fi;
       else
         export MNE_SKIP_TESTING_DATASET_TESTS=true;

doc/python_reference.rst

@@ -166,6 +166,7 @@ Functions:
    read_source_spaces
    read_surface
    read_trans
+   read_tri
    save_stc_as_volume
    write_labels_to_annot
    write_bem_solution

mne/__init__.py

@@ -57,7 +57,7 @@
                               spatio_temporal_tris_connectivity,
                               spatio_temporal_dist_connectivity,
                               save_stc_as_volume, extract_label_time_course)
-from .surface import (read_surface, write_surface, decimate_surface,
+from .surface import (read_surface, write_surface, decimate_surface, read_tri,
                       read_morph_map, get_head_surf, get_meg_helmet_surf)
 from .source_space import (read_source_spaces, vertex_to_mni,
                            write_source_spaces, setup_source_space,

mne/bem.py

@@ -1013,7 +1013,7 @@ def make_watershed_bem(subject, subjects_dir=None, overwrite=False,
     -----
     .. versionadded:: 0.10
     """
-    from .surface import read_surface
+    from .surface import read_surface, write_surface, _read_surface_geom
     from .viz.misc import plot_bem
     env, mri_dir = _prepare_env(subject, subjects_dir,
                                 requires_freesurfer=True,
@@ -1063,14 +1063,20 @@ def make_watershed_bem(subject, subjects_dir=None, overwrite=False,
     run_subprocess(cmd, env=env, stdout=sys.stdout, stderr=sys.stderr)
 
     if op.isfile(T1_mgz):
-        # XXX : do this with python code
+        new_info = _extract_volume_info(T1_mgz)
+        if new_info is None:
+            warn('nibabel is required to replace the volume info. Volume info'
+                 'not updated in the written surface.')
+            new_info = dict()
         surfs = ['brain', 'inner_skull', 'outer_skull', 'outer_skin']
         for s in surfs:
             surf_ws_out = op.join(ws_dir, '%s_%s_surface' % (subject, s))
-            cmd = ['mne_convert_surface', '--surf', surf_ws_out, '--mghmri',
-                   T1_mgz, '--surfout', s, "--replacegeom"]
-            run_subprocess(cmd, env=env, stdout=sys.stdout, stderr=sys.stderr)
 
+            surf, volume_info = _read_surface_geom(surf_ws_out,
+                                                   read_metadata=True)
+            volume_info.update(new_info)  # replace volume info, 'head' stays
+
+            write_surface(s, surf['rr'], surf['tris'], volume_info=volume_info)
             # Create symbolic links
             surf_out = op.join(bem_dir, '%s.surf' % s)
             if not overwrite and op.exists(surf_out):
@@ -1114,6 +1120,28 @@ def make_watershed_bem(subject, subjects_dir=None, overwrite=False,
     logger.info('Created %s\n\nComplete.' % (fname_head,))
 
 
+def _extract_volume_info(mgz, raise_error=True):
+    """Helper for extracting volume info from a mgz file."""
+    try:
+        import nibabel as nib
+    except ImportError:
+        return  # warning raised elsewhere
+    header = nib.load(mgz).header
+    new_info = dict()
+    version = header['version']
+    if version == 1:
+        version = '%s  # volume info valid' % version
+    else:
+        raise ValueError('Volume info invalid.')
+    new_info['valid'] = version
+    new_info['filename'] = mgz
+    new_info['volume'] = header['dims'][:3]
+    new_info['voxelsize'] = header['delta']
+    new_info['xras'], new_info['yras'], new_info['zras'] = header['Mdc'].T
+    new_info['cras'] = header['Pxyz_c']
+    return new_info
+
+
 # ############################################################################
 # Read
 
@@ -1642,7 +1670,7 @@ def convert_flash_mris(subject, flash30=True, convert=True, unwarp=False,
 
 @verbose
 def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None,
-                   verbose=None):
+                   flash_path=None, verbose=None):
     """Create 3-Layer BEM model from prepared flash MRI images
 
     Parameters
@@ -1655,30 +1683,40 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None,
         Show surfaces to visually inspect all three BEM surfaces (recommended).
     subjects_dir : string, or None
         Path to SUBJECTS_DIR if it is not set in the environment.
+    flash_path : str | None
+        Path to the flash images. If None (default), mri/flash/parameter_maps
+        within the subject reconstruction is used.
+
+        .. versionadded:: 0.13.0
+
     verbose : bool, str, int, or None
         If not None, override default verbose level (see mne.verbose).
 
     Notes
     -----
-    This program assumes that FreeSurfer and MNE are installed and
-    sourced properly.
+    This program assumes that FreeSurfer is installed and sourced properly.
 
     This function extracts the BEM surfaces (outer skull, inner skull, and
     outer skin) from multiecho FLASH MRI data with spin angles of 5 and 30
     degrees, in mgz format.
 
-    This function assumes that the flash images are available in the
-    folder mri/bem/flash within the freesurfer subject reconstruction.
-
     See Also
     --------
     convert_flash_mris
     """
     from .viz.misc import plot_bem
+    from .surface import write_surface, read_tri
+
+    is_test = os.environ.get('MNE_SKIP_FS_FLASH_CALL', False)
+
     env, mri_dir, bem_dir = _prepare_env(subject, subjects_dir,
                                          requires_freesurfer=True,
-                                         requires_mne=True)
+                                         requires_mne=False)
 
+    if flash_path is None:
+        flash_path = op.join(mri_dir, 'flash', 'parameter_maps')
+    else:
+        flash_path = op.abspath(flash_path)
     curdir = os.getcwd()
     subjects_dir = env['SUBJECTS_DIR']
 
@@ -1690,22 +1728,25 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None,
                                        op.join(bem_dir, 'flash')))
     # Step 4 : Register with MPRAGE
     logger.info("\n---- Registering flash 5 with MPRAGE ----")
-    if not op.exists('flash5_reg.mgz'):
+    flash5 = op.join(flash_path, 'flash5.mgz')
+    flash5_reg = op.join(flash_path, 'flash5_reg.mgz')
+    if not op.exists(flash5_reg):
         if op.exists(op.join(mri_dir, 'T1.mgz')):
             ref_volume = op.join(mri_dir, 'T1.mgz')
         else:
             ref_volume = op.join(mri_dir, 'T1')
-        cmd = ['fsl_rigid_register', '-r', ref_volume, '-i', 'flash5.mgz',
-               '-o', 'flash5_reg.mgz']
+        cmd = ['fsl_rigid_register', '-r', ref_volume, '-i', flash5,
+               '-o', flash5_reg]
         run_subprocess(cmd, env=env, stdout=sys.stdout, stderr=sys.stderr)
     else:
         logger.info("Registered flash 5 image is already there")
     # Step 5a : Convert flash5 into COR
     logger.info("\n---- Converting flash5 volume into COR format ----")
     shutil.rmtree(op.join(mri_dir, 'flash5'), ignore_errors=True)
     os.makedirs(op.join(mri_dir, 'flash5'))
-    cmd = ['mri_convert', 'flash5_reg.mgz', op.join(mri_dir, 'flash5')]
-    run_subprocess(cmd, env=env, stdout=sys.stdout, stderr=sys.stderr)
+    if not is_test:  # CIs don't have freesurfer, skipped when testing.
+        cmd = ['mri_convert', flash5_reg, op.join(mri_dir, 'flash5')]
+        run_subprocess(cmd, env=env, stdout=sys.stdout, stderr=sys.stderr)
     # Step 5b and c : Convert the mgz volumes into COR
     os.chdir(mri_dir)
     convert_T1 = False
@@ -1733,9 +1774,11 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None,
     else:
         logger.info("Brain volume is already in COR format")
     # Finally ready to go
-    logger.info("\n---- Creating the BEM surfaces ----")
-    cmd = ['mri_make_bem_surfaces', subject]
-    run_subprocess(cmd, env=env, stdout=sys.stdout, stderr=sys.stderr)
+    if not is_test:  # CIs don't have freesurfer, skipped when testing.
+        logger.info("\n---- Creating the BEM surfaces ----")
+        cmd = ['mri_make_bem_surfaces', subject]
+        run_subprocess(cmd, env=env, stdout=sys.stdout, stderr=sys.stderr)
+
     logger.info("\n---- Converting the tri files into surf files ----")
     os.chdir(bem_dir)
     if not op.exists('flash'):
@@ -1744,11 +1787,16 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None,
     surfs = ['inner_skull', 'outer_skull', 'outer_skin']
     for surf in surfs:
         shutil.move(op.join(bem_dir, surf + '.tri'), surf + '.tri')
-        cmd = ['mne_convert_surface', '--tri', surf + '.tri', '--surfout',
-               surf + '.surf', '--swap', '--mghmri',
-               op.join(subjects_dir, subject, 'mri', 'flash', 'parameter_maps',
-                       'flash5_reg.mgz')]
-        run_subprocess(cmd, env=env, stdout=sys.stdout, stderr=sys.stderr)
+
+        nodes, tris = read_tri(surf + '.tri', swap=True)
+        vol_info = _extract_volume_info(flash5_reg)
+        if vol_info is None:
+            warn('nibabel is required to update the volume info. Volume info '
+                 'omitted from the written surface.')
+        else:
+            vol_info['head'] = np.array([20])
+        write_surface(surf + '.surf', nodes, tris, volume_info=vol_info)
+
     # Cleanup section
     logger.info("\n---- Cleaning up ----")
     os.chdir(bem_dir)

mne/commands/mne_flash_bem.py

@@ -87,7 +87,7 @@ def run():
     convert_flash_mris(subject=subject, subjects_dir=subjects_dir,
                        flash30=flash30, convert=convert, unwarp=unwarp)
     make_flash_bem(subject=subject, subjects_dir=subjects_dir,
-                   overwrite=overwrite, show=show)
+                   overwrite=overwrite, show=show, flash_path='.')
 
 is_main = (__name__ == '__main__')
 if is_main:

mne/commands/tests/test_commands.py

@@ -212,7 +212,6 @@ def test_watershed_bem():
 
 
 @ultra_slow_test
-@requires_mne
 @requires_freesurfer
 @sample.requires_sample_data
 def test_flash_bem():

mne/datasets/utils.py

@@ -167,7 +167,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True,
     path = _get_path(path, key, name)
     # To update the testing or misc dataset, push commits, then make a new
     # release on GitHub. Then update the "releases" variable:
-    releases = dict(testing='0.23', misc='0.1')
+    releases = dict(testing='0.24', misc='0.1')
     # And also update the "hashes['testing']" variable below.
 
     # To update any other dataset, update the data archive itself (upload
@@ -211,7 +211,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True,
         sample='1d5da3a809fded1ef5734444ab5bf857',
         somato='f3e3a8441477bb5bacae1d0c6e0964fb',
         spm='f61041e3f3f2ba0def8a2ca71592cc41',
-        testing='7ddb41eef97bb2ae6bdb371dc3b5d53d',
+        testing='566d9e97421972c874fe8a76a837c6ed'
     )
     folder_origs = dict(  # not listed means None
         misc='mne-misc-data-%s' % releases['misc'],

mne/fixes.py

@@ -1242,3 +1242,61 @@ def within_tol(x, y, atol, rtol):
     isclose = _isclose
 else:
     isclose = np.isclose
+
+
+def _read_volume_info(fobj):
+    """An implementation of nibabel.freesurfer.io._read_volume_info, since old
+    versions of nibabel (<=2.1.0) don't have it.
+    """
+    volume_info = dict()
+    head = np.fromfile(fobj, '>i4', 1)
+    if not np.array_equal(head, [20]):  # Read two bytes more
+        head = np.concatenate([head, np.fromfile(fobj, '>i4', 2)])
+        if not np.array_equal(head, [2, 0, 20]):
+            warnings.warn("Unknown extension code.")
+            return volume_info
+
+    volume_info['head'] = head
+    for key in ['valid', 'filename', 'volume', 'voxelsize', 'xras', 'yras',
+                'zras', 'cras']:
+        pair = fobj.readline().decode('utf-8').split('=')
+        if pair[0].strip() != key or len(pair) != 2:
+            raise IOError('Error parsing volume info.')
+        if key in ('valid', 'filename'):
+            volume_info[key] = pair[1].strip()
+        elif key == 'volume':
+            volume_info[key] = np.array(pair[1].split()).astype(int)
+        else:
+            volume_info[key] = np.array(pair[1].split()).astype(float)
+    # Ignore the rest
+    return volume_info
+
+
+def _serialize_volume_info(volume_info):
+    """An implementation of nibabel.freesurfer.io._serialize_volume_info, since
+    old versions of nibabel (<=2.1.0) don't have it."""
+    keys = ['head', 'valid', 'filename', 'volume', 'voxelsize', 'xras', 'yras',
+            'zras', 'cras']
+    diff = set(volume_info.keys()).difference(keys)
+    if len(diff) > 0:
+        raise ValueError('Invalid volume info: %s.' % diff.pop())
+
+    strings = list()
+    for key in keys:
+        if key == 'head':
+            if not (np.array_equal(volume_info[key], [20]) or np.array_equal(
+                    volume_info[key], [2, 0, 20])):
+                warnings.warn("Unknown extension code.")
+            strings.append(np.array(volume_info[key], dtype='>i4').tostring())
+        elif key in ('valid', 'filename'):
+            val = volume_info[key]
+            strings.append('{0} = {1}\n'.format(key, val).encode('utf-8'))
+        elif key == 'volume':
+            val = volume_info[key]
+            strings.append('{0} = {1} {2} {3}\n'.format(
+                key, val[0], val[1], val[2]).encode('utf-8'))
+        else:
+            val = volume_info[key]
+            strings.append('{0} = {1:0.10g} {2:0.10g} {3:0.10g}\n'.format(
+                key.ljust(6), val[0], val[1], val[2]).encode('utf-8'))
+    return b''.join(strings)