#ifndef SDF
#define SDF

#include "utils.hlsl"

static const float3 fAklemanVectors_[19] = {
	{1.000, 0.000, 0.000},
	{0.000, 1.000, 0.000},
	{0.000, 0.000, 1.000},

	{0.577, 0.577, 0.577},
	{-0.577, 0.577, 0.577},
	{0.577, -0.577, 0.577},
	{0.577, 0.577, -0.577},

	{0.000, 0.357, 0.934},
	{0.000, -0.357, 0.934},
	{0.934, 0.000, 0.357},
	{-0.934, 0.000, 0.357},
	{0.357, 0.934, 0.000},
	{-0.357, 0.934, 0.000},

	{0.000, 0.851, 0.526},
	{0.000, -0.851, 0.526},
	{0.526, 0.000, 0.851},
	{-0.526, 0.000, 0.851},
	{0.851, 0.526, 0.000},
	{-0.851, 0.526, 0.000}
};


static const float3 fAklemanVectors[19] = {
	normalize(float3(1, 0, 0)),
	normalize(float3(0, 1, 0)),
	normalize(float3(0, 0, 1)),

	normalize(float3(1, 1, 1 )),
	normalize(float3(-1, 1, 1)),
	normalize(float3(1, -1, 1)),
	normalize(float3(1, 1, -1)),

	normalize(float3(0, 1, PHI+1)),
	normalize(float3(0, -1, PHI+1)),
	normalize(float3(PHI+1, 0, 1)),
	normalize(float3(-PHI-1, 0, 1)),
	normalize(float3(1, PHI+1, 0)),
	normalize(float3(-1, PHI+1, 0)),

	normalize(float3(0, PHI, 1)),
	normalize(float3(0, -PHI, 1)),
	normalize(float3(1, 0, PHI)),
	normalize(float3(-1, 0, PHI)),
	normalize(float3(PHI, 1, 0)),
	normalize(float3(-PHI, 1, 0))
};

float3 fOpBend(float3 p) {
	if (p.y > 0) {
		p.xy = float2(length(p.xy), atan2f(p.y, p.x));
	} else if (p.x < 0) {
		p.xy *= -1;
		p.y += 2 * asin(1);
	}
	return p;
}

float3 fOpBend(float3 p, float k, float s0) {
	if (abs(k) < 0.01)
		return p;
	float sk = sign(k);
	float r = 1/k;
	float2 d = sk * float2(p.x - s0, r - p.y);
	float a = atan2f(d.x, d.y);
	return float3(a*r + s0, -length(d)*sk + r, p.z);
}

float3 fOpBend(float3 p, float angle, float s0, float sBegin, float sEnd) {
	float k = angle / (sEnd - sBegin);
	if (abs(k) < 0.0001)
		return p;
	float sk = sign(k);
	float ak = abs(k);
	float r = 1/k;
	float2 d = sk * float2(p.x - s0, r - p.y);
	float a = atan2f(d.x, d.y);
	float2 boundary = (float2(sBegin, sEnd) - s0) * k;
	float b = clamp(a, min(boundary.x, boundary.y), max(boundary.x, boundary.y));
	float s = b*r + s0;
	if (abs(a - b) > 0) {
		float2 c = float2(cos(b), sin(b));
		return float3(-sk * (c * d.x * float2(-1,1) + c.yx * d.y) + float2(s, r), p.z);
	}
	return float3(s, -length(d)*sk + r, p.z);
}

#define fOpR45(p) p = 0.7071 * p + 0.7071 * float2(p.y, -p.x)

float fOpIntersectD(float d1, float d2, float r) {
	return (-d1 < r) && (-d2 < r) ? (d1 + r + d2) * sqrt(0.5) : max(d1,d2);
}

float fOpDivideD(float d1, float d2, float r) {
	return fOpIntersectD(d1, -d2, r);
}

float fOpCombineD(float d1, float d2, float r) {
	return (d1 < r) && (d2 < r) ? (d1 - r + d2) * sqrt(0.5) : min(d1,d2);
}

float3 fOpVoxelize(float3 p, float s) {
	return floor(p/s + 0.5) * s;
}

float fTooth(float3 p) {
	// x^4 + y^4 + z^4 - x^2 + y^2 + z^2
	float d = dot(pow(p, 4), 1) - dot(pow(p, 2), 1);

	// d/dx = 4x^3 - 2x
	// d/dy = 4y^3 - 2y
	// d/dy = 4z^3 - 2z
	float3 g = 4*p*p*p - 2*p;
	return d / length(g);
}

float fTeardrop(float3 p) {
	float d = p.y*p.y*p.y*p.y + p.y*p.y*p.y + p.x*p.x + p.z*p.z;
	float3 g = float3(2*p.x, -4*p.y*p.y*p.y - 3*p.y*p.y, 2*p.z);
	return d / length(g);
}

float fKusnerSchmitt(float3 p) {
	float d = (p.x*p.x + 3) * (p.y*p.y + 3) * (p.z*p.z + 3)  - 32 * (p.x*p.y*p.z + 1);

	// d/dx = 2 x (y^2+3) (z^2+3) - 32yz
	// d/dy = 2 y (x^2+3) (z^2+3) - 32xz
	// d/dy = 2 z (x^2+3) (y^2+3) - 32xy
	float3 g = 2 * p.xyz * (p.yxx * p.yxx + 3) * (p.zzy * p.zzy + 3) - 32 * p.yxx * p.zzy;
	return d / length(g);
}

float fSteinerRoman(float3 p) {
	//x^2 + y^2 + z^2 - 1)^2 - ((z - 1)^2 - 2*x^2)*((z + 1)^2 - 2*y^2)
	float d =  (p.x*p.x + p.y*p.y + p.z*p.z - 1) * (p.x*p.x + p.y*p.y + p.z*p.z - 1) -
		((p.z - 1) * (p.z - 1) - 2*p.x*p.x) * ((p.z + 1) * (p.z + 1) - 2*p.y*p.y);

	float3 g = 4 * float3(p.x, p.y, 1) * 
		float3(
			p.x*p.x - p.y*p.y + 2*p.z*(p.z + 1),
			-p.x*p.x + p.y*p.y + 2*p.z*(p.z - 1),
			p.x*p.x*(2*p.z + 1) + p.y*p.y*(2*p.z - 1)
		);

	return d / length(g);
}

float fTetrahedralSymmetry(float3 p) {
	float d  = (p.x*p.x + p.y*p.y + p.z*p.z) * (p.x*p.x + p.y*p.y + p.z*p.z) +
		8*p.x*p.y*p.z - 10*(p.x*p.x + p.y*p.y + p.z*p.z) + 25;
	float3 g = 4 * float3(
			p.x*p.x*p.x + p.x*(p.y*p.y + p.z*p.z - 5) + 2*p.y*p.z,
			p.y*(p.x*p.x + p.z*p.z - 5) + 2*p.x*p.z + p.y*p.y*p.y,
			p.x*p.x*p.z + 2*p.x*p.y + p.z*(p.y*p.y + p.z*p.z - 5)
		);
	return d / max(16, length(g));
}

float fHyperbolicOctahedron(float3 p) {
	float3 p23 = pow(p*p, 1.0/3.0);
	float d = p23.x + p23.y + p23.z - 1;
	float3 g = 2*p / (3*p23*p23);
	return d / max(3, length(g));
}

float fMcMullenK3(float3 p) {
	float d = (p.x*p.x + 1) * (p.y*p.y + 1) * (p.z*p.z + 1) +
		8*p.x*p.y*p.z - 2;
	float3 g = 2*p*(p.yxx*p.yxx + 1)*(p.zzy*p.zzy + 1) +
		8*p.yxx*p.zzy;
	return d / max(1, length(g));
}

float fAkleman(float3 p, float r, float e, int begin, int end) {
	float d = 0;
	[unroll] for (int i = begin; i <= end; ++i)
		d += pow(abs(dot(p, fAklemanVectors[i])), e);
	return pow(d, 1/e) - r;
}

float fAkleman(float3 p, float r, int begin, int end) {
	float d = 0;
	[unroll] for (int i = begin; i <= end; ++i)
		d = max(d, abs(dot(p, fAklemanVectors[i])));
	return d - r;
}

float fPowAkleman(float3 p, float e, int begin, int end) {
	float d = 0;
	[unroll] for (int i = begin; i <= end; ++i)
		d += pow(abs(dot(p, fAklemanVectors[i])), e);
	return d;
}

float fPyramid(float3 p, float r, float e) {
	float d = 0;
	[unroll] for (int i = 3; i <= 6; ++i)
		d += pow(abs(dot(p, fAklemanVectors[i])), e);
	d += pow(max(0, -dot(p - float3(0,0.95,0), float3(0,1,0))), e);
	return pow(d, 1/e) - r;
}

float fCrystal(float3 p, float r, float e) {
	float d = 0;
	[unroll] for (int i = 0; i < 6; ++i) {
		float w0 = (PI/3 * i), w1 = w0 + PI/6;
		d += pow(max(0, -dot(p, normalize(float3(cos(w0), 0.75, sin(w0))))), e);
		d += pow(max(0, +dot(p, normalize(float3(cos(w1), 1.5, sin(w1))))), e);
	}
	d += pow(max(0, dot(p + float3(0, 0.5, 0), float3(0,1,0))), e);
	return pow(d, 1/e) - r;
}

float fCrystal(float3 p, float r) {
	float d = 0;
	[unroll] for (int i = 0; i < 6; ++i) {
		float w0 = (PI/3 * i), w1 = w0 + PI/6;
		d = max(d, -dot(p, normalize(float3(cos(w0), 0.75, sin(w0)))));
		d = max(d, dot(p, normalize(float3(cos(w1), 1.5, sin(w1)))));
	}
	d = max(d, dot(p + float3(0, 0.5, 0), float3(0,1,0)));
	return d - r;
}

float fPentagonalTrapezohedron(float3 p, float r, float e) {
	float w = 2 * PI / 5;
	float3 v[] = {
		normalize((float3(sin(w * 0), 0.85, cos(w * 0)))),
		normalize((float3(sin(w * 1), 0.85, cos(w * 1)))),
		normalize((float3(sin(w * 2), 0.85, cos(w * 2)))),
		normalize((float3(sin(w * 3), 0.85, cos(w * 3)))),
		normalize((float3(sin(w * 4), 0.85, cos(w * 4)))),
	};

	float d = 0;
	[unroll] for (int i = 0; i < 5; ++i)
		d += pow(abs(dot(p, v[i])), e);
	return pow(d, 1/e) - r;
}

float fPentagonalTrapezohedron(float3 p, float r) {
	float w = 2 * PI / 5;
	float3 v[] = {
		normalize((float3(sin(w * 0), 0.85, cos(w * 0)))),
		normalize((float3(sin(w * 1), 0.85, cos(w * 1)))),
		normalize((float3(sin(w * 2), 0.85, cos(w * 2)))),
		normalize((float3(sin(w * 3), 0.85, cos(w * 3)))),
		normalize((float3(sin(w * 4), 0.85, cos(w * 4)))),
	};

	float d = 0;
	[unroll] for (int i = 0; i < 5; ++i)
		d = max(d, abs(dot(p, v[i])));
	return d - r;
}


float fDodecahedron(float3 p, float r, float e) {
	return fAkleman(p, r, e, 13, 18);
}

float fOctahedron(float3 p, float r, float e) {
	return fAkleman(p, r, e, 3, 6);
}

float fIcosahedron(float3 p, float r, float e) {
	return fAkleman(p, r, e, 3, 12);
}

float fTruncatedOctahedron(float3 p, float r, float e) {
	return fAkleman(p, r, e, 0, 6);
}

float fTruncatedIcosahedron(float3 p, float r, float e) {
	return fAkleman(p, r, e, 3, 18);
}

float fOctahedron(float3 p, float r) {
	return fAkleman(p, r, 3, 6);
}

float fDodecahedron(float3 p, float r) {
	return fAkleman(p, r, 13, 18);
}

float fIcosahedron(float3 p, float r) {
	return fAkleman(p, r, 3, 12);
}

float fTruncatedOctahedron(float3 p, float r) {
	return fAkleman(p, r, 0, 6);
}

float fTruncatedIcosahedron(float3 p, float r) {
	return fAkleman(p, r, 3, 18);
}

float fBox(float3 p, float3 b){
	float3 d = abs(p) - b;
	return length(max(d, 0)) + min(max(max(d.x, d.y), d.z), 0);
}

float fBoxNonEuclidean(float3 p, float3 b) {
	float3 d = abs(p) - b;
	return max(max(d.x, d.y), d.z);
}

float fCylinder(float3 p, float r) {
	return length(p.xz) - r;
}

float fLineSegment(float3 p, float3 a, float3 b) {
	float3 ab = b - a;
	float t = saturate(dot(p - a, ab) / dot(ab, ab));
	return length((ab*t + a) - p);
}

float fTorus(float3 p, float rInner, float rOuter) {
	return length(float2(length(p.xz) - rInner, p.y)) - rOuter;
}

float fCircle(float3 p, float r) {
	float l = length(p.xz) - r;
	return length(float2(p.y, l));
}

float fDisc(float3 p, float r) {
	float l = length(p.xz) - r;
	return l < 0 ? abs(p.y) : length(float2(p.y, l));
}

float fCone(float3 p, float2 d) {
	return dot(d, float2(length(p.xz), p.y));
}

float fHexagon(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);
}

float fHelix(float3 p, float f, float r) {
	float s = 0.5 * f / PI;
	float x = PI * s;
	float a = atan2f(p.z, p.y) * s;
	float b = mod(p.x, 2*x);
	a = abs(a-b);
	a = a > x ? 2*x - a : a;
	return length(float2(length(p.yz) - r, a));
}

float fDiscNonEuclidean(float3 p, float r) {
	return max(length(p.xz) - r, abs(p.y));
}

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

float fStar(float3 p, float r) {
	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;
	
	return p.y - r;
}

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

	float b = max(max(max(
		dot(p, normalize(float3(1, 1, 1))),
		dot(p.xz, normalize(float2(PHI+1, 1)))),
		dot(p.yx, normalize(float2(1, PHI)))),
		dot(p.xz, normalize(float2(1, PHI))));

	float l = length(p);
	
	return l - 1.6125 + 0.25 * (smoothstepQuintic(0.0, 0.343, sin(acos(b/l-0.0075))));
	//return l - 1.5 + 1.5 * saturate(1-b/l); // Bluming
	return l - 1.5 - 0.2 * (1.5 / 2)* cos(min(sqrt(1.01 - b / l)*(PI / 0.25), PI));
}

float fHG(float3 p, float depthFactor = 0.5) {
	float hgUnit = 1.0 / 8.0;
	float depth = hgUnit * depthFactor;
	float d = INFINITY;
	d = min(d, fBox(p - float3(0, 0, 0), float3(hgUnit * 4, hgUnit, depth)));
	d = min(d, fBox(p - float3(0, -hgUnit * 2.5, 0), float3(hgUnit, hgUnit * 3.5, depth)));
	p.x = abs(p.x);
	d = min(d, fBox(p - float3(hgUnit * 7.5, 0, 0), float3(hgUnit * 0.5, hgUnit * 8, depth)));
	d = min(d, fBox(p - float3(hgUnit * 3, hgUnit * 2.5, 0), float3(hgUnit, hgUnit * 3.5, depth)));
	p.y = abs(p.y);
	d = min(d, fBox(p - float3(0, hgUnit * 7.5, 0), float3(hgUnit * 8, hgUnit * 0.5, depth)));
	d = min(d, fBox(p - float3(0, hgUnit * 3, 0), float3(hgUnit * 4, hgUnit, depth)));
	return d;
}

float fGlass(float3 p) {
	float l = length(p.xz);
	float3 q = float3(l, p.y, 0);
	q = fOpBend(q.yxz, PI/3, -0.25, -0.25, 0.0).yxz;
	q = fOpBend(q.yxz, -PI/2, 0.0, 0.0, 1.5).yxz;
	q = fOpBend(q.yxz, -PI/2, -0.75, -0.75, -1.5).yxz;
	return fBox(q, float3(0.05-abs(q.y)*0.02, 1.5, 0.1)) * 0.8;
}

float fGummibaer(float3 p) {
	if (length(p) > 1.41)
		return length(p) - 1.4;

	const float4 data[] = {
		// Koerper
		float4(1.0, 0.0, 0.0, 2.5),
		float4(-1.0, 0.0, 0.0, 2.5),
		float4(1.0, 2.0, 0.0, 2.5),
		float4(-1.0, 2.0, 0.0, 2.5),
		float4(1.0, 4.0, 0.0, 1.0),
		float4(-1.0, 4.0, 0.0, 1.0),
		
		// Arme / Beine
		float4(1.5, -0.0, 1.5, 1.2),
		float4(-1.5, -0.0, 1.5, 1.2),
		float4(1.55, 2.0, 1.5, 0.0),
		float4(-1.5, 2.0, 1.5, 0.0),
		
		// Ohren
		float4(1.3, 4.8, -0.2, 0.3),
		float4(-1.3, 4.8, -0.2, 0.3),

		// Gesicht
		float4(-0.0, 4.0, 1.0, 0.0)
	};

	float globalScale = 0.1;
	p /= globalScale;

	p.y += 4.8 / 2;

	float b = 0;
	for (int i = 0; i < 13; ++i)
		b += 0.5 - atan(length(p - 3 * data[i].xyz) -  data[i].w) / PI;

	return ((tan((-clamp(b, 0.0001, 0.78) - 0.5) * PI) + 1)) * globalScale;
}

#endif