diff --git a/nitransforms/base.py b/nitransforms/base.py
index 420fc3ce..6c1eba17 100644
--- a/nitransforms/base.py
+++ b/nitransforms/base.py
@@ -193,4 +193,4 @@ def to_filename(self, filename, fmt='X5'):
             root = out_file.create_group('/0')
             self._to_hdf5(root)
 
-        return filename
\ No newline at end of file
+        return filename
diff --git a/nitransforms/conftest.py b/nitransforms/conftest.py
new file mode 100644
index 00000000..4bdbd723
--- /dev/null
+++ b/nitransforms/conftest.py
@@ -0,0 +1,27 @@
+"""py.test configuration."""
+import os
+from pathlib import Path
+import numpy as np
+import nibabel as nb
+import pytest
+import tempfile
+
+
+@pytest.fixture(autouse=True)
+def doctest_autoimport(doctest_namespace):
+    """Make available some fundamental modules to doctest modules."""
+    doctest_namespace['np'] = np
+    doctest_namespace['nb'] = nb
+    doctest_namespace['os'] = os
+    doctest_namespace['Path'] = Path
+    doctest_namespace['datadir'] = os.path.join(os.path.dirname(__file__), 'tests/data')
+
+    tmpdir = tempfile.TemporaryDirectory()
+    doctest_namespace['tmpdir'] = tmpdir.name
+
+    testdata = np.zeros((11, 11, 11), dtype='uint8')
+    nifti_fname = str(Path(tmpdir.name) / 'test.nii.gz')
+    nb.Nifti1Image(testdata, np.eye(4)).to_filename(nifti_fname)
+    doctest_namespace['testfile'] = nifti_fname
+    yield
+    tmpdir.cleanup()
diff --git a/nitransforms/linear.py b/nitransforms/linear.py
index 9f926ee5..9e0d379c 100644
--- a/nitransforms/linear.py
+++ b/nitransforms/linear.py
@@ -30,6 +30,7 @@ class Affine(TransformBase):
     def __init__(self, matrix=None, reference=None):
         """
         Initialize a linear transform.
+
         Parameters
         ----------
         matrix : ndarray
@@ -37,14 +38,16 @@ def __init__(self, matrix=None, reference=None):
             coordinates**, mapping coordinates from *reference* space
             into *moving* space.
             This matrix should be provided in homogeneous coordinates.
+
         Examples
         --------
         >>> xfm = Affine([[1, 0, 0, 4], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]])
         >>> xfm.matrix  # doctest: +NORMALIZE_WHITESPACE
-        array([[1, 0, 0, 4],
-               [0, 1, 0, 0],
-               [0, 0, 1, 0],
-               [0, 0, 0, 1]])
+        array([[[1, 0, 0, 4],
+                [0, 1, 0, 0],
+                [0, 0, 1, 0],
+                [0, 0, 0, 1]]])
+
         """
         super(Affine, self).__init__()
         if matrix is None:
@@ -71,6 +74,7 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True,
                  output_dtype=None):
         """
         Resample the moving image in reference space.
+
         Parameters
         ----------
         moving : `spatialimage`
@@ -92,17 +96,27 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True,
             slightly blurred if *order > 1*, unless the input is prefiltered,
             i.e. it is the result of calling the spline filter on the original
             input.
+
         Returns
         -------
         moved_image : `spatialimage`
             The moving imaged after resampling to reference space.
+
         Examples
         --------
-        >>> import nibabel as nib
         >>> xfm = Affine([[1, 0, 0, 4], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]])
-        >>> ref = nib.load('image.nii.gz')
+        >>> ref = nb.load(testfile)
         >>> xfm.reference = ref
-        >>> xfm.resample(ref, order=0)
+        >>> refdata = ref.get_fdata()
+        >>> np.allclose(refdata, 0)
+        True
+
+        >>> refdata[5, 5, 5] = 1  # Set a one in the middle voxel
+        >>> moving = nb.Nifti1Image(refdata, ref.affine, ref.header)
+        >>> resampled = xfm.resample(moving, order=0).get_fdata()
+        >>> resampled[1, 5, 5]
+        1.0
+
         """
         if output_dtype is None:
             output_dtype = moving.header.get_data_dtype()
@@ -323,4 +337,4 @@ def _fsl_aff_adapt(space):
     if np.linalg.det(aff) > 0:
         swp[0, 0] = -1.0
         swp[0, 3] = (space.shape[0] - 1) * zooms[0]
-    return swp, np.diag(zooms)
\ No newline at end of file
+    return swp, np.diag(zooms)
diff --git a/nitransforms/nonlinear.py b/nitransforms/nonlinear.py
index 90aa9d0b..c9cf543e 100644
--- a/nitransforms/nonlinear.py
+++ b/nitransforms/nonlinear.py
@@ -102,18 +102,25 @@ def _cache_moving(self, moving):
     def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True,
                  output_dtype=None):
         """
-        Resample the `moving` image applying the deformation field.
+        Resample the ``moving`` image applying the deformation field.
+
         Examples
         --------
-        >>> import numpy as np
-        >>> import nibabel as nb
-        >>> ref = nb.load('t1_weighted.nii.gz')
+        >>> ref = nb.load(testfile)
+        >>> refdata = ref.get_fdata()
+        >>> np.allclose(refdata, 0)
+        True
+
+        >>> refdata[5, 5, 5] = 1  # Set a one in the middle voxel
+        >>> moving = nb.Nifti1Image(refdata, ref.affine, ref.header)
         >>> field = np.zeros(tuple(list(ref.shape) + [3]))
         >>> field[..., 0] = 4.0
         >>> fieldimg = nb.Nifti1Image(field, ref.affine, ref.header)
-        >>> xfm = nb.transform.DeformationFieldTransform(fieldimg)
-        >>> new = xfm.resample(ref)
-        >>> new.to_filename('deffield.nii.gz')
+         >>> xfm = DeformationFieldTransform(fieldimg)
+        >>> resampled = xfm.resample(moving, order=0).get_fdata()
+        >>> resampled[1, 5, 5]
+        1.0
+
         """
         self._cache_moving(moving)
         return super(DeformationFieldTransform, self).resample(
@@ -222,20 +229,19 @@ def map_voxel(self, index, moving=None):
     def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True,
                  output_dtype=None):
         """
-        Resample the `moving` image applying the deformation field.
+        Resample the ``moving`` image applying the deformation field.
+
         Examples
         --------
-        >>> import numpy as np
-        >>> import nibabel as nb
-        >>> ref = nb.load('t1_weighted.nii.gz')
+        >>> ref = nb.load(os.path.join(datadir, 'someones_anatomy.nii.gz'))
         >>> coeffs = np.zeros((6, 6, 6, 3))
         >>> coeffs[2, 2, 2, ...] = [10.0, -20.0, 0]
         >>> aff = ref.affine
-        >>> aff[:3, :3] = aff.dot(np.eye(3) * np.array(ref.header.get_zooms()[:3] / 6.0)
+        >>> aff[:3, :3] = aff[:3, :3].dot(np.eye(3) * np.array(ref.header.get_zooms()[:3]) / 6.0)
         >>> coeffsimg = nb.Nifti1Image(coeffs, ref.affine, ref.header)
-        >>> xfm = nb.transform.BSplineFieldTransform(ref, coeffsimg)
-        >>> new = xfm.resample(ref)
-        >>> new.to_filename('deffield.nii.gz')
+        >>> xfm = BSplineFieldTransform(ref, coeffsimg)
+        >>> new = xfm.resample(ref)  # doctest: +SKIP
+
         """
         self._cache_moving()
         return super(BSplineFieldTransform, self).resample(
@@ -243,4 +249,4 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True,
 
 
 def _pprint(inlist):
-    return 'x'.join(['%d' % s for s in inlist])
\ No newline at end of file
+    return 'x'.join(['%d' % s for s in inlist])