diff --git a/fsl/fsleyes/displaycontext/tensoropts.py b/fsl/fsleyes/displaycontext/tensoropts.py
index f8455697db7028cc88d3c1f3ca4181d19e78a3f4..e1ab15c25680d1908e197036d07961a69d2120ed 100644
--- a/fsl/fsleyes/displaycontext/tensoropts.py
+++ b/fsl/fsleyes/displaycontext/tensoropts.py
@@ -29,6 +29,7 @@ class TensorOpts(vectoropts.VectorOpts):
 
     # Tensor ellipsoid resolution
     tensorResolution = props.Int(default=10)
+
     
     def __init__(self, *args, **kwargs):
         
diff --git a/fsl/fsleyes/gl/gl21/gltensor_funcs.py b/fsl/fsleyes/gl/gl21/gltensor_funcs.py
index 08a66538971583f6b7cdc261bdd687229feff610..babd552002c44025d342f61a8c49bda8cd9554e0 100644
--- a/fsl/fsleyes/gl/gl21/gltensor_funcs.py
+++ b/fsl/fsleyes/gl/gl21/gltensor_funcs.py
@@ -12,6 +12,7 @@ import numpy                          as np
 import OpenGL.GL                      as gl
 import OpenGL.raw.GL._types           as gltypes
 import OpenGL.GL.ARB.instanced_arrays as arbia
+import OpenGL.GL.ARB.draw_instanced   as arbdi
 
 import fsl.utils.transform      as transform
 import fsl.fsleyes.gl.resources as glresources
@@ -63,6 +64,7 @@ def init(self):
         setattr(self, '{}Texture'.format(name), tex)
 
     self.vertexBuffer = gl.glGenBuffers(1)
+    self.indexBuffer  = gl.glGenBuffers(1)
     self.voxelBuffer  = gl.glGenBuffers(1)
 
     self.shaders = None
@@ -73,6 +75,7 @@ def init(self):
 
 def destroy(self):
     gl.glDeleteBuffers(1, gltypes.GLuint(self.vertexBuffer))
+    gl.glDeleteBuffers(1, gltypes.GLuint(self.indexBuffer))
     gl.glDeleteBuffers(1, gltypes.GLuint(self.voxelBuffer))
     gl.glDeleteProgram(self.shaders)
     
@@ -189,14 +192,23 @@ def updateShaderState(self):
     # shader will transform these vertices
     # into the tensor ellipsoid for each
     # voxel.
-    vertices = glroutines.unitSphere(resolution).ravel('C')
- 
-    self.nVertices = len(vertices) / 3
+    vertices, indices = glroutines.unitSphere(resolution)
+    
+    self.nVertices = len(indices)
+    vertices       = vertices.ravel('C')
+
+    gl.glUseProgram(0)    
 
     gl.glBindBuffer(gl.GL_ARRAY_BUFFER, self.vertexBuffer)
     gl.glBufferData(
         gl.GL_ARRAY_BUFFER, vertices.nbytes, vertices, gl.GL_STATIC_DRAW)
-    gl.glUseProgram(0)
+
+    gl.glBindBuffer(gl.GL_ELEMENT_ARRAY_BUFFER, self.indexBuffer)
+    gl.glBufferData(
+        gl.GL_ELEMENT_ARRAY_BUFFER, indices.nbytes, indices, gl.GL_STATIC_DRAW)
+
+    gl.glBindBuffer(gl.GL_ARRAY_BUFFER,         0)
+    gl.glBindBuffer(gl.GL_ELEMENT_ARRAY_BUFFER, 0)
 
 
 def preDraw(self):
@@ -239,20 +251,23 @@ def draw(self, zpos, xform=None):
 
     voxels  = transform.transform(voxels, d2vMat)
     nVoxels = len(voxels)
+    
+    voxels  = np.array(voxels, dtype=np.float32).ravel('C')
 
-    voxels = np.array(voxels, dtype=np.float32).ravel('C')
-
+    # Copy the voxel coordinates to the voxel buffer
     gl.glBindBuffer(gl.GL_ARRAY_BUFFER, self.voxelBuffer)
     gl.glBufferData(
         gl.GL_ARRAY_BUFFER, voxels.nbytes, voxels, gl.GL_STATIC_DRAW)
     gl.glVertexAttribPointer(
         self.voxelPos, 3, gl.GL_FLOAT, gl.GL_FALSE, 0, None)
+
+    # Use one set of voxel coordinates for every sphere drawn
     arbia.glVertexAttribDivisorARB(self.voxelPos, 1)
 
     # Bind the vertex buffer
     gl.glBindBuffer(gl.GL_ARRAY_BUFFER, self.vertexBuffer)
     gl.glVertexAttribPointer(
-        self.vertexPos, 3, gl.GL_FLOAT, gl.GL_FALSE, 0, None) 
+        self.vertexPos, 3, gl.GL_FLOAT, gl.GL_FALSE, 0, None)
 
     if xform is None: xform = v2dMat
     else:             xform = transform.concat(v2dMat, xform)
@@ -260,12 +275,17 @@ def draw(self, zpos, xform=None):
     xform = np.array(xform, dtype=np.float32).ravel('C') 
     gl.glUniformMatrix4fv(self.voxToDisplayMatPos, 1, False, xform)
 
-    gl.glDrawArraysInstanced(gl.GL_QUADS, 0, self.nVertices, nVoxels)
+    # And the vertex index buffer
+    gl.glBindBuffer(gl.GL_ELEMENT_ARRAY_BUFFER, self.indexBuffer)
+
+    arbdi.glDrawElementsInstancedARB(
+        gl.GL_QUADS, self.nVertices, gl.GL_UNSIGNED_INT, None, nVoxels)
 
 
 def postDraw(self):
     gl.glUseProgram(0)
-    gl.glBindBuffer(gl.GL_ARRAY_BUFFER, 0)
+    gl.glBindBuffer(gl.GL_ARRAY_BUFFER,         0)
+    gl.glBindBuffer(gl.GL_ELEMENT_ARRAY_BUFFER, 0)
     self.v1Texture.unbindTexture()
     self.v2Texture.unbindTexture()
     self.v3Texture.unbindTexture()
diff --git a/fsl/fsleyes/gl/routines.py b/fsl/fsleyes/gl/routines.py
index 61bc7738c08c3aa430fd9cfa76756b76864a92dd..2ec6374e9ca8ddb818c66b3f6d9601353371558c 100644
--- a/fsl/fsleyes/gl/routines.py
+++ b/fsl/fsleyes/gl/routines.py
@@ -648,14 +648,74 @@ def unitSphere(res):
 
     :arg res: Resolution - the number of angles to sample.
 
+    :returns: A tuple comprising:
+    
+              - a ``numpy.float32`` array of size ``((res + 1)**2, 3)``
+                containing a set of ``(x, y, z)`` vertices which define
+                the ellipsoid surface.
+     
+              - A ``numpy.uint32`` array of size ``(4 * res * res)``
+                containing a list of indices into the vertex array,
+                defining a vertex ordering that can be used to draw
+                the ellipsoid using the OpenGL ``GL_QUAD`` primitive type.
+    """
+
+    ustep = np.pi       / res
+    vstep = np.pi * 2.0 / res
+
+    # All angles to be sampled
+    u = np.linspace(-np.pi / 2, np.pi / 2 + ustep, res + 1, dtype=np.float32)
+    v = np.linspace(-np.pi,     np.pi     + vstep, res + 1, dtype=np.float32)
+
+    cosu = np.cos(u)
+    cosv = np.cos(v)
+    sinu = np.sin(u)
+    sinv = np.sin(v) 
+
+    cucv = np.outer(cosu, cosv)
+    cusv = np.outer(cosu, sinv)
+
+    vertices = np.zeros(((res + 1) ** 2, 3), dtype=np.float32)
+
+    # All x coordinates are of the form cos(u) * cos(v),
+    # y coordinates are of the form cos(u) * sin(v),
+    # and z coordinates of the form sin(u).
+    vertices[:, 0] = cucv.flatten()
+    vertices[:, 1] = cusv.flatten()
+    vertices[:, 2] = sinu.repeat(res + 1)
+
+    # Generate a list of indices which join the
+    # vertices so they can be used to draw the
+    # sphere as GL_QUADs.
+    # 
+    # The vertex locations for each quad follow
+    # this pattern:
+    # 
+    #  1. (u,         v)
+    #  2. (u + ustep, v)
+    #  3. (u + ustep, v + vstep)
+    #  4. (u,         v + vstep)
+    nquads   = res * res
+    quadIdxs = np.array([0, res + 1, res + 2, 1], dtype=np.uint32)
+
+    indices  = np.tile(quadIdxs, nquads)
+    indices += np.arange(nquads, dtype=np.uint32).repeat(4)
+    indices += np.arange(res,    dtype=np.uint32).repeat(4 * res)
+    
+    return vertices, indices
+
+
+def fullUnitSphere(res):
+    """Generates a unit sphere in the same way as :func:`unitSphere`, but
+    returns all vertices, instead of the unique vertices and an index array.
+
+    :arg res: Resolution - the number of angles to sample.
+
     :returns: A ``numpy.float32`` array of size ``(4 * res**2, 3)`` containing
               the ``(x, y, z)`` vertices which can be used to draw a unit
               sphere (using the ``GL_QUAD`` primitive type).
     """
 
-    # TODO Eliminate duplicate vertices,
-    #      and return an index array
-
     ustep = np.pi       / res
     vstep = np.pi * 2.0 / res