ENH: Prototype X5 io module and x5<->flirt converter program.

fsl/scripts/fsl_convert_x5.py

+#!/usr/bin/env python
+#
+# fsl_convert_x5.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+
+import os.path as op
+import            sys
+import            shutil
+import            logging
+import            argparse
+
+import fsl.data.image      as fslimage
+import fsl.utils.transform as transform
+
+
+log = logging.getLogger(__name__)
+
+
+def parseArgs(args):
+
+    parser     = argparse.ArgumentParser('fsl_convert_x5')
+    subparsers = parser.add_subparsers(dest='ctype')
+    flirt      = subparsers.add_parser('flirt')
+
+    flirt.add_argument('input')
+    flirt.add_argument('output')
+    flirt.add_argument('-s', '--source')
+    flirt.add_argument('-r', '--reference')
+    flirt.add_argument('-if', '--input_format',  choices=('x5', 'mat'))
+    flirt.add_argument('-of', '--output_format', choices=('x5', 'mat'))
+
+    args = parser.parse_args(args)
+
+    def getfmt(fname):
+        ext = op.splitext(fname)[1]
+        if ext not in ('.x5', '.mat'):
+            raise argparse.ArgumentError('Could not infer format from '
+                                         'filename: {}'.format(args.input))
+        return ext[1:]
+
+    if args.ctype == 'flirt':
+        if args.input_format  is None: args.input_format  = getfmt(args.input)
+        if args.output_format is None: args.output_format = getfmt(args.output)
+
+    return args
+
+
+def flirtToX5(args):
+    src   = fslimage.Image(args.source)
+    ref   = fslimage.Image(args.reference)
+    xform = transform.readFlirt(args.input)
+    transform.writeFlirtX5(args.output, xform, src, ref)
+
+
+def X5ToFlirt(args):
+    xform, src, ref = transform.readFlirtX5(args.input)
+    xform           = transform.toFlirt(xform, src, ref, 'world', 'world')
+    transform.writeFlirt(xform, args.output)
+
+
+def main(args=None):
+
+    if args is None:
+        args = sys.argv[1:]
+
+    args   = parseArgs(args)
+    ctype  = args.ctype
+
+    if ctype == 'flirt':
+        infmt  = args.input_format
+        outfmt = args.output_format
+
+        if   (infmt, outfmt) == ('x5', 'mat'): X5ToFlirt(args)
+        elif (infmt, outfmt) == ('mat', 'x5'): flirtToX5(args)
+        else: shutil.copy(args.input, args.output)
+
+
+if __name__ == '__main__':
+    sys.exit(main())

fsl/utils/transform/__init__.py

@@ -25,8 +25,16 @@ from .affine import (  # noqa
     transform,
     transformNormal,
     rmsdev)
+
 from .flirt  import (  # noqa
+    readFlirt,
+    writeFlirt,
     fromFlirt,
     toFlirt,
     flirtMatrixToSform,
     sformToFlirtMatrix)
+
+from .x5 import (  # noqa
+    readFlirtX5,
+    writeFlirtX5
+)

fsl/utils/transform/x5.py

+#!/usr/bin/env python
+#
+# x5.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+"""
+
+Functions for reading/writing linear/non-linear FSL transformations from/to
+BIDS X5 files.
+"""
+
+import json
+
+import numpy        as np
+import numpy.linalg as npla
+import nibabel      as nib
+import h5py
+
+from . import flirt
+
+
+def _writeLinearTransform(group, xform):
+    group.attrs['Type'] = 'linear'
+    group.create_dataset('Transform', data=xform)
+    group.create_dataset('Inverse',   data=npla.inv(xform))
+
+
+def _readLinearTransform(group):
+    if group.attrs['Type'] != 'linear':
+        raise ValueError('Not a linear transform')
+    return np.array(group['Transform'])
+
+
+def _writeLinearMapping(group, img):
+    group.attrs['Type']   = 'image'
+    group.attrs['Size']   = img.shape[ :3]
+    group.attrs['Scales'] = img.pixdim[:3]
+
+    mapping = group.create_group('Mapping')
+    _writeLinearTransform(mapping, img.getAffine('voxel', 'world'))
+
+def _readLinearMapping(group):
+
+    import fsl.data.image as fslimage
+
+    if group.attrs['Type'] != 'image':
+        raise ValueError('Not an image mapping')
+
+    shape  = group.attrs['Size']
+    pixdim = group.attrs['Scales']
+    xform  = _readLinearTransform(group['Mapping'])
+
+    hdr = nib.Nifti2Header()
+    hdr.set_data_shape(shape)
+    hdr.set_zooms(     pixdim)
+    hdr.set_sform(     xform, 'aligned')
+    return fslimage.Nifti(hdr)
+
+
+def writeFlirtX5(fname, xform, src, ref):
+    """
+    """
+
+    xform = flirt.fromFlirt(xform, src, ref, 'world', 'world')
+
+    with h5py.File(fname, 'w') as f:
+        f.attrs['Format']   = 'X5'
+        f.attrs['Version']  = '0.0.1'
+        f.attrs['Metadata'] = json.dumps({'software' : 'flirt'})
+
+        _writeLinearTransform(f, xform)
+
+        from_ = f.create_group('/From')
+        to    = f.create_group('/To')
+
+        _writeLinearMapping(from_, src)
+        _writeLinearMapping(to,    ref)
+
+
+def readFlirtX5(fname):
+    """
+    """
+    with h5py.File(fname, 'r') as f:
+        xform = _readLinearTransform(f['/'])
+        src   = _readLinearMapping(  f['/From'])
+        ref   = _readLinearMapping(  f['/To'])
+
+    return xform, src, ref