diff --git a/fsl/fsleyes/gl/gl21/gltensor_funcs.py b/fsl/fsleyes/gl/gl21/gltensor_funcs.py
index 4cc8667f4998fec890d668fd90c6ba780ec539f6..33419e57fabf84104af5c75a0078dc4624a76e8d 100644
--- a/fsl/fsleyes/gl/gl21/gltensor_funcs.py
+++ b/fsl/fsleyes/gl/gl21/gltensor_funcs.py
@@ -9,6 +9,7 @@ import logging
import numpy as np
+import numpy.linalg as npla
import OpenGL.GL as gl
import OpenGL.raw.GL._types as gltypes
import OpenGL.GL.ARB.instanced_arrays as arbia
@@ -104,7 +105,8 @@ def compileShaders(self):
'v1ValXform', 'v2ValXform', 'v3ValXform',
'l1ValXform', 'l2ValXform', 'l3ValXform',
'voxToDisplayMat', 'imageShape', 'resolution',
- 'lighting', 'lightDir']
+ 'lighting', 'lightPos', 'normalMatrix',
+ 'eigValNorm', 'zax']
vertAtts = ['voxel', 'vertex']
@@ -185,8 +187,12 @@ def updateShaderState(self):
lighting = 1 if opts.lighting else 0
useSpline = 0
+ l1 = self.image. L1()
+ eigValNorm = 0.5 / abs(l1.data).max()
+
gl.glUniform3fv(self.imageShapePos, 1, imageShape)
gl.glUniform1f( self.resolutionPos, resolution)
+ gl.glUniform1f( self.eigValNormPos, eigValNorm)
gl.glUniform1f( self.lightingPos, lighting)
gl.glUniform1f( self.modThresholdPos, modThreshold)
gl.glUniform1f( self.useSplinePos, useSpline)
@@ -220,21 +226,26 @@ def preDraw(self):
"""
gl.glUseProgram(self.shaders)
- # Define the light direction in
+ # Define the light position in
# the world coordinate system
- lightDir = np.array([1, 1, 1], dtype=np.float32)
- lightDir[self.zax] = 3
+ lightPos = np.array([1, 1, -1], dtype=np.float32)
+
+ lightPos[self.zax] *= 3
- # Transform it into the display
- # coordinate system
- mvMat = gl.glGetFloatv(gl.GL_MODELVIEW_MATRIX)
- lightDir = np.dot(mvMat[:3, :3], lightDir)
+ # Transform the light position into
+ # the display coordinate system,
+ # and normalise to unit length
+ mvMat = gl.glGetFloatv(gl.GL_MODELVIEW_MATRIX)[:3, :3]
+ lightPos = np.dot(mvMat, lightPos)
+ lightPos /= np.sqrt(np.sum(lightPos ** 2))
- # Normalise to unit length
- ldLen = np.sqrt(np.sum(lightDir ** 2))
- lightDir /= ldLen
+ # Calculate a transformation matrix for
+ # normal vectors - T(I(MV matrix))
+ normalMatrix = npla.inv(mvMat).T
- gl.glUniform3fv(self.lightDirPos, 1, lightDir)
+ gl.glUniform1f( self.zaxPos, self.zax)
+ gl.glUniform3fv( self.lightPosPos, 1, lightPos)
+ gl.glUniformMatrix3fv(self.normalMatrixPos, 1, False, normalMatrix)
self.v1Texture.bindTexture(gl.GL_TEXTURE5)
self.v2Texture.bindTexture(gl.GL_TEXTURE6)
diff --git a/fsl/fsleyes/gl/gl21/gltensor_vert.glsl b/fsl/fsleyes/gl/gl21/gltensor_vert.glsl
index 3a7994be53a6a3a0440c29c457a9b36a24a3579d..3775bcab1ccb59bcfa33d476acd0af11f7fb64e2 100644
--- a/fsl/fsleyes/gl/gl21/gltensor_vert.glsl
+++ b/fsl/fsleyes/gl/gl21/gltensor_vert.glsl
@@ -5,6 +5,11 @@
*/
#version 120
+/*
+ * Textures containing the eigenvectors (v1, v2, v3)
+ * and eigenvalues (l1, l2, l3) of the diffusion
+ * tensor matrix.
+ */
uniform sampler3D v1Texture;
uniform sampler3D v2Texture;
uniform sampler3D v3Texture;
@@ -12,6 +17,11 @@ uniform sampler3D l1Texture;
uniform sampler3D l2Texture;
uniform sampler3D l3Texture;
+/*
+ * Transforms (scales/offsets) for transforming from
+ * data in the above textures to their original data
+ * range.
+ */
uniform mat4 v1ValXform;
uniform mat4 v2ValXform;
uniform mat4 v3ValXform;
@@ -19,22 +29,83 @@ uniform mat4 l1ValXform;
uniform mat4 l2ValXform;
uniform mat4 l3ValXform;
+/*
+ * Transformation matrix which transforms voxel
+ * coordinates into the display coordinate system.
+ */
uniform mat4 voxToDisplayMat;
+/*
+ * Transformation matrix which transforms normal
+ * vectors. This should be set to the transpose
+ * of the inverse of the model-view matrix (see
+ * http://www.scratchapixel.com/lessons/\
+ * mathematics-physics-for-computer-graphics/\
+ * geometry/transforming-normals for a good
+ * explanation).
+ */
+uniform mat3 normalMatrix;
+
+/*
+ * Image shape (x, y, z).
+ */
uniform vec3 imageShape;
-uniform float resolution;
+/*
+ * Scaling factor for eigenvalues - this controls
+ * the maximum size that any tensor ellipsoid is
+ * drawn. Set this to:
+ *
+ * 0.5 / max(abs(l1))
+ *
+ * to make the maximum ellipsoid circumference the
+ * size of one voxel.
+ */
+uniform float eigValNorm;
+/*
+ * Enable/disable a simple directional lighting model.
+ */
uniform bool lighting;
-uniform vec3 lightDir;
+/*
+ * The display coordinate system axis which is the depth
+ * axis. Required for some horrible lighting related hackiness.
+ */
+uniform float zax;
+
+/*
+ * Position of the directional light - must be specified
+ * in the view coordinate system..
+ */
+uniform vec3 lightPos;
+
+/*
+ * Voxel corresponding to the current vertex.
+ */
+attribute vec3 voxel;
-attribute vec3 voxel;
-attribute vec3 vertex;
+/*
+ * The current vertex on a unit sphere. The vertex
+ * will be transformed into an ellipsoid using the
+ * tensor matrix eigen-decomposition.
+ */
+attribute vec3 vertex;
+/*
+ * Voxel coordinate passed through to the fragment shader.
+ */
varying vec3 fragVoxCoord;
+
+/*
+ * Texture coordinate passed through to the fragment shader.
+ */
varying vec3 fragTexCoord;
+/*
+ * Multiplicative colour factor passed through to the
+ * fragment shader, used for lighting.
+ */
varying vec4 fragColourFactor;
@@ -42,7 +113,6 @@ void main(void) {
// Lookup the tensor parameters from the textures
vec3 texCoord = (voxel + 0.5) / imageShape;
- vec3 light;
vec3 v1 = texture3D(v1Texture, texCoord).xyz;
vec3 v2 = texture3D(v2Texture, texCoord).xyz;
@@ -76,54 +146,103 @@ void main(void) {
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;
+ // Scale the ellipsoid by a constant
+ // factor (calculated in gltensor_funcs)
+ pos.x *= l1 * eigValNorm;
+ pos.y *= l2 * eigValNorm;
+ pos.z *= l3 * eigValNorm;
pos = eigvecs * pos;
pos += voxel;
+
+ // Apply lighting if it is enabled
+ vec3 light;
if (lighting) {
// Calculate the normal
// vector for this vertex.
vec3 norm = vertex;
- norm.x /= l1;
- norm.y /= l2;
- norm.z /= l3;
- // The matrix made up of the eigenvectors
- // should be orthonormal, so can be used
- // to transform the vertex surface normals.
+ // Divide by the ellipsoid radii
+ // to get the normal of a point
+ // on the surface
+ norm.x /= l1;
+ norm.y /= l2;
+ norm.z /= l3;
+
+ // Transform the normal vectors by eigvecs
+ // to rotate them into the voxel coordinate
+ // system, and then transform them by the
+ // normal matrix to put them into display
+ // space.
+ //
+ // The eigvecs matrix should be orthonormal,
+ // so can be used to transform the vertex
+ // surface normals. The normalMatrix is
+ // calculated for us in gltensor_funcs.
//
// [ 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 = eigvecs * norm;
+ norm = normalMatrix * norm;
norm = normalize(norm);
+ // I honestly have no idea why this is necessary.
+ // There is some weird interaction going on in the
+ // transformation of the normal vectors from unit
+ // sphere space to ellipsoid space (through the
+ // eigvecs matrix), and then on to display space
+ // (through the normalMatrix).
+ //
+ // If our depth axis (zax) is y or z (1 or 2),
+ // and the eigvecs matrix has a negative determinant,
+ // we need to flip the z value of the normal vector.
+ // But if zax is 0, and the eigvecs matrix has a
+ // positive determinant, we need to flip the y value
+ // of the normal vector.
+ //
+ // If we don't do this, the lighting on the ellipsoids
+ // appears to be coming from behind them, instead of in
+ // front. I don't know what the hell is going on here.
+ float det = eigvecs[0][0] * eigvecs[1][1] * eigvecs[2][2] +
+ eigvecs[0][1] * eigvecs[1][2] * eigvecs[2][0] +
+ eigvecs[1][0] * eigvecs[2][1] * eigvecs[0][2] -
+ eigvecs[2][0] * eigvecs[1][1] * eigvecs[0][2] -
+ eigvecs[0][0] * eigvecs[1][2] * eigvecs[2][1] -
+ eigvecs[0][1] * eigvecs[1][0] * eigvecs[2][2];
+
+ // Avert your eyes.
+ if (det < 0) {
+ if (zax == 0)
+ norm.y = -norm.y;
+ }
+ else if (zax > 0) {
+ norm.z = -norm.z;
+ }
+
// Calculate the diffuse and
- // ambient light components
- float ambient = 0.1;
- float diffuse = clamp(dot(norm, -lightDir), 0, 1);
+ // ambient light components
+ float ambient = 0.3;
+ float angle = dot(norm, -lightPos);
+ float diffuse = max(sign(angle) * angle * angle, 0);
- diffuse = (diffuse + ambient);
+ diffuse = diffuse + ambient;
light = vec3(diffuse, diffuse, diffuse);
}
-
+
+ // If lighting is not enabled, the
+ // fragment colour is not modified.
else
light = vec3(1, 1, 1);
- // Transform from voxels into display
+ // Transform the vertex from the
+ // voxel coordinate system into
+ // the display coordinate system.
gl_Position = gl_ModelViewProjectionMatrix *
voxToDisplayMat *
vec4(pos, 1);
+ // Send the voxel and texture coordinates, and
+ // the colour scaling factor to the fragment shader.
fragVoxCoord = floor(voxel + 0.5);
fragTexCoord = texCoord;
fragColourFactor = vec4(light, 1);
diff --git a/fsl/fsleyes/gl/routines.py b/fsl/fsleyes/gl/routines.py
index 72ed37a4bedc7a3e26c982593f5585910f575ad4..f3d7ec1b3755ec8e48278f29ed912ec2d6c6cbdd 100644
--- a/fsl/fsleyes/gl/routines.py
+++ b/fsl/fsleyes/gl/routines.py
@@ -76,10 +76,10 @@ def show2D(xax, yax, width, height, lo, hi):
# TODO There's got to be a more generic way
# to perform this rotation. This will break
# if I add functionality allowing the user
- # to specifty the x/y axes on initialisation.
+ # to specifty the x/y axes on initialisation.
if zax == 0:
- gl.glRotatef(-90, 1, 0, 0)
- gl.glRotatef(-90, 0, 0, 1)
+ gl.glRotatef(270, 1, 0, 0)
+ gl.glRotatef(270, 0, 0, 1)
elif zax == 1:
gl.glRotatef(270, 1, 0, 0)