fix: resolved conflict

Author: satra

bin/nib-ls

@@ -13,7 +13,7 @@ Output a summary table for neuroimaging files (resolution, dimensionality, etc.)
 from __future__ import division, print_function, absolute_import
 
 __author__ = 'Yaroslav Halchenko'
-__copyright__ = 'Copyright (c) 2011-2012 Yaroslav Halchenko ' \
+__copyright__ = 'Copyright (c) 2011-2015 Yaroslav Halchenko ' \
                 'and NiBabel contributors'
 __license__ = 'MIT'
 
@@ -141,6 +141,10 @@ def get_opt_parser():
                dest="verbose", default=0,
                help="Make more noise.  Could be specified multiple times"),
 
+        Option("-H", "--header-fields",
+               dest="header_fields", default='',
+               help="Header fields (comma separated) to be printed as well (if present)"),
+
         Option("-s", "--stats",
                action="store_true", dest='stats', default=False,
                help="Output basic data statistics"),
@@ -180,6 +184,23 @@ def proc_file(f, opts):
     else:
         row += [ '' ]
 
+    if opts.header_fields:
+        # signals "all fields"
+        if opts.header_fields == 'all':
+            # TODO: might vary across file types, thus prior sensing
+            # would be needed
+            header_fields = h.keys()
+        else:
+            header_fields = opts.header_fields.split(',')
+
+        for f in header_fields:
+            if not f: # skip empty
+                continue
+            try:
+                row += [str(h[f])]
+            except (KeyError, ValueError):
+                row += [ 'error' ]
+
     try:
         if (hasattr(h, 'get_qform') and hasattr(h, 'get_sform')
             and (h.get_qform() != h.get_sform()).any()):

nibabel/affines.py

@@ -8,6 +8,13 @@
 from .externals.six.moves import reduce
 
 
+class AffineError(ValueError):
+    """ Errors in calculating or using affines """
+    # Inherits from ValueError to keep compatibility with ValueError previously
+    # raised in append_diag
+    pass
+
+
 def apply_affine(aff, pts):
     """ Apply affine matrix `aff` to points `pts`
 
@@ -213,7 +220,7 @@ def append_diag(aff, steps, starts=()):
     if len(starts) == 0:
         starts = np.zeros(n_steps, dtype=steps.dtype)
     elif len(starts) != n_steps:
-        raise ValueError('Steps should have same length as starts')
+        raise AffineError('Steps should have same length as starts')
     old_n_out, old_n_in = aff.shape[0] - 1, aff.shape[1] - 1
     # make new affine
     aff_plus = np.zeros((old_n_out + n_steps + 1,
@@ -247,3 +254,45 @@ def dot_reduce(*args):
         ...  arg[N-2].dot(arg[N-1])))...``
     """
     return reduce(lambda x, y: np.dot(y, x), args[::-1])
+
+
+def voxel_sizes(affine):
+    r""" Return voxel size for each input axis given `affine`
+
+    The `affine` is the mapping between array (voxel) coordinates and mm
+    (world) coordinates.
+
+    The voxel size for the first voxel (array) axis is the distance moved in
+    world coordinates when moving one unit along the first voxel (array) axis.
+    This is the distance between the world coordinate of voxel (0, 0, 0) and
+    the world coordinate of voxel (1, 0, 0).  The world coordinate vector of
+    voxel coordinate vector (0, 0, 0) is given by ``v0 = affine.dot((0, 0, 0,
+    1)[:3]``.  The world coordinate vector of voxel vector (1, 0, 0) is
+    ``v1_ax1 = affine.dot((1, 0, 0, 1))[:3]``.  The final 1 in the voxel
+    vectors and the ``[:3]`` at the end are because the affine works on
+    homogenous coodinates.  The translations part of the affine is ``trans =
+    affine[:3, 3]``, and the rotations, zooms and shearing part of the affine
+    is ``rzs = affine[:3, :3]``. Because of the final 1 in the input voxel
+    vector, ``v0 == rzs.dot((0, 0, 0)) + trans``, and ``v1_ax1 == rzs.dot((1,
+    0, 0)) + trans``, and the difference vector is ``rzs.dot((0, 0, 0)) -
+    rzs.dot((1, 0, 0)) == rzs.dot((1, 0, 0)) == rzs[:, 0]``.  The distance
+    vectors in world coordinates between (0, 0, 0) and (1, 0, 0), (0, 1, 0),
+    (0, 0, 1) are given by ``rzs.dot(np.eye(3)) = rzs``.  The voxel sizes are
+    the Euclidean lengths of the distance vectors.  So, the voxel sizes are
+    the Euclidean lengths of the columns of the affine (excluding the last row
+    and column of the affine).
+
+    Parameters
+    ----------
+    affine : 2D array-like
+        Affine transformation array.  Usually shape (4, 4), but can be any 2D
+        array.
+
+    Returns
+    -------
+    vox_sizes : 1D array
+        Voxel sizes for each input axis of affine.  Usually 1D array length 3,
+        but in general has length (N-1) where input `affine` is shape (M, N).
+    """
+    top_left = affine[:-1, :-1]
+    return np.sqrt(np.sum(top_left ** 2, axis=0))

nibabel/freesurfer/io.py

@@ -6,6 +6,7 @@
 
 
 from .. externals.six.moves import xrange
+from ..openers import Opener
 
 
 def _fread3(fobj):
@@ -160,23 +161,50 @@ def read_morph_data(filepath):
     return curv
 
 
-def write_morph_data(filename, values):
-    '''
-    '''
-    with open(filename, 'wb') as f:
+def write_morph_data(file_like, values, fnum=0):
+    """Write Freesurfer morphometry data `values` to file-like `file_like`
 
-      # magic number
-      np.array([255], dtype='>u1').tofile(f)
-      np.array([255], dtype='>u1').tofile(f)
-      np.array([255], dtype='>u1').tofile(f)
+    Equivalent to FreeSurfer's `write_curv.m`_
 
-      # vertices number and two un-used int4
-      np.array([len(values)], dtype='>i4').tofile(f)
-      np.array([0], dtype='>i4').tofile(f)
-      np.array([1], dtype='>i4').tofile(f)
+    See also:
+    http://www.grahamwideman.com/gw/brain/fs/surfacefileformats.htm#CurvNew
 
-      # now the data
-      np.array(values, dtype='>f4').tofile(f)
+    .. _write_curv.m: \
+    https://github.com/neurodebian/freesurfer/blob/debian-sloppy/matlab/write_curv.m
+
+    Parameters
+    ----------
+    file_like : file-like
+        String containing path of file to be written, or file-like object, open
+        in binary write (`'wb'` mode, implementing the `write` method)
+    values : array-like
+        Surface morphometry values
+
+        Shape must be (N,), (N, 1), (1, N) or (N, 1, 1)
+    fnum : int, optional
+        Number of faces in the associated surface
+    """
+    magic_bytes = np.array([255, 255, 255], dtype=np.uint8)
+
+    vector = np.asarray(values)
+    vnum = np.prod(vector.shape)
+    if vector.shape not in ((vnum,), (vnum, 1), (1, vnum), (vnum, 1, 1)):
+        raise ValueError("Invalid shape: argument values must be a vector")
+
+    i4info = np.iinfo('i4')
+    if vnum > i4info.max:
+        raise ValueError("Too many values for morphometry file")
+    if not i4info.min <= fnum <= i4info.max:
+        raise ValueError("Argument fnum must be between {0} and {1}".format(
+                         i4info.min, i4info.max))
+
+    with Opener(file_like, 'wb') as fobj:
+        fobj.write(magic_bytes)
+
+        # vertex count, face count (unused), vals per vertex (only 1 supported)
+        fobj.write(np.array([vnum, fnum, 1], dtype='>i4'))
+
+        fobj.write(vector.astype('>f4'))
 
 
 def read_annot(filepath, orig_ids=False):

nibabel/freesurfer/tests/test_io.py

@@ -10,12 +10,13 @@
 
 from nose.tools import assert_true
 import numpy as np
-from numpy.testing import assert_equal, dec
+from numpy.testing import assert_equal, assert_raises, dec
 
 from .. import (read_geometry, read_morph_data, read_annot, read_label,
-                write_geometry, write_annot)
+                write_geometry, write_morph_data, write_annot)
 
 from ...tests.nibabel_data import get_nibabel_data
+from ...fileslice import strided_scalar
 
 
 DATA_SDIR = 'fsaverage'
@@ -92,6 +93,33 @@ def test_morph_data():
     curv = read_morph_data(curv_path)
     assert_true(-1.0 < curv.min() < 0)
     assert_true(0 < curv.max() < 1.0)
+    with InTemporaryDirectory():
+        new_path = 'test'
+        write_morph_data(new_path, curv)
+        curv2 = read_morph_data(new_path)
+        assert_equal(curv2, curv)
+
+
+def test_write_morph_data():
+    """Test write_morph_data edge cases"""
+    values = np.arange(20, dtype='>f4')
+    okay_shapes = [(20,), (20, 1), (20, 1, 1), (1, 20)]
+    bad_shapes = [(10, 2), (1, 1, 20, 1, 1)]
+    big_num = np.iinfo('i4').max + 1
+    with InTemporaryDirectory():
+        for shape in okay_shapes:
+            write_morph_data('test.curv', values.reshape(shape))
+            # Check ordering is preserved, regardless of shape
+            assert_equal(values, read_morph_data('test.curv'))
+        assert_raises(ValueError, write_morph_data, 'test.curv',
+                      np.zeros(shape), big_num)
+        # Windows 32-bit overflows Python int
+        if np.dtype(np.int) != np.dtype(np.int32):
+            assert_raises(ValueError, write_morph_data, 'test.curv',
+                          strided_scalar((big_num,)))
+        for shape in bad_shapes:
+            assert_raises(ValueError, write_morph_data, 'test.curv',
+                          values.reshape(shape))
 
 
 @freesurfer_test