Volume rendering artefacts when looking along the z-plane

Question

I've been trying to get some WebGL volume rendering up and running. I've managed to tile 3D data to a 2D texture (tiling both laterally and along the RGB channels). And it now displays some synthetic data (a sphere).

When viewed along the x,y planes (i.e. looking face on at the tiles), everything seems ok.

But as soon as you start viewing rays that are cast entering the cube from a different side, there is a weird haze effect which blocks out the volume being rendered, like this

and when you look directly along the z planes (i.e. edge on to the tiles) then things get really messed up

Has anyone seen this before? The complete code is available here

Here is the business end of my shader - the part that looks samples the tiled array and does the ray marching.

varying vec3 worldSpaceCoords;
varying vec4 projectedCoords;

uniform sampler2D firstPassTexture, dataTexture; //i.e. tex and cubeTex
uniform float steps;
uniform float alphaCorrection;


const int MAX_STEPS = 512;

const vec3 dataShape = vec3(99, 93, 330); // really this should be passed not hard coded
const vec2 textureShape = vec2(1024, 1024);

vec3 getDatumColor(float datum) {
    vec3 color = vec3(255, 255, 255);
    if (datum == 9999.9999){
        color = vec3(255, 0, 0);
    }
    return color;
}

float getDatumAlpha(float datum) {
    return datum * alphaCorrection;
}

float getDatum(sampler2D tex, vec3 pos, vec3 dataShape, vec2 nTiles, float tilesPerLayer, vec2 tileDim, float thisTileN){
    float zTile = floor(thisTileN/tilesPerLayer);
    float yTile = floor((thisTileN - (zTile * tilesPerLayer)) / nTiles.x);
    float xTile = mod((thisTileN - (zTile * tilesPerLayer)), nTiles.x);

    vec2 thisPoint = vec2(xTile+pos.x, yTile+pos.y) * tileDim;

    float datum;
    if (zTile == 0.0){
        datum = texture2D(tex, thisPoint).r;
    }else if (zTile == 1.0){
        datum = texture2D(tex, thisPoint).g;
    }else if (zTile == 2.0){
        datum = texture2D(tex, thisPoint).b;
    }

    return datum;
}

float sampleAs3DTexture(sampler2D tex, vec3 pos, vec3 dataShape, vec2 texShape) {
    /* 
    A function to reference a 2D RGBA texture which contains tiles 3D array data.

    Tiling goes column, row, channel

    Args:
        * tex: texture of tiled data
        * pos: position of the datum
        * dataShape: the x,y,z shape of the data which has been tiled
        * texShape: the x,y dims of the tiles texture
    */

    vec2 fracNTiles = texShape.xy / dataShape.xy;
    vec2 nTiles = vec2(floor(fracNTiles.x), floor(fracNTiles.y));
    float tilesPerLayer = nTiles.x * nTiles.y;
    vec2 tileDim = vec2(1.0, 1.0) / fracNTiles;
    float thisTileN = floor((dataShape.z-1.0) * pos.z);
    float thisTileNp1 = min(thisTileN+1.0, dataShape.z);

    float datumN = getDatum(tex, pos, dataShape, nTiles, tilesPerLayer, tileDim, thisTileN);
    float datumNp1 = getDatum(tex, pos, dataShape, nTiles, tilesPerLayer, tileDim, thisTileNp1);

    float zDiff = mod((dataShape.z-1.0) * pos.z, 1.0);

    return ((1.0 - zDiff) * datumN) + (zDiff * datumNp1);
}

vec4 getRGBAfromDataTex(sampler2D tex, vec3 pos, vec3 dataShape, vec2 texShape){
    float datum = sampleAs3DTexture(tex, pos, dataShape, texShape);
    vec3 color = getDatumColor(datum);
    float alpha = getDatumAlpha(datum);

    return vec4(color.xyz, alpha);
}

// max 2d size is 4096 x 4096

void main( void ) {
    //Transform the coordinates it from [-1;1] to [0;1]
    vec2 firstPassTexCoord = vec2(((projectedCoords.x / projectedCoords.w) + 1.0 ) / 2.0,
                    ((projectedCoords.y / projectedCoords.w) + 1.0 ) / 2.0 );

    //The back position is the world space position stored in the texture.
    vec3 backPos = texture2D(firstPassTexture, firstPassTexCoord).xyz;

    //The front position is the world space position of the second render pass.
    vec3 frontPos = worldSpaceCoords;

    //The direction from the front position to back position.
    vec3 dir = backPos - frontPos;

    float rayLength = length(dir);

    //Calculate how long to increment in each step.
    float delta = 1.0 / steps;

    //The increment in each direction for each step.
    vec3 deltaDirection = normalize(dir) * delta;
    float deltaDirectionLength = length(deltaDirection);

    //Start the ray casting from the front position.
    vec3 currentPosition = frontPos;

    //The color accumulator.
    vec3 accumulatedColor = vec3(0.0);

    //The alpha value accumulated so far.
    float accumulatedAlpha = 0.0;

    //How long has the ray travelled so far.
    float accumulatedLength = 0.0;

    //vec4 dataSample;
    vec4 dataSample;

    float alphaSample;

    //Perform the ray marching iterations
    for(int i = 0; i < MAX_STEPS; i++){
        //Get the voxel intensity value from the 3D texture.    
        dataSample = getRGBAfromDataTex(dataTexture, currentPosition, dataShape, textureShape);

        //Perform the composition.
        accumulatedColor += (1.0 - accumulatedAlpha) * dataSample.xyz * dataSample.a;
        //accumulatedColor += dataSample;

        //Store the alpha accumulated so far.
        accumulatedAlpha += dataSample.a;

        //Advance the ray.
        currentPosition += deltaDirection;
        accumulatedLength += deltaDirectionLength;

        //If the length traversed is more than the ray length, or if the alpha accumulated reaches 1.0 then exit.
        if(accumulatedLength >= rayLength || accumulatedAlpha >= 1.0 )
            break;
    }
    vec4 fragColor = vec4(accumulatedColor.x, accumulatedColor.y, accumulatedColor.z, accumulatedAlpha);
    gl_FragColor = fragColor;
}

Thanks for any help everyone!