fsl/tests/test_scripts/test_imln.py

+#!/usr/bin/env python
+#
+# test_imln.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+import os
+import os.path as op
+from unittest import mock
+
+import pytest
+
+from   fsl.utils.tempdir import tempdir
+import fsl.utils.path    as     fslpath
+import fsl.scripts.imln  as     imln
+
+from fsl.tests import touch
+
+
+def test_usage():
+    assert imln.main([]) != 0
+    with mock.patch('sys.argv', []):
+        assert imln.main() != 0
+
+
+def test_imln():
+
+    # (files, command, expected links)
+    tests = [
+        ('a.nii',    'a.nii    link.nii',    'link.nii'),
+        ('a.nii',    'a        link',        'link.nii'),
+        ('a.nii',    'a.nii    link',        'link.nii'),
+        ('a.nii',    'a        link.nii',    'link.nii'),
+        ('a.nii',    'a        link.nii.gz', 'link.nii'),
+        ('a.nii.gz', 'a.nii.gz link.nii.gz', 'link.nii.gz'),
+        ('a.nii.gz', 'a        link',        'link.nii.gz'),
+        ('a.nii.gz', 'a.nii.gz link',        'link.nii.gz'),
+        ('a.nii.gz', 'a        link.nii.gz', 'link.nii.gz'),
+        ('a.nii.gz', 'a        link.nii',    'link.nii.gz'),
+
+        ('a.img a.hdr', 'a      link',     'link.img link.hdr'),
+        ('a.img a.hdr', 'a      link.img', 'link.img link.hdr'),
+        ('a.img a.hdr', 'a      link.hdr', 'link.img link.hdr'),
+        ('a.img a.hdr', 'a.img  link',     'link.img link.hdr'),
+        ('a.img a.hdr', 'a.hdr  link',     'link.img link.hdr'),
+        ('a.img a.hdr', 'a.img  link.img', 'link.img link.hdr'),
+        ('a.img a.hdr', 'a.hdr  link.hdr', 'link.img link.hdr'),
+        ('a.img a.hdr', 'a.img  link.hdr', 'link.img link.hdr'),
+        ('a.img a.hdr', 'a.hdr  link.img', 'link.img link.hdr'),
+    ]
+
+
+    for files, command, explinks in tests:
+        with tempdir():
+            files    = files.split()
+            command  = command.split()
+            explinks = explinks.split()
+
+            for f in files:
+                touch(f)
+
+            assert imln.main(command) == 0
+
+            assert sorted(os.listdir('.')) == sorted(files + explinks)
+
+            for f, l in zip(sorted(files), sorted(explinks)):
+                assert op.islink(l)
+                assert op.isfile(f) and not op.islink(f)
+                assert op.realpath(l) == op.abspath(f)
+
+
+    # subdirs - imln currently only
+    # works with absolute paths (we
+    # make all paths absolute below)
+    tests = [
+        ('dir/a.nii',           'dir/a dir/link', 'dir/link.nii'),
+        ('dir/a.img dir/a.hdr', 'dir/a dir/link', 'dir/link.img dir/link.hdr'),
+        ('dir/a.nii',           'dir/a link',     'link.nii'),
+        ('dir/a.img dir/a.hdr', 'dir/a link',     'link.img link.hdr'),
+        ('a.nii',               'a     dir/link', 'dir/link.nii'),
+        ('a.img a.hdr',         'a     dir/link', 'dir/link.img dir/link.hdr'),
+    ]
+    for files, command, explinks in tests:
+        with tempdir():
+            files    = files.split()
+            command  = [op.abspath(c) for c in command.split()]
+            explinks = explinks.split()
+
+            os.mkdir('dir')
+
+            for f in files:
+                touch(f)
+
+            assert imln.main(command) == 0
+
+            for f, l in zip(sorted(files), sorted(explinks)):
+                assert op.islink(l)
+                assert op.isfile(f) and not op.islink(f)
+                assert op.realpath(l) == op.abspath(f)
+
+    # error cases
+    # (files, commnad)
+    tests = [
+        ('a.img',             'a link'),
+        ('a.nii a.img a.hdr', 'a link'),
+    ]
+    for files, command in tests:
+        with tempdir():
+            files   = files.split()
+            command = command.split()
+
+            for f in files:
+                touch(f)
+
+            assert imln.main(command) != 0
+            assert sorted(os.listdir('.')) == sorted(files)

tests/test_immv_imcp.py → fsl/tests/test_scripts/test_immv_imcp.py

 #!/usr/bin/env python
-#
-# test_immv_imcp.py -
-#
-# Author: Paul McCarthy <pauldmccarthy@gmail.com>
-#
-
-from __future__ import print_function
 
 
 
@@ -13,101 +6,25 @@ import os.path    as op
 import itertools  as it
 import subprocess as sp
 import               os
+import               glob
+import               gzip
 import               shutil
 import               tempfile
 
+from unittest import mock
+
 import               pytest
 
+import numpy as np
 import nibabel as nib
 
-import fsl.utils.imcp   as imcp
-import fsl.scripts.imcp as imcp_script
-import fsl.scripts.immv as immv_script
-import fsl.data.image   as fslimage
-
-from . import make_random_image
-from . import make_dummy_file
-from . import looks_like_image
-from . import cleardir
-
-
-real_print = print
-
-def print(*args, **kwargs):
-    pass
-
-
-def makeImage(filename):
-    return hash(make_random_image(filename).get_data().tobytes())
-
-
-def checkImageHash(filename, datahash):
-    """Checks that the given NIFTI image matches the given hash.
-    """
-
-    img = nib.load(filename)
-    assert hash(img.get_data().tobytes()) == datahash
-
-
-def checkFilesToExpect(files, outdir, outputType, datahashes):
-
-    exts = {
-        'NIFTI'      : ['.nii'],
-        'NIFTI_PAIR' : ['.hdr', '.img'],
-        'NIFTI_GZ'   : ['.nii.gz'],
-        ''           : ['.nii.gz'],
-    }.get(outputType, None)
-
-    allFiles = []
-
-    if isinstance(files, str):
-        files = files.split()
+from   fsl.utils.tempdir import tempdir
+import fsl.scripts.imcp      as imcp_script
+import fsl.scripts.immv      as immv_script
+import fsl.data.image        as fslimage
 
-    for f in files:
-
-        f, fe = fslimage.splitExt(f)
-        fexts = exts
-
-        if fexts is None:
-            fexts = {
-                '.img'    : ['.hdr', '.img'],
-                '.hdr'    : ['.hdr', '.img'],
-                '.nii'    : ['.nii'],
-                '.nii.gz' : ['.nii.gz']
-            }.get(fe, [])
-        # filename already has a different extension
-        elif fe != '' and fe not in fexts:
-            fexts = [fe]
-
-        for e in fexts:
-
-            expected = op.join(outdir, f + e)
-
-            allFiles.append(expected)
-
-            print('  ', expected)
-
-            assert op.exists(expected)
-
-    allThatExist = os.listdir(outdir)
-    allThatExist = [f for f in allThatExist if op.isfile(op.join(outdir, f))]
-
-    assert len(allThatExist) == len(allFiles)
-
-    for i, f in enumerate(files):
-        f = fslimage.addExt(op.join(outdir, f), mustExist=True)
-
-        if isinstance(datahashes, list):
-            if len(datahashes) > len(files):
-                diff = len(datahashes) - len(files)
-                h    = datahashes[i + diff]
-
-            else:
-                h = datahashes[i]
-        else:
-            h = datahashes[op.basename(f)]
-
-        checkImageHash(f, h)
+from fsl.tests import cleardir
+from fsl.tests.test_immv_imcp import makeImage, checkImageHash, checkFilesToExpect
 
 
 def test_imcp_script_shouldPass(move=False):
@@ -313,80 +230,67 @@ def test_imcp_script_shouldPass(move=False):
 
         for outputType, reldir in it.product(outputTypes, reldirs):
 
-            os.environ['FSLOUTPUTTYPE'] = outputType
-
-            for files_to_create, imcp_args, files_to_expect in tests:
-
-                imageHashes = []
-
-                print()
-                print('files_to_create: ', files_to_create)
-                print('imcp_args:       ', imcp_args)
-                print('files_to_expect: ', files_to_expect)
+            with mock.patch.dict(os.environ, {'FSLOUTPUTTYPE' : outputType}):
 
-                for i, fname in enumerate(files_to_create.split()):
-                    imageHashes.append(makeImage(op.join(indir, fname)))
+                for files_to_create, imcp_args, files_to_expect in tests:
 
-                imcp_args = imcp_args.split()
+                    imageHashes = []
 
-                tindir  = indir
-                toutdir = outdir
+                    for i, fname in enumerate(files_to_create.split()):
+                        imageHashes.append(makeImage(op.join(indir, fname)))
 
-                if   reldir == 'neutral': reldir = startdir
-                elif reldir == 'indir':   reldir = tindir
-                elif reldir == 'outdir':  reldir = toutdir
-                elif reldir == 'samedir':
-                    reldir  = tindir
-                    toutdir = tindir
+                    imcp_args = imcp_args.split()
 
-                    if not move:
+                    tindir  = indir
+                    toutdir = outdir
 
-                        infiles = os.listdir(tindir)
-
-                        files_to_expect = files_to_expect +  ' ' + \
-                                          ' '.join(infiles)
+                    if   reldir == 'neutral': reldir = startdir
+                    elif reldir == 'indir':   reldir = tindir
+                    elif reldir == 'outdir':  reldir = toutdir
+                    elif reldir == 'samedir':
+                        reldir  = tindir
+                        toutdir = tindir
 
-                        for inf in infiles:
-                            img     = nib.load(op.join(tindir, inf),
-                                               mmap=False)
-                            imghash = hash(img.get_data().tobytes())
-                            img = None
-                            imageHashes.append(imghash)
+                        if not move:
 
-                print('adj files_to_expect: ', files_to_expect)
+                            infiles = os.listdir(tindir)
 
-                os.chdir(reldir)
+                            files_to_expect = files_to_expect +  ' ' + \
+                                              ' '.join(infiles)
 
-                imcp_args[:-1] = [op.join(tindir, a) for a in imcp_args[:-1]]
-                imcp_args[ -1] =  op.join(toutdir, imcp_args[-1])
+                            for inf in infiles:
+                                img     = nib.load(op.join(tindir, inf),
+                                                   mmap=False)
+                                imghash = hash(np.asanyarray(img.dataobj).tobytes())
+                                img = None
+                                imageHashes.append(imghash)
 
-                for i, a in enumerate(imcp_args):
-                    if op.splitdrive(a)[0] == op.splitdrive(reldir)[0]:
-                        imcp_args[i] = op.relpath(a, reldir)
+                    os.chdir(reldir)
 
-                print('indir before:    ', os.listdir(tindir))
-                print('outdir before:   ', os.listdir(toutdir))
+                    imcp_args[:-1] = [op.join(tindir, a) for a in imcp_args[:-1]]
+                    imcp_args[ -1] =  op.join(toutdir, imcp_args[-1])
 
-                if move: result = immv_script.main(imcp_args)
-                else:    result = imcp_script.main(imcp_args)
+                    for i, a in enumerate(imcp_args):
+                        if op.splitdrive(a)[0] == op.splitdrive(reldir)[0]:
+                            imcp_args[i] = op.relpath(a, reldir)
 
-                print('indir after:     ', os.listdir(tindir))
-                print('outdir after:    ', os.listdir(toutdir))
+                    if move: result = immv_script.main(imcp_args)
+                    else:    result = imcp_script.main(imcp_args)
 
-                assert result == 0
+                    assert result == 0
 
-                checkFilesToExpect(
-                    files_to_expect, toutdir, outputType, imageHashes)
+                    checkFilesToExpect(
+                        files_to_expect, toutdir, outputType, imageHashes)
 
-                # too hard if indir == outdir
-                if move and tindir != toutdir:
-                    infiles = os.listdir(tindir)
-                    infiles = [f for f in infiles if op.isfile(f)]
-                    infiles = [f for f in infiles if op.isfile(f)]
-                    assert len(infiles) == 0
+                    # too hard if indir == outdir
+                    if move and tindir != toutdir:
+                        infiles = os.listdir(tindir)
+                        infiles = [f for f in infiles if op.isfile(f)]
+                        infiles = [f for f in infiles if op.isfile(f)]
+                        assert len(infiles) == 0
 
-                cleardir(indir)
-                cleardir(outdir)
+                    cleardir(indir)
+                    cleardir(outdir)
 
     finally:
         os.chdir(startdir)
@@ -476,18 +380,9 @@ def test_imcp_script_shouldFail(move=False):
                     cmd = cmd.replace('indir', indir).replace('outdir', outdir)
                     sp.call(cmd.split())
 
-            print('calling {} {}'.format('immv' if move else 'imcp',
-                                         ' '.join(imcp_args)))
-
-            print('indir before:   {}'.format(os.listdir(indir)))
-            print('out dir before: {}'.format(os.listdir(outdir)))
-
             if move: result = immv_script.main(imcp_args)
             else:    result = imcp_script.main(imcp_args)
 
-            print('indir after:   {}'.format(os.listdir(indir)))
-            print('out dir after: {}'.format(os.listdir(outdir)))
-
             assert result != 0
 
             sp.call('chmod u+rwx {}'.format(indir) .split())
@@ -509,7 +404,6 @@ def test_imcp_script_shouldFail(move=False):
     else:    assert imcp_script.main(['wa']) != 0
 
 
-
 def test_immv_script_shouldPass():
     test_imcp_script_shouldPass(move=True)
 
@@ -519,214 +413,123 @@ def test_immv_script_shouldFail():
     test_imcp_script_shouldFail(move=True)
 
 
-def test_imcp_shouldPass(move=False):
+def test_imcp_badExt():
+    with tempdir():
 
+        ihash  = makeImage('file.nii.gz')
+        result = imcp_script.main(['file.nii', 'dest'])
 
-    #
-    # if not useDefaultExt:
-    #
-    #      imcp src.img dest     -> dest.img
-    #      imcp src.img dest.nii -> dest.nii
-    #
-    # if defaultExt:
-    #      imgp src.img dest     -> dest.nii.gz
-    #      imgp src.img dest.nii -> dest.nii.gz
+        assert result == 0
+        assert op.exists('dest.nii.gz')
+        checkImageHash('dest.nii.gz', ihash)
 
-    #
-    # (files_to_create,
-    #    [( imcp_args, files_which_should_exist),
-    #     ...
-    #    ]
-    # )
-    shouldPass = [
-        ('file.img', [
-            ('file     file',     'file.img'),
-            ('file     file.img', 'file.img'),
-            ('file     file.hdr', 'file.img'),
-            ('file     .',        'file.img'),
-            ('file.img file',     'file.img'),
-            ('file.img file.img', 'file.img'),
-            ('file.img file.hdr', 'file.img'),
-            ('file.img .',        'file.img'),
-            ('file.hdr file',     'file.img'),
-            ('file.hdr file.img', 'file.img'),
-            ('file.hdr file.hdr', 'file.img'),
-            ('file.hdr .',        'file.img'),
-        ]),
-
-        ('file.img file.blob', [
-            ('file     file',     'file.img'),
-            ('file     file.img', 'file.img'),
-            ('file     file.hdr', 'file.img'),
-            ('file     .',        'file.img'),
-            ('file.img file',     'file.img'),
-            ('file.img file.img', 'file.img'),
-            ('file.img file.hdr', 'file.img'),
-            ('file.img .',        'file.img'),
-            ('file.hdr file',     'file.img'),
-            ('file.hdr file.img', 'file.img'),
-            ('file.hdr file.hdr', 'file.img'),
-            ('file.hdr .',        'file.img'),
-        ]),
-
-
-        ('file.img file.nii', [
-            ('file.img file',     'file.img'),
-            ('file.img file.img', 'file.img'),
-            ('file.img file.hdr', 'file.img'),
-            ('file.img .',        'file.img'),
-            ('file.hdr file',     'file.img'),
-            ('file.hdr file.img', 'file.img'),
-            ('file.hdr file.hdr', 'file.img'),
-            ('file.hdr .',        'file.img'),
-            ('file.nii file',     'file.nii'),
-            ('file.nii file.nii', 'file.nii'),
-
-            # TODO both img/nii files
-        ]),
-
-
-        ('file.nii', [
-            ('file     file',     'file.nii'),
-            ('file     file.nii', 'file.nii'),
-            ('file     .',        'file.nii'),
-            ('file.nii file',     'file.nii'),
-            ('file.nii file.nii', 'file.nii'),
-            ('file.nii .',        'file.nii'),
-        ]),
-
-        ('file.nii.gz', [
-            ('file        file',        'file.nii.gz'),
-            ('file        file.nii.gz', 'file.nii.gz'),
-            ('file        .',           'file.nii.gz'),
-            ('file.nii.gz file',        'file.nii.gz'),
-            ('file.nii.gz file.nii.gz', 'file.nii.gz'),
-            ('file.nii.gz .',           'file.nii.gz'),
-        ]),
-
-
-        ('file.nii file.blob', [
-            ('file     file',     'file.nii'),
-            ('file     file.nii', 'file.nii'),
-            ('file     .',        'file.nii'),
-            ('file.nii file',     'file.nii'),
-            ('file.nii file.nii', 'file.nii'),
-            ('file.nii .',        'file.nii'),
-        ]),
-
-
-        ('file.nii file.nii.gz', [
-            ('file.nii    file',        'file.nii'),
-            ('file.nii    file.nii',    'file.nii'),
-            ('file.nii    .',           'file.nii'),
-            ('file.nii.gz file',        'file.nii.gz'),
-            ('file.nii.gz file.nii.gz', 'file.nii.gz'),
-            ('file.nii.gz .',           'file.nii.gz'),
-
-            # TODO both
-        ]),
-
-        ('file.img file.nii file.nii.gz', [
-            ('file.img file.nii file.nii.gz .', 'file.img file.nii file.nii.gz'),
-            ('file.img                      .', 'file.img'),
-            ('file.img                      .', 'file.img'),
-            ('file.nii                      .', 'file.nii'),
-            ('file.nii file.nii.gz          .', 'file.nii file.nii.gz'),
-        ]),
-
-
-        ('001.img 002.img 003.img', [
-
-            ('001     002     003     .', '001.img 002.img 003.img'),
-            ('001.img 002.img 003.img .', '001.img 002.img 003.img'),
-            ('001.hdr 002.hdr 003.hdr .', '001.img 002.img 003.img'),
-
-            ('001.img 002     003     .', '001.img 002.img 003.img'),
-            ('001.hdr 002     003     .', '001.img 002.img 003.img'),
-
-            ('001.img 002.hdr 003.img .', '001.img 002.img 003.img'),
-            ('001.hdr 002.img 003.hdr .', '001.img 002.img 003.img'),
-
-            ('001     003             .', '001.img 003.img'),
-            ('001.img 003.img         .', '001.img 003.img'),
-            ('001.hdr 003.hdr         .', '001.img 003.img'),
-
-            ('001.img 003             .', '001.img 003.img'),
-            ('001.hdr 003             .', '001.img 003.img'),
-
-            ('001.img 003.img         .', '001.img 003.img'),
-            ('001.hdr 003.hdr         .', '001.img 003.img'),
-        ]),
-    ]
 
+def test_immv_badExt():
+    with tempdir():
 
-    indir  = tempfile.mkdtemp()
-    outdir = tempfile.mkdtemp()
+        ihash  = makeImage('file.nii.gz')
+        result = immv_script.main(['file.nii', 'dest'])
 
-    try:
+        assert result == 0
+        assert op.exists('dest.nii.gz')
+        checkImageHash('dest.nii.gz', ihash)
 
-        for files_to_create, tests in shouldPass:
 
-            files_to_create = files_to_create.split()
 
-            for imcp_args, should_exist in tests:
+def _make_file(prefix, ftype, dtype):
 
-                should_exist    = should_exist.split()
-                imcp_args       = imcp_args.split()
-                imcp_srcs       = imcp_args[:-1]
-                imcp_dest       = imcp_args[ -1]
+    mapping = {
+        fslimage.FileType.NIFTI          : (nib.Nifti1Image,  'nii'),
+        fslimage.FileType.NIFTI2         : (nib.Nifti2Image,  'nii'),
+        fslimage.FileType.ANALYZE        : (nib.AnalyzeImage, 'img'),
+        fslimage.FileType.NIFTI_PAIR     : (nib.Nifti1Pair,   'img'),
+        fslimage.FileType.NIFTI2_PAIR    : (nib.Nifti2Pair,   'img'),
+        fslimage.FileType.ANALYZE_GZ     : (nib.AnalyzeImage, 'img.gz'),
+        fslimage.FileType.NIFTI_GZ       : (nib.Nifti1Image,  'nii.gz'),
+        fslimage.FileType.NIFTI2_GZ      : (nib.Nifti2Image,  'nii.gz'),
+        fslimage.FileType.NIFTI_PAIR_GZ  : (nib.Nifti1Pair,   'img.gz'),
+        fslimage.FileType.NIFTI2_PAIR_GZ : (nib.Nifti2Pair,   'img.gz'),
+    }
 
-                hashes = {}
-                for fn in files_to_create:
-                    if looks_like_image(fn):
-                        hashes[fn] = makeImage(op.join(indir, fn))
-                    else:
-                        hashes[fn] = make_dummy_file(op.join(indir, fn))
+    if np.issubdtype(dtype, np.complex64):
+        data = np.random.random((20, 20, 20)).astype(np.float32) + \
+               np.random.random((20, 20, 20)).astype(np.float32) * 1j
+    else:
+        data = np.random.random((20, 20, 20)).astype(dtype)
 
-                print()
-                print('files_to_create: ', files_to_create)
-                print('imcp_srcs:       ', imcp_srcs)
-                print('imcp_dest:       ', imcp_dest)
-                print('should_exist:    ', should_exist)
-                print('indir:           ', os.listdir(indir))
+    cls, suffix = mapping[ftype]
+    filename    = f'{prefix}.{suffix}'
 
-                for src in imcp_srcs:
+    cls(data, None).to_filename(filename)
 
-                    print('  src: {}'.format(src))
+    return filename
 
-                    src = op.join(indir, src)
 
-                    if move: imcp.immv(src, op.join(outdir, imcp_dest), overwrite=True)
-                    else:    imcp.imcp(src, op.join(outdir, imcp_dest), overwrite=True)
+def _is_gzip(filename):
+    try:
+        with gzip.GzipFile(filename, 'rb') as f:
+            f.read()
+        return True
+    except Exception:
+        return False
 
 
-                print('indir after:     ', os.listdir(indir))
-                print('outdir after:    ', os.listdir(outdir))
+def _is_pair(imgfile):
+    prefix, suffix = fslimage.splitExt(imgfile)
+    hdrfile        = imgfile
+    if   suffix == '.hdr':    imgfile = f'{prefix}.img'
+    elif suffix == '.hdr.gz': imgfile = f'{prefix}.img.gz'
+    elif suffix == '.img':    hdrfile = f'{prefix}.hdr'
+    elif suffix == '.img.gz': hdrfile = f'{prefix}.hdr.gz'
+    return op.exists(imgfile) and op.exists(hdrfile)
 
-                # check file contents
-                checkFilesToExpect(should_exist,
-                                   outdir,
-                                   None,
-                                   hashes)
 
-                # If move, check that
-                # input files are gone
-                if move:
-                    for f in should_exist:
-                         assert not op.exists(op.join(indir, f))
+def _is_analyze(filename):
+    img = nib.load(filename)
+    return  _is_pair(filename)                                     and \
+            isinstance(img, (nib.AnalyzeImage, nib.AnalyzeHeader)) and \
+        not isinstance(img, (nib.Nifti1Image,  nib.Nifti1Header, nib.Nifti1Pair))
 
-                for f in os.listdir(indir):
-                    try:    os.remove(op.join(indir,  f))
-                    except: pass
 
-                for f in os.listdir(outdir):
-                    os.remove(op.join(outdir, f))
+def _is_nifti1(filename):
+    img = nib.load(filename)
+    return  isinstance(img, (nib.Nifti1Image, nib.Nifti1Header, nib.Nifti1Pair)) and \
+        not isinstance(img, (nib.Nifti2Image, nib.Nifti2Header, nib.Nifti2Pair))
 
 
-    finally:
-        shutil.rmtree(indir)
-        shutil.rmtree(outdir)
+def _is_nifti2(filename):
+    img = nib.load(filename)
+    return isinstance(img, (nib.Nifti2Image, nib.Nifti2Header, nib.Nifti2Pair))
+
+
+def _check_file(prefix, expftype, dtype):
+    filename = fslimage.addExt(prefix)
+    if   'NIFTI2'  in expftype.name: assert _is_nifti2( filename)
+    elif 'NIFTI'   in expftype.name: assert _is_nifti1( filename)
+    elif 'ANALYZE' in expftype.name: assert _is_analyze(filename)
+    if   'GZ'      in expftype.name: assert _is_gzip(   filename)
+    if   'PAIR'    in expftype.name: assert _is_pair(   filename)
+
+    img = nib.load(filename)
+    assert np.issubdtype(img.get_data_dtype(), dtype)
+
+
+def test_imcp_script_correct_output_type(move=False):
+
+    # only testing dtypes supported by ANALYZE
+    dtypes = [np.uint8, np.int16, np.float32, np.float64, np.complex64]
+
+    # from, to
+    for from_, to_, dtype in it.product(fslimage.FileType, fslimage.FileType, dtypes):
+        with tempdir():
+
+            fname = f'f{from_.name}'
+            tname = f't{to_.name}'
+
+            _make_file(fname, from_, dtype)
 
+            with mock.patch.dict(os.environ, {'FSLOUTPUTTYPE' : to_.name}):
+                if move: immv_script.main([fname, tname])
+                else:    imcp_script.main([fname, tname])
 
-def test_immv_shouldPass():
-    test_imcp_shouldPass(move=True)
+            _check_file(tname, to_, dtype)

fsl/tests/test_scripts/test_imrm.py

+#!/usr/bin/env python
+#
+# test_imrm.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+import os
+
+from fsl.utils.tempdir import tempdir
+import fsl.scripts.imrm as imrm
+
+from fsl.tests import touch
+
+def test_imrm_usage():
+    assert imrm.main([]) != 0
+
+
+def test_imrm():
+    # (files present, command, expected)
+    tests = [
+        ('a.nii',       'a',       ''),
+        ('a.nii.gz',    'a',       ''),
+        ('a.img a.hdr', 'a',       ''),
+        ('a.img',       'a',       ''),
+        ('a.hdr',       'a',       ''),
+        ('a.nii b.nii', 'a',       'b.nii'),
+        ('a.nii b.nii', 'a b',     ''),
+        ('a.nii b.nii', 'a b.nii', ''),
+
+        # suffix doesn't have to be correct
+        ('a.nii.gz', 'a.nii',  ''),
+
+        # files don't exist -> no problem
+        ('a.nii', 'b',  'a.nii'),
+
+        # wildcards
+        ('img_a.nii img_b.nii', 'img_?',     ''),
+        ('img_a.nii img_b.nii', 'img_?.nii', ''),
+        ('img_a.nii img_b.nii', 'img_*',     ''),
+        ('img_a.nii img_b.nii', 'img_*.nii', ''),
+
+        ('img_a.nii img_b.nii img_cc.nii', 'img_?',     'img_cc.nii'),
+        ('img_a.nii img_b.nii img_cc.nii', 'img_?.nii', 'img_cc.nii'),
+        ('img_a.nii img_b.nii img_cc.nii', 'img_*',     ''),
+        ('img_a.nii img_b.nii img_cc.nii', 'img_*.nii', ''),
+    ]
+
+    for files, command, expected in tests:
+        with tempdir():
+
+            for f in files.split():
+                touch(f)
+
+            print('files',    files)
+            print('command',  command)
+            print('expected', expected)
+
+            ret = imrm.main(('imrm ' + command).split())
+
+            assert ret == 0
+            assert sorted(os.listdir()) == sorted(expected.split())

fsl/tests/test_scripts/test_imtest.py

+#!/usr/bin/env python
+#
+# test_imtest.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+import os
+import os.path as op
+
+import fsl.utils.path     as fslpath
+from   fsl.utils.tempdir  import tempdir
+
+import fsl.scripts.imtest as imtest
+
+from fsl.tests import CaptureStdout, touch
+
+
+def test_wrongargs():
+    cap = CaptureStdout()
+    with cap:
+        assert imtest.main([]) == 0
+    assert cap.stdout.strip() == '0'
+
+
+def test_imtest():
+    # (files, input, expected)
+    tests = [
+        ('a.nii', 'a',        '1'),
+        ('a.nii', 'a.nii',    '1'),
+        ('a.nii', 'a.nii.gz', '1'), # imtest is suffix-agnostic
+
+        ('a.img a.hdr', 'a',     '1'),
+        ('a.img a.hdr', 'a.img', '1'),
+        ('a.img a.hdr', 'a.hdr', '1'),
+        ('a.img',       'a',     '0'),
+        ('a.img',       'a.img', '0'),
+        ('a.img',       'a.hdr', '0'),
+        ('a.hdr',       'a',     '0'),
+        ('a.hdr',       'a.img', '0'),
+        ('a.hdr',       'a.hdr', '0'),
+
+        ('dir/a.nii',           'dir/a', '1'),
+        ('dir/a.img dir/a.hdr', 'dir/a', '1'),
+    ]
+
+    for files, input, expected in tests:
+        with tempdir():
+            for f in files.split():
+                dirname = op.dirname(f)
+                if dirname != '':
+                    os.makedirs(dirname, exist_ok=True)
+                touch(f)
+
+            cap = CaptureStdout()
+            with cap:
+                assert imtest.main([input]) == 0
+
+            assert cap.stdout.strip() == expected
+
+    # test that sym-links are
+    # followed correctly
+    with tempdir():
+        touch('image.nii.gz')
+        os.symlink('image.nii.gz', 'link.nii.gz')
+        cap = CaptureStdout()
+        with cap:
+            assert imtest.main(['link']) == 0
+        assert cap.stdout.strip() == '1'
+
+    # sym-links in sub-directories
+    # (old imtest would not work
+    # in this scenario)
+    with tempdir():
+        os.mkdir('subdir')
+        impath = op.join('subdir', 'image.nii.gz')
+        lnpath = op.join('subdir', 'link.nii.gz')
+        touch(impath)
+        os.symlink('image.nii.gz', lnpath)
+        cap = CaptureStdout()
+        with cap:
+            assert imtest.main([lnpath]) == 0
+        assert cap.stdout.strip() == '1'

fsl/tests/test_scripts/test_remove_ext.py

+#!/usr/bin/env python
+#
+# test_remove_ext.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+
+import fsl.scripts.remove_ext as remove_ext
+
+from fsl.tests import CaptureStdout
+
+def test_usage():
+    assert remove_ext.main([]) != 0
+
+def test_remove_ext():
+    # (input, expected output)
+    tests = [
+        ('a',                 'a'),
+        ('a.nii',             'a'),
+        ('a.nii.gz',          'a'),
+        ('a.txt',             'a.txt'),
+        ('a.nii b.img c.hdr', 'a b c'),
+        ('a.nii b.img b.hdr', 'a b b'),
+        ('a b.img c.txt',     'a b c.txt'),
+        ('a.nii.gz b c.mnc',  'a b c'),
+    ]
+
+    for input, expected in tests:
+
+        cap = CaptureStdout()
+        with cap:
+            ret = remove_ext.main(input.split())
+
+        assert ret == 0
+
+        got = cap.stdout.split()
+        assert sorted(got) == sorted(expected.split())

fsl/tests/test_scripts/test_resample_image.py

+#!/usr/bin/env python
+
+
+import numpy as np
+
+import pytest
+
+import fsl.scripts.resample_image as resample_image
+
+
+import fsl.transform.affine as affine
+from fsl.utils.tempdir import tempdir
+from fsl.data.image    import Image
+
+from fsl.tests import make_random_image
+
+
+
+def test_resample_image_shape():
+    with tempdir():
+        img = Image(make_random_image('image.nii.gz', dims=(10, 10, 10)))
+        resample_image.main('image resampled -s 20,20,20'.split())
+        res = Image('resampled')
+
+        expv2w = affine.concat(
+            img.voxToWorldMat,
+            affine.scaleOffsetXform([0.5, 0.5, 0.5], 0))
+
+        assert np.all(np.isclose(res.shape, (20, 20, 20)))
+        assert np.all(np.isclose(res.pixdim, (0.5, 0.5, 0.5)))
+        assert np.all(np.isclose(res.voxToWorldMat, expv2w))
+        assert np.all(np.isclose(
+            np.array(affine.axisBounds(res.shape, res.voxToWorldMat)) - 0.25,
+                     affine.axisBounds(img.shape, img.voxToWorldMat)))
+
+        resample_image.main('image resampled -s 20,20,20 -o corner'.split())
+        res = Image('resampled')
+        assert np.all(np.isclose(
+            affine.axisBounds(res.shape, res.voxToWorldMat),
+            affine.axisBounds(img.shape, img.voxToWorldMat)))
+
+
+def test_resample_image_shape_4D():
+    with tempdir():
+
+        # Can specify three dims
+        img = Image(make_random_image('image.nii.gz', dims=(10, 10, 10, 10)))
+        resample_image.main('image resampled -s 20,20,20'.split())
+        res = Image('resampled')
+
+        assert np.all(np.isclose(res.shape, (20, 20, 20, 10)))
+        assert np.all(np.isclose(res.pixdim, (0.5, 0.5, 0.5, 1)))
+
+        # Or resample along the higher dims
+        resample_image.main('image resampled -s 20,20,20,20'.split())
+        res = Image('resampled')
+        assert np.all(np.isclose(res.shape, (20, 20, 20, 20)))
+        assert np.all(np.isclose(res.pixdim, (0.5, 0.5, 0.5, 0.5)))
+
+
+def test_resample_image_dim():
+    with tempdir():
+        img = Image(make_random_image('image.nii.gz', dims=(10, 10, 10)))
+
+        resample_image.main('image resampled -d 0.5,0.5,0.5'.split())
+
+        res = Image('resampled')
+        expv2w = affine.concat(
+            img.voxToWorldMat,
+            affine.scaleOffsetXform([0.5, 0.5, 0.5], 0))
+
+        assert np.all(np.isclose(res.shape, (20, 20, 20)))
+        assert np.all(np.isclose(res.pixdim, (0.5, 0.5, 0.5)))
+        assert np.all(np.isclose(res.voxToWorldMat, expv2w))
+
+
+def test_resample_image_ref():
+    with tempdir():
+        img = Image(make_random_image('image.nii.gz', dims=(10, 10, 10)))
+        ref = Image(make_random_image('ref.nii.gz',   dims=(20, 20, 20),
+                                      pixdims=(0.5, 0.5, 0.5)))
+
+        resample_image.main('image resampled -r ref'.split())
+
+        res    = Image('resampled')
+        expv2w = ref.voxToWorldMat
+
+        assert np.all(np.isclose(res.shape, (20, 20, 20)))
+        assert np.all(np.isclose(res.pixdim, (0.5, 0.5, 0.5)))
+        assert np.all(np.isclose(res.voxToWorldMat, expv2w))
+
+        # 3D / 4D
+        img = Image(make_random_image('image.nii.gz', dims=(10, 10, 10)))
+        ref = Image(make_random_image('ref.nii.gz',   dims=(20, 20, 20, 20),
+                                      pixdims=(0.5, 0.5, 0.5, 1)))
+
+        resample_image.main('image resampled -r ref'.split())
+        res    = Image('resampled')
+        assert np.all(np.isclose(res.shape, (20, 20, 20)))
+        assert np.all(np.isclose(res.pixdim, (0.5, 0.5, 0.5)))
+
+        # 4D / 3D
+        img = Image(make_random_image('image.nii.gz', dims=(10, 10, 10, 10)))
+        ref = Image(make_random_image('ref.nii.gz',   dims=(20, 20, 20),
+                                      pixdims=(0.5, 0.5, 0.5)))
+
+        resample_image.main('image resampled -r ref'.split())
+        res    = Image('resampled')
+        assert np.all(np.isclose(res.shape, (20, 20, 20, 10)))
+        assert np.all(np.isclose(res.pixdim, (0.5, 0.5, 0.5, 1)))
+
+        # 4D / 4D - no resampling along fourth dim
+        img = Image(make_random_image('image.nii.gz', dims=(10, 10, 10, 10)))
+        ref = Image(make_random_image('ref.nii.gz',   dims=(20, 20, 20, 20),
+                                      pixdims=(0.5, 0.5, 0.5, 1)))
+
+        resample_image.main('image resampled -r ref'.split())
+        res    = Image('resampled')
+        assert np.all(np.isclose(res.shape, (20, 20, 20, 10)))
+        assert np.all(np.isclose(res.pixdim, (0.5, 0.5, 0.5, 1)))
+
+
+def test_resample_image_bad_options():
+    with tempdir():
+        img = Image(make_random_image('image.nii.gz', dims=(10, 10, 10)))
+
+        # No args - should print help and exit(0)
+        with pytest.raises(SystemExit) as e:
+            resample_image.main([])
+        assert e.value.code == 0
+
+        with pytest.raises(SystemExit) as e:
+            resample_image.main('image resampled -d 0.5,0.5,0.5 '
+                                '-s 20,20,20'.split())
+        assert e.value.code != 0
+
+        with pytest.raises(SystemExit) as e:
+            resample_image.main('image resampled -s 20,20'.split())
+        assert e.value.code != 0
+
+        with pytest.raises(SystemExit) as e:
+            resample_image.main('image resampled -s 20,20,20,20'.split())
+        assert e.value.code != 0

fsl/tests/test_scripts/test_vest2text.py

+#!/usr/bin/env python
+#
+# test_vest2text.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+import textwrap as tw
+import numpy as np
+
+import fsl.data.vest         as fslvest
+import fsl.scripts.Text2Vest as Text2Vest
+import fsl.scripts.Vest2Text as Vest2Text
+
+from fsl.tests import tempdir
+
+
+def test_usage():
+    assert Vest2Text.main([]) == 0
+    assert Text2Vest.main([]) == 0
+
+
+def test_Vest2Text():
+    with tempdir():
+        data = np.random.random((20, 10))
+        vest = fslvest.generateVest(data)
+
+        with open('data.vest', 'wt') as f:
+            f.write(vest)
+
+        assert Vest2Text.main(['data.vest', 'data.txt']) == 0
+
+        got = np.loadtxt('data.txt')
+
+        assert np.all(np.isclose(data, got))
+
+
+def test_Text2Vest():
+    with tempdir():
+        data = np.random.random((20, 10))
+
+        np.savetxt('data.txt', data)
+
+        assert Text2Vest.main(['data.txt', 'data.vest']) == 0
+
+        got = fslvest.loadVestFile('data.vest', ignoreHeader=False)
+
+        assert np.all(np.isclose(data, got))
+
+        # make sure that 1D files are treated correctly (fsl/fslpy!387)
+        with open('colvec.txt', 'wt') as f:
+            f.write('1\n2\n3\n')
+        with open('rowvec.txt', 'wt') as f:
+            f.write('1 2 3\n')
+
+        colexp = tw.dedent("""
+        /NumWaves 1
+        /NumPoints 3
+        /Matrix
+        1.000000000000
+        2.000000000000
+        3.000000000000
+        """).strip()
+        rowexp = tw.dedent("""
+        /NumWaves 3
+        /NumPoints 1
+        /Matrix
+        1.000000000000 2.000000000000 3.000000000000
+        """).strip()
+
+        assert Text2Vest.main(['colvec.txt', 'colvec.vest']) == 0
+        assert Text2Vest.main(['rowvec.txt', 'rowvec.vest']) == 0
+
+        assert open('colvec.vest', 'rt').read().strip() == colexp
+        assert open('rowvec.vest', 'rt').read().strip() == rowexp

tests/test_settings.py → fsl/tests/test_settings.py

@@ -12,16 +12,11 @@ import            pickle
 import            textwrap
 import            tempfile
 
-# python 3
-try:
-    import unittest.mock as mock
-# python 2
-except:
-    import mock
+import unittest.mock as mock
 
 import pytest
 
-import tests
+import fsl.tests          as tests
 import fsl.utils.settings as settings
 import fsl.utils.tempdir  as tempdir
 

tests/test_tempdir.py → fsl/tests/test_tempdir.py

@@ -6,6 +6,7 @@
 #
 
 
+import glob
 import os
 import os.path as op
 
@@ -53,13 +54,27 @@ def test_tempdir_changeto():
     assert op.realpath(os.getcwd()) == cwd
 
 
+def test_tempdir_prefix():
+    with tempdir.tempdir(prefix='mytempdirtest') as tdir:
+        assert op.basename(tdir).startswith('mytempdirtest')
+
+    with tempdir.tempdir() as parent:
+        with tempdir.tempdir(prefix='mytempdirtest', root='.') as tdir:
+            assert list(glob.glob(op.join(parent, 'mytempdirtest*'))) == [tdir]
+
+
 def test_tempdir_override():
     with tempdir.tempdir() as parent:
+        override = op.abspath((op.join('override', 'directory')))
+
+        os.makedirs(override)
 
         # tempdir should not create/change to
         # a new temp directory, but should
         # stay in the override directory
-        with tempdir.tempdir(override=parent):
-            assert op.realpath(os.getcwd()) == op.realpath(parent)
+        with tempdir.tempdir(override=override):
+            assert op.realpath(os.getcwd()) == op.realpath(override)
         # override should not be deleted
-        assert op.exists(parent)
+        assert op.exists(override)
+        # and we should be moved back into the original cwd
+        assert op.realpath(os.getcwd()) == op.realpath(parent)

tests/test_transform.py → fsl/tests/test_transform/test_affine.py

@@ -6,8 +6,6 @@
 #
 
 
-from __future__ import division
-
 import                 random
 import                 glob
 import os.path      as op
@@ -15,12 +13,9 @@ import itertools    as it
 import numpy        as np
 import numpy.linalg as npla
 
-import six
-
 import pytest
 
-import fsl.utils.transform as transform
-import fsl.data.image      as fslimage
+import fsl.transform.affine  as affine
 
 
 datadir = op.join(op.dirname(__file__), 'testdata')
@@ -40,8 +35,7 @@ def readlines(filename):
         #
         # Pass it [bytes, bytes, ...], and it works
         # fine.
-        if six.PY3:
-            lines = [l.encode('ascii') for l in lines]
+        lines = [l.encode('ascii') for l in lines]
 
     return lines
 
@@ -57,7 +51,7 @@ def test_invert():
 
         x      = testdata[i * 4:i * 4 + 4, 0:4]
         invx   = testdata[i * 4:i * 4 + 4, 4:8]
-        result = transform.invert(x)
+        result = affine.invert(x)
 
         assert np.all(np.isclose(invx, result))
 
@@ -87,7 +81,7 @@ def test_concat():
 
     for inputs, expected in tests:
 
-        result = transform.concat(*inputs)
+        result = affine.concat(*inputs)
 
         assert np.all(np.isclose(result, expected))
 
@@ -95,7 +89,7 @@ def test_concat():
 def test_veclength(seed):
 
     def l(v):
-        v = np.array(v, copy=False).reshape((-1, 3))
+        v = np.asarray(v).reshape((-1, 3))
         x = v[:, 0]
         y = v[:, 1]
         z = v[:, 2]
@@ -109,10 +103,10 @@ def test_veclength(seed):
         vtype = random.choice((list, tuple, np.array))
         v     = vtype(v)
 
-        assert np.isclose(transform.veclength(v), l(v))
+        assert np.isclose(affine.veclength(v), l(v))
 
     # Multiple vectors in parallel
-    result   = transform.veclength(vectors)
+    result   = affine.veclength(vectors)
     expected = l(vectors)
     assert np.all(np.isclose(result, expected))
 
@@ -122,8 +116,8 @@ def test_normalise(seed):
     vectors = -100 + 200 * np.random.random((200, 3))
 
     def parallel(v1, v2):
-        v1 = v1 / transform.veclength(v1)
-        v2 = v2 / transform.veclength(v2)
+        v1 = v1 / affine.veclength(v1)
+        v2 = v2 / affine.veclength(v2)
 
         return np.isclose(np.dot(v1, v2), 1)
 
@@ -131,16 +125,16 @@ def test_normalise(seed):
 
         vtype = random.choice((list, tuple, np.array))
         v     = vtype(v)
-        vn    = transform.normalise(v)
-        vl    = transform.veclength(vn)
+        vn    = affine.normalise(v)
+        vl    = affine.veclength(vn)
 
         assert np.isclose(vl, 1.0)
         assert parallel(v, vn)
 
     # normalise should also be able
     # to do multiple vectors at once
-    results = transform.normalise(vectors)
-    lengths = transform.veclength(results)
+    results = affine.normalise(vectors)
+    lengths = affine.veclength(results)
     pars    = np.zeros(200)
     for i in range(200):
 
@@ -153,6 +147,24 @@ def test_normalise(seed):
     assert np.all(pars)
 
 
+def test_flip():
+
+    shape        = [5, 5, 5]
+    xform        = np.diag([2, 2, 2, 1])
+    xform[:3, 3] = [-4.5, -4.5, -4.5]
+
+    allaxes = list(it.chain(
+        it.combinations([0, 1, 2], 1),
+        it.combinations([0, 1, 2], 2),
+        it.combinations([0, 1, 2], 3)))
+
+    for axes in allaxes:
+        flipped  = affine.flip(shape, xform, *axes)
+        flo, fhi = affine.axisBounds(shape, flipped, boundary=None)
+        assert np.all(np.isclose(flo, [-5, -5, -5]))
+        assert np.all(np.isclose(fhi, [ 5,  5,  5]))
+
+
 def test_scaleOffsetXform():
 
     # Test numerically
@@ -172,7 +184,7 @@ def test_scaleOffsetXform():
         expected = [[float(v) for v in l.split()] for l in expected]
         expected = np.array(expected)
 
-        result = transform.scaleOffsetXform(scales, offsets)
+        result = affine.scaleOffsetXform(scales, offsets)
 
         assert np.all(np.isclose(result, expected))
 
@@ -208,13 +220,24 @@ def test_scaleOffsetXform():
 
     for (scale, expected) in stests:
         expected = np.array(expected).reshape(4, 4)
-        result   = transform.scaleOffsetXform(scale, 0)
+        result   = affine.scaleOffsetXform(scale, 0)
         assert np.all(np.isclose(result, expected))
     for (offset, expected) in otests:
         expected = np.array(expected).reshape(4, 4)
-        result   = transform.scaleOffsetXform(1, offset)
+        result   = affine.scaleOffsetXform(1, offset)
         assert np.all(np.isclose(result, expected))
 
+    # both args are optional
+    assert np.all(np.isclose(
+        affine.scaleOffsetXform(), np.eye(4)))
+    assert np.all(np.isclose(
+        affine.scaleOffsetXform([1, 2, 3]), np.diag([1, 2, 3, 1])))
+
+    exp = np.eye(4)
+    exp[:3, 3] = [1, 2, 3]
+    assert np.all(np.isclose(
+        affine.scaleOffsetXform(None, [1, 2, 3]), exp))
+
 
 def test_compose_and_decompose():
 
@@ -227,8 +250,10 @@ def test_compose_and_decompose():
         xform                      = lines[i * 4: i * 4 + 4]
         xform                      = np.genfromtxt(xform)
 
-        scales, offsets, rotations = transform.decompose(xform)
-        result = transform.compose(scales, offsets, rotations)
+        scales, offsets, rotations, shears = affine.decompose(
+            xform, shears=True)
+
+        result = affine.compose(scales, offsets, rotations, shears=shears)
 
         assert np.all(np.isclose(xform, result, atol=1e-5))
 
@@ -236,22 +261,24 @@ def test_compose_and_decompose():
         # different rotation origin, but we test
         # explicitly passing the origin for
         # completeness
-        scales, offsets, rotations = transform.decompose(xform)
-        result = transform.compose(scales, offsets, rotations, [0, 0, 0])
+        scales, offsets, rotations, shears = affine.decompose(
+            xform, shears=True)
+        result = affine.compose(
+            scales, offsets, rotations, origin=[0, 0, 0], shears=shears)
 
         assert np.all(np.isclose(xform, result, atol=1e-5))
 
     # compose should also accept a rotation matrix
     rots = [np.pi / 5, np.pi / 4, np.pi / 3]
-    rmat  = transform.axisAnglesToRotMat(*rots)
-    xform = transform.compose([1, 1, 1], [0, 0, 0], rmat)
+    rmat  = affine.axisAnglesToRotMat(*rots)
+    xform = affine.compose([1, 1, 1], [0, 0, 0], rmat)
 
     # And the angles flag should cause decompose
     # to return the rotation matrix, instead of
     # the axis angls
-    sc,   of,   rot   = transform.decompose(xform)
-    scat, ofat, rotat = transform.decompose(xform, angles=True)
-    scaf, ofaf, rotaf = transform.decompose(xform, angles=False)
+    sc,   of,   rot   = affine.decompose(xform)
+    scat, ofat, rotat = affine.decompose(xform, angles=True)
+    scaf, ofaf, rotaf = affine.decompose(xform, angles=False)
 
     sc,   of,   rot   = np.array(sc),   np.array(of),   np.array(rot)
     scat, ofat, rotat = np.array(scat), np.array(ofat), np.array(rotat)
@@ -268,6 +295,27 @@ def test_compose_and_decompose():
     assert np.all(np.isclose(rotat, rots))
     assert np.all(np.isclose(rotaf, rmat))
 
+    # decompose should accept a 3x3
+    # affine, and return translations of 0
+    affine.decompose(xform[:3, :3])
+    sc,   of,   rot   = affine.decompose(xform[:3, :3])
+    sc,   of,   rot   = np.array(sc), np.array(of), np.array(rot)
+    assert np.all(np.isclose(sc,    [1, 1, 1]))
+    assert np.all(np.isclose(of,    [0, 0, 0]))
+    assert np.all(np.isclose(rot,   rots))
+
+
+def test_compose_scaleAtOrigin():
+
+    for _ in range(100):
+        origin = np.random.randint(-50, 50, 3)
+        rots   = -np.pi + 2  * np.pi * np.random.random(3)
+        scales = -10    + 20 *         np.random.random(3)
+        xform  = affine.compose(scales, [0, 0, 0], rots, origin=origin,
+                                scaleAtOrigin=True)
+
+        assert np.all(np.isclose(affine.transform(origin, xform), origin))
+
 
 def test_rotMatToAxisAngles(seed):
 
@@ -280,8 +328,8 @@ def test_rotMatToAxisAngles(seed):
                 -pi2 + 2 * pi2 * np.random.random(),
                 -pi  + 2 * pi  * np.random.random()]
 
-        rmat    = transform.axisAnglesToRotMat(*rots)
-        gotrots = transform.rotMatToAxisAngles(rmat)
+        rmat    = affine.axisAnglesToRotMat(*rots)
+        gotrots = affine.rotMatToAxisAngles(rmat)
 
         assert np.all(np.isclose(rots, gotrots))
 
@@ -300,9 +348,9 @@ def test_rotMatToAffine(seed):
         if np.random.random() < 0.5: origin = None
         else:                        origin = np.random.random(3)
 
-        rmat   = transform.axisAnglesToRotMat(*rots)
-        mataff = transform.rotMatToAffine(rmat, origin)
-        rotaff = transform.rotMatToAffine(rots, origin)
+        rmat   = affine.axisAnglesToRotMat(*rots)
+        mataff = affine.rotMatToAffine(rmat, origin)
+        rotaff = affine.rotMatToAffine(rots, origin)
 
         exp         = np.eye(4)
         exp[:3, :3] = rmat
@@ -339,7 +387,7 @@ def test_axisBounds():
         shape, origin, boundary, xform, expected = readTest(i)
 
         for axes in allAxes:
-            result = transform.axisBounds(shape,
+            result = affine.axisBounds(shape,
                                           xform,
                                           axes=axes,
                                           origin=origin,
@@ -361,14 +409,31 @@ def test_axisBounds():
     # US-spelling
     assert np.all(np.isclose(
         expected,
-        transform.axisBounds(
+        affine.axisBounds(
             shape, xform, origin='center', boundary=boundary)))
 
     # Bad origin/boundary values
     with pytest.raises(ValueError):
-        transform.axisBounds(shape, xform, origin='Blag', boundary=boundary)
+        affine.axisBounds(shape, xform, origin='Blag', boundary=boundary)
     with pytest.raises(ValueError):
-        transform.axisBounds(shape, xform, origin=origin, boundary='Blufu')
+        affine.axisBounds(shape, xform, origin=origin, boundary='Blufu')
+
+
+def test_mergeBounds():
+
+    bounds = [[[ 0,  0, 0], [10, 10, 10]],
+              [[ 0, -2, 5], [11,  8,  7]],
+              [[ 4,  6, 3], [12,  9,  8]],
+              [[-1,  3, 2], [ 9, 11, 21]]]
+    expect =  ((-1, -2, 0), (12, 11, 21))
+    assert affine.mergeBounds(*bounds) == expect
+
+    bounds = [[[ 0, 10], [ 0, 10], [0, 10]],
+              [[ 0, 11], [-2,  8], [5,  7]],
+              [[ 4, 12], [ 6,  9], [3,  8]],
+              [[-1,  9], [ 3, 11], [2, 21]]]
+    expect =  ((-1, 12), (-2, 11), (0, 21))
+    assert affine.mergeBounds(*bounds) == expect
 
 
 def test_transform():
@@ -381,7 +446,7 @@ def test_transform():
         mask = np.array([[1, 0, 0, 1],
                          [0, 1, 0, 1],
                          [0, 0, 1, 1],
-                         [0, 0, 0, 1]], dtype=np.bool)
+                         [0, 0, 0, 1]], dtype=bool)
 
         return np.all((xform != 0) == mask)
 
@@ -402,7 +467,7 @@ def test_transform():
         lines    = readlines(testfile)
         xform    = np.genfromtxt(lines[:4])
         expected = np.genfromtxt(lines[ 4:])
-        result   = transform.transform(testcoords, xform)
+        result   = affine.transform(testcoords, xform)
 
         assert np.all(np.isclose(expected, result))
 
@@ -412,7 +477,7 @@ def test_transform():
         for axes in allAxes:
             atestcoords = testcoords[:, axes]
             aexpected   = expected[  :, axes]
-            aresult     = transform.transform(atestcoords, xform, axes=axes)
+            aresult     = affine.transform(atestcoords, xform, axes=axes)
 
             assert np.all(np.isclose(aexpected, aresult))
 
@@ -423,26 +488,26 @@ def test_transform():
     coords    = badcoords[:, :3]
 
     with pytest.raises(IndexError):
-        transform.transform(coords, badxform)
+        affine.transform(coords, badxform)
 
     with pytest.raises(ValueError):
-        transform.transform(badcoords, xform)
+        affine.transform(badcoords, xform)
 
     with pytest.raises(ValueError):
-        transform.transform(badcoords.reshape(5, 2, 4), xform)
+        affine.transform(badcoords.reshape(5, 2, 4), xform)
 
     with pytest.raises(ValueError):
-        transform.transform(badcoords.reshape(5, 2, 4), xform, axes=1)
+        affine.transform(badcoords.reshape(5, 2, 4), xform, axes=1)
 
     with pytest.raises(ValueError):
-        transform.transform(badcoords[:, (1, 2, 3)], xform, axes=[1, 2])
+        affine.transform(badcoords[:, (1, 2, 3)], xform, axes=[1, 2])
 
 
 def test_transform_vector(seed):
 
     # Some transform with a
     # translation component
-    xform = transform.compose([1, 2, 3],
+    xform = affine.compose([1, 2, 3],
                               [5, 10, 15],
                               [np.pi / 2, np.pi / 2, 0])
 
@@ -453,9 +518,9 @@ def test_transform_vector(seed):
         vecExpected = np.dot(xform, list(v) + [0])[:3]
         ptExpected  = np.dot(xform, list(v) + [1])[:3]
 
-        vecResult   = transform.transform(v, xform,         vector=True)
-        vec33Result = transform.transform(v, xform[:3, :3], vector=True)
-        ptResult    = transform.transform(v, xform,         vector=False)
+        vecResult   = affine.transform(v, xform,         vector=True)
+        vec33Result = affine.transform(v, xform[:3, :3], vector=True)
+        ptResult    = affine.transform(v, xform,         vector=False)
 
         assert np.all(np.isclose(vecExpected, vecResult))
         assert np.all(np.isclose(vecExpected, vec33Result))
@@ -478,83 +543,33 @@ def test_transformNormal(seed):
         rotations = -np.pi + np.random.random(3) * 2 * np.pi
         origin    = -100   + np.random.random(3) * 200
 
-        xform = transform.compose(scales,
+        xform = affine.compose(scales,
                                   offsets,
                                   rotations,
                                   origin)
 
         expected = tn(n, xform)
-        result   = transform.transformNormal(n, xform)
+        result   = affine.transformNormal(n, xform)
 
         assert np.all(np.isclose(expected, result))
 
 
-def test_flirtMatrixToSform():
-
-    testfile = op.join(datadir, 'test_transform_test_flirtMatrixToSform.txt')
-    lines    = readlines(testfile)
-    ntests   = len(lines) // 18
-
-    for i in range(ntests):
-        tlines    = lines[i * 18: i * 18 + 18]
-        srcShape  = [int(  w) for w in tlines[0].split()]
-        srcXform  = np.genfromtxt(tlines[1:5])
-        refShape  = [int(  w) for w in tlines[5].split()]
-        refXform  = np.genfromtxt(tlines[6:10])
-        flirtMat  = np.genfromtxt(tlines[10:14])
-        expected  = np.genfromtxt(tlines[14:18])
-
-        srcImg = fslimage.Image(np.zeros(srcShape), xform=srcXform)
-        refImg = fslimage.Image(np.zeros(refShape), xform=refXform)
-
-        result = transform.flirtMatrixToSform(flirtMat, srcImg, refImg)
-
-        assert np.all(np.isclose(result, expected))
-
-
-def test_sformToFlirtMatrix():
-    testfile = op.join(datadir, 'test_transform_test_flirtMatrixToSform.txt')
-    lines    = readlines(testfile)
-    ntests   = len(lines) // 18
-
-    for i in range(ntests):
-        tlines      = lines[i * 18: i * 18 + 18]
-        srcShape    = [int(  w) for w in tlines[0].split()]
-        srcXform    = np.genfromtxt(tlines[1:5])
-        refShape    = [int(  w) for w in tlines[5].split()]
-        refXform    = np.genfromtxt(tlines[6:10])
-        expected    = np.genfromtxt(tlines[10:14])
-        srcXformOvr = np.genfromtxt(tlines[14:18])
-
-        srcImg1 = fslimage.Image(np.zeros(srcShape), xform=srcXform)
-        srcImg2 = fslimage.Image(np.zeros(srcShape), xform=srcXform)
-        refImg  = fslimage.Image(np.zeros(refShape), xform=refXform)
-
-        srcImg2.voxToWorldMat = srcXformOvr
-
-        result1 = transform.sformToFlirtMatrix(srcImg1, refImg, srcXformOvr)
-        result2 = transform.sformToFlirtMatrix(srcImg2, refImg)
-
-        assert np.all(np.isclose(result1, expected))
-        assert np.all(np.isclose(result2, expected))
-
-
 def test_rmsdev():
 
     t1 = np.eye(4)
-    t2 = transform.scaleOffsetXform([1, 1, 1], [2, 0, 0])
+    t2 = affine.scaleOffsetXform([1, 1, 1], [2, 0, 0])
 
-    assert np.isclose(transform.rmsdev(t1, t2), 2)
-    assert np.isclose(transform.rmsdev(t1, t2, R=2), 2)
-    assert np.isclose(transform.rmsdev(t1, t2, R=2, xc=(1, 1, 1)), 2)
+    assert np.isclose(affine.rmsdev(t1, t2), 2)
+    assert np.isclose(affine.rmsdev(t1, t2, R=2), 2)
+    assert np.isclose(affine.rmsdev(t1, t2, R=2, xc=(1, 1, 1)), 2)
 
     t1       = np.eye(3)
     lastdist = 0
 
     for i in range(1, 11):
         rot    = np.pi * i / 10.0
-        t2     = transform.axisAnglesToRotMat(rot, 0, 0)
-        result = transform.rmsdev(t1, t2)
+        t2     = affine.axisAnglesToRotMat(rot, 0, 0)
+        result = affine.rmsdev(t1, t2)
 
         assert result > lastdist
 
@@ -562,9 +577,54 @@ def test_rmsdev():
 
     for i in range(11, 20):
         rot    = np.pi * i / 10.0
-        t2     = transform.axisAnglesToRotMat(rot, 0, 0)
-        result = transform.rmsdev(t1, t2)
+        t2     = affine.axisAnglesToRotMat(rot, 0, 0)
+        result = affine.rmsdev(t1, t2)
 
         assert result < lastdist
 
         lastdist = result
+
+
+def test_rescale():
+
+    with pytest.raises(ValueError):
+        affine.rescale((5, 5), (10, 10, 10))
+
+    assert np.all(affine.rescale((5, 5),       (5, 5))       == np.eye(3))
+    assert np.all(affine.rescale((5, 5, 5),    (5, 5, 5))    == np.eye(4))
+    assert np.all(affine.rescale((5, 5, 5, 5), (5, 5, 5, 5)) == np.eye(5))
+
+    # (old shape, new shape, origin, expect)
+    tests = [
+        ((5, 5), (10, 10), 'centre', np.array([[0.5, 0,    0],
+                                               [0,   0.5,  0],
+                                               [0,   0,    1]])),
+        ((5, 5), (10, 10), 'corner', np.array([[0.5, 0,   -0.25],
+                                               [0,   0.5, -0.25],
+                                               [0,   0,    1]])),
+        ((5, 5, 5), (10, 10, 10), 'centre', np.array([[0.5, 0,    0,   0],
+                                                      [0,   0.5,  0,   0],
+                                                      [0,   0,    0.5, 0],
+                                                      [0,   0,    0,   1]])),
+        ((5, 5, 5), (10, 10, 10), 'corner', np.array([[0.5, 0,    0,   -0.25],
+                                                      [0,   0.5,  0,   -0.25],
+                                                      [0,   0,    0.5, -0.25],
+                                                      [0,   0,    0,    1]])),
+        ((5, 5, 5, 5), (10, 10, 10, 10), 'centre', np.array([[0.5, 0,    0,   0,   0],
+                                                             [0,   0.5,  0,   0,   0],
+                                                             [0,   0,    0.5, 0,   0],
+                                                             [0,   0,    0,   0.5, 0],
+                                                             [0,   0,    0,   0,   1]])),
+        ((5, 5, 5, 5), (10, 10, 10, 10), 'corner', np.array([[0.5, 0,    0,   0,   -0.25],
+                                                             [0,   0.5,  0,   0,   -0.25],
+                                                             [0,   0,    0.5, 0,   -0.25],
+                                                             [0,   0,    0,   0.5, -0.25],
+                                                             [0,   0,    0,   0,   1]])),
+    ]
+
+    for oldshape, newshape, origin, expect in tests:
+
+        got = affine.rescale(oldshape, newshape, origin)
+        assert np.all(np.isclose(got, expect))
+        got = affine.rescale(newshape, oldshape, origin)
+        assert np.all(np.isclose(got, affine.invert(expect)))

fsl/tests/test_transform/test_flirt.py

+#!/usr/bin/env python
+#
+# test_transform_flirt.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+import itertools as it
+import os.path as op
+
+import numpy as np
+
+import fsl.data.image       as fslimage
+import fsl.transform.flirt  as flirt
+import fsl.transform.affine as affine
+import fsl.utils.tempdir    as tempdir
+
+from fsl.tests.test_transform.test_affine import readlines
+
+from fsl.tests import make_random_image
+
+
+datadir = op.join(op.dirname(__file__), 'testdata')
+
+
+def test_read_write():
+    with tempdir.tempdir():
+        aff = np.random.random((4, 4))
+        flirt.writeFlirt(aff, 'aff.mat')
+        got = flirt.readFlirt('aff.mat')
+        assert np.all(np.isclose(aff, got))
+
+
+def test_fromFlirt():
+
+    src = affine.compose(np.random.randint( 1,  5,  3),
+                         np.random.randint(-20, 20, 3),
+                         np.random.random(3) - 0.5)
+    ref = affine.compose(np.random.randint( 1,  5,  3),
+                         np.random.randint(-20, 20, 3),
+                         np.random.random(3) - 0.5)
+
+    src      = fslimage.Image(make_random_image(xform=src))
+    ref      = fslimage.Image(make_random_image(xform=ref))
+    flirtmat = affine.concat(ref.getAffine('world', 'fsl'),
+                             src.getAffine('fsl',   'world'))
+
+
+    spaces = it.permutations(('voxel', 'fsl', 'world'), 2)
+
+    for from_, to in spaces:
+        got = flirt.fromFlirt(flirtmat, src, ref, from_, to)
+        exp = affine.concat(ref.getAffine('fsl', to),
+                            flirtmat,
+                            src.getAffine(from_, 'fsl'))
+
+        assert np.all(np.isclose(got, exp))
+
+
+def test_toFlirt():
+
+    src = affine.compose(np.random.randint( 1,  5,  3),
+                         np.random.randint(-20, 20, 3),
+                         np.random.random(3) - 0.5)
+    ref = affine.compose(np.random.randint( 1,  5,  3),
+                         np.random.randint(-20, 20, 3),
+                         np.random.random(3) - 0.5)
+
+    src      = fslimage.Image(make_random_image(xform=src))
+    ref      = fslimage.Image(make_random_image(xform=ref))
+    flirtmat = affine.concat(ref.getAffine('world', 'fsl'),
+                             src.getAffine('fsl',   'world'))
+
+
+    spaces = it.permutations(('voxel', 'fsl', 'world'), 2)
+
+    for from_, to in spaces:
+        xform = affine.concat(ref.getAffine('world', to),
+                              src.getAffine(from_, 'world'))
+        got = flirt.toFlirt(xform, src, ref, from_, to)
+
+        assert np.all(np.isclose(got, flirtmat))
+
+
+def test_flirtMatrixToSform():
+
+    testfile = op.join(datadir, 'test_transform_test_flirtMatrixToSform.txt')
+    lines    = readlines(testfile)
+    ntests   = len(lines) // 18
+
+    for i in range(ntests):
+        tlines    = lines[i * 18: i * 18 + 18]
+        srcShape  = [int(  w) for w in tlines[0].split()]
+        srcXform  = np.genfromtxt(tlines[1:5])
+        refShape  = [int(  w) for w in tlines[5].split()]
+        refXform  = np.genfromtxt(tlines[6:10])
+        flirtMat  = np.genfromtxt(tlines[10:14])
+        expected  = np.genfromtxt(tlines[14:18])
+
+        srcImg = fslimage.Image(np.zeros(srcShape), xform=srcXform)
+        refImg = fslimage.Image(np.zeros(refShape), xform=refXform)
+
+        result = flirt.flirtMatrixToSform(flirtMat, srcImg, refImg)
+
+        assert np.all(np.isclose(result, expected))
+
+
+def test_sformToFlirtMatrix():
+    testfile = op.join(datadir, 'test_transform_test_flirtMatrixToSform.txt')
+    lines    = readlines(testfile)
+    ntests   = len(lines) // 18
+
+    for i in range(ntests):
+        tlines      = lines[i * 18: i * 18 + 18]
+        srcShape    = [int(  w) for w in tlines[0].split()]
+        srcXform    = np.genfromtxt(tlines[1:5])
+        refShape    = [int(  w) for w in tlines[5].split()]
+        refXform    = np.genfromtxt(tlines[6:10])
+        expected    = np.genfromtxt(tlines[10:14])
+        srcXformOvr = np.genfromtxt(tlines[14:18])
+
+        srcImg1 = fslimage.Image(np.zeros(srcShape), xform=srcXform)
+        srcImg2 = fslimage.Image(np.zeros(srcShape), xform=srcXform)
+        refImg  = fslimage.Image(np.zeros(refShape), xform=refXform)
+
+        srcImg2.voxToWorldMat = srcXformOvr
+
+        result1 = flirt.sformToFlirtMatrix(srcImg1, refImg, srcXformOvr)
+        result2 = flirt.sformToFlirtMatrix(srcImg2, refImg)
+
+        assert np.all(np.isclose(result1, expected))
+        assert np.all(np.isclose(result2, expected))

fsl/tests/test_transform/test_fnirt.py

+#!/usr/bin/env python
+#
+# test_fnirt.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+import os.path as op
+import itertools as it
+
+import numpy as np
+import nibabel as nib
+
+import pytest
+
+import fsl.data.image          as fslimage
+import fsl.utils.tempdir       as tempdir
+import fsl.data.constants      as constants
+import fsl.transform.affine    as affine
+import fsl.transform.nonlinear as nonlinear
+import fsl.transform.fnirt     as fnirt
+
+from fsl.tests.test_transform.test_nonlinear import _random_affine_field
+
+
+datadir = op.join(op.dirname(__file__), 'testdata', 'nonlinear')
+
+
+def test_readFnirt():
+
+    src  = op.join(datadir, 'src')
+    ref  = op.join(datadir, 'ref')
+    coef = op.join(datadir, 'coefficientfield')
+    disp = op.join(datadir, 'displacementfield')
+
+    src  = fslimage.Image(src)
+    ref  = fslimage.Image(ref)
+    coef = fnirt.readFnirt(coef, src, ref)
+    disp = fnirt.readFnirt(disp, src, ref)
+
+    with pytest.raises(ValueError):
+        fnirt.readFnirt(src, src, ref)
+
+    assert isinstance(coef, nonlinear.CoefficientField)
+    assert isinstance(disp, nonlinear.DeformationField)
+
+    assert coef.src.sameSpace(src)
+    assert coef.ref.sameSpace(ref)
+    assert disp.src.sameSpace(src)
+    assert disp.ref.sameSpace(ref)
+    assert coef.srcSpace == 'fsl'
+    assert coef.refSpace == 'fsl'
+    assert disp.srcSpace == 'fsl'
+    assert disp.refSpace == 'fsl'
+
+
+def test_readFnirt_defType_intent():
+    src  = op.join(datadir, 'src.nii.gz')
+    ref  = op.join(datadir, 'ref.nii.gz')
+    coef = op.join(datadir, 'coefficientfield.nii.gz')
+    disp = op.join(datadir, 'displacementfield.nii.gz')
+
+    src  = fslimage.Image(src)
+    ref  = fslimage.Image(ref)
+
+    field = fnirt.readFnirt(disp, src, ref, defType='absolute')
+    assert field.deformationType == 'absolute'
+    field = fnirt.readFnirt(disp, src, ref, defType='relative')
+    assert field.deformationType == 'relative'
+
+    img = nib.load(coef)
+    img.header['intent_code'] = 0
+    with tempdir.tempdir():
+        img.to_filename('field.nii.gz')
+
+        with pytest.raises(ValueError):
+            fnirt.readFnirt('field', src, ref)
+
+        field = fnirt.readFnirt(
+            'field', src, ref,
+            intent=constants.FSL_CUBIC_SPLINE_COEFFICIENTS)
+        assert isinstance(field, nonlinear.CoefficientField)
+
+        field = fnirt.readFnirt(
+            'field', src, ref,
+            intent=constants.FSL_FNIRT_DISPLACEMENT_FIELD)
+        assert isinstance(field, nonlinear.DeformationField)
+
+
+def test_toFnirt():
+
+    def check(got, exp):
+        tol = dict(atol=1e-5, rtol=1e-5)
+        assert np.all(np.isclose(got.data, exp.data, **tol))
+        assert got.src.sameSpace(exp.src)
+        assert got.ref.sameSpace(exp.ref)
+        assert got.srcSpace == 'fsl'
+        assert got.refSpace == 'fsl'
+
+    basefield, xform = _random_affine_field()
+    src = basefield.src
+    ref = basefield.ref
+
+    spaces = it.permutations(('voxel', 'fsl', 'world'), 2)
+
+    for from_, to in spaces:
+        field = nonlinear.convertDeformationSpace(basefield, from_, to)
+        got = fnirt.toFnirt(field)
+        check(got, basefield)
+
+    src  = fslimage.Image(op.join(datadir, 'src'))
+    ref  = fslimage.Image(op.join(datadir, 'ref'))
+    coef = fnirt.readFnirt(op.join(datadir, 'coefficientfield'),  src, ref)
+    got  = fnirt.toFnirt(coef)
+    check(got, coef)
+
+
+def test_fromFnirt():
+
+    basefield, basexform = _random_affine_field()
+    src = basefield.src
+    ref = basefield.ref
+    spaces = list(it.permutations(('voxel', 'fsl', 'world'), 2))
+
+    for from_, to in spaces:
+
+        got = fnirt.fromFnirt(basefield, from_, to)
+
+        assert got.srcSpace == to
+        assert got.refSpace == from_
+
+        coords = [np.random.randint(0, basefield.shape[0], 5),
+                  np.random.randint(0, basefield.shape[1], 5),
+                  np.random.randint(0, basefield.shape[2], 5)]
+        coords = np.array(coords).T
+
+        coords = affine.transform(coords, ref.getAffine('voxel', from_))
+
+        aff = affine.concat(src.getAffine('fsl', to),
+                            basexform,
+                            ref.getAffine(from_, 'fsl'))
+
+        got = got.transform(coords)
+        exp = affine.transform(coords, aff)
+
+        enan = np.isnan(exp)
+        gnan = np.isnan(got)
+
+        assert np.all(np.isclose(enan, gnan))
+        assert np.all(np.isclose(exp[~enan], got[~gnan]))
+
+    # Converting from a FNIRT coefficient field
+    src  = fslimage.Image(op.join(datadir, 'src'))
+    ref  = fslimage.Image(op.join(datadir, 'ref'))
+    coef = fnirt.readFnirt(op.join(datadir, 'coefficientfield'),  src, ref)
+    disp = fnirt.readFnirt(op.join(datadir, 'displacementfield'), src, ref)
+
+    for from_, to in spaces:
+
+        cnv = fnirt.fromFnirt(coef, from_, to)
+        exp = nonlinear.convertDeformationSpace(disp, from_, to)
+        tol = dict(atol=1e-5, rtol=1e-5)
+        assert np.all(np.isclose(cnv.data, exp.data, **tol))

fsl/tests/test_transform/test_nonlinear.py

+#!/usr/bin/env python
+
+import itertools as it
+import os.path   as op
+
+import numpy   as np
+
+import fsl.data.image           as fslimage
+import fsl.utils.image.resample as resample
+import fsl.utils.image.roi      as roi
+import fsl.transform.affine     as affine
+import fsl.transform.nonlinear  as nonlinear
+import fsl.transform.fnirt      as fnirt
+
+
+datadir = op.join(op.dirname(__file__), 'testdata')
+
+
+def _random_image():
+    vx, vy, vz = np.random.randint(10, 50, 3)
+    dx, dy, dz = np.random.randint( 1, 10, 3)
+    data       = (np.random.random((vx, vy, vz)) - 0.5) * 10
+    aff        = affine.compose(
+        (dx, dy, dz),
+        np.random.randint(1, 100, 3),
+        np.random.random(3) * np.pi / 2)
+
+    return fslimage.Image(data, xform=aff)
+
+
+def _random_field():
+
+    src        = _random_image()
+    vx, vy, vz = np.random.randint(10, 50, 3)
+    dx, dy, dz = np.random.randint( 1, 10, 3)
+
+    field = (np.random.random((vx, vy, vz, 3)) - 0.5) * 10
+    aff   = affine.compose(
+        (dx, dy, dz),
+        np.random.randint(1, 100, 3),
+        np.random.random(3) * np.pi / 2)
+
+    return nonlinear.DeformationField(field, src=src, xform=aff)
+
+
+def _affine_field(src, ref, xform, srcSpace, refSpace, shape=None, fv2w=None):
+
+    if shape is None: shape = ref.shape[:3]
+    if fv2w  is None: fv2w  = ref.getAffine('voxel', 'world')
+
+    rx, ry, rz = np.meshgrid(np.arange(shape[0]),
+                             np.arange(shape[1]),
+                             np.arange(shape[2]), indexing='ij')
+
+    rvoxels  = np.vstack((rx.flatten(), ry.flatten(), rz.flatten())).T
+    f2r      = affine.concat(ref.getAffine('world', refSpace), fv2w)
+    rcoords  = affine.transform(rvoxels, f2r)
+    scoords  = affine.transform(rcoords, xform)
+
+    field    = np.zeros(list(shape[:3]) + [3])
+    field[:] = (scoords - rcoords).reshape(*it.chain(shape, [3]))
+    field    = nonlinear.DeformationField(field, src, ref,
+                                          srcSpace=srcSpace,
+                                          refSpace=refSpace,
+                                          xform=fv2w,
+                                          header=ref.header,
+                                          defType='relative')
+    return field
+
+
+def _random_affine_field():
+
+    src = _random_image()
+    ref = _random_image()
+
+    # our test field just encodes an affine
+    xform = affine.compose(
+        np.random.randint(2, 5, 3),
+        np.random.randint(1, 10, 3),
+        np.random.random(3))
+
+    rx, ry, rz = np.meshgrid(np.arange(ref.shape[0]),
+                             np.arange(ref.shape[1]),
+                             np.arange(ref.shape[2]), indexing='ij')
+
+    rvoxels  = np.vstack((rx.flatten(), ry.flatten(), rz.flatten())).T
+    rcoords  = affine.transform(rvoxels, ref.getAffine('voxel', 'fsl'))
+    scoords  = affine.transform(rcoords, xform)
+
+    field    = np.zeros(list(ref.shape[:3]) + [3])
+    field[:] = (scoords - rcoords).reshape(*it.chain(ref.shape, [3]))
+    field    = nonlinear.DeformationField(field, src, ref,
+                                          header=ref.header,
+                                          defType='relative')
+    return field, xform
+
+
+def _field_coords(field):
+    vx, vy, vz = field.shape[ :3]
+    coords     = np.meshgrid(np.arange(vx),
+                             np.arange(vy),
+                             np.arange(vz), indexing='ij')
+    coords = np.array(coords).transpose((1, 2, 3, 0))
+    return affine.transform(
+        coords.reshape(-1, 3),
+        field.getAffine('voxel', 'fsl')).reshape(field.shape)
+
+
+def test_detectDeformationType():
+    relfield = _random_field()
+    coords   = _field_coords(relfield)
+    absfield = nonlinear.DeformationField(
+        relfield.data + coords,
+        src=relfield.src,
+        xform=relfield.voxToWorldMat)
+
+    assert nonlinear.detectDeformationType(relfield) == 'relative'
+    assert nonlinear.detectDeformationType(absfield) == 'absolute'
+
+
+def test_convertDeformationType():
+
+    relfield = _random_field()
+    coords   = _field_coords(relfield)
+    absfield = nonlinear.DeformationField(
+        relfield.data + coords,
+        src=relfield.src,
+        xform=relfield.voxToWorldMat)
+
+    gotconvrel1 = nonlinear.convertDeformationType(relfield)
+    gotconvabs1 = nonlinear.convertDeformationType(absfield)
+    gotconvrel2 = nonlinear.convertDeformationType(relfield, 'absolute')
+    gotconvabs2 = nonlinear.convertDeformationType(absfield, 'relative')
+
+    tol = dict(atol=1e-3, rtol=1e-3)
+
+    assert np.all(np.isclose(gotconvrel1, absfield.data, **tol))
+    assert np.all(np.isclose(gotconvabs1, relfield.data, **tol))
+    assert np.all(np.isclose(gotconvrel2, absfield.data, **tol))
+    assert np.all(np.isclose(gotconvabs2, relfield.data, **tol))
+
+
+def test_convertDeformationSpace():
+
+    basefield, xform = _random_affine_field()
+    src              = basefield.src
+    ref              = basefield.ref
+
+    # generate reference fsl->fsl coordinate mappings
+
+    # For each combination of srcspace->tospace
+    # Generate random coordinates, check that
+    # displacements are correct
+    spaces = ['fsl', 'voxel', 'world']
+    spaces = list(it.combinations_with_replacement(spaces, 2))
+    spaces = spaces + [(r, s) for s, r in spaces]
+    spaces = list(set(spaces))
+
+    for from_, to in spaces:
+
+        refcoords = [np.random.randint(0, basefield.shape[0], 5),
+                     np.random.randint(0, basefield.shape[1], 5),
+                     np.random.randint(0, basefield.shape[2], 5)]
+        refcoords = np.array(refcoords, dtype=int).T
+        refcoords = affine.transform(refcoords, ref.getAffine('voxel', 'fsl'))
+        srccoords = basefield.transform(refcoords)
+
+        field   = nonlinear.convertDeformationSpace(basefield, from_, to)
+        premat  = ref.getAffine('fsl', from_)
+        postmat = src.getAffine('fsl', to)
+
+        input  = affine.transform(refcoords, premat)
+        expect = affine.transform(srccoords, postmat)
+
+        got  = field.transform(input)
+        enan = np.isnan(expect)
+        gnan = np.isnan(got)
+
+        assert np.all(np.isclose(enan, gnan))
+        assert np.all(np.isclose(expect[~enan], got[~gnan]))
+
+
+def test_DeformationField_transform():
+
+    relfield, xform = _random_affine_field()
+    src             = relfield.src
+    ref             = relfield.ref
+
+    rx, ry, rz = np.meshgrid(np.arange(ref.shape[0]),
+                             np.arange(ref.shape[1]),
+                             np.arange(ref.shape[2]), indexing='ij')
+    rvoxels  = np.vstack((rx.flatten(), ry.flatten(), rz.flatten())).T
+    rcoords  = affine.transform(rvoxels, ref.getAffine('voxel', 'fsl'))
+    scoords  = affine.transform(rcoords, xform)
+    svoxels  = affine.transform(scoords, src.getAffine('fsl', 'voxel'))
+
+    absfield    = np.zeros(list(ref.shape[:3]) + [3])
+    absfield[:] = scoords.reshape(*it.chain(ref.shape, [3]))
+    absfield    = nonlinear.DeformationField(absfield, src, ref,
+                                             header=ref.header,
+                                             defType='absolute')
+
+    got = relfield.transform(rcoords)
+    assert np.all(np.isclose(got, scoords))
+    got = absfield.transform(rcoords)
+    assert np.all(np.isclose(got, scoords))
+
+    # test single set of coords
+    got = absfield.transform(rcoords[0])
+    assert np.all(np.isclose(got, scoords[0]))
+
+    got = relfield.transform(rvoxels, from_='voxel', to='voxel')
+    assert np.all(np.isclose(got, svoxels))
+    got = absfield.transform(rvoxels, from_='voxel', to='voxel')
+    assert np.all(np.isclose(got, svoxels))
+
+    # test out of bounds are returned as nan
+    rvoxels = np.array([[-1, -1, -1],
+                        [ 0,  0,  0]])
+    rcoords  = affine.transform(rvoxels, ref.getAffine('voxel', 'fsl'))
+    scoords  = affine.transform(rcoords, xform)
+    svoxels  = affine.transform(scoords, src.getAffine('fsl', 'voxel'))
+
+    got = relfield.transform(rcoords)
+    assert np.all(np.isnan(got[0, :]))
+    assert np.all(np.isclose(got[1, :], scoords[1, :]))
+    got = absfield.transform(rcoords)
+    assert np.all(np.isnan(got[0, :]))
+    assert np.all(np.isclose(got[1, :], scoords[1, :]))
+
+
+def test_CoefficientField_displacements():
+
+    nldir = op.join(datadir, 'nonlinear')
+    src   = op.join(nldir, 'src.nii.gz')
+    ref   = op.join(nldir, 'ref.nii.gz')
+    cf    = op.join(nldir, 'coefficientfield.nii.gz')
+    df    = op.join(nldir, 'displacementfield_no_premat.nii.gz')
+
+    src = fslimage.Image(src)
+    ref = fslimage.Image(ref)
+    cf  = fnirt.readFnirt(cf, src, ref)
+    df  = fnirt.readFnirt(df, src, ref)
+
+    ix, iy, iz = ref.shape[:3]
+    x,  y,  z  = np.meshgrid(np.arange(ix),
+                             np.arange(iy),
+                             np.arange(iz), indexing='ij')
+    x          = x.flatten()
+    y          = y.flatten()
+    z          = z.flatten()
+    xyz        = np.vstack((x, y, z)).T
+
+    disps = cf.displacements(xyz)
+    disps = disps.reshape(df.shape)
+
+    tol = dict(atol=1e-5, rtol=1e-5)
+    assert np.all(np.isclose(disps, df.data, **tol))
+
+
+def test_CoefficientField_transform():
+    nldir = op.join(datadir, 'nonlinear')
+    src   = op.join(nldir, 'src.nii.gz')
+    ref   = op.join(nldir, 'ref.nii.gz')
+    cf    = op.join(nldir, 'coefficientfield.nii.gz')
+    df    = op.join(nldir, 'displacementfield.nii.gz')
+    dfnp  = op.join(nldir, 'displacementfield_no_premat.nii.gz')
+
+    src  = fslimage.Image(src)
+    ref  = fslimage.Image(ref)
+    cf   = fnirt.readFnirt(cf,   src, ref)
+    df   = fnirt.readFnirt(df,   src, ref)
+    dfnp = fnirt.readFnirt(dfnp, src, ref)
+
+    spaces = ['fsl', 'voxel', 'world']
+    spaces = list(it.combinations_with_replacement(spaces, 2))
+    spaces = spaces + [(r, s) for s, r in spaces]
+    spaces = list(set(spaces))
+
+    rx, ry, rz = np.meshgrid(np.arange(ref.shape[0]),
+                             np.arange(ref.shape[1]),
+                             np.arange(ref.shape[2]), indexing='ij')
+    rvoxels  = np.vstack((rx.flatten(), ry.flatten(), rz.flatten())).T
+
+    refcoords = {
+        'voxel' : rvoxels,
+        'fsl'   : affine.transform(rvoxels, ref.getAffine('voxel', 'fsl')),
+        'world' : affine.transform(rvoxels, ref.getAffine('voxel', 'world'))
+    }
+
+    srccoords = refcoords['fsl'] + df.data.reshape(-1, 3)
+    srccoords = {
+        'voxel' : affine.transform(srccoords, src.getAffine('fsl', 'voxel')),
+        'fsl'   : srccoords,
+        'world' : affine.transform(srccoords, src.getAffine('fsl', 'world'))
+    }
+
+    srccoordsnp = refcoords['fsl'] + dfnp.data.reshape(-1, 3)
+    srccoordsnp = {
+        'voxel' : affine.transform(srccoordsnp, src.getAffine('fsl', 'voxel')),
+        'fsl'   : srccoordsnp,
+        'world' : affine.transform(srccoordsnp, src.getAffine('fsl', 'world'))
+    }
+
+    tol = dict(atol=1e-5, rtol=1e-5)
+    for srcspace, refspace in spaces:
+        got   = cf.transform(refcoords[refspace], refspace, srcspace)
+        gotnp = cf.transform(refcoords[refspace], refspace, srcspace,
+                             premat=False)
+        assert np.all(np.isclose(got,   srccoords[  srcspace], **tol))
+        assert np.all(np.isclose(gotnp, srccoordsnp[srcspace], **tol))
+
+
+def test_coefficientField_transform_altref():
+
+    # test coordinates (manually determined).
+    # original ref image is 2mm isotropic,
+    # resampled is 1mm. Each tuple contains:
+    #
+    # (src, ref2mm, ref1mm)
+    coords = [
+        ((18.414, 26.579, 25.599), (11, 19, 11), (22, 38, 22)),
+        ((14.727, 22.480, 20.340), ( 8, 17,  8), (16, 34, 16)),
+        ((19.932, 75.616, 27.747), (11, 45,  5), (22, 90, 10))
+    ]
+
+    nldir  = op.join(datadir, 'nonlinear')
+    src    = op.join(nldir, 'src.nii.gz')
+    ref    = op.join(nldir, 'ref.nii.gz')
+    cf     = op.join(nldir, 'coefficientfield.nii.gz')
+
+    src      = fslimage.Image(src)
+    ref2mm   = fslimage.Image(ref)
+    ref1mm   = ref2mm.adjust((1, 1, 1))
+    cfref2mm = fnirt.readFnirt(cf, src, ref2mm)
+    cfref1mm = fnirt.readFnirt(cf, src, ref1mm)
+
+    for srcc, ref2mmc, ref1mmc in coords:
+        ref2mmc = cfref2mm.transform(ref2mmc, 'voxel', 'voxel')
+        ref1mmc = cfref1mm.transform(ref1mmc, 'voxel', 'voxel')
+
+        assert np.all(np.isclose(ref2mmc, srcc, 1e-4))
+        assert np.all(np.isclose(ref1mmc, srcc, 1e-4))
+
+
+def test_coefficientFieldToDeformationField():
+
+    nldir = op.join(datadir, 'nonlinear')
+    src   = op.join(nldir, 'src.nii.gz')
+    ref   = op.join(nldir, 'ref.nii.gz')
+    cf    = op.join(nldir, 'coefficientfield.nii.gz')
+    df    = op.join(nldir, 'displacementfield.nii.gz')
+    dfnp  = op.join(nldir, 'displacementfield_no_premat.nii.gz')
+
+    src   = fslimage.Image(src)
+    ref   = fslimage.Image(ref)
+    cf    = fnirt.readFnirt(cf,   src, ref)
+    rdf   = fnirt.readFnirt(df,   src, ref)
+    rdfnp = fnirt.readFnirt(dfnp, src, ref)
+    adf   = nonlinear.convertDeformationType(rdf)
+    adfnp = nonlinear.convertDeformationType(rdfnp)
+
+    rcnv   = nonlinear.coefficientFieldToDeformationField(cf)
+    acnv   = nonlinear.coefficientFieldToDeformationField(cf,
+                                                          defType='absolute')
+    acnvnp = nonlinear.coefficientFieldToDeformationField(cf,
+                                                          defType='absolute',
+                                                          premat=False)
+    rcnvnp = nonlinear.coefficientFieldToDeformationField(cf,
+                                                          premat=False)
+
+    tol = dict(atol=1e-5, rtol=1e-5)
+    assert np.all(np.isclose(rcnv.data,   rdf  .data, **tol))
+    assert np.all(np.isclose(acnv.data,   adf  .data, **tol))
+    assert np.all(np.isclose(rcnvnp.data, rdfnp.data, **tol))
+    assert np.all(np.isclose(acnvnp.data, adfnp.data, **tol))
+
+
+def test_applyDeformation():
+
+    src2ref = affine.compose(
+        np.random.randint(2, 5, 3),
+        np.random.randint(1, 10, 3),
+        np.random.random(3))
+    ref2src = affine.invert(src2ref)
+
+    srcdata = np.random.randint(1, 65536, (10, 10, 10), dtype=np.int32)
+    refdata = np.random.randint(1, 65536, (10, 10, 10), dtype=np.int32)
+
+    src   = fslimage.Image(srcdata)
+    ref   = fslimage.Image(refdata, xform=src2ref)
+    field = _affine_field(src, ref, ref2src, 'world', 'world')
+
+    expect, xf = resample.resampleToReference(
+        src, ref, matrix=src2ref, order=1, mode='nearest')
+    result = nonlinear.applyDeformation(
+        src, field, order=1, mode='nearest')
+
+    assert np.all(np.isclose(expect, result))
+
+
+def test_applyDeformation_altsrc():
+
+    src2ref = affine.compose(
+        np.random.randint(2, 5, 3),
+        np.random.randint(1, 10, 3),
+        [0, 0, 0])
+    ref2src = affine.invert(src2ref)
+
+    srcdata = np.random.randint(1, 65536, (10, 10, 10), dtype=np.int32)
+    refdata = np.random.randint(1, 65536, (10, 10, 10), dtype=np.int32)
+
+    src   = fslimage.Image(srcdata)
+    ref   = fslimage.Image(refdata, xform=src2ref)
+    field = _affine_field(src, ref, ref2src, 'world', 'world')
+
+    # First try a down-sampled version
+    # of the original source image
+    altsrc, xf = resample.resample(src, (5, 5, 5), origin='corner')
+    altsrc     = fslimage.Image(altsrc, xform=xf, header=src.header)
+    expect, xf = resample.resampleToReference(
+        altsrc, ref, matrix=src2ref, order=1, mode='nearest')
+    result = nonlinear.applyDeformation(
+        altsrc, field, order=1, mode='nearest')
+    assert np.all(np.isclose(expect, result))
+
+    # Now try a down-sampled ROI
+    # of the original source image
+    altsrc     = roi.roi(src, [(2, 9), (2, 9), (2, 9)])
+    altsrc, xf = resample.resample(altsrc, (4, 4, 4))
+    altsrc     = fslimage.Image(altsrc, xform=xf, header=src.header)
+    expect, xf = resample.resampleToReference(
+        altsrc, ref, matrix=src2ref, order=1, mode='nearest')
+    result = nonlinear.applyDeformation(
+        altsrc, field, order=1, mode='nearest')
+    assert np.all(np.isclose(expect, result))
+
+    # down-sampled and offset ROI
+    # of the original source image
+    altsrc     = roi.roi(src, [(-5, 8), (-5, 8), (-5, 8)])
+    altsrc, xf = resample.resample(altsrc, (6, 6, 6))
+    altsrc     = fslimage.Image(altsrc, xform=xf, header=src.header)
+    expect, xf = resample.resampleToReference(
+        altsrc, ref, matrix=src2ref, order=1, mode='nearest')
+    result = nonlinear.applyDeformation(
+        altsrc, field, order=1, mode='nearest')
+    assert np.all(np.isclose(expect, result))
+
+
+def test_applyDeformation_premat():
+
+    src2ref = affine.compose(
+        np.random.randint(2, 5, 3),
+        np.random.randint(1, 10, 3),
+        [0, 0, 0])
+    ref2src = affine.invert(src2ref)
+
+    srcdata = np.random.randint(1, 65536, (10, 10, 10), dtype=np.int32)
+    refdata = np.random.randint(1, 65536, (10, 10, 10), dtype=np.int32)
+
+    src   = fslimage.Image(srcdata)
+    ref   = fslimage.Image(refdata, xform=src2ref)
+    field = _affine_field(src, ref, ref2src, 'world', 'world')
+
+    # First try a down-sampled version
+    # of the original source image
+    altsrc, xf = resample.resample(src, (5, 5, 5), origin='corner')
+    altsrc     = fslimage.Image(altsrc, xform=xf, header=src.header)
+    expect, xf = resample.resampleToReference(
+        altsrc, ref, matrix=src2ref, order=1, mode='nearest')
+    premat = affine.concat(src   .getAffine('world', 'voxel'),
+                           altsrc.getAffine('voxel', 'world'))
+    result = nonlinear.applyDeformation(
+        altsrc, field, order=1, mode='nearest', premat=premat)
+    assert np.all(np.isclose(expect, result))
+
+    # Now try a down-sampled ROI
+    # of the original source image
+    altsrc     = roi.roi(src, [(2, 9), (2, 9), (2, 9)])
+    altsrc, xf = resample.resample(altsrc, (4, 4, 4))
+    altsrc     = fslimage.Image(altsrc, xform=xf, header=src.header)
+    expect, xf = resample.resampleToReference(
+        altsrc, ref, matrix=src2ref, order=1, mode='nearest')
+    premat = affine.concat(src   .getAffine('world', 'voxel'),
+                           altsrc.getAffine('voxel', 'world'))
+    result = nonlinear.applyDeformation(
+        altsrc, field, order=1, mode='nearest', premat=premat)
+    assert np.all(np.isclose(expect, result))
+
+    # down-sampled and offset ROI
+    # of the original source image
+    altsrc     = roi.roi(src, [(-5, 8), (-5, 8), (-5, 8)])
+    altsrc, xf = resample.resample(altsrc, (6, 6, 6))
+    altsrc     = fslimage.Image(altsrc, xform=xf, header=src.header)
+    expect, xf = resample.resampleToReference(
+        altsrc, ref, matrix=src2ref, order=1, mode='nearest')
+    premat = affine.concat(src   .getAffine('world', 'voxel'),
+                           altsrc.getAffine('voxel', 'world'))
+    result = nonlinear.applyDeformation(
+        altsrc, field, order=1, mode='nearest', premat=premat)
+    assert np.all(np.isclose(expect, result))
+
+
+def test_applyDeformation_altref():
+    src2ref = affine.compose(
+        np.random.randint(2, 5, 3),
+        np.random.randint(1, 10, 3),
+        np.random.random(3))
+    ref2src = affine.invert(src2ref)
+
+    srcdata = np.random.randint(1, 65536, (10, 10, 10), dtype=np.int32)
+    refdata = np.random.randint(1, 65536, (10, 10, 10), dtype=np.int32)
+
+    src   = fslimage.Image(srcdata)
+    ref   = fslimage.Image(refdata, xform=src2ref)
+    field = _affine_field(src, ref, ref2src, 'world', 'world')
+
+    altrefxform = affine.concat(
+        src2ref,
+        affine.scaleOffsetXform([1, 1, 1], [5, 0, 0]))
+
+    altref = fslimage.Image(refdata, xform=altrefxform)
+
+    expect, xf = resample.resampleToReference(
+        src, altref, matrix=src2ref, order=1, mode='constant', cval=0)
+    result = nonlinear.applyDeformation(
+        src, field, ref=altref, order=1, mode='constant', cval=0)
+
+    # boundary voxels can get truncated
+    # (4 is the altref-ref overlap boundary)
+    expect[4, :, :] = 0
+    result[4, :, :] = 0
+    expect = expect[1:-1, 1:-1, 1:-1]
+    result = result[1:-1, 1:-1, 1:-1]
+
+    assert np.all(np.isclose(expect, result))
+
+
+# test when reference/field
+# are not voxel-aligned
+def test_applyDeformation_worldAligned():
+    refv2w   = affine.scaleOffsetXform([1, 1, 1], [10,   10,   10])
+    fieldv2w = affine.scaleOffsetXform([2, 2, 2], [10.5, 10.5, 10.5])
+    src2ref  = refv2w
+    ref2src  = affine.invert(src2ref)
+
+    srcdata = np.random.randint(1, 65536, (10, 10, 10), dtype=np.int32)
+
+    src   = fslimage.Image(srcdata)
+    ref   = fslimage.Image(srcdata, xform=src2ref)
+    field = _affine_field(src, ref, ref2src, 'world', 'world',
+                          shape=(5, 5, 5), fv2w=fieldv2w)
+
+    field = nonlinear.DeformationField(
+        nonlinear.convertDeformationType(field, 'absolute'),
+        header=field.header,
+        src=src,
+        ref=ref,
+        srcSpace='world',
+        refSpace='world',
+        defType='absolute')
+
+    expect, xf = resample.resampleToReference(
+        src, ref, matrix=src2ref, order=1, mode='constant', cval=0)
+    result = nonlinear.applyDeformation(
+        src, field, order=1, mode='constant', cval=0)
+
+    expect = expect[1:-1, 1:-1, 1:-1]
+    result = result[1:-1, 1:-1, 1:-1]
+
+    assert np.all(np.isclose(expect, result))

fsl/tests/test_transform/test_x5.py

+#!/usr/bin/env python
+#
+# test_x5.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+
+import os.path as op
+import numpy as np
+
+import pytest
+
+import h5py
+
+import fsl.data.image          as fslimage
+import fsl.utils.tempdir       as tempdir
+import fsl.transform.affine    as affine
+import fsl.transform.fnirt     as fnirt
+import fsl.transform.nonlinear as nonlinear
+import fsl.transform.x5        as x5
+
+from fsl.tests import make_random_image
+
+
+def _check_metadata(group):
+    assert group.attrs['Format']  == x5.X5_FORMAT
+    assert group.attrs['Version'] == x5.X5_VERSION
+
+
+def _check_affine(group, xform):
+    assert group.attrs['Type'] == 'affine'
+    gotxform = np.array(group['Matrix'])
+    assert np.all(np.isclose(gotxform, xform))
+
+
+def _check_space(group, img):
+    assert group.attrs['Type'] == 'image'
+    assert np.all(np.isclose(group.attrs['Size'],   img.shape[ :3]))
+    assert np.all(np.isclose(group.attrs['Scales'], img.pixdim[:3]))
+    _check_affine(group['Mapping'], img.voxToWorldMat)
+
+
+def _check_deformation(group, field):
+    assert group.attrs['Type'] == 'deformation'
+    assert group.attrs['SubType'] == field.deformationType
+    xform = np.array(group['Matrix'])
+    assert np.all(np.isclose(xform, field.data))
+    _check_affine(group['Mapping'], field.voxToWorldMat)
+
+
+def test_readWriteLinearX5():
+    with tempdir.tempdir():
+        make_random_image('src.nii')
+        make_random_image('ref.nii')
+        xform = affine.compose(
+            np.random.randint(1, 5, 3),
+            np.random.randint(-10, 10, 3),
+            -np.pi / 4 + np.random.random(3) * np.pi / 2)
+
+        src = fslimage.Image('src.nii')
+        ref = fslimage.Image('ref.nii')
+
+        x5.writeLinearX5('linear.x5', xform, src, ref)
+
+        gotxform, gotsrc, gotref = x5.readLinearX5('linear.x5')
+        assert np.all(np.isclose(gotxform, xform))
+        assert gotsrc.sameSpace(src)
+        assert gotref.sameSpace(ref)
+
+        with h5py.File('linear.x5', 'r') as f:
+            _check_metadata(f)
+            assert f.attrs['Type'] == 'linear'
+            _check_affine(f['/Transform'], xform)
+            _check_space( f['/A'],         src)
+            _check_space( f['/B'],         ref)
+
+
+def test_readWriteNonLinearX5():
+    datadir = op.join(op.dirname(__file__), 'testdata', 'nonlinear')
+    dffile  = op.join(datadir, 'displacementfield.nii.gz')
+    srcfile = op.join(datadir, 'src.nii.gz')
+    reffile = op.join(datadir, 'ref.nii.gz')
+
+    src     = fslimage.Image(srcfile)
+    ref     = fslimage.Image(reffile)
+    dfield  = fnirt.readFnirt(dffile, src, ref)
+    wdfield = nonlinear.convertDeformationSpace(dfield, 'world', 'world')
+
+    with tempdir.tempdir():
+
+        # field must be world->world
+        with pytest.raises(x5.X5Error):
+            x5.writeNonLinearX5('nonlinear.x5', dfield)
+
+        x5.writeNonLinearX5('nonlinear.x5', wdfield)
+
+        gotdfield = x5.readNonLinearX5('nonlinear.x5')
+
+        assert gotdfield.src.sameSpace(src)
+        assert gotdfield.ref.sameSpace(ref)
+        assert gotdfield.srcSpace == wdfield.srcSpace
+        assert gotdfield.refSpace == wdfield.refSpace
+        assert gotdfield.deformationType == wdfield.deformationType
+        assert np.all(np.isclose(gotdfield.data, wdfield.data))
+
+        with h5py.File('nonlinear.x5', 'r') as f:
+            assert f.attrs['Type'] == 'nonlinear'
+            _check_metadata(f)
+            _check_deformation(f['/Transform'], wdfield)
+            _check_space(      f['/A'],         ref)
+            _check_space(      f['/B'],         src)