diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 019e6e6eff6f74ab7c951eeb9f32b435507fb3fb..6c8862394e6cc4b415c5ed8e2208bcab94674c19 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -22,6 +22,8 @@ Fixed
 
 
 * Improved the algorithm used by the :func:`.mesh.needsFixing` function.
+* The :meth:`.fslmaths.run` method now accepts :attr:`.wrappers.LOAD` as an
+  output specification.
 
 
 Deprecated
diff --git a/fsl/wrappers/__init__.py b/fsl/wrappers/__init__.py
index b37413ed3678a5fbd84f9056bdae2d52f76a64dd..453313f2ca14e592732cf8866cabd78065413d29 100644
--- a/fsl/wrappers/__init__.py
+++ b/fsl/wrappers/__init__.py
@@ -63,11 +63,11 @@ if we want to FLIRT two images and get the result, we can do this::
 Similarly, we can run a ``fslmaths`` command on in-memory images::
 
     import nibabel as nib
-    from fsl.wrappers import fslmaths, LOAD
+    from fsl.wrappers import fslmaths
 
     image  = nib.load('image.nii')
     mask   = nib.load('mask.nii')
-    output = fslmaths(image).mas(mask).bin().run(LOAD)
+    output = fslmaths(image).mas(mask).bin().run()
 
 
 If you are *writing* wrapper functions, take a look at the
diff --git a/fsl/wrappers/fslmaths.py b/fsl/wrappers/fslmaths.py
index 5fa790ab46d6ef472874c4ab63e160a04c59ed7f..2bfd7129183739b79ae7af7b150f417c862e9244 100644
--- a/fsl/wrappers/fslmaths.py
+++ b/fsl/wrappers/fslmaths.py
@@ -13,7 +13,28 @@ from . import wrapperutils as wutils
 
 
 class fslmaths(object):
-    """Perform mathematical manipulation of images."""
+    """Perform mathematical manipulation of images.
+
+    ``fslmaths`` is unlike the other FSL wrapper tools in that it provides an
+    object-oriented method-chaining interface, which is hopefully easier to
+    use than constructing a ``fslmaths`` command-line call. For example, the
+    following call to the ``fslmaths`` wrapper function::
+
+        fslmaths('input.nii').thr(0.25).mul(-1).run('output.nii')
+
+    will be translated into the following command-line call::
+
+        fslmaths input.nii -thr 0.25 -mul -1 output.nii
+
+    The ``fslmaths`` wrapper function can also be used with in-memory
+    images. If no output file name is passed to the :meth:`run` method, the
+    result is loaded into memory and returned as a ``nibabel`` image.  For
+    example::
+
+        import nibabel as nib
+        input  = nib.load('input.nii')
+        output = fslmaths(input).thr(0.25).mul(-1).run()
+    """
 
     def __init__(self, input):
         """Constructor."""
@@ -138,20 +159,24 @@ class fslmaths(object):
         return self
 
     def run(self, output=None):
-        """Save output of operations to image."""
+        """Save output of operations to image. Set ``output`` to a filename to have
+        the result saved to file, or omit ``output`` entirely to have the
+        result returned as a ``nibabel`` image.
+        """
 
         cmd = ['fslmaths', self.__input] + self.__args
 
-        if output is None: cmd += [wutils.LOAD]
-        else:              cmd += [output]
+        if output is None:
+            output = wutils.LOAD
 
+        cmd   += [output]
         result = self.__run(*cmd)
 
         # if output is LOADed, there
         # will only be one entry in
         # the result dict.
-        if output is None: return list(result.values())[0]
-        else:              return result
+        if output == wutils.LOAD: return list(result.values())[0]
+        else:                     return result
 
     @wutils.fileOrImage()
     @wutils.fslwrapper
diff --git a/tests/__init__.py b/tests/__init__.py
index 2ed387e0baab5e22feae602f190801997b199be2..24587392fc29de81bea199cb8e991b6012cba22e 100644
--- a/tests/__init__.py
+++ b/tests/__init__.py
@@ -50,7 +50,10 @@ def mockFSLDIR(**kwargs):
                 if not op.isdir(subdir):
                     os.makedirs(subdir)
                 for fname in files:
-                    touch(op.join(subdir, fname))
+                    fname = op.join(subdir, fname)
+                    touch(fname)
+                    if subdir == bindir:
+                        os.chmod(fname, 0o755)
             fslplatform.fsldir = fsldir
             fslplatform.fsldevdir = None
 
diff --git a/tests/test_wrappers.py b/tests/test_wrappers.py
index fc0050a31b3ecde9a4e343927917eb8fdf78f51c..9f4819af0580ac61061e6e12dfa7306be9bf779e 100644
--- a/tests/test_wrappers.py
+++ b/tests/test_wrappers.py
@@ -5,14 +5,19 @@
 # Author: Paul McCarthy <pauldmccarthy@gmail.com>
 #
 
+import              os
 import os.path   as op
 import itertools as it
+import textwrap  as tw
+
+import numpy as np
 
 import fsl.wrappers                       as fw
 import fsl.utils.assertions               as asrt
 import fsl.utils.run                      as run
+from fsl.utils.tempdir import tempdir
 
-from . import mockFSLDIR
+from . import mockFSLDIR, make_random_image
 
 
 def checkResult(cmd, base, args, stripdir=None):
@@ -272,7 +277,23 @@ def test_fslmaths():
 
         assert result.output[0] == expected
 
-        # TODO test LOAD output
+    # test LOAD output
+    with tempdir() as td, mockFSLDIR(bin=('fslmaths',)) as fsldir:
+        expect = make_random_image(op.join(td, 'output.nii.gz'))
+
+        with open(op.join(fsldir, 'bin', 'fslmaths'), 'wt') as f:
+            f.write(tw.dedent("""
+            #!/usr/bin/env python
+            import sys
+            import shutil
+            shutil.copy('{}', sys.argv[2])
+            """.format(op.join(td, 'output.nii.gz'))).strip())
+            os.chmod(op.join(fsldir, 'bin', 'fslmaths'), 0o755)
+
+        got = fw.fslmaths('input').run()
+        assert np.all(expect.dataobj[:] == got.dataobj[:])
+        got = fw.fslmaths('input').run(fw.LOAD)
+        assert np.all(expect.dataobj[:] == got.dataobj[:])
 
 def test_fast():
     with asrt.disabled(), run.dryrun(), mockFSLDIR(bin=('fast',)) as fsldir: