more accurate Fresnel reflectance computations involving conductors

include/mitsuba/core/spectrum.h

@@ -493,6 +493,14 @@ public:
 		return value;
 	}
 
+	/// Component-wise square root
+	inline TSpectrum safe_sqrt() const {
+		TSpectrum value;
+		for (int i=0; i<N; i++)
+			value.s[i] = math::safe_sqrt(s[i]);
+		return value;
+	}
+
 	/// Component-wise exponentation
 	inline TSpectrum exp() const {
 		TSpectrum value;

include/mitsuba/core/util.h

@@ -507,22 +507,99 @@ inline Float fresnelDielectricExt(Float cosThetaI, Float eta) { Float cosThetaT;
 	return fresnelDielectricExt(cosThetaI, cosThetaT, eta); }
 
 /**
- * \brief Calculates the unpolarized fresnel reflection coefficient
+ * \brief Calculates the unpolarized Fresnel reflection coefficient
- * at a planar interface between vacuum and a conductor.
+ * at a planar interface having a complex-valued relative index of
+ * refraction (approximate scalar version)
+ *
+ * The implementation of this function relies on a simplified expression
+ * that becomes increasingly accurate as k grows.
+ *
+ * The name of this function is a slight misnomer, since it supports
+ * the general case of a complex-valued relative index of refraction
+ * (rather than being restricted to conductors)
  *
  * \param cosThetaI
  * 		Cosine of the angle between the normal and the incident ray
  * \param eta
- * 		Real refractive index (wavelength-dependent)
+ * 		Relative refractive index (real component)
  * \param k
- * 		Imaginary refractive index (wavelength-dependent)
+ * 		Relative refractive index (imaginary component)
  * \ingroup libpython
  */
-extern MTS_EXPORT_CORE Spectrum fresnelConductor(Float cosThetaI,
+extern MTS_EXPORT_CORE Float fresnelConductorApprox(Float cosThetaI,
+		Float eta, Float k);
+
+/**
+ * \brief Calculates the unpolarized Fresnel reflection coefficient
+ * at a planar interface having a complex-valued relative index of
+ * refraction (approximate vectorized version)
+ *
+ * The implementation of this function relies on a simplified expression
+ * that becomes increasingly accurate as k grows.
+ *
+ * The name of this function is a slight misnomer, since it supports
+ * the general case of a complex-valued relative index of refraction
+ * (rather than being restricted to conductors)
+ *
+ * \param cosThetaI
+ * 		Cosine of the angle between the normal and the incident ray
+ * \param eta
+ * 		Relative refractive index (real component)
+ * \param k
+ * 		Relative refractive index (imaginary component)
+ * \ingroup libpython
+ */
+extern MTS_EXPORT_CORE Spectrum fresnelConductorApprox(Float cosThetaI,
 		const Spectrum &eta, const Spectrum &k);
 
 /**
- * \brief Calculates the diffuse unpolarized fresnel reflectance of
+ * \brief Calculates the unpolarized Fresnel reflection coefficient
+ * at a planar interface having a complex-valued relative index of
+ * refraction (accurate scalar version)
+ *
+ * The implementation of this function computes the exact unpolarized
+ * Fresnel reflectance for a complex index of refraction change.
+ *
+ * The name of this function is a slight misnomer, since it supports
+ * the general case of a complex-valued relative index of refraction
+ * (rather than being restricted to conductors)
+ *
+ * \param cosThetaI
+ * 		Cosine of the angle between the normal and the incident ray
+ * \param eta
+ * 		Relative refractive index (real component)
+ * \param k
+ * 		Relative refractive index (imaginary component)
+ * \ingroup libpython
+ */
+extern MTS_EXPORT_CORE Float fresnelConductorExact(Float cosThetaI,
+		Float eta, Float k);
+
+/**
+ * \brief Calculates the unpolarized Fresnel reflection coefficient
+ * at a planar interface having a complex-valued relative index of
+ * refraction (accurate vectorized version)
+ *
+ * The implementation of this function computes the exact unpolarized
+ * Fresnel reflectance for a complex index of refraction change.
+ *
+ * The name of this function is a slight misnomer, since it supports
+ * the general case of a complex-valued relative index of refraction
+ * (rather than being restricted to conductors)
+ *
+ * \param cosThetaI
+ * 		Cosine of the angle between the normal and the incident ray
+ * \param eta
+ * 		Relative refractive index (real component)
+ * \param k
+ * 		Relative refractive index (imaginary component)
+ * \ingroup libpython
+ */
+extern MTS_EXPORT_CORE Spectrum fresnelConductorExact(Float cosThetaI,
+		const Spectrum &eta, const Spectrum &k);
+
+/**
+ * \brief Calculates the diffuse unpolarized Fresnel reflectance of
  * a dielectric material (sometimes referred to as "Fdr").
  *
  * This value quantifies what fraction of diffuse incident illumination

src/bsdfs/conductor.cpp

@@ -20,6 +20,7 @@
 #include <mitsuba/core/fresolver.h>
 #include <mitsuba/hw/basicshader.h>
 #include <boost/algorithm/string.hpp>
+#include "ior.h"
 
 MTS_NAMESPACE_BEGIN
 
@@ -29,11 +30,12 @@ MTS_NAMESPACE_BEGIN
  * \parameters{
  *     \parameter{material}{\String}{Name of a material preset, see
  *           \tblref{conductor-iors}.\!\default{\texttt{Cu} / copper}}
- *     \parameter{eta}{\Spectrum}{Real part of the material's index
+ *     \parameter{eta, k}{\Spectrum}{Real and imaginary components of the material's index of
- *           of refraction \default{based on the value of \texttt{material}}}
+ *             refraction \default{based on the value of \texttt{material}}}
- *     \parameter{k}{\Spectrum}{Imaginary part of the material's index of
+ *     \parameter{extEta}{\Float\Or\String}{
- *             refraction, also known as absorption coefficient.
+ *           Real-valued index of refraction of the surrounding dielectric,
- *             \default{based on the value of \texttt{material}}}
+ *           or a material name of a dielectric \default{\code{air}}
+ *     }
  *     \parameter{specular\showbreak Reflectance}{\Spectrum\Or\Texture}{Optional
  *         factor that can be used to modulate the specular reflection component. Note
  *         that for physical realism, this parameter should never be touched. \default{1.0}}
@@ -154,21 +156,23 @@ public:
 		m_specularReflectance = new ConstantSpectrumTexture(
 			props.getSpectrum("specularReflectance", Spectrum(1.0f)));
 
-		std::string material = props.getString("material", "Cu");
+		std::string materialName = props.getString("material", "Cu");
 
-		Spectrum materialEta, materialK;
+		Spectrum intEta, intK;
-		if (boost::to_lower_copy(material) == "none") {
+		if (boost::to_lower_copy(materialName) == "none") {
-			materialEta = Spectrum(0.0f);
+			intEta = Spectrum(0.0f);
-			materialK = Spectrum(1.0f);
+			intK = Spectrum(1.0f);
 		} else {
-			materialEta.fromContinuousSpectrum(InterpolatedSpectrum(
+			intEta.fromContinuousSpectrum(InterpolatedSpectrum(
-				fResolver->resolve("data/ior/" + material + ".eta.spd")));
+				fResolver->resolve("data/ior/" + materialName + ".eta.spd")));
-			materialK.fromContinuousSpectrum(InterpolatedSpectrum(
+			intK.fromContinuousSpectrum(InterpolatedSpectrum(
-				fResolver->resolve("data/ior/" + material + ".k.spd")));
+				fResolver->resolve("data/ior/" + materialName + ".k.spd")));
 		}
 
-		m_eta = props.getSpectrum("eta", materialEta);
+		Float extEta = lookupIOR(props, "extEta", "air");
-		m_k = props.getSpectrum("k", materialK);
+
+		m_eta = props.getSpectrum("eta", intEta) / extEta;
+		m_k   = props.getSpectrum("k", intK) / extEta;
 	}
 
 	SmoothConductor(Stream *stream, InstanceManager *manager)
@@ -229,7 +233,7 @@ public:
 			return Spectrum(0.0f);
 
 		return m_specularReflectance->eval(bRec.its) *
-			fresnelConductor(Frame::cosTheta(bRec.wi), m_eta, m_k);
+			fresnelConductorExact(Frame::cosTheta(bRec.wi), m_eta, m_k);
 	}
 
 	Float pdf(const BSDFSamplingRecord &bRec, EMeasure measure) const {
@@ -260,7 +264,7 @@ public:
 		bRec.eta = 1.0f;
 
 		return m_specularReflectance->eval(bRec.its) *
-			fresnelConductor(Frame::cosTheta(bRec.wi), m_eta, m_k);
+			fresnelConductorExact(Frame::cosTheta(bRec.wi), m_eta, m_k);
 	}
 
 	Spectrum sample(BSDFSamplingRecord &bRec, Float &pdf, const Point2 &sample) const {
@@ -277,7 +281,7 @@ public:
 		pdf = 1;
 
 		return m_specularReflectance->eval(bRec.its) *
-			fresnelConductor(Frame::cosTheta(bRec.wi), m_eta, m_k);
+			fresnelConductorExact(Frame::cosTheta(bRec.wi), m_eta, m_k);
 	}
 
 	Float getRoughness(const Intersection &its, int component) const {
@@ -320,7 +324,7 @@ public:
 		m_specularReflectanceShader = renderer->registerShaderForResource(m_specularReflectance.get());
 
 		/* Compute the reflectance at perpendicular incidence */
-		m_R0 = fresnelConductor(1.0f, eta, k);
+		m_R0 = fresnelConductorExact(1.0f, eta, k);
 
 		m_alpha = 0.4f;
 	}

src/bsdfs/roughconductor.cpp

@@ -20,6 +20,7 @@
 #include <mitsuba/render/bsdf.h>
 #include <mitsuba/hw/basicshader.h>
 #include "microfacet.h"
+#include "ior.h"
 
 MTS_NAMESPACE_BEGIN
 
@@ -32,15 +33,15 @@ MTS_NAMESPACE_BEGIN
  *          used to model the surface roughness.
  *       \begin{enumerate}[(i)]
  *           \item \code{beckmann}: Physically-based distribution derived from
- *               Gaussian random surfaces. This is the default.
+ *               Gaussian random surfaces. This is the default.\vspace{-1mm}
  *           \item \code{ggx}: New distribution proposed by
  *              Walter et al. \cite{Walter07Microfacet}, which is meant to better handle
  *              the long tails observed in measurements of ground surfaces.
- *              Renderings with this distribution may converge slowly.
+ *              Renderings with this distribution may converge slowly.\vspace{-1mm}
  *           \item \code{phong}: Classical $\cos^p\theta$ distribution.
  *              Due to the underlying microfacet theory,
  *              the use of this distribution here leads to more realistic
- *              behavior than the separately available \pluginref{phong} plugin.
+ *              behavior than the separately available \pluginref{phong} plugin.\vspace{-1mm}
  *           \item \code{as}: Anisotropic Phong-style microfacet distribution proposed by
  *              Ashikhmin and Shirley \cite{Ashikhmin2005Anisotropic}.\vspace{-3mm}
  *       \end{enumerate}
@@ -59,11 +60,12 @@ MTS_NAMESPACE_BEGIN
  *     }
  *     \parameter{material}{\String}{Name of a material preset, see
  *           \tblref{conductor-iors}.\!\default{\texttt{Cu} / copper}}
- *     \parameter{eta}{\Spectrum}{Real part of the material's index
+ *     \parameter{eta, k}{\Spectrum}{Real and imaginary components of the material's index of
- *           of refraction \default{based on the value of \texttt{material}}}
+ *             refraction \default{based on the value of \texttt{material}}}
- *     \parameter{k}{\Spectrum}{Imaginary part of the material's index of
+ *     \parameter{extEta}{\Float\Or\String}{
- *             refraction (the absorption coefficient).
+ *           Real-valued index of refraction of the surrounding dielectric,
- *             \default{based on \texttt{material}}}
+ *           or a material name of a dielectric \default{\code{air}}
+ *     }
  *     \parameter{specular\showbreak Reflectance}{\Spectrum\Or\Texture}{Optional
  *         factor that can be used to modulate the specular reflection component. Note
  *         that for physical realism, this parameter should never be touched. \default{1.0}}
@@ -158,21 +160,23 @@ public:
 		m_specularReflectance = new ConstantSpectrumTexture(
 			props.getSpectrum("specularReflectance", Spectrum(1.0f)));
 
-		std::string material = props.getString("material", "Cu");
+		std::string materialName = props.getString("material", "Cu");
-		Spectrum materialEta, materialK;
 
-		if (boost::to_lower_copy(material) == "none") {
+		Spectrum intEta, intK;
-			materialEta = Spectrum(0.0f);
+		if (boost::to_lower_copy(materialName) == "none") {
-			materialK = Spectrum(1.0f);
+			intEta = Spectrum(0.0f);
+			intK = Spectrum(1.0f);
 		} else {
-			materialEta.fromContinuousSpectrum(InterpolatedSpectrum(
+			intEta.fromContinuousSpectrum(InterpolatedSpectrum(
-				fResolver->resolve("data/ior/" + material + ".eta.spd")));
+				fResolver->resolve("data/ior/" + materialName + ".eta.spd")));
-			materialK.fromContinuousSpectrum(InterpolatedSpectrum(
+			intK.fromContinuousSpectrum(InterpolatedSpectrum(
-				fResolver->resolve("data/ior/" + material + ".k.spd")));
+				fResolver->resolve("data/ior/" + materialName + ".k.spd")));
 		}
 
-		m_eta = props.getSpectrum("eta", materialEta);
+		Float extEta = lookupIOR(props, "extEta", "air");
-		m_k = props.getSpectrum("k", materialK);
+
+		m_eta = props.getSpectrum("eta", intEta) / extEta;
+		m_k   = props.getSpectrum("k", intK) / extEta;
 
 		m_distribution = MicrofacetDistribution(
 			props.getString("distribution", "beckmann")
@@ -262,7 +266,7 @@ public:
 			return Spectrum(0.0f);
 
 		/* Fresnel factor */
-		const Spectrum F = fresnelConductor(dot(bRec.wi, H), m_eta, m_k);
+		const Spectrum F = fresnelConductorExact(dot(bRec.wi, H), m_eta, m_k);
 
 		/* Smith's shadow-masking function */
 		const Float G = m_distribution.G(bRec.wi, bRec.wo, H, alphaU, alphaV);
@@ -324,7 +328,7 @@ public:
 		if (Frame::cosTheta(bRec.wo) <= 0)
 			return Spectrum(0.0f);
 
-		const Spectrum F = fresnelConductor(dot(bRec.wi, m),
+		const Spectrum F = fresnelConductorExact(dot(bRec.wi, m),
 				m_eta, m_k);
 
 		Float numerator = m_distribution.eval(m, alphaU, alphaV)
@@ -367,7 +371,7 @@ public:
 		if (Frame::cosTheta(bRec.wo) <= 0)
 			return Spectrum(0.0f);
 
-		const Spectrum F = fresnelConductor(dot(bRec.wi, m),
+		const Spectrum F = fresnelConductorExact(dot(bRec.wi, m),
 				m_eta, m_k);
 
 		Float numerator = m_distribution.eval(m, alphaU, alphaV)
@@ -460,7 +464,7 @@ public:
 		m_alphaVShader = renderer->registerShaderForResource(m_alphaV.get());
 
 		/* Compute the reflectance at perpendicular incidence */
-		m_R0 = fresnelConductor(1.0f, eta, k);
+		m_R0 = fresnelConductorExact(1.0f, eta, k);
 	}
 
 	bool isComplete() const {

src/libcore/util.cpp

@@ -586,18 +586,84 @@ Float fresnelDielectricExt(Float cosThetaI_, Float &cosThetaT_, Float eta) {
 	return 0.5f * (Rs * Rs + Rp * Rp);
 }
 
-Spectrum fresnelConductor(Float cosThetaI, const Spectrum &eta, const Spectrum &k) {
+Float fresnelConductorApprox(Float cosThetaI, Float eta, Float k) {
-	Spectrum tmp = (eta*eta + k*k) * (cosThetaI * cosThetaI);
+	Float cosThetaI2 = cosThetaI*cosThetaI;
 
-	Spectrum rParl2 = (tmp - (eta * (2.0f * cosThetaI)) + Spectrum(1.0f))
+	Float tmp = (eta*eta + k*k) * cosThetaI2;
-					/ (tmp + (eta * (2.0f * cosThetaI)) + Spectrum(1.0f));
+
+	Float Rp2 = (tmp - (eta * (2 * cosThetaI)) + 1)
+	          / (tmp + (eta * (2 * cosThetaI)) + 1);
+
+	Float tmpF = eta*eta + k*k;
+
+	Float Rs2 = (tmpF - (eta * (2 * cosThetaI)) + cosThetaI2) /
+	            (tmpF + (eta * (2 * cosThetaI)) + cosThetaI2);
+
+	return 0.5f * (Rp2 + Rs2);
+}
+
+Spectrum fresnelConductorApprox(Float cosThetaI, const Spectrum &eta, const Spectrum &k) {
+	Float cosThetaI2 = cosThetaI*cosThetaI;
+
+	Spectrum tmp = (eta*eta + k*k) * cosThetaI2;
+
+	Spectrum Rp2 = (tmp - (eta * (2 * cosThetaI)) + Spectrum(1.0f))
+	             / (tmp + (eta * (2 * cosThetaI)) + Spectrum(1.0f));
 
 	Spectrum tmpF = eta*eta + k*k;
 
-	Spectrum rPerp2 = (tmpF - (eta * (2.0f * cosThetaI)) + Spectrum(cosThetaI*cosThetaI)) /
+	Spectrum Rs2 = (tmpF - (eta * (2 * cosThetaI)) + Spectrum(cosThetaI2)) /
-					  (tmpF + (eta * (2.0f * cosThetaI)) + Spectrum(cosThetaI*cosThetaI));
+	               (tmpF + (eta * (2 * cosThetaI)) + Spectrum(cosThetaI2));
 
-	return (rParl2 + rPerp2) / 2.0f;
+	return 0.5f * (Rp2 + Rs2);
+}
+
+Float fresnelConductorExact(Float cosThetaI, Float eta, Float k) {
+	/* Modified from "Optics" by K.D. Moeller, University Science Books, 1988 */
+
+	Float cosThetaI2 = cosThetaI*cosThetaI,
+	      sinThetaI2 = 1-cosThetaI2,
+		  sinThetaI4 = sinThetaI2*sinThetaI2;
+
+	Float temp1 = eta*eta - k*k - sinThetaI2,
+	      a2pb2 = math::safe_sqrt(temp1*temp1 + 4*k*k*eta*eta),
+	      a     = math::safe_sqrt(0.5f * (a2pb2 + temp1));
+
+	Float term1 = a2pb2 + cosThetaI2,
+	      term2 = 2*a*cosThetaI;
+
+	Float Rs2 = (term1 - term2) / (term1 + term2);
+
+	Float term3 = a2pb2*cosThetaI2 + sinThetaI4,
+	      term4 = term2*sinThetaI2;
+
+	Float Rp2 = Rs2 * (term3 - term4) / (term3 + term4);
+
+	return 0.5f * (Rp2 + Rs2);
+}
+
+Spectrum fresnelConductorExact(Float cosThetaI, const Spectrum &eta, const Spectrum &k) {
+	/* Modified from "Optics" by K.D. Moeller, University Science Books, 1988 */
+
+	Float cosThetaI2 = cosThetaI*cosThetaI,
+	      sinThetaI2 = 1-cosThetaI2,
+		  sinThetaI4 = sinThetaI2*sinThetaI2;
+
+	Spectrum temp1 = eta*eta - k*k - Spectrum(sinThetaI2),
+	         a2pb2 = (temp1*temp1 + k*k*eta*eta*4).safe_sqrt(),
+	         a     = ((a2pb2 + temp1) * 0.5f).safe_sqrt();
+
+	Spectrum term1 = a2pb2 + Spectrum(cosThetaI2),
+	         term2 = a*(2*cosThetaI);
+
+	Spectrum Rs2 = (term1 - term2) / (term1 + term2);
+
+	Spectrum term3 = a2pb2*cosThetaI2 + Spectrum(sinThetaI4),
+	         term4 = term2*sinThetaI2;
+
+	Spectrum Rp2 = Rs2 * (term3 - term4) / (term3 + term4);
+
+	return 0.5f * (Rp2 + Rs2);
 }
 
 Vector reflect(const Vector &wi, const Normal &n) {

src/libpython/core.cpp

@@ -1261,11 +1261,19 @@ void export_core() {
 		.staticmethod("glOrthographic")
 		.staticmethod("fromFrame");
 
+	Float (*fresnelConductorApprox1)(Float, Float, Float) = &fresnelConductorApprox;
+	Float (*fresnelConductorExact1)(Float, Float, Float) = &fresnelConductorExact;
+	Spectrum (*fresnelConductorApprox2)(Float, const Spectrum &, const Spectrum &) = &fresnelConductorApprox;
+	Spectrum (*fresnelConductorExact2)(Float, const Spectrum &, const Spectrum &) = &fresnelConductorExact;
+
 	/* Functions from utility.h */
 	bp::def("fresnelDielectric", &fresnelDielectric);
 	bp::def("fresnelDielectricExt", &fresnelDielectricExt1);
 	bp::def("fresnelDielectricExt", &fresnelDielectricExt2);
-	bp::def("fresnelConductor", &fresnelConductor, BP_RETURN_VALUE);
+	bp::def("fresnelConductorApprox", fresnelConductorApprox1, BP_RETURN_VALUE);
+	bp::def("fresnelConductorApprox", fresnelConductorApprox2, BP_RETURN_VALUE);
+	bp::def("fresnelConductorExact", fresnelConductorExact1, BP_RETURN_VALUE);
+	bp::def("fresnelConductorExact", fresnelConductorExact2, BP_RETURN_VALUE);
 	bp::def("fresnelDiffuseReflectance", &fresnelDiffuseReflectance);
 	bp::def("reflect", &reflect);
 	bp::def("refract", &refract1);

src/mtsgui/mainwindow.cpp

@@ -1478,9 +1478,6 @@ void MainWindow::on_actionExportImage_triggered() {
 }
 
 void MainWindow::onExportDialogClose(int reason) {
-	int currentIndex = ui->tabBar->currentIndex();
-	SceneContext *ctx = m_context[currentIndex];
-
 	QSettings settings;
 	QFileDialog *dialog = static_cast<QFileDialog *>(sender());
 	m_currentChild = NULL;