added a HK GLSL shader

data/tests/test_bsdf.xml

@@ -53,6 +53,18 @@
 		<bsdf type="diffuse"/>
 	</bsdf>
 
+	<!-- Test the Hanrahan-Krueger model with an
+		 isotropic phase function -->
+	<bsdf type="hk"/>
+	
+	<!-- Test the Hanrahan-Krueger model with a 
+		 forward-scattering phase function -->
+	<bsdf type="hk">
+		<phase type="hg">
+			<float name="g" value="0.8"/>
+		</phase>
+	</bsdf>
+
 	<!-- Test the rough glass model with the 
 		 Beckmann microfacet distribution -->
 	<bsdf type="roughdielectric">
@@ -137,18 +149,6 @@
 		<bsdf type="diffuse"/>
 	</bsdf>
 
-	<!-- Test the Hanrahan-Krueger model with an
-		 isotropic phase function -->
-	<bsdf type="hk"/>
-	
-	<!-- Test the Hanrahan-Krueger model with a 
-		 forward-scattering phase function -->
-	<bsdf type="hk">
-		<phase type="hg">
-			<float name="g" value="0.8"/>
-		</phase>
-	</bsdf>
-
 	<!-- Test the smooth coating over a diffuse base material -->
 	<bsdf type="coating">
 		<float name="intIOR" value="1.5"/>

src/bsdfs/hk.cpp

@@ -185,7 +185,6 @@ public:
 			for (int i = 0; i < SPECTRUM_SAMPLES; i++)
 				albedo[i] = sigmaT[i] > 0 ? (sigmaS[i]/sigmaT[i]) : (Float) 0;
 
-
 			const Float cosThetaI = Frame::cosTheta(bRec.wi),
 				        cosThetaO = Frame::cosTheta(bRec.wo),
 				        dp = cosThetaI*cosThetaO;
@@ -202,7 +201,7 @@ public:
 				const Float phaseVal = m_phase->eval(pRec);
 
 				result = albedo * (phaseVal*cosThetaI/(cosThetaI+cosThetaO)) *
-					(Spectrum(1.0f)-((-tauD/std::abs(cosThetaI))+(-tauD/std::abs(cosThetaO))).exp());
+					(Spectrum(1.0f)-((-1.0f/std::abs(cosThetaI)-1.0f/std::abs(cosThetaO)) * tauD).exp());
 			}
 
 			/* ==================================================================== */
@@ -378,6 +377,8 @@ public:
 			<< "]";
 		return oss.str();
 	}
+	
+	Shader *createShader(Renderer *renderer) const;
 
 	MTS_DECLARE_CLASS()
 private:
@@ -387,6 +388,71 @@ private:
 	Float m_thickness;
 };
 
+
+// ================ Hardware shader implementation ================ 
+
+/**
+ * This is a relatively approximate GLSL shader for the HK model.
+ * It assumes that the layer is infinitely thick (i.e. there is no
+ * transmission) and that the phase function is isotropic
+ */
+class HanrahanKruegerShader : public Shader {
+public:
+	HanrahanKruegerShader(Renderer *renderer, const Texture *sigmaS, const Texture *sigmaA) 
+		: Shader(renderer, EBSDFShader), m_sigmaS(sigmaS), m_sigmaA(sigmaA) {
+		m_sigmaSShader = renderer->registerShaderForResource(m_sigmaS.get());
+		m_sigmaAShader = renderer->registerShaderForResource(m_sigmaA.get());
+	}
+
+	bool isComplete() const {
+		return m_sigmaSShader.get() != NULL
+			&& m_sigmaAShader.get() != NULL;
+	}
+
+	void cleanup(Renderer *renderer) {
+		renderer->unregisterShaderForResource(m_sigmaS.get());
+		renderer->unregisterShaderForResource(m_sigmaA.get());
+	}
+
+	void putDependencies(std::vector<Shader *> &deps) {
+		deps.push_back(m_sigmaSShader.get());
+		deps.push_back(m_sigmaAShader.get());
+	}
+
+	void generateCode(std::ostringstream &oss,
+			const std::string &evalName,
+			const std::vector<std::string> &depNames) const {
+		oss << "vec3 " << evalName << "(vec2 uv, vec3 wi, vec3 wo) {" << endl
+			<< "    vec3 sigmaS = " << depNames[0] << "(uv);" << endl
+			<< "    vec3 sigmaA = " << depNames[1] << "(uv);" << endl
+			<< "    vec3 albedo = sigmaS/(sigmaS + sigmaA);" << endl
+			<< "    float cosThetaI = abs(cosTheta(wi));" << endl
+			<< "    float cosThetaO = abs(cosTheta(wo));" << endl
+			<< "    return albedo * (0.079577*cosThetaI*cosThetaO/(cosThetaI + cosThetaO));" << endl
+			<< "}" << endl
+			<< endl
+			<< "vec3 " << evalName << "_diffuse(vec2 uv, vec3 wi, vec3 wo) {" << endl
+			<< "    vec3 sigmaS = " << depNames[0] << "(uv);" << endl
+			<< "    vec3 sigmaA = " << depNames[1] << "(uv);" << endl
+			<< "    vec3 albedo = sigmaS/(sigmaS + sigmaA);" << endl
+			<< "    float cosThetaO = abs(cosTheta(wo));" << endl
+			<< "    return albedo * 0.079577 * cosThetaO;" << endl
+			<< "}" << endl;
+	}
+
+	MTS_DECLARE_CLASS()
+private:
+	ref<const Texture> m_sigmaS;
+	ref<const Texture> m_sigmaA;
+	ref<Shader> m_sigmaSShader;
+	ref<Shader> m_sigmaAShader;
+};
+
+Shader *HanrahanKrueger::createShader(Renderer *renderer) const { 
+	return new HanrahanKruegerShader(renderer, m_sigmaS.get(), m_sigmaA.get());
+}
+
+MTS_IMPLEMENT_CLASS(HanrahanKruegerShader, false, Shader)
 MTS_IMPLEMENT_CLASS_S(HanrahanKrueger, false, BSDF)
 MTS_EXPORT_PLUGIN(HanrahanKrueger, "Hanrahan-Krueger BSDF");
 MTS_NAMESPACE_END

src/bsdfs/microfacet.h

@@ -367,6 +367,10 @@ public:
 							+ alphaV * sinPhiM*sinPhiM;
 					pdf = std::sqrt((alphaU + 1) * (alphaV + 1))
 						* INV_TWOPI * std::pow(cosThetaM, exponent);
+						
+					/* Prevent potential numerical issues in other stages of the model */
+					if (pdf < 1e-20f)
+						pdf = 0;
 					
 					return Vector(
 						sinThetaM * cosPhiM,
@@ -379,6 +383,10 @@ public:
 				SLog(EError, "Invalid distribution function!");
 		}
 
+		/* Prevent potential numerical issues in other stages of the model */
+		if (pdf < 1e-20f)
+			pdf = 0;
+					
 		const Float sinThetaM = std::sqrt(
 			std::max((Float) 0, 1 - cosThetaM*cosThetaM));
 		Float phiM = (2.0f * M_PI) * sample.y;

src/bsdfs/roughconductor.cpp

@@ -302,6 +302,9 @@ public:
 		const Normal m = m_distribution.sample(sample, 
 			alphaU, alphaV, microfacetPDF);
 
+		if (microfacetPDF == 0)
+			return Spectrum(0.0f);
+
 		/* Perfect specular reflection based on the microsurface normal */
 		bRec.wo = reflect(bRec.wi, m);
 		bRec.sampledComponent = 0;
@@ -341,6 +344,9 @@ public:
 		const Normal m = m_distribution.sample(sample, 
 			alphaU, alphaV, pdf);
 
+		if (pdf == 0)
+			return Spectrum(0.0f);
+
 		/* Perfect specular reflection based on the microsurface normal */
 		bRec.wo = reflect(bRec.wi, m);
 		bRec.sampledComponent = 0;

src/bsdfs/roughdielectric.cpp

@@ -464,6 +464,9 @@ public:
 		const Normal m = m_distribution.sample(sample,
 				sampleAlphaU, sampleAlphaV, microfacetPDF);
 
+		if (microfacetPDF == 0)
+			return Spectrum(0.0f);
+
 		Float F = fresnel(dot(bRec.wi, m), m_extIOR, m_intIOR),
 			  numerator = 1.0f;
 
@@ -551,11 +554,8 @@ public:
 		const Normal m = m_distribution.sample(sample,
 				sampleAlphaU, sampleAlphaV, microfacetPDF);
 
-#if 1
+		if (microfacetPDF == 0)
-		Float ref = m_distribution.pdf(m, sampleAlphaU, sampleAlphaV);
+			return Spectrum(0.0f);
-		if (std::abs(ref-microfacetPDF)/ref > Epsilon)
-			cout << "OOPS! ref=" << ref << ", got=" << microfacetPDF << endl;
-#endif
 
 		pdf = microfacetPDF;