From ddf9db1a929cda46d925db01bd5df8e35ab5b3fc Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Mon, 6 Jun 2016 09:54:41 -0400 Subject: [PATCH 1/6] NF: Handle metadata in read/write_geometry --- nibabel/freesurfer/io.py | 63 +++++++++++++++++++++++++++++++++++++--- 1 file changed, 59 insertions(+), 4 deletions(-) diff --git a/nibabel/freesurfer/io.py b/nibabel/freesurfer/io.py index 2ac5592090..08568c6c88 100644 --- a/nibabel/freesurfer/io.py +++ b/nibabel/freesurfer/io.py @@ -3,6 +3,7 @@ import numpy as np import getpass import time +from collections import OrderedDict from ..externals.six.moves import xrange @@ -44,13 +45,17 @@ def _fread3_many(fobj, n): return (b1 << 16) + (b2 << 8) + b3 -def read_geometry(filepath): +def read_geometry(filepath, read_metadata=False, read_stamp=False): """Read a triangular format Freesurfer surface mesh. Parameters ---------- filepath : str Path to surface file + read_metadata : bool + Read metadata as key-value pairs + read_stamp : bool + Return the comment from the file Returns ------- @@ -58,7 +63,14 @@ def read_geometry(filepath): nvtx x 3 array of vertex (x, y, z) coordinates faces : numpy array nfaces x 3 array of defining mesh triangles + volume_info : dict-like + If read_metadata is true, key-value pairs found in the geometry file + create_stamp : str + If read_stamp is true, the comment added by the program that saved + the file """ + volume_info = None + with open(filepath, "rb") as fobj: magic = _fread3(fobj) if magic == 16777215: # Quad file @@ -86,20 +98,46 @@ def read_geometry(filepath): nface += 1 elif magic == 16777214: # Triangle file - fobj.readline() # create_stamp + create_stamp = fobj.readline().rstrip(b'\n').decode('utf-8') fobj.readline() vnum = np.fromfile(fobj, ">i4", 1)[0] 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) + + extra = fobj.read() if read_metadata else b'' + if extra: + if extra[:4] != b'\x00\x00\x00\x14': + warnings.warn("Unknown extension code.") + volume_info = OrderedDict() + try: + for line in extra[4:].split(b'\n'): + key, val = map(bytes.strip, line.split(b'=', 1)) + key = key.decode('utf-8') + if key in ('voxelsize', 'xras', 'yras', 'zras'): + val = np.fromstring(val, sep=' ') + elif key == 'volume': + val = np.fromstring(val, sep=' ', dtype=np.uint) + volume_info[key] = val + except ValueError: + raise ValueError("Error parsing volume info") + else: raise ValueError("File does not appear to be a Freesurfer surface") coords = coords.astype(np.float) # XXX: due to mayavi bug on mac 32bits - return coords, faces + ret = (coords, faces) + if read_metadata: + ret += (volume_info,) + if read_stamp: + ret += (create_stamp,) -def write_geometry(filepath, coords, faces, create_stamp=None): + return ret + + +def write_geometry(filepath, coords, faces, create_stamp=None, + volume_info=None): """Write a triangular format Freesurfer surface mesh. Parameters @@ -112,6 +150,8 @@ def write_geometry(filepath, coords, faces, create_stamp=None): nfaces x 3 array of defining mesh triangles create_stamp : str User/time stamp (default: "created by on ") + volume_info : dict-like or None + Key-value pairs to encode at the end of the file """ magic_bytes = np.array([255, 255, 254], dtype=np.uint8) @@ -119,6 +159,18 @@ def write_geometry(filepath, coords, faces, create_stamp=None): create_stamp = "created by %s on %s" % (getpass.getuser(), time.ctime()) + postlude = b'' + if volume_info is not None: + postlude = [b'\x00\x00\x00\x14'] + for key, val in volume_info.items(): + if key in ('voxelsize', 'xras', 'yras', 'zras'): + val = '{:.3f} {:.3f} {:.3f}'.format(*val) + elif key == 'volume': + val = '{:d} {:d} {:d}'.format(*val) + key = key.ljust(6) + postlude.append('{} = {}'.format(key, val).encode('utf-8')) + postlude = b'\n'.join(postlude) + with open(filepath, 'wb') as fobj: magic_bytes.tofile(fobj) fobj.write(("%s\n\n" % create_stamp).encode('utf-8')) @@ -129,6 +181,9 @@ def write_geometry(filepath, coords, faces, create_stamp=None): coords.astype('>f4').reshape(-1).tofile(fobj) faces.astype('>i4').reshape(-1).tofile(fobj) + # Add volume info, if given + fobj.write(postlude) + def read_morph_data(filepath): """Read a Freesurfer morphometry data file. From 63fd6c979cf27f51ad36a125973b60d8e40c7826 Mon Sep 17 00:00:00 2001 From: Alexandre Gramfort Date: Tue, 7 Jun 2016 14:42:50 +0200 Subject: [PATCH 2/6] add test + misc --- nibabel/freesurfer/io.py | 40 +++++++++++++++++++---------- nibabel/freesurfer/tests/test_io.py | 25 +++++++++++++----- 2 files changed, 44 insertions(+), 21 deletions(-) diff --git a/nibabel/freesurfer/io.py b/nibabel/freesurfer/io.py index 08568c6c88..7ab00d52ba 100644 --- a/nibabel/freesurfer/io.py +++ b/nibabel/freesurfer/io.py @@ -1,5 +1,6 @@ from __future__ import division, print_function, absolute_import +import warnings import numpy as np import getpass import time @@ -107,20 +108,31 @@ def read_geometry(filepath, read_metadata=False, read_stamp=False): extra = fobj.read() if read_metadata else b'' if extra: + volume_info = OrderedDict() + if extra[:4] != b'\x00\x00\x00\x14': warnings.warn("Unknown extension code.") - volume_info = OrderedDict() - try: - for line in extra[4:].split(b'\n'): - key, val = map(bytes.strip, line.split(b'=', 1)) - key = key.decode('utf-8') - if key in ('voxelsize', 'xras', 'yras', 'zras'): - val = np.fromstring(val, sep=' ') - elif key == 'volume': - val = np.fromstring(val, sep=' ', dtype=np.uint) - volume_info[key] = val - except ValueError: - raise ValueError("Error parsing volume info") + else: + try: + for line in extra[4:].split(b'\n'): + if len(line) == 0: + continue + key, val = map(bytes.strip, line.split(b'=', 1)) + print(key, val) + key = key.decode('utf-8') + if key in ('voxelsize', 'xras', 'yras', 'zras', 'cras'): + val = np.fromstring(val, sep=' ') + val = val.astype(np.float) + elif key == 'volume': + val = np.fromstring(val, sep=' ', dtype=np.uint) + val = val.astype(np.int) + volume_info[key] = val + except ValueError: + raise ValueError("Error parsing volume info") + + if len(volume_info) == 0: + warnings.warn("Volume geometry info is either " + "not contained or not valid.") else: raise ValueError("File does not appear to be a Freesurfer surface") @@ -160,10 +172,10 @@ def write_geometry(filepath, coords, faces, create_stamp=None, time.ctime()) postlude = b'' - if volume_info is not None: + if volume_info is not None and len(volume_info) > 0: postlude = [b'\x00\x00\x00\x14'] for key, val in volume_info.items(): - if key in ('voxelsize', 'xras', 'yras', 'zras'): + if key in ('voxelsize', 'xras', 'yras', 'zras', 'cras'): val = '{:.3f} {:.3f} {:.3f}'.format(*val) elif key == 'volume': val = '{:d} {:d} {:d}'.format(*val) diff --git a/nibabel/freesurfer/tests/test_io.py b/nibabel/freesurfer/tests/test_io.py index 69728f57a4..ce26518e52 100644 --- a/nibabel/freesurfer/tests/test_io.py +++ b/nibabel/freesurfer/tests/test_io.py @@ -4,6 +4,7 @@ import getpass import time import hashlib +import warnings from ...tmpdirs import InTemporaryDirectory @@ -17,7 +18,7 @@ from ...tests.nibabel_data import get_nibabel_data from ...fileslice import strided_scalar - +from ...testing import clear_and_catch_warnings DATA_SDIR = 'fsaverage' @@ -55,10 +56,17 @@ def test_geometry(): assert_equal(coords.shape[0], faces.max() + 1) # Test quad with sphere - surf_path = pjoin(data_path, "surf", "%s.%s" % ("lh", "sphere")) - coords, faces = read_geometry(surf_path) - assert_equal(0, faces.min()) - assert_equal(coords.shape[0], faces.max() + 1) + with clear_and_catch_warnings() as w: + warnings.filterwarnings('always', category=DeprecationWarning) + surf_path = pjoin(data_path, "surf", "%s.%s" % ("lh", "sphere")) + coords, faces, volume_info, create_stamp = \ + read_geometry(surf_path, read_metadata=True, read_stamp=True) + assert_equal(0, faces.min()) + assert_equal(coords.shape[0], faces.max() + 1) + assert_equal(0, len(volume_info)) + assert_equal(u'created by greve on Thu Jun 8 19:17:51 2006', + create_stamp) + assert_equal(len(w), 2) # Test equivalence of freesurfer- and nibabel-generated triangular files # with respect to read_geometry() @@ -66,10 +74,13 @@ def test_geometry(): surf_path = 'test' create_stamp = "created by %s on %s" % (getpass.getuser(), time.ctime()) - write_geometry(surf_path, coords, faces, create_stamp) + volume_info['cras'] = np.array([1., 2., 3.]) + write_geometry(surf_path, coords, faces, create_stamp, volume_info) - coords2, faces2 = read_geometry(surf_path) + coords2, faces2, volume_info2 = \ + read_geometry(surf_path, read_metadata=True) + assert_equal(volume_info2['cras'], volume_info['cras']) with open(surf_path, 'rb') as fobj: np.fromfile(fobj, ">u1", 3) read_create_stamp = fobj.readline().decode().rstrip('\n') From 995dc758620ed9e47348a934c779e9726213f019 Mon Sep 17 00:00:00 2001 From: Alexandre Gramfort Date: Tue, 7 Jun 2016 22:46:59 +0200 Subject: [PATCH 3/6] address comments --- nibabel/freesurfer/io.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/nibabel/freesurfer/io.py b/nibabel/freesurfer/io.py index 7ab00d52ba..3583545109 100644 --- a/nibabel/freesurfer/io.py +++ b/nibabel/freesurfer/io.py @@ -4,9 +4,8 @@ import numpy as np import getpass import time -from collections import OrderedDict - +from ..externals import OrderedDict from ..externals.six.moves import xrange from ..openers import Opener @@ -176,7 +175,7 @@ def write_geometry(filepath, coords, faces, create_stamp=None, postlude = [b'\x00\x00\x00\x14'] for key, val in volume_info.items(): if key in ('voxelsize', 'xras', 'yras', 'zras', 'cras'): - val = '{:.3f} {:.3f} {:.3f}'.format(*val) + val = '{:.4f} {:.4f} {:.4f}'.format(*val) elif key == 'volume': val = '{:d} {:d} {:d}'.format(*val) key = key.ljust(6) From c41faadc5d8f73b9484c29476364795d9c87346e Mon Sep 17 00:00:00 2001 From: Alexandre Gramfort Date: Wed, 8 Jun 2016 08:43:48 +0200 Subject: [PATCH 4/6] address comments --- nibabel/freesurfer/io.py | 68 ++++++++++++++++++++++------------------ 1 file changed, 37 insertions(+), 31 deletions(-) diff --git a/nibabel/freesurfer/io.py b/nibabel/freesurfer/io.py index 3583545109..f0c76da9e3 100644 --- a/nibabel/freesurfer/io.py +++ b/nibabel/freesurfer/io.py @@ -45,6 +45,38 @@ def _fread3_many(fobj, n): return (b1 << 16) + (b2 << 8) + b3 +def _read_volume_info(extra): + volume_info = OrderedDict() + + if extra is None: + return volume_info + + if extra[:4] != b'\x00\x00\x00\x14': + warnings.warn("Unknown extension code.") + else: + try: + for line in extra[4:].split(b'\n'): + if len(line) == 0: + continue + key, val = map(bytes.strip, line.split(b'=', 1)) + key = key.decode('utf-8') + if key in ('voxelsize', 'xras', 'yras', 'zras', 'cras'): + val = np.fromstring(val, sep=' ') + val = val.astype(np.float) + elif key == 'volume': + val = np.fromstring(val, sep=' ', dtype=np.uint) + val = val.astype(np.int) + volume_info[key] = val + except ValueError: + raise ValueError("Error parsing volume info") + + if len(volume_info) == 0: + warnings.warn("Volume geometry info is either " + "not contained or not valid.") + + return volume_info + + def read_geometry(filepath, read_metadata=False, read_stamp=False): """Read a triangular format Freesurfer surface mesh. @@ -63,7 +95,7 @@ def read_geometry(filepath, read_metadata=False, read_stamp=False): nvtx x 3 array of vertex (x, y, z) coordinates faces : numpy array nfaces x 3 array of defining mesh triangles - volume_info : dict-like + volume_info : OrderedDict If read_metadata is true, key-value pairs found in the geometry file create_stamp : str If read_stamp is true, the comment added by the program that saved @@ -106,33 +138,7 @@ def read_geometry(filepath, read_metadata=False, read_stamp=False): faces = np.fromfile(fobj, ">i4", fnum * 3).reshape(fnum, 3) extra = fobj.read() if read_metadata else b'' - if extra: - volume_info = OrderedDict() - - if extra[:4] != b'\x00\x00\x00\x14': - warnings.warn("Unknown extension code.") - else: - try: - for line in extra[4:].split(b'\n'): - if len(line) == 0: - continue - key, val = map(bytes.strip, line.split(b'=', 1)) - print(key, val) - key = key.decode('utf-8') - if key in ('voxelsize', 'xras', 'yras', 'zras', 'cras'): - val = np.fromstring(val, sep=' ') - val = val.astype(np.float) - elif key == 'volume': - val = np.fromstring(val, sep=' ', dtype=np.uint) - val = val.astype(np.int) - volume_info[key] = val - except ValueError: - raise ValueError("Error parsing volume info") - - if len(volume_info) == 0: - warnings.warn("Volume geometry info is either " - "not contained or not valid.") - + volume_info = _read_volume_info(extra) else: raise ValueError("File does not appear to be a Freesurfer surface") @@ -175,11 +181,11 @@ def write_geometry(filepath, coords, faces, create_stamp=None, postlude = [b'\x00\x00\x00\x14'] for key, val in volume_info.items(): if key in ('voxelsize', 'xras', 'yras', 'zras', 'cras'): - val = '{:.4f} {:.4f} {:.4f}'.format(*val) + val = '{0:.4f} {1:.4f} {2:.4f}'.format(*val) elif key == 'volume': - val = '{:d} {:d} {:d}'.format(*val) + val = '{0:d} {1:d} {2:d}'.format(*val) key = key.ljust(6) - postlude.append('{} = {}'.format(key, val).encode('utf-8')) + postlude.append('{0} = {1}'.format(key, val).encode('utf-8')) postlude = b'\n'.join(postlude) with open(filepath, 'wb') as fobj: From 89203bbb6526c794063c7ff24792596436768f83 Mon Sep 17 00:00:00 2001 From: Alexandre Gramfort Date: Wed, 8 Jun 2016 09:41:58 +0200 Subject: [PATCH 5/6] don't count warnings --- nibabel/freesurfer/tests/test_io.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/nibabel/freesurfer/tests/test_io.py b/nibabel/freesurfer/tests/test_io.py index ce26518e52..f75ad88112 100644 --- a/nibabel/freesurfer/tests/test_io.py +++ b/nibabel/freesurfer/tests/test_io.py @@ -56,7 +56,7 @@ def test_geometry(): assert_equal(coords.shape[0], faces.max() + 1) # Test quad with sphere - with clear_and_catch_warnings() as w: + with clear_and_catch_warnings(): warnings.filterwarnings('always', category=DeprecationWarning) surf_path = pjoin(data_path, "surf", "%s.%s" % ("lh", "sphere")) coords, faces, volume_info, create_stamp = \ @@ -66,7 +66,6 @@ def test_geometry(): assert_equal(0, len(volume_info)) assert_equal(u'created by greve on Thu Jun 8 19:17:51 2006', create_stamp) - assert_equal(len(w), 2) # Test equivalence of freesurfer- and nibabel-generated triangular files # with respect to read_geometry() From 45055cd385adc7cb5fb5d2f29771f4bf25e5bc4f Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Fri, 10 Jun 2016 09:18:27 -0400 Subject: [PATCH 6/6] FIX: Minor fixes --- nibabel/freesurfer/io.py | 26 ++++++++++++-------------- nibabel/freesurfer/tests/test_io.py | 23 ++++++++++++++++------- 2 files changed, 28 insertions(+), 21 deletions(-) diff --git a/nibabel/freesurfer/io.py b/nibabel/freesurfer/io.py index f0c76da9e3..34c79719bf 100644 --- a/nibabel/freesurfer/io.py +++ b/nibabel/freesurfer/io.py @@ -48,11 +48,12 @@ def _fread3_many(fobj, n): def _read_volume_info(extra): volume_info = OrderedDict() - if extra is None: - return volume_info - - if extra[:4] != b'\x00\x00\x00\x14': - warnings.warn("Unknown extension code.") + if extra is None: # we pass b'' when user does not request info + pass + elif len(extra) == 0: + warnings.warn('No volume information contained in the file') + elif extra[:4] != b'\x00\x00\x00\x14': + warnings.warn("Unknown extension code for volume info: %r" % extra[:4]) else: try: for line in extra[4:].split(b'\n'): @@ -67,13 +68,11 @@ def _read_volume_info(extra): val = np.fromstring(val, sep=' ', dtype=np.uint) val = val.astype(np.int) volume_info[key] = val - except ValueError: - raise ValueError("Error parsing volume info") - - if len(volume_info) == 0: - warnings.warn("Volume geometry info is either " - "not contained or not valid.") - + except Exception as exp: + raise ValueError("Error parsing volume info: %s" % str(exp)) + if len(volume_info) == 0: + warnings.warn("Volume geometry info is either " + "not contained or not valid.") return volume_info @@ -136,8 +135,7 @@ def read_geometry(filepath, read_metadata=False, read_stamp=False): 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) - - extra = fobj.read() if read_metadata else b'' + extra = fobj.read() if read_metadata else None volume_info = _read_volume_info(extra) else: raise ValueError("File does not appear to be a Freesurfer surface") diff --git a/nibabel/freesurfer/tests/test_io.py b/nibabel/freesurfer/tests/test_io.py index f75ad88112..57a1bf5934 100644 --- a/nibabel/freesurfer/tests/test_io.py +++ b/nibabel/freesurfer/tests/test_io.py @@ -56,16 +56,17 @@ def test_geometry(): assert_equal(coords.shape[0], faces.max() + 1) # Test quad with sphere - with clear_and_catch_warnings(): + surf_path = pjoin(data_path, "surf", "%s.%s" % ("lh", "sphere")) + with clear_and_catch_warnings() as w: warnings.filterwarnings('always', category=DeprecationWarning) - surf_path = pjoin(data_path, "surf", "%s.%s" % ("lh", "sphere")) coords, faces, volume_info, create_stamp = \ read_geometry(surf_path, read_metadata=True, read_stamp=True) - assert_equal(0, faces.min()) - assert_equal(coords.shape[0], faces.max() + 1) - assert_equal(0, len(volume_info)) - assert_equal(u'created by greve on Thu Jun 8 19:17:51 2006', - create_stamp) + assert_true(any('extension code' in str(ww.message) for ww in w)) + assert_equal(0, faces.min()) + assert_equal(coords.shape[0], faces.max() + 1) + assert_equal(0, len(volume_info)) + assert_equal(u'created by greve on Thu Jun 8 19:17:51 2006', + create_stamp) # Test equivalence of freesurfer- and nibabel-generated triangular files # with respect to read_geometry() @@ -84,6 +85,14 @@ def test_geometry(): np.fromfile(fobj, ">u1", 3) read_create_stamp = fobj.readline().decode().rstrip('\n') + # now write an incomplete file + write_geometry(surf_path, coords, faces) + with clear_and_catch_warnings() as w: + warnings.filterwarnings('always', category=DeprecationWarning) + read_geometry(surf_path, read_metadata=True) + assert_true(any('volume information contained' in str(ww.message) + for ww in w)) + assert_equal(create_stamp, read_create_stamp) np.testing.assert_array_equal(coords, coords2)