GGX importance sampling for specular

projects/path_tracer/shaders/path_tracer.comp

@@ -243,12 +243,12 @@ vec3 computeSpecularBRDF(vec3 f0, float r, float NoV, float NoL, float NoH){
 
 // ---- pathtracing and main ----
 
+// distributions: https://link.springer.com/content/pdf/10.1007/978-1-4842-4427-2_16.pdf
 float cosineDistributionPDF(float NoL){
 	return NoL / pi;
 }
 
-// https://link.springer.com/content/pdf/10.1007/978-1-4842-4427-2_16.pdf
-vec3 cosineDistribution(vec2 xi){
+vec3 sampleCosineDistribution(vec2 xi){
 	float phi = 2 * pi * xi.y;
 	return vec3(
 		sqrt(xi.x) * cos(phi),
@@ -256,6 +256,37 @@ vec3 cosineDistribution(vec2 xi){
 		sqrt(xi.x) * sin(phi));
 }
 
+float uniformDistributionPDF(){
+	return 1.f / (2 * pi);
+}
+
+vec3 sampleUniformDistribution(vec2 xi){
+	float phi = 2 * pi * xi.y;
+	return vec3(
+		sqrt(xi.x) * cos(phi),
+		1 - xi.x,
+		sqrt(xi.x) * sin(phi));
+}
+
+float ggxDistributionPDF(float r, float NoH){
+	return distributionGGX(r, NoH) * NoH;
+}
+
+float ggxDistributionPDFReflected(float r, float NoH, float NoV){
+	float jacobian = 0.25 / NoV;
+	return ggxDistributionPDF(r, NoH) * jacobian;
+}
+
+vec3 sampleGGXDistribution(vec2 xi, float r){
+	float phi 		= 2 * pi * xi.y;
+	float cosTheta 	= sqrt((1 - xi.x) / ((r*r - 1) * xi.x + 1));
+	float sinTheta	= sqrt(1 - cosTheta*cosTheta);
+	return vec3(
+		cos(phi) * sinTheta,
+		cosTheta,
+		sin(phi) * sinTheta);
+}
+
 vec3 sampleTangentToWorldSpace(vec3 tangentSpaceSample, vec3 N){
 	Basis tangentBasis = buildBasisAroundNormal(N);
 	return
@@ -269,7 +300,7 @@ vec3 castRay(Ray ray) {
     vec3   	throughput 	= vec3(1);
     vec3	color   	= vec3(0);
     
-	const int maxDepth = 10;
+	const int maxDepth = 3;
     for(int i = 0; i < maxDepth; i++){
 
         Intersection intersection = sceneIntersect(ray);
@@ -278,11 +309,38 @@ vec3 castRay(Ray ray) {
         vec3 brdf 			= vec3(1);
 		
 		// V is where the ray came from and will lead back to the camera (over multiple bounces)
-		vec3 	V 			= -normalize(ray.direction);
-		vec3 	R 			= reflect(-V, intersection.N);
+		vec3 	V 				= -normalize(ray.direction);
+		vec3 	R 				= reflect(-V, intersection.N);
+		float 	NoV				= max(dot(intersection.N, V), 0);
+		
+		float 	r 				= 0.15;
+		float 	kd 				= 1 - intersection.material.ks;
+		bool 	sampleDiffuse 	= random().x < kd;
 		
-		vec3 sampleTangentSpace = cosineDistribution(random());
-        ray.direction   		= sampleTangentToWorldSpace(sampleTangentSpace, intersection.N);
+		vec3 	sampleTangentSpace;
+		float 	pdf;
+		if(sampleDiffuse){
+			sampleTangentSpace 	= sampleCosineDistribution(random());
+			ray.direction   	= sampleTangentToWorldSpace(sampleTangentSpace, intersection.N);
+			
+			float NoL 			= max(dot(intersection.N, ray.direction), 0);
+			pdf 				= cosineDistributionPDF(NoL);
+		}
+		else{			
+			#define IMPORTANCE
+			
+			#ifdef IMPORTANCE
+			sampleTangentSpace 	= sampleGGXDistribution(random(), r);
+			ray.direction   	= sampleTangentToWorldSpace(sampleTangentSpace, R);		
+			vec3 	L 			= normalize(ray.direction);
+			pdf 				= ggxDistributionPDFReflected(r, max(sampleTangentSpace.y, 0.01), max(dot(intersection.N, V), 0.01));
+			#else
+			sampleTangentSpace 	= sampleUniformDistribution(random());
+			ray.direction   	= sampleTangentToWorldSpace(sampleTangentSpace, intersection.N);		
+			pdf 				= uniformDistributionPDF();
+			#endif
+		}
+        
         ray.origin				= biasHitPosition(intersection.pos, ray.direction, intersection.N);
 		
 		// L is where the ray is going, as that is the direction where light will from
@@ -290,17 +348,14 @@ vec3 castRay(Ray ray) {
 		vec3 	H			= normalize(L + V);
 		
         float 	NoL 		= max(dot(intersection.N, L), 0);
-		float 	NoV			= max(dot(intersection.N, V), 0);
 		float	NoH			= max(dot(intersection.N, H), 0);
 		
         if(intersection.hit){
-            brdf 		= computeDiffuseBRDF(intersection.material);
-			float 	kd 	= 1 - intersection.material.ks;
-			
+            vec3 	diffuseBRDF 	= computeDiffuseBRDF(intersection.material);
+
 			vec3 	f0  			= vec3(0.04);
-			float 	r 				= 0.15;
 			vec3 	specularBRDF 	= computeSpecularBRDF(f0, r, NoV, NoL, NoH);
-			brdf 					= mix(brdf, specularBRDF, intersection.material.ks);
+			brdf 					= mix(diffuseBRDF, specularBRDF, intersection.material.ks);
 			
 			
 			hitLighting = intersection.material.emission * max(sign(NoV), 0);	// objects only emit in direction of normal
@@ -310,7 +365,7 @@ vec3 castRay(Ray ray) {
         }
         
         color       	+= hitLighting * throughput;
-        throughput  	*= brdf * NoL / max(cosineDistributionPDF(NoL), 0.0001);
+        throughput  	*= brdf * NoL / max(pdf, 0.0001);
         
         if(!intersection.hit)
             break;