#include "World/DescentSDF.h"
#include "World/DescentWorldDef.h"
#include "Core/Tuning/DescentTune_World.h"

namespace WTune = Descent::Tune::World;

namespace Descent::World
{
	FDescentSDF::FDescentSDF(const FDescentWorldDef& InDef)
		: Def(InDef)
		, NoiseA(WTune::NoiseSeedA)
		, NoiseBObj(WTune::NoiseSeedB)
	{
	}

	// ---------------------------------------------------------------------------------------------------------
	// Primitive distance functions (port of 03_world.js 129 - 165 with the axes of section 0)
	// ---------------------------------------------------------------------------------------------------------

	double FDescentSDF::SegDistR(const FSegment& S, double X, double Y, double Z)
	{
		const double Bx = S.B.X - S.A.X, By = S.B.Y - S.A.Y, Bz = S.B.Z - S.A.Z;
		const double Px = X - S.A.X, Py = Y - S.A.Y, Pz = Z - S.A.Z;
		double T = (Px * Bx + Py * By + Pz * Bz) / (Bx * Bx + By * By + Bz * Bz);
		T = T < 0.0 ? 0.0 : (T > 1.0 ? 1.0 : T);
		const double Dx = Px - Bx * T, Dy = Py - By * T, Dz = Pz - Bz * T;
		return FMath::Sqrt(Dx * Dx + Dy * Dy + Dz * Dz) - (S.Ra + (S.Rb - S.Ra) * T);
	}

	double FDescentSDF::SegDist3(const FVector& A, const FVector& B, double X, double Y, double Z)
	{
		const double Bx = B.X - A.X, By = B.Y - A.Y, Bz = B.Z - A.Z;
		const double Px = X - A.X, Py = Y - A.Y, Pz = Z - A.Z;
		double T = (Px * Bx + Py * By + Pz * Bz) / (Bx * Bx + By * By + Bz * Bz);
		T = T < 0.0 ? 0.0 : (T > 1.0 ? 1.0 : T);
		const double Dx = Px - Bx * T, Dy = Py - By * T, Dz = Pz - Bz * T;
		return FMath::Sqrt(Dx * Dx + Dy * Dy + Dz * Dz);
	}

	double FDescentSDF::WellSDF(const FWell& W, double X, double Y, double Z)
	{
		return FMath::Max3(FMath::Sqrt((X - W.X) * (X - W.X) + (Y - W.Y) * (Y - W.Y)) - W.R, Z - W.Z1, W.Z0 - Z);
	}

	double FDescentSDF::SlotSDF(const FShelf& S, double X, double Y)
	{
		const double Dx = X - S.X, Dy = Y - S.Y;
		const double A = S.Axis.X * Dx + S.Axis.Y * Dy;       // along the slot axis
		const double O = -S.Axis.Y * Dx + S.Axis.X * Dy;      // across
		const double Lim = S.Half - S.W;
		const double Al = FMath::Max(-Lim, FMath::Min(Lim, A));
		return FMath::Sqrt((A - Al) * (A - Al) + O * O) - S.W;
	}

	double FDescentSDF::ShelfRock(const FShelf& S, double X, double Y, double Z, double Noise)
	{
		const double Slab = S.T + Noise - FMath::Abs(Z - S.Zc);
		const double Clip = S.Rclip - FMath::Sqrt((X - S.X) * (X - S.X) + (Y - S.Y) * (Y - S.Y));
		double V = FMath::Min3(Slab, SlotSDF(S, X, Y), Clip);
		for (int32 K = 0; K + 1 < S.Crack.Num(); ++K)
		{
			V = FMath::Min(V, SegDist3(S.Crack[K], S.Crack[K + 1], X, Y, Z) - S.CrackR);
		}
		return V;
	}

	double FDescentSDF::PinchRock(const FPinch& P, double X, double Y, double Z)
	{
		const double Dx = X - P.C.X, Dy = Y - P.C.Y, Dz = Z - P.C.Z;
		const double Along = Dx * P.Dir.X + Dy * P.Dir.Y + Dz * P.Dir.Z;
		const double U = Dx * P.Side.X + Dy * P.Side.Y + Dz * P.Side.Z;
		const double V = Dx * P.Up.X + Dy * P.Up.Y + Dz * P.Up.Z;
		const double Ea = U / P.Ax, Eb = V / P.Ay;
		const double El = (FMath::Sqrt(Ea * Ea + Eb * Eb) - 1.0) * FMath::Min(P.Ax, P.Ay);
		return FMath::Min3(P.Half - FMath::Abs(Along), El, (P.TubeR + 14.0) - FMath::Sqrt(U * U + V * V));
	}

	double FDescentSDF::RingRock(const FRing& R, const FWell& W, double X, double Y, double Z)
	{
		const double H = FMath::Sqrt((X - W.X) * (X - W.X) + (Y - W.Y) * (Y - W.Y));
		return FMath::Min3(R.Half - FMath::Abs(Z - R.Zc), H - R.RInner, (W.R + 12.0) - H);
	}

	double FDescentSDF::AbyssSDF(double X, double Y, double Z)
	{
		const double Rz = Descent::Lerp(WTune::AbyssRadiusTop, WTune::AbyssRadiusBottom, Descent::Smooth(WTune::AbyssTopZ, WTune::AbyssBottomZ, Z));
		const double H = FMath::Sqrt((X - WTune::AbyssX) * (X - WTune::AbyssX) + (Y - WTune::AbyssY) * (Y - WTune::AbyssY));
		// Z + 1400 is positive above the shaft's top, -1640 - Z positive under its bottom: both close it.
		return FMath::Max3(H - Rz - 5.0 * FMath::Sin(0.05 * Z + 0.04 * Y), Z - WTune::AbyssTopZ, WTune::AbyssBottomZ - Z);
	}

	double FDescentSDF::BoxSDF(const FBox& B, double X, double Y, double Z)
	{
		const FVector C = (B.Min + B.Max) * 0.5, H = (B.Max - B.Min) * 0.5;
		const double Qx = FMath::Abs(X - C.X) - H.X, Qy = FMath::Abs(Y - C.Y) - H.Y, Qz = FMath::Abs(Z - C.Z) - H.Z;
		const double Ox = FMath::Max(Qx, 0.0), Oy = FMath::Max(Qy, 0.0), Oz = FMath::Max(Qz, 0.0);
		return FMath::Sqrt(Ox * Ox + Oy * Oy + Oz * Oz) + FMath::Min(FMath::Max3(Qx, Qy, Qz), 0.0);
	}

	double FDescentSDF::CanyonCentre(double X)
	{
		return 30.0 * FMath::Sin(0.006 * X) + 12.0 * FMath::Sin(0.017 * X);
	}

	double FDescentSDF::CanyonFloor(double X)
	{
		return Descent::Lerp(-120.0, -260.0, Descent::Smooth(-100.0, 800.0, X));
	}

	double FDescentSDF::SpongeSDF(const FSpongeZone& Zone, double X, double Y, double Z) const
	{
		// Two ridged noises intersected give worm-like tunnels 5 - 9 m wide (section 2.2).
		const double NA = NoiseA(0.016 * X + 37.0, 0.016 * Y + 37.0, 0.016 * Z + 37.0);
		const double NB = NoiseBObj(0.016 * X - 53.0, 0.016 * Y - 53.0, 0.016 * Z - 53.0);
		double V = (FMath::Max(FMath::Abs(NA), FMath::Abs(NB)) - WTune::SpongeNoiseCut) * WTune::SpongeScale;
		V = FMath::Max(V, BoxSDF(Zone.Box, X, Y, Z) + WTune::SpongeBoxGrow);
		if (V < 0.0)
		{
			// A procedural void never opens a bypass around a choke: masked within 30 m of a well or gallery (section 2.9).
			const double D = Def.CarrierRouteDistance(FVector(X, Y, Z));
			V = FMath::Max(V, WTune::NearCarrierRouteM - D);
		}
		return V;
	}

	// ---------------------------------------------------------------------------------------------------------
	// The density function (section 2.2)
	// ---------------------------------------------------------------------------------------------------------

	double FDescentSDF::GenRaw(const FVector& P, const FPrimList& L) const
	{
		const double X = P.X, Y = P.Y, Z = P.Z;
		double Open = 1e9, Obig = 1e9;

		for (const FSegment* S : L.Segs)
		{
			const double D = SegDistR(*S, X, Y, Z);
			if (S->bBig)
			{
				if (D < Obig) { Obig = D; }
			}
			else if (D < Open)
			{
				Open = D;
			}
		}

		if (L.bCanyon)
		{
			const double Fl = CanyonFloor(X) + 4.0 * FMath::Sin(0.07 * Y - 0.02 * X) + 3.0 * FMath::Sin(0.045 * X);
			const double Hw = WTune::CanyonHalfWidth0 + WTune::CanyonWallSlope * FMath::Max(0.0, Z - Fl);
			double Cyn = FMath::Max(FMath::Max(FMath::Abs(Y - CanyonCentre(X)) - Hw, Fl - Z), FMath::Max(X - WTune::CanyonX1, WTune::CanyonX0 - X));
			// Everything above the rim is open sea (and below it the world is rock out to infinity: the sea floor).
			Cyn = FMath::Min(Cyn, (WTune::RimZ + WTune::RimWobble * FMath::Sin(0.03 * Y) * FMath::Cos(0.025 * X)) - Z);
			Open = Descent::SMin(Open, Cyn, WTune::CanyonBlendK);
		}

		for (const FChamber* C : L.Chambers)
		{
			const double Dx = (X - C->C.X) / C->R.X, Dy = (Y - C->C.Y) / C->R.Y, Dz = (Z - C->C.Z) / C->R.Z;
			const double MinR = FMath::Min3(C->R.X, C->R.Y, C->R.Z);
			Open = Descent::SMin(Open, (FMath::Sqrt(Dx * Dx + Dy * Dy + Dz * Dz) - 1.0) * MinR, C->K);
		}
		for (const FWell* W : L.Wells)
		{
			Obig = Descent::SMin(Obig, WellSDF(*W, X, Y, Z), WTune::BigWellBlendK);
		}
		if (L.bAbyss)
		{
			Open = Descent::SMin(Open, AbyssSDF(X, Y, Z), WTune::AbyssBlendK);
		}
		for (const FPore* Pr : L.Pores)
		{
			const double Dx = X - Pr->C.X, Dy = Y - Pr->C.Y, Dz = Z - Pr->C.Z;
			Open = Descent::SMin(Open, FMath::Sqrt(Dx * Dx + Dy * Dy + Dz * Dz) - Pr->R, WTune::PoreBlendK);
		}
		for (const FSpongeZone* Zone : L.Sponges)
		{
			Open = Descent::SMin(Open, SpongeSDF(*Zone, X, Y, Z), WTune::SpongeBlendK);
		}

		const double N1 = NoiseA(0.02 * X, 0.02 * Y, 0.024 * Z);
		const double N2 = NoiseA(0.07 * X + 11.0, 0.07 * Y, 0.07 * Z);
		double Rock = Open + WTune::AmpN1 * N1 + WTune::AmpN2 * N2;
		if (Z < WTune::DeepZ)
		{
			Rock += WTune::AmpDeepB * NoiseBObj(0.05 * X, 0.05 * Y, 0.05 * Z);
		}
		if (Obig < 1e8)
		{
			// Carrier voids: low noise, so their clearances are trustworthy (section 2.8).
			const double Rb = Obig + WTune::AmpBigN1 * N1 + WTune::AmpBigN2 * N2;
			Rock = (Rock < 1e8) ? Descent::SMin(Rock, Rb, WTune::BigBlendK) : Rb;
		}

		// The chokes: rock added back with max.
		for (const FShelf* S : L.Shelves)
		{
			const double Noise = 1.4 * NoiseBObj(0.09 * X, 0.09 * Y, 0.09 * Z) + 0.8 * N2;
			const double V = ShelfRock(*S, X, Y, Z, Noise);
			if (V > Rock) { Rock = V; }
		}
		for (const FPinch* Pn : L.Pinches)
		{
			const double V = PinchRock(*Pn, X, Y, Z) + 0.7 * N2;
			if (V > Rock) { Rock = V; }
		}
		for (const FRing* R : L.Rings)
		{
			const FWell& W = Def.Wells[R->Well];
			const double V = RingRock(*R, W, X, Y, Z) + 0.7 * N2;
			if (V > Rock) { Rock = V; }
		}
		return Rock;
	}

	double FDescentSDF::Gen(const FVector& P, const FPrimList& L) const
	{
		const double V = GenRaw(P, L);
		return V > WTune::ClampD ? WTune::ClampD : (V < -WTune::ClampD ? -WTune::ClampD : V);
	}

	void FDescentSDF::Gather(const FVector& Centre, double Radius, FPrimList& Out) const
	{
		Def.GatherInto(Centre, Radius, Out);
	}

	double FDescentSDF::GenRawSlow(const FVector& Pm) const
	{
		FPrimList L;
		Gather(Pm, 24.0, L);
		return GenRaw(Pm, L);
	}

	FVector FDescentSDF::GradientRaw(const FVector& P, const FPrimList& L, double H) const
	{
		const FVector G(GenRaw(P + FVector(H, 0, 0), L) - GenRaw(P - FVector(H, 0, 0), L),
			GenRaw(P + FVector(0, H, 0), L) - GenRaw(P - FVector(0, H, 0), L),
			GenRaw(P + FVector(0, 0, H), L) - GenRaw(P - FVector(0, 0, H), L));
		const double Len = G.Size();
		return Len > 1e-9 ? G / Len : FVector(0, 0, 1);
	}
}