#version 430


in vec4 posG;
in vec3 normalG;

in float triangleIdG;

layout(location = 0) out vec4 frag;
// layout(location = 1) out vec4 frag2;

#define PI 3.1415926

uniform float g_time;


uniform mat4 modelViewMatrix;
uniform mat4 modelViewInvMatrix;
uniform mat4 viewInvMatrix;
uniform mat4 viewMatrix;
uniform mat4 projectionMatrix;
uniform mat4 projectionInvMatrix;

uniform float g_windowWidth = 1280.0;
uniform float g_windowHeight = 720.0;

vec4 rotateXY(vec4 p, float a) {
  vec4 r = p;
  r.x = cos(a)*p.x - sin(a)*p.y;
  r.y = sin(a)*p.x + cos(a)*p.y;
  return r;
}

vec3 rotateXY3(vec3 p, float a) {
  return rotateXY(vec4(p, 0.0), a).xyz;
}

// google glsl rand gave this, thanks and credit flies to
// http://stackoverflow.com/questions/4200224/random-noise-functions-for-glsl
float rand(vec2 co){
return fract(sin(dot(co.xy ,vec2(12.9898,78.233))) * 43758.5453);
}

const float zFar = 1.0;
const float zNear = 0.0;

float getPointDist(float z) {
    float clipA = zFar / (zFar - zNear);
    float clipB = zFar*zNear / (zNear - zFar);
    return clipB/(z-clipA);
}
float getPointZ(float d) {
    float clipA = zFar / (zFar - zNear);
    float clipB = zFar*zNear / (zNear - zFar);
    return (clipB + d*clipA)/d;
}

vec4 CalcEyeFromWindow(in vec3 windowSpace) {
    vec3 ndcPos;
    vec4 viewport = vec4(0.0, 0.0, g_windowWidth, g_windowHeight);
    ndcPos.xy = ((2.0 * windowSpace.xy) - (2.0 * viewport.xy)) / (viewport.zw) - 1;
//    ndcPos.z = (2.0 * windowSpace.z - zNear - zFar) / (zFar - zNear);
    ndcPos.z = (2.0 * windowSpace.z - zNear - zFar) / (zFar - zNear);
    vec4 clipPos;
    clipPos.w = projectionMatrix[3][2]/(ndcPos.z-(projectionMatrix[2][2]/projectionMatrix[2][3]));
    clipPos.xyz = ndcPos * clipPos.w;
    return projectionInvMatrix * clipPos;
}

float atanSafe(float y, float x) {
    float ret=0.0;
    if (x!=0.0) {
        if (x>0.0) {
            ret=atan(y/x);
        } else	{
            ret=atan(y/x)+3.141592;
        }
    } else {
        if (y>=0.0) {
            ret=0.5*3.141592;
        } else {
            ret=-0.5*3.141592;
        }
    }
    return ret;
}

// from pouet raymarching thread by las of mercury
float perlin(vec3 p) {
    vec3 i = floor(p);
    vec4 a = dot(i, vec3(1., 57., 21.)) + vec4(0., 57., 21., 78.);
    vec3 f = cos((p-i)*3.14159265)*(-.5)+.5;
    a = mix(sin(cos(a)*a),sin(cos(1.+a)*(1.+a)), f.x);
    a.xy = mix(a.xz, a.yw, f.y);
    return mix(a.x, a.y, f.z);
}

float turb(vec3 c) {
        float r=0.0;
        float s=0.5;
        c*=1.0;
        for (int i=0;i<8;i++) {
                //if (i<2) continue;
                r+=s*(perlin(c)*1.0+0.0);
                c*=3.0;
                s*=0.55;
        }
        return r;
}

layout(binding=0) uniform sampler2D texVoxelPos;
layout(binding=1) uniform sampler2D texVoxelNor;


layout(binding=0, offset=0) uniform atomic_uint ac;

layout(binding=0, r32i) uniform iimage2D voxelHP;
layout(binding=1, r32i) uniform iimage2D voxelNP;
layout(binding=2, rgba32f) uniform image2D voxPosComp;
layout(binding=3, rgba32f) uniform image2D voxNorComp;
layout(binding=4, rgba32f) uniform image2D voxVel;

uniform float g_voxelDim;

void main() {

  vec3 normal = normalG;
  vec4 diffuse = vec4(1.0);

  ivec2 mp = ivec2(gl_FragCoord);
  vec4 s = imageLoad(voxelHP, mp);
  if (s.r < 0.5) {
    return;
  }

  float vd2 = g_voxelDim*g_voxelDim;
  vec3 pos3D = vec3(0.0);
  vec2 pd2 = vec2(mp)/vd2;
  pos3D.xy = fract(pd2)*2.0-vec2(1.0);

  vec2 tk2 = floor(pd2);
  pos3D.z = ((tk2.x+tk2.y*g_voxelDim)/vd2)*2.0-1.0;

// debug
  frag.rgb = vec3(0.0);
  frag.rgb = vec3(pos3D); // +vec3(1.0);

// real deal
  //pos3D = vec3(1.0);
  uint vd3 = int(vd2*g_voxelDim);
  uint counter = atomicCounterIncrement(ac);
  uint cy = counter/vd3;
  uint cx = counter-uint(cy*vd3);
  ivec2 acp = ivec2(cx,cy);

  uint hp = uint(imageLoad(voxelHP, mp).r);

  vec3 vp = vec3(0.0);
  vec3 vn = vec3(0.0);
  float count = 0.0;
  while (hp > 0) {
      uint vcy = hp/vd3;
      uint vcx = hp-uint(vcy*vd3);
      ivec2 vcp = ivec2(vcx,vcy);
      vp += texelFetch(texVoxelPos, vcp, 0).xyz;
      vn += texelFetch(texVoxelNor, vcp, 0).xyz;
      count += 1.0;
      hp = uint(imageLoad(voxelNP, vcp).r);
  }
  if (count > 0.0) {
      vp /= count;
      vn /= count;
  }
  imageStore(voxVel, mp, vec4(vn, 0.0));

  imageStore(voxPosComp, acp, vec4(pos3D+vp*2.0/(g_voxelDim*g_voxelDim), 0.0));
  imageStore(voxNorComp, acp, vec4(vn, 0.0));


  return;

}

