diff --git a/nitransforms/conftest.py b/nitransforms/conftest.py
index 4bdbd723..c93daea9 100644
--- a/nitransforms/conftest.py
+++ b/nitransforms/conftest.py
@@ -5,6 +5,11 @@
 import nibabel as nb
 import pytest
 import tempfile
+import pkg_resources
+
+SOMEONES_ANATOMY = pkg_resources.resource_filename(
+    'nitransforms', 'tests/data/someones_anatomy.nii.gz')
+_data = None
 
 
 @pytest.fixture(autouse=True)
@@ -25,3 +30,41 @@ def doctest_autoimport(doctest_namespace):
     doctest_namespace['testfile'] = nifti_fname
     yield
     tmpdir.cleanup()
+
+
+@pytest.fixture
+def data_path():
+    """Return the test data folder."""
+    return os.path.join(os.path.dirname(__file__), 'tests/data')
+
+
+@pytest.fixture
+def get_data():
+    """Generate data in the requested orientation."""
+    global _data
+
+    if _data is not None:
+        return _data
+
+    img = nb.load(SOMEONES_ANATOMY)
+    imgaff = img.affine
+
+    _data = {'RAS': img}
+    newaff = imgaff.copy()
+    newaff[0, 0] *= -1.0
+    newaff[0, 3] = imgaff.dot(np.hstack((np.array(img.shape[:3]) - 1, 1.0)))[0]
+    _data['LAS'] = nb.Nifti1Image(np.flip(img.get_fdata(), 0), newaff, img.header)
+    newaff = imgaff.copy()
+    newaff[0, 0] *= -1.0
+    newaff[1, 1] *= -1.0
+    newaff[:2, 3] = imgaff.dot(np.hstack((np.array(img.shape[:3]) - 1, 1.0)))[:2]
+    _data['LPS'] = nb.Nifti1Image(np.flip(np.flip(img.get_fdata(), 0), 1), newaff, img.header)
+    A = nb.volumeutils.shape_zoom_affine(img.shape, img.header.get_zooms(), x_flip=False)
+    R = nb.affines.from_matvec(nb.eulerangles.euler2mat(x=0.09, y=0.001, z=0.001))
+    newaff = R.dot(A)
+    oblique_img = nb.Nifti1Image(img.get_fdata(), newaff, img.header)
+    oblique_img.header.set_qform(newaff, 1)
+    oblique_img.header.set_sform(newaff, 1)
+    _data['oblique'] = oblique_img
+
+    return _data
diff --git a/nitransforms/tests/data/affine-LAS-itk.tfm b/nitransforms/tests/data/affine-LAS-itk.tfm
deleted file mode 120000
index 849a5e16..00000000
--- a/nitransforms/tests/data/affine-LAS-itk.tfm
+++ /dev/null
@@ -1 +0,0 @@
-affine-RAS-itk.tfm
\ No newline at end of file
diff --git a/nitransforms/tests/data/affine-LAS.itk.tfm b/nitransforms/tests/data/affine-LAS.itk.tfm
new file mode 120000
index 00000000..629b9bf2
--- /dev/null
+++ b/nitransforms/tests/data/affine-LAS.itk.tfm
@@ -0,0 +1 @@
+affine-RAS.itk.tfm
\ No newline at end of file
diff --git a/nitransforms/tests/data/affine-LPS-itk.tfm b/nitransforms/tests/data/affine-LPS.itk.tfm
similarity index 100%
rename from nitransforms/tests/data/affine-LPS-itk.tfm
rename to nitransforms/tests/data/affine-LPS.itk.tfm
diff --git a/nitransforms/tests/data/affine-RAS-itk.tfm b/nitransforms/tests/data/affine-RAS.itk.tfm
similarity index 100%
rename from nitransforms/tests/data/affine-RAS-itk.tfm
rename to nitransforms/tests/data/affine-RAS.itk.tfm
diff --git a/nitransforms/tests/data/affine-oblique-itk.tfm b/nitransforms/tests/data/affine-oblique.itk.tfm
similarity index 100%
rename from nitransforms/tests/data/affine-oblique-itk.tfm
rename to nitransforms/tests/data/affine-oblique.itk.tfm
diff --git a/nitransforms/tests/test_transform.py b/nitransforms/tests/test_transform.py
index 86f404d3..cdad85d2 100644
--- a/nitransforms/tests/test_transform.py
+++ b/nitransforms/tests/test_transform.py
@@ -2,73 +2,135 @@
 # vi: set ft=python sts=4 ts=4 sw=4 et:
 """Tests of the transform module."""
 import os
-import numpy as np
-from numpy.testing import assert_array_equal, assert_almost_equal, \
-    assert_array_almost_equal
 import pytest
+import numpy as np
+from subprocess import check_call
 
-from nibabel.loadsave import load as loadimg
-from nibabel.nifti1 import Nifti1Image
+import nibabel as nb
 from nibabel.eulerangles import euler2mat
 from nibabel.affines import from_matvec
-from ..patched import shape_zoom_affine
-from .. import linear as nbl
-from nibabel.testing import (assert_equal, assert_not_equal, assert_true,
-                       assert_false, assert_raises,
-                       suppress_warnings, assert_dt_equal)
 from nibabel.tmpdirs import InTemporaryDirectory
+from .. import linear as nbl
+from .utils import assert_affines_by_filename
 
-data_path = os.path.join(os.path.dirname(__file__), 'data')
-SOMEONES_ANATOMY = os.path.join(data_path, 'someones_anatomy.nii.gz')
+TESTS_BORDER_TOLERANCE = 0.05
+APPLY_LINEAR_CMD = {
+    'fsl': """\
+flirt -setbackground 0 -interp nearestneighbour -in {moving} -ref {reference} \
+-applyxfm -init {transform} -out resampled.nii.gz\
+""".format,
+    'itk': """\
+antsApplyTransforms -d 3 -r {reference} -i {moving} \
+-o resampled.nii.gz -n NearestNeighbor -t {transform} --float\
+""".format,
+    'afni': """\
+3dAllineate -base {reference} -input {moving} \
+-prefix resampled.nii.gz -1Dmatrix_apply {transform} -final NN\
+""".format,
+}
 
 
-@pytest.mark.parametrize('image_orientation', ['RAS', 'LAS', 'LPS', 'oblique'])
-def test_affines_save(image_orientation):
-    """Check implementation of exporting affines to formats."""
-    # Generate test transform
-    img = loadimg(SOMEONES_ANATOMY)
-    imgaff = img.affine
-
-    if image_orientation == 'LAS':
-        newaff = imgaff.copy()
-        newaff[0, 0] *= -1.0
-        newaff[0, 3] = imgaff.dot(np.hstack((np.array(img.shape[:3]) - 1, 1.0)))[0]
-        img = Nifti1Image(np.flip(img.get_fdata(), 0), newaff, img.header)
-    elif image_orientation == 'LPS':
-        newaff = imgaff.copy()
-        newaff[0, 0] *= -1.0
-        newaff[1, 1] *= -1.0
-        newaff[:2, 3] = imgaff.dot(np.hstack((np.array(img.shape[:3]) - 1, 1.0)))[:2]
-        img = Nifti1Image(np.flip(np.flip(img.get_fdata(), 0), 1), newaff, img.header)
-    elif image_orientation == 'oblique':
-        A = shape_zoom_affine(img.shape, img.header.get_zooms(), x_flip=False)
-        R = from_matvec(euler2mat(x=0.09, y=0.001, z=0.001))
-        newaff = R.dot(A)
-        img = Nifti1Image(img.get_fdata(), newaff, img.header)
-        img.header.set_qform(newaff, 1)
-        img.header.set_sform(newaff, 1)
+@pytest.mark.parametrize('image_orientation', [
+    'RAS', 'LAS', 'LPS',  # 'oblique',
+])
+@pytest.mark.parametrize('sw_tool', ['itk', 'fsl', 'afni'])
+def test_linear_load(tmpdir, data_path, get_data, image_orientation, sw_tool):
+    """Check implementation of loading affines from formats."""
+    tmpdir.chdir()
 
+    img = get_data[image_orientation]
+    img.to_filename('img.nii.gz')
+
+    # Generate test transform
     T = from_matvec(euler2mat(x=0.9, y=0.001, z=0.001), [4.0, 2.0, -1.0])
+    xfm = nbl.Affine(T)
+    xfm.reference = img
+
+    ext = ''
+    if sw_tool == 'itk':
+        ext = '.tfm'
+
+    fname = 'affine-%s.%s%s' % (image_orientation, sw_tool, ext)
+    xfm_fname = os.path.join(data_path, fname)
+
+    if sw_tool == 'fsl':
+        with pytest.raises("ValueError"):
+            loaded = nbl.load(xfm_fname, fmt=fname.split('.')[-1])
+        with pytest.raises("ValueError"):
+            loaded = nbl.load(xfm_fname, fmt=fname.split('.')[-1],
+                              reference='img.nii.gz')
+        with pytest.raises("ValueError"):
+            loaded = nbl.load(xfm_fname, fmt=fname.split('.')[-1],
+                              moving='img.nii.gz')
+
+        loaded = nbl.load(xfm_fname, fmt=fname.split('.')[-1],
+                          moving='img.nii.gz', reference='img.nii.gz')
+    if sw_tool == 'afni':
+        with pytest.raises("ValueError"):
+            loaded = nbl.load(xfm_fname, fmt=fname.split('.')[-1])
+
+        loaded = nbl.load(xfm_fname, fmt=fname.split('.')[-1],
+                          reference='img.nii.gz')
 
+    assert loaded == xfm
+
+
+@pytest.mark.parametrize('image_orientation', [
+    'RAS', 'LAS', 'LPS',  # 'oblique',
+])
+@pytest.mark.parametrize('sw_tool', ['itk', 'fsl', 'afni'])
+def test_linear_save(data_path, get_data, image_orientation, sw_tool):
+    """Check implementation of exporting affines to formats."""
+    img = get_data[image_orientation]
+    # Generate test transform
+    T = from_matvec(euler2mat(x=0.9, y=0.001, z=0.001), [4.0, 2.0, -1.0])
     xfm = nbl.Affine(T)
     xfm.reference = img
 
-    itk = nbl.load(os.path.join(data_path, 'affine-%s-itk.tfm' % image_orientation),
-                   fmt='itk')
-    fsl = np.loadtxt(os.path.join(data_path, 'affine-%s.fsl' % image_orientation))
-    afni = np.loadtxt(os.path.join(data_path, 'affine-%s.afni' % image_orientation))
+    ext = ''
+    if sw_tool == 'itk':
+        ext = '.tfm'
 
     with InTemporaryDirectory():
-        xfm.to_filename('M.tfm', fmt='itk')
-        xfm.to_filename('M.fsl', fmt='fsl')
-        xfm.to_filename('M.afni', fmt='afni')
+        xfm_fname1 = 'M.%s%s' % (sw_tool, ext)
+        xfm.to_filename(xfm_fname1, fmt=sw_tool)
+
+        xfm_fname2 = os.path.join(
+            data_path, 'affine-%s.%s%s' % (image_orientation, sw_tool, ext))
+        assert_affines_by_filename(xfm_fname1, xfm_fname2)
+
+
+@pytest.mark.parametrize('image_orientation', [
+    'RAS', 'LAS', 'LPS',  # 'oblique',
+])
+@pytest.mark.parametrize('sw_tool', ['itk', 'fsl', 'afni'])
+def test_apply_linear_transform(tmpdir, data_path, get_data, image_orientation, sw_tool):
+    """Check implementation of exporting affines to formats."""
+    tmpdir.chdir()
+
+    img = get_data[image_orientation]
+    # Generate test transform
+    T = from_matvec(euler2mat(x=0.9, y=0.001, z=0.001), [4.0, 2.0, -1.0])
+    xfm = nbl.Affine(T)
+    xfm.reference = img
+
+    ext = ''
+    if sw_tool == 'itk':
+        ext = '.tfm'
 
-        nb_itk = nbl.load('M.tfm', fmt='itk')
-        nb_fsl = np.loadtxt('M.fsl')
-        nb_afni = np.loadtxt('M.afni')
+    img.to_filename('img.nii.gz')
+    xfm_fname = 'M.%s%s' % (sw_tool, ext)
+    xfm.to_filename(xfm_fname, fmt=sw_tool)
 
-    assert_equal(itk, nb_itk)
-    assert_almost_equal(fsl, nb_fsl)
-    assert_almost_equal(afni, nb_afni)
+    cmd = APPLY_LINEAR_CMD[sw_tool](
+        transform=os.path.abspath(xfm_fname),
+        reference=os.path.abspath('img.nii.gz'),
+        moving=os.path.abspath('img.nii.gz'))
+    exit_code = check_call([cmd], shell=True)
+    assert exit_code == 0
+    sw_moved = nb.load('resampled.nii.gz')
 
-# Create version not aligned to canonical
\ No newline at end of file
+    nt_moved = xfm.resample(img, order=0)
+    diff = sw_moved.get_fdata() - nt_moved.get_fdata()
+    # A certain tolerance is necessary because of resampling at borders
+    assert (np.abs(diff) > 1e-3).sum() / diff.size < TESTS_BORDER_TOLERANCE
diff --git a/nitransforms/tests/utils.py b/nitransforms/tests/utils.py
new file mode 100644
index 00000000..3bea4357
--- /dev/null
+++ b/nitransforms/tests/utils.py
@@ -0,0 +1,21 @@
+"""Utilities for testing."""
+from pathlib import Path
+import numpy as np
+
+from .. import linear as nbl
+
+
+def assert_affines_by_filename(affine1, affine2):
+    """Check affines by filename."""
+    affine1 = Path(affine1)
+    affine2 = Path(affine2)
+    assert affine1.suffix == affine2.suffix, 'affines of different type'
+
+    if affine1.suffix.endswith('.tfm'):  # An ITK transform
+        xfm1 = nbl.load(str(affine1), fmt='itk')
+        xfm2 = nbl.load(str(affine2), fmt='itk')
+        assert xfm1 == xfm2
+    else:
+        xfm1 = np.loadtxt(str(affine1))
+        xfm2 = np.loadtxt(str(affine2))
+        np.testing.assert_almost_equal(xfm1, xfm2)