cbuffer cb0 : register(b0) {
	float time;
	float dummy0, dummy1, dummy2;

	float3 viewPosition;
	float dummy3;

	float3 viewDirection;
	float dummy4;

	float3 viewUp;
	float dummy5;
}

//#include "ogl2dx.hlsl"
#include "utils.hlsl"
#include "lighting.hlsl"
#include "material.hlsl"
#include "sdf.hlsl"

RWTexture2D<float4> out0:register(u0);
SamplerState sampler0:register(s0);
Texture2D tex0:register(t0);

//---
// Settings
//---
#define NUM_ST_ITERATIONS 80
#define NUM_ST_ITERATIONS_REFLECTIONS 80
#define NUM_ST_ITERATIONS_REFRACTIONS 80

#define USE_SOR
#define SOR_OMEGA 1.3

#define USE_ERROR_SMOOTHING
#define NUM_ERROR_SMOOTHING_ITERATIONS 3
#define _USE_ERROR_SMOOTHING_ALONG_NORMAL

#define NUM_REFRACTIONS 8
#define NUM_REFLECTIONS 2
//#define NUM_REFRACTIONS 0
//#define NUM_REFLECTIONS 0

#define FOV 45

static const float3 lightPosition = float3(-0.5,5.5,5);
static const float lightSize = 0.75;
static const float3 lightColor = 10*float3(1,0.8,0.7);

#define GRADIENT_EPSILON_THRESHOLD 0.00045
#define GRADIENT_EPSILON_MINIMUM 0.0004
#define POSITION_EPSILON_THRESHOLD 0.0045
#define POSITION_EPSILON_MINIMUM 0.004

float FJewels(float3 p, bool evaluateMaterial, bool evaluateGradient) {

	float dGround = p.y + 1.65;
	float dOuterSphere = -fSphere(p, 10);
	
	float w = 2 * PI/3;
	float wo =  0*time*0.3;

	float3 pLasBlob = p - 1.85 * float3(sin(w*2+wo), 0, cos(w*2+wo));
	float3 pIcosahedron = p - 1.85 * float3(sin(w*1+wo), 0.1, cos(w*1+wo));
	float3 pCrystal = mul(rotX(0.6+0.05), mul(rotX(-PI/5), (p - 2 * float3(sin(w*3+wo), 0, cos(w*3+wo))).xzy).zyx);

	//pLasBlob = mul(rotX(time*0.25 + pLasBlob.x*sin(time)*0.125), pLasBlob);

	float dLasBlob = (fLasBlob(pLasBlob));
	//if (abs(dLasBlob) < 0.02)
	//	dLasBlob += fpn(p*150)*0.0045;
	//dLasBlob *= evaluateGradient ? 1 : 0.67;
	
	float dIcosahedron = fIcosahedron(pIcosahedron, 1.4, 30);
	dIcosahedron *= (evaluateGradient ? 1.025 : 1);

	float dCrystal = fCrystal(pCrystal, 1);
	float dCrystalSoft = fCrystal(pCrystal, 1, 40);
	dCrystalSoft *= (evaluateGradient ? 1.025 : 1);
	dCrystal = lerp(dCrystal, dCrystalSoft, 0.4);

	float d = INFINITY;
	minMaterial(d, dGround, evaluateMaterial, makeDiffuseMaterial(float3(0.5,0.5,0.5)));
	minMaterial(d, dOuterSphere, evaluateMaterial, makeDiffuseMaterial(float3(0.5,0.5,0.5)));
	minMaterial(d, dCrystalSoft, evaluateMaterial, makeWaterMaterial(float3(0.1, 1.9, 0.5), 1.5, 0.0));
	minMaterial(d, dLasBlob, evaluateMaterial, makeWaterMaterial(float3(0.1, 0.5, 1.9), 1.1, 0.2));
	minMaterial(d, dIcosahedron, evaluateMaterial, makeWaterMaterial(float3(1.9, 0.5, 0.1), 1.3));
	return d;
}

float fLaserEmitter(float3 p, float dScene, bool evaluateMaterial){
	float radius = length(p.xz);

	float h = radius - 1.5;
	h = max(h, p.y);

	if (h > 0.1)
		return h;

	float shaft = radius - 1;
	shaft = max(shaft, 0.8-radius);
	shaft = fOpDivideD(shaft, -p.y, 0.05);
	//shaft = max(shaft, p.y);

	float angle = 2*3.1415926/12;
	float angle2 = 2*3.1415926/48;
	float at = atan2f(p.z,p.x);
	float a = at + angle*0.5;
	float a2 = at + angle2*0.5;
	a = mod(a, angle) - angle*0.5;
	a2 = mod(a2, angle2) - angle2*0.5;
	p.xz = float2(cos(a), sin(a)) * radius;
	float3 p2 = float3(cos(a2) * radius, p.y , sin(a2) * radius);

	float3 q = p;

	p += float3(-0.9, 0.3, 0.0);
	p.z = abs(p.z);

	p.x -= 0.2;
	p.z -= 0.3;

	shaft = max(shaft, 0.2 - length(p.xz));

	q += float3(-1, 10+2, 0);
	float ring = fBox(q, float3(0.35, 10+0.2, 0.4));
	ring = fOpCombineD(ring, fBox(q,float3(0.5, 10+0.3, 0.2)), 0.05);
	//ring = min(ring, fBox(q,float3(0.5, 10+0.3, 0.2)));

	q.x -= 0.1;
	p2.x -= 0.8;
	shaft = max(shaft, 0.03 - length(p2.xz));
	shaft = min(shaft, ring);

	p += float3(-0.2,  0.6, 0.25);

	//if (p.x < p.y*2.5)
	//	fOpR45(p.xy);//p.xy = p.yx;
	p = fOpBend(p.yxz, -PI/4, 0, 0, 0.25).yxz;

	float pipe = length(p.xz) - 0.047;
	pipe = max(pipe,-0.2 - q.x);
	minMaterial(dScene, shaft, evaluateMaterial, makeDiffuseMaterial(float3(0.8, 0.9, 1)));

	//shaft = fOpCombineD(shaft, pipe, 0.01);
	//shaft = min(shaft, pipe);
	minMaterialTransmissive(shaft, pipe, evaluateMaterial, makeWaterMaterial(8*float3(0.1, 0.5, 1.9), 1.3, 0.25));

	return shaft;// * (evaluateGradient ? 1 : 0.9);
}

float FScreenSpaceMetric(float3 p, bool evaluateMaterial, bool evaluateGradient) {
	float dGround = p.y + 2.05;
	float dOuterSphere = -fSphere(p, 20);
	
	float dLasBlob =  fLasBlob(p) * (evaluateGradient ? 1 : 0.67);
	//float dSphere = fSphere(p, 1);
	float dBox = fBox(p - float3(0,0,0), float3(1,1,1)) - 0.0;
	float dSchale = max(max(fSphere(p - float3(0, 1, 0), 2), -fSphere(p - float3(0, 1.1, 0), 1.95)), p.y - 1.7) * (evaluateGradient ? 1 : 0.9);

	float dIcosahedron = fIcosahedron(p, 1.4, 30) * (evaluateGradient ? 1 : 0.9);
	float dStar = fStar(p, 1);
	float d = INFINITY;
	minMaterial(d, dGround, evaluateMaterial, makeDiffuseMaterial(3*float3(0.5,0.5,0.5)));
	minMaterial(d, dOuterSphere, evaluateMaterial, makeDiffuseMaterial(3*float3(0.5,0.5,0.5)));
	minMaterial(d, dSchale, evaluateMaterial, makeErrorTestMaterial(p*float3(0.25,0.25,200)));
	return d;
}

float FPinkBlob(float3 p, bool evaluateMaterial, bool evaluateGradient) {
	float dGround = p.y + 2.05;
	float dOuterSphere = -fSphere(p, 10);

	float dLasBlob =  fLasBlob(p);
	//if (abs(dLasBlob) < 0.02)
	//	dLasBlob += fpn(p*150)*0.0045;
	//dLasBlob *= evaluateGradient ? 1 : 1;

	float d = INFINITY;
	minMaterial(d, dGround, evaluateMaterial, makeDiffuseMaterial(float3(0.5,0.5,0.5)));
	minMaterial(d, dOuterSphere, evaluateMaterial, makeDiffuseMaterial(float3(0.5,0.5,0.5)));
	minMaterial(d, dLasBlob, evaluateMaterial, makeWaterMaterial(float3(0.1, 1.9, 0.5), 1.1));
	return d;
}

float F(float3 p, bool evaluateMaterial, bool evaluateGradient) {
	//return FSmoothErrorMaterial(p, evaluateMaterial, evaluateGradient);
	//return FSmoothErrorSphere(p, evaluateMaterial, evaluateGradient);
	//return FScreenSpaceMetric(p, evaluateMaterial, evaluateGradient);
	//return FJewels(p, evaluateMaterial, evaluateGradient);
	//return FPinkBlob(p, evaluateMaterial, evaluateGradient);
	//return FCity(p, evaluateMaterial, evaluateGradient);

	//if (!evaluateGradient)
	//	p = fOpVoxelize(p, 0.015);

	float d = INFINITY;
	float dOuterBox = -fBoxNonEuclidean(p, float3(10, 10, 10));
	minMaterial(d, dOuterBox, evaluateMaterial, makeDiffuseMaterial(float3(0.5, 0.5, 0.5)));

	//p = fOpBend(p.yzx, PI / 4 * sin(0.1 * time * PI), 0, -0.5, 0.5).zxy;
	//p = fOpBend(p.yxz, PI / 4 * cos(0.1 * time * PI), 0, -0.5, 0.5).yxz;

	//float dHG = fAkleman(p, 1, 13, 17);// fHG(p);
//	float dHG = fHexagon(p, float2(1, 1));// = fSphere(p, 2);//fPentagonalTrapezohedron(p, 1, 20);// 
	//p = fOpBend(p.yzx, PI / 4 * sin(0.1 * time * PI), 0, -0.5, 0.5).zxy;
	//float dX = fHelix(p.yxz + float3(0*0.5*time, 0, 0), 0.6, 0.6) - 0.2;
	float dHG = fLasBlob(p);

	minMaterialTransmissive(d, dHG, evaluateMaterial, makeWaterMaterial(1 * float3(0.0, 0.3333, 1.0), 1.1));

	//float omega = 2*PI/8;
	//float theta = atan2f(p.z,p.x);
	//float phi = mod(omega*0.5 + theta, omega) - omega*0.5;
	//p.xz = float2(cos(phi), sin(phi)) * length(p.xz);
	//p.x -= 3;
	//p = fOpBend(p.yxz, PI/3, -2, -2, 2).yxz;
	//
	//float dSpiralA = fHelix(p.yxz + float3(0*0.5*time, 0, 0), 0.6, 0.6);//fGlass(p);
	//dSpiralA = dSpiralA - 0.175;//max(dSpiralA - 0.175, -dSpiralA + 0.125);
	//float dCylinder = fCylinder(p, 0.4);//fSpiral(p.yxz + float3(0.5, 0, 0), 0.5, 0.5) - 0.05;
	//
	//float phi = atan2f(p.x, p.z) + PI;
	//float x = abs(frac(12*phi/(2*PI) - p.y*5) - 0.5);
	//x = smoothstep(0.24, 0.26, x) * 0.5 + 0.5;
	//
	//minMaterial(d, dCylinder, evaluateMaterial, makeEmissiveMaterial(100 * x * float3(1, 0.1, 0.01)));

	//if (abs(dSpiralA) < 0.0025)
	//	dSpiralA += fpn(p*500) * 0.000625;	

	//minMaterialTransmissive(d, dSpiralA + 0.1, evaluateMaterial, makeWaterMaterial(2*float3(8, 8., 1.), 1/1.1));
	//minMaterialTransmissive(d, dSpiralA, evaluateMaterial, makeWaterMaterial(10*float3(0.1, 0.1, 0.1), 1.1));

	return d;
}

float fMaterial(float3 p) {
	return F(p, true, false);
}

float fGradient(float3 p) {
	return F(p, false, true);
}

float f(float3 p) {
	return F(p, false, false);
}

float3 gradient(float3 p, float epsilon) {
	const float3 e[] = {
		{+epsilon, 0, 0}, {-epsilon, 0, 0},
		{0, +epsilon, 0}, {0, -epsilon, 0},
		{0, 0, +epsilon}, {0, 0, -epsilon}
	};

	float g[6];
	for (int i = 0; i < 6; ++i)
		g[i] = fGradient(p + e[i]);
	return float3(g[0]-g[1], g[2]-g[3], g[4]-g[5]);
}

float hessianSample(float3 p, float3 e, float h) {
	return (f(p + h*e) - 2 * f(p) + f(p - h*e)) / (h*h);
}

float hessianSample(float3 p, float3 ei, float3 ej, float h) {
	return (f(p + h*ei + h*ej) - f(p + h*ei - h*ej) + f(p - h*ei - h*ej) - f(p - h*ei + h*ej)) / (4*h*h);
}

float3x3 hessian(float3 p, float epsilon) {
	float3 X = float3(1, 0, 0);
	float3 Y = float3(0, 1, 0);
	float3 Z = float3(0, 0, 1);

	float3 D = float3(
		hessianSample(p, X, epsilon),
		hessianSample(p, Y, epsilon),
		hessianSample(p, Z, epsilon)
	);

	float3 O = float3(
		hessianSample(p, X, Y, epsilon),
		hessianSample(p, Y, Z, epsilon),
		hessianSample(p, Z, X, epsilon)
	);

	return float3x3(
		D.x, O.x, O.z,
		O.x, D.y, O.y,
		O.z, O.y, D.z
	);
}

float gaussianCurvature(float3 p, float3 gradient, float epsilon) {
	float3x3 H = hessian(p, epsilon);
	return dot(mul(cofactorMatrix(H), gradient), gradient);
}

float meanCurvature(float3 p, float3 gradient, float epsilon) {
	float3x3 H = hessian(p, epsilon);
	float traceH = H[0][0] + H[1][1] + H[2][2];
	return -dot(mul(H, gradient), gradient) + traceH;
	
}

float specularOcclusion(float3 p, float3 d, float dotNV, float stepSize, float numSteps) {
	float a = 1, e = 1;
	for (float i = 1; i <= numSteps; ++i)
		a -= (i * dotNV * stepSize - abs(f(p + d * i * stepSize)))/(e*=2);
	return saturate(a);
}

float ambientOcclusion(float3 p, float3 d, float stepSize, float numSteps) {
	float a = 1, e = 1;
	for (float i = 1; i <= numSteps; ++i)
		a -= (i * stepSize - abs(f(p + d * i * stepSize)))/(e*=2);
	return saturate(a);
}

float3 computeMaterialColor(float3 p, float3 d, float3 n, float gradientError) {
	float ao = ambientOcclusion(p, n, 0.3, 4);

	float3 r = reflect(d, n);
	float dotNV = saturate(-dot(n, d));
	float so =  specularOcclusion(p, r, dotNV, 0.1, 4); so *= so;

	float roughness = saturate(max(material.roughness, gradientError));
	float metallic = material.metallic;
	float f0 = lerp(0.04, 1.0, metallic);

	float3 specular = lerp(1, material.color, metallic) * lerp(1, ao, roughness) * so;
	float3 albedo = lerp(material.color, 0, metallic) * ao;

	float3 color = brdfSphereLight(p, n, -d, r, lightPosition, lightSize, albedo, specular, roughness, f0) * lightColor;
	color += material.emissive;
	return color;
}

float3 getNormalAndGradientError(float3 p, float error, out float gradientError) {
	float gradientEpsilon = max(abs(error), GRADIENT_EPSILON_MINIMUM);
	float3 g = gradient(p, gradientEpsilon);
	float gradientLength = length(g);
	gradientError = (1-saturate((gradientLength/(2*gradientEpsilon))));
	return g / gradientLength;
}

void smoothErrorAlongDirection(float3 o, float3 d, float functionSign, float pixelRadius, inout float3 p, inout float error) {
	for (int i = 0; i < NUM_ERROR_SMOOTHING_ITERATIONS; ++i) {
		p -= functionSign * d * (error - f(p));
		error = pixelRadius * functionSign * distance(o, p);
	}
}

void smoothErrorAlongNormal(float3 o, float3 n, float functionSign, float pixelRadius, inout float3 p, inout float error) {
	for (int i = 0; i < NUM_ERROR_SMOOTHING_ITERATIONS; ++i) {
		p += n * (error - f(p));
		error = pixelRadius * functionSign * distance(o, p);
	}
}

#define ST_T(X) (X).x
#define ST_SIGN(X) ((X).y < 0 ? -1 : +1)
#define ST_ABS_RADIUS(X) abs((X).y)
#define ST_RADIUS(X) (X).y
#define ST_HIT(X) ((X).x != INFINITY)

float2 sphereTrace(float3 o, float3 d, float tmax, float pixelRadius, int numIterations, bool forceHit) {
	float tmin = 0.0001;
	float omega = SOR_OMEGA, t = tmin, previousRadius = 0, stepLength = 0;
	float functionSign = f(o) < 0 ? -1 : +1;
	float2 intersection = {0, INFINITY};

	for (int i = 0; i < numIterations; ++i) {
		float signedRadius = functionSign * f(d*t + o);
		float radius = abs(signedRadius);

#ifdef USE_SOR
		bool sorFail = omega > 1 && radius + previousRadius < stepLength;
		if (sorFail) {
			stepLength = stepLength - omega * stepLength;
			omega = 1;
		} else {
			stepLength = signedRadius * omega;
		}
#else
		stepLength = signedRadius;
#endif

		previousRadius = radius;
		float screenSpaceError = radius / t;
#ifdef USE_SOR
		if (!sorFail && screenSpaceError < intersection.y)
			intersection = float2(t, screenSpaceError);
		if (!sorFail && screenSpaceError < pixelRadius || t > tmax)
			break;
#else
		if (screenSpaceError < intersection.y)
			intersection = float2(t, screenSpaceError);
		if (screenSpaceError < pixelRadius || t > tmax)
			break;
#endif
		t += stepLength;
	}

	return float2((t > tmax || intersection.y > pixelRadius) && !forceHit ? INFINITY : intersection.x,
		functionSign * max(intersection.x * intersection.y, POSITION_EPSILON_MINIMUM));
}

float3 traceReflections(float3 o, float3 d, float3 n, float2 h, float pixelRadius) {
	bool reflective = true;
	float3 pathThroughput = 1;
	float3 a = 0;

	[fastopt] for (int i = 0; i < NUM_REFLECTIONS && reflective; ++i) {
		pathThroughput *= saturate(1 - transmittance(d, n, 1, material.eta));

		d = reflect(d, n);
		o += n * ST_ABS_RADIUS(h);

		h = sphereTrace(o, d, 30, pixelRadius, NUM_ST_ITERATIONS_REFLECTIONS, true);
		if (!ST_HIT(h))
			break;

		float3 p = d * ST_T(h) + o;
		float error = pixelRadius * ST_SIGN(h) * ST_T(h);

#if defined(USE_ERROR_SMOOTHING) && !defined(USE_ERROR_SMOOTHING_ALONG_NORMAL)
		smoothErrorAlongDirection(o, d, ST_SIGN(h), pixelRadius, p, error);
#endif

		float gradientError;
		n = getNormalAndGradientError(p, error, gradientError);

#if defined(USE_ERROR_SMOOTHING) && defined(USE_ERROR_SMOOTHING_ALONG_NORMAL)
		smoothErrorAlongNormal(o, n, ST_SIGN(h), pixelRadius, p, error);
#endif

		ST_RADIUS(h) = ST_SIGN(h) * max(abs(fMaterial(p)), POSITION_EPSILON_MINIMUM);

		reflective = isReflectiveMaterial(material);
		n *= ST_SIGN(h);

		if (!reflective)
			a = pathThroughput * computeMaterialColor(p, d, n, gradientError);
		o = p;
	}

	return a;
}

float3 traceRefractions(float3 o, float3 d, float3 n, float2 h, float pixelRadius) {
	bool refractive = true;
	float3 pathThroughput = 1;
	float3 a = 0;

	[fastopt] for (int i = 0; i < NUM_REFRACTIONS && refractive; ++i) {
		float eta1 = 1.0, eta2 = material.eta;
		if (ST_SIGN(h) < 0)
			eta1 = material.eta, eta2 = 1.0;

		pathThroughput *= saturate(transmittance(d, n, eta1, eta2));
		float3 nd = refract(d, n, eta1/eta2);
		bool tir = abs(dot(nd, nd)) == 0;
		d = tir ? reflect(d, n) : nd;
		o -=  2 * n * (tir ? -1 : 1) * ST_ABS_RADIUS(h);
		//if (tir) return float3(1,0,0);

		h = sphereTrace(o, d, 30, pixelRadius, NUM_ST_ITERATIONS_REFRACTIONS, true);
		if (!ST_HIT(h))
			break;

		float3 p = d * ST_T(h) + o;
		float error = pixelRadius * ST_SIGN(h) * ST_T(h);

#if defined(USE_ERROR_SMOOTHING) && !defined(USE_ERROR_SMOOTHING_ALONG_NORMAL)
		smoothErrorAlongDirection(o, d, ST_SIGN(h), pixelRadius, p, error);
#endif

		float gradientError;
		n = getNormalAndGradientError(p, abs(error), gradientError);

		if (ST_SIGN(h) < 0)
			pathThroughput *= exp(-ST_T(h) * material.attenuation.xyz);

#if defined(USE_ERROR_SMOOTHING) && defined(USE_ERROR_SMOOTHING_ALONG_NORMAL)
		smoothErrorAlongNormal(o, n, ST_SIGN(h), pixelRadius, p, error);
#endif

		//if (isnan(dot(pathThroughput, pathThroughput)) && isinf(dot(pathThroughput, pathThroughput)))
		//	return float3(1, 0, 1);

		ST_RADIUS(h) = ST_SIGN(h) * max(abs(fMaterial(p)), POSITION_EPSILON_MINIMUM);

		refractive = isRefractiveMaterial(material);
		n = dot(d, n) > 0 ? -n : +n;

		if (!refractive || (refractive && i <= 1)) {
			a += computeMaterialColor(p, d, n, gradientError) * pathThroughput;
		}
		o = p;
	}

	return a;
}

float3 trace(float3 o, float3 d, float pixelRadius) {
	float2 h = sphereTrace(o, d, 30, pixelRadius, NUM_ST_ITERATIONS, true);

	float3 p = d * ST_T(h) + o;
	//if (ST_T(h) == INFINITY)
	//	return float3(0, 1, 1);
	float error = pixelRadius * ST_SIGN(h) * ST_T(h);

#if defined(USE_ERROR_SMOOTHING) && !defined(USE_ERROR_SMOOTHING_ALONG_NORMAL)
	smoothErrorAlongDirection(o, d, ST_SIGN(h), pixelRadius, p, error);
#endif

	float gradientError;
	float3 n = getNormalAndGradientError(p, error, gradientError) * ST_SIGN(h);

#if defined(USE_ERROR_SMOOTHING) && defined(USE_ERROR_SMOOTHING_ALONG_NORMAL)
	smoothErrorAlongNormal(o, n, ST_SIGN(h), pixelRadius, p, error);
#endif

	ST_RADIUS(h) = ST_SIGN(h) * max(abs(fMaterial(p)), POSITION_EPSILON_MINIMUM);

	float3 color = computeMaterialColor(p, d, n, gradientError);

	Material backupMaterial = material;
	if (isRefractiveMaterial(backupMaterial)) {
		color += traceRefractions(p, d, dot(d, n) > 0 ? -n : +n, h, pixelRadius);
	}

	if (isReflectiveMaterial(backupMaterial)) {
		color += traceReflections(p, d, n, h, pixelRadius);
	}

	if (ST_SIGN(h) < 0)
		color *= exp(-ST_T(h) * backupMaterial.attenuation.xyz);

	//float gradientEpsilon = thresholdedIdentity(abs(error), GRADIENT_EPSILON_THRESHOLD, GRADIENT_EPSILON_MINIMUM);
	//float gCurvature = gaussianCurvature(p, n, gradientEpsilon * 2);
	//float mCurvature = meanCurvature(p, n, gradientEpsilon * 2);
	//
	//float k1 = mCurvature + sqrt(mCurvature * mCurvature - gCurvature);
	//float k2 = 2 * mCurvature - k1;
	//float curvature = k2;

	//color = abs(curvature) * (curvature < 0 ? float3(1, 0, 0) : float3(0, 1, 0));
	//color = normalize(n) * 0.5 + 0.5;
	//color = abs(gradientError);
	return color;
}

[numthreads(16, 16, 1)] void csRT(int3 id:SV_DispatchThreadID) {
	float2 resolution;
	out0.GetDimensions(resolution.x, resolution.y);
	float2 tc = (id.xy + 0.5) / resolution;
	float2 sc = (tc*2 - 1)  * float2(resolution.x/resolution.y, -1.0);

	float3 u = cross(viewDirection, viewUp);
	float3 v = viewUp;
	float3 w = -viewDirection;

	float3 o = viewPosition;
	float3 d = normalize(mul(transpose(float3x3(u, v, w)), float3(sc/2, -0.5/tan((.5*3.1415)*FOV/180.))));

	float pixelRadius = tan((.5*3.1415)*FOV/180.)/resolution.y;
	float3 color = trace(o, d, pixelRadius);

	float4 result = float4(color, 1);

	if (isnan(dot(result, result)) || isinf(dot(result, result))) {
		out0[id.xy] = float4(0,1,0,1);
		return;
	}


	out0[id.xy] = result;
}