Added an experiment for hallucinated ZH3 irradiance evaluation for L1 SH

Source/Probulator/ExperimentSH.h

@@ -107,4 +107,35 @@ class ExperimentSHL1Geomerics : public Experiment
 		});
 	}
 };
+
+class ExperimentSHL1ZH3: public Experiment
+{
+public:
+	void run(SharedData& data) override
+	{
+		SphericalHarmonicsL1RGB shRadiance = {};
+
+		const ivec2 imageSize = data.m_outputSize;
+		data.m_directionImage.forPixels2D([&](const vec3& direction, ivec2 pixelPos)
+		{
+			float texelArea = latLongTexelArea(pixelPos, imageSize);
+			vec3 radiance = (vec3)data.m_radianceImage.at(pixelPos);
+			shAddWeighted(shRadiance, shEvaluateL1(direction), radiance * texelArea);
+		});
+
+		m_radianceImage = Image(data.m_outputSize);
+		m_irradianceImage = Image(data.m_outputSize);
+
+		data.m_directionImage.forPixels2D([&](const vec3& direction, ivec2 pixelPos)
+		{
+			SphericalHarmonicsL1 directionSh = shEvaluateL1(direction);
+
+			vec3 sampleSh = max(vec3(0.0f), shDot(shRadiance, directionSh));
+			m_radianceImage.at(pixelPos) = vec4(sampleSh, 1.0f);
+
+			vec3 sampleIrradianceSh = max(vec3(0.0f), shEvaluateDiffuseL1ZH3Hallucinate(shRadiance, direction) / pi);
+			m_irradianceImage.at(pixelPos) = vec4(sampleIrradianceSh, 1.0f);
+		});
+	}
+};
 }

Source/Probulator/Experiments.cpp

@@ -49,6 +49,8 @@ void addAllExperiments(ExperimentList& experiments)
 
     addExperiment<ExperimentSHL1Geomerics>(experiments, "Spherical Harmonics L1 [Geomerics]", "SHL1G");
 
+    addExperiment<ExperimentSHL1ZH3>(experiments, "Spherical Harmonics L1 [ZH3 Hallucinate]", "SHL1ZH3");
+
     addExperiment<ExperimentSH<1>>(experiments, "Spherical Harmonics L1", "SHL1");
     addExperiment<ExperimentSH<2>>(experiments, "Spherical Harmonics L2", "SHL2");
     addExperiment<ExperimentSH<3>>(experiments, "Spherical Harmonics L3", "SHL3");

Source/Probulator/SphericalHarmonics.h

@@ -248,6 +248,26 @@ namespace Probulator
 		return shEvaluateDiffuse<T, 2>(sh, direction);
 	}
 
+	inline vec3 shEvaluateDiffuseL1ZH3Hallucinate(const SphericalHarmonicsL1RGB& sh, const vec3& n)
+	{
+		// From "ZH3: Quadratic Zonal Harmonics" - https://torust.me/ZH3.pdf
+
+		const vec3 lumCoefficients = vec3(0.2126f, 0.7152f, 0.0722f);
+		const vec3 zonalAxis = normalize(vec3(glm::dot(sh[3], lumCoefficients), glm::dot(sh[1], lumCoefficients), glm::dot(sh[2], lumCoefficients)));
+
+		const vec3 ratio = vec3(glm::abs(glm::dot(vec3(sh[3].x, sh[1].x, sh[2].x), zonalAxis)),
+								glm::abs(glm::dot(vec3(sh[3].y, sh[1].y, sh[2].y), zonalAxis)),
+								glm::abs(glm::dot(vec3(sh[3].z, sh[1].z, sh[2].z), zonalAxis))) / sh[0];
+		const vec3 zonalL2Coeff = sh[0] * (0.08f * ratio + 0.6f * ratio * ratio);
+
+		const float fZ = glm::dot(zonalAxis, n);
+		const float zhDir = sqrt(5.0f / (16.0f * pi)) * (3.0f * fZ * fZ - 1.0f);
+
+		const vec3 baseIrradiance = shEvaluateDiffuseL1(sh, n);
+
+		return baseIrradiance + (pi * 0.25f * zonalL2Coeff * zhDir);
+	}
+
 	template <size_t L>
 	float shFindWindowingLambda(const SphericalHarmonicsT<float, L>& sh, float maxLaplacian)
 	{