Feature/shot noise

Source/cielim/Actors/CameraModel.cpp

@@ -68,6 +68,7 @@ ACameraModel::ACameraModel()
 	this->CameraParams.QuECurveG = FVector3f::One();
 	this->CameraParams.QuECurveB = FVector3f::One();
 	this->CameraParams.CorrectionFactor = 1.0f;
+	this->CameraParams.bEnableShotNoise = false;
 	this->CameraParams.FullWellCapacity = 50000.0f;
 	this->CameraParams.Gamma = 2.2f;
 }
@@ -124,6 +125,8 @@ void ACameraModel::SetCameraParameters(const cielimMessage::CielimMessage &Cieli
 		if (SensorModel.exposuretime() > 0.0f)
 			this->CameraParams.ExposureTime = CameraModel.sensormodel().exposuretime();
 
+		this->CameraParams.bEnableShotNoise = CameraModel.sensormodel().shotnoise();
+
 		if (SensorModel.fullwellcapacity() > 0.0f)
 			this->CameraParams.FullWellCapacity = CameraModel.sensormodel().fullwellcapacity();
 

Source/cielim/Actors/CameraModel.h

@@ -36,6 +36,7 @@ struct FCameraParams
 	FVector3f QuECurveG;
 	FVector3f QuECurveB;
 	float CorrectionFactor;
+	bool bEnableShotNoise;
 	float FullWellCapacity;
 	float Gamma;
 };

Source/cielim/CielimSceneViewExtension.cpp

@@ -88,6 +88,7 @@ void FCielimSceneViewExtension::PrePostProcessPass_RenderThread(FRDGBuilder &Gra
 		CameraParams.QuECurveG = FVector3f::One();
 		CameraParams.QuECurveB = FVector3f::One();
 		CameraParams.CorrectionFactor = 1.0f;
+		CameraParams.bEnableShotNoise = false;
 		CameraParams.FullWellCapacity = 50000.0f;
 		CameraParams.Gamma = 2.2f;
 
@@ -165,6 +166,8 @@ void FCielimSceneViewExtension::QuETonemapPass(FRDGBuilder &GraphBuilder, const
 	QuEParams->QuECurveB = FVector4f(CameraParams.QuECurveB, 1.0f);
 	QuEParams->SimpsonFactor = FMath::Abs(CameraParams.Wavelength1 - CameraParams.Wavelength3) / 6.0f;
 	QuEParams->CorrectionFactor = CameraParams.CorrectionFactor;
+	QuEParams->CurrentTime = static_cast<uint32>(FDateTime::UtcNow().ToUnixTimestamp());
+	QuEParams->EnableShotNoise = static_cast<uint32>(CameraParams.bEnableShotNoise);
 	QuEParams->InvFullWellCapacity = 1.0f / FMath::Max(CameraParams.FullWellCapacity, 1e-6);
 	QuEParams->RenderTargets[0] = FRenderTargetBinding(TextureOut, ERenderTargetLoadAction::EClear);
 

Source/cielim/Shaders/QuETonemap.h

@@ -31,6 +31,8 @@ class FQuETonemap : public FGlobalShader
 	SHADER_PARAMETER(FVector4f, QuECurveB)
 	SHADER_PARAMETER(float, SimpsonFactor)
 	SHADER_PARAMETER(float, CorrectionFactor)
+	SHADER_PARAMETER(uint32, CurrentTime)
+	SHADER_PARAMETER(uint32, EnableShotNoise)
 	SHADER_PARAMETER(float, InvFullWellCapacity)
 	RENDER_TARGET_BINDING_SLOTS()
 	END_SHADER_PARAMETER_STRUCT()

Source/cielim/Shaders/QuETonemap.usf

@@ -22,8 +22,66 @@ float4 QuECurveG;
 float4 QuECurveB;
 float SimpsonFactor;
 float CorrectionFactor;
+uint CurrentTime;
+uint EnableShotNoise; // 0 or 1
 float InvFullWellCapacity; // This is in # of electrons (inverse)
 
+// RNG is a version of the process described in:
+// https://developer.nvidia.com/gpugems/gpugems3/part-vi-gpu-computing/chapter-37-efficient-random-number-generation-and-application
+
+// https://www.reedbeta.com/blog/hash-functions-for-gpu-rendering/
+uint PcgHash(uint Input)
+{
+	uint State = Input * 747796405u + 2891336453u;
+	uint Word = (State >> ((State >> 28u) + 4u) ^ State) * 277803737u;
+	return (Word >> 22u) ^ Word;
+}
+
+uint CombinedTausworthe(inout uint Z, int S1, int S2, int S3, uint M)
+{
+	uint b = (Z << S1 ^ Z) >> S2;
+	return Z = (Z & M) << S3 ^ b;
+}
+
+float UniformRNG(uint Z1, uint Z2, uint Z3, inout uint Z4)
+{
+	uint Step1 = CombinedTausworthe(Z1, 13, 19, 12, 4294967294u);
+	uint Step2 = CombinedTausworthe(Z2, 2, 25, 4, 4294967288u);
+	uint Step3 = CombinedTausworthe(Z3, 3, 11, 17, 4294967280u);
+	Z4 = Z4 * 1664525 + 1013904223u;
+	return 2.3283064365387e-10 * float(Step1 ^ Step2 ^ Step3 ^ Z4);
+}
+
+float2 BoxMuller(float U1, float U2)
+{
+	U1 = max(U1, 1e-6);
+	float R = sqrt(-2 * log(U1));
+	float Theta = 6.28318 * U2;
+	return float2(R * sin(Theta), R * cos(Theta));
+}
+
+float3 ComputeShotNoise(uint Seed, float3 NumElectrons)
+{
+	uint Z1 = PcgHash(Seed + 1);
+	uint Z2 = PcgHash(Seed + 2);
+	uint Z3 = PcgHash(Seed + 3);
+	uint Z4 = PcgHash(Seed + 4);
+
+	float U1 = UniformRNG(Z1, Z3, Z4, Z2);
+	float U2 = UniformRNG(Z2, Z3, Z1, Z4);
+	float U3 = UniformRNG(Z3, Z1, Z2, Z4);
+	float U4 = UniformRNG(Z4, Z3, Z1, Z2);
+
+	float2 RandGaussian1 = BoxMuller(U1, U2);
+	float2 RandGaussian2 = BoxMuller(U3, U4);
+
+	float ShotNoiseR = RandGaussian1.x * sqrt(NumElectrons.r);
+	float ShotNoiseG = RandGaussian1.y * sqrt(NumElectrons.g);
+	float ShotNoiseB = RandGaussian2.x * sqrt(NumElectrons.b);
+
+	return float3(ShotNoiseR, ShotNoiseG, ShotNoiseB);
+}
+
 float4 MainPS(float4 UV : TEXCOORD0, float4 Position : SV_Position) : SV_Target0
 {
 	// This is in W/m^2/nm/str
@@ -55,9 +113,19 @@ float4 MainPS(float4 UV : TEXCOORD0, float4 Position : SV_Position) : SV_Target0
 	NumElectronsG *= CorrectionFactor;
 	NumElectronsB *= CorrectionFactor;
 
-	float ColorR = saturate(NumElectronsR * InvFullWellCapacity);
-	float ColorG = saturate(NumElectronsG * InvFullWellCapacity);
-	float ColorB = saturate(NumElectronsB * InvFullWellCapacity);
+	// Compute shot noise with Gaussian approximation and add to electron count
+
+	uint Seed = (uint(Position.x) * 1973u) ^ (uint(Position.y) * 9277u) ^ (CurrentTime * 2663u);
+
+	float3 ShotNoise = EnableShotNoise * ComputeShotNoise(Seed, float3(NumElectronsR, NumElectronsG, NumElectronsB));
+
+	NumElectronsR += ShotNoise.r;
+	NumElectronsG += ShotNoise.g;
+	NumElectronsB += ShotNoise.b;
+
+	// Set final output color to NumElectrons / FWC clamped to [0-1]
+
+	float3 Color = saturate(float3(NumElectronsR, NumElectronsG, NumElectronsB) * InvFullWellCapacity);
 
-	return float4(ColorR, ColorG, ColorB, TotalRadiance.a);
+	return float4(Color, TotalRadiance.a);
 }

Source/cielim/Shaders/ReadNoise.usf

@@ -25,7 +25,6 @@ uint PcgHash(uint Input)
 	return (Word >> 22u) ^ Word;
 }
 
-
 uint CombinedTausworthe(inout uint Z, int S1, int S2, int S3, uint M)
 {
 	uint b = (Z << S1 ^ Z) >> S2;
@@ -41,7 +40,6 @@ float UniformRNG(uint Z1, uint Z2, uint Z3, inout uint Z4)
 	return 2.3283064365387e-10 * float(Step1 ^ Step2 ^ Step3 ^ Z4);
 }
 
-
 float2 BoxMuller(float U1, float U2)
 {
 	U1 = max(U1, 1e-6);
@@ -69,7 +67,7 @@ float4 MainPS(float4 UV : TEXCOORD0, float4 Position : SV_Position) : SV_Target0
 	float2 RandGaussian1 = BoxMuller(U1, U2) * ReadNoiseSigma;
 	float2 RandGaussian2 = BoxMuller(U3, U4) * ReadNoiseSigma;
 
-	float4 Noise = float4(saturate(RandGaussian1.x), saturate(RandGaussian1.y), saturate(RandGaussian2.x), 0.0f);
+	float4 ReadNoise = float4(RandGaussian1.x, RandGaussian1.y, RandGaussian2.x, 0.0f);
 
-	return Color + Noise;
+	return saturate(Color + ReadNoise);
 }