Except slight issues, cook-torrance-brdf is now working and well-documented.

assets/shaders/raytracer/datastructure/mesh_bvh.glsl

@@ -25,7 +25,7 @@ bool traverseBVH(const in Mesh mesh, const in Ray local_ray, const in bool use_f
 
 	while (hit_scene)
 	{
-		if (mesh.nodes[current_node].type == 0) //Inner node.
+		while (mesh.nodes[current_node].type == 0) //Inner node.
 		{
 			int id_left = int(mesh.nodes[current_node].left_idx);
 			bool hit_left = local_ray.intersectsBounds(getNodeBounds(id_left, mesh), max_distance);
@@ -41,7 +41,6 @@ bool traverseBVH(const in Mesh mesh, const in Ray local_ray, const in bool use_f
 				bitstack = bitstack << 1;
 				current_node = id_left;
 				bitstack = bitstack | 1;
-				continue;
 			}
 			// only left hit
 			else if (hit_left && !hit_right)
@@ -49,7 +48,6 @@ bool traverseBVH(const in Mesh mesh, const in Ray local_ray, const in bool use_f
 				// Not branching here, the other sibling-check won't be needed here.
 				bitstack = bitstack << 1;
 				current_node = id_left;
-				continue;
 			}
 			// only right hit
 			else if (!hit_left && hit_right)
@@ -57,11 +55,11 @@ bool traverseBVH(const in Mesh mesh, const in Ray local_ray, const in bool use_f
 				// Not branching here, the other sibling-check won't be needed here.
 				bitstack = bitstack << 1;
 				current_node = id_right;
-				continue;
 			}
+			else break;
 		}
-		else
-		{
+
+		if (mesh.nodes[current_node].type == 1) {
 			//Is leaf
 			// intersect ray with primitives.
 			//shorten ray if closer intersection found.

assets/shaders/raytracer/properties/bsdf.glsl

@@ -10,6 +10,12 @@
 #include <screenshader/gbuffer/cook_torrance.glsl>
 #include <raytracer/properties/ggx.glsl>
 
+// BSDF IDs
+const uint eDiffuse = 0;
+const uint eReflect = 1;
+const uint eTransmit = 2;
+const uint eEmit = 3;
+
 struct BSDFResult
 {
 	// A newly generated ray for shading
@@ -64,54 +70,81 @@ SampledMaterial getMaterialParameters(const in Material material, const in Verte
 BSDFResult computeBSDF(const in Material material, const in vec2 random, const in Vertex vertex, const in Ray ray_in)
 {
 	const SampledMaterial sampled_material = material.getMaterialParameters(vertex);
+
+	// Importance sampling: compute an importance sample {phi, theta} and generate a local (z-up) microsurface normal, as well as an according normal distribution
 	HemisphereSample hemisphere_sample = ggxSample(random, sampled_material.roughness);
 
+	// inverse incoming direction
 	vec3 incoming = normalize(-ray_in.direction);
 
-	bool entering = dot(incoming, vertex.normal.xyz) > 0.f;
-	vec3 micro_normal = normalize(toWorld(hemisphere_sample.micro_normal, vertex.normal.xyz));
+	// Turn around normal if looking at it from the backside
+	vec3 normal = faceforward(vertex.normal.xyz, -incoming, vertex.normal.xyz);
+
+	// Transform the computed local microsurface normal to world space along the hit normal
+	vec3 micro_normal = normalize(toWorld(hemisphere_sample.micro_normal, normal));
+
+	// Now we can generate a direction of the ray being reflected at the microsurface normal
 	vec3 outgoing = normalize(2 * clamp(dot(incoming, micro_normal), 0, 1) * micro_normal - incoming);
 
-	float geometry = ggxGeometry(incoming, outgoing, vertex.normal.xyz, micro_normal, sampled_material.roughness);
+	// Compute geometric attenuation via GGX
+	float geometry = ggxGeometry(incoming, outgoing, normal, micro_normal, sampled_material.roughness);
+
+	// TODO: what exactly is causing the slight offset in color on metallic smooth surfaces?
+	float arbitrary_correction_factor = mix(1.f, 0.1f, sampled_material.metallic * (1-sampled_material.roughness));
 
-	float normal_response = normalResponse(sampled_material.ior, sampled_material.metallic);
-	float fresnel = fresnelSchlick(clamp(dot(incoming, micro_normal), 0, 1), normal_response);
+	// Schlick's fresnel approximation with the material's normal response
+	float normal_response = arbitrary_correction_factor * normalResponse(sampled_material.ior, sampled_material.metallic);
+	float fresnel = arbitrary_correction_factor * fresnelSchlick(clamp(dot(incoming, micro_normal), 0, 1), normal_response);
+
+	// Reflections of metallic materials are tinted with the material base color.
 	vec3 fresnel_tint = fresnelTint(fresnel, sampled_material.metallic, sampled_material.diffuse.rgb);
 
-	float n_dot_in = clamp(dot(vertex.normal.xyz, incoming), 0, 1);
-	float n_dot_out = clamp(dot(vertex.normal.xyz, outgoing), 0, 1);
+	float n_dot_in = clamp(dot(normal, incoming), 0, 1);
+	float n_dot_out = clamp(dot(normal, outgoing), 0, 1);
 	float m_dot_out = clamp(dot(micro_normal, outgoing), 0, 1);
 
+	// In this implementation, it is defined that the importance sample generated is cos(theta) instead of the usual approach of only theta.
 	float cos_theta = hemisphere_sample.importance.y;
-	float sin_theta = sqrt(1-(cos_theta*cos_theta));
-	float sin_aa = sqrt(1-(m_dot_out*m_dot_out));
 
-	float jacobian = 1 / (4 * clamp(dot(outgoing, micro_normal), 0, 1));
-	float prop_m = hemisphere_sample.distribution * cos_theta;
-	float prop_o = prop_m * jacobian;
+	// Compute the microsurface normal probability and with that the final probability density function
+	// Clamp up denominator to not divide by zero.
+	float jacobian = 1 / clamp(4 * m_dot_out, 0.001f, 1);
+	float propability_m = hemisphere_sample.distribution * cos_theta;
+	float pdf = propability_m * jacobian;
+
+	// Same here... clamp up denominator to not divide by zero
+	float brdf_denominator = clamp((4 * n_dot_in * n_dot_out), 0.001f, 1.f);
 
-	vec3 brdf = clamp(fresnel_tint * geometry * hemisphere_sample.distribution / clamp((4 * n_dot_in * n_dot_out), 0.1f, 1.f), 0, 1);
+	// Finally, build the brdf
+	vec3 brdf = fresnel_tint * geometry * hemisphere_sample.distribution / brdf_denominator;
 
 
 	BSDFResult result;
+	//Use the pixel coordinates of the incoming ray.
 	result.generated_ray = ray_in;
-	result.evaluation = (1-sampled_material.metallic)*sampled_material.diffuse.rgb * ONE_OVER_PI;
 
-	if(sampled_material.metallic == 1 || fresnel*PI > random.x)
+	// russian roulette via K_s = fresnel => K_d = 1-fresnel
+	if(sampled_material.metallic == 1 || fresnel > random.x)
 	{
+		// Use as brdf sample
 		result.generated_ray.direction = outgoing;
 		result.radiance = brdf;
-		result.probability_density = hemisphere_sample.distribution * jacobian * PI/2;
+		result.probability_density = pdf;
+		result.evaluation = vec3(0); // No diffuse if reflecting
+		result.bsdf_id = eReflect;
 	}
 	else
 	{
-		result.generated_ray.direction = normalize(toWorld(randCosineHemisphere(random.x, random.y), vertex.normal.xyz));
+		// Sample lambertian diffuse
+		result.evaluation = (1-fresnel)*(1-sampled_material.metallic)*sampled_material.diffuse.rgb * ONE_OVER_PI;
+		result.generated_ray.direction = normalize(toWorld(randCosineHemisphere(random.x, random.y), normal));
 		result.probability_density = abs(dot(vec4(result.generated_ray.direction, 0), vertex.normal)) * ONE_OVER_PI;
-		result.radiance = (1-sampled_material.metallic)*sampled_material.diffuse.rgb * ONE_OVER_PI;
+		result.radiance = (1-fresnel)*(1-sampled_material.metallic)*sampled_material.diffuse.rgb * ONE_OVER_PI;
+		result.bsdf_id = eDiffuse;
 	}
 
-	result.generated_ray.origin = vertex.position.xyz + 1e-5f * result.generated_ray.direction;
-	result.bsdf_id = 0;
+	// offset by newly generated direction
+	result.generated_ray.origin = vertex.position.xyz + 1e-8f * result.generated_ray.direction;
 	return result;
 }
 
@@ -129,6 +162,8 @@ BSDFResult computeBSDF(const in Material material, const in vec2 random, const i
 
 
 
+
+
 subroutine vec3 bsdfSampler(const in Material material, const in vec2 random, const in Vertex vertex, const in Ray ray_in, inout Ray ray_out, out vec3 light_influence, out float pdf, out uint bsdf_id);
 
 subroutine(bsdfSampler)

assets/shaders/raytracer/properties/ggx.glsl

@@ -5,13 +5,20 @@
 
 #define GGX_PI 3.14159265359f
 
+// Tidy up the importance sampling result values into a return struct
 struct HemisphereSample
 {
+  // The generated microsurface normal
   vec3 micro_normal;
+
+  // The importance sample {phi, theta}
   vec2 importance;
+
+  // The computed normal distribution function result for the generated microsurface normal
   float distribution;
 };
 
+// called F0 in several papers, it helps computing the fresnel via Schlick's approximation
 float normalResponse(float ior, float metallic)
 {
   float response = abs((1-ior) / (1+ior));
@@ -19,16 +26,21 @@ float normalResponse(float ior, float metallic)
   return mix(response, 1.f, metallic);
 }
 
-float fresnelSchlick(float view_dot_half, const in float normal_response)
+// Schlick's fresnel approximation
+// @param cos_theta_view The cosine of the incident ray and the microsurface normal. Calculated by dot(w_i, m).
+// @param normal_response A precalculated material normal response value. see float normalResponse(ior, metallic)
+float fresnelSchlick(float cos_theta_view, const in float normal_response)
 {
-    return normal_response + (1-normal_response) * pow(1 - view_dot_half, 5);
+    return normal_response + (1-normal_response) * pow(1 - cos_theta_view, 5);
 }
 
+// Thints the one-channel fresnel value depending on the metalness input value
 vec3 fresnelTint(float fresnel, float metallic, const in vec3 material_color)
 {
   return mix(vec3(1), material_color, metallic) * fresnel;
 }
 
+// Computes an importance sample {phi, theta} from any random sample {(0..1], (0..1]}.
 vec2 ggxImportanceSample(const in vec2 random_sample, const in float roughness)
 {
   float phi = random_sample.y;
@@ -37,6 +49,7 @@ vec2 ggxImportanceSample(const in vec2 random_sample, const in float roughness)
   return vec2(phi, theta);
 }
 
+// With a given random sample, this function will generate  an importance sample (output via the importance parameter), as well as an accordingly sampled z-up local microsurface normal
 vec3 ggxImportantHemisphereSample(const in vec2 random_sample, const in float roughness, out vec2 importance)
 {
   importance = ggxImportanceSample(random_sample, roughness);
@@ -46,6 +59,25 @@ vec3 ggxImportantHemisphereSample(const in vec2 random_sample, const in float ro
   return vec3(sinTheta * cos(phi), sinTheta * sin(phi), cosTheta);
 }
 
+// The funny GGX chi function
+float ggxChi(float v)
+{
+  return v > 0 ? 1 : 0;
+}
+
+// Calculates a normaal distribution for an input cos(theta) value (theta = angle between normal and microsurface normal).
+float ggxDistribution(float cosT, float roughness)
+{
+  float normal_dot_half = cosT;
+  float roughness2 = roughness * roughness;
+
+  float normal_dot_half2 = normal_dot_half * normal_dot_half;
+  float denominator = normal_dot_half2 * roughness2 + (1-normal_dot_half2);
+
+  return (ggxChi(normal_dot_half) * roughness2) / (GGX_PI * denominator * denominator);
+}
+
+// Half of the full ggx geometry attenuation term
 float ggxPartialGeometry(const in vec3 vector, const in vec3 normal, const in vec3 half_view, float roughness)
 {
   float eye_dot_half2 = clamp(dot(vector, half_view), 0, 1);
@@ -56,6 +88,7 @@ float ggxPartialGeometry(const in vec3 vector, const in vec3 normal, const in ve
   return (chi*2) / (1 + sqrt(1 + roughness * roughness * tan2));
 }
 
+// The full ggx geometric attenuation term.
 float ggxGeometry(const in vec3 vector_in, const in vec3 vector_out, const in vec3 normal, const in vec3 half_view, float roughness)
 {
   float geom_in = ggxPartialGeometry(vector_in, normal, half_view, roughness);
@@ -63,6 +96,7 @@ float ggxGeometry(const in vec3 vector_in, const in vec3 vector_out, const in ve
   return clamp(geom_in * geom_out, 0, 1);
 }
 
+// Samples an importance sample, a z-up local microsurface normal and the fitting normal distribution
 HemisphereSample ggxSample(const in vec2 random, const in float roughness)
 {
   HemisphereSample hemisphere_sample;
@@ -71,31 +105,7 @@ HemisphereSample ggxSample(const in vec2 random, const in float roughness)
   return hemisphere_sample;
 }
 
-float ggxDist(const in vec3 normal, const in vec3 vector, float roughness, float cosT)
-{
-  roughness = clamp(roughness, 1e-5f, 1.f);
-  float normal_dot_half = cosT;
-  float roughness2 = roughness * roughness;
-  float normal_dot_half2 = normal_dot_half * normal_dot_half;
-  float tan2 = (1-normal_dot_half2) / normal_dot_half2;
-  float mult = roughness2 + tan2;
-  mult *= mult;
-  float denominator = GGX_PI * normal_dot_half2 * normal_dot_half2 * mult;
-
-  float d = (ggxChi(dot(normal, vector)) * roughness2) / denominator;
-
-  if(isinf(d)) return 1;//cosT == 0 ? 1 : 0;
-  return d;
-}
-
-float ggxGeom(const in vec3 normal, const in vec3 vector, const in vec3 micro_normal, float roughness)
-{
-  float eye_dot_half2 = clamp(dot(vector, micro_normal), 0, 1);
-  float chi = ggxChi(eye_dot_half2 / clamp(dot(vector, normal), 0, 1));
-  eye_dot_half2 = eye_dot_half2 * eye_dot_half2;
-
-  float tan2 = (1 - eye_dot_half2) / eye_dot_half2;
-  return (chi*2) / (1 + sqrt(1 + roughness * roughness * tan2));
-}
+//Clean up defines
+#undef GGX_PI
 
 #endif //!INCLUDE_GGX_GLSL

assets/shaders/screenshader/gbuffer/cook_torrance.glsl

@@ -10,12 +10,12 @@ vec3 ggxHalfView(const in vec3 to_light, const in vec3 to_eye)
   return normalize(to_light + to_eye);
 }
 
-float ggxChi(float v)
+float ggxChi2(float v)
 {
   return v > 0 ? 1 : 0;
 }
 
-float ggxDistribution(float cosT, float roughness)
+float ggxDistribution2(float cosT, float roughness)
 {
   float normal_dot_half = cosT;
   float roughness2 = roughness * roughness;
@@ -23,13 +23,13 @@ float ggxDistribution(float cosT, float roughness)
   float normal_dot_half2 = normal_dot_half * normal_dot_half;
   float denominator = normal_dot_half2 * roughness2 + (1-normal_dot_half2);
 
-  return (ggxChi(normal_dot_half) * roughness2) / (GGX_PI * denominator * denominator);
+  return (ggxChi2(normal_dot_half) * roughness2) / (GGX_PI * denominator * denominator);
 }
 
 float ggxPartialGeometrys(const in vec3 to_eye, const in vec3 normal, const in vec3 half_view, float roughness)
 {
   float eye_dot_half2 = clamp(dot(to_eye, half_view), 0, 1);
-  float chi = ggxChi(eye_dot_half2 / clamp(dot(to_eye, normal), 0, 1));
+  float chi = ggxChi2(eye_dot_half2 / clamp(dot(to_eye, normal), 0, 1));
   eye_dot_half2 = eye_dot_half2 * eye_dot_half2;
 
   float tan2 = (1 - eye_dot_half2) / eye_dot_half2;