
// [Eberly2014] GPGPU Programming for Games and Science
float fastacos(float x)
{
    float res = -0.156583 * abs(x) + PI*0.5;
    res *= sqrt(1.0 - abs(x));
    return x >= 0 ? res : PI - res;
}

float integrateArc(float h1, float h2) {
    float cosN = 1.0;//cos(n);
    float sinN = 0.0;//sin(n);
    return 0.25 * (-cos(2.0 * h1) + cosN + 2.0 * h1 * sinN - cos(2.0 * h2) + cosN + 2.0 * h2 * sinN);
}

float IntegrateArc(float h1, float h2, float n) {
    float cosN = cos(n);
    float sinN = sin(n);
    return 0.25 * (-cos(2.0 * h1 - n) + cosN + 2.0 * h1 * sinN - cos(2.0 * h2 - n) + cosN + 2.0 * h2 * sinN);
}

vec3 multiBounceApprox(float ao, vec3 albedo) {
	vec3 a = 2.0404*albedo - 0.3324;
	vec3 b = 4.7951*albedo - 0.6417;
	vec3 c = 2.7552*albedo + 0.6903;
	return a*ao*ao*ao - b*ao*ao + c*ao;
}

float intersectPlane(vec3 origin, vec3 direction, vec3 point, vec3 normal) {
    return clamp(dot(point - origin, normal) / dot(direction, normal), -1.0, 9991999.0);
}

vec3 getInitialHorizonDirection(ivec2 icoord, vec3 viewPos, vec3 normal, vec2 pixelStep) {
	vec3 dir = normalize(getDepthPoint(downsampleAmount*(icoord + pixelStep*0.5)/resolution, 1.0));
	float dist = intersectPlane(vec3(0.0), dir, viewPos, normal);
	vec3 hitPos = dir*dist;
	return normalize(hitPos - viewPos);
}

float searchMaximumHorizonAngle(ivec2 icoord, vec3 viewPos, vec3 normal, vec2 pixelStep, float dither) {

	 float maxRadius = length(viewPos);

	// horizonDirection = getInitialHorizonDirection(icoord, viewPos, normal, pixelStep);

	float thetaMin = 0.0;//PI*0.5; // Minimum theta is maximum horizon angle
	float thetaPrev = thetaMin;
	ivec2 sampleCoord = ivec2(icoord + dither*pixelStep);
	for(int i = 0; i < 3; i++, sampleCoord += ivec2(pixelStep)) {
		
		float sampleDepth = texelFetch(depthtex1, sampleCoord*downsampleAmount, 0).r;

		if(clamp(sampleCoord, ivec2(0), ivec2(resolution*downsampleFactor-1)) != sampleCoord || sampleDepth == 1.0) continue;

		vec3 samplePoint = getDepthPoint(downsampleAmount*sampleCoord/resolution,sampleDepth );

		float falloff = clamp((maxRadius - length(samplePoint - viewPos)) / maxRadius, 0.0, 1.0);

		vec3 dirToSamplePoint = normalize(samplePoint - viewPos);
		float thetaSample = (dot(dirToSamplePoint, normal));//fastacos


		if(thetaSample >= thetaMin) {
			// thetaMin = min(thetaMin, thetaSample);
			thetaMin = mix(thetaMin, thetaSample, falloff);
			// horizonDirection = dirToSamplePoint;
		} else {
			// thetaMin = mix(thetaMin, thetaSample, falloff);
		}
		// else thetaMin = mix(thetaPrev, thetaSample, 0.2);

		thetaPrev = thetaSample;
		// pixelStep *= 2.0;
	}
	return fastacos(thetaMin);
}




float calculateGTAO(ivec2 icoord, vec3 viewPos, vec3 normal, inout float seed, out vec3 bentNormal) {
	float sum = 0.0;
	bentNormal = vec3(0.0);

	const float stepSize = 12.0*length(resolution)/length(vec2(1920.0, 1080.0));//100.0/length(viewPos);//resolution.x / 64.0;

	const int N = 2;
	float sliceRot =  PI/float(N) * hash1(seed);// float((frameCounter % 16))/16.0;

	float dither = hash1(seed);

	vec3 viewDir = -normalize(viewPos);

	for(int sliceIdx = 0; sliceIdx < N; sliceIdx++) {
		
		vec2 sliceDir = vec2(cos(sliceRot), sin(sliceRot));

		// vec3 toDir = normalize(getDepthPoint(downsampleAmount*(icoord + sliceDir)/resolution, 1.0));
		// vec3 planeNormal = normalize(cross(viewDir, -toDir));
		// vec3 projectedNormal = normal - planeNormal * dot(normal, planeNormal);


		// vec3 projectedDir = normalize(normalize(toDir) + viewDir);
    	// float n = fastacos(dot(-projectedDir, normalize(projectedNormal))) - PI_HALF;

		vec3 bitang = normalize(cross(vec3(sliceDir, 0.0), viewDir));
		vec3 tang = cross(viewDir, bitang);
		vec3 nx = normal - bitang * dot(normal, bitang);

		float nnx		= sqrt(dot(nx, nx));
    	float invnnx	= 1.0 / (nnx + 1e-6);	
		float cosxi		= dot(nx, tang) * invnnx;
		float gamma		= fastacos(cosxi) - PI_HALF;
		float cos_gamma	= dot(nx, viewDir) * invnnx;
		float sin_gamma = -2.0 * cosxi;

		// vec3 dir1, dir2;
		float h1a = -searchMaximumHorizonAngle(icoord, viewPos, normal, sliceDir*stepSize, dither);
		float h2a = searchMaximumHorizonAngle(icoord, viewPos, normal, -sliceDir*stepSize, dither);
		 
		float theta1 = gamma + max(h1a - gamma, -PI_HALF);//, dir1);
		float theta2 = gamma + min(h2a - gamma, PI_HALF);//, dir2);

		float bentAngle = (theta1+theta2)*0.5;

		sum += (0.5 * cos_gamma + 0.5 * (theta1 + theta2) * sin_gamma - 0.25 * (cos(2.0 * theta1 - gamma) + cos(2.0 * theta2 - gamma)))*nnx;
		// sum += IntegrateArc(theta1, theta2, n)* (nnx);
		bentNormal += viewDir * cos(bentAngle) - tang * sin(bentAngle);

		sliceRot += PI/float(N);
		// sliceRot =  PI * hash1(seed);
	}
	bentNormal = normalize(normalize(bentNormal) - viewDir * 0.5);
	bentNormal = viewToWorld(bentNormal);
	return clamp01(sum   / float(N));
}