diff --git a/fsl/data/featanalysis.py b/fsl/data/featanalysis.py
index 39e40c9122c163d68adbfe0d308bd9146ef37306..e95ff8ae1317eeac8e8905b4c9e35aa78ee465bb 100644
--- a/fsl/data/featanalysis.py
+++ b/fsl/data/featanalysis.py
@@ -427,48 +427,49 @@ def loadClusterResults(featdir, settings, contrast):
         # whitespace, and discarding
         # empty lines
         lines = f.readlines()
-        lines = [l.strip() for l in lines]
-        lines = [l         for l in lines if l != '']
-
-        # the first line should contain column
-        # names, and each other line should
-        # contain the data for one cluster
-        colNames     = lines[0]
-        clusterLines = lines[1:]
-
-        # each line should be tab-separated
-        colNames     =  colNames.split('\t')
-        clusterLines = [cl      .split('\t') for cl in clusterLines]
-
-        # Turn each cluster line into a
-        # Cluster instance. An error will
-        # be raised if the columm names
-        # are unrecognised (i.e. not in
-        # the colmap above), or if the
-        # file is poorly formed.
-        clusters = [Cluster(**dict(zip(colNames, cl))) for cl in clusterLines]
-
-        # Make sure all coordinates are in voxels -
-        # for first level analyses, the coordXform
-        # will be an identity transform (the coords
-        # are already in voxels). But for higher
-        # level, the coords are in mm, and need to
-        # be transformed to voxels.
-        for c in clusters:
-
-            zmax    = [c.zmaxx,    c.zmaxy,    c.zmaxz]
-            zcog    = [c.zcogx,    c.zcogy,    c.zcogz]
-            copemax = [c.copemaxx, c.copemaxy, c.copemaxz]
-
-            zmax    = transform.transform([zmax],    coordXform)[0].round()
-            zcog    = transform.transform([zcog],    coordXform)[0].round()
-            copemax = transform.transform([copemax], coordXform)[0].round()
-
-            c.zmaxx,   c.zmaxy,    c.zmaxz    = zmax
-            c.zcogx,   c.zcogy,    c.zcogz    = zcog
-            c.copemax, c.copemaxy, c.copemaxz = copemax
-
-        return clusters
+
+    lines = [line.strip() for line in lines]
+    lines = [line         for line in lines if line != '']
+
+    # the first line should contain column
+    # names, and each other line should
+    # contain the data for one cluster
+    colNames     = lines[0]
+    clusterLines = lines[1:]
+
+    # each line should be tab-separated
+    colNames     =  colNames.split('\t')
+    clusterLines = [cl      .split('\t') for cl in clusterLines]
+
+    # Turn each cluster line into a
+    # Cluster instance. An error will
+    # be raised if the columm names
+    # are unrecognised (i.e. not in
+    # the colmap above), or if the
+    # file is poorly formed.
+    clusters = [Cluster(**dict(zip(colNames, cl))) for cl in clusterLines]
+
+    # Make sure all coordinates are in voxels -
+    # for first level analyses, the coordXform
+    # will be an identity transform (the coords
+    # are already in voxels). But for higher
+    # level, the coords are in mm, and need to
+    # be transformed to voxels.
+    for c in clusters:
+
+        zmax    = [c.zmaxx,    c.zmaxy,    c.zmaxz]
+        zcog    = [c.zcogx,    c.zcogy,    c.zcogz]
+        copemax = [c.copemaxx, c.copemaxy, c.copemaxz]
+
+        zmax    = transform.transform([zmax],    coordXform)[0].round()
+        zcog    = transform.transform([zcog],    coordXform)[0].round()
+        copemax = transform.transform([copemax], coordXform)[0].round()
+
+        c.zmaxx,   c.zmaxy,    c.zmaxz    = zmax
+        c.zcogx,   c.zcogy,    c.zcogz    = zcog
+        c.copemax, c.copemaxy, c.copemaxz = copemax
+
+    return clusters
 
 
 def getDataFile(featdir):
diff --git a/fsl/data/featdesign.py b/fsl/data/featdesign.py
index e28f038c6bb3dfb864689aff8eb99f5e15a5fdac..68e7e398ed6d066c506f0365482e07daeb1d21f9 100644
--- a/fsl/data/featdesign.py
+++ b/fsl/data/featdesign.py
@@ -169,30 +169,15 @@ class FEATFSFDesign(object):
         if level == 1: getEVs = getFirstLevelEVs
         else:          getEVs = getHigherLevelEVs
 
-        self.__settings = collections.OrderedDict(settings.items())
-        self.__design   = np.array(designMatrix)
-        self.__numEVs   = self.__design.shape[1]
-        self.__evs      = getEVs(featDir, self.__settings, self.__design)
+        self.__loadVoxEVs = loadVoxelwiseEVs
+        self.__settings   = collections.OrderedDict(settings.items())
+        self.__design     = np.array(designMatrix)
+        self.__numEVs     = self.__design.shape[1]
+        self.__evs        = getEVs(featDir, self.__settings, self.__design)
 
         if len(self.__evs) != self.__numEVs:
             raise FSFError('Number of EVs does not match design.mat')
 
-        # Load the voxelwise images now,
-        # so they're ready to be used by
-        # the getDesign method.
-        for ev in self.__evs:
-
-            if not isinstance(ev, (VoxelwiseEV, VoxelwiseConfoundEV)):
-                continue
-
-            ev.image = None
-
-            # The path to some voxelwise
-            # EVs may not be present -
-            # see the VoxelwisEV class.
-            if loadVoxelwiseEVs and (ev.filename is not None):
-                ev.image = fslimage.Image(ev.filename)
-
 
     def getEVs(self):
         """Returns a list containing the :class:`EV` instances that represent
@@ -224,7 +209,7 @@ class FEATFSFDesign(object):
             if not isinstance(ev, (VoxelwiseEV, VoxelwiseConfoundEV)):
                 continue
 
-            if ev.image is None:
+            if (not self.__loadVoxEVs) or (ev.filename is None):
                 log.warning('Voxel EV image missing '
                             'for ev {}'.format(ev.index))
                 continue
@@ -300,14 +285,15 @@ class VoxelwiseEV(NormalEV):
 
     ============ ======================================================
     ``filename`` Path to the image file containing the data for this EV
+    ``image``    Reference to the :class:`.Image` object
     ============ ======================================================
 
     .. note:: The file for voxelwise EVs in a higher level analysis are not
               copied into the FEAT directory, so if the user has removed them,
               or moved the .gfeat directory, the file path here will not be
               valid. Therefore, a ``VoxelwiseEV`` will test to see if the
-              file exists, and will set the ``filename`` attribute to ``None``
-              it it does not exist.
+              file exists, and will set the ``filename`` and ``image``
+              attributes to ``None`` it it does not exist.
     """
 
     def __init__(self, realIdx, origIdx, title, filename):
@@ -330,6 +316,27 @@ class VoxelwiseEV(NormalEV):
                         'exist: {}'.format(filename))
             self.filename = None
 
+        self.__image = None
+
+
+    def __del__(self):
+        """Clears any reference to the voxelwise EV image. """
+        self.__image = None
+
+
+    @property
+    def image(self):
+        """Returns the :class:`.Image` containing the voxelwise EV data. """
+
+        if self.filename is None:
+            return None
+
+        if self.__image is not None:
+            return self.__image
+
+        self.__image = fslimage.Image(self.filename, mmap=False)
+        return self.__image
+
 
 class ConfoundEV(EV):
     """Class representing a confound EV.