From 90ebbab0879c7a603829be976372039b301dcf3c Mon Sep 17 00:00:00 2001 From: jaeilepp Date: Tue, 31 May 2016 13:52:54 +0200 Subject: [PATCH 01/18] [WIP] convert surface using python. --- mne/bem.py | 111 ++++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 105 insertions(+), 6 deletions(-) diff --git a/mne/bem.py b/mne/bem.py index 79ac99a69dc..b7d867b3fde 100644 --- a/mne/bem.py +++ b/mne/bem.py @@ -16,7 +16,8 @@ from .fixes import partial from .utils import verbose, logger, run_subprocess, get_subjects_dir, warn -from .transforms import _ensure_trans, apply_trans +from .transforms import (_ensure_trans, apply_trans, Transform, + transform_surface_to) from .io import Info from .io.constants import FIFF from .io.write import (start_file, start_block, write_float, write_int, @@ -1063,13 +1064,11 @@ 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 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) + + convert_surface(subject, subjects_dir, surf_ws_out, s, T1_mgz) # Create symbolic links surf_out = op.join(bem_dir, '%s.surf' % s) @@ -1108,12 +1107,112 @@ def make_watershed_bem(subject, subjects_dir=None, overwrite=False, # Show computed BEM surfaces if show: - plot_bem(subject=subject, subjects_dir=subjects_dir, + plot_bem(subject=subject, sub1jects_dir=subjects_dir, orientation='coronal', slices=None, show=True) logger.info('Created %s\n\nComplete.' % (fname_head,)) +def convert_surface(fname_in, fname_out, T1_mgz, useLAS=False): + """Convert surface. + + Parameters + ---------- + fname_in : str + Path to surface in FreeSurfer surf format. + fname_out : str + File to write. + T1_mgz : str + Full path to T1.mgz file. + useLAS : + Flip x axis. + + Returns + ------- + surf : dict + Converted surface. + """ + from .source_space import _get_mgz_header + from .surface import _read_surface_geom, write_surface + header = _get_mgz_header(T1_mgz) + surf_in = _read_surface_geom(fname_in) + surf_in['coord_frame'] = FIFF.FIFFV_COORD_MRI + ex, ey, ez, _ = header['ras2vox'][:3].T + rot = np.zeros((3, 3)) + + ex /= np.linalg.norm(ex) + ey /= np.linalg.norm(ey) + ez /= np.linalg.norm(ez) + right = ant = sup = ex + + if (abs(ey[0]) > abs(right[0])): + right = ey + if (abs(ez[0]) > abs(right[0])): + right = ez + if (abs(ey[1]) > abs(ant[1])): + ant = ey + if (abs(ez[1]) > abs(ant[1])): + ant = ez + if (abs(ey[2]) > abs(sup[2])): + sup = ey + if (abs(ez[2]) > abs(sup[2])): + sup = ez + + if np.all(right == ant) or np.all(right == sup) or np.all(ant == sup): + raise ValueError("Cannot decide the RAS directions.") + + if right[0] < 0: + right *= -1. + if ant[1] < 0: + ant *= -1. + if sup[2] < 0: + sup *= -1. + for j in range(3): + if useLAS: + rot[j][0] = -right[j] + else: + rot[j][0] = right[j] + rot[j][1] = ant[j] + rot[j][2] = sup[j] + + corners = np.zeros((8, 3), dtype=float) + corners[1][0] = header['dims'][0] - 1.0 + corners[2][0] = header['dims'][0] - 1.0 + corners[2][1] = header['dims'][1] - 1.0 + corners[3][1] = header['dims'][1] - 1.0 + for j in range(4, 8): + corners[j][0] = corners[j - 4][0] + corners[j][1] = corners[j - 4][1] + corners[j][2] = header['dims'][2] - 1.0 + + r0 = np.zeros(3, dtype=float) + for j in range(8): + corners[j] = apply_trans(header['ras2vox'], corners[j], move=True) + r0 += corners[j] + r0 /= 8.0 + + dir = (-1.0, -1.0, -1.0) + maxdot = -1. + LPI_RPI = None + for j in range(8): + corners[j] -= r0 + thisdot = np.dot(corners[j], dir) / (np.linalg.norm(corners[j]) * + np.linalg.norm(dir)) + if thisdot > maxdot: + LPI_RPI = corners[j] + maxdot = thisdot + if LPI_RPI is None: + raise ValueError("Could not determine the corner") + + move = LPI_RPI + r0 + trans = np.vstack((np.vstack((rot, move)).T, np.zeros(4))) + trans[3][3] = 1. + trans = Transform(FIFF.FIFFV_COORD_UNKNOWN, FIFF.FIFFV_COORD_MRI, trans) + surf_out = transform_surface_to(surf_in, FIFF.FIFFV_COORD_MRI, trans) + write_surface(fname_out, surf_out['rr'], surf_out['tris']) + return surf_out + + # ############################################################################ # Read From 534cceba435b59b1e53d71d775123cf9488b7d01 Mon Sep 17 00:00:00 2001 From: jaeilepp Date: Fri, 3 Jun 2016 15:50:09 +0200 Subject: [PATCH 02/18] Remade the whole thing. --- mne/bem.py | 160 +++++++++---------------------------------------- mne/surface.py | 40 +++++++++++++ 2 files changed, 68 insertions(+), 132 deletions(-) diff --git a/mne/bem.py b/mne/bem.py index b7d867b3fde..502274e7693 100644 --- a/mne/bem.py +++ b/mne/bem.py @@ -16,8 +16,7 @@ from .fixes import partial from .utils import verbose, logger, run_subprocess, get_subjects_dir, warn -from .transforms import (_ensure_trans, apply_trans, Transform, - transform_surface_to) +from .transforms import _ensure_trans, apply_trans from .io import Info from .io.constants import FIFF from .io.write import (start_file, start_block, write_float, write_int, @@ -1014,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,31 +1062,29 @@ def make_watershed_bem(subject, subjects_dir=None, overwrite=False, os.makedirs(op.join(ws_dir, 'ws')) run_subprocess(cmd, env=env, stdout=sys.stdout, stderr=sys.stderr) - if op.isfile(T1_mgz): - surfs = ['brain', 'inner_skull', 'outer_skull', 'outer_skin'] - for s in surfs: - surf_ws_out = op.join(ws_dir, '%s_%s_surface' % (subject, s)) + surfs = ['brain', 'inner_skull', 'outer_skull', 'outer_skin'] + for s in surfs: + surf_ws_out = op.join(ws_dir, '%s_%s_surface' % (subject, s)) - convert_surface(subject, subjects_dir, surf_ws_out, s, T1_mgz) - - # Create symbolic links - surf_out = op.join(bem_dir, '%s.surf' % s) - if not overwrite and op.exists(surf_out): - skip_symlink = True - else: - if op.exists(surf_out): - os.remove(surf_out) - _symlink(surf_ws_out, surf_out) - skip_symlink = False - - if skip_symlink: - logger.info("Unable to create all symbolic links to .surf files " - "in bem folder. Use --overwrite option to recreate " - "them.") - dest = op.join(bem_dir, 'watershed') + surf = _read_surface_geom(surf_ws_out) + write_surface(s, surf['rr'], surf['tris']) + # Create symbolic links + surf_out = op.join(bem_dir, '%s.surf' % s) + if not overwrite and op.exists(surf_out): + skip_symlink = True else: - logger.info("Symbolic links to .surf files created in bem folder") - dest = bem_dir + if op.exists(surf_out): + os.remove(surf_out) + _symlink(surf_ws_out, surf_out) + skip_symlink = False + + if skip_symlink: + logger.info("Unable to create all symbolic links to .surf files in " + "bem folder. Use --overwrite option to recreate them.") + dest = op.join(bem_dir, 'watershed') + else: + logger.info("Symbolic links to .surf files created in bem folder") + dest = bem_dir logger.info("\nThank you for waiting.\nThe BEM triangulations for this " "subject are now available at:\n%s." % dest) @@ -1107,112 +1104,12 @@ def make_watershed_bem(subject, subjects_dir=None, overwrite=False, # Show computed BEM surfaces if show: - plot_bem(subject=subject, sub1jects_dir=subjects_dir, + plot_bem(subject=subject, subjects_dir=subjects_dir, orientation='coronal', slices=None, show=True) logger.info('Created %s\n\nComplete.' % (fname_head,)) -def convert_surface(fname_in, fname_out, T1_mgz, useLAS=False): - """Convert surface. - - Parameters - ---------- - fname_in : str - Path to surface in FreeSurfer surf format. - fname_out : str - File to write. - T1_mgz : str - Full path to T1.mgz file. - useLAS : - Flip x axis. - - Returns - ------- - surf : dict - Converted surface. - """ - from .source_space import _get_mgz_header - from .surface import _read_surface_geom, write_surface - header = _get_mgz_header(T1_mgz) - surf_in = _read_surface_geom(fname_in) - surf_in['coord_frame'] = FIFF.FIFFV_COORD_MRI - ex, ey, ez, _ = header['ras2vox'][:3].T - rot = np.zeros((3, 3)) - - ex /= np.linalg.norm(ex) - ey /= np.linalg.norm(ey) - ez /= np.linalg.norm(ez) - right = ant = sup = ex - - if (abs(ey[0]) > abs(right[0])): - right = ey - if (abs(ez[0]) > abs(right[0])): - right = ez - if (abs(ey[1]) > abs(ant[1])): - ant = ey - if (abs(ez[1]) > abs(ant[1])): - ant = ez - if (abs(ey[2]) > abs(sup[2])): - sup = ey - if (abs(ez[2]) > abs(sup[2])): - sup = ez - - if np.all(right == ant) or np.all(right == sup) or np.all(ant == sup): - raise ValueError("Cannot decide the RAS directions.") - - if right[0] < 0: - right *= -1. - if ant[1] < 0: - ant *= -1. - if sup[2] < 0: - sup *= -1. - for j in range(3): - if useLAS: - rot[j][0] = -right[j] - else: - rot[j][0] = right[j] - rot[j][1] = ant[j] - rot[j][2] = sup[j] - - corners = np.zeros((8, 3), dtype=float) - corners[1][0] = header['dims'][0] - 1.0 - corners[2][0] = header['dims'][0] - 1.0 - corners[2][1] = header['dims'][1] - 1.0 - corners[3][1] = header['dims'][1] - 1.0 - for j in range(4, 8): - corners[j][0] = corners[j - 4][0] - corners[j][1] = corners[j - 4][1] - corners[j][2] = header['dims'][2] - 1.0 - - r0 = np.zeros(3, dtype=float) - for j in range(8): - corners[j] = apply_trans(header['ras2vox'], corners[j], move=True) - r0 += corners[j] - r0 /= 8.0 - - dir = (-1.0, -1.0, -1.0) - maxdot = -1. - LPI_RPI = None - for j in range(8): - corners[j] -= r0 - thisdot = np.dot(corners[j], dir) / (np.linalg.norm(corners[j]) * - np.linalg.norm(dir)) - if thisdot > maxdot: - LPI_RPI = corners[j] - maxdot = thisdot - if LPI_RPI is None: - raise ValueError("Could not determine the corner") - - move = LPI_RPI + r0 - trans = np.vstack((np.vstack((rot, move)).T, np.zeros(4))) - trans[3][3] = 1. - trans = Transform(FIFF.FIFFV_COORD_UNKNOWN, FIFF.FIFFV_COORD_MRI, trans) - surf_out = transform_surface_to(surf_in, FIFF.FIFFV_COORD_MRI, trans) - write_surface(fname_out, surf_out['rr'], surf_out['tris']) - return surf_out - - # ############################################################################ # Read @@ -1774,6 +1671,7 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, convert_flash_mris """ from .viz.misc import plot_bem + from .surface import write_surface, _load_ascii_surface env, mri_dir, bem_dir = _prepare_env(subject, subjects_dir, requires_freesurfer=True, requires_mne=True) @@ -1843,11 +1741,9 @@ 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) + surf_out = _load_ascii_surface(surf + '.tri', swap=True) + write_surface(surf + '.surf', surf_out[0], surf_out[1]) + # Cleanup section logger.info("\n---- Cleaning up ----") os.chdir(bem_dir) diff --git a/mne/surface.py b/mne/surface.py index a80cec34160..fe09b2bf8c7 100644 --- a/mne/surface.py +++ b/mne/surface.py @@ -1113,3 +1113,43 @@ def mesh_dist(tris, vert): axis=1)) dist_matrix = csr_matrix((dist, (edges.row, edges.col)), shape=edges.shape) return dist_matrix + + +def _load_ascii_surface(filepath, swap=False): + """Function for reading triangle definitions from an ascii file. + Parameters + ---------- + fname_in : str + Path to surface ASCII file (ending with '.tri'). + swap : bool + Assume the ASCII file vertex ordering is clockwise instead of + counterclockwise. + Returns + ------- + surf : tuple (nodes, tris) + The surface.""" + with open(filepath, "r") as fid: + lines = fid.readlines() + n_nodes = int(lines[0]) + n_tris = int(lines[n_nodes + 1]) + n_items = len(lines[1].split()) + if n_items in [3, 6, 14, 17]: + inds = range(3) + elif n_items in [4, 7]: + inds = range(1, 4) + else: + raise IOError('Unrecognized format of data.') + nodes = np.array([np.array([float(v) for v in l.split()])[inds] + for l in lines[1:n_nodes + 1]]) + tris = np.array([np.array([int(v) for v in l.split()])[inds] + for l in lines[n_nodes + 2:n_nodes + 2 + n_tris]]) + if swap: + tris[:, [2, 1]] = tris[:, [1, 2]] + tris -= 1 + logger.info('Loaded surface from %s with %s nodes and %s triangles.' % + (filepath, n_nodes, n_tris)) + if n_items in [3, 4]: + logger.info('Node normals were not included in the source file.') + else: + warn('Node normals were not read.') + return (nodes, tris) From f6e6beb51a1aeaa39c89c9aafb8bbbb2f7cc6890 Mon Sep 17 00:00:00 2001 From: jaeilepp Date: Mon, 6 Jun 2016 15:43:15 +0200 Subject: [PATCH 03/18] Fixes. --- mne/bem.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/mne/bem.py b/mne/bem.py index 502274e7693..5d820171f80 100644 --- a/mne/bem.py +++ b/mne/bem.py @@ -1663,8 +1663,8 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, 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. + This function assumes that the flash images are available in the current + folder. See Also -------- @@ -1769,13 +1769,14 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, os.remove(surf) _symlink(op.join('flash', surf), op.join(surf)) skip_symlink = False - if skip_symlink: - logger.info("Unable to create all symbolic links to .surf files " - "in bem folder. Use --overwrite option to recreate them.") - dest = op.join(bem_dir, 'flash') - else: - logger.info("Symbolic links to .surf files created in bem folder") - dest = bem_dir + if skip_symlink: + logger.info("Unable to create all symbolic links to .surf files " + "in bem folder. Use --overwrite option to recreate " + "them.") + dest = op.join(bem_dir, 'flash') + else: + logger.info("Symbolic links to .surf files created in bem folder") + dest = bem_dir logger.info("\nThank you for waiting.\nThe BEM triangulations for this " "subject are now available at:\n%s.\nWe hope the BEM meshes " "created will facilitate your MEG and EEG data analyses." From 63a628de71f90330d2f38256f3302bcf36004a1a Mon Sep 17 00:00:00 2001 From: jaeilepp Date: Mon, 20 Jun 2016 14:37:23 +0200 Subject: [PATCH 04/18] Read/write volume info --- mne/bem.py | 54 +++++++++++++++++++++++++++++++------------ mne/fixes.py | 58 ++++++++++++++++++++++++++++++++++++++++++++++ mne/surface.py | 62 ++++++++++++++++++++++++++++++++++++++++---------- 3 files changed, 147 insertions(+), 27 deletions(-) diff --git a/mne/bem.py b/mne/bem.py index 5d820171f80..ac062828bbe 100644 --- a/mne/bem.py +++ b/mne/bem.py @@ -1062,21 +1062,26 @@ def make_watershed_bem(subject, subjects_dir=None, overwrite=False, os.makedirs(op.join(ws_dir, 'ws')) run_subprocess(cmd, env=env, stdout=sys.stdout, stderr=sys.stderr) - surfs = ['brain', 'inner_skull', 'outer_skull', 'outer_skin'] - for s in surfs: - surf_ws_out = op.join(ws_dir, '%s_%s_surface' % (subject, s)) - - surf = _read_surface_geom(surf_ws_out) - write_surface(s, surf['rr'], surf['tris']) - # Create symbolic links - surf_out = op.join(bem_dir, '%s.surf' % s) - if not overwrite and op.exists(surf_out): - skip_symlink = True - else: - if op.exists(surf_out): - os.remove(surf_out) - _symlink(surf_ws_out, surf_out) - skip_symlink = False + if op.isfile(T1_mgz): + new_info = _read_volume_info(T1_mgz) + surfs = ['brain', 'inner_skull', 'outer_skull', 'outer_skin'] + for s in surfs: + surf_ws_out = op.join(ws_dir, '%s_%s_surface' % (subject, s)) + + surf, volume_info = _read_surface_geom(surf_ws_out, + read_metadata=True) + volume_info.update(new_info) # replace geometry, '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): + skip_symlink = True + else: + if op.exists(surf_out): + os.remove(surf_out) + _symlink(surf_ws_out, surf_out) + skip_symlink = False if skip_symlink: logger.info("Unable to create all symbolic links to .surf files in " @@ -1110,6 +1115,25 @@ def make_watershed_bem(subject, subjects_dir=None, overwrite=False, logger.info('Created %s\n\nComplete.' % (fname_head,)) +def _read_volume_info(T1): + """Helper for extracting volume info from T1.mgz.""" + import nibabel as nib + header = nib.load(T1).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'] = T1 + 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 diff --git a/mne/fixes.py b/mne/fixes.py index 09651a3ea6f..a8a8d53b143 100644 --- a/mne/fixes.py +++ b/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) diff --git a/mne/surface.py b/mne/surface.py index fe09b2bf8c7..07340e4987e 100644 --- a/mne/surface.py +++ b/mne/surface.py @@ -9,22 +9,24 @@ import sys from struct import pack from glob import glob +from distutils.version import LooseVersion import numpy as np from scipy.sparse import coo_matrix, csr_matrix, eye as speye +import nibabel as nib from .bem import read_bem_surfaces from .io.constants import FIFF from .io.open import fiff_open from .io.tree import dir_tree_find from .io.tag import find_tag -from .io.write import (write_int, start_file, end_block, - start_block, end_file, write_string, - write_float_sparse_rcs) +from .io.write import (write_int, start_file, end_block, start_block, end_file, + write_string, write_float_sparse_rcs) from .channels.channels import _get_meg_system from .transforms import transform_surface_to from .utils import logger, verbose, get_subjects_dir, warn from .externals.six import string_types +from .fixes import _read_volume_info, _serialize_volume_info ############################################################################### @@ -406,7 +408,7 @@ def read_curvature(filepath): @verbose -def read_surface(fname, verbose=None): +def read_surface(fname, read_metadata=False, verbose=None): """Load a Freesurfer surface mesh in triangular format Parameters @@ -428,6 +430,9 @@ def read_surface(fname, verbose=None): -------- write_surface """ + if LooseVersion(nib.__version__) > LooseVersion('2.1.0'): + return nib.freesurfer.read_geometry(fname, read_metadata=read_metadata) + TRIANGLE_MAGIC = 16777214 QUAD_MAGIC = 16777215 NEW_QUAD_MAGIC = 16777213 @@ -462,6 +467,8 @@ def read_surface(fname, verbose=None): fnum = np.fromfile(fobj, ">i4", 1)[0] coords = np.fromfile(fobj, ">f4", vnum * 3).reshape(vnum, 3) faces = np.fromfile(fobj, ">i4", fnum * 3).reshape(fnum, 3) + if read_metadata: + volume_info = _read_volume_info(fobj) else: raise ValueError("%s does not appear to be a Freesurfer surface" % fname) @@ -469,19 +476,25 @@ def read_surface(fname, verbose=None): % (create_stamp.strip(), len(coords), len(faces))) coords = coords.astype(np.float) # XXX: due to mayavi bug on mac 32bits - return coords, faces + + ret = (coords, faces) + if read_metadata: + if len(volume_info) == 0: + warn('No volume information contained in the file') + ret += (volume_info,) + return ret @verbose -def _read_surface_geom(fname, patch_stats=True, norm_rr=False, verbose=None): +def _read_surface_geom(fname, patch_stats=True, norm_rr=False, + read_metadata=False, verbose=None): """Load the surface as dict, optionally add the geometry information""" # based on mne_load_surface_geom() in mne_surface_io.c if isinstance(fname, string_types): - rr, tris = read_surface(fname) # mne_read_triangle_file() - nvert = len(rr) - ntri = len(tris) - s = dict(rr=rr, tris=tris, use_tris=tris, ntri=ntri, - np=nvert) + ret = read_surface(fname, read_metadata=read_metadata) + nvert = len(ret[0]) + ntri = len(ret[1]) + s = dict(rr=ret[0], tris=ret[1], use_tris=ret[1], ntri=ntri, np=nvert) elif isinstance(fname, dict): s = fname else: @@ -490,6 +503,8 @@ def _read_surface_geom(fname, patch_stats=True, norm_rr=False, verbose=None): s = _complete_surface_info(s) if norm_rr is True: _normalize_vectors(s['rr']) + if read_metadata: + return s, ret[2] return s @@ -671,7 +686,7 @@ def _create_surf_spacing(surf, hemi, subject, stype, sval, ico_surf, return surf -def write_surface(fname, coords, faces, create_stamp=''): +def write_surface(fname, coords, faces, create_stamp='', volume_info=None): """Write a triangular Freesurfer surface mesh Accepts the same data format as is returned by read_surface(). @@ -688,6 +703,20 @@ def write_surface(fname, coords, faces, create_stamp=''): create_stamp : str Comment that is written to the beginning of the file. Can not contain line breaks. + volume_info : dict-like or None + Key-value pairs to encode at the end of the file. + Valid keys: + * 'head' : array of int + * 'valid' : str + * 'filename' : str + * 'volume' : array of int, shape (3,) + * 'voxelsize' : array of float, shape (3,) + * 'xras' : array of float, shape (3,) + * 'yras' : array of float, shape (3,) + * 'zras' : array of float, shape (3,) + * 'cras' : array of float, shape (3,) + + .. versionadded:: 0.13.0 See Also -------- @@ -696,6 +725,11 @@ def write_surface(fname, coords, faces, create_stamp=''): if len(create_stamp.splitlines()) > 1: raise ValueError("create_stamp can only contain one line") + if LooseVersion(nib.__version__) > LooseVersion('2.1.0'): + nib.freesurfer.io.write_geometry(fname, coords, faces, + create_stamp=create_stamp, + volume_info=volume_info) + return with open(fname, 'wb') as fid: fid.write(pack('>3B', 255, 255, 254)) strs = ['%s\n' % create_stamp, '\n'] @@ -707,6 +741,10 @@ def write_surface(fname, coords, faces, create_stamp=''): fid.write(np.array(coords, dtype='>f4').tostring()) fid.write(np.array(faces, dtype='>i4').tostring()) + # Add volume info, if given + if volume_info is not None and len(volume_info) > 0: + fid.write(_serialize_volume_info(volume_info)) + ############################################################################### # Decimation From 19909944b32b25e4485d2db1c424415d0450f5cd Mon Sep 17 00:00:00 2001 From: jaeilepp Date: Mon, 20 Jun 2016 15:54:58 +0200 Subject: [PATCH 05/18] Tests --- mne/bem.py | 2 +- mne/surface.py | 24 +++++++++++++++++++++--- mne/tests/test_surface.py | 9 ++++++--- 3 files changed, 28 insertions(+), 7 deletions(-) diff --git a/mne/bem.py b/mne/bem.py index ac062828bbe..ba8760129fc 100644 --- a/mne/bem.py +++ b/mne/bem.py @@ -1070,7 +1070,7 @@ def make_watershed_bem(subject, subjects_dir=None, overwrite=False, surf, volume_info = _read_surface_geom(surf_ws_out, read_metadata=True) - volume_info.update(new_info) # replace geometry, 'head' stays + volume_info.update(new_info) # replace volume info, 'head' stays write_surface(s, surf['rr'], surf['tris'], volume_info=volume_info) # Create symbolic links diff --git a/mne/surface.py b/mne/surface.py index 07340e4987e..4b724cfff33 100644 --- a/mne/surface.py +++ b/mne/surface.py @@ -417,6 +417,20 @@ def read_surface(fname, read_metadata=False, verbose=None): The name of the file containing the surface. verbose : bool, str, int, or None If not None, override default verbose level (see mne.verbose). + read_metadata : bool + Read metadata as key-value pairs. + Valid keys: + * 'head' : array of int + * 'valid' : str + * 'filename' : str + * 'volume' : array of int, shape (3,) + * 'voxelsize' : array of float, shape (3,) + * 'xras' : array of float, shape (3,) + * 'yras' : array of float, shape (3,) + * 'zras' : array of float, shape (3,) + * 'cras' : array of float, shape (3,) + + .. versionadded:: 0.13.0 Returns ------- @@ -425,14 +439,18 @@ def read_surface(fname, read_metadata=False, verbose=None): tris : int array, shape=(n_faces, 3) Triangulation (each line contains indexes for three points which together form a face). - + volume_info : dict-like + If read_metadata is true, key-value pairs found in the geometry file. See Also -------- write_surface """ - if LooseVersion(nib.__version__) > LooseVersion('2.1.0'): - return nib.freesurfer.read_geometry(fname, read_metadata=read_metadata) + # XXX: Tests fail here due to numerical error. + # if LooseVersion(nib.__version__) > LooseVersion('2.1.0'): + # return nib.freesurfer.read_geometry(fname, + # read_metadata=read_metadata) + volume_info = dict() TRIANGLE_MAGIC = 16777214 QUAD_MAGIC = 16777215 NEW_QUAD_MAGIC = 16777213 diff --git a/mne/tests/test_surface.py b/mne/tests/test_surface.py index 7394882eefd..f8cba010be9 100644 --- a/mne/tests/test_surface.py +++ b/mne/tests/test_surface.py @@ -16,6 +16,7 @@ from mne.utils import _TempDir, requires_mayavi, run_tests_if_main, slow_test from mne.io import read_info from mne.transforms import _get_trans +from mne.io.meas_info import _is_equal_dict data_path = testing.data_path(download=False) subjects_dir = op.join(data_path, 'subjects') @@ -130,11 +131,13 @@ def test_io_surface(): fname_tri = op.join(data_path, 'subjects', 'fsaverage', 'surf', 'lh.inflated') for fname in (fname_quad, fname_tri): - pts, tri = read_surface(fname) - write_surface(op.join(tempdir, 'tmp'), pts, tri) - c_pts, c_tri = read_surface(op.join(tempdir, 'tmp')) + pts, tri, vol_info = read_surface(fname, read_metadata=True) + write_surface(op.join(tempdir, 'tmp'), pts, tri, volume_info=vol_info) + c_pts, c_tri, c_vol_info = read_surface(op.join(tempdir, 'tmp'), + read_metadata=True) assert_array_equal(pts, c_pts) assert_array_equal(tri, c_tri) + _is_equal_dict([vol_info, c_vol_info]) @testing.requires_testing_data From c5f6db14ef8ff2fb2fc3ea8af87678e83d1987de Mon Sep 17 00:00:00 2001 From: jaeilepp Date: Tue, 21 Jun 2016 10:42:45 +0200 Subject: [PATCH 06/18] Fixes. --- mne/bem.py | 35 ++++++++++++++++++----------------- mne/surface.py | 18 +++++++++++------- 2 files changed, 29 insertions(+), 24 deletions(-) diff --git a/mne/bem.py b/mne/bem.py index ba8760129fc..c4d03a8af28 100644 --- a/mne/bem.py +++ b/mne/bem.py @@ -1083,13 +1083,14 @@ def make_watershed_bem(subject, subjects_dir=None, overwrite=False, _symlink(surf_ws_out, surf_out) skip_symlink = False - if skip_symlink: - logger.info("Unable to create all symbolic links to .surf files in " - "bem folder. Use --overwrite option to recreate them.") - dest = op.join(bem_dir, 'watershed') - else: - logger.info("Symbolic links to .surf files created in bem folder") - dest = bem_dir + if skip_symlink: + logger.info("Unable to create all symbolic links to .surf files " + "in bem folder. Use --overwrite option to recreate " + "them.") + dest = op.join(bem_dir, 'watershed') + else: + logger.info("Symbolic links to .surf files created in bem folder") + dest = bem_dir logger.info("\nThank you for waiting.\nThe BEM triangulations for this " "subject are now available at:\n%s." % dest) @@ -1765,8 +1766,8 @@ 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') - surf_out = _load_ascii_surface(surf + '.tri', swap=True) - write_surface(surf + '.surf', surf_out[0], surf_out[1]) + nodes, tris = _load_ascii_surface(surf + '.tri', swap=True) + write_surface(surf + '.surf', nodes, tris) # Cleanup section logger.info("\n---- Cleaning up ----") @@ -1793,14 +1794,14 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, os.remove(surf) _symlink(op.join('flash', surf), op.join(surf)) skip_symlink = False - if skip_symlink: - logger.info("Unable to create all symbolic links to .surf files " - "in bem folder. Use --overwrite option to recreate " - "them.") - dest = op.join(bem_dir, 'flash') - else: - logger.info("Symbolic links to .surf files created in bem folder") - dest = bem_dir + if skip_symlink: + logger.info("Unable to create all symbolic links to .surf files " + "in bem folder. Use --overwrite option to recreate " + "them.") + dest = op.join(bem_dir, 'flash') + else: + logger.info("Symbolic links to .surf files created in bem folder") + dest = bem_dir logger.info("\nThank you for waiting.\nThe BEM triangulations for this " "subject are now available at:\n%s.\nWe hope the BEM meshes " "created will facilitate your MEG and EEG data analyses." diff --git a/mne/surface.py b/mne/surface.py index 4b724cfff33..9aa654a8253 100644 --- a/mne/surface.py +++ b/mne/surface.py @@ -13,7 +13,6 @@ import numpy as np from scipy.sparse import coo_matrix, csr_matrix, eye as speye -import nibabel as nib from .bem import read_bem_surfaces from .io.constants import FIFF @@ -415,8 +414,6 @@ def read_surface(fname, read_metadata=False, verbose=None): ---------- fname : str The name of the file containing the surface. - verbose : bool, str, int, or None - If not None, override default verbose level (see mne.verbose). read_metadata : bool Read metadata as key-value pairs. Valid keys: @@ -432,6 +429,9 @@ def read_surface(fname, read_metadata=False, verbose=None): .. versionadded:: 0.13.0 + verbose : bool, str, int, or None + If not None, override default verbose level (see mne.verbose). + Returns ------- rr : array, shape=(n_vertices, 3) @@ -445,10 +445,13 @@ def read_surface(fname, read_metadata=False, verbose=None): -------- write_surface """ - # XXX: Tests fail here due to numerical error. - # if LooseVersion(nib.__version__) > LooseVersion('2.1.0'): - # return nib.freesurfer.read_geometry(fname, - # read_metadata=read_metadata) + import nibabel as nib + if LooseVersion(nib.__version__) > LooseVersion('2.1.0'): + ret = nib.freesurfer.read_geometry(fname, read_metadata=read_metadata) + coords = ret[0].astype(np.float) # XXX: due to mayavi bug on mac 32b + if read_metadata: + return coords, ret[1], ret[2] + return coords, ret[1] volume_info = dict() TRIANGLE_MAGIC = 16777214 @@ -740,6 +743,7 @@ def write_surface(fname, coords, faces, create_stamp='', volume_info=None): -------- read_surface """ + import nibabel as nib if len(create_stamp.splitlines()) > 1: raise ValueError("create_stamp can only contain one line") From e5cfb1e48eaf1e13d5f0c745a9120129674fe7d8 Mon Sep 17 00:00:00 2001 From: jaeilepp Date: Tue, 21 Jun 2016 11:47:42 +0200 Subject: [PATCH 07/18] Volume info for tri surface. --- mne/bem.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/mne/bem.py b/mne/bem.py index c4d03a8af28..f1a9f3f3340 100644 --- a/mne/bem.py +++ b/mne/bem.py @@ -1767,7 +1767,11 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, for surf in surfs: shutil.move(op.join(bem_dir, surf + '.tri'), surf + '.tri') nodes, tris = _load_ascii_surface(surf + '.tri', swap=True) - write_surface(surf + '.surf', nodes, tris) + vol_info = _read_volume_info(op.join(subjects_dir, subject, 'mri', + 'flash', 'parameter_maps', + 'flash5_reg.mgz')) + vol_info['head'] = np.array([20]) + write_surface(surf + '.surf', nodes, tris, volume_info=vol_info) # Cleanup section logger.info("\n---- Cleaning up ----") @@ -1796,8 +1800,7 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, skip_symlink = False if skip_symlink: logger.info("Unable to create all symbolic links to .surf files " - "in bem folder. Use --overwrite option to recreate " - "them.") + "in bem folder. Use --overwrite option to recreate them.") dest = op.join(bem_dir, 'flash') else: logger.info("Symbolic links to .surf files created in bem folder") From 9a14a56860fb2728af60c18fcab6061403f79a68 Mon Sep 17 00:00:00 2001 From: jaeilepp Date: Tue, 21 Jun 2016 13:36:51 +0200 Subject: [PATCH 08/18] Checks for nibabel. --- mne/bem.py | 31 +++++++++++++++++++++---------- mne/surface.py | 28 ++++++++++++++++------------ 2 files changed, 37 insertions(+), 22 deletions(-) diff --git a/mne/bem.py b/mne/bem.py index f1a9f3f3340..a611e67bc9a 100644 --- a/mne/bem.py +++ b/mne/bem.py @@ -1063,7 +1063,11 @@ 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): - new_info = _read_volume_info(T1_mgz) + 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)) @@ -1116,10 +1120,13 @@ def make_watershed_bem(subject, subjects_dir=None, overwrite=False, logger.info('Created %s\n\nComplete.' % (fname_head,)) -def _read_volume_info(T1): - """Helper for extracting volume info from T1.mgz.""" - import nibabel as nib - header = nib.load(T1).header +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: @@ -1127,7 +1134,7 @@ def _read_volume_info(T1): else: raise ValueError('Volume info invalid.') new_info['valid'] = version - new_info['filename'] = T1 + 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 @@ -1767,10 +1774,14 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, for surf in surfs: shutil.move(op.join(bem_dir, surf + '.tri'), surf + '.tri') nodes, tris = _load_ascii_surface(surf + '.tri', swap=True) - vol_info = _read_volume_info(op.join(subjects_dir, subject, 'mri', - 'flash', 'parameter_maps', - 'flash5_reg.mgz')) - vol_info['head'] = np.array([20]) + vol_info = _extract_volume_info(op.join(subjects_dir, subject, 'mri', + 'flash', 'parameter_maps', + 'flash5_reg.mgz')) + 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 diff --git a/mne/surface.py b/mne/surface.py index 9aa654a8253..1bcacf7add2 100644 --- a/mne/surface.py +++ b/mne/surface.py @@ -445,13 +445,13 @@ def read_surface(fname, read_metadata=False, verbose=None): -------- write_surface """ - import nibabel as nib - if LooseVersion(nib.__version__) > LooseVersion('2.1.0'): - ret = nib.freesurfer.read_geometry(fname, read_metadata=read_metadata) - coords = ret[0].astype(np.float) # XXX: due to mayavi bug on mac 32b - if read_metadata: - return coords, ret[1], ret[2] - return coords, ret[1] + try: + import nibabel as nib + has_nibabel = True + except ImportError: + has_nibabel = False + if has_nibabel and LooseVersion(nib.__version__) > LooseVersion('2.1.0'): + return nib.freesurfer.read_geometry(fname, read_metadata=read_metadata) volume_info = dict() TRIANGLE_MAGIC = 16777214 @@ -743,15 +743,19 @@ def write_surface(fname, coords, faces, create_stamp='', volume_info=None): -------- read_surface """ - import nibabel as nib - if len(create_stamp.splitlines()) > 1: - raise ValueError("create_stamp can only contain one line") - - if LooseVersion(nib.__version__) > LooseVersion('2.1.0'): + try: + import nibabel as nib + has_nibabel = True + except ImportError: + has_nibabel = False + if has_nibabel and LooseVersion(nib.__version__) > LooseVersion('2.1.0'): nib.freesurfer.io.write_geometry(fname, coords, faces, create_stamp=create_stamp, volume_info=volume_info) return + if len(create_stamp.splitlines()) > 1: + raise ValueError("create_stamp can only contain one line") + with open(fname, 'wb') as fid: fid.write(pack('>3B', 255, 255, 254)) strs = ['%s\n' % create_stamp, '\n'] From 2757b9ebfc384586a263c881650015899be4e5bc Mon Sep 17 00:00:00 2001 From: jaeilepp Date: Tue, 21 Jun 2016 15:22:25 +0200 Subject: [PATCH 09/18] Added flash_path to make_flash_bem. --- mne/bem.py | 24 +++++++++++++++--------- mne/commands/tests/test_commands.py | 1 - 2 files changed, 15 insertions(+), 10 deletions(-) diff --git a/mne/bem.py b/mne/bem.py index a611e67bc9a..2c77eab24e2 100644 --- a/mne/bem.py +++ b/mne/bem.py @@ -1670,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 @@ -1683,6 +1683,9 @@ 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), the current folder is + used. verbose : bool, str, int, or None If not None, override default verbose level (see mne.verbose). @@ -1695,9 +1698,6 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, 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 current - folder. - See Also -------- convert_flash_mris @@ -1719,13 +1719,19 @@ 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'): + if flash_path is None: + flash5 = 'flash5.mgz' + flash5_reg = 'flash5_reg.mgz' + else: + 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") @@ -1733,7 +1739,7 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, 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')] + 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) @@ -1776,7 +1782,7 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, nodes, tris = _load_ascii_surface(surf + '.tri', swap=True) vol_info = _extract_volume_info(op.join(subjects_dir, subject, 'mri', 'flash', 'parameter_maps', - 'flash5_reg.mgz')) + flash5_reg)) if vol_info is None: warn('nibabel is required to update the volume info. Volume info ' 'omitted from the written surface.') diff --git a/mne/commands/tests/test_commands.py b/mne/commands/tests/test_commands.py index 180017f723b..4d86a49479c 100644 --- a/mne/commands/tests/test_commands.py +++ b/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(): From 17a070767a56aae4b9fdcab9cbc1ad293ec69c77 Mon Sep 17 00:00:00 2001 From: jaeilepp Date: Wed, 22 Jun 2016 13:22:10 +0200 Subject: [PATCH 10/18] Test for make_flash_bem. --- mne/bem.py | 37 ++++++++++++++++++++--------------- mne/commands/mne_flash_bem.py | 2 +- mne/datasets/utils.py | 4 ++-- mne/surface.py | 2 ++ mne/tests/test_bem.py | 23 +++++++++++++++++++++- 5 files changed, 48 insertions(+), 20 deletions(-) diff --git a/mne/bem.py b/mne/bem.py index 2c77eab24e2..d6242f388e8 100644 --- a/mne/bem.py +++ b/mne/bem.py @@ -1684,8 +1684,11 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, 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), the current folder is - used. + 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). @@ -1704,10 +1707,15 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, """ from .viz.misc import plot_bem from .surface import write_surface, _load_ascii_surface + + # Hack to enable testing with travis. + istest = os.environ.get('FREESURFER_HOME', '').endswith('MNE-testing-data') env, mri_dir, bem_dir = _prepare_env(subject, subjects_dir, requires_freesurfer=True, requires_mne=True) + if flash_path is None: + flash_path = op.join(mri_dir, 'flash', 'parameter_maps') curdir = os.getcwd() subjects_dir = env['SUBJECTS_DIR'] @@ -1719,12 +1727,8 @@ 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 flash_path is None: - flash5 = 'flash5.mgz' - flash5_reg = 'flash5_reg.mgz' - else: - flash5 = op.join(flash_path, 'flash5.mgz') - flash5_reg = op.join(flash_path, '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') @@ -1739,8 +1743,9 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, 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, op.join(mri_dir, 'flash5')] - run_subprocess(cmd, env=env, stdout=sys.stdout, stderr=sys.stderr) + if not istest: + 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 @@ -1768,9 +1773,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 istest: + 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'): @@ -1780,9 +1787,7 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, for surf in surfs: shutil.move(op.join(bem_dir, surf + '.tri'), surf + '.tri') nodes, tris = _load_ascii_surface(surf + '.tri', swap=True) - vol_info = _extract_volume_info(op.join(subjects_dir, subject, 'mri', - 'flash', 'parameter_maps', - flash5_reg)) + 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.') diff --git a/mne/commands/mne_flash_bem.py b/mne/commands/mne_flash_bem.py index bf65be4195d..cc52010fee7 100644 --- a/mne/commands/mne_flash_bem.py +++ b/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: diff --git a/mne/datasets/utils.py b/mne/datasets/utils.py index db47a77039e..27c9f41fd14 100644 --- a/mne/datasets/utils.py +++ b/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='24c7f0cd1ed94f359a3fa5478a4c42ce' ) folder_origs = dict( # not listed means None misc='mne-misc-data-%s' % releases['misc'], diff --git a/mne/surface.py b/mne/surface.py index 1bcacf7add2..c98ad79d02d 100644 --- a/mne/surface.py +++ b/mne/surface.py @@ -1181,6 +1181,7 @@ def mesh_dist(tris, vert): def _load_ascii_surface(filepath, swap=False): """Function for reading triangle definitions from an ascii file. + Parameters ---------- fname_in : str @@ -1188,6 +1189,7 @@ def _load_ascii_surface(filepath, swap=False): swap : bool Assume the ASCII file vertex ordering is clockwise instead of counterclockwise. + Returns ------- surf : tuple (nodes, tris) diff --git a/mne/tests/test_bem.py b/mne/tests/test_bem.py index 4c8b1e56fb1..5815b58cd42 100644 --- a/mne/tests/test_bem.py +++ b/mne/tests/test_bem.py @@ -20,9 +20,13 @@ from mne.utils import run_tests_if_main, _TempDir, slow_test, catch_logging from mne.bem import (_ico_downsample, _get_ico_map, _order_surfaces, _assert_complete_surface, _assert_inside, - _check_surface_size, _bem_find_surface) + _check_surface_size, _bem_find_surface, make_flash_bem) +from mne.surface import read_surface from mne.io import read_info +import matplotlib +matplotlib.use('Agg') # for testing don't use X server + warnings.simplefilter('always') fname_raw = op.join(op.dirname(__file__), '..', 'io', 'tests', 'data', @@ -333,4 +337,21 @@ def test_fit_sphere_to_headshape(): assert_raises(TypeError, fit_sphere_to_headshape, 1, units='m') +@testing.requires_testing_data +def test_make_flash_bem(): + """Test computing bem from flash images.""" + import matplotlib.pyplot as plt + flash_path = op.join(subjects_dir, 'sample', 'mri', 'flash') + # This function deletes some files at the end. + make_flash_bem('sample', overwrite=True, subjects_dir=subjects_dir, + flash_path=flash_path) + plt.close('all') + inner_skull = op.join(subjects_dir, 'sample', 'bem', 'inner_skull.surf') + outer_skull = op.join(subjects_dir, 'sample', 'bem', 'outer_skull.surf') + outer_skin = op.join(subjects_dir, 'sample', 'bem', 'outer_skin.surf') + for surf in (inner_skull, outer_skull, outer_skin): + coords, faces = read_surface(surf) + assert_equal(0, faces.min()) + assert_equal(coords.shape[0], faces.max() + 1) + run_tests_if_main() From a42e1c7368be92456ef97fde047971900a376b71 Mon Sep 17 00:00:00 2001 From: jaeilepp Date: Wed, 22 Jun 2016 14:07:06 +0200 Subject: [PATCH 11/18] Fix to testing. --- mne/bem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mne/bem.py b/mne/bem.py index d6242f388e8..08063662dca 100644 --- a/mne/bem.py +++ b/mne/bem.py @@ -1711,7 +1711,7 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, # Hack to enable testing with travis. istest = os.environ.get('FREESURFER_HOME', '').endswith('MNE-testing-data') env, mri_dir, bem_dir = _prepare_env(subject, subjects_dir, - requires_freesurfer=True, + requires_freesurfer=not istest, requires_mne=True) if flash_path is None: From 8931e10a52af4ba0ac01e648cfc89d9956d7885b Mon Sep 17 00:00:00 2001 From: jaeilepp Date: Thu, 23 Jun 2016 08:12:15 +0200 Subject: [PATCH 12/18] Fix. --- mne/bem.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/mne/bem.py b/mne/bem.py index 08063662dca..202213a8e3b 100644 --- a/mne/bem.py +++ b/mne/bem.py @@ -1709,9 +1709,10 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, from .surface import write_surface, _load_ascii_surface # Hack to enable testing with travis. - istest = os.environ.get('FREESURFER_HOME', '').endswith('MNE-testing-data') + is_test = not os.environ.get('FREESURFER_HOME', '').endswith( + 'MNE-testing-data') env, mri_dir, bem_dir = _prepare_env(subject, subjects_dir, - requires_freesurfer=not istest, + requires_freesurfer=not is_test, requires_mne=True) if flash_path is None: @@ -1743,7 +1744,7 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, 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')) - if not istest: + if not is_test: 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 @@ -1773,7 +1774,7 @@ 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 - if not istest: + if not is_test: logger.info("\n---- Creating the BEM surfaces ----") cmd = ['mri_make_bem_surfaces', subject] run_subprocess(cmd, env=env, stdout=sys.stdout, stderr=sys.stderr) From d891b9f89f88c0451b17531a683f2197d635f596 Mon Sep 17 00:00:00 2001 From: jaeilepp Date: Thu, 23 Jun 2016 08:59:21 +0200 Subject: [PATCH 13/18] Better fix. --- mne/bem.py | 6 +++--- mne/tests/test_bem.py | 4 +++- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/mne/bem.py b/mne/bem.py index 202213a8e3b..de9666849d9 100644 --- a/mne/bem.py +++ b/mne/bem.py @@ -1709,11 +1709,11 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, from .surface import write_surface, _load_ascii_surface # Hack to enable testing with travis. - is_test = not os.environ.get('FREESURFER_HOME', '').endswith( + is_test = os.environ.get('FREESURFER_HOME', '').endswith( 'MNE-testing-data') env, mri_dir, bem_dir = _prepare_env(subject, subjects_dir, - requires_freesurfer=not is_test, - requires_mne=True) + requires_freesurfer=True, + requires_mne=False) if flash_path is None: flash_path = op.join(mri_dir, 'flash', 'parameter_maps') diff --git a/mne/tests/test_bem.py b/mne/tests/test_bem.py index 5815b58cd42..93e4dcb7555 100644 --- a/mne/tests/test_bem.py +++ b/mne/tests/test_bem.py @@ -17,7 +17,8 @@ from mne.io.constants import FIFF from mne.transforms import translation from mne.datasets import testing -from mne.utils import run_tests_if_main, _TempDir, slow_test, catch_logging +from mne.utils import (run_tests_if_main, _TempDir, slow_test, catch_logging, + requires_freesurfer) from mne.bem import (_ico_downsample, _get_ico_map, _order_surfaces, _assert_complete_surface, _assert_inside, _check_surface_size, _bem_find_surface, make_flash_bem) @@ -337,6 +338,7 @@ def test_fit_sphere_to_headshape(): assert_raises(TypeError, fit_sphere_to_headshape, 1, units='m') +@requires_freesurfer @testing.requires_testing_data def test_make_flash_bem(): """Test computing bem from flash images.""" From 443160fe34058717aece37718782d6e73227c5f7 Mon Sep 17 00:00:00 2001 From: jaeilepp Date: Thu, 23 Jun 2016 10:11:06 +0200 Subject: [PATCH 14/18] Fix doc. --- mne/bem.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/mne/bem.py b/mne/bem.py index de9666849d9..8f0b803008c 100644 --- a/mne/bem.py +++ b/mne/bem.py @@ -1694,8 +1694,7 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, 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 From 7d65468e925df5da446284a23073859b91e21218 Mon Sep 17 00:00:00 2001 From: jaeilepp Date: Fri, 24 Jun 2016 08:56:44 +0200 Subject: [PATCH 15/18] Testing. --- mne/bem.py | 1 + mne/surface.py | 3 ++- mne/tests/test_bem.py | 36 ++++++++++++++++++++++++++---------- 3 files changed, 29 insertions(+), 11 deletions(-) diff --git a/mne/bem.py b/mne/bem.py index 8f0b803008c..df9cd74e09b 100644 --- a/mne/bem.py +++ b/mne/bem.py @@ -1786,6 +1786,7 @@ 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') + nodes, tris = _load_ascii_surface(surf + '.tri', swap=True) vol_info = _extract_volume_info(flash5_reg) if vol_info is None: diff --git a/mne/surface.py b/mne/surface.py index c98ad79d02d..efa88fc34a6 100644 --- a/mne/surface.py +++ b/mne/surface.py @@ -1193,7 +1193,8 @@ def _load_ascii_surface(filepath, swap=False): Returns ------- surf : tuple (nodes, tris) - The surface.""" + The surface. + """ with open(filepath, "r") as fid: lines = fid.readlines() n_nodes = int(lines[0]) diff --git a/mne/tests/test_bem.py b/mne/tests/test_bem.py index 93e4dcb7555..d230744bb0a 100644 --- a/mne/tests/test_bem.py +++ b/mne/tests/test_bem.py @@ -343,17 +343,33 @@ def test_fit_sphere_to_headshape(): def test_make_flash_bem(): """Test computing bem from flash images.""" import matplotlib.pyplot as plt + from shutil import copy + tmp = _TempDir() + bemdir = op.join(subjects_dir, 'sample', 'bem') flash_path = op.join(subjects_dir, 'sample', 'mri', 'flash') - # This function deletes some files at the end. - make_flash_bem('sample', overwrite=True, subjects_dir=subjects_dir, - flash_path=flash_path) + + for surf in ('inner_skull', 'outer_skull', 'outer_skin'): + copy(op.join(bemdir, surf + '.surf'), tmp) + copy(op.join(bemdir, surf + '.tri'), tmp) + copy(op.join(bemdir, 'inner_skull_tmp.tri'), tmp) + + # This function deletes the tri files at the end. + try: + make_flash_bem('sample', overwrite=True, subjects_dir=subjects_dir, + flash_path=flash_path) + for surf in ('inner_skull', 'outer_skull', 'outer_skin'): + coords, faces = read_surface(op.join(bemdir, surf + '.surf')) + coords_c, faces_c = read_surface(op.join(tmp, surf + '.surf')) + assert_equal(0, faces.min()) + assert_equal(coords.shape[0], faces.max() + 1) + assert_allclose(coords, coords_c) + assert_allclose(faces, faces_c) + finally: + for surf in ('inner_skull', 'outer_skull', 'outer_skin'): + copy(op.join(tmp, surf + '.tri'), bemdir) # return deleted tri + copy(op.join(tmp, surf + '.surf'), bemdir) # return moved surf + copy(op.join(tmp, 'inner_skull_tmp.tri'), bemdir) plt.close('all') - inner_skull = op.join(subjects_dir, 'sample', 'bem', 'inner_skull.surf') - outer_skull = op.join(subjects_dir, 'sample', 'bem', 'outer_skull.surf') - outer_skin = op.join(subjects_dir, 'sample', 'bem', 'outer_skin.surf') - for surf in (inner_skull, outer_skull, outer_skin): - coords, faces = read_surface(surf) - assert_equal(0, faces.min()) - assert_equal(coords.shape[0], faces.max() + 1) + run_tests_if_main() From c4f07ae303290109460c925a62483ee62e22e9d5 Mon Sep 17 00:00:00 2001 From: jaeilepp Date: Fri, 24 Jun 2016 16:01:48 +0200 Subject: [PATCH 16/18] Fixes. --- .travis.yml | 1 + appveyor.yml | 1 + mne/bem.py | 7 ++++--- mne/tests/test_bem.py | 6 ++++-- mne/utils.py | 1 + 5 files changed, 11 insertions(+), 5 deletions(-) diff --git a/.travis.yml b/.travis.yml index 9d57618236e..565dbf79ae4 100644 --- a/.travis.yml +++ b/.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_TESTING=1; fi; else export MNE_SKIP_TESTING_DATASET_TESTS=true; diff --git a/appveyor.yml b/appveyor.yml index 3953842fffc..09e71a5fea6 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -27,6 +27,7 @@ install: - "pip install nose-timer nibabel nitime" - "python setup.py develop" - "SET MNE_SKIP_NETWORK_TESTS=1" + - "SET MNE_TESTING=1" - "SET MNE_FORCE_SERIAL=true" # otherwise joblib will bomb - "SET MNE_LOGGING_LEVEL=warning" - "python -c \"import mne; mne.datasets.testing.data_path()\"" diff --git a/mne/bem.py b/mne/bem.py index df9cd74e09b..b1c06b48cc3 100644 --- a/mne/bem.py +++ b/mne/bem.py @@ -1707,15 +1707,16 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, from .viz.misc import plot_bem from .surface import write_surface, _load_ascii_surface - # Hack to enable testing with travis. - is_test = os.environ.get('FREESURFER_HOME', '').endswith( - 'MNE-testing-data') + is_test = os.environ.get('MNE_TESTING', False) # to enable travis testing + env, mri_dir, bem_dir = _prepare_env(subject, subjects_dir, requires_freesurfer=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'] diff --git a/mne/tests/test_bem.py b/mne/tests/test_bem.py index d230744bb0a..3edd7b96650 100644 --- a/mne/tests/test_bem.py +++ b/mne/tests/test_bem.py @@ -3,7 +3,9 @@ # License: BSD 3 clause from copy import deepcopy +from os import remove import os.path as op +from shutil import copy import warnings import numpy as np @@ -343,7 +345,6 @@ def test_fit_sphere_to_headshape(): def test_make_flash_bem(): """Test computing bem from flash images.""" import matplotlib.pyplot as plt - from shutil import copy tmp = _TempDir() bemdir = op.join(subjects_dir, 'sample', 'bem') flash_path = op.join(subjects_dir, 'sample', 'mri', 'flash') @@ -357,7 +358,7 @@ def test_make_flash_bem(): try: make_flash_bem('sample', overwrite=True, subjects_dir=subjects_dir, flash_path=flash_path) - for surf in ('inner_skull', 'outer_skull', 'outer_skin'): + for surf in ('inner_skull', 'outer_skull'): coords, faces = read_surface(op.join(bemdir, surf + '.surf')) coords_c, faces_c = read_surface(op.join(tmp, surf + '.surf')) assert_equal(0, faces.min()) @@ -366,6 +367,7 @@ def test_make_flash_bem(): assert_allclose(faces, faces_c) finally: for surf in ('inner_skull', 'outer_skull', 'outer_skin'): + remove(op.join(bemdir, surf + '.surf')) # delete symlinks copy(op.join(tmp, surf + '.tri'), bemdir) # return deleted tri copy(op.join(tmp, surf + '.surf'), bemdir) # return moved surf copy(op.join(tmp, 'inner_skull_tmp.tri'), bemdir) diff --git a/mne/utils.py b/mne/utils.py index 05c0d6885db..da2e0e5b993 100644 --- a/mne/utils.py +++ b/mne/utils.py @@ -1206,6 +1206,7 @@ def set_memmap_min_size(memmap_min_size): 'MNE_SKIP_TESTING_DATASET_TESTS', 'MNE_STIM_CHANNEL', 'MNE_USE_CUDA', + 'MNE_TESTING', 'SUBJECTS_DIR', ) From 1723797bd6c408920dce1e090fec038358c52a16 Mon Sep 17 00:00:00 2001 From: jaeilepp Date: Tue, 28 Jun 2016 09:08:11 +0200 Subject: [PATCH 17/18] Address comments. Updated dataset. --- .travis.yml | 2 +- appveyor.yml | 2 +- mne/bem.py | 6 +++--- mne/datasets/utils.py | 2 +- mne/surface.py | 19 +++++++++++-------- mne/tests/test_bem.py | 5 ++++- mne/tests/test_surface.py | 2 +- mne/utils.py | 2 +- 8 files changed, 23 insertions(+), 17 deletions(-) diff --git a/.travis.yml b/.travis.yml index 565dbf79ae4..89b0833f738 100644 --- a/.travis.yml +++ b/.travis.yml @@ -104,7 +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_TESTING=1; + export MNE_SKIP_FLASH_CALL=1; fi; else export MNE_SKIP_TESTING_DATASET_TESTS=true; diff --git a/appveyor.yml b/appveyor.yml index 09e71a5fea6..c14e4a5a735 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -27,7 +27,7 @@ install: - "pip install nose-timer nibabel nitime" - "python setup.py develop" - "SET MNE_SKIP_NETWORK_TESTS=1" - - "SET MNE_TESTING=1" + - "SET MNE_SKIP_FLASH_CALL=1" - "SET MNE_FORCE_SERIAL=true" # otherwise joblib will bomb - "SET MNE_LOGGING_LEVEL=warning" - "python -c \"import mne; mne.datasets.testing.data_path()\"" diff --git a/mne/bem.py b/mne/bem.py index b1c06b48cc3..281acbc1d7c 100644 --- a/mne/bem.py +++ b/mne/bem.py @@ -1705,9 +1705,9 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, convert_flash_mris """ from .viz.misc import plot_bem - from .surface import write_surface, _load_ascii_surface + from .surface import write_surface, read_tri - is_test = os.environ.get('MNE_TESTING', False) # to enable travis testing + is_test = os.environ.get('MNE_SKIP_FLASH_CALL', False) # to enable testing env, mri_dir, bem_dir = _prepare_env(subject, subjects_dir, requires_freesurfer=True, @@ -1788,7 +1788,7 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, for surf in surfs: shutil.move(op.join(bem_dir, surf + '.tri'), surf + '.tri') - nodes, tris = _load_ascii_surface(surf + '.tri', swap=True) + 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 ' diff --git a/mne/datasets/utils.py b/mne/datasets/utils.py index 27c9f41fd14..1ee00ff07ed 100644 --- a/mne/datasets/utils.py +++ b/mne/datasets/utils.py @@ -211,7 +211,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True, sample='1d5da3a809fded1ef5734444ab5bf857', somato='f3e3a8441477bb5bacae1d0c6e0964fb', spm='f61041e3f3f2ba0def8a2ca71592cc41', - testing='24c7f0cd1ed94f359a3fa5478a4c42ce' + testing='566d9e97421972c874fe8a76a837c6ed' ) folder_origs = dict( # not listed means None misc='mne-misc-data-%s' % releases['misc'], diff --git a/mne/surface.py b/mne/surface.py index efa88fc34a6..cc6fd5e1657 100644 --- a/mne/surface.py +++ b/mne/surface.py @@ -1179,7 +1179,7 @@ def mesh_dist(tris, vert): return dist_matrix -def _load_ascii_surface(filepath, swap=False): +def read_tri(fname_in, swap=False): """Function for reading triangle definitions from an ascii file. Parameters @@ -1192,10 +1192,13 @@ def _load_ascii_surface(filepath, swap=False): Returns ------- - surf : tuple (nodes, tris) - The surface. + rr : array, shape=(n_vertices, 3) + Coordinate points. + tris : int array, shape=(n_faces, 3) + Triangulation (each line contains indexes for three points which + together form a face). """ - with open(filepath, "r") as fid: + with open(fname_in, "r") as fid: lines = fid.readlines() n_nodes = int(lines[0]) n_tris = int(lines[n_nodes + 1]) @@ -1206,17 +1209,17 @@ def _load_ascii_surface(filepath, swap=False): inds = range(1, 4) else: raise IOError('Unrecognized format of data.') - nodes = np.array([np.array([float(v) for v in l.split()])[inds] - for l in lines[1:n_nodes + 1]]) + rr = np.array([np.array([float(v) for v in l.split()])[inds] + for l in lines[1:n_nodes + 1]]) tris = np.array([np.array([int(v) for v in l.split()])[inds] for l in lines[n_nodes + 2:n_nodes + 2 + n_tris]]) if swap: tris[:, [2, 1]] = tris[:, [1, 2]] tris -= 1 logger.info('Loaded surface from %s with %s nodes and %s triangles.' % - (filepath, n_nodes, n_tris)) + (fname_in, n_nodes, n_tris)) if n_items in [3, 4]: logger.info('Node normals were not included in the source file.') else: warn('Node normals were not read.') - return (nodes, tris) + return (rr, tris) diff --git a/mne/tests/test_bem.py b/mne/tests/test_bem.py index 3edd7b96650..3bf1e0dd7d8 100644 --- a/mne/tests/test_bem.py +++ b/mne/tests/test_bem.py @@ -353,13 +353,15 @@ def test_make_flash_bem(): copy(op.join(bemdir, surf + '.surf'), tmp) copy(op.join(bemdir, surf + '.tri'), tmp) copy(op.join(bemdir, 'inner_skull_tmp.tri'), tmp) + copy(op.join(bemdir, 'outer_skin_from_testing.surf'), tmp) # This function deletes the tri files at the end. try: make_flash_bem('sample', overwrite=True, subjects_dir=subjects_dir, flash_path=flash_path) - for surf in ('inner_skull', 'outer_skull'): + for surf in ('inner_skull', 'outer_skull', 'outer_skin'): coords, faces = read_surface(op.join(bemdir, surf + '.surf')) + surf = 'outer_skin_from_testing' if surf == 'outer_skin' else surf coords_c, faces_c = read_surface(op.join(tmp, surf + '.surf')) assert_equal(0, faces.min()) assert_equal(coords.shape[0], faces.max() + 1) @@ -371,6 +373,7 @@ def test_make_flash_bem(): copy(op.join(tmp, surf + '.tri'), bemdir) # return deleted tri copy(op.join(tmp, surf + '.surf'), bemdir) # return moved surf copy(op.join(tmp, 'inner_skull_tmp.tri'), bemdir) + copy(op.join(tmp, 'outer_skin_from_testing.surf'), bemdir) plt.close('all') diff --git a/mne/tests/test_surface.py b/mne/tests/test_surface.py index f8cba010be9..2ff52f4e265 100644 --- a/mne/tests/test_surface.py +++ b/mne/tests/test_surface.py @@ -137,7 +137,7 @@ def test_io_surface(): read_metadata=True) assert_array_equal(pts, c_pts) assert_array_equal(tri, c_tri) - _is_equal_dict([vol_info, c_vol_info]) + assert_true(_is_equal_dict([vol_info, c_vol_info])) @testing.requires_testing_data diff --git a/mne/utils.py b/mne/utils.py index da2e0e5b993..107e181dec1 100644 --- a/mne/utils.py +++ b/mne/utils.py @@ -1206,7 +1206,7 @@ def set_memmap_min_size(memmap_min_size): 'MNE_SKIP_TESTING_DATASET_TESTS', 'MNE_STIM_CHANNEL', 'MNE_USE_CUDA', - 'MNE_TESTING', + 'MNE_SKIP_FLASH_CALL', 'SUBJECTS_DIR', ) From 0196b28f9ce36070d9ffd3082dd19693037ff045 Mon Sep 17 00:00:00 2001 From: jaeilepp Date: Wed, 29 Jun 2016 08:12:35 +0200 Subject: [PATCH 18/18] Address comments. --- .travis.yml | 2 +- appveyor.yml | 1 - doc/python_reference.rst | 1 + mne/__init__.py | 2 +- mne/bem.py | 6 +++--- mne/surface.py | 22 ++++++++++++++++++---- mne/utils.py | 2 +- 7 files changed, 25 insertions(+), 11 deletions(-) diff --git a/.travis.yml b/.travis.yml index 89b0833f738..38aa6fc2676 100644 --- a/.travis.yml +++ b/.travis.yml @@ -104,7 +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_FLASH_CALL=1; + export MNE_SKIP_FS_FLASH_CALL=1; fi; else export MNE_SKIP_TESTING_DATASET_TESTS=true; diff --git a/appveyor.yml b/appveyor.yml index c14e4a5a735..3953842fffc 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -27,7 +27,6 @@ install: - "pip install nose-timer nibabel nitime" - "python setup.py develop" - "SET MNE_SKIP_NETWORK_TESTS=1" - - "SET MNE_SKIP_FLASH_CALL=1" - "SET MNE_FORCE_SERIAL=true" # otherwise joblib will bomb - "SET MNE_LOGGING_LEVEL=warning" - "python -c \"import mne; mne.datasets.testing.data_path()\"" diff --git a/doc/python_reference.rst b/doc/python_reference.rst index 15f97941537..32c33ac1043 100644 --- a/doc/python_reference.rst +++ b/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 diff --git a/mne/__init__.py b/mne/__init__.py index ba128b34703..e7d8d12ef5b 100644 --- a/mne/__init__.py +++ b/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, diff --git a/mne/bem.py b/mne/bem.py index 281acbc1d7c..b9cee310946 100644 --- a/mne/bem.py +++ b/mne/bem.py @@ -1707,7 +1707,7 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, from .viz.misc import plot_bem from .surface import write_surface, read_tri - is_test = os.environ.get('MNE_SKIP_FLASH_CALL', False) # to enable testing + is_test = os.environ.get('MNE_SKIP_FS_FLASH_CALL', False) env, mri_dir, bem_dir = _prepare_env(subject, subjects_dir, requires_freesurfer=True, @@ -1744,7 +1744,7 @@ def make_flash_bem(subject, overwrite=False, show=True, subjects_dir=None, 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')) - if not is_test: + 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 @@ -1774,7 +1774,7 @@ 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 - if not is_test: + 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) diff --git a/mne/surface.py b/mne/surface.py index cc6fd5e1657..82dc966995c 100644 --- a/mne/surface.py +++ b/mne/surface.py @@ -437,13 +437,14 @@ def read_surface(fname, read_metadata=False, verbose=None): rr : array, shape=(n_vertices, 3) Coordinate points. tris : int array, shape=(n_faces, 3) - Triangulation (each line contains indexes for three points which + Triangulation (each line contains indices for three points which together form a face). volume_info : dict-like If read_metadata is true, key-value pairs found in the geometry file. See Also -------- write_surface + read_tri """ try: import nibabel as nib @@ -719,7 +720,7 @@ def write_surface(fname, coords, faces, create_stamp='', volume_info=None): coords : array, shape=(n_vertices, 3) Coordinate points. faces : int array, shape=(n_faces, 3) - Triangulation (each line contains indexes for three points which + Triangulation (each line contains indices for three points which together form a face). create_stamp : str Comment that is written to the beginning of the file. Can not contain @@ -742,6 +743,7 @@ def write_surface(fname, coords, faces, create_stamp='', volume_info=None): See Also -------- read_surface + read_tri """ try: import nibabel as nib @@ -1179,7 +1181,8 @@ def mesh_dist(tris, vert): return dist_matrix -def read_tri(fname_in, swap=False): +@verbose +def read_tri(fname_in, swap=False, verbose=None): """Function for reading triangle definitions from an ascii file. Parameters @@ -1189,14 +1192,25 @@ def read_tri(fname_in, swap=False): swap : bool Assume the ASCII file vertex ordering is clockwise instead of counterclockwise. + verbose : bool, str, int, or None + If not None, override default verbose level (see mne.verbose). Returns ------- rr : array, shape=(n_vertices, 3) Coordinate points. tris : int array, shape=(n_faces, 3) - Triangulation (each line contains indexes for three points which + Triangulation (each line contains indices for three points which together form a face). + + Notes + ----- + .. versionadded:: 0.13.0 + + See Also + -------- + read_surface + write_surface """ with open(fname_in, "r") as fid: lines = fid.readlines() diff --git a/mne/utils.py b/mne/utils.py index 107e181dec1..2dececc1425 100644 --- a/mne/utils.py +++ b/mne/utils.py @@ -1206,7 +1206,7 @@ def set_memmap_min_size(memmap_min_size): 'MNE_SKIP_TESTING_DATASET_TESTS', 'MNE_STIM_CHANNEL', 'MNE_USE_CUDA', - 'MNE_SKIP_FLASH_CALL', + 'MNE_SKIP_FS_FLASH_CALL', 'SUBJECTS_DIR', )