added a function for computing the diffuse fresnel reflectance (aka Fdr)

include/mitsuba/core/util.h

@@ -463,6 +463,20 @@ extern MTS_EXPORT_CORE Point2 squareToStdNormal(const Point2 &sample);
 extern MTS_EXPORT_CORE Float fresnelDielectric(Float cosThetaI, 
 		Float cosThetaT, Float etaI, Float etaT);
 
+/**
+ * \brief Calculates the unpolarized fresnel reflection coefficient on
+ * an interface to a conductor.
+ *
+ * \param cosThetaI
+ * 		Cosine of the angle between the normal and the incident ray
+ * \param eta
+ * 		Real refractive index (wavelength-dependent)
+ * \param k
+ * 		Imaginary refractive index (wavelength-dependent)
+ */
+extern MTS_EXPORT_CORE Spectrum fresnelConductor(Float cosThetaI, 
+		const Spectrum &eta, const Spectrum &k);
+
 /**
  * \brief Calculates the unpolarized fresnel reflection coefficient for a 
  * dielectric material. Handles incidence from either sides.
@@ -479,7 +493,7 @@ extern MTS_EXPORT_CORE Float fresnel(Float cosThetaI, Float extIOR,
 
 /**
  * \brief Calculates the diffuse unpolarized fresnel reflectance of
- * a dielectric material. 
+ * a dielectric material (sometimes referred to as "Fdr"). 
  *
  * This value quantifies what fraction of completely diffuse incident 
  * illumination will be reflected by a dielectric material on average.
@@ -495,21 +509,8 @@ extern MTS_EXPORT_CORE Float fresnel(Float cosThetaI, Float extIOR,
  *      a wider range of refraction coefficients, but at a cost
  *      in terms of performance.
  */
-extern MTS_EXPORT_CORE Float diffuseFresnelReflectance(Float eta, bool fast = false);
-
-/**
- * \brief Calculates the unpolarized fresnel reflection coefficient on
- * an interface to a conductor.
- *
- * \param cosThetaI
- * 		Cosine of the angle between the normal and the incident ray
- * \param eta
- * 		Real refractive index (wavelength-dependent)
- * \param k
- * 		Imaginary refractive index (wavelength-dependent)
- */
-extern MTS_EXPORT_CORE Spectrum fresnelConductor(Float cosThetaI, 
-		const Spectrum &eta, const Spectrum &k);
+extern MTS_EXPORT_CORE Float fresnelDiffuseReflectance(
+	Float eta, bool fast = false);
 
 /*! @} */
 

src/bsdfs/dipolebrdf.cpp

@@ -78,8 +78,6 @@ public:
 		if (eta > 1) {
 			return  -1.440f / (eta * eta) + 0.710f / eta + 0.668f + 0.0636f * eta;
 		} else if (eta < 1) {
-			/* Average reflectance due to mismatched indices of refraction
-			 * at the boundary - [Egan et al. 1973] */
 			return  -0.4399f + 0.7099f / eta - 0.3319f / (eta * eta)
 					+ 0.0636f / (eta * eta * eta);
 		} else {
@@ -107,11 +105,10 @@ public:
 		m_usesRayDifferentials = m_sigmaS->usesRayDifferentials()
 			|| m_sigmaA->usesRayDifferentials();
 
-		/* relative index of refraction */
-		const Float eta = m_intIOR / m_extIOR;
+		/* Numerically approximate the diffuse Fresnel reflectance */
+		const Float Fdr = fresnelDiffuse(m_extIOR / m_intIOR, false);
 
-		/* Approximate DipoleBRDF boundary condition term */
-		const Float Fdr = evalFdr(eta);
+		/* Compute the extrapolation distance */
 		m_A = (1 + Fdr) / (1 - Fdr);
 
 		BSDF::configure();

src/libcore/util.cpp

@@ -17,6 +17,8 @@
 */
 
 #include <mitsuba/core/random.h>
+#include <mitsuba/core/quad.h>
+#include <boost/bind.hpp>
 #include <stdarg.h>
 #include <iomanip>
 #include <errno.h>
@@ -584,7 +586,6 @@ void latinHypercube(Random *random, Float *dest, size_t nSamples, size_t nDim) {
 	}
 }
 
-
 Vector sphericalDirection(Float theta, Float phi) {
 	Float sinTheta, cosTheta, sinPhi, cosPhi;
 
@@ -727,6 +728,9 @@ Float lanczosSinc(Float t, Float tau) {
 	by Paul S. Heckbert. */
 Float fresnelDielectric(Float cosThetaI, Float cosThetaT, 
 						Float etaI, Float etaT) {
+	if (etaI == etaT)
+		return 0.0f;
+
 	Float Rs = (etaI * cosThetaI - etaT * cosThetaT)
 			/ (etaI * cosThetaI + etaT * cosThetaT);
 	Float Rp = (etaT * cosThetaI - etaI * cosThetaT)
@@ -772,6 +776,63 @@ Float fresnel(Float cosThetaI, Float extIOR, Float intIOR) {
 		cosThetaT, etaI, etaT);
 }
 
+static Float fresnelDiffuseIntegrand(Float eta, Float xi) {
+	if (eta > 1)
+		return fresnel(std::sqrt(xi), 1, eta);
+	else
+		return fresnel(std::sqrt(xi), 1/eta, 1);
+}
+
+Float fresnelDiffuseReflectance(Float eta, bool fast) {
+	if (fast) {
+		/* Fast mode: the following code approximates the
+		 * diffuse Frensel reflectance for the eta<1 and 
+		 * eta>1 cases. An evalution of the accuracy led
+		 * to the following scheme, which cherry-picks
+		 * fits from two papers where they are best.
+		 */
+		if (eta < 1) {
+			/* Fit by Egan and Hilgeman (1973). Works
+			   reasonably well for "normal" IOR values (<2).
+	
+			   Max rel. error in 1.0 - 1.5 : 0.1%
+			   Max rel. error in 1.5 - 2   : 0.6%
+			   Max rel. error in 2.0 - 5   : 9.5%
+			*/
+			return -1.4399f * (eta * eta) 
+				  + 0.7099f * eta 
+				  + 0.6681f 
+				  + 0.0636f / eta;
+		} else {
+			/* Fit by d'Eon and Irving (2011)
+			 *
+			 * Maintains a good accuracy even for 
+			 * unrealistic IOR values.
+			 *
+			 * Max rel. error in 1.0 - 2.0   : 0.1%
+			 * Max rel. error in 2.0 - 10.0  : 0.2%
+			 */
+			Float invEta = 1.0f / eta,
+				  invEta2 = invEta*invEta,
+				  invEta3 = invEta2*invEta,
+				  invEta4 = invEta3*invEta,
+				  invEta5 = invEta4*invEta;
+
+			return 0.919317f - 3.4793f * invEta 
+				 + 6.75335f * invEta2
+				 - 7.80989f * invEta3 
+				 + 4.98554f * invEta4 
+				 - 1.36881f * invEta5;
+		}
+	} else {
+		GaussLobattoIntegrator quad(1024, 0, 1e-5f);
+		return quad.integrate(
+			boost::bind(&fresnelDiffuseIntegrand, eta, _1), 0, 1);
+	}
+
+	return 0.0f;
+}
+
 Float radicalInverse(int b, size_t i) {
 	Float invB = (Float) 1 / (Float) b;
 	Float x = 0.0f, f = invB;

src/libpython/core.cpp

@@ -1061,6 +1061,12 @@ void export_core() {
 		.staticmethod("glOrthographic")
 		.staticmethod("fromFrame");
 
+	/* Functions from utility.h */
+	bp::def("fresnel", &fresnel);
+	bp::def("fresnelDielectric", &fresnelDielectric);
+	bp::def("fresnelConductor", &fresnelConductor, BP_RETURN_VALUE);
+	bp::def("fresnelDiffuseReflectance", &fresnelDiffuseReflectance);
+
 	bp::detail::current_scope = oldScope;
 }