Improve shadow cascade quality

AnKi/Collision/Functions.h

@@ -9,6 +9,7 @@
 #include <AnKi/Collision/Plane.h>
 #include <AnKi/Collision/Ray.h>
 #include <AnKi/Collision/Aabb.h>
+#include <AnKi/Util/WeakArray.h>
 
 namespace anki {
 
@@ -225,7 +226,10 @@ void extractClipPlanes(const Mat4& mvp, Array<Plane, 6>& planes);
 void extractClipPlane(const Mat4& mvp, FrustumPlaneType id, Plane& plane);
 
 /// Compute the edges of the far plane of a frustum
-void computeEdgesOfFrustum(F32 far, F32 fovX, F32 fovY, Vec4 points[4]);
+void computeEdgesOfFrustum(F32 far, F32 fovX, F32 fovY, Vec3 points[4]);
+
+/// Welzl's algorithm that computes a compact bounding sphere given a point cloud.
+Sphere computeBoundingSphere(ConstWeakArray<Vec3> points);
 
 /// @}
 

AnKi/Collision/FunctionsMisc.cpp

@@ -4,11 +4,15 @@
 // http://www.anki3d.org/LICENSE
 
 #include <AnKi/Collision/Functions.h>
+#include <AnKi/Collision/Sphere.h>
 
 namespace anki {
 
 void extractClipPlane(const Mat4& mvp, FrustumPlaneType id, Plane& plane)
 {
+	// This code extracts the planes assuming the projection matrices are DX-like right-handed (eg D3DXMatrixPerspectiveFovRH)
+	// See "Fast Extraction of Viewing Frustum Planes from the WorldView-Projection Matrix" paper for more info
+
 #define ANKI_CASE(i, a, op, b) \
 	case i: \
 	{ \
@@ -21,7 +25,14 @@ void extractClipPlane(const Mat4& mvp, FrustumPlaneType id, Plane& plane)
 
 	switch(id)
 	{
-		ANKI_CASE(FrustumPlaneType::kNear, 3, +, 2)
+	case FrustumPlaneType::kNear:
+	{
+		const Vec4 planeEqationCoefs = mvp.getRow(2);
+		const Vec4 n = planeEqationCoefs.xyz0();
+		const F32 len = n.getLength();
+		plane = Plane(n / len, -planeEqationCoefs.w() / len);
+		break;
+	}
 		ANKI_CASE(FrustumPlaneType::kFar, 3, -, 2)
 		ANKI_CASE(FrustumPlaneType::kLeft, 3, +, 0)
 		ANKI_CASE(FrustumPlaneType::kRight, 3, -, 0)
@@ -42,17 +53,152 @@ void extractClipPlanes(const Mat4& mvp, Array<Plane, 6>& planes)
 	}
 }
 
-void computeEdgesOfFrustum(F32 far, F32 fovX, F32 fovY, Vec4 points[4])
+void computeEdgesOfFrustum(F32 far, F32 fovX, F32 fovY, Vec3 points[4])
 {
 	// This came from unprojecting. It works, don't touch it
 	const F32 x = far * tan(fovY / 2.0f) * fovX / fovY;
 	const F32 y = tan(fovY / 2.0f) * far;
 	const F32 z = -far;
 
-	points[0] = Vec4(x, y, z, 0.0f); // top right
-	points[1] = Vec4(-x, y, z, 0.0f); // top left
-	points[2] = Vec4(-x, -y, z, 0.0f); // bot left
-	points[3] = Vec4(x, -y, z, 0.0f); // bot right
+	points[0] = Vec3(x, y, z); // top right
+	points[1] = Vec3(-x, y, z); // top left
+	points[2] = Vec3(-x, -y, z); // bot left
+	points[3] = Vec3(x, -y, z); // bot right
+}
+
+Vec4 computeBoundingSphereRecursive(WeakArray<const Vec3*> pPoints, U32 begin, U32 p, U32 b)
+{
+	Vec4 sphere;
+
+	switch(b)
+	{
+	case 0:
+		sphere = Vec4(0.0f, 0.0f, 0.0f, -1.0f);
+		break;
+	case 1:
+		sphere = Vec4(*pPoints[begin - 1], kEpsilonf);
+		break;
+	case 2:
+	{
+		const Vec3 O = *pPoints[begin - 1];
+		const Vec3 A = *pPoints[begin - 2];
+
+		const Vec3 a = A - O;
+
+		const Vec3 o = 0.5f * a;
+
+		const F32 radius = o.getLength() + kEpsilonf;
+		const Vec3 center = O + o;
+
+		sphere = Vec4(center, radius);
+		break;
+	}
+	case 3:
+	{
+		const Vec3 O = *pPoints[begin - 1];
+		const Vec3 A = *pPoints[begin - 2];
+		const Vec3 B = *pPoints[begin - 3];
+
+		const Vec3 a = A - O;
+		const Vec3 b = B - O;
+
+		const Vec3 acrossb = a.cross(b);
+		const F32 denominator = 2.0f * (acrossb.dot(acrossb));
+
+		Vec3 o = b.dot(b) * acrossb.cross(a);
+		o += a.dot(a) * b.cross(acrossb);
+		o /= denominator;
+
+		const F32 radius = o.getLength() + kEpsilonf;
+		const Vec3 center = O + o;
+
+		sphere = Vec4(center, radius);
+
+		return sphere;
+	}
+#if 0
+	// There is this case as well but it fails if points are coplanar so avoid it
+	case 4:
+	{
+		const Vec3 O = *pPoints[begin - 1];
+		const Vec3 A = *pPoints[begin - 2];
+		const Vec3 B = *pPoints[begin - 3];
+		const Vec3 C = *pPoints[begin - 4];
+
+		const Vec3 a = A - O;
+		const Vec3 b = B - O;
+		const Vec3 c = C - O;
+
+		auto compDet = [](F32 m11, F32 m12, F32 m13, F32 m21, F32 m22, F32 m23, F32 m31, F32 m32, F32 m33) {
+			return m11 * (m22 * m33 - m32 * m23) - m21 * (m12 * m33 - m32 * m13) + m31 * (m12 * m23 - m22 * m13);
+		};
+
+		const F32 denominator = 2.0f * compDet(a.x(), a.y(), a.z(), b.x(), b.y(), b.z(), c.x(), c.y(), c.z());
+
+		Vec3 o = c.dot(c) * a.cross(b);
+		o += b.dot(b) * c.cross(a);
+		o += a.dot(a) * b.cross(c);
+		o /= denominator;
+
+		const F32 radius = o.getLength() + kEpsilonf;
+		const Vec3 center = O + o;
+
+		sphere = Vec4(center, radius);
+
+		return sphere;
+	}
+#endif
+	default:
+		ANKI_ASSERT(0);
+	}
+
+	for(U32 i = 0; i < p; i++)
+	{
+		const F32 distSq = (sphere.xyz() - *pPoints[begin + i]).getLengthSquared();
+		const F32 radiusSq = sphere.w() * sphere.w();
+
+		if(distSq > radiusSq)
+		{
+			for(U32 j = i; j > 0; j--)
+			{
+				const Vec3* T = pPoints[begin + j];
+				pPoints[begin + j] = pPoints[begin + j - 1];
+				pPoints[begin + j - 1] = T;
+			}
+
+			sphere = computeBoundingSphereRecursive(pPoints, begin + 1, i, b + 1);
+		}
+	}
+
+	return sphere;
+}
+
+Sphere computeBoundingSphere(ConstWeakArray<Vec3> points)
+{
+	ANKI_ASSERT(points.getSize() >= 3);
+
+	DynamicArray<const Vec3*> pPointsDyn;
+	Array<const Vec3*, 8> pPointsArr;
+	WeakArray<const Vec3*> pPoints;
+
+	if(points.getSize() > pPointsArr.getSize())
+	{
+		pPointsDyn.resize(points.getSize());
+		pPoints = pPointsDyn;
+	}
+	else
+	{
+		pPoints = {&pPointsArr[0], points.getSize()};
+	}
+
+	for(U32 i = 0; i < points.getSize(); ++i)
+	{
+		pPoints[i] = &points[i];
+	}
+
+	const Vec4 sphere = computeBoundingSphereRecursive(pPoints, 0, pPoints.getSize(), 0);
+
+	return Sphere(sphere.xyz(), sphere.w());
 }
 
 } // end namespace anki

AnKi/Math/Mat.h

@@ -1240,6 +1240,7 @@ public:
 #endif
 
 	/// Calculate a perspective projection matrix. The z is mapped in [0, 1] range just like DX and Vulkan.
+	/// Same as D3DXMatrixPerspectiveFovRH
 	ANKI_ENABLE_METHOD(kTColumnCount == 4 && kTRowCount == 4)
 	[[nodiscard]] static TMat calculatePerspectiveProjectionMatrix(T fovX, T fovY, T near, T far)
 	{
@@ -1269,6 +1270,7 @@ public:
 	}
 
 	/// Calculate an orthographic projection matrix. The z is mapped in [0, 1] range just like DX and Vulkan.
+	/// Same as D3DXMatrixOrthoOffCenterRH.
 	ANKI_ENABLE_METHOD(kTColumnCount == 4 && kTRowCount == 4)
 	[[nodiscard]] static TMat calculateOrthographicProjectionMatrix(T right, T left, T top, T bottom, T near, T far)
 	{

AnKi/Scene/Components/LightComponent.cpp

@@ -161,7 +161,7 @@ Error LightComponent::update(SceneComponentUpdateInfo& info, Bool& updated)
 		gpuLight.m_direction = -m_worldTransform.getRotation().getZAxis();
 		gpuLight.m_outerCos = cos(m_spot.m_outerAngle / 2.0f);
 
-		Array<Vec4, 4> points;
+		Array<Vec3, 4> points;
 		computeEdgesOfFrustum(m_spot.m_distance, m_spot.m_outerAngle, m_spot.m_outerAngle, &points[0]);
 		for(U32 i = 0; i < 4; ++i)
 		{
@@ -223,40 +223,33 @@ void LightComponent::computeCascadeFrustums(const Frustum& primaryFrustum, Const
 
 		// Compute a sphere per cascade
 		Array<Sphere, kMaxShadowCascades> boundingSpheres;
+		Array<Vec3, 4> prevFarPlaneEdges;
 		for(U32 cascade = 0; cascade < shadowCascadeCount; ++cascade)
 		{
-			// Compute the center of the sphere
-			//           ^ z
-			//           |
-			// ----------|---------- A(a, -f)
-			//  \        |        /
-			//   \       |       /
-			//    \    C(0,z)   /
-			//     \     |     /
-			//      \    |    /
-			//       \---|---/ B(b, -n)
-			//        \  |  /
-			//         \ | /
-			//           v
-			// --------------------------> x
-			//           |
-			// The square distance of A-C is equal to B-C. Solve the equation to find the z.
-			const F32 f = cascadeDistances[cascade]; // Cascade far
-			const F32 n = (cascade == 0) ? primaryFrustum.getNear() : cascadeDistances[cascade - 1]; // Cascade near
-			const F32 a = f * tan(fovY / 2.0f) * fovX / fovY;
-			const F32 b = n * tan(fovY / 2.0f) * fovX / fovY;
-			const F32 z = (b * b + n * n - a * a - f * f) / (2.0f * (f - n));
-			ANKI_ASSERT(absolute((Vec2(a, -f) - Vec2(0, z)).getLength() - (Vec2(b, -n) - Vec2(0, z)).getLength()) <= kEpsilonf * 100.0f);
-
-			Vec3 C(0.0f, 0.0f, z); // Sphere center
-
-			// Compute the radius of the sphere
-			const Vec3 A(a, tan(fovY / 2.0f) * f, -f);
-			const F32 r = (A - C).getLength();
-
-			// Set the sphere
-			boundingSpheres[cascade].setRadius(r);
-			boundingSpheres[cascade].setCenter(primaryFrustum.getWorldTransform().transform(C));
+			if(cascade == 0)
+			{
+				Array<Vec3, 5> edgePoints;
+				edgePoints[0] = Vec3(0.0f);
+
+				computeEdgesOfFrustum(cascadeDistances[cascade], fovX, fovY, &edgePoints[1]);
+
+				boundingSpheres[cascade] = computeBoundingSphere(edgePoints);
+
+				memcpy(&prevFarPlaneEdges[0], &edgePoints[1], sizeof(prevFarPlaneEdges));
+			}
+			else
+			{
+				Array<Vec3, 8> edgePoints;
+
+				computeEdgesOfFrustum(cascadeDistances[cascade], fovX, fovY, &edgePoints[0]);
+				memcpy(&edgePoints[4], &prevFarPlaneEdges[0], sizeof(prevFarPlaneEdges));
+
+				boundingSpheres[cascade] = computeBoundingSphere(edgePoints);
+
+				memcpy(&prevFarPlaneEdges[0], &edgePoints[0], sizeof(prevFarPlaneEdges));
+			}
+
+			boundingSpheres[cascade].setCenter(primaryFrustum.getWorldTransform().transform(boundingSpheres[cascade].getCenter()));
 		}
 
 		// Compute the matrices

AnKi/Scene/Frustum.cpp

@@ -122,10 +122,10 @@ Bool Frustum::update()
 		if(m_frustumType == FrustumType::kPerspective)
 		{
 			m_perspective.m_edgesW[0] = m_worldTransform.getOrigin();
-			m_perspective.m_edgesW[1] = m_worldTransform.transform(m_perspective.m_edgesL[0]);
-			m_perspective.m_edgesW[2] = m_worldTransform.transform(m_perspective.m_edgesL[1]);
-			m_perspective.m_edgesW[3] = m_worldTransform.transform(m_perspective.m_edgesL[2]);
-			m_perspective.m_edgesW[4] = m_worldTransform.transform(m_perspective.m_edgesL[3]);
+			m_perspective.m_edgesW[1] = m_worldTransform.transform(m_perspective.m_edgesL[0].xyz0());
+			m_perspective.m_edgesW[2] = m_worldTransform.transform(m_perspective.m_edgesL[1].xyz0());
+			m_perspective.m_edgesW[3] = m_worldTransform.transform(m_perspective.m_edgesL[2].xyz0());
+			m_perspective.m_edgesW[4] = m_worldTransform.transform(m_perspective.m_edgesL[3].xyz0());
 
 			m_perspective.m_hull = ConvexHullShape(&m_perspective.m_edgesW[0], m_perspective.m_edgesW.getSize());
 		}

AnKi/Scene/Frustum.h

@@ -239,7 +239,7 @@ private:
 		F32 m_fovX;
 		F32 m_fovY;
 		Array<Vec4, 5> m_edgesW;
-		Array<Vec4, 4> m_edgesL; ///< Don't need the eye point.
+		Array<Vec3, 4> m_edgesL; ///< Don't need the eye point.
 		ConvexHullShape m_hull;
 	};
 

AnKi/Shaders/Include/Common.h

@@ -359,6 +359,36 @@ Vec4 mul(Mat4 m, Vec4 v)
 	return Vec4(a, b, c, d);
 }
 
+Mat4 mul(Mat4 a_, Mat4 b_)
+{
+	const Vec4 a[4] = {a_.m_row0, a_.m_row1, a_.m_row2, a_.m_row3};
+	const Vec4 b[4] = {b_.m_row0, b_.m_row1, b_.m_row2, b_.m_row3};
+	Vec4 c[4];
+
+	[unroll] for(U32 i = 0; i < 4; i++)
+	{
+		Vec4 t1, t2;
+
+		t1 = a[i][0];
+		t2 = b[0] * t1;
+		t1 = a[i][1];
+		t2 += b[1] * t1;
+		t1 = a[i][2];
+		t2 += b[2] * t1;
+		t1 = a[i][3];
+		t2 += b[3] * t1;
+
+		c[i] = t2;
+	}
+
+	Mat4 o;
+	o.m_row0 = c[0];
+	o.m_row1 = c[1];
+	o.m_row2 = c[2];
+	o.m_row3 = c[3];
+	return o;
+}
+
 Vec3 mul(Mat3x4 m, Vec4 v)
 {
 	const F32 a = dot(m.m_row0, v);

AnKi/Shaders/Include/MiscRendererTypes.h

@@ -124,12 +124,4 @@ struct VolumetricLightingUniforms
 	F32 m_maxZSplitsToProcessf;
 };
 
-struct HzbUniforms
-{
-	Mat4 m_reprojectionMatrix; ///< For the main camera.
-	Mat4 m_invertedViewProjectionMatrix; ///< NDC to world for the main camera.
-	Mat4 m_projectionMatrix;
-	Mat4 m_shadowCascadeViewProjectionMatrices[kMaxShadowCascades];
-};
-
 ANKI_END_NAMESPACE

AnKi/Shaders/Intellisense.hlsl

@@ -132,6 +132,12 @@ struct Texture2D
 template<typename T>
 using RWTexture2D = Texture2D<T>;
 
+template<typename T>
+using Texture3D = Texture2D<T>;
+
+template<typename T>
+using RWTexture3D = Texture2D<T>;
+
 template<typename T>
 struct StructuredBuffer
 {
@@ -140,6 +146,12 @@ struct StructuredBuffer
 	void GetDimensions(U32& length, U32& stride);
 };
 
+template<typename T>
+struct Buffer
+{
+	T& operator[](U32 index);
+};
+
 template<typename T>
 using RWStructuredBuffer = StructuredBuffer<T>;
 
@@ -181,6 +193,12 @@ T saturate(T a);
 template<typename T>
 float dot(T a, T b);
 
+U32 asuint(float f);
+
+F32 asfloat(U32 u);
+
+U32 NonUniformResourceIndex(U32 x);
+
 // Atomics
 
 template<typename T>
@@ -189,6 +207,18 @@ void InterlockedAdd(T dest, T value, T& originalValue);
 template<typename T>
 void InterlockedAdd(T dest, T value);
 
+template<typename T>
+void InterlockedMin(T dest, T value, T& originalValue);
+
+template<typename T>
+void InterlockedMin(T dest, T value);
+
+template<typename T>
+void InterlockedMax(T dest, T value, T& originalValue);
+
+template<typename T>
+void InterlockedMax(T dest, T value);
+
 template<typename T>
 void InterlockedExchange(T dest, T value, T& originalValue);
 

AnKi/Util/Tracer.h

@@ -147,7 +147,6 @@ public:
 private:
 	const char* m_name;
 	TracerEventHandle m_handle;
-	Tracer* m_tracer;
 };
 
 #	define ANKI_TRACE_SCOPED_EVENT(name_) TracerScopedEvent _tse##name_(ANKI_STRINGIZE(ANKI_CONCATENATE(t, name_)))