UnrealEngineUWP/Engine/Shaders/Private/ParticipatingMediaCommon.ush

/*=============================================================================
	ParticipatingMediaCommon.ush
=============================================================================*/

#pragma once

#include "Common.ush"

#define PARTICIPATING_MEDIA_MIN_MFP_METER		0.000000000001f
#define PARTICIPATING_MEDIA_MIN_EXTINCTION		0.000000000001f
#define PARTICIPATING_MEDIA_MIN_TRANSMITTANCE	0.000000000001f



//---------------------------------------------
// Participating media representation
//---------------------------------------------

struct FParticipatingMedia
{
	float3 ScatteringCoef;	// sigma_s (1/meter)
	float3 AbsorptionCoef;	// sigma_a (1/meter)
	float3 ExtinctionCoef;	// sigma_t (1/meter)
	float3 MeanFreePath;	// (meter)
	float3 Albedo;			// Represent sigma_s / sigma_t (unitless)
	float3 BaseColor;		// Represent the reflectance albedo resulting from single+multiple scattering assuming an isotropic phase function (unitless)
};

// GetBaseColorFromAlbedo: returns the color of the participating media resulting from the single+multiple subsurface scattering. 
// GetAlbedoFromBaseColor is the inverse transform.
// [Kulla and Conty 2017, "Revisiting Physically Based Shading at Imageworks"] https://blog.selfshadow.com/publications/s2017-shading-course/imageworks/s2017_pbs_imageworks_slides_v2.pdf, slide 62 
// [d'Eon, "A Hitchhiker's guide to multiple scattering"] or http://www.eugenedeon.com/wp-content/uploads/2016/09/hitchhikers_v0.1.3.pdf
float3 GetBaseColorFromAlbedo(const float3 Albedo, const float g = 0.0f)
{
	const float3 s = sqrt((1 - Albedo) / (1.0f - Albedo * g));
	const float3 BaseColor = ((1.0f - s) * (1 - 0.139 * s)) / (1.0f + 1.17 * s);
	return BaseColor;
}
float3 GetAlbedoFromBaseColor(const float3 BaseColor, const float g = 0.0f)
{
	const float3 s = 4.09712 + 4.20863 * BaseColor - sqrt(9.59217 + 41.6808 * BaseColor + 17.7126 * BaseColor * BaseColor);
	const float3 Albedo = (1.0f - s * s) / (1.0f - g * s * s);
	return Albedo;
}

// The mean free path must be in meters
FParticipatingMedia CreateMediumFromAlbedoMFP(float3 Albedo, float3 MeanFreePathMeters)
{
	FParticipatingMedia PM = (FParticipatingMedia)0;
	PM.Albedo = Albedo;
	PM.BaseColor = GetBaseColorFromAlbedo(Albedo);
	PM.MeanFreePath = MeanFreePathMeters;
	PM.ExtinctionCoef = 1.0f / max(PARTICIPATING_MEDIA_MIN_MFP_METER, PM.MeanFreePath);
	PM.ScatteringCoef = PM.Albedo * PM.ExtinctionCoef;
	PM.AbsorptionCoef = max(0.0f, PM.ExtinctionCoef - PM.ScatteringCoef);
	return PM;
}

// The mean free path must be in meters
FParticipatingMedia CreateMediumFromBaseColorMFP(float3 BaseColor, float3 MeanFreePathMeters)
{
	FParticipatingMedia PM = (FParticipatingMedia)0;
	PM.Albedo = GetAlbedoFromBaseColor(BaseColor);
	PM.BaseColor = BaseColor;
	PM.MeanFreePath = MeanFreePathMeters;
	PM.ExtinctionCoef = 1.0f / max(PARTICIPATING_MEDIA_MIN_MFP_METER, PM.MeanFreePath);
	PM.ScatteringCoef = PM.Albedo * PM.ExtinctionCoef;
	PM.AbsorptionCoef = max(0.0f, PM.ExtinctionCoef - PM.ScatteringCoef);
	return PM;
}



//---------------------------------------------
// Phase functions
//---------------------------------------------

float IsotropicPhase()
{
	return 1.0f / (4.0f * PI);
}

// Follows PBRT convention http://www.pbr-book.org/3ed-2018/Volume_Scattering/Phase_Functions.html#PhaseHG
float HenyeyGreensteinPhase(float G, float CosTheta)
{
	// Reference implementation (i.e. not schlick approximation). 
	// See http://www.pbr-book.org/3ed-2018/Volume_Scattering/Phase_Functions.html
	float Numer = 1.0f - G * G;
	float Denom = 1.0f + G * G + 2.0f * G * CosTheta;
	return Numer / (4.0f * PI * Denom * sqrt(Denom));
}

float RayleighPhase(float CosTheta)
{
	float Factor = 3.0f / (16.0f * PI);
	return Factor * (1.0f + CosTheta * CosTheta);
}

// Schlick phase function approximating henyey-greenstein
float SchlickPhaseFromK(float K, float CosTheta)
{
	const float SchlickPhaseFactor = 1.0f + K * CosTheta;
	const float PhaseValue = (1.0f - K * K) / (4.0f * PI * SchlickPhaseFactor * SchlickPhaseFactor);
	return PhaseValue;
}
float SchlickPhase(float G, float CosTheta)
{
	const float K = 1.55f * G - 0.55f * G * G * G;
	return SchlickPhaseFromK(K, CosTheta);
}

// Follows PBRT convention http://www.pbr-book.org/3ed-2018/Light_Transport_II_Volume_Rendering/Sampling_Volume_Scattering.html#SamplingPhaseFunctions
float HenyeyGreensteinPhaseInvertCDF(float E, float G)
{
	// The naive expression requires a special case for G==0, rearranging terms a bit
	// gives the following result which does not require any special cases
	float t0 = (1.0 - G) + G * E;
	float t1 = (1.0 - E) + E * E;
	float t2 = t1 + (G * E) * t0;
	float t3 = (2.0 * E - 1.0) - G * G;
	float Num = t3 + (2.0 * G) * t2;
	float Den = t0 + G * E;
	return Num / (Den * Den);
}

// Follows PBRT convention http://www.pbr-book.org/3ed-2018/Light_Transport_II_Volume_Rendering/Sampling_Volume_Scattering.html#SamplingPhaseFunctions
float4 ImportanceSampleHenyeyGreensteinPhase(float2 E, float G)
{
	float Phi = 2.0f * PI * E.x;
	float CosTheta = HenyeyGreensteinPhaseInvertCDF(E.y, G);
	float SinTheta = sqrt(max(0.0f, 1.0f - CosTheta * CosTheta));

	float3 H = float3(SinTheta * sin(Phi), SinTheta * cos(Phi), CosTheta);

	return float4(H, HenyeyGreensteinPhase(G, CosTheta));
}

// Returns CosTheta given a random value in [0,1)
float RayleighPhaseInvertCdf(float E)
{
	// [Frisvad, 2011] "Importance sampling the Rayleigh phase function"
	// Optical Society of America. Journal A: Optics, Image Science, and Vision
	float Z = E * 4.0 - 2.0;
	float InvZ = sqrt(Z * Z + 1.0);
	float u = pow(Z + InvZ, 1.0 / 3.0); // cbrt
	return u - rcp(u);
}

float4 ImportanceSampleRayleigh(float2 E)
{
	float Phi = 2.0f * PI * E.x;
	float CosTheta = RayleighPhaseInvertCdf(E.y);
	float SinTheta = sqrt(max(0.0f, 1.0f - CosTheta * CosTheta));
	float3 H = float3(SinTheta * sin(Phi), SinTheta * cos(Phi), CosTheta);
	return float4(H, RayleighPhase(CosTheta));
}


//---------------------------------------------
// Utilities
//---------------------------------------------

float3 TransmittanceToExtinction(in float3 TransmittanceColor, in float ThicknessMeters)
{
	// TransmittanceColor	= exp(-Extinction * Thickness)
	// Extinction			= -log(TransmittanceColor) / Thickness
	return -log(clamp(TransmittanceColor, PARTICIPATING_MEDIA_MIN_TRANSMITTANCE, 1.0f)) / max(PARTICIPATING_MEDIA_MIN_MFP_METER, ThicknessMeters);
}

float3 TransmittanceToMeanFreePath(in float3 TransmittanceColor, in float ThicknessMeters)
{
	return 1.0f / max(PARTICIPATING_MEDIA_MIN_EXTINCTION, TransmittanceToExtinction(TransmittanceColor, ThicknessMeters));
}

float3 ExtinctionToTransmittance(in float3 Extinction, in float ThicknessMeters)
{
	return exp(-Extinction * ThicknessMeters);
}

//---------------------------------------------
// Lighting evaluation function for an Isotropic participating media a slab of 1 meter
//---------------------------------------------

float3 IsotropicMediumSlabDirectionalAlbedoFade(float3 BaseColor, float3 MFP)
{
	float3 Fade;
	// Ensure the scattering fade out to 0 as the BaseColor approach 0. This Starts when BaseColor is below 10%.
	const float BaseColorFadesOutBelowPercentage = 10.0f;
	Fade = saturate(BaseColor * BaseColorFadesOutBelowPercentage);
	// And also fade out BaseColor to zero from MFP = 20meters (the last measured MFP for which we have fit the curve) to a large MFP=1000 meters
	const float FitLastMeasuredSampleMFP = 20.0f;
	const float AlbedoIsZeroForMFP = 1000.0f;
	Fade*= saturate(1.0f - (MFP - FitLastMeasuredSampleMFP) / (AlbedoIsZeroForMFP - FitLastMeasuredSampleMFP));

	return Fade;
}

float3 IsotropicMediumSlabPunctualDirectionalAlbedo(FParticipatingMedia PM)
{
	// Our measurement are only valid for a MFP of 0.01 meter so we clamp input MFP to that.
	// This is relatively ok since all computation are done for a slab of 1 meter.
	const float3 MFP = max(0.01f, PM.MeanFreePath);

	const float3 EvaluateForBaseColor1	= 0.0855674 / (0.237742 + (MFP + ((0.0310849 - MFP) / (1.95492 * MFP + 2.07238))));
	const float3 EvaluateForBaseColor01 = 0.0167964 / (0.541037 * (pow(1.17902, (-4.33046) / MFP) * (-0.294969 + MFP)) + 0.797592);

	// Lerp between the two extreme BaseColor curves. Measurement range was [0.1f, 1.0f].
	float3 FinalEvaluate = lerp(EvaluateForBaseColor01, EvaluateForBaseColor1, (PM.BaseColor - 0.1f) / (1.0f - 0.1f));
	return FinalEvaluate * IsotropicMediumSlabDirectionalAlbedoFade(PM.BaseColor, MFP);
}

float3 IsotropicMediumSlabEnvDirectionalAlbedo(FParticipatingMedia PM)
{
	// Our measurement are only valid for a MFP of 0.01 meter so we clamp input MFP to that.
	// This is relatively ok since all computation are done for a slab of 1 meter.
	const float3 MFP = max(0.01f, PM.MeanFreePath);

	const float3 EvaluateForBaseColor1	= 0.00231881 + (0.51379 / (pow(MFP, 1.03577) + 0.510465)); 
	const float3 EvaluateForBaseColor01 = 0.189167 / (1.55597 + (MFP + pow(0.182843, 0.0666775 + MFP)));

	// Lerp between the two extreme BaseColor curves. Measurement range was [0.1f, 1.0f].
	float3 FinalEvaluate = lerp(EvaluateForBaseColor01, EvaluateForBaseColor1, (PM.BaseColor - 0.1f) / (1.0f - 0.1f));
	return FinalEvaluate * IsotropicMediumSlabDirectionalAlbedoFade(PM.BaseColor, MFP);
}

float3 IsotropicMediumSlabTransmittance(FParticipatingMedia PM, float SlabThickness, float NoV)
{
	const float3 SafeExtinctionThreshold = 0.000001f;
	const float3 SafeExtinctionCoefficients = max(SafeExtinctionThreshold, PM.ExtinctionCoef);

	const float PathLength = SlabThickness / max(0.0001f, abs(NoV));
	const float3 SafePathSegmentTransmittance = exp(-SafeExtinctionCoefficients * PathLength);

	return SafePathSegmentTransmittance;
}