Binary scene voxelization using imageStore - problem

Question

This problem has been haunting me for some time but I was busy with other things. I can't finish my GI technique (specifically occlusion step) without solving it.

I made a binary scene voxelization in a single pass. The result is stored in a RGBA32UI 3D texture. It is three way so I set three ortho matrices each for X, Y, Z direction. These matrices are centered in camera position and form frustum in cube shape. And in geometry shader choose axis based on triangle normal.

First, what works - without binary voxelization (but it's a waste of VRAM). Here is fragment shader code of voxelization pass:

//gridDim - gridDimensions in voxels
ivec3 temp = ivec3( gl_FragCoord.xy, gridDim.z * gl_FragCoord.z) ;
ivec3 texcoord[3];

//correction of texture coord based on chosen axis and ortho matrix 
//(f_axis - 0, 1 or 2 from geometry shader)
texcoord[0] = temp;

texcoord[1].x = int(gridDim.z) - temp.x;
texcoord[1].y = temp.z;
texcoord[1].z = temp.y;

texcoord[2].x = temp.z;
texcoord[2].y = temp.y;
texcoord[2].z = int(gridDim.z) - temp.x;

imageStore(geometryGrid, texcoord[f_axis], uvec4(1u));

Drawing pass is done using instancing with cubes. The voxels are quite big because more of them results in a big framerate drop:

//vpVoxels is ortho matrix, worldPosCenter - center of each cube from vertex shader
vec4 geomCoords = vpVoxels * vec4(worldPosCenter, 1.0);
geomCoords.xyz /= geomCoords.w;
geomCoords.xyz = geomCoords.xyz * 0.5 + 0.5;

ivec3 index = ivec3(geomCoords.xyz * gridDim);

uvec4 voxel2 = texelFetch(geometryTexture, index, 0);
if(voxel2.r == 0u)
    discard;

normalColor = vec4(normal0, 1.0);
fragColor = vec4(1.0, 0.0, 0.0, 1.0);

And it works:

But if I change it to the binary version:

//input - position 1-128, output: coded binary position as 1 in uvec4 format
    uvec4 encodeBinaryVoxel(int position)
    {
        // check for positioning to R channel
        position = max(position, 0);
        uint voxelId1 = (position > 32 || position == 0)? 0u : 1u << uint(position - 1);

        // check for positioning to G channel
        position = max(position - 32, 0);
        uint voxelId2 = (position > 32 || position == 0)? 0u : 1u << uint(position - 1);

        // check for positioning to B channel
        position = max(position - 32, 0);
        uint voxelId3 = (position > 32 || position == 0)? 0u : 1u << uint(position - 1);

        // check for positioning to A channel
        position = max(position - 32, 0);
        uint voxelId4 = (position > 32 || position == 0)? 0u : 1u << uint(position - 1);

        return uvec4(voxelId1, voxelId2, voxelId3, voxelId4);
    }

void main()
{
    ivec3 temp = ivec3( gl_FragCoord.xy, gridDim.z * gl_FragCoord.z) ;
    ivec3 texcoord[3];

    texcoord[0] = temp;

    texcoord[1].x = int(gridDim.z) - temp.x;
    texcoord[1].y = temp.z;
    texcoord[1].z = temp.y;

    texcoord[2].x = temp.z;
    texcoord[2].y = temp.y;
    texcoord[2].z = int(gridDim.z) - temp.x;

    //get z index of texture, each z slice has room for 128 bits (32*4)
    int indexZ = texcoord[f_axis].z / 128; 
    //get position in slice
    int bitZ = texcoord[f_axis].z % 128 + 1; 
    ivec3 index = ivec3(texcoord[f_axis].xy, indexZ);
    uvec4 current = imageLoad(geometryGrid, index);
    //encode position as binary uvec4
    uvec4 bit = encodeBinaryVoxel(bitZ); 
    //or operation with stored value
    imageStore(geometryGrid, index, bit | current); 
}

And read:

vec4 geomCoords = vpVoxels * vec4(worldPosCenter, 1.0);
geomCoords.xyz /= geomCoords.w;
geomCoords.xyz = geomCoords.xyz * 0.5 + 0.5;

ivec3 index = ivec3(geomCoords.xyz * gridDim);
int indexZ = index.z / 128;
int bitZ = index.z % 128 + 1;
uvec4 binZ = encodeBinaryVoxel(bitZ);

uvec4 voxel = texelFetch(geometryTexture, ivec3(index.xy, indexZ), 0);
uvec4 occupy = binZ & voxel;
if((occupy.r | occupy.g | occupy.b | occupy.a) == 0u)
    discard;

normalColor = vec4(normal0, 1.0);
fragColor = vec4(1.0, 0.0, 0.0, 1.0);

The result:

As you see in the binary version it flickers and I have failed to find out the logic behind this (yet). However, voxel positions seem to be fine. I can render these voxels also by raymarching through the texture and this flickers too.

EDIT

I've made a simplest voxelization using only vertex + fragment and one ortho matrix in X direction (then it is incorrect but whatever), and store value only in first slice (my gridDim is less than 128 so it fits) and it still flickers:

ivec3 index = ivec3( gl_FragCoord.xy, gridDim.z * gl_FragCoord.z) ;

uvec4 bit = encodeBinaryVoxel(index.z + 1);

uvec4 current = imageLoad(geometryGrid, ivec3(index.xy, 0));

imageStore(geometryGrid, ivec3(index.xy, 0), bit | current);

Answer 1

In the binary version you're doing an imageLoad and an imageStore on the same texture to modify it. This is a classic race condition: the GPU executes multiple threads of the shader in parallel, and there's no guarantee for those load-modify-store sequences to happen atomically. The operations from different threads can be interleaved, resulting in wrong data stored.

To fix it, you should use imageAtomicOr to combine the bits you want to write with the bits already in the buffer without a race condition. For more information, see this OpenGL wiki page.