diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 0347a5df351c952d247feb8efaa3416d70d62c36..912d98a07979aa01b3feb25ff2b4cf8cd355ac3e 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -12,6 +12,20 @@ Added
 
 * New tensor conversion routines in the :mod:`.dtifit` module (Michiel
   Cottaar).
+* New :func:`.makeWriteable` function which ensures that a ``numpy.array`` is
+  writeable, and creates a copy if necessary
+
+
+Changed
+^^^^^^^
+
+
+* The :class:`.GiftiMesh` class no longer creates copies of the mesh
+  vertex/index arrays. This means that, with ``numpy>=1.16`` these arrays
+  will be flagged as read-only.
+* The :class:`.Mesh` class handles vertex data sets requiring different
+  triangle unwinding orders, at the cost of potentially having to store
+  two copies of the mesh indices.
 
 
 Fixed
diff --git a/fsl/data/__init__.py b/fsl/data/__init__.py
index c1890b375baadbf3a78f1ac8437145040cab8570..3c5afc3f2d09c49b5c1e5be18d20b7b5eab0cc56 100644
--- a/fsl/data/__init__.py
+++ b/fsl/data/__init__.py
@@ -9,4 +9,5 @@ models, constants, and other data-like things used throughout ``fslpy``.
 """
 
 
-from .utils import guessType  # noqa
+from .utils import (guessType,  # noqa
+                    makeWriteable)
diff --git a/fsl/data/gifti.py b/fsl/data/gifti.py
index abd41f7d43e468136c08be7e418ec8cbdb0bda9f..aa8b20ed7b8cedf7a168c901458f0aa85e8d770e 100644
--- a/fsl/data/gifti.py
+++ b/fsl/data/gifti.py
@@ -211,8 +211,8 @@ def loadGiftiMesh(filename):
         raise ValueError('{}: GIFTI surface files must contain '
                          'at least one pointset array'.format(filename))
 
-    vertices = [np.array(ps.data) for ps in pointsets]
-    indices  = np.array(triangles[0].data)
+    vertices = [ps.data for ps in pointsets]
+    indices  = triangles[0].data
 
     if len(vdata) == 0: vdata = None
     else:               vdata = prepareGiftiVertexData(vdata, filename)
diff --git a/fsl/data/mesh.py b/fsl/data/mesh.py
index f5ca20adea3747ca53945755f61c85dc28c539d1..e581c96037108376d1ba551b46f32d9b8720c0df 100644
--- a/fsl/data/mesh.py
+++ b/fsl/data/mesh.py
@@ -98,6 +98,12 @@ class Mesh(notifier.Notifier, meta.Meta):
        selectedVertices
        vertexSets
 
+    .. note:: Internally the ``Mesh`` class may store two versions of the
+              triangles, with opposite unwinding orders. It keeps track of the
+              required triangle unwinding order for each vertex set, so that
+              the :meth:`indices` method will return the appropriate copy for
+              the currently selected vertex set.
+
 
     **Vertex data**
 
@@ -183,19 +189,17 @@ class Mesh(notifier.Notifier, meta.Meta):
 
         self.__name       = name
         self.__dataSource = dataSource
-        self.__indices    = np.asarray(indices).reshape((-1, 3))
-        self.__nvertices  = self.__indices.max() + 1
-
-        # This attribute is used to store
-        # the currently selected vertex set,
-        # used as a kety into all of the
-        # dictionaries below.
-        self.__selected = None
-
-        # Flag used to keep track of whether
-        # the triangle winding order has been
-        # "fixed" - see the addVertices method.
-        self.__fixed = False
+        self.__nvertices  = indices.max() + 1
+        self.__selected   = None
+
+        # We potentially store two copies of
+        # the indices, with opposite unwinding
+        # orders. The vindices dict stores refs
+        # to one or the other for each vertex
+        # set.
+        self.__indices      = np.asarray(indices).reshape((-1, 3))
+        self.__fixedIndices = None
+        self.__vindices     = collections.OrderedDict()
 
         # All of these are populated
         # in the addVertices method
@@ -287,7 +291,7 @@ class Mesh(notifier.Notifier, meta.Meta):
     @property
     def indices(self):
         """The ``(M, 3)`` triangles of this mesh. """
-        return self.__indices
+        return self.__vindices[self.__selected]
 
 
     @property
@@ -297,7 +301,7 @@ class Mesh(notifier.Notifier, meta.Meta):
         """
 
         selected = self.__selected
-        indices  = self.__indices
+        indices  = self.__vindices[selected]
         vertices = self.__vertices[selected]
         fnormals = self.__faceNormals.get(selected, None)
 
@@ -315,7 +319,7 @@ class Mesh(notifier.Notifier, meta.Meta):
         """
 
         selected = self.__selected
-        indices  = self.__indices
+        indices  = self.__vindices[selected]
         vertices = self.__vertices[selected]
         vnormals = self.__vertNormals.get(selected, None)
 
@@ -334,10 +338,8 @@ class Mesh(notifier.Notifier, meta.Meta):
 
             ``((xlow, ylow, zlow), (xhigh, yhigh, zhigh))``
         """
-
         lo = self.__loBounds[self.__selected]
         hi = self.__hiBounds[self.__selected]
-
         return lo, hi
 
 
@@ -411,26 +413,26 @@ class Mesh(notifier.Notifier, meta.Meta):
                                                     self.nvertices))
 
         self.__vertices[key] = vertices
+        self.__vindices[key] = self.__indices
         self.__loBounds[key] = lo
         self.__hiBounds[key] = hi
 
         if select:
             self.vertices = key
 
-        # indices already fixed?
-        if fixWinding and (not self.__fixed):
-            indices      = self.indices
-            normals      = self.normals
-            needsFix     = needsFixing(vertices, indices, normals, lo, hi)
-            self.__fixed = True
+        if fixWinding:
+            indices  = self.__indices
+            normals  = self.normals
+            needsFix = needsFixing(vertices, indices, normals, lo, hi)
 
             # See needsFixing documentation
             if needsFix:
 
-                indices[:, [1, 2]] = indices[:, [2, 1]]
+                if self.__fixedIndices is None:
+                    self.__fixedIndices = indices[:, [0, 2, 1]]
 
-                for k, fn in self.__faceNormals.items():
-                    self.__faceNormals[k] = fn * -1
+                self.__vindices[   key] = self.__fixedIndices
+                self.__faceNormals[key] = normals * -1
 
         return vertices
 
diff --git a/fsl/data/utils.py b/fsl/data/utils.py
index e0fae151e05c02d5ecfe6289e5e1779ee9276a61..d902dc18a2b77c753c7d43f09d70398ebb25b71a 100644
--- a/fsl/data/utils.py
+++ b/fsl/data/utils.py
@@ -10,7 +10,7 @@ with the data types defined in the :mod:`fsl.data` package.
 
 
 import os.path as op
-
+import numpy   as np
 
 def guessType(path):
     """A convenience function which, given the name of a file or directory,
@@ -77,3 +77,14 @@ def guessType(path):
     # Otherwise, I don't
     # know what to do
     return None, path
+
+
+def makeWriteable(array):
+    """Updates the given ``numpy.array`` so that it is writeable. If this
+    is not possible, a copy is created and returned.
+    """
+    try:
+        array.flags['WRITEABLE'] = True
+    except ValueError:
+        array = np.array(array)
+    return array
diff --git a/tests/test_fsl_data_utils.py b/tests/test_fsl_data_utils.py
index 3fd02a8f948d921fda54c0ca9481a83154d371fc..1ef030bd228dbc263a0b390329a646f162e186a3 100644
--- a/tests/test_fsl_data_utils.py
+++ b/tests/test_fsl_data_utils.py
@@ -9,6 +9,8 @@ import            shutil
 import            os
 import os.path as op
 
+import numpy   as np
+
 import fsl.utils.tempdir        as tempdir
 import fsl.data.utils           as dutils
 
@@ -105,3 +107,22 @@ def test_guessType():
         asrt('norecognise.txt', None)
         os.remove('norecognise')
         os.remove('norecognise.txt')
+
+
+def test_makeWriteable():
+    robuf = bytes(    b'\01\02\03\04')
+    wbuf  = bytearray(b'\01\02\03\04')
+
+    roarr = np.ndarray((4,), dtype=np.uint8, buffer=robuf)
+    warr  = np.ndarray((4,), dtype=np.uint8, buffer=wbuf)
+
+    warr.flags['WRITEABLE'] = False
+
+    rocopy = dutils.makeWriteable(roarr)
+    wcopy  = dutils.makeWriteable(warr)
+
+    assert rocopy.base is not roarr.base
+    assert wcopy .base is     warr .base
+
+    rocopy[1] = 100
+    wcopy[ 1] = 100
diff --git a/tests/test_mesh.py b/tests/test_mesh.py
index 8b97bfa1c180cd301af462a2f067aff02c04e4b1..20ff3cdb772878e1ece865f6fc98f7b813fa1824 100644
--- a/tests/test_mesh.py
+++ b/tests/test_mesh.py
@@ -417,3 +417,29 @@ def test_planeIntersection():
     assert lines.shape  == (0, 2, 3)
     assert faces.shape  == (0, )
     assert dists.shape  == (0, 2, 3)
+
+
+def test_mesh_different_winding_orders():
+
+    verts1    =  CUBE_VERTICES
+    verts2    = -CUBE_VERTICES
+    tris      =  CUBE_TRIANGLES_CCW
+    trisfixed =  CUBE_TRIANGLES_CW
+
+    mnofix = fslmesh.Mesh(tris)
+    mfix   = fslmesh.Mesh(tris)
+
+    mnofix.addVertices(verts1, key='v1', fixWinding=False)
+    mnofix.addVertices(verts2, key='v2', fixWinding=False)
+    mfix  .addVertices(verts1, key='v1', fixWinding=True)
+    mfix  .addVertices(verts2, key='v2', fixWinding=True)
+
+    mnofix.vertices = 'v1'
+    assert np.all(mnofix.indices == tris)
+    mnofix.vertices = 'v2'
+    assert np.all(mnofix.indices == tris)
+
+    mfix.vertices = 'v1'
+    assert np.all(mfix.indices == tris)
+    mfix.vertices = 'v2'
+    assert np.all(mfix.indices == trisfixed)