#define PI 3.14159265
#define T [numthreads(16, 16, 1)]

RWTexture2D<float4> out0:register(u0);
RWTexture2D<float4> out1:register(u1);

Texture2D tex0: register(t0);
Texture2D tex1: register(t1);
Texture2D tex2: register(t2);
Texture2D tex3: register(t3);

SamplerState sampler0 : register(s0);

static uint rndSeed = 0;

uint hash(uint x) {
	return x += x << 10, x ^= x >> 6, x += x << 3, x ^= x >> 11, x += x << 15, x;
}

float rnd() {
	return asfloat((rndSeed = hash(rndSeed) & 0x007fffff) | 0x3f800000) - 1;
}

float srnd() {
	return asfloat((rndSeed = hash(rndSeed) & 0x007fffff) | 0x40000000) - 3;
}

float rnd(uint s) {
	return asfloat((s & 0x007fffff) | 0x3f800000) - 1;
}

float srnd(uint s) {
	return asfloat((s & 0x007fffff) | 0x40000000) - 3;
}

float expf(float x) {
	return asfloat(int(12102203 * x + 1064866805));
}

float2 mod(float2 x, float y) {
	return x - y * floor(x / y);
}

float pulse(float center, float width, float x) {
	float t = abs(x - center);
	if (t > width) return 0;
	t /= width;
	return 1 - t*t*(3 - 2 * t);
}

float fhex(float2 p, float s) {
	return p = abs(p), max(p.x + p.y*0.57735, p.y*1.1547) - s;
}

//---
// Bokeh
//---
static const float bkRadius = 32;
static const float bkFocus = 1.0;
static const float bkAperature = 8;
static const float bkBorderStrength = 0;

float bkCircleOfConfusion(float depth) {
	return clamp(bkAperature * abs(depth - bkFocus) / bkFocus, 1, bkRadius);
}

T void csBokehPreprocess(int3 id:SV_DispatchThreadID) {
	float2 resolution;
	out0.GetDimensions(resolution.x, resolution.y);

	float2 tc = (id.xy) / resolution;

	// Compute and premultiply CoC
	float4 s = tex0.SampleLevel(sampler0, tc, 0);
	float coc = bkCircleOfConfusion(s.w); 

	// Alternative CoC computations
	//float coc = bkCircleOfConfusion(tex0.SampleLevel(sampler0, tc, 2).w);
	
	//float4 depth = tex0.GatherAlpha(sampler0, tc);
	//float coc = (
	//		bkCircleOfConfusion(depth.x) +
	//		bkCircleOfConfusion(depth.y) +
	//		bkCircleOfConfusion(depth.z) +
	//		bkCircleOfConfusion(depth.w)
	//	) * 0.25;

	out0[id.xy] = float4(s.xyz, 1) * coc;
}

float4 bkSample(float2 tc, float2 d, Texture2D tex) {
	// Fetch reference CoC
	float referenceCoc = tex.SampleLevel(sampler0, tc, 0).w;

	float sigma = 0;
	float4 a = 0;
	for (float i = 0.5; i < bkRadius; ++i) {
		// Fetch sample float4(RGB,1)*coc
		float4 s = tex.SampleLevel(sampler0, d * i + tc, 0);

		// Compute and accumulate weighting factor
		float weight = saturate(referenceCoc - i) * (saturate(s.w - i) +  bkBorderStrength * step(bkRadius,i+1));
		sigma += weight;

		// Accumulate weighted contribution
		a += s * weight;
	}

	return a / sigma;
}

T void csBokehFirst(int3 id:SV_DispatchThreadID) {
	float2 resolution;
	out0.GetDimensions(resolution.x, resolution.y);

	float2 tc = (id.xy + 0.5) / resolution;

	float4 sA = bkSample(tc, float2(        0,   1) / resolution, tex0);
	float4 sB = bkSample(tc, float2(-.8660254, -.5) / resolution, tex0);

	out0[id.xy] = sA;
	out1[id.xy] = (sA + sB) * 0.5;
}

T void csBokehSecond(int3 id:SV_DispatchThreadID) {
	float2 resolution;
	out0.GetDimensions(resolution.x, resolution.y);

	float2 tc = (id.xy + 0.5) / resolution;

	float4 sA = bkSample(tc, float2(-.8660254, -.5) / resolution, tex0);
	float4 sB = bkSample(tc, float2(+.8660254, -.5) / resolution, tex1);

	out0[id.xy] = (sA += sB * 2) / sA.w;
}

T void csScatterBokeh(int3 id:SV_DispatchThreadID) {
	
}

//---
// Radial and circumferential blur
//---
static const float rcfMipmapFactor = 1.5;
static const float rcfStepFactor = 5;
static const float rcfStrength = 100;
static const float rcfCircumFerentialStrength = 1.0;
static const float rcfRadius = 5;
static const float2 rcfCenter = {0.5, 0.5};

float rcfRadialFalloff(float radius, float amount) {
	return radius * (1 + amount*radius*0.01);
}

T void csRadialCircumferentialBlur(int3 id:SV_DispatchThreadID) {
	float2 resolution;
	out0.GetDimensions(resolution.x, resolution.y);

	float2 tc = (id.xy + 0.5) / resolution;

	// Direction from rcfCenter to the current tc, aspect ratio corrected
	float2 d = (tc - rcfCenter) * float2(resolution.x/resolution.y, 1);

	// Compute convolution step size
	float delta = length(d);
	float distortedDelta = rcfRadialFalloff(delta, rcfStrength);
	float stepSize = abs(delta - distortedDelta) * rcfStepFactor;

	float radius = rcfRadius;
	// Circumferential blur
	if (rcfCircumFerentialStrength > 0) {
		d = float2(d.y, -d.x);
		stepSize *= rcfCircumFerentialStrength;
		radius = floor(rcfRadius * rcfCircumFerentialStrength);
	}

	// Early exit
	if (radius <= 0.5 || stepSize <= 1/resolution.x) {
		out0[id.xy] = tex0.SampleLevel(sampler0, tc, 0);
		return;
	}

	// Adjust direction
	d = normalize(d) * stepSize / resolution;

	// Convolution
	float mipLevel = sqrt(stepSize) * rcfMipmapFactor;
	float sigma = radius/2;
	float twoSigmaSquared = 2 * sigma * sigma;
	float4 a = 0;
	for (float i = -radius; i <= radius; ++i) {
		float weight = exp(-i*i/twoSigmaSquared);
		a +=  tex0.SampleLevel(sampler0, i * d + tc, mipLevel) * weight;
	}

	out0[id.xy] = a / sqrt(twoSigmaSquared * PI);
}

//---
// Small lens dirt
//---


//---
// Lens dirt
//---
static const float ldMinRadius = 16;
static const float ldMaxRadius = 32;
static const float ldCellWidth = 128;
static const float ldBorderWidthHalf = 1.5;
static const float ldSeed = 2;
static const float ldNumLayers = 4;

T void csLensDirt(int3 id:SV_DispatchThreadID) {
	float2 tc = id.xy + .5;

	rndSeed = ldSeed;
	
	float4 a = 0;
	for (int i = 0; i < ldNumLayers; ++i) {
		// Modulo repeat, find cell
		float2 p = mod(tc, ldCellWidth) - ldCellWidth / 2;
		float2 c = floor(tc / ldCellWidth);

		// Backup seed
		uint seedBackup = rndSeed;

		// Modifiy seed with cell coordinates
		rndSeed += c.x + c.y*ldCellWidth + i*ldCellWidth*ldCellWidth;

		// Determine radius and jitter
		float r = lerp(ldMinRadius, ldMaxRadius, rnd());
		p += float2(srnd(), srnd()) * (ldCellWidth / 2 - r);

		// Shapes
		// Hexagon
		//p = abs(p);
		//float d = max(p.y + p.x*0.57735, p.x*1.1547);

		// Circle
		//float d = length(p);

		// Pentagon
		float phi = 0.145;
		float2 n[] = {
			{ sin(0*PI/5 + phi), cos(0*PI/5 + phi) },
			{ sin(2*PI/5 + phi), cos(2*PI/5 + phi) },
			{ sin(4*PI/5 + phi), cos(4*PI/5 + phi) },
			{ sin(6*PI/5 + phi), cos(6*PI/5 + phi) },
			{ sin(8*PI/5 + phi), cos(8*PI/5 + phi) }
		};
		// Difference between circumcircle and incircle radii
		float delta = r - (1 + sqrt(5)) / 4 * r;
		float d = max(max(max(max(dot(p, n[0]), dot(p, n[1])), dot(p, n[2])), dot(p, n[3])), dot(p, n[4]));
		// Mix between sphere and pentagon
		d = lerp(length(p), d+ delta, 0.5);

		// Random color
		a += (float4(rnd(), rnd(), rnd(), 1) + 0.5) * (smoothstep(r, r - ldBorderWidthHalf * 2, d) + smoothstep(ldBorderWidthHalf, 0, abs(d - r + ldBorderWidthHalf * 2)));
		//a += float4(c, 0, 1)*0.1;

		// Restore seed
		rndSeed = seedBackup;

		// Displace grid
		tc += ldCellWidth / ldNumLayers;
	}

	out0[id.xy] = float4(a.xyz / ldNumLayers, 1);
}

//---
// Blur
//---
static const float bRadius = 160;
static const float bSigma = bRadius / 3;
static const float bTwoSigmaSquared = 2 * bSigma * bSigma;

void bBlur(int3 id, float2 direction) {
	float2 resolution;
	out0.GetDimensions(resolution.x, resolution.y);

	float2 sourceResolution;
	tex0.GetDimensions(sourceResolution.x, sourceResolution.y);

	float mipLevel = log2(sourceResolution.x / resolution.x);

	float2 tc = (id.xy + 0.5) / resolution;
	float2 d = direction / resolution;

	float4 a = 0;
	for (float i = -bRadius; i <= bRadius; ++i) {
		float weight = expf(-i*i/bTwoSigmaSquared);
		a += tex0.SampleLevel(sampler0, i * d + tc, mipLevel) * weight;
	}

	out0[id.xy] = a / sqrt(bTwoSigmaSquared * PI);
}

T void csBlurH(int3 id:SV_DispatchThreadID) { bBlur(id, float2(1,0)); }
T void csBlurV(int3 id:SV_DispatchThreadID) { bBlur(id, float2(0,1)); }

//---
// Streaks
//---
static const float stRadius = 160;

T void csStreaks(int3 id:SV_DispatchThreadID) {
	float2 resolution;
	out0.GetDimensions(resolution.x, resolution.y);
	
	float2 sourceResolution;
	tex0.GetDimensions(sourceResolution.x, sourceResolution.y);

	float mipLevel = log2(sourceResolution.x / resolution.x);

	float2 tc = (id.xy + 0.5) / resolution;
	float2 d = float2(1, 0) / resolution;

	// Convolution with parabolic tent filter
	float4 a = 0;
	for (float i = -stRadius+1; i < stRadius; ++i) {
		float weight = (stRadius - abs(i)) * (stRadius - abs(i));
		a += tex0.SampleLevel(sampler0, i * d + tc, mipLevel) * weight;
	}

	out0[id.xy] = a * 1.5 / (stRadius * stRadius * stRadius);
}

//---
// Merge
//---
T void csMerge(int3 id:SV_DispatchThreadID) {

	out0[id.xy] = 1;
}


/*
#define ST_NUMTHREADS_X 16
#define ST_NUMTHREADS_Y 16

#define ST_STRIDE (ST_NUMTHREADS_X + stRadius * 2)
groupshared float4 sData[ST_STRIDE * ST_NUMTHREADS_Y];
#define ST_ACCESS(X,Y) sData[clamp(X+stRadius, 0, ST_STRIDE-1) + ST_STRIDE * Y]
[numthreads(ST_NUMTHREADS_X, ST_NUMTHREADS_Y, 1)]
void csStreaksTest(int3 id:SV_DispatchThreadID, int3 gid : SV_GroupThreadID) {
	float2 resolution;

	tex0.GetDimensions(resolution.x, resolution.y);
	float2 d = normalize(float2(1, 0)) / resolution;

	float2 tc = (id.xy + 0.5) / resolution;
	ST_ACCESS(gid.x, gid.y) = tex0.SampleLevel(sampler0, tc, 0);

	int N = (stRadius+ST_NUMTHREADS_X-1) / ST_NUMTHREADS_X;
	int M = stRadius / N;
	if (gid.x < stRadius)
	for (int i = 1 ; i <= N; ++i) {
		ST_ACCESS(gid.x + M*i, gid.y) = tex0.SampleLevel(sampler0, tc + d * M*i, 0);
		ST_ACCESS(gid.x - M*i, gid.y) = tex0.SampleLevel(sampler0, tc - d * M*i, 0);
	}
		
	GroupMemoryBarrier();

	float4 color = 0;
	int j = ST_STRIDE * gid.y;
	////int k = gid.x;
	//for (float i = 0; i < stRadius; ++i) {
	//	int ii = (i+k) % ST_STRIDE;
	//	float weight = (stRadius - abs(ii)) * (stRadius - abs(ii));
	//	//color += ST_ACCESS(gid.x+i, gid.y) * weight;
	//	color += sData[j++] * weight;
	//}

	for (float i = -stRadius+1; i < stRadius; ++i) {
		float weight = (stRadius - abs(i)) * (stRadius - abs(i));
		//color += ST_ACCESS(gid.x+i, gid.y) * weight;
		color += sData[j++] * weight;
	}
	
	out0[id.xy] = color * 1.5 / (stRadius * stRadius * stRadius);
}

*/
//---
// Test 1
//---
static float3 d, p, n;
static float h, i, j, k;
static float d0, d1, d2, d3, dOut;


float ff(float3 p){
	d0 = 30 - length(p);
	d1 = length(p - float3(0, 1, 0)) - 0.02;
	d2 = p.y;
	float3 q = p;
	q.x = abs(q.x) - 2;
	q.y -= 1;
	q.y = abs(q.y) - 0.5;
	//
	q.z = frac(q.z);
	d3 = length(q) - 0.15;
	return dOut = min(min(d0, d1), min(d2, d3));
}

T void csRaytrace(int3 id:SV_DispatchThreadID) {
	//out0[id.xy] = float4(id.y % 16, id.x % 16,  0, 1)*0.1;
	p = float3(0, .5, 4);
	float3 q = p;

	d = mul(float3x3(.002, 0, -1.279, 0, -.002, 0.719, 0, 0, 1), float3(id.xy, 1));
	d = normalize(float3(d.xy, -1));

	for (h = .001, i = k = 1; k + 2 > i + h; i += h)
		k = ff(p), p += d*k;

	float4 color = 0.25;
	if (dOut == d0)
		color = float4(0, 0, 0, 1); // background
	if (dOut == d1)
		color = float4(1, 1, 0, 1)*100; // sphere
	if (dOut == d2)
		color = float4(0.5, 0, 0, 1); // ground
	if (dOut == d3)
		color = float4(1, 1, 0, 1); // sphere


	n = normalize(float3(ff(float3(h, 0, 0) + p), ff(float3(0, h, 0) + p), ff(float3(0, 0, h) + p)) - ff(p));

	for (h = .1, i = j = k = 1; i < 6; i += 1)
		j -= (i*h - ff(p + n*i*h)) / (k *= 2);

	float depth = dot(float3(0, 0, -1), d*length(p - q));
	//float coc =  bkCircleOfConfusion(depth);
	float4 c =  float4(color.xyz * (j*0.5 + (n*0.5 + 0.5)*0.5), depth);
		//c.w = 1;
	out0[id.xy] = c;
}
//---- RM-


T void cs2(int3 id:SV_DispatchThreadID) {
	out0[id.xy] = float4(1, id.xy / float2(1280, 720), 1);
	if (abs(id.y-359.5) <= 4 && abs(id.x-639.5) <= 100)
		out0[id.xy] = 1;
}


T void cs__(int3 id:SV_DispatchThreadID) {
	out0[id.xy] = float4((id.xy + 0.5) / float2(1280, 720), 0, 1);

}

//---
// Test 2
//---
static const float t = 100.0;

#define VP float3(.947556,4.117568,-1.34871)
#define LP (VP+float3(0.16591,-.983647,-0.06044))
#define UP float3(0.944249,0.179789,-.275808)

//#define VP vViewPosition.xyz
//uniform float4 vViewPosition;
//#define LP (VP+vViewDirection.xyz)
//uniform float4 vViewDirection;
//#define UP (vViewUp.xyz)
//uniform float4 vViewUp;

static const float2 r = float2(1280, 720);
static const float3 viewPosition = VP;
static const float3 viewDirection = normalize(LP - VP);
static const float3 viewUp = UP;
static const float fovOffset = 0.4;
static const float lensDistortion0 = -0.125;
static const float lensDistortion1 = 0;
static const float lensDistortion2 = +0.125;
static const float lensDistortion3 = 0;
static const float cam_shake = 0;

#define uTime t

float pn(float3 p) {
	p *= 1.5;
	float3 i = floor(p);
	float4 a = dot(i, float3(1, 57, 21)) + float4(0, 57, 21, 78);
	float3 f = cos((p - i)*acos(-1.))*(-.5) + .5;
	a = lerp(sin(cos(a)*a), sin(cos(1 + a)*(1 + a)), f.x);
	a.xy = lerp(a.xz, a.yw, f.y);
	return lerp(a.x, a.y, f.z);
}

float fpn(float3 p) {
	float r;
	p.xz = cos(1.) * p.xz + float2(-sin(1.), sin(1.)) * p.zx;
	r = pn(p*.06125)*.5;
	p.xy = cos(1.) * p.xy + float2(-sin(1.), sin(1.)) * p.yx;
	r += pn(p*.125)*.25;
	p.zy = cos(1.) * p.zy + float2(-sin(1.), sin(1.)) * p.yz;
	r += pn(p*.25)*.125;
	return r;
}

// Globals
float3 modXZ(float3 p, float m) {
	return p.xz = mod(p.xz + m / 2, m) - m / 2, p;
}

float max3(float3 v) {
	return max(max(v.x, v.y), v.z);
}

float fbox(float3 p, float3 b) {
	float3 d = abs(p) - b;
		return length(max(d, 0)) + max3(min(d, 0));
}

float fsphere(float3 p, float r) {
	return length(p) - r;
}

float fblob(float3 p) {
	p = abs(p);
	if (p.x <= max(p.z, p.y)) p = p.yzx;
	if (p.x <= max(p.z, p.y)) p = p.yzx;

	float b = max(max(max(
		dot(p, float3(.577, .577, .577)),
		dot(p.xz, float2(.934, .357))),
		dot(p.yx, float2(.526, .851))),
		dot(p.xz, float2(.526, .851))
		), l = length(p);

	return l - 2.5 - 0.3*(2.5 / 2)* cos(min(sqrt(1.01 - b / l)*(PI / 0.25), PI));
}

float ftest(float3 p) {
#define phi (.5*(1+sqrt(5)))
	float3 n0 = normalize(float3(-phi, phi - 1, 1));
	float3 n1 = normalize(float3(1, -phi, phi + 1));
	p = abs(p);
	p -= 2 * max(0, dot(p, n0))*n0;
	p -= 2 * max(0, dot(p, n1))*n1;
	p = abs(p);
	p -= 2 * max(0, dot(p, n0))*n0;
	p -= 2 * max(0, dot(p, n1))*n1;
	p = abs(p);
	return p.y - 1.5;
}

float fhex(float3 p, float2 h) {
	float3 q = abs(p);
	return max(q.y - h.y, max(q.z*0.866025 + q.x*0.5, q.x) - h.x);
}

float3x3 rotX(float a) {
	return float3x3 (1, 0, 0, 0, cos(a), -sin(a), 0, sin(a), cos(a));
}

// Hammer anfang
#define Ro(p, a) p=cos((a))*p+sin((a))*float2(p.y, -p.x)
float combine_d2(float d1, float d2, float r) {
	float result;
	if ((d1 < r) && (d2 < r)) {
		result = (d1 - r + d2)*sqrt(0.5);
	} else {
		result = min(d1, d2);
	}
	return result;
}
float combine_d(float d1, float d2, float r) {
	float result;
	if ((d1 < r) && (d2 < r)) {
		result = (d1 - r + d2)*sqrt(0.5);
	} else {
		result = min(d1, d2);
	}
	return result;
}

float cylinder(float3 p, float r, float c) {
	float d = length(p.xz) - r;
	d = max(d, -p.y - c);
	d = max(d, p.y - c);
	return d;
}

float hammer(float3 p) {
	//float h = box(p-float3(0,1,0),float3(2,10,3));
	//if (h > 1) return h; 

	float S = 0.5;

	p.y -= 9 * S;
	float3 q = p*(1 - 1 / (0.1 + sqrt(length(p)))*0.3 / S); //bauchiger machen
	float d = fbox(q, S*float3(1, 1, 1.4)) - 0.02*S;

	//return d;
	float2 pxy = abs(p.xy);
	Ro(pxy, -PI*0.5);
	d = max(d, pxy.x + .8);
	d = combine_d(d, fbox(q, S*float3(0.7, 0.7, 1.75)) - 0.02*S, S*0.3);
	d = combine_d2(d, cylinder(p + S*float3(0, 0.9, 0), 0.4*S, 1.6*S), S*0.1);
	d = combine_d2(d, cylinder(p + S*float3(0, 6, 0), 0.5*S, 7.3*S), S*0.1);
	d = combine_d2(d, cylinder(p + S*float3(0, 11.6, 0), 0.4*S, 1.7*S), S*0.1);
	return d;
}
//Hammer ende

static float s, s0, s1, s2, s3;
float f(float3 p) {
	// GroundPlane
	s0 = min(p.y + 3, -length(p) + 100);
	s1 = s2 = s3 = 100000;

	//s1 = fhex(p-float3(sin(acos(-1)*4./3),0.0,cos(acos(-1)*4./3)), float2(1.5, 0.3));
	//s2 = fhex(p-float3(sin(acos(-1)*2./3),0.7,cos(acos(-1)*2./3)), float2(1.5, 0.3));
	//s3 = fhex(p-float3(sin(acos(-1)*0./3),1.4,cos(acos(-1)*0/3)), float2(1.5, 0.3));

	s2 = fblob(p);//fblob(p);
	p.yx = p.xy;
	p.z += 5;
	Ro(p.yz, .5);
	s3 = hammer(p);
	return s = min(min(s0, s1), min(s2, s3));
}
float3 skydome(float3 direction) {
	direction = normalize(direction);
	float3 col = 0.4 * float3(1, 0.8, 0.5)* max(0.0, +fpn(direction * 13))*3.5;
	col += 0.3 * float3(0.8, 0.9, 1.0)* max(0.0, +fpn(direction * 8 + 100))*2.7;
	col += float3(0.5, 0.6, 1.0)* max(0.0, fpn(direction * 20 + 20) - 0.4) * 10;
	float3 p = direction;

		//col += float3(1.0,0.6,0.15)* pow(max(0.0,dot(normalize(p), normalize(float3(0.0,1,-8)))),600.0)*400; //sun
	col = pow(col, float3(2, 2, 2));
	return col*1.8;
}


float3 G(float3 p) {
	return normalize(float3(f(float3(.001, 0, 0) + p), f(float3(0, .001, 0) + p), f(float3(0, 0, .001) + p)) - f(p));
}

float2 I(float3 p, float3 d, float k, float e, float m) {
	float3 q = p;
		float i = 0, j = 0, l = sign(f(p));
	while (i < 1 && abs(s) >= j*e&&j < m)
		f(p), j += abs(s), p += l*d*s, i += k;
	return float2(2 * l*max(abs(s), 0.025)*step(i, 1)*step(j, m), distance(q, p));
}

float A(float3 p, float3 d) {
	float h = .5, i = 1, j = 1, k = 1;
	while (i < 8)
		j -= (i*h - f(p + d*i*abs(h))) / (k *= 2), i++;
	return j;
}

#define R0 0.18
float fresnel(float3 d, float3 n) {
	float x = 1. - abs(dot(n, d));
	return R0 + (1. - R0) * x*x*x*x;
}

float4 K() {
	return
		//(s==s3?float4(1.0,0.8,0.3,1): 
		(s == s3 ? float4(1.0, 0.8, 0.3, 1) :
		(s == s2 ? float4(0.8, 0.3, 0.1, 1)*1.2 :
		(s == s1 ? float4(0.3, 0.6, 0.5, 1) :
		(s == s0 ? float4(1, 1, 1, 1) :
		(s == 10 ? float4(0, 1, 1, 1) :
		float4(1, 1, 1, 1))))));
}

float4 K(float3 p, float3 d, float3 n) {
	float3 l = normalize(float3(2, 20, 17) +
		10 * float3(cos(t * 2), 0, sin(t * 2)) - p);
	float3 l2 = normalize(float3(-10, 20, 7));
	float3 l3 = normalize(float3(1, 1, 1.3));

	float4 c = K();
	if (s == s0) {
		c *= (A(p, n)*.5 + .5);
		c *= 1.3 / pow(log(length(p)), 1.3);
		c += ((0.0 + 0.5*fpn(p*float3(2, 2, 2)))*1.125) + .02*fpn(p.zxy * 1000);
		if (c.x > 0.5) c.x = 1 - c.x;
		if (c.y > 0.5) c.y = 1 - c.y;
		if (c.z > 0.5) c.z = 1 - c.z;
		//c = float4(1) - c;
	} else if (s == s1 || s == s2 || s == s3) {

		c = c*.5 + .5;
		c *= max(max(
			pow(max(abs(dot(n, normalize(l - p))), 0), 248.) * 10.0,
			pow(max(abs(dot(n, normalize(l2 - p))), 0), 248.) * 10.0),
			pow(max(abs(dot(n, normalize(l3 - p))), 0), 648.) * 200.0);
		c.xyz += skydome(reflect(d, n)) * 0.5;
	} else {
		c = c*dot(n, l);
	}

	return c;
}

float4 B(float3 p, float3 d, float3 n, float2 h) {
	float4 a = float4(1, 1, 1, 1), c = float4(0, 0, 0, 0);
	for (int i = 0; i < 8 && (f(p) == s1 || s == s2 || s == s3); ++i) {
		float3 x = d;
		a *= 1 - fresnel(d, n);
		d = refract(d, n * sign(h.x), lerp(1.05, 1 / 1.05, sign(h.x) * .5 + .5));
		if (dot(d, d) == 0)d = reflect(x, n), h.x *= -1.;
		p -= n * h.x, h = I(p, d, 1 / 64.0, 0.005, 200), p += d*h.y;
		if (h.x == 0) {

			break;
		}
		float J = 1.0;
		if (s == s3)
			J = 1.0;
		a *= lerp(exp(-h.y * J * 3 * (1 - K())), float4(1, 1, 1, 1), sign(h.x) * .5 + .5);
		n = G(p), c += K(p, d, n) * a;

	}
	return c;
}

float4 R(float3 p, float3 d, float3 n, float2 h) {
	float4 a = float4(1, 1, 1, 1), c = float4(0, 0, 0, 0);
	for (int i = 0; i < 6 && (f(p) == s1 || s == s2 || s == s3); ++i) {
		a *= fresnel(d, n);
		d = reflect(d, n);
		p += n * h.x, h = I(p, d, 1 / 64.0, 0.01, 200), p += d*h.y;
		if (h.x == 0.0) break;
		n = G(p), c += K(p, d, n) * a;
	}
	return c;
}

T void csTestScene(int3 id:SV_DispatchThreadID) {
	float3 up = viewUp;
	float3 right = -cross(viewDirection, viewUp);
	float3 p = viewPosition;
	float3 d = normalize(float3((id.xy + 0.5 - float2(0.5*r)) / r.y, -1.0 + fovOffset)); //fov is fixed
	float r = length(d.xy);
	d.xy = normalize(d.xy)*(r + r*r*lensDistortion0 + pow(r, 4)*lensDistortion1 * 2 + pow(r, 6)*lensDistortion2 * 8 + pow(r, 8)*lensDistortion3 * 16);
	d.xy *= -1;
	d = mul(transpose(float3x3(right, viewUp, -viewDirection)), d);
	d = normalize(d);

	float2 h = I(p, d, 1 / 128., 0.001, 200);
		//if(h.x==0){o0=float4(.3,0.4,.8,1);return;}
	p += d*h.y;
	float3 n = G(p);
	float4 c = K(p, d, n);
	c += B(p, d, n, h);
	c += R(p, d, n, h);
	float depth = dot(viewDirection, d*h.y);
	float coc = bkCircleOfConfusion(depth);
	out0[id.xy] = float4(c.xyz, h.y);
}