Binary scene voxelization using imageStore - problem

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;

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

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

texcoord.x = temp.z;
texcoord.y = temp.y;
texcoord.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)

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;

texcoord = temp;

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

texcoord.x = temp.z;
texcoord.y = temp.y;
texcoord.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);
}

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)

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);