Tensor lighting implemented, but not working yet. Can't figure out what

I'm doing wrong.

fsl/fsleyes/gl/gl21/gltensor_funcs.py

@@ -104,7 +104,7 @@ def compileShaders(self):
                     'v1ValXform',      'v2ValXform', 'v3ValXform',
                     'l1ValXform',      'l2ValXform', 'l3ValXform',
                     'voxToDisplayMat', 'imageShape', 'resolution',
-                    'lighting']
+                    'lighting',        'lightDir']
 
     vertAtts     = ['voxel', 'vertex']
 
@@ -220,11 +220,21 @@ def preDraw(self):
     """
     gl.glUseProgram(self.shaders)
 
-    if self.displayOpts.lighting:
-        gl.glEnable(gl.GL_LIGHTING)
-        gl.glEnable(gl.GL_LIGHT0)
-    
-        # gl.glLight(gl.GL_LIGHT0, gl.GL_)
+    # Define the light direction in
+    # the world coordinate system
+    lightDir = np.array([1, 1, 1], dtype=np.float32)
+    lightDir[self.zax] = 3
+
+    # Transform it into the display
+    # coordinate system
+    mvMat     = gl.glGetFloatv(gl.GL_MODELVIEW_MATRIX)
+    lightDir  = np.dot(mvMat[:3, :3], lightDir)
+
+    # Normalise to unit length
+    ldLen     = np.sqrt(np.sum(lightDir ** 2))
+    lightDir /= ldLen
+
+    gl.glUniform3fv(self.lightDirPos, 1, lightDir)
     
     self.v1Texture.bindTexture(gl.GL_TEXTURE5)
     self.v2Texture.bindTexture(gl.GL_TEXTURE6)
@@ -236,6 +246,9 @@ def preDraw(self):
     gl.glEnableVertexAttribArray(self.voxelPos)
     gl.glEnableVertexAttribArray(self.vertexPos)
 
+    gl.glEnable(gl.GL_CULL_FACE)
+    gl.glCullFace(gl.GL_BACK)
+    
 
 def draw(self, zpos, xform=None):
 
@@ -296,10 +309,8 @@ def draw(self, zpos, xform=None):
 def postDraw(self):
     gl.glUseProgram(0)
 
-    if self.displayOpts.lighting: 
-        gl.glDisable(gl.GL_LIGHT0)
-        gl.glDisable(gl.GL_LIGHTING)
-        
+    gl.glDisable(gl.GL_CULL_FACE)
+    
     gl.glBindBuffer(gl.GL_ARRAY_BUFFER,         0)
     gl.glBindBuffer(gl.GL_ELEMENT_ARRAY_BUFFER, 0)
     

fsl/fsleyes/gl/gl21/gltensor_vert.glsl

@@ -27,6 +27,8 @@ uniform float resolution;
 
 uniform bool lighting;
 
+uniform vec3 lightDir;
+
 attribute vec3  voxel;
 attribute vec3  vertex;
 
@@ -63,36 +65,66 @@ void main(void) {
   // shift it by the voxel coordinates
   // to transform it in to the image
   // voxel coordinate system.
-  vec3 pos  = vertex;
+  //
+  // See these web pages for details
+  // on drawing ellipsoids:
+  //
+  //   - http://archive.gamedev.net/archive/reference/  \
+  //     programming/features/superquadric/index.html
+  //
+  //   - paulbourke.net/geometry/circlesphere/
+  vec3 pos     = vertex;
+  mat3 eigvecs = mat3(v1, v2, v3);
 
   // TODO scale by the range of l1/l2/l3
   pos.x *= l1 * 150;
   pos.y *= l2 * 150;
   pos.z *= l3 * 150;
-  
-  pos  = mat3(v1, v2, v3) * pos;
-  pos += voxel;
+  pos    = eigvecs * pos;
+  pos   += voxel;
 
   if (lighting) {
-    
+
+    // Calculate the normal
+    // vector for this vertex.
     vec3 norm = vertex;
-    norm.x   /= l2;
+    norm.x   /= l1;
     norm.y   /= l2;
     norm.z   /= l3;
-    norm      = mat3(v1, v2, v3) * pos;
-    
-    light     = vec3(1, 1, 1);
+
+    // The matrix made up of the eigenvectors
+    // should be orthonormal, so can be used
+    // to transform the vertex surface normals.
+    //
+    // [ orthonormal -> eigvecs = T(I(eigvecs)) ]
+    norm      = eigvecs * norm;
+
+    // GLSL 1.20 calculates the normal
+    // matrix for us; we will have to
+    // calculate it ourselves if we move
+    // to a more modern version of GLSL.
+    norm     = gl_Normal * norm;
+
+    norm = normalize(norm);
+
+    // Calculate the diffuse and
+    // ambient light components
+    float ambient = 0.1;
+    float diffuse = clamp(dot(norm, -lightDir), 0, 1);
+
+    diffuse = (diffuse + ambient);
+    light   = vec3(diffuse, diffuse, diffuse);
   }
   
   else
-    light     = vec3(1, 1, 1);
+    light = vec3(1, 1, 1);
 
   // Transform from voxels into display
   gl_Position = gl_ModelViewProjectionMatrix *
                 voxToDisplayMat              *
                 vec4(pos, 1);
 
-  fragVoxCoord     = voxel;
+  fragVoxCoord     = floor(voxel + 0.5);
   fragTexCoord     = texCoord;
   fragColourFactor = vec4(light, 1);
 }

fsl/fsleyes/gl/routines.py

@@ -83,9 +83,6 @@ def show2D(xax, yax, width, height, lo, hi):
     elif zax == 1:
         gl.glRotatef(270, 1, 0, 0)
 
-    gl.glDepthMask(gl.GL_FALSE)
-    gl.glDisable(gl.GL_DEPTH_TEST)
-
 
 def calculateSamplePoints(shape, resolution, xform, xax, yax, origin='centre'):
     """Calculates a uniform grid of points, in the display coordinate system
@@ -650,22 +647,23 @@ def unitSphere(res):
 
     :returns: A tuple comprising:
     
-              - a ``numpy.float32`` array of size ``((res + 1)**2, 3)``
+              - a ``numpy.float32`` array of size ``(res**2, 3)``
                 containing a set of ``(x, y, z)`` vertices which define
                 the ellipsoid surface.
      
-              - A ``numpy.uint32`` array of size ``(4 * res * res)``
+              - A ``numpy.uint32`` array of size ``(4 * (res - 1)**2)``
                 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
+
+    .. todo:: Generate indices to use with ``GL_TRIANGLES`` instead of
+              ``GL_QUADS``.
+    """
 
     # 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)
+    u = np.linspace(-np.pi / 2, np.pi / 2, res, dtype=np.float32)
+    v = np.linspace(-np.pi,     np.pi,     res, dtype=np.float32)
 
     cosu = np.cos(u)
     cosv = np.cos(v)
@@ -675,14 +673,14 @@ def unitSphere(res):
     cucv = np.outer(cosu, cosv)
     cusv = np.outer(cosu, sinv)
 
-    vertices = np.zeros(((res + 1) ** 2, 3), dtype=np.float32)
+    vertices = np.zeros((res ** 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)
+    vertices[:, 2] = sinu.repeat(res)
 
     # Generate a list of indices which join the
     # vertices so they can be used to draw the
@@ -695,12 +693,12 @@ def unitSphere(res):
     #  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)
+    nquads   = (res - 1) ** 2
+    quadIdxs = np.array([0, res, res + 1, 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)
+    indices += np.arange(nquads,  dtype=np.uint32).repeat(4)
+    indices += np.arange(res - 1, dtype=np.uint32).repeat(4 * (res - 1))
     
     return vertices, indices
 
@@ -711,23 +709,20 @@ def fullUnitSphere(res):
 
     :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).
+    :returns: A ``numpy.float32`` array of size ``(4 * (res - 1)**2, 3)``
+              containing the ``(x, y, z)`` vertices which can be used to draw
+              a unit sphere (using the ``GL_QUADS`` primitive type).
     """
 
-    ustep = np.pi       / res
-    vstep = np.pi * 2.0 / res
-
-    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)
+    u = np.linspace(-np.pi / 2, np.pi / 2, res, dtype=np.float32)
+    v = np.linspace(-np.pi,     np.pi,     res, dtype=np.float32)
 
     cosu = np.cos(u)
     cosv = np.cos(v)
     sinu = np.sin(u)
     sinv = np.sin(v) 
 
-    vertices = np.zeros((res * res * 4, 3), dtype=np.float32)
+    vertices = np.zeros(((res - 1) * (res - 1) * 4, 3), dtype=np.float32)
 
     cucv   = np.outer(cosu[:-1], cosv[:-1]).flatten()
     cusv   = np.outer(cosu[:-1], sinv[:-1]).flatten()
@@ -738,8 +733,8 @@ def fullUnitSphere(res):
     cucv1  = np.outer(cosu[:-1], cosv[1:]) .flatten()
     cusv1  = np.outer(cosu[:-1], sinv[1:]) .flatten()
     
-    su     = np.repeat(sinu[:-1], res)
-    s1u    = np.repeat(sinu[1:],  res)
+    su     = np.repeat(sinu[:-1], res - 1)
+    s1u    = np.repeat(sinu[1:],  res - 1)
 
     vertices.T[:,  ::4] = [cucv,   cusv,   su]
     vertices.T[:, 1::4] = [cu1cv,  cu1sv,  s1u]