diff --git a/fsl/wrappers/__init__.py b/fsl/wrappers/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..3bf0f2a63670655824cf4c47aed717615fb1748d
--- /dev/null
+++ b/fsl/wrappers/__init__.py
@@ -0,0 +1,32 @@
+#!/usr/bin/env python
+#
+# __init__.py - Wrappers for FSL command-line tools.
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+"""This package contains wrappers for various FSL command line tools, allowing
+them to be called from Python.
+"""
+
+
+from .bet      import (bet,)            # noqa
+from .eddy     import (eddy_cuda,       # noqa
+                       topup)
+from .flirt    import (flirt,           # noqa
+                       invxfm,
+                       applyxfm,
+                       concatxfm,
+                       mcflirt)
+from .fnirt    import (fnirt,           # noqa
+                       applywarp,
+                       invwarp,
+                       convertwarp)
+from .fslmaths import (fslmaths,)       # noqa
+from .fugue    import (fugue,           # noqa
+                       sigloss)
+from .melodic  import (melodic,         # noqa
+                       fsl_regfilt)
+from .misc     import (fslreorient2std, # noqa
+                       fslroi,
+                       slicer,
+                       cluster)
diff --git a/fsl/wrappers/bet.py b/fsl/wrappers/bet.py
new file mode 100644
index 0000000000000000000000000000000000000000..58c4aab3b78aa64c8fa9eb73c18e3c05b9b6e1ad
--- /dev/null
+++ b/fsl/wrappers/bet.py
@@ -0,0 +1,46 @@
+#!/usr/bin/env python
+#
+# bet.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+
+import fsl.utils.assertions as asrt
+import fsl.utils.run        as run
+from . import wrapperutils  as wutils
+
+
+@wutils.fileOrImage('input', 'output')
+def bet(input, output, **kwargs):
+    """Delete non-brain tissue from an image of the whole head.
+
+    :arg input:         Required
+    :arg output:        Required
+    :arg mask:
+    :arg robust:
+    :arg fracintensity:
+    :arg seg:
+
+    Refer to the ``bet`` command-line help for details on all arguments.
+    """
+
+    asrt.assertIsNifti(input)
+
+    argmap = {
+        'mask'          : 'm',
+        'robust'        : 'R',
+        'fracintensity' : 'f',
+        'seg'           : 'n',
+    }
+
+    valmap = {
+        'm' : wutils.SHOW_IF_TRUE,
+        'R' : wutils.SHOW_IF_TRUE,
+        'n' : wutils.HIDE_IF_TRUE,
+    }
+
+    cmd  = ['bet', input, output]
+    cmd += wutils.applyArgStyle('-', argmap, valmap, **kwargs)
+
+    return run.runfsl(*cmd)
diff --git a/fsl/wrappers/eddy.py b/fsl/wrappers/eddy.py
new file mode 100644
index 0000000000000000000000000000000000000000..1081edb6d244eea27e153e5fc38fea21dcdf5f28
--- /dev/null
+++ b/fsl/wrappers/eddy.py
@@ -0,0 +1,113 @@
+#!/usr/bin/env python
+#
+# eddy.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+
+def eddy_cuda(imain, mask, index, acqp, bvecs, bvals, out, very_verbose=False,
+              niter=None, fwhm=None, s2v_niter=None, mporder=None, nvoxhp=None,
+              slspec=None, b0_only=False, topup=None, field=None,
+              field_mat=None, debug=None, s2v_fwhm=None, interp=None,
+              dont_mask_output=False, s2v_interp=None, ref_scan_no=None,
+              data_is_shelled=False, estimate_move_by_susceptibility=False,
+              mbs_niter=None, mbs_lambda=None, mbs_ksp=None, s2v_lambda=None,
+              cnr_maps=False, residuals=False):
+    """Correct eddy current-induced distortions and subject movements."""
+    asrt.assertFileExists(imain, mask, index, acqp, bvecs, bvals)
+    asrt.assertIsNifti(imain, mask)
+
+    assert not (topup and field), "topup and field arguments are incompatible"
+
+    out = img.splitext(out)[0]
+
+    opts = " --imain={0} --mask={1} --index={2} --bvals={3} " \
+           "--bvecs={4} --acqp={5} --out={6}".format(imain, mask, index, bvals,
+                                                     bvecs, acqp, out)
+
+    cmd = op.join(os.getenv('DHCP_EDDY_PATH', ''), 'eddy_cuda')
+    # cmd = 'eddy_cuda'
+    cmd = cmd + opts
+
+    if very_verbose:
+        cmd += " --very_verbose"
+    if estimate_move_by_susceptibility:
+        cmd += " --estimate_move_by_susceptibility"
+    if data_is_shelled:
+        cmd += " --data_is_shelled"
+    if mbs_niter is not None:
+        cmd += " --mbs_niter={0}".format(mbs_niter)
+    if mbs_lambda is not None:
+        cmd += " --mbs_lambda={0}".format(mbs_lambda)
+    if mbs_ksp is not None:
+        cmd += " --mbs_ksp={0}".format(mbs_ksp)
+    if niter is not None:
+        cmd += " --niter={0}".format(niter)
+    if fwhm is not None:
+        cmd += " --fwhm={0}".format(fwhm)
+    if s2v_fwhm is not None:
+        cmd += " --s2v_fwhm={0}".format(s2v_fwhm)
+    if s2v_niter is not None:
+        cmd += " --s2v_niter={0}".format(s2v_niter)
+    if s2v_interp is not None:
+        cmd += " --s2v_interp={0}".format(s2v_interp)
+    if interp is not None:
+        cmd += " --interp={0}".format(interp)
+    if mporder is not None:
+        cmd += " --mporder={0}".format(mporder)
+    if nvoxhp is not None:
+        cmd += " --nvoxhp={0}".format(nvoxhp)
+    if slspec is not None:
+        cmd += " --slspec={0}".format(slspec)
+    if topup is not None:
+        cmd += " --topup={0}".format(topup)
+    if field is not None:
+        field = img.splitext(field)[0]
+        cmd += " --field={0}".format(field)
+    if b0_only:
+        cmd += " --b0_only"
+    if field_mat is not None:
+        cmd += " --field_mat={0}".format(field_mat)
+    if debug is not None:
+        cmd += " --debug={0}".format(debug)
+    if dont_mask_output:
+        cmd += " --dont_mask_output"
+    if ref_scan_no is not None:
+        cmd += " --ref_scan_no={0}".format(ref_scan_no)
+    if s2v_lambda is not None:
+        cmd += " --s2v_lambda={0}".format(s2v_lambda)
+    if cnr_maps:
+        cmd += " --cnr_maps"
+    if residuals:
+        cmd += " --residuals"
+
+    shops.run(cmd)
+
+
+
+
+def topup(imain, datain, config=None, fout=None, iout=None, out=None,
+          verbose=False, subsamp=None, logout=None):
+    """Estimate and correct susceptibility induced distortions."""
+    asrt.assertFileExists(imain, datain)
+    asrt.assertIsNifti(imain)
+
+    cmd = "topup --imain={0} --datain={1}".format(imain, datain)
+
+    if config is not None:
+        cmd += " --config={0}".format(config)
+    if fout is not None:
+        cmd += " --fout={0}".format(fout)
+    if iout is not None:
+        cmd += " --iout={0}".format(iout)
+    if out is not None:
+        cmd += " --out={0}".format(out)
+    if subsamp is not None:
+        cmd += " --subsamp={0}".format(subsamp)
+    if logout is not None:
+        cmd += " --logout={0}".format(logout)
+    if verbose:
+        cmd += " -v"
+
+    shops.run(cmd)
diff --git a/fsl/wrappers/flirt.py b/fsl/wrappers/flirt.py
new file mode 100644
index 0000000000000000000000000000000000000000..9247ce4557e74bf4b9103ace52c3c279739646af
--- /dev/null
+++ b/fsl/wrappers/flirt.py
@@ -0,0 +1,124 @@
+#!/usr/bin/env python
+#
+# flirt.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+"""
+
+.. autosummary::
+   :nosignatures:
+
+   flirt
+   invxfm
+   applyxfm
+   concatxfm
+   mcflirt
+"""
+
+
+import fsl.utils.run        as run
+import fsl.utils.assertions as asrt
+import fsl.data.image       as fslimage
+
+
+def flirt(src, ref, out=None, omat=None, dof=None, cost=None, wmseg=None,
+          init=None, schedule=None, echospacing=None, pedir=None,
+          fieldmap=None, fieldmapmask=None, bbrslope=None, bbrtype=None,
+          interp=None, refweight=None, applyisoxfm=None, usesqform=False,
+          nosearch=False, verbose=0):
+    """FLIRT (FMRIB's Linear Image Registration Tool)."""
+    asrt.assertIsNifti(src, ref)
+    asrt.assertFileExists(src, ref)
+
+    cmd = "flirt -in {0} -ref {1}".format(src, ref)
+
+    if out is not None:
+        asrt.assertIsNifti(out)
+        cmd += " -out {0}".format(out)
+    if omat is not None:
+        cmd += " -omat {0}".format(omat)
+    if dof is not None:
+        cmd += " -dof {0}".format(dof)
+    if cost is not None:
+        cmd += " -cost {0}".format(cost)
+    if wmseg is not None:
+        asrt.assertIsNifti(wmseg)
+        cmd += " -wmseg {0}".format(wmseg)
+    if init is not None:
+        cmd += " -init {0}".format(init)
+    if schedule is not None:
+        cmd += " -schedule {0}".format(schedule)
+    if echospacing is not None:
+        cmd += " -echospacing {0}".format(echospacing)
+    if pedir is not None:
+        cmd += " -pedir {0}".format(pedir)
+    if fieldmap is not None:
+        cmd += " -fieldmap {0}".format(fieldmap)
+    if fieldmapmask is not None:
+        cmd += " -fieldmapmask {0}".format(fieldmapmask)
+    if bbrslope is not None:
+        cmd += " -bbrslope {0}".format(bbrslope)
+    if bbrtype is not None:
+        cmd += " -bbrtype {0}".format(bbrtype)
+    if interp is not None:
+        cmd += " -interp {0}".format(interp)
+    if refweight is not None:
+        asrt.assertIsNifti(refweight)
+        cmd += " -refweight {0}".format(refweight)
+    if applyisoxfm is not None:
+        cmd += " -applyisoxfm {0}".format(applyisoxfm)
+    if verbose is not None:
+        cmd += " -verbose {0}".format(verbose)
+    if usesqform:
+        cmd += " -usesqform"
+    if nosearch:
+        cmd += " -nosearch"
+
+    return run.runfsl(cmd)
+
+
+def invxfm(inmat, omat):
+    """Tool for inverting FSL transformation matrices."""
+    asrt.assertFileExists(inmat)
+
+    cmd = "convert_xfm -omat {0} -inverse {1}".format(omat, inmat)
+    return run.runfsl(cmd)
+
+
+def applyxfm(src, ref, mat, out, interp='spline'):
+    """Tool for applying FSL transformation matrices."""
+    asrt.assertFileExists(src, ref)
+    asrt.assertIsNifti(src, ref)
+
+    cmd = "flirt -init {0} -in {1} -ref {2} -applyxfm -out {3} -interp {4}"
+    return run.runfsl(cmd.format(mat, src, ref, out, interp))
+
+
+def concatxfm(inmat1, inmat2, outmat):
+    """Tool to concatenate two FSL transformation matrices."""
+    asrt.assertFileExists(inmat1, inmat2)
+
+    cmd = "convert_xfm -omat {0} -concat {1} {2}"
+    return run.runfsl(cmd.format(outmat, inmat2, inmat1))
+
+
+def mcflirt(infile, outfile, reffile=None, spline_final=True, plots=True,
+            mats=True, refvol=None):
+    """Rigid-body motion correction using mcflirt."""
+
+    outfile = fslimage.removeExt(outfile)
+    cmd = "mcflirt -in {0} -out {1} -rmsrel -rmsabs".format(infile, outfile)
+
+    if reffile is not None:
+        cmd += " -reffile {0}".format(reffile)
+    if refvol is not None:
+        cmd += " -refvol {0}".format(refvol)
+    if spline_final:
+        cmd += " -spline_final"
+    if plots:
+        cmd += " -plots"
+    if mats:
+        cmd += " -mats"
+
+    return run.runfsl(cmd)
diff --git a/fsl/wrappers/fnirt.py b/fsl/wrappers/fnirt.py
new file mode 100644
index 0000000000000000000000000000000000000000..17c42b036eb181b3d4df7f8daf98641055dd94b7
--- /dev/null
+++ b/fsl/wrappers/fnirt.py
@@ -0,0 +1,200 @@
+#!/usr/bin/env python
+#
+# fnirt.py - FNIRT wrapper functions.
+#
+# Author: Sean Fitzgibbon <sean.fitzgibbon@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+"""This module provides wrapper functions for the FSL `FNIRT
+<https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FNIRT/>`_ tool, used for performing
+non-linear registration of 3D images.
+
+.. autosummary::
+   :nosignatures:
+
+   fnirt
+   applywarp
+   invwarp
+   convertwarp
+"""
+
+
+import os.path as op
+import            glob
+
+
+import fsl.utils.run        as run
+import fsl.utils.assertions as asrt
+
+
+def fnirt(src, ref, aff=None, imprefm=None, impinm=None, applyrefmask=None,
+          applyinmask=None, subsamp=None, miter=None, infwhm=None,
+          reffwhm=None, lmbda=None, estint=None, warpres=None, ssqlambda=None,
+          regmod=None, intmod=None, intorder=None, biasres=None,
+          biaslambda=None, refderiv=None, cout=None, intout=None, refout=None,
+          iout=None, interp=None, inwarp=None, minmet=None, verbose=False,
+          intin=None, jout=None):
+    """Do nonlinear image registration."""
+    cmd = 'fnirt --in={0} --ref={1}'.format(src, ref)
+
+    if aff is not None:
+        cmd += " --aff={0}".format(aff)
+    if imprefm is not None:
+        cmd += " --imprefm={0}".format(imprefm)
+    if impinm is not None:
+        cmd += " --impinm={0}".format(impinm)
+    if applyrefmask is not None:
+        cmd += " --applyrefmask={0}".format(applyrefmask)
+    if applyinmask is not None:
+        cmd += " --applyinmask={0}".format(applyinmask)
+    if subsamp is not None:
+        cmd += " --subsamp={0}".format(subsamp)
+    if miter is not None:
+        cmd += " --miter={0}".format(miter)
+    if infwhm is not None:
+        cmd += " --infwhm={0}".format(infwhm)
+    if reffwhm is not None:
+        cmd += " --reffwhm={0}".format(reffwhm)
+    if lmbda is not None:
+        cmd += " --lambda={0}".format(lmbda)
+    if estint is not None:
+        cmd += " --estint={0}".format(estint)
+    if warpres is not None:
+        cmd += " --warpres={0}".format(warpres)
+    if ssqlambda is not None:
+        cmd += " --ssqlambda={0}".format(ssqlambda)
+    if regmod is not None:
+        cmd += " --regmod={0}".format(regmod)
+    if intmod is not None:
+        cmd += " --intmod={0}".format(intmod)
+    if intorder is not None:
+        cmd += " --intorder={0}".format(intorder)
+    if biasres is not None:
+        cmd += " --biasres={0}".format(biasres)
+    if biaslambda is not None:
+        cmd += " --biaslambda={0}".format(biaslambda)
+    if refderiv is not None:
+        cmd += " --refderiv={0}".format(refderiv)
+    if cout is not None:
+        cmd += " --cout={0}".format(cout)
+    if intout is not None:
+        cmd += " --intout={0}".format(intout)
+    if refout is not None:
+        cmd += " --refout={0}".format(refout)
+    if iout is not None:
+        cmd += " --iout={0}".format(iout)
+    if interp is not None:
+        cmd += " --interp={0}".format(interp)
+    if inwarp is not None:
+        cmd += " --inwarp={0}".format(inwarp)
+    if minmet is not None:
+        cmd += " --minmet={0}".format(minmet)
+    if intin is not None:
+        cmd += " --intin={0}".format(intin)
+    if jout is not None:
+        cmd += " --jout={0}".format(jout)
+    if verbose:
+        cmd += " --verbose"
+
+    return run.runfsl(cmd)
+
+
+def applywarp(src, ref, out, warp=None, premat=None, prematdir=None,
+              postmat=None, postmatdir=None, interp="spline",
+              paddingsize=None, abs=False, rel=False):
+    """Tool for applying FSL warps.
+
+    The ``prematdir`` and ``postmatdir`` arguments can be used when warping a
+    4D image. You can specify a directory containing multiple affines, named
+    ``MAT_????`` (as output by e.g. ``mcflirt``). Each file will be applied to
+    one volume (in order) of the image.
+    """
+    assert (warp or premat or postmat or prematdir or postmatdir), \
+        "either a warp or mat (premat, postmat or prematdir) must be supplied"
+    assert not (premat and prematdir), \
+        "cannot use premat and prematdir arguments together"
+    assert not (postmat and postmatdir), \
+        "cannot use postmat and postmatdir arguments together"
+
+    def catmats(matdir, out):
+        """Concatenate FSL trasformations files into a single file."""
+        mats = sorted(glob.glob(op.join(matdir, "MAT_????")))
+        with open(out, 'w') as outfile:
+            for fname in mats:
+                with open(fname) as infile:
+                    outfile.write(infile.read())
+
+    cmd = "--in={0} --ref={1} --out={2} --interp={3}".format(src, ref, out,
+                                                             interp)
+    cmd = "applywarp " + cmd
+
+    if prematdir:
+        premat = op.join(prematdir, 'allmats.txt')
+        catmats(prematdir, premat)
+    if postmatdir:
+        postmat = op.join(postmatdir, 'allmats.txt')
+        catmats(postmatdir, postmat)
+    if warp:
+        cmd += " --warp={0}".format(warp)
+    if premat:
+        cmd += " --premat={0}".format(premat)
+    if postmat:
+        cmd += " --postmat={0}".format(postmat)
+    if paddingsize:
+        cmd += " --paddingsize={0}".format(paddingsize)
+    if abs:
+        cmd += " --abs"
+    if rel:
+        cmd += " --rel"
+
+    return run.runfsl(cmd)
+
+
+def invwarp(inwarp, ref, outwarp):
+    """Tool for inverting FSL warps."""
+
+    asrt.assertFileExists(inwarp, ref)
+    asrt.assertIsNifti(inwarp, ref, outwarp)
+
+    cmd  = 'invwarp'
+    cmd += ' --warp={}'.format(inwarp)
+    cmd += ' --out={}'.format(outwarp)
+    cmd += ' --ref={}'.format(ref)
+
+    return run.runfsl(cmd)
+
+
+def convertwarp(out, ref, warp1=None, warp2=None, premat=None, midmat=None,
+                postmat=None, shiftmap=None, shiftdir=None, absout=False,
+                abs=False, rel=False, relout=False):
+    """Tool for converting FSL warps."""
+
+    assert (warp1 or warp2 or premat or midmat or postmat), \
+        "either a warp (warp1 or warp2) or mat (premat, midmat, or " + \
+        "postmat) must be supplied"
+
+    cmd = "convertwarp --ref={0} --out={1}".format(ref, out)
+    if warp1:
+        cmd = cmd + " --warp1={0}".format(warp1)
+    if warp2:
+        cmd = cmd + " --warp2={0}".format(warp2)
+    if premat:
+        cmd = cmd + " --premat={0}".format(premat)
+    if midmat:
+        cmd = cmd + " --midmat={0}".format(midmat)
+    if postmat:
+        cmd = cmd + " --postmat={0}".format(postmat)
+    if shiftmap:
+        cmd = cmd + " --shiftmap={0}".format(shiftmap)
+    if shiftdir:
+        cmd = cmd + " --shiftdir={0}".format(shiftdir)
+    if absout:
+        cmd = cmd + " --absout"
+    if relout:
+        cmd = cmd + " --relout"
+    if abs:
+        cmd = cmd + " --abs"
+    if rel:
+        cmd = cmd + " --rel"
+
+    return run.runfsl(cmd)
diff --git a/fsl/wrappers/fslmaths.py b/fsl/wrappers/fslmaths.py
new file mode 100644
index 0000000000000000000000000000000000000000..614d0c68e22f17b07335209536b3bc1ed434ef66
--- /dev/null
+++ b/fsl/wrappers/fslmaths.py
@@ -0,0 +1,155 @@
+#!/usr/bin/env python
+#
+# fslmaths.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+
+from fsl.utils.run import runfsl
+
+
+class fslmaths(object):
+    """Perform mathematical manipulation of images."""
+
+    def __init__(self, input):
+        """Constructor."""
+        self.inputImage = input
+        self.outputImage = None
+        self.operations = []
+
+    def abs(self):
+        """Absolute value."""
+        self.operations.append("-abs")
+        return self
+
+    def bin(self):
+        """Use (current image>0) to binarise."""
+        self.operations.append("-bin")
+        return self
+
+    def binv(self):
+        """Binarise and invert (binarisation and logical inversion)."""
+        self.operations.append("-binv")
+        return self
+
+    def recip(self):
+        """Reciprocal (1/current image)."""
+        self.operations.append("-recip")
+        return self
+
+    def Tmean(self):
+        """Mean across time."""
+        self.operations.append("-Tmean")
+        return self
+
+    def Tstd(self):
+        """Standard deviation across time."""
+        self.operations.append("-Tstd")
+        return self
+
+    def Tmin(self):
+        """Min across time."""
+        self.operations.append("-Tmin")
+        return self
+
+    def Tmax(self):
+        """Max across time."""
+        self.operations.append("-Tmax")
+        return self
+
+    def fillh(self):
+        """fill holes in a binary mask (holes are internal - i.e. do not touch
+        the edge of the FOV)."""
+        self.operations.append("-fillh")
+        return self
+
+    def ero(self, repeat=1):
+        """Erode by zeroing non-zero voxels when zero voxels in kernel."""
+        for i in range(repeat):
+            self.operations.append("-ero")
+        return self
+
+    def dilM(self, repeat=1):
+        """Mean Dilation of non-zero voxels."""
+        for i in range(repeat):
+            self.operations.append("-dilM")
+        return self
+
+    def dilF(self, repeat=1):
+        """Maximum filtering of all voxels."""
+        for i in range(repeat):
+            self.operations.append("-dilF")
+        return self
+
+    def add(self, input):
+        """Add input to current image."""
+        self.operations.append("-add {0}".format(input))
+        return self
+
+    def sub(self, input):
+        """Subtract input from current image."""
+        self.operations.append("-sub {0}".format(input))
+        return self
+
+    def mul(self, input):
+        """Multiply current image by input."""
+        self.operations.append("-mul {0}".format(input))
+        return self
+
+    def div(self, input):
+        """Divide current image by input."""
+        self.operations.append("-div {0}".format(input))
+        return self
+
+    def mas(self, image):
+        """Use image (>0) to mask current image."""
+        self.operations.append("-mas {0}".format(image))
+        return self
+
+    def rem(self, input):
+        """Divide current image by following input and take remainder."""
+        self.operations.append("-rem {0}".format(input))
+        return self
+
+    def thr(self, input):
+        """use input number to threshold current image (zero < input)."""
+        self.operations.append("-thr {0}".format(input))
+        return self
+
+    def uthr(self, input):
+        """use input number to upper-threshold current image (zero
+        anything above the number)."""
+        self.operations.append("-uthr {0}".format(input))
+        return self
+
+    def inm(self, input):
+        """Intensity normalisation (per 3D volume mean)"""
+        self.operations.append("-inm {0}".format(input))
+        return self
+
+    def bptf(self, hp_sigma, lp_sigma):
+        """Bandpass temporal filtering; nonlinear highpass and Gaussian linear
+        lowpass (with sigmas in volumes, not seconds); set either sigma<0 to
+        skip that filter."""
+        self.operations.append("-bptf {0} {1}".format(hp_sigma, lp_sigma))
+        return self
+
+    def toString(self):
+        """Generate fslmaths command string."""
+        cmd = "fslmaths {0} ".format(self.inputImage)
+        for oper in self.operations:
+            cmd = cmd + oper + " "
+        cmd = cmd + self.outputImage
+        return cmd
+
+    def run(self, output=None):
+        """Save output of operations to image."""
+        if output:
+            self.outputImage = output
+        else:
+            self.outputImage = self.inputImage
+
+        runfsl(self.toString())
+
+        return self.outputImage
diff --git a/fsl/wrappers/fugue.py b/fsl/wrappers/fugue.py
new file mode 100644
index 0000000000000000000000000000000000000000..298493beb08f457c9b8d78f3572de793ec707a1a
--- /dev/null
+++ b/fsl/wrappers/fugue.py
@@ -0,0 +1,42 @@
+#!/usr/bin/env python
+#
+# fugue.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+"""
+"""
+
+
+import fsl.utils.run as run
+
+
+def fugue(**kwargs):
+    """FMRIB's Utility for Geometric Unwarping of EPIs."""
+
+    cmd = "fugue"
+
+    if kwargs.pop('unmaskshift', False):
+        cmd += " --unmaskshift"
+    if kwargs.pop('despike', False):
+        cmd += " --despike"
+    if kwargs.pop('unmaskfmap', False):
+        cmd += " --unmaskfmap"
+
+    cmd += ' '.join(['--{}={}'.format(k, v) for k, v in kwargs.items()])
+
+    return run.runfsl(cmd)
+
+
+def sigloss(input, output, te=None, slicedir=None, mask=None):
+    """Estimate signal loss from a field map (in rad/s)."""
+    cmd = "sigloss -i {0} -s {1}".format(input, output)
+
+    if te is not None:
+        cmd += " --te={0}".format(te)
+    if slicedir is not None:
+        cmd += " --slicedir={0}".format(slicedir)
+    if mask is not None:
+        cmd += " --mask={0}".format(mask)
+
+    return run.runfsl(cmd)
diff --git a/fsl/wrappers/melodic.py b/fsl/wrappers/melodic.py
new file mode 100644
index 0000000000000000000000000000000000000000..82bf28844540bb10d4051d48ac36eaac351d85f4
--- /dev/null
+++ b/fsl/wrappers/melodic.py
@@ -0,0 +1,57 @@
+#!/usr/bin/env python
+#
+# melodic.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+
+import fsl.utils.run        as run
+import fsl.utils.assertions as asrt
+
+
+def melodic(input, outdir, dim=None, tr=None, mmthresh=None, report=True,
+            prefix=None, nomask=False, updatemask=False, nobet=False,
+            mask=None):
+    """Multivariate Exploratory Linear Optimised ICA."""
+    asrt.assertIsNifti(input)
+
+    cmd = "melodic -i {0} -v --Oall --outdir={1}".format(input, outdir)
+
+    if mmthresh:
+        cmd += " --mmthresh={0}".format(mmthresh)
+    if dim:
+        cmd += " -d -{0}".format(dim)
+    if report:
+        cmd += " --report"
+    if tr:
+        cmd += " --tr={0}".format(tr)
+    if nomask:
+        cmd += " --nomask"
+    if updatemask:
+        cmd += " --update_mask"
+    if nobet:
+        cmd += " --nobet"
+    if prefix:
+        cmd = prefix + " " + cmd
+    if mask is not None:
+        cmd += " --mask={0}".format(mask)
+    return run.runfsl(cmd)
+
+
+def fsl_regfilt(infile, outfile, mix, ics):
+    """Data de-noising by regression.
+
+    Data de-noising by regressing out part of a design matrix
+    using simple OLS regression on 4D images
+    """
+    asrt.assertIsNifti(infile, outfile)
+
+    icstr = '"'
+    for i in range(0, len(ics) - 1):
+        icstr = icstr + '{0},'.format(ics[i] + 1)
+    icstr = icstr + '{0}"'.format(ics[-1] + 1)
+
+    cmd = "fsl_regfilt -i {0} -o {1} -d {2} -f {3}".format(infile, outfile,
+                                                           mix, icstr)
+    return run.runfsl(cmd)
diff --git a/fsl/wrappers/misc.py b/fsl/wrappers/misc.py
new file mode 100644
index 0000000000000000000000000000000000000000..94df91ac704a42a2f384ceb10f4760ca411458c4
--- /dev/null
+++ b/fsl/wrappers/misc.py
@@ -0,0 +1,81 @@
+#!/usr/bin/env python
+#
+# misc.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+
+import fsl.utils.run        as run
+import fsl.utils.assertions as asrt
+
+
+def fslreorient2std(input, output):
+    """reorient to match the approx. orientation of the standard (MNI152)."""
+
+    asrt.assertIsNifti(input, output)
+    asrt.assertFileExists(input)
+
+    cmd = 'fslreorient2std {0} {1}'.format(input, output)
+    return run.runfsl(cmd)
+
+
+def fslroi(input, output, xmin=None, xsize=None, ymin=None, ysize=None,
+           zmin=None, zsize=None, tmin=None, tsize=None):
+    """Extract region of interest (ROI) from an image."""
+    assert ((tmin is not None and tsize is not None) or
+            (xmin is not None and xsize is not None and
+            ymin is not None and ysize is not None and
+            zmin is not None and zsize is not None)), \
+        "either time min/size or x/y/z min/size must be provided"
+
+    cmd = "fslroi {0} {1}".format(input, output)
+
+    if xmin is not None:
+        cmd += " {0} {1} {2} {3} {4} {5}".format(xmin, xsize, ymin, ysize,
+                                                 zmin, zsize)
+    if tmin is not None:
+        cmd += " {0} {1}".format(tmin, tsize)
+
+    return run.runfsl(cmd)
+
+
+def slicer(input, input2=None, label=None, lut=None, intensity=None,
+           edgethreshold=None, x=None, y=None, z=None):
+
+    cmd = "slicer {0}".format(input)
+
+    if input2 is not None:
+        cmd += " {0}".format(input2)
+    if label is not None:
+        cmd += " -L {0}".format(label)
+    if lut is not None:
+        cmd += " -l {0}".format(lut)
+    if intensity is not None:
+        cmd += " -i {0} {1}".format(intensity[0], intensity[1])
+    if edgethreshold is not None:
+        cmd += " -e {0}".format(edgethreshold)
+    if x is not None:
+        cmd += " -x {0} {1}".format(x[0], x[1])
+    if y is not None:
+        cmd += " -y {0} {1}".format(y[0], y[1])
+    if z is not None:
+        cmd += " -z {0} {1}".format(z[0], z[1])
+
+    return run.runfsl(cmd)
+
+
+def cluster(infile, thresh, oindex=None, no_table=False):
+    """
+    Form clusters, report information about clusters and/or perform
+    cluster-based inference.
+    """
+    cmd = "cluster --in={0} --thresh={1}".format(infile, thresh)
+
+    if oindex is not None:
+        cmd += " -o {0}".format(oindex)
+
+    if no_table:
+        cmd += " --no_table"
+
+    return run.runfsl(cmd)
diff --git a/fsl/wrappers/wrapperutils.py b/fsl/wrappers/wrapperutils.py
new file mode 100644
index 0000000000000000000000000000000000000000..629845bf2df840db9614af71849bef16c70c7621
--- /dev/null
+++ b/fsl/wrappers/wrapperutils.py
@@ -0,0 +1,230 @@
+#!/usr/bin/env python
+#
+# wrapperutils.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+
+import os.path as op
+import            os
+import            inspect
+import            tempfile
+import            warnings
+import            functools
+import            collections
+
+import            six
+import nibabel as nib
+
+import fsl.utils.tempdir as tempdir
+import fsl.data.image    as fslimage
+
+
+class _BooleanFlag(object):
+    def __init__(self, show):
+        self.show = show
+    def __eq__(self, other):
+        return type(other) == type(self) and self.show == other.show
+
+
+SHOW_IF_TRUE = _BooleanFlag(True)
+HIDE_IF_TRUE = _BooleanFlag(False)
+
+
+def applyArgStyle(style, argmap=None, valmap=None, **kwargs):
+
+    def fmtarg(arg, style):
+        if   style in ('-',  '-='):  arg =  '-{}'.format(arg)
+        elif style in ('--', '--='): arg = '--{}'.format(arg)
+        return arg
+
+    def fmtval(val, style=None):
+        if     isinstance(val, collections.Sequence) and \
+           not isinstance(val, six.string_types):
+            return ' '.join([str(v) for v in val])
+        else:
+            return str(val)
+
+    if style not in ('-', '--', '-=', '--='):
+        raise ValueError('Invalid style: {}'.format(style))
+
+    args = []
+
+    for k, v in kwargs.items():
+
+        k    = argmap.get(k, k)
+        mapv = valmap.get(k, fmtval(v, style))
+        k    = fmtarg(k, style)
+
+        if mapv in (SHOW_IF_TRUE, HIDE_IF_TRUE):
+            if v == mapv.show:
+                args.append(k)
+        elif '=' in style:
+            args.append('{}={}'.format(k, mapv))
+        else:
+            args.extend((k, mapv))
+
+    return args
+
+
+def required(*reqargs):
+    """Decorator which makes sure that all specified keyword arguments are
+    present before calling the decorated function.
+    """
+    def decorator(func):
+        def wrapper(**kwargs):
+            for reqarg in reqargs:
+                assert reqarg in kwargs
+            return func(**kwargs)
+        return wrapper
+    return decorator
+
+
+def argsToKwargs(func, args):
+    """Given a function, and a sequence of positional arguments destined
+    for that function, converts the positional arguments into a dict
+    of keyword arguments. Used by the :class:`_FileOrImage` class.
+    """
+    # getfullargspec is the only way to get the names
+    # of positional arguments in Python 2.x. It is
+    # deprecated in python 3.5, but not in python 3.6.
+    with warnings.catch_warnings():
+        warnings.filterwarnings('ignore', category=DeprecationWarning)
+        spec = inspect.getfullargspec(func)
+
+    kwargs = collections.OrderedDict()
+    for name, val in zip(spec.args, args):
+        kwargs[name] = val
+
+    return kwargs
+
+
+RETURN = object()
+"""
+"""
+
+
+class _FileOrImage(object):
+    """
+
+    Inputs:
+      - In-memory nibabel images loaded from a file. The image is replaced with
+        its file name.
+
+      - In-memory nibabel images. The image is saved to a temporary file, and
+        replaced with the temporary file's name. The file is deleted after the
+        function has returned.
+
+    Outputs:
+      - File name:  The file name is passed straight through to the function.
+      - ``RETURN``: A temporary file name is passed to the function. After the
+        function has completed, the image is loaded into memory and the
+        temporary file is deleted. The image is returned from the function
+        call.
+    """
+
+
+    def __init__(self, *imgargs):
+        """
+        """
+        self.__imgargs = imgargs
+
+
+    def __call__(self, func):
+        """
+        """
+        return functools.partial(self.__wrapper, func)
+
+
+    def __wrapper(self, func, *args, **kwargs):
+        """
+        """
+
+        kwargs.update(argsToKwargs(func, args))
+
+        # Create a tempdir to store any temporary
+        # input/output images, but don't change
+        # into it, as file paths passed to the
+        # function may be relative.
+        with tempdir.tempdir(changeto=False) as td:
+
+            kwargs, infiles, outfiles = self.__prepareArgs(td, kwargs)
+
+            # Call the function
+            result  = func(**kwargs)
+
+            # Load the output images that
+            # were specified as RETURN
+            outimgs = []
+            for of in outfiles:
+
+                # output file didn't get created
+                if not op.exists(of):
+                    oi = None
+
+                # load the file, and create
+                # an in-memory copy (the file
+                # is going to get deleted)
+                else:
+                    oi = nib.load(of)
+                    oi = nib.nifti1.Nifti1Image(oi.get_data(), None, oi.header)
+
+                outimgs.append(oi)
+
+            return tuple([result] + outimgs)
+
+
+    def __prepareArgs(self, workdir, kwargs):
+        """
+        """
+
+        kwargs   = dict(kwargs)
+        infiles  = []
+        outfiles = []
+
+        for imgarg in self.__imgargs:
+
+            img = kwargs.get(imgarg, None)
+
+            # Not specified, nothing to do
+            if img is None:
+                continue
+
+            # This is an input image which has
+            # been specified as an in-memory
+            # nibabel image. if the image has
+            # a backing file, replace the image
+            # object with the file name.
+            # Otherwise, save the image out to
+            # a temporary file, and replace the
+            # image with the file name.
+            if isinstance(img, nib.nifti1.Nifti1Image):
+                imgfile = img.get_filename()
+
+                # in-memory image - we have
+                # to save it out to a file
+                if imgfile is None:
+
+                    hd, imgfile = tempfile.mkstemp(fslimage.defaultExt())
+
+                    os.close(hd)
+                    img.to_filename(imgfile)
+                    infiles.append(imgfile)
+
+                # replace the image with its
+                # file name
+                kwargs[img] = imgfile
+
+            # This is an output image, and the
+            # caller has requested that it be
+            # returned from the function call
+            # as an in-memory image.
+            if img == RETURN:
+                kwargs[imgarg] = '{}.nii.gz'.format(imgarg)
+                outfiles.append(imgarg)
+
+        return kwargs, infiles, outfiles
+
+
+fileOrImage = _FileOrImage