From d02d4ebedf1e3dae68710f02a4c3ea5c882c8d22 Mon Sep 17 00:00:00 2001
From: Wenzel Jakob <wenzel@inf.ethz.ch>
Date: Thu, 31 Jul 2014 22:33:00 +0200
Subject: [PATCH] ported all microfacet models to a new visible normal sampling
 implementation

---
 doc/gendoc.py                 |   4 +-
 doc/main.bib                  |  21 +
 src/bsdfs/microfacet.h        | 883 ++++++++++++++++++++--------------
 src/bsdfs/roughcoating.cpp    | 114 +++--
 src/bsdfs/roughconductor.cpp  | 305 ++++++------
 src/bsdfs/roughdielectric.cpp | 352 +++++++-------
 src/bsdfs/roughplastic.cpp    | 147 +++---
 src/tests/test_microfacet.cpp | 176 +++++++
 8 files changed, 1218 insertions(+), 784 deletions(-)
 create mode 100644 src/tests/test_microfacet.cpp

diff --git a/doc/gendoc.py b/doc/gendoc.py
index d83aed88..da7a9fc8 100755
--- a/doc/gendoc.py
+++ b/doc/gendoc.py
@@ -78,8 +78,8 @@ def texify(texfile):
 	else:
 		check_call(['pdflatex', texfile])
 		check_call(['bibtex',   texfile.replace('.tex', '.aux')])
-		check_call(['pdflatex', texfile])
-		check_call(['pdflatex', texfile])
+		#check_call(['pdflatex', texfile])
+		#check_call(['pdflatex', texfile])
 
 os.chdir(os.path.dirname(os.path.abspath(__file__)))
 with open('plugins_generated.tex', 'w') as f:
diff --git a/doc/main.bib b/doc/main.bib
index 094f0075..4adf4cce 100644
--- a/doc/main.bib
+++ b/doc/main.bib
@@ -563,3 +563,24 @@
 	year = 2008,
 	pages = {183--190}
 }
+
+@article{Trowbridge19975Average,
+	author = {T. S. Trowbridge and K. P. Reitz},
+	journal = {J. Opt. Soc. Am.},
+	number = {5},
+	pages = {531--536},
+	publisher = {OSA},
+	title = {Average irregularity representation of a rough surface for ray reflection},
+	volume = {65},
+	month = {May},
+	year = {1975}
+}
+
+@article{Heitz1014Importance,
+  author = {Heitz, Eric and D'Eon, Eugene},
+  title = {Importance Sampling Microfacet-Based BSDFs using the Distribution of Visible Normals},
+  journal = {Computer Graphics Forum},
+  publisher = {Blackwell Publishing},
+  year = {2014},
+  month = Jun
+}
diff --git a/src/bsdfs/microfacet.h b/src/bsdfs/microfacet.h
index a83df449..2c8dc419 100644
--- a/src/bsdfs/microfacet.h
+++ b/src/bsdfs/microfacet.h
@@ -19,16 +19,28 @@
 #if !defined(__MICROFACET_H)
 #define __MICROFACET_H
 
-#include <mitsuba/core/timer.h>
+#include <mitsuba/mitsuba.h>
+#include <mitsuba/core/frame.h>
+#include <mitsuba/core/properties.h>
 #include <boost/algorithm/string.hpp>
-#include <boost/bind.hpp>
 
 MTS_NAMESPACE_BEGIN
 
 /**
- * Implements the microfacet distributions discussed in
- * "Microfacet Models for Refraction through Rough Surfaces"
- * by Bruce Walter, Stephen R. Marschner, Hongsong Li, and Kenneth E. Torrance
+ * \brief Implementation of the Beckman and GGX / Trowbridge-Reitz microfacet
+ * distributions and various useful sampling routines
+ *
+ * Based on the papers
+ *
+ *   "Microfacet Models for Refraction through Rough Surfaces"
+ *    by Bruce Walter, Stephen R. Marschner, Hongsong Li, and Kenneth E. Torrance
+ *
+ * and
+ *
+ *   "Importance Sampling Microfacet-Based BSDFs using the Distribution of Visible Normals"
+ *    by Eric Heitz and Eugene D'Eon
+ *
+ *  The visible normal sampling code was provided by Eric Heitz and Eugene D'Eon.
  */
 class MicrofacetDistribution {
 public:
@@ -36,264 +48,332 @@ public:
 	enum EType {
 		/// Beckmann distribution derived from Gaussian random surfaces
 		EBeckmann         = 0,
-		/// Long-tailed distribution proposed by Walter et al.
+
+		/// GGX: Long-tailed distribution for very rough surfaces (aka. Trowbridge-Reitz distr.)
 		EGGX              = 1,
-		/// Classical Phong distribution
-		EPhong            = 2,
-		/// Anisotropic distribution by Ashikhmin and Shirley
-		EAshikhminShirley = 3
+
+		/// Phong distribution (with the anisotropic extension by Ashikhmin and Shirley)
+		EPhong            = 2
 	};
 
-	/// Create a microfacet distribution of the specified type
-	MicrofacetDistribution(EType type = EBeckmann)
-		: m_type(type) { }
+	/**
+	 * Create an isotropic microfacet distribution of the specified type
+	 *
+	 * \param type
+	 *     The desired type of microfacet distribution
+	 * \param alpha
+	 *     The surface roughness
+	 */
+	inline MicrofacetDistribution(EType type, Float alpha, bool sampleVisible = true)
+		: m_type(type), m_alphaU(alpha), m_alphaV(alpha), m_sampleVisible(sampleVisible),
+	      m_exponentU(0.0f), m_exponentV(0.0f) {
+		if (m_type == EPhong)
+			computePhongExponent();
+	}
 
 	/**
-	 * \brief Create a microfacet distribution of the specified name
-	 * (ggx/phong/beckmann/as)
+	 * Create an anisotropic microfacet distribution of the specified type
+	 *
+	 * \param type
+	 *     The desired type of microfacet distribution
+	 * \param alphaU
+	 *     The surface roughness in the tangent direction
+	 * \param alphaV
+	 *     The surface roughness in the bitangent direction
 	 */
-	MicrofacetDistribution(const std::string &name) : m_type(EBeckmann) {
-		std::string distr = boost::to_lower_copy(name);
+	inline MicrofacetDistribution(EType type, Float alphaU, Float alphaV, bool sampleVisible = true)
+		: m_type(type), m_alphaU(alphaU), m_alphaV(alphaV), m_sampleVisible(sampleVisible),
+	      m_exponentU(0.0f), m_exponentV(0.0f) {
+		if (m_type == EPhong)
+			computePhongExponent();
+	}
 
-		if (distr == "beckmann")
-			m_type = EBeckmann;
-		else if (distr == "phong")
-			m_type = EPhong;
-		else if (distr == "ggx")
-			m_type = EGGX;
-		else if (distr == "as")
-			m_type = EAshikhminShirley;
-		else
-			SLog(EError, "Specified an invalid distribution \"%s\", must be "
-				"\"beckmann\", \"phong\", \"ggx\", or \"as\"!", distr.c_str());
+	/**
+	 * \brief Create a microfacet distribution from a Property data
+	 * structure
+	 */
+	MicrofacetDistribution(const Properties &props, EType type = EBeckmann,
+		Float alphaU = 0.1f, Float alphaV = 0.1f, bool sampleVisible = true)
+		: m_type(type), m_alphaU(alphaU), m_alphaV(alphaV), m_exponentU(0.0f),
+		  m_exponentV(0.0f) {
+
+		if (props.hasProperty("distribution")) {
+			std::string distr = boost::to_lower_copy(props.getString("distribution"));
+			if (distr == "beckmann")
+				m_type = EBeckmann;
+			else if (distr == "ggx")
+				m_type = EGGX;
+			else if (distr == "phong" || distr == "as")
+				m_type = EPhong;
+			else
+				SLog(EError, "Specified an invalid distribution \"%s\", must be "
+					"\"beckmann\", \"ggx\", or \"phong\"/\"as\"!", distr.c_str());
+		}
+
+		if (props.hasProperty("alpha")) {
+			m_alphaU = m_alphaV = props.getFloat("alpha");
+			if (props.hasProperty("alphaU") || props.hasProperty("alphaV"))
+				SLog(EError, "Microfacet model: please specify either 'alpha' or 'alphaU'/'alphaV'.");
+		} else if (props.hasProperty("alphaU") || props.hasProperty("alphaV")) {
+			if (!props.hasProperty("alphaU") || !props.hasProperty("alphaV"))
+				SLog(EError, "Microfacet model: both 'alphaU' and 'alphaV' must be specified.");
+			if (props.hasProperty("alpha"))
+				SLog(EError, "Microfacet model: please specify either 'alpha' or 'alphaU'/'alphaV'.");
+			m_alphaU = props.getFloat("alphaU");
+			m_alphaV = props.getFloat("alphaV");
+		}
+
+		m_sampleVisible = props.getBoolean("sampleVisible", sampleVisible);
+
+		/* Visible normal sampling is not supported for the Phong / Ashikhmin-Shirley distribution */
+		if (m_type == EPhong) {
+			m_sampleVisible = false;
+			computePhongExponent();
+		}
 	}
 
 	/// Return the distribution type
 	inline EType getType() const { return m_type; }
 
+	/// Return the roughness (isotropic case)
+	inline Float getAlpha() const { return m_alphaU; }
+
+	/// Return the roughness along the tangent direction
+	inline Float getAlphaU() const { return m_alphaU; }
+
+	/// Return the roughness along the bitangent direction
+	inline Float getAlphaV() const { return m_alphaV; }
+
+	/// Return the Phong exponent (isotropic case)
+	inline Float getExponent() const { return m_exponentU; }
+
+	/// Return the Phong exponent along the tangent direction
+	inline Float getExponentU() const { return m_exponentU; }
+
+	/// Return the Phong exponent along the bitangent direction
+	inline Float getExponentV() const { return m_exponentV; }
+
+	/// Return whether or not only visible normals are sampled?
+	inline bool getSampleVisible() const { return m_sampleVisible; }
+
 	/// Is this an anisotropic microfacet distribution?
-	bool isAnisotropic() const {
-		return m_type == EAshikhminShirley;
+	inline bool isAnisotropic() const { return m_alphaU != m_alphaV; }
+
+	/// Is this an anisotropic microfacet distribution?
+	inline bool isIsotropic() const { return m_alphaU == m_alphaV; }
+
+	/// Scale the roughness values by some constant
+	inline void scaleAlpha(Float value) {
+		m_alphaU *= value;
+		m_alphaV *= value;
+		if (m_type == EPhong)
+			computePhongExponent();
 	}
 
 	/**
-	 * \brief Convert the roughness values so that they behave similarly to the
-	 * Beckmann distribution.
+	 * \brief Evaluate the microfacet distribution function
 	 *
-	 * Also clamps to the minimal exponent to to avoid numerical issues
-	 * (For lower roughness values, please switch to the smooth BSDF variants)
+	 * \param m
+	 *     The microfacet normal
 	 */
-	Float transformRoughness(Float value) const {
-		#if defined(SINGLE_PRECISION)
-			value = std::max(value, (Float) 1e-3f);
-		#else
-			value = std::max(value, (Float) 1e-5f);
-		#endif
-		if (m_type == EPhong || m_type == EAshikhminShirley)
-			value = std::max(2 / (value * value) - 2, (Float) 0.1f);
-		return value;
-	}
-
-	/**
-	 * \brief Implements the microfacet distribution function D
-	 *
-	 * \param m The microsurface normal
-	 * \param alpha  The surface roughness
-	 */
-	inline Float eval(const Vector &m, Float alpha) const {
-		return eval(m, alpha, alpha);
-	}
-
-	/**
-	 * \brief Implements the microfacet distribution function D
-	 *
-	 * \param m The microsurface normal
-	 * \param alphaU  The surface roughness in the tangent direction
-	 * \param alphaV  The surface roughness in the bitangent direction
-	 */
-	Float eval(const Vector &m, Float alphaU, Float alphaV) const {
+	inline Float eval(const Vector &m) const {
 		if (Frame::cosTheta(m) <= 0)
 			return 0.0f;
 
+		Float cosTheta2 = Frame::cosTheta2(m);
+		Float beckmannExponent = ((m.x*m.x) / (m_alphaU * m_alphaU)
+				+ (m.y*m.y) / (m_alphaV * m_alphaV)) / cosTheta2;
+
 		Float result;
 		switch (m_type) {
 			case EBeckmann: {
-					/* Beckmann distribution function for Gaussian random surfaces */
-					const Float exponent  = Frame::tanTheta2(m) / (alphaU*alphaU);
-					const Float cosTheta2 = Frame::cosTheta2(m);
-
-					result = math::fastexp(-exponent) /
-						(M_PI * alphaU*alphaU*cosTheta2*cosTheta2);
+					/* Beckmann distribution function for Gaussian random surfaces - [Walter 2005] evaluation */
+					result = math::fastexp(-beckmannExponent) /
+						(M_PI * m_alphaU * m_alphaV * cosTheta2 * cosTheta2);
 				}
 				break;
 
 			case EGGX: {
-					/* Empirical GGX distribution function for rough surfaces */
-					const Float tanTheta2 = Frame::tanTheta2(m),
-						        cosTheta2 = Frame::cosTheta2(m);
-
-					const Float root = alphaU / (cosTheta2 *
-								(alphaU*alphaU + tanTheta2));
-
-					result = INV_PI * (root * root);
+					/* GGX / Trowbridge-Reitz distribution function for rough surfaces */
+					Float root = ((Float) 1 + beckmannExponent) * cosTheta2;
+					result = (Float) 1 / (M_PI * m_alphaU * m_alphaV * root * root);
 				}
 				break;
 
 			case EPhong: {
-					/* Phong distribution function */
-					result = (alphaU + 2) * INV_TWOPI
-							* std::pow(Frame::cosTheta(m), alphaU);
+					/* Isotropic case: Phong distribution. Anisotropic case: Ashikhmin-Shirley distribution */
+					Float exponent = interpolatePhongExponent(m);
+					result = std::sqrt((m_exponentU + 2) * (m_exponentV + 2))
+						* INV_TWOPI * std::pow(Frame::cosTheta(m), exponent);
 				}
 				break;
 
-			case EAshikhminShirley: {
-					const Float cosTheta = Frame::cosTheta(m);
-					const Float ds = 1 - cosTheta * cosTheta;
-					if (ds < 0)
-						return 0.0f;
-					const Float exponent = (alphaU * m.x * m.x
-							+ alphaV * m.y * m.y) / ds;
-					result = std::sqrt((alphaU + 2) * (alphaV + 2))
-						* INV_TWOPI * std::pow(cosTheta, exponent);
-				}
-				break;
 
 			default:
-				SLog(EError, "Invalid distribution function!");
-				return 0.0f;
+				SLog(EError, "Invalid distribution type!");
+				return -1;
 		}
 
 		/* Prevent potential numerical issues in other stages of the model */
-		if (result < 1e-20f)
+		if (result * Frame::cosTheta(m) < 1e-20f)
 			result = 0;
 
 		return result;
 	}
 
 	/**
-	 * \brief Returns the density function associated with
-	 * the \ref sample() function.
-	 * \param m The microsurface normal
-	 * \param alpha  The surface roughness
+	 * \brief Wrapper function which calls \ref sampleAll() or \ref sampleVisible()
+	 * depending on the parameters of this class
 	 */
-	inline Float pdf(const Vector &m, Float alpha) const {
-		return pdf(m, alpha, alpha);
+	inline Normal sample(const Vector &wi, const Point2 &sample, Float &pdf) const {
+		Normal m;
+		if (m_sampleVisible) {
+			m = sampleVisible(wi, sample);
+			pdf = pdfVisible(wi, m);
+		} else {
+			m = sampleAll(sample, pdf);
+		}
+		return m;
 	}
 
 	/**
-	 * \brief Returns the density function associated with
-	 * the \ref sample() function.
-	 * \param m The microsurface normal
-	 * \param alphaU  The surface roughness in the tangent direction
-	 * \param alphaV  The surface roughness in the bitangent direction
+	 * \brief Wrapper function which calls \ref sampleAll() or \ref sampleVisible()
+	 * depending on the parameters of this class
 	 */
-	Float pdf(const Vector &m, Float alphaU, Float alphaV) const {
-		/* Usually, this is just D(m) * cos(theta_M) */
-		if (m_type != EAshikhminShirley)
-			return eval(m, alphaU, alphaV) * Frame::cosTheta(m);
-
-		/* For the Ashikhmin-Shirley model, the sampling density
-		   does not include the cos(theta_M) factor, and the
-		   normalization is slightly different than in eval(). */
-		const Float cosTheta = Frame::cosTheta(m);
-		const Float ds = 1 - cosTheta * cosTheta;
-		if (ds < 0)
-			return 0.0f;
-		const Float exponent = (alphaU * m.x * m.x
-				+ alphaV * m.y * m.y) / ds;
-		Float result = std::sqrt((alphaU + 1) * (alphaV + 1))
-			* INV_TWOPI * std::pow(cosTheta, exponent);
-
-		/* Prevent potential numerical issues in other stages of the model */
-		if (result < 1e-20f)
-			result = 0;
-
-		return result;
+	inline Normal sample(const Vector &wi, const Point2 &sample) const {
+		Normal m;
+		if (m_sampleVisible) {
+			m = sampleVisible(wi, sample);
+		} else {
+			Float pdf;
+			m = sampleAll(sample, pdf);
+		}
+		return m;
 	}
 
-	/// Helper routine: sample the first quadrant of the A&S distribution
-	void sampleFirstQuadrant(Float alphaU, Float alphaV, Float u1, Float u2,
-			Float &phi, Float &cosTheta) const {
-		if (alphaU == alphaV)
-			phi = M_PI * u1 * 0.5f;
+	/**
+	 * \brief Wrapper function which calls \ref pdfAll() or \ref pdfVisible()
+	 * depending on the parameters of this class
+	 */
+	inline Float pdf(const Vector &wi, const Vector &m) const {
+		if (m_sampleVisible)
+			return pdfVisible(wi, m);
 		else
-			phi = std::atan(
-				std::sqrt((alphaU + 1.0f) / (alphaV + 1.0f)) *
-				std::tan(M_PI * u1 * 0.5f));
-		const Float cosPhi = std::cos(phi), sinPhi = std::sin(phi);
-		cosTheta = std::pow(u2, 1.0f /
-			(alphaU * cosPhi * cosPhi + alphaV * sinPhi * sinPhi + 1.0f));
+			return pdfAll(m);
 	}
 
 	/**
-	 * \brief Draw a sample from the microsurface normal distribution
+	 * \brief Draw a sample from the microfacet normal distribution
+	 * (including *all* normals) and return the associated
+	 * probability density
 	 *
-	 * \param sample  A uniformly distributed 2D sample
-	 * \param alpha  The surface roughness
+	 * \param sample
+	 *    A uniformly distributed 2D sample
+	 * \param pdf
+	 *    The probability density wrt. solid angles
 	 */
-	inline Normal sample(const Point2 &sample, Float alpha) const {
-		return MicrofacetDistribution::sample(sample, alpha, alpha);
-	}
-
-	/**
-	 * \brief Draw a sample from the microsurface normal distribution
-	 *
-	 * \param sample  A uniformly distributed 2D sample
-	 * \param alphaU  The surface roughness in the tangent direction
-	 * \param alphaV  The surface roughness in the bitangent direction
-	 */
-	Normal sample(const Point2 &sample, Float alphaU, Float alphaV) const {
+	inline Normal sampleAll(const Point2 &sample, Float &pdf) const {
 		/* The azimuthal component is always selected
 		   uniformly regardless of the distribution */
-		Float cosThetaM = 0.0f, phiM = (2.0f * M_PI) * sample.y;
+		Float cosThetaM = 0.0f;
+		Float sinPhiM, cosPhiM;
+		Float alphaSqr;
 
 		switch (m_type) {
 			case EBeckmann: {
-					Float tanThetaMSqr = -alphaU*alphaU * math::fastlog(1.0f - sample.x);
-					cosThetaM = 1.0f / std::sqrt(1 + tanThetaMSqr);
+					/* Beckmann distribution function for Gaussian random surfaces */
+					if (isIsotropic()) {
+						/* Sample phi component (isotropic case) */
+						math::sincos((2.0f * M_PI) * sample.y, &sinPhiM, &cosPhiM);
+
+						alphaSqr = m_alphaU * m_alphaU;
+					} else {
+						/* Sample phi component (anisotropic case) */
+						Float phiM = std::atan(m_alphaV / m_alphaU *
+							std::tan(M_PI + 2*M_PI*sample.y)) + M_PI * std::floor(2*sample.y + 0.5f);
+						math::sincos(phiM, &sinPhiM, &cosPhiM);
+
+						Float cosSc = cosPhiM / m_alphaU, sinSc = sinPhiM / m_alphaV;
+						alphaSqr = 1.0f / (cosSc*cosSc + sinSc*sinSc);
+					}
+
+					/* Sample theta component */
+					Float tanThetaMSqr = alphaSqr * -math::fastlog(1.0f - sample.x);
+					cosThetaM = 1.0f / std::sqrt(1.0f + tanThetaMSqr);
+
+					/* Compute probability density of the sampled position */
+					pdf = (1.0f - sample.x) / (M_PI*m_alphaU*m_alphaV*cosThetaM*cosThetaM*cosThetaM);
 				}
 				break;
 
 			case EGGX: {
-					Float tanThetaMSqr = alphaU * alphaU * sample.x / (1.0f - sample.x);
-					cosThetaM = 1.0f / std::sqrt(1 + tanThetaMSqr);
+					/* GGX / Trowbridge-Reitz distribution function for rough surfaces */
+					if (isIsotropic()) {
+						/* Sample phi component (isotropic case) */
+						math::sincos((2.0f * M_PI) * sample.y, &sinPhiM, &cosPhiM);
+
+						/* Sample theta component */
+						alphaSqr = m_alphaU*m_alphaU;
+					} else {
+						/* Sample phi component (anisotropic case) */
+						Float phiM = std::atan(m_alphaV / m_alphaU *
+							std::tan(M_PI + 2*M_PI*sample.y)) + M_PI * std::floor(2*sample.y + 0.5f);
+						math::sincos(phiM, &sinPhiM, &cosPhiM);
+
+						Float cosSc = cosPhiM / m_alphaU, sinSc = sinPhiM / m_alphaV;
+						alphaSqr = 1.0f / (cosSc*cosSc + sinSc*sinSc);
+					}
+
+					/* Sample theta component */
+					Float tanThetaMSqr = alphaSqr * sample.x / (1.0f - sample.x);
+					cosThetaM = 1.0f / std::sqrt(1.0f + tanThetaMSqr);
+
+					/* Compute probability density of the sampled position */
+					Float temp = 1+tanThetaMSqr/alphaSqr;
+					pdf = INV_PI / (m_alphaU*m_alphaV*cosThetaM*cosThetaM*cosThetaM*temp*temp);
 				}
 				break;
 
 			case EPhong: {
-					cosThetaM = std::pow(sample.x, 1 / (alphaU + 2));
+					Float phiM;
+					Float exponent;
+					if (isIsotropic()) {
+						phiM = (2.0f * M_PI) * sample.y;
+						exponent = m_exponentU;
+					} else {
+						/* Sampling method based on code from PBRT */
+						if (sample.y < 0.25f) {
+							sampleFirstQuadrant(4 * sample.y, phiM, exponent);
+						} else if (sample.y < 0.5f) {
+							sampleFirstQuadrant(4 * (0.5f - sample.y), phiM, exponent);
+							phiM = M_PI - phiM;
+						} else if (sample.y < 0.75f) {
+							sampleFirstQuadrant(4 * (sample.y - 0.5f), phiM, exponent);
+							phiM += M_PI;
+						} else {
+							sampleFirstQuadrant(4 * (1 - sample.y), phiM, exponent);
+							phiM = 2 * M_PI - phiM;
+						}
+					}
+					math::sincos(phiM, &sinPhiM, &cosPhiM);
+					cosThetaM = std::pow(sample.x, 1.0f / (exponent + 2.0f));
+					pdf = std::sqrt((m_exponentU + 2.0f) * (m_exponentV + 2.0f))
+						* INV_TWOPI * std::pow(cosThetaM, exponent + 1.0f);
 				}
 				break;
 
-			case EAshikhminShirley: {
-					/* Sampling method based on code from PBRT */
-					if (sample.x < 0.25f) {
-						sampleFirstQuadrant(alphaU, alphaV,
-							4 * sample.x, sample.y, phiM, cosThetaM);
-					} else if (sample.x < 0.5f) {
-						sampleFirstQuadrant(alphaU, alphaV,
-							4 * (0.5f - sample.x), sample.y, phiM, cosThetaM);
-						phiM = M_PI - phiM;
-					} else if (sample.x < 0.75f) {
-						sampleFirstQuadrant(alphaU, alphaV,
-							4 * (sample.x - 0.5f), sample.y, phiM, cosThetaM);
-						phiM += M_PI;
-					} else {
-						sampleFirstQuadrant(alphaU, alphaV,
-							4 * (1 - sample.x), sample.y, phiM, cosThetaM);
-						phiM = 2 * M_PI - phiM;
-					}
-				}
-				break;
 			default:
-				SLog(EError, "Invalid distribution function!");
+				SLog(EError, "Invalid distribution type!");
+				pdf = -1;
+				return Vector(-1);
 		}
 
-		const Float sinThetaM = std::sqrt(
-			std::max((Float) 0, 1 - cosThetaM*cosThetaM));
+		/* Prevent potential numerical issues in other stages of the model */
+		if (pdf < 1e-20f)
+			pdf = 0;
 
-		Float sinPhiM, cosPhiM;
-		math::sincos(phiM, &sinPhiM, &cosPhiM);
+		Float sinThetaM = std::sqrt(
+			std::max((Float) 0, 1 - cosThetaM*cosThetaM));
 
 		return Vector(
 			sinThetaM * cosPhiM,
@@ -302,224 +382,327 @@ public:
 		);
 	}
 
+	/**
+	 * \brief Returns the density function associated with
+	 * the \ref sampleAll() function.
+	 *
+	 * \param m
+	 *     The microfacet normal
+	 */
+	inline Float pdfAll(const Vector &m) const {
+		/* PDF is just D(m) * cos(theta_M) */
+		return eval(m) * Frame::cosTheta(m);
+	}
+
 
 	/**
-	 * \brief Draw a sample from the microsurface normal distribution
+	 * \brief Draw a sample from the distribution of visible normals
 	 * and return the associated probability density
 	 *
-	 * \param sample  A uniformly distributed 2D sample
-	 * \param alphaU  The surface roughness in the tangent direction
-	 * \param alphaV  The surface roughness in the bitangent direction
-	 * \param pdf     The probability density wrt. solid angles
+	 * \param _wi
+	 *    A reference direction that defines the set of visible normals
+	 * \param sample
+	 *    A uniformly distributed 2D sample
+	 * \param pdf
+	 *    The probability density wrt. solid angles
 	 */
-	Normal sample(const Point2 &sample, Float alphaU, Float alphaV, Float &pdf) const {
-		/* The azimuthal component is always selected
-		   uniformly regardless of the distribution */
-		Float cosThetaM = 0.0f;
+	inline Normal sampleVisible(const Vector &_wi, const Point2 &sample) const {
+		/* Step 1: stretch wi */
+		Vector wi = normalize(Vector(
+			m_alphaU * _wi.x,
+			m_alphaV * _wi.y,
+			_wi.z
+		));
 
-		switch (m_type) {
-			case EBeckmann: {
-					Float tanThetaMSqr = -alphaU*alphaU * math::fastlog(1.0f - sample.x);
-					cosThetaM = 1.0f / std::sqrt(1 + tanThetaMSqr);
-					Float cosThetaM2 = cosThetaM * cosThetaM,
-						  cosThetaM3 = cosThetaM2 * cosThetaM;
-					pdf = (1.0f - sample.x) / (M_PI * alphaU*alphaU * cosThetaM3);
-				}
-				break;
-
-			case EGGX: {
-					Float alphaUSqr = alphaU * alphaU;
-					Float tanThetaMSqr = alphaUSqr * sample.x / (1.0f - sample.x);
-					cosThetaM = 1.0f / std::sqrt(1 + tanThetaMSqr);
-
-					Float cosThetaM2 = cosThetaM * cosThetaM,
-						  cosThetaM3 = cosThetaM2 * cosThetaM,
-						  temp = alphaUSqr + tanThetaMSqr;
-
-					pdf = INV_PI * alphaUSqr / (cosThetaM3 * temp * temp);
-				}
-				break;
-
-			case EPhong: {
-					Float exponent = 1 / (alphaU + 2);
-					cosThetaM = std::pow(sample.x, exponent);
-					pdf = (alphaU + 2) * INV_TWOPI * std::pow(sample.x, (alphaU+1) * exponent);
-				}
-				break;
-
-			case EAshikhminShirley: {
-					Float phiM;
-
-					/* Sampling method based on code from PBRT */
-					if (sample.x < 0.25f) {
-						sampleFirstQuadrant(alphaU, alphaV,
-							4 * sample.x, sample.y, phiM, cosThetaM);
-					} else if (sample.x < 0.5f) {
-						sampleFirstQuadrant(alphaU, alphaV,
-							4 * (0.5f - sample.x), sample.y, phiM, cosThetaM);
-						phiM = M_PI - phiM;
-					} else if (sample.x < 0.75f) {
-						sampleFirstQuadrant(alphaU, alphaV,
-							4 * (sample.x - 0.5f), sample.y, phiM, cosThetaM);
-						phiM += M_PI;
-					} else {
-						sampleFirstQuadrant(alphaU, alphaV,
-							4 * (1 - sample.x), sample.y, phiM, cosThetaM);
-						phiM = 2 * M_PI - phiM;
-					}
-					const Float sinThetaM = std::sqrt(
-						    std::max((Float) 0, 1 - cosThetaM*cosThetaM)),
-						  sinPhiM = std::sin(phiM),
-						  cosPhiM = std::cos(phiM);
-
-					const Float exponent = alphaU * cosPhiM*cosPhiM
-							+ 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,
-						sinThetaM * sinPhiM,
-						cosThetaM
-					);
-				}
-				break;
-			default:
-				pdf = 0;
-				SLog(EError, "Invalid distribution function!");
+		/* Get polar coordinates */
+		Float theta = 0, phi = 0;
+		if (wi.z < (Float) 0.99999) {
+			theta = std::acos(wi.z);
+			phi = std::atan2(wi.y, wi.x);
 		}
+		Float sinPhi, cosPhi;
+		math::sincos(phi, &sinPhi, &cosPhi);
 
-		/* Prevent potential numerical issues in other stages of the model */
-		if (pdf < 1e-20f)
-			pdf = 0;
+		/* Step 2: simulate P22_{wi}(slope.x, slope.y, 1, 1) */
+		Vector2 slope = sampleVisible11(theta, sample);
 
-		const Float sinThetaM = std::sqrt(
-			std::max((Float) 0, 1 - cosThetaM*cosThetaM));
-		Float phiM = (2.0f * M_PI) * sample.y;
-		return Vector(
-			sinThetaM * std::cos(phiM),
-			sinThetaM * std::sin(phiM),
-			cosThetaM
+		/* Step 3: rotate */
+		slope = Vector2(
+			cosPhi * slope.x - sinPhi * slope.y,
+			sinPhi * slope.x + cosPhi * slope.y);
+
+		/* Step 4: unstretch */
+		slope.x *= m_alphaU;
+		slope.y *= m_alphaV;
+
+		/* Step 5: compute normal */
+		Float normalization = (Float) 1 / std::sqrt(slope.x*slope.x
+				+ slope.y*slope.y + (Float) 1.0);
+
+		return Normal(
+			-slope.x * normalization,
+			-slope.y * normalization,
+			normalization
 		);
 	}
 
+	/// Implements the probability density of the function \ref sampleVisible()
+	Float pdfVisible(const Vector &wi, const Vector &m) const {
+		return smithG1(wi, m) * absDot(wi, m) * eval(m) / std::abs(Frame::cosTheta(wi));
+	}
+
 	/**
 	 * \brief Smith's shadowing-masking function G1 for each
 	 * of the supported microfacet distributions
 	 *
-	 * \param v An arbitrary direction
-	 * \param m The microsurface normal
-	 * \param alpha The surface roughness
+	 * \param v
+	 *     An arbitrary direction
+	 * \param m
+	 *     The microfacet normal
 	 */
-	Float smithG1(const Vector &v, const Vector &m, Float alpha) const {
-		const Float tanTheta = std::abs(Frame::tanTheta(v));
-
-		/* perpendicular incidence -- no shadowing/masking */
-		if (tanTheta == 0.0f)
-			return 1.0f;
-
-		/* Can't see the back side from the front and vice versa */
+	Float smithG1(const Vector &v, const Vector &m) const {
+		/* Ensure consistent orientation (can't see the back
+		   of the microfacet from the front and vice versa) */
 		if (dot(v, m) * Frame::cosTheta(v) <= 0)
 			return 0.0f;
 
+		/* Perpendicular incidence -- no shadowing/masking */
+		Float tanTheta = std::abs(Frame::tanTheta(v));
+		if (tanTheta == 0.0f)
+			return 1.0f;
+
+		Float alpha = projectRoughness(v);
 		switch (m_type) {
 			case EPhong:
 			case EBeckmann: {
-					Float a;
-					/* Approximation recommended by Bruce Walter: Use
-					   the Beckmann shadowing-masking function with
-					   specially chosen roughness value */
-					if (m_type != EBeckmann)
-						a = std::sqrt(0.5f * alpha + 1) / tanTheta;
-					else
-						a = 1.0f / (alpha * tanTheta);
-
+					Float a = 1.0f / (alpha * tanTheta);
 					if (a >= 1.6f)
 						return 1.0f;
 
 					/* Use a fast and accurate (<0.35% rel. error) rational
 					   approximation to the shadowing-masking function */
-					const Float aSqr = a * a;
+					Float aSqr = a*a;
 					return (3.535f * a + 2.181f * aSqr)
 						 / (1.0f + 2.276f * a + 2.577f * aSqr);
 				}
 				break;
 
 			case EGGX: {
-					const Float root = alpha * tanTheta;
+					Float root = alpha * tanTheta;
 					return 2.0f / (1.0f + std::sqrt(1.0f + root*root));
 				}
 				break;
 
 			default:
-				SLog(EError, "Invalid distribution function!");
-				return 0.0f;
+				SLog(EError, "Invalid distribution type!");
+				return -1.0f;
 		}
 	}
 
 	/**
-	 * \brief Shadow-masking function for each of the supported
-	 * microfacet distributions
-	 *
-	 * \param wi The incident direction
-	 * \param wo The exitant direction
-	 * \param m The microsurface normal
-	 * \param alpha The surface roughness
+	 * \brief Separable shadow-masking function based on Smith's
+	 * one-dimensional masking model
 	 */
-	inline Float G(const Vector &wi, const Vector &wo, const Vector &m, Float alpha) const {
-		return G(wi, wo, m, alpha, alpha);
+	inline Float G(const Vector &wi, const Vector &wo, const Vector &m) const {
+		return smithG1(wi, m) * smithG1(wo, m);
 	}
 
-	/**
-	 * \brief Shadow-masking function for each of the supported
-	 * microfacet distributions
-	 *
-	 * \param wi The incident direction
-	 * \param wo The exitant direction
-	 * \param m The microsurface normal
-	 * \param alphaU  The surface roughness in the tangent direction
-	 * \param alphaV  The surface roughness in the bitangent direction
-	 */
-	Float G(const Vector &wi, const Vector &wo, const Vector &m, Float alphaU, Float alphaV) const {
-		if (m_type != EAshikhminShirley) {
-			return smithG1(wi, m, alphaU)
-				 * smithG1(wo, m, alphaU);
-		} else {
-			/* Can't see the back side from the front and vice versa */
-			if (dot(wi, m) * Frame::cosTheta(wi) <= 0 ||
-				dot(wo, m) * Frame::cosTheta(wo) <= 0)
-				return 0.0f;
-
-			/* Infinite groove shadowing/masking */
-			const Float nDotM  = Frame::cosTheta(m),
-						nDotWo = Frame::cosTheta(wo),
-						nDotWi = Frame::cosTheta(wi),
-						woDotM = dot(wo, m),
-						wiDotM = dot(wi, m);
-
-			return std::min((Float) 1,
-				std::min(std::abs(2 * nDotM * nDotWo / woDotM),
-						 std::abs(2 * nDotM * nDotWi / wiDotM)));
-		}
-	}
-
-	std::string toString() const {
-		switch (m_type) {
+	/// Return a string representation of the name of a distribution
+	inline static std::string distributionName(EType type) {
+		switch (type) {
 			case EBeckmann: return "beckmann"; break;
-			case EPhong: return "phong"; break;
 			case EGGX: return "ggx"; break;
-			case EAshikhminShirley: return "as"; break;
-			default:
-				SLog(EError, "Invalid distribution function");
-				return "";
+			case EPhong: return "phong"; break;
+			default: return "invalid"; break;
 		}
 	}
+
+	/// Return a string representation of the contents of this instance
+	std::string toString() const {
+		return formatString("MicrofacetDistribution[type=\"%s\", alphaU=%f, alphaV=%f]",
+			distributionName(m_type).c_str(), m_alphaU, m_alphaV);
+	}
+protected:
+	/// Compute the effective roughness projected on direction \c v
+	inline Float projectRoughness(const Vector &v) const {
+		Float invSinTheta2 = 1 / Frame::sinTheta2(v);
+
+		if (isIsotropic() || invSinTheta2 <= 0)
+			return m_alphaU;
+
+		Float cosPhi2 = v.x * v.x * invSinTheta2;
+		Float sinPhi2 = v.y * v.y * invSinTheta2;
+
+		return std::sqrt(cosPhi2 * m_alphaU * m_alphaU + sinPhi2 * m_alphaV * m_alphaV);
+	}
+
+	/// Compute the interpolated roughness for the Phong model
+	inline Float interpolatePhongExponent(const Vector &v) const {
+		Float invSinTheta2 = 1 / Frame::sinTheta2(v);
+
+		if (isIsotropic() || invSinTheta2 <= 0)
+			return m_exponentU;
+
+		Float cosPhi2 = v.x * v.x * invSinTheta2;
+		Float sinPhi2 = v.y * v.y * invSinTheta2;
+
+		return m_exponentU * cosPhi2 + m_exponentV * sinPhi2;
+	}
+
+	/**
+	 * \brief Visible normal sampling code for the alpha=1 case
+	 *
+	 * Source: supplemental material of "Importance Sampling
+	 * Microfacet-Based BSDFs using the Distribution of Visible Normals"
+	 */
+	Vector2 sampleVisible11(Float thetaI, Point2 sample) const {
+		const Float SQRT_PI_INV = 1 / std::sqrt(M_PI);
+		Vector2 slope;
+
+		switch (m_type) {
+			case EBeckmann: {
+					/* Special case (normal incidence) */
+					if (thetaI < 1e-4f) {
+						Float sinPhi, cosPhi;
+						Float r = std::sqrt(-math::fastlog(1.0f-sample.x));
+						math::sincos(2 * M_PI * sample.y, &sinPhi, &cosPhi);
+						return Vector2(r * cosPhi, r * sinPhi);
+					}
+
+					/* The original inversion routine from the paper contained
+					   discontinuities, which causes issues for QMC integration
+					   and techniques like Kelemen-style MLT. The following code
+					   performs a numerical inversion with better behavior */
+					Float tanThetaI = std::tan(thetaI);
+					Float cotThetaI = 1 / tanThetaI;
+
+					/* Search interval -- everything is parameterized
+					   in the erf() domain */
+					Float a = -1, c = math::erf(cotThetaI);
+					Float sample_x = std::max(sample.x, (Float) 1e-6f);
+
+					/* Start with a good initial guess */
+					//Float b = (1-sample_x) * a + sample_x * c;
+
+					/* We can do better (inverse of an approximation computed in Mathematica) */
+					Float fit = 1 + thetaI*(-0.876f + thetaI * (0.4265f - 0.0594f*thetaI));
+					Float b = c - (1+c) * std::pow(1-sample_x, fit);
+
+					/* Normalization factor for the CDF */
+					Float normalization = 1 / (1 + c + SQRT_PI_INV*
+						tanThetaI*std::exp(-cotThetaI*cotThetaI));
+
+					int it = 0;
+					while (++it < 10) {
+						/* Bisection criterion -- the oddly-looking
+						   boolean expression are intentional to check
+						   for NaNs at little additional cost */
+						if (!(b >= a && b <= c))
+							b = 0.5f * (a + c);
+
+						/* Evaluate the CDF and its derivative
+						   (i.e. the density function) */
+						Float invErf = math::erfinv(b);
+						Float value = normalization*(1 + b + SQRT_PI_INV*
+							tanThetaI*std::exp(-invErf*invErf)) - sample_x;
+						Float derivative = normalization * (1
+							- invErf*tanThetaI);
+
+						if (std::abs(value) < 1e-5f)
+							break;
+
+						/* Update bisection intervals */
+						if (value > 0)
+							c = b;
+						else
+							a = b;
+
+						b -= value / derivative;
+					}
+
+					/* Now convert back into a slope value */
+					slope.x = math::erfinv(b);
+
+					/* Simulate Y component */
+					slope.y = math::erfinv(2.0f*std::max(sample.y, (Float) 1e-6f) - 1.0f);
+				};
+				break;
+
+			case EGGX: {
+					/* Special case (normal incidence) */
+					if (thetaI < 1e-4f) {
+						Float sinPhi, cosPhi;
+						Float r = math::safe_sqrt(sample.x / (1 - sample.x));
+						math::sincos(2 * M_PI * sample.y, &sinPhi, &cosPhi);
+						return Vector2(r * cosPhi, r * sinPhi);
+					}
+
+					/* Precomputations */
+					Float tanThetaI = std::tan(thetaI);
+					Float a = 1 / tanThetaI;
+					Float G1 = 2.0f / (1.0f + math::safe_sqrt(1.0f + 1.0f / (a*a)));
+
+					/* Simulate X component */
+					Float A = 2.0f * sample.x / G1 - 1.0f;
+					if (std::abs(A) == 1)
+						A -= math::signum(A)*Epsilon;
+					Float tmp = 1.0f / (A*A - 1.0f);
+					Float B = tanThetaI;
+					Float D = math::safe_sqrt(B*B*tmp*tmp - (A*A - B*B) * tmp);
+					Float slope_x_1 = B * tmp - D;
+					Float slope_x_2 = B * tmp + D;
+					slope.x = (A < 0.0f || slope_x_2 > 1.0f / tanThetaI) ? slope_x_1 : slope_x_2;
+
+					/* Simulate Y component */
+					Float S;
+					if (sample.y > 0.5f) {
+						S = 1.0f;
+						sample.y = 2.0f * (sample.y - 0.5f);
+					} else {
+						S = -1.0f;
+						sample.y = 2.0f * (0.5f - sample.y);
+					}
+
+					/* Improved fit */
+					Float z =
+						(sample.y * (sample.y * (sample.y * (-(Float) 0.365728915865723) + (Float) 0.790235037209296) -
+							(Float) 0.424965825137544) + (Float) 0.000152998850436920) /
+						(sample.y * (sample.y * (sample.y * (sample.y * (Float) 0.169507819808272 - (Float) 0.397203533833404) -
+							(Float) 0.232500544458471) + (Float) 1) - (Float) 0.539825872510702);
+
+					slope.y = S * z * std::sqrt(1.0f + slope.x*slope.x);
+				};
+				break;
+
+			default:
+				SLog(EError, "Invalid distribution type!");
+				return Vector2(-1);
+		};
+		return slope;
+	}
+
+
+	/// Helper routine: convert from Beckmann-style roughness values to Phong exponents (Walter et al.)
+	void computePhongExponent() {
+		m_exponentU = std::max(2.0f / (m_alphaU * m_alphaU) - 2.0f, (Float) 0.0f);
+		m_exponentV = std::max(2.0f / (m_alphaV * m_alphaV) - 2.0f, (Float) 0.0f);
+	}
+
+	/// Helper routine: sample the azimuthal part of the first quadrant of the A&S distribution
+	void sampleFirstQuadrant(Float u1, Float &phi, Float &exponent) const {
+		Float cosPhi, sinPhi;
+		phi = std::atan(
+				std::sqrt((m_exponentU + 2.0f) / (m_exponentV + 2.0f)) *
+				std::tan(M_PI * u1 * 0.5f));
+		math::sincos(phi, &sinPhi, &cosPhi);
+		/* Return the interpolated roughness */
+		exponent = m_exponentU * cosPhi * cosPhi + m_exponentV * sinPhi * sinPhi;
+	}
 protected:
 	EType m_type;
+	Float m_alphaU, m_alphaV;
+	bool m_sampleVisible;
+	Float m_exponentU, m_exponentV;
 };
 
 MTS_NAMESPACE_END
diff --git a/src/bsdfs/roughcoating.cpp b/src/bsdfs/roughcoating.cpp
index 3baa690e..a4817cdb 100644
--- a/src/bsdfs/roughcoating.cpp
+++ b/src/bsdfs/roughcoating.cpp
@@ -31,17 +31,20 @@ MTS_NAMESPACE_BEGIN
  *     \parameter{distribution}{\String}{
  *          Specifies the type of microfacet normal distribution
  *          used to model the surface roughness.
+ *          \vspace{-1mm}
  *       \begin{enumerate}[(i)]
  *           \item \code{beckmann}: Physically-based distribution derived from
- *               Gaussian random surfaces. This is the default.
- *           \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.
- *           \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.
+ *               Gaussian random surfaces. This is the default.\vspace{-1.5mm}
+ *           \item \code{ggx}: The GGX \cite{Walter07Microfacet} distribution (also known as
+ *               Trowbridge-Reitz \cite{Trowbridge19975Average} distribution)
+ *               was designed to better approximate the long tails observed in measurements
+ *               of ground surfaces, which are not modeled by the Beckmann distribution.
+ *           \vspace{-1.5mm}
+ *           \item \code{phong}: Classical Phong distribution.
+ *              In most cases, the \code{ggx} and \code{beckmann} distributions
+ *              should be preferred, since they provide better importance sampling
+ *              and accurate shadowing/masking computations.
+ *              \vspace{-4mm}
  *       \end{enumerate}
  *     }
  *     \parameter{alpha}{\Float\Or\Texture}{
@@ -50,6 +53,11 @@ MTS_NAMESPACE_BEGIN
  *         \emph{root mean square} (RMS) slope of the microfacets.
  *         \default{0.1}.
  *     }
+ *     \parameter{sampleVisible}{\Boolean}{
+ *	       Enables an improved importance sampling technique. Refer to
+ *	       pages \pageref{plg:roughconductor} and \pageref{sec:visiblenormal-sampling}
+ *	       for details. \default{\code{true}}
+ *     }
  *     \parameter{intIOR}{\Float\Or\String}{Interior index of refraction specified
  *      numerically or using a known material name. \default{\texttt{bk7} / 1.5046}}
  *     \parameter{extIOR}{\Float\Or\String}{Exterior index of refraction specified
@@ -90,9 +98,10 @@ MTS_NAMESPACE_BEGIN
  * loss for materials that reflect much of their energy near or below the critical
  * angle (i.e. diffuse or very rough materials).
  *
- * The implementation here is influenced by the paper
+ * The implementation here is motivated by the paper
  * ``Arbitrarily Layered Micro-Facet Surfaces'' by Weidlich and
- * Wilkie \cite{Weidlich2007Arbitrarily}.
+ * Wilkie \cite{Weidlich2007Arbitrarily}, though the implementation
+ * works differently.
  */
 class RoughCoating : public BSDF {
 public:
@@ -127,25 +136,23 @@ public:
 		m_specularReflectance = new ConstantSpectrumTexture(
 			props.getSpectrum("specularReflectance", Spectrum(1.0f)));
 
-		m_distribution = MicrofacetDistribution(
-			props.getString("distribution", "beckmann")
-		);
+		MicrofacetDistribution distr(props);
+		m_type = distr.getType();
+		m_sampleVisible = distr.getSampleVisible();
 
-		if (m_distribution.isAnisotropic())
-			Log(EError, "The 'roughcoating' plugin currently does not support "
+		if (distr.isAnisotropic())
+			Log(EError, "The 'roughplastic' plugin currently does not support "
 				"anisotropic microfacet distributions!");
 
-		m_alpha = new ConstantFloatTexture(
-			props.getFloat("alpha", 0.1f));
+		m_alpha = new ConstantFloatTexture(distr.getAlpha());
 
 		m_specularSamplingWeight = 0.0f;
 	}
 
 	RoughCoating(Stream *stream, InstanceManager *manager)
 	 : BSDF(stream, manager) {
-		m_distribution = MicrofacetDistribution(
-			(MicrofacetDistribution::EType) stream->readUInt()
-		);
+		m_type = (MicrofacetDistribution::EType) stream->readUInt();
+		m_sampleVisible = stream->readBool();
 		m_nested = static_cast<BSDF *>(manager->getInstance(stream));
 		m_sigmaA = static_cast<Texture *>(manager->getInstance(stream));
 		m_specularReflectance = static_cast<Texture *>(manager->getInstance(stream));
@@ -160,7 +167,8 @@ public:
 	void serialize(Stream *stream, InstanceManager *manager) const {
 		BSDF::serialize(stream, manager);
 
-		stream->writeUInt((uint32_t) m_distribution.getType());
+		stream->writeUInt((uint32_t) m_type);
+		stream->writeBool(m_sampleVisible);
 		manager->serialize(stream, m_nested.get());
 		manager->serialize(stream, m_sigmaA.get());
 		manager->serialize(stream, m_specularReflectance.get());
@@ -200,8 +208,7 @@ public:
 		if (!m_roughTransmittance.get()) {
 			/* Load precomputed data used to compute the rough
 			   transmittance through the dielectric interface */
-			m_roughTransmittance = new RoughTransmittance(
-				m_distribution.getType());
+			m_roughTransmittance = new RoughTransmittance(m_type);
 
 			m_roughTransmittance->checkEta(m_eta);
 			m_roughTransmittance->checkAlpha(m_alpha->getMinimum().average());
@@ -254,9 +261,13 @@ public:
 			&& (bRec.component == -1 || bRec.component == (int) m_components.size()-1)
 			&& measure == ESolidAngle;
 
-		/* Evaluate the roughness texture */
-		Float alpha = m_alpha->eval(bRec.its).average();
-		Float alphaT = m_distribution.transformRoughness(alpha);
+		/* Construct the microfacet distribution matching the
+		   roughness values at the current surface position. */
+		MicrofacetDistribution distr(
+			m_type,
+			m_alpha->eval(bRec.its).average(),
+			m_sampleVisible
+		);
 
 		Spectrum result(0.0f);
 		if (hasSpecular && Frame::cosTheta(bRec.wo) * Frame::cosTheta(bRec.wi) > 0) {
@@ -264,14 +275,14 @@ public:
 			const Vector H = normalize(bRec.wo+bRec.wi)
 				* math::signum(Frame::cosTheta(bRec.wo));
 
-			/* Evaluate the microsurface normal distribution */
-			const Float D = m_distribution.eval(H, alphaT);
+			/* Evaluate the microfacet normal distribution */
+			const Float D = distr.eval(H);
 
 			/* Fresnel term */
 			const Float F = fresnelDielectricExt(absDot(bRec.wi, H), m_eta);
 
 			/* Smith's shadow-masking function */
-			const Float G = m_distribution.G(bRec.wi, bRec.wo, H, alphaT);
+			const Float G = distr.G(bRec.wi, bRec.wo, H);
 
 			/* Calculate the specular reflection component */
 			Float value = F * D * G /
@@ -286,8 +297,8 @@ public:
 			bRecInt.wo = refractTo(EInterior, bRec.wo);
 
 			Spectrum nestedResult = m_nested->eval(bRecInt, measure) *
-				m_roughTransmittance->eval(std::abs(Frame::cosTheta(bRec.wi)), alpha) *
-				m_roughTransmittance->eval(std::abs(Frame::cosTheta(bRec.wo)), alpha);
+				m_roughTransmittance->eval(std::abs(Frame::cosTheta(bRec.wi)), distr.getAlpha()) *
+				m_roughTransmittance->eval(std::abs(Frame::cosTheta(bRec.wo)), distr.getAlpha());
 
 			Spectrum sigmaA = m_sigmaA->eval(bRec.its) * m_thickness;
 			if (!sigmaA.isZero())
@@ -319,15 +330,19 @@ public:
 		const Vector H = normalize(bRec.wo+bRec.wi)
 				* math::signum(Frame::cosTheta(bRec.wo));
 
-		/* Evaluate the roughness texture */
-		Float alpha = m_alpha->eval(bRec.its).average();
-		Float alphaT = m_distribution.transformRoughness(alpha);
+		/* Construct the microfacet distribution matching the
+		   roughness values at the current surface position. */
+		MicrofacetDistribution distr(
+			m_type,
+			m_alpha->eval(bRec.its).average(),
+			m_sampleVisible
+		);
 
 		Float probNested, probSpecular;
 		if (hasSpecular && hasNested) {
 			/* Find the probability of sampling the specular component */
 			probSpecular = 1-m_roughTransmittance->eval(
-				std::abs(Frame::cosTheta(bRec.wi)), alpha);
+				std::abs(Frame::cosTheta(bRec.wi)), distr.getAlpha());
 
 			/* Reallocate samples */
 			probSpecular = (probSpecular*m_specularSamplingWeight) /
@@ -344,8 +359,8 @@ public:
 			/* Jacobian of the half-direction mapping */
 			const Float dwh_dwo = 1.0f / (4.0f * absDot(bRec.wo, H));
 
-			/* Evaluate the microsurface normal distribution */
-			const Float prob = m_distribution.pdf(H, alphaT);
+			/* Evaluate the microfacet model sampling density function */
+			const Float prob = distr.pdf(bRec.wi, H);
 
 			result = prob * dwh_dwo * probSpecular;
 		}
@@ -377,14 +392,19 @@ public:
 		bool choseSpecular = hasSpecular;
 		Point2 sample(_sample);
 
-		/* Evaluate the roughness texture */
-		Float alpha = m_alpha->eval(bRec.its).average();
-		Float alphaT = m_distribution.transformRoughness(alpha);
+		/* Construct the microfacet distribution matching the
+		   roughness values at the current surface position. */
+		MicrofacetDistribution distr(
+			m_type,
+			m_alpha->eval(bRec.its).average(),
+			m_sampleVisible
+		);
 
 		Float probSpecular;
 		if (hasSpecular && hasNested) {
 			/* Find the probability of sampling the diffuse component */
-			probSpecular = 1 - m_roughTransmittance->eval(std::abs(Frame::cosTheta(bRec.wi)), alpha);
+			probSpecular = 1 - m_roughTransmittance->eval(
+				std::abs(Frame::cosTheta(bRec.wi)), distr.getAlpha());
 
 			/* Reallocate samples */
 			probSpecular = (probSpecular*m_specularSamplingWeight) /
@@ -400,8 +420,8 @@ public:
 		}
 
 		if (choseSpecular) {
-			/* Perfect specular reflection based on the microsurface normal */
-			Normal m = m_distribution.sample(sample, alphaT);
+			/* Perfect specular reflection based on the microfacet normal */
+			Normal m = distr.sample(bRec.wi, sample);
 			bRec.wo = reflect(bRec.wi, m);
 			bRec.sampledComponent = (int) m_components.size() - 1;
 			bRec.sampledType = EGlossyReflection;
@@ -464,7 +484,8 @@ public:
 		std::ostringstream oss;
 		oss << "RoughCoating[" << endl
 			<< "  id = \"" << getID() << "\"," << endl
-			<< "  distribution = " << m_distribution.toString() << "," << endl
+			<< "  distribution = " << MicrofacetDistribution::distributionName(m_type) << "," << endl
+			<< "  sampleVisible = " << m_sampleVisible << "," << endl
 			<< "  alpha = " << indent(m_alpha->toString()) << "," << endl
 			<< "  sigmaA = " << indent(m_sigmaA->toString()) << "," << endl
 			<< "  specularReflectance = " << indent(m_specularReflectance->toString()) << "," << endl
@@ -480,7 +501,7 @@ public:
 
 	MTS_DECLARE_CLASS()
 private:
-	MicrofacetDistribution m_distribution;
+	MicrofacetDistribution::EType m_type;
 	ref<RoughTransmittance> m_roughTransmittance;
 	ref<Texture> m_sigmaA;
 	ref<Texture> m_alpha;
@@ -489,6 +510,7 @@ private:
 	Float m_eta, m_invEta;
 	Float m_specularSamplingWeight;
 	Float m_thickness;
+	bool m_sampleVisible;
 };
 
 /**
diff --git a/src/bsdfs/roughconductor.cpp b/src/bsdfs/roughconductor.cpp
index 76c73aab..2146c60b 100644
--- a/src/bsdfs/roughconductor.cpp
+++ b/src/bsdfs/roughconductor.cpp
@@ -31,32 +31,30 @@ MTS_NAMESPACE_BEGIN
  *     \parameter{distribution}{\String}{
  *          Specifies the type of microfacet normal distribution
  *          used to model the surface roughness.
+ *          \vspace{-1mm}
  *       \begin{enumerate}[(i)]
  *           \item \code{beckmann}: Physically-based distribution derived from
- *               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.\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.\vspace{-1mm}
- *           \item \code{as}: Anisotropic Phong-style microfacet distribution proposed by
- *              Ashikhmin and Shirley \cite{Ashikhmin2005Anisotropic}.\vspace{-3mm}
+ *               Gaussian random surfaces. This is the default.\vspace{-1.5mm}
+ *           \item \code{ggx}: The GGX \cite{Walter07Microfacet} distribution (also known as
+ *               Trowbridge-Reitz \cite{Trowbridge19975Average} distribution)
+ *               was designed to better approximate the long tails observed in measurements
+ *               of ground surfaces, which are not modeled by the Beckmann distribution.
+ *           \vspace{-1.5mm}
+ *           \item \code{phong}: Anisotropic Phong distribution by
+ *              Ashikhmin and Shirley \cite{Ashikhmin2005Anisotropic}.
+ *              In most cases, the \code{ggx} and \code{beckmann} distributions
+ *              should be preferred, since they provide better importance sampling
+ *              and accurate shadowing/masking computations.
+ *              \vspace{-4mm}
  *       \end{enumerate}
  *     }
- *     \parameter{alpha}{\Float\Or\Texture}{
- *         Specifies the roughness of the unresolved surface micro-geometry.
- *         When the Beckmann distribution is used, this parameter is equal to the
- *         \emph{root mean square} (RMS) slope of the microfacets. This
- *         parameter is only valid when \texttt{distribution=beckmann/phong/ggx}.
- *         \default{0.1}.
- *     }
- *     \parameter{alphaU, alphaV}{\Float\Or\Texture}{
- *         Specifies the anisotropic roughness values along the tangent and
- *         bitangent directions. These parameter are only valid when
- *         \texttt{distribution=as}. \default{0.1}.
+ *     \parameter{alpha, alphaU, alphaV}{\Float\Or\Texture}{
+ *         Specifies the roughness of the unresolved surface micro-geometry
+ *         along the tangent and bitangent directions. When the Beckmann
+ *         distribution is used, this parameter is equal to the
+ *         \emph{root mean square} (RMS) slope of the microfacets.
+ *         \code{alpha} is a convenience parameter to initialize both
+ *         \code{alphaU} and \code{alphaV} to the same value. \default{0.1}.
  *     }
  *     \parameter{material}{\String}{Name of a material preset, see
  *           \tblref{conductor-iors}.\!\default{\texttt{Cu} / copper}}
@@ -66,30 +64,36 @@ MTS_NAMESPACE_BEGIN
  *           Real-valued index of refraction of the surrounding dielectric,
  *           or a material name of a dielectric \default{\code{air}}
  *     }
+ *     \parameter{sampleVisible}{\Boolean}{
+ *         Enables a sampling technique proposed by Heitz and D'Eon~\cite{Heitz1014Importance},
+ *         which focuses computation on the visible parts of the microfacet normal
+ *         distribution, considerably reducing variance in some cases.
+ *         \default{\code{true}, i.e. use visible normal sampling}
+ *     }
  *     \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}}
  * }
- * \vspace{4mm}
+ * \vspace{3mm}
  * This plugin implements a realistic microfacet scattering model for rendering
  * rough conducting materials, such as metals. It can be interpreted as a fancy
  * version of the Cook-Torrance model and should be preferred over
- * heuristic models like \pluginref{phong} and \pluginref{ward} when possible.
+ * heuristic models like \pluginref{phong} and \pluginref{ward} if possible.
  * \renderings{
  *     \rendering{Rough copper (Beckmann, $\alpha=0.1$)}
  *     	   {bsdf_roughconductor_copper.jpg}
- *     \rendering{Vertically brushed aluminium (Ashikhmin-Shirley,
+ *     \rendering{Vertically brushed aluminium (Anisotropic Phong,
  *         $\alpha_u=0.05,\ \alpha_v=0.3$), see
  *         \lstref{roughconductor-aluminium}}
  *         {bsdf_roughconductor_anisotropic_aluminium.jpg}
  * }
  *
- * Microfacet theory describes rough
- * surfaces as an arrangement of unresolved and ideally specular facets, whose
- * normal directions are given by a specially chosen \emph{microfacet distribution}.
- * By accounting for shadowing and masking effects between these facets, it is
- * possible to reproduce the important off-specular reflections peaks observed
- * in real-world measurements of such materials.
+ * Microfacet theory describes rough surfaces as an arrangement of unresolved
+ * and ideally specular facets, whose normal directions are given by a
+ * specially chosen \emph{microfacet distribution}. By accounting for shadowing
+ * and masking effects between these facets, it is possible to reproduce the
+ * important off-specular reflections peaks observed in real-world measurements
+ * of such materials.
  *
  * This plugin is essentially the ``roughened'' equivalent of the (smooth) plugin
  * \pluginref{conductor}. For very low values of $\alpha$, the two will
@@ -98,7 +102,7 @@ MTS_NAMESPACE_BEGIN
  *
  * The implementation is based on the paper ``Microfacet Models
  * for Refraction through Rough Surfaces'' by Walter et al.
- * \cite{Walter07Microfacet}. It supports several different types of microfacet
+ * \cite{Walter07Microfacet}. It supports three different types of microfacet
  * distributions and has a texturable roughness parameter.
  * To facilitate the tedious task of specifying spectrally-varying index of
  * refraction information, this plugin can access a set of measured materials
@@ -109,41 +113,48 @@ MTS_NAMESPACE_BEGIN
  * 100% reflecting mirror.
  *
  * When no parameters are given, the plugin activates the default settings,
- * which describe copper with a light amount of roughness modeled using a
+ * which describe copper with a medium amount of roughness modeled using a
  * Beckmann distribution.
  *
- * To get an intuition about the effect of the surface roughness
- * parameter $\alpha$, consider the following approximate classification:
- * a value of $\alpha=0.001-0.01$ corresponds to a material
- * with slight imperfections on an
- * otherwise smooth surface finish, $\alpha=0.1$ is relatively rough,
+ * To get an intuition about the effect of the surface roughness parameter
+ * $\alpha$, consider the following approximate classification: a value of
+ * $\alpha=0.001-0.01$ corresponds to a material with slight imperfections
+ * on an otherwise smooth surface finish, $\alpha=0.1$ is relatively rough,
  * and $\alpha=0.3-0.7$ is \emph{extremely} rough (e.g. an etched or ground
  * finish). Values significantly above that are probably not too realistic.
  * \vspace{4mm}
  * \begin{xml}[caption={A material definition for brushed aluminium}, label=lst:roughconductor-aluminium]
  * <bsdf type="roughconductor">
  *     <string name="material" value="Al"/>
- *     <string name="distribution" value="as"/>
+ *     <string name="distribution" value="phong"/>
  *     <float name="alphaU" value="0.05"/>
  *     <float name="alphaV" value="0.3"/>
  * </bsdf>
  * \end{xml}
  *
  * \subsubsection*{Technical details}
- * When rendering with the Ashikhmin-Shirley or Phong microfacet
- * distributions, a conversion is used to turn the specified
- * $\alpha$ roughness value into the exponents of these distributions.
- * This is done in a way, such that the different
- * distributions all produce a similar appearance for the same value of
- * $\alpha$.
+ * All microfacet distributions allow the specification of two distinct
+ * roughness values along the tangent and bitangent directions. This can be
+ * used to provide a material with a ``brushed'' appearance. The alignment
+ * of the anisotropy will follow the UV parameterization of the underlying
+ * mesh. This means that such an anisotropic material cannot be applied to
+ * triangle meshes that are missing texture coordinates.
  *
- * The Ashikhmin-Shirley microfacet distribution allows the specification
- * of two distinct roughness values along the tangent and bitangent
- * directions. This can be used to provide a material with a ``brushed''
- * appearance. The alignment of the anisotropy will follow the UV
- * parameterization of the underlying mesh in this case. This also means that
- * such an anisotropic material cannot be applied to triangle meshes that
- * are missing texture coordinates.
+ * \label{sec:visiblenormal-sampling}
+ * Since Mitsuba 0.5.1, this plugin uses a new importance sampling technique
+ * contributed by Eric Heitz and Eugene D'Eon, which restricts the sampling
+ * domain to the set of visible (unmasked) microfacet normals. The previous
+ * approach of sampling all normals is still available and can be enabled
+ * by setting \code{sampleVisible} to \code{false}.
+ * Note that this new method is only available for the \code{beckmann} and
+ * \code{ggx} microfacet distributions. When the \code{phong} distribution
+ * is selected, the parameter has no effect.
+ *
+ * When rendering with the Phong microfacet distribution, a conversion is
+ * used to turn the specified Beckmann-equivalent $\alpha$ roughness value
+ * into the exponent parameter of this distribution. This is done in a way,
+ * such that the same value $\alpha$ will produce a similar appearance across
+ * different microfacet distributions.
  *
  * When using this plugin, you should ideally compile Mitsuba with support for
  * spectral rendering to get the most accurate results. While it also works
@@ -178,26 +189,21 @@ public:
 		m_eta = props.getSpectrum("eta", intEta) / extEta;
 		m_k   = props.getSpectrum("k", intK) / extEta;
 
-		m_distribution = MicrofacetDistribution(
-			props.getString("distribution", "beckmann")
-		);
+		MicrofacetDistribution distr(props);
+		m_type = distr.getType();
+		m_sampleVisible = distr.getSampleVisible();
 
-		Float alpha = props.getFloat("alpha", 0.1f),
-			  alphaU = props.getFloat("alphaU", alpha),
-			  alphaV = props.getFloat("alphaV", alpha);
-
-		m_alphaU = new ConstantFloatTexture(alphaU);
-		if (alphaU == alphaV)
+		m_alphaU = new ConstantFloatTexture(distr.getAlphaU());
+		if (distr.getAlphaU() == distr.getAlphaV())
 			m_alphaV = m_alphaU;
 		else
-			m_alphaV = new ConstantFloatTexture(alphaV);
+			m_alphaV = new ConstantFloatTexture(distr.getAlphaV());
 	}
 
 	RoughConductor(Stream *stream, InstanceManager *manager)
 	 : BSDF(stream, manager) {
-		m_distribution = MicrofacetDistribution(
-			(MicrofacetDistribution::EType) stream->readUInt()
-		);
+		m_type = (MicrofacetDistribution::EType) stream->readUInt();
+		m_sampleVisible = stream->readBool();
 		m_alphaU = static_cast<Texture *>(manager->getInstance(stream));
 		m_alphaV = static_cast<Texture *>(manager->getInstance(stream));
 		m_specularReflectance = static_cast<Texture *>(manager->getInstance(stream));
@@ -207,17 +213,23 @@ public:
 		configure();
 	}
 
+	void serialize(Stream *stream, InstanceManager *manager) const {
+		BSDF::serialize(stream, manager);
+
+		stream->writeUInt((uint32_t) m_type);
+		stream->writeBool(m_sampleVisible);
+		manager->serialize(stream, m_alphaU.get());
+		manager->serialize(stream, m_alphaV.get());
+		manager->serialize(stream, m_specularReflectance.get());
+		m_eta.serialize(stream);
+		m_k.serialize(stream);
+	}
+
 	void configure() {
 		unsigned int extraFlags = 0;
-		if (m_alphaU != m_alphaV) {
+		if (m_alphaU != m_alphaV)
 			extraFlags |= EAnisotropic;
-			if (m_distribution.getType() !=
-				MicrofacetDistribution::EAshikhminShirley)
-				Log(EError, "Different roughness values along the tangent and "
-						"bitangent directions are only supported when using the "
-						"anisotropic Ashikhmin-Shirley microfacet distribution "
-						"(named \"as\")");
-		}
+
 		if (!m_alphaU->isConstant() || !m_alphaV->isConstant() ||
 			!m_specularReflectance->isConstant())
 			extraFlags |= ESpatiallyVarying;
@@ -254,27 +266,31 @@ public:
 		/* Calculate the reflection half-vector */
 		Vector H = normalize(bRec.wo+bRec.wi);
 
-		/* Evaluate the roughness */
-		Float alphaU = m_distribution.transformRoughness(
-					m_alphaU->eval(bRec.its).average()),
-			  alphaV = m_distribution.transformRoughness(
-					m_alphaV->eval(bRec.its).average());
+		/* Construct the microfacet distribution matching the
+		   roughness values at the current surface position. */
+		MicrofacetDistribution distr(
+			m_type,
+			m_alphaU->eval(bRec.its).average(),
+			m_alphaV->eval(bRec.its).average(),
+			m_sampleVisible
+		);
 
-		/* Evaluate the microsurface normal distribution */
-		const Float D = m_distribution.eval(H, alphaU, alphaV);
+		/* Evaluate the microfacet normal distribution */
+		const Float D = distr.eval(H);
 		if (D == 0)
 			return Spectrum(0.0f);
 
 		/* Fresnel factor */
-		const Spectrum F = fresnelConductorExact(dot(bRec.wi, H), m_eta, m_k);
+		const Spectrum F = fresnelConductorExact(dot(bRec.wi, H), m_eta, m_k) *
+			m_specularReflectance->eval(bRec.its);
 
 		/* Smith's shadow-masking function */
-		const Float G = m_distribution.G(bRec.wi, bRec.wo, H, alphaU, alphaV);
+		const Float G = distr.G(bRec.wi, bRec.wo, H);
 
 		/* Calculate the total amount of reflection */
-		Float value = D * G / (4.0f * Frame::cosTheta(bRec.wi));
+		Float model = D * G / (4.0f * Frame::cosTheta(bRec.wi));
 
-		return m_specularReflectance->eval(bRec.its) * F * value;
+		return F * model;
 	}
 
 	Float pdf(const BSDFSamplingRecord &bRec, EMeasure measure) const {
@@ -288,14 +304,20 @@ public:
 		/* Calculate the reflection half-vector */
 		Vector H = normalize(bRec.wo+bRec.wi);
 
-		/* Evaluate the roughness */
-		Float alphaU = m_distribution.transformRoughness(
-					m_alphaU->eval(bRec.its).average()),
-			  alphaV = m_distribution.transformRoughness(
-					m_alphaV->eval(bRec.its).average());
+		/* Construct the microfacet distribution matching the
+		   roughness values at the current surface position. */
+		MicrofacetDistribution distr(
+			m_type,
+			m_alphaU->eval(bRec.its).average(),
+			m_alphaV->eval(bRec.its).average(),
+			m_sampleVisible
+		);
 
-		return m_distribution.pdf(H, alphaU, alphaV)
-			/ (4 * absDot(bRec.wo, H));
+		if (m_sampleVisible)
+			return distr.eval(H) * distr.smithG1(bRec.wi, H)
+				/ (4.0f * Frame::cosTheta(bRec.wi));
+		else
+			return distr.pdf(bRec.wi, H) / (4 * absDot(bRec.wo, H));
 	}
 
 	Spectrum sample(BSDFSamplingRecord &bRec, const Point2 &sample) const {
@@ -304,21 +326,23 @@ public:
 			!(bRec.typeMask & EGlossyReflection)))
 			return Spectrum(0.0f);
 
-		/* Evaluate the roughness */
-		Float alphaU = m_distribution.transformRoughness(
-					m_alphaU->eval(bRec.its).average()),
-			  alphaV = m_distribution.transformRoughness(
-					m_alphaV->eval(bRec.its).average());
+		/* Construct the microfacet distribution matching the
+		   roughness values at the current surface position. */
+		MicrofacetDistribution distr(
+			m_type,
+			m_alphaU->eval(bRec.its).average(),
+			m_alphaV->eval(bRec.its).average(),
+			m_sampleVisible
+		);
 
-		/* Sample M, the microsurface normal */
-		Float microfacetPDF;
-		const Normal m = m_distribution.sample(sample,
-			alphaU, alphaV, microfacetPDF);
+		/* Sample M, the microfacet normal */
+		Float pdf;
+		Normal m = distr.sample(bRec.wi, sample, pdf);
 
-		if (microfacetPDF == 0)
+		if (pdf == 0)
 			return Spectrum(0.0f);
 
-		/* Perfect specular reflection based on the microsurface normal */
+		/* Perfect specular reflection based on the microfacet normal */
 		bRec.wo = reflect(bRec.wi, m);
 		bRec.eta = 1.0f;
 		bRec.sampledComponent = 0;
@@ -328,18 +352,18 @@ public:
 		if (Frame::cosTheta(bRec.wo) <= 0)
 			return Spectrum(0.0f);
 
-		const Spectrum F = fresnelConductorExact(dot(bRec.wi, m),
-				m_eta, m_k);
+		Spectrum F = fresnelConductorExact(dot(bRec.wi, m),
+			m_eta, m_k) * m_specularReflectance->eval(bRec.its);
 
-		Float numerator = m_distribution.eval(m, alphaU, alphaV)
-			* m_distribution.G(bRec.wi, bRec.wo, m, alphaU, alphaV)
-			* dot(bRec.wi, m);
+		Float weight;
+		if (m_sampleVisible) {
+			weight = distr.smithG1(bRec.wo, m);
+		} else {
+			weight = distr.eval(m) * distr.G(bRec.wi, bRec.wo, m)
+				* dot(bRec.wi, m) / (pdf * Frame::cosTheta(bRec.wi));
+		}
 
-		Float denominator = microfacetPDF
-			* Frame::cosTheta(bRec.wi);
-
-		return m_specularReflectance->eval(bRec.its) * F
-				* (numerator / denominator);
+		return F * weight;
 	}
 
 	Spectrum sample(BSDFSamplingRecord &bRec, Float &pdf, const Point2 &sample) const {
@@ -348,20 +372,22 @@ public:
 			!(bRec.typeMask & EGlossyReflection)))
 			return Spectrum(0.0f);
 
-		/* Evaluate the roughness */
-		Float alphaU = m_distribution.transformRoughness(
-					m_alphaU->eval(bRec.its).average()),
-			  alphaV = m_distribution.transformRoughness(
-					m_alphaV->eval(bRec.its).average());
+		/* Construct the microfacet distribution matching the
+		   roughness values at the current surface position. */
+		MicrofacetDistribution distr(
+			m_type,
+			m_alphaU->eval(bRec.its).average(),
+			m_alphaV->eval(bRec.its).average(),
+			m_sampleVisible
+		);
 
-		/* Sample M, the microsurface normal */
-		const Normal m = m_distribution.sample(sample,
-			alphaU, alphaV, pdf);
+		/* Sample M, the microfacet normal */
+		Normal m = distr.sample(bRec.wi, sample, pdf);
 
 		if (pdf == 0)
 			return Spectrum(0.0f);
 
-		/* Perfect specular reflection based on the microsurface normal */
+		/* Perfect specular reflection based on the microfacet normal */
 		bRec.wo = reflect(bRec.wi, m);
 		bRec.eta = 1.0f;
 		bRec.sampledComponent = 0;
@@ -371,23 +397,23 @@ public:
 		if (Frame::cosTheta(bRec.wo) <= 0)
 			return Spectrum(0.0f);
 
-		const Spectrum F = fresnelConductorExact(dot(bRec.wi, m),
-				m_eta, m_k);
+		Spectrum F = fresnelConductorExact(dot(bRec.wi, m),
+			m_eta, m_k) * m_specularReflectance->eval(bRec.its);
 
-		Float numerator = m_distribution.eval(m, alphaU, alphaV)
-			* m_distribution.G(bRec.wi, bRec.wo, m, alphaU, alphaV)
-			* dot(bRec.wi, m);
-
-		Float denominator = pdf * Frame::cosTheta(bRec.wi);
+		Float weight;
+		if (m_sampleVisible) {
+			weight = distr.smithG1(bRec.wo, m);
+		} else {
+			weight = distr.eval(m) * distr.G(bRec.wi, bRec.wo, m)
+				* dot(bRec.wi, m) / (pdf * Frame::cosTheta(bRec.wi));
+		}
 
 		/* Jacobian of the half-direction mapping */
 		pdf /= 4.0f * dot(bRec.wo, m);
 
-		return m_specularReflectance->eval(bRec.its) * F
-				* (numerator / denominator);
+		return F * weight;
 	}
 
-
 	void addChild(const std::string &name, ConfigurableObject *child) {
 		if (child->getClass()->derivesFrom(MTS_CLASS(Texture))) {
 			if (name == "alpha")
@@ -405,17 +431,6 @@ public:
 		}
 	}
 
-	void serialize(Stream *stream, InstanceManager *manager) const {
-		BSDF::serialize(stream, manager);
-
-		stream->writeUInt((uint32_t) m_distribution.getType());
-		manager->serialize(stream, m_alphaU.get());
-		manager->serialize(stream, m_alphaV.get());
-		manager->serialize(stream, m_specularReflectance.get());
-		m_eta.serialize(stream);
-		m_k.serialize(stream);
-	}
-
 	Float getRoughness(const Intersection &its, int component) const {
 		return 0.5f * (m_alphaU->eval(its).average()
 			+ m_alphaV->eval(its).average());
@@ -425,7 +440,8 @@ public:
 		std::ostringstream oss;
 		oss << "RoughConductor[" << endl
 			<< "  id = \"" << getID() << "\"," << endl
-			<< "  distribution = " << m_distribution.toString() << "," << endl
+			<< "  distribution = " << MicrofacetDistribution::distributionName(m_type) << "," << endl
+			<< "  sampleVisible = " << m_sampleVisible << "," << endl
 			<< "  alphaU = " << indent(m_alphaU->toString()) << "," << endl
 			<< "  alphaV = " << indent(m_alphaV->toString()) << "," << endl
 			<< "  specularReflectance = " << indent(m_specularReflectance->toString()) << "," << endl
@@ -439,9 +455,10 @@ public:
 
 	MTS_DECLARE_CLASS()
 private:
-	MicrofacetDistribution m_distribution;
+	MicrofacetDistribution::EType m_type;
 	ref<Texture> m_specularReflectance;
 	ref<Texture> m_alphaU, m_alphaV;
+	bool m_sampleVisible;
 	Spectrum m_eta, m_k;
 };
 
diff --git a/src/bsdfs/roughdielectric.cpp b/src/bsdfs/roughdielectric.cpp
index be0e286f..028cdbda 100644
--- a/src/bsdfs/roughdielectric.cpp
+++ b/src/bsdfs/roughdielectric.cpp
@@ -24,12 +24,6 @@
 
 MTS_NAMESPACE_BEGIN
 
-/* Suggestion by Bruce Walter: sample the model using a slightly
-   wider density function. This in practice limits the importance
-   weights to values <= 4.
-*/
-#define ENLARGE_LOBE_TRICK 1
-
 /*!\plugin{roughdielectric}{Rough dielectric material}
  * \order{5}
  * \icon{bsdf_roughdielectric}
@@ -37,42 +31,44 @@ MTS_NAMESPACE_BEGIN
  *     \parameter{distribution}{\String}{
  *          Specifies the type of microfacet normal distribution
  *          used to model the surface roughness.
+ *          \vspace{-1mm}
  *       \begin{enumerate}[(i)]
  *           \item \code{beckmann}: Physically-based distribution derived from
- *               Gaussian random surfaces. This is the default.
- *           \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.
- *           \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.
- *           \item \code{as}: Anisotropic Phong-style microfacet distribution proposed by
- *              Ashikhmin and Shirley \cite{Ashikhmin2005Anisotropic}.\vspace{-3mm}
+ *               Gaussian random surfaces. This is the default.\vspace{-1.5mm}
+ *           \item \code{ggx}: The GGX \cite{Walter07Microfacet} distribution (also known as
+ *               Trowbridge-Reitz \cite{Trowbridge19975Average} distribution)
+ *               was designed to better approximate the long tails observed in measurements
+ *               of ground surfaces, which are not modeled by the Beckmann distribution.
+ *           \vspace{-1.5mm}
+ *           \item \code{phong}: Anisotropic Phong distribution by
+ *              Ashikhmin and Shirley \cite{Ashikhmin2005Anisotropic}.
+ *              In most cases, the \code{ggx} and \code{beckmann} distributions
+ *              should be preferred, since they provide better importance sampling
+ *              and accurate shadowing/masking computations.
+ *              \vspace{-4mm}
  *       \end{enumerate}
  *     }
- *     \parameter{alpha}{\Float\Or\Texture}{
- *         Specifies the roughness of the unresolved surface micro-geometry.
- *         When the Beckmann distribution is used, this parameter is equal to the
- *         \emph{root mean square} (RMS) slope of the microfacets. This
- *         parameter is only valid when \texttt{distribution=beckmann/phong/ggx}.
- *         \default{0.1}.
- *     }
- *     \parameter{alphaU, alphaV}{\Float\Or\Texture}{
- *         Specifies the anisotropic roughness values along the tangent and
- *         bitangent directions. These parameter are only valid when
- *         \texttt{distribution=as}. \default{0.1}.
+ *     \parameter{alpha, alphaU, alphaV}{\Float\Or\Texture}{
+ *         Specifies the roughness of the unresolved surface micro-geometry
+ *         along the tangent and bitangent directions. When the Beckmann
+ *         distribution is used, this parameter is equal to the
+ *         \emph{root mean square} (RMS) slope of the microfacets.
+ *         \code{alpha} is a convenience parameter to initialize both
+ *         \code{alphaU} and \code{alphaV} to the same value. \default{0.1}.
  *     }
  *     \parameter{intIOR}{\Float\Or\String}{Interior index of refraction specified
- *      numerically or using a known material name. \default{\texttt{bk7} / 1.5046}}
+ *         numerically or using a known material name. \default{\texttt{bk7} / 1.5046}}
  *     \parameter{extIOR}{\Float\Or\String}{Exterior index of refraction specified
- *      numerically or using a known material name. \default{\texttt{air} / 1.000277}}
- *     \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}}
- *     \parameter{specular\showbreak Transmittance}{\Spectrum\Or\Texture}{Optional
- *         factor that can be used to modulate the specular transmission component. Note
+ *         numerically or using a known material name. \default{\texttt{air} / 1.000277}}
+ *     \parameter{sampleVisible}{\Boolean}{
+ *         Enables a sampling technique proposed by Heitz and D'Eon~\cite{Heitz1014Importance},
+ *         which focuses computation on the visible parts of the microfacet normal
+ *         distribution, considerably reducing variance in some cases.
+ *         \default{\code{true}, i.e. use visible normal sampling}
+ *     }
+ *     \parameter{specular\showbreak Reflectance,\newline
+ *         specular\showbreak Transmittance}{\Spectrum\Or\Texture}{Optional
+ *         factor that can be used to modulate the specular reflection/transmission component. Note
  *         that for physical realism, this parameter should never be touched. \default{1.0}}
  * }\vspace{4mm}
  *
@@ -98,7 +94,7 @@ MTS_NAMESPACE_BEGIN
  *
  * The implementation is based on the paper ``Microfacet Models
  * for Refraction through Rough Surfaces'' by Walter et al.
- * \cite{Walter07Microfacet}. It supports several different types of microfacet
+ * \cite{Walter07Microfacet}. It supports three different types of microfacet
  * distributions and has a texturable roughness parameter. Exterior and
  * interior IOR values can be specified independently, where ``exterior''
  * refers to the side that contains the surface normal. Similar to the
@@ -109,13 +105,12 @@ MTS_NAMESPACE_BEGIN
  * glass BK7/air interface with a light amount of roughness modeled using a
  * Beckmann distribution.
  *
- * To get an intuition about the effect of the surface roughness
- * parameter $\alpha$, consider the following approximate classification:
- * a value of $\alpha=0.001-0.01$ corresponds to a material
- * with slight imperfections on an
- * otherwise smooth surface finish, $\alpha=0.1$ is relatively rough,
+ * To get an intuition about the effect of the surface roughness parameter
+ * $\alpha$, consider the following approximate classification: a value of
+ * $\alpha=0.001-0.01$ corresponds to a material with slight imperfections
+ * on an otherwise smooth surface finish, $\alpha=0.1$ is relatively rough,
  * and $\alpha=0.3-0.7$ is \emph{extremely} rough (e.g. an etched or ground
- * finish).
+ * finish). Values significantly above that are probably not too realistic.
  *
  * Please note that when using this plugin, it is crucial that the scene contains
  * meaningful and mutually compatible index of refraction changes---see
@@ -127,20 +122,33 @@ MTS_NAMESPACE_BEGIN
  * converge slowly.
  *
  * \subsubsection*{Technical details}
- * When rendering with the Ashikhmin-Shirley or Phong microfacet
- * distributions, a conversion is used to turn the specified
- * $\alpha$ roughness value into the exponents of these distributions.
- * This is done in a way, such that the different
- * distributions all produce a similar appearance for the same value of
- * $\alpha$.
+ * All microfacet distributions allow the specification of two distinct
+ * roughness values along the tangent and bitangent directions. This can be
+ * used to provide a material with a ``brushed'' appearance. The alignment
+ * of the anisotropy will follow the UV parameterization of the underlying
+ * mesh. This means that such an anisotropic material cannot be applied to
+ * triangle meshes that are missing texture coordinates.
  *
- * The Ashikhmin-Shirley microfacet distribution allows the specification
- * of two distinct roughness values along the tangent and bitangent
- * directions. This can be used to provide a material with a ``brushed''
- * appearance. The alignment of the anisotropy will follow the UV
- * parameterization of the underlying mesh in this case. This also means that
- * such an anisotropic material cannot be applied to triangle meshes that
- * are missing texture coordinates.\newpage
+ * Since Mitsuba 0.5.1, this plugin uses a new importance sampling technique
+ * contributed by Eric Heitz and Eugene D'Eon, which restricts the sampling
+ * domain to the set of visible (unmasked) microfacet normals. The previous
+ * approach of sampling all normals is still available and can be enabled
+ * by setting \code{sampleVisible} to \code{false}.
+ * Note that this new method is only available for the \code{beckmann} and
+ * \code{ggx} microfacet distributions. When the \code{phong} distribution
+ * is selected, the parameter has no effect.
+ *
+ * When rendering with the Phong microfacet distribution, a conversion is
+ * used to turn the specified Beckmann-equivalent $\alpha$ roughness value
+ * into the exponent parameter of this distribution. This is done in a way,
+ * such that the same value $\alpha$ will produce a similar appearance across
+ * different microfacet distributions.
+ *
+ * When rendering with the Phong microfacet distribution, a conversion is
+ * used to turn the specified Beckmann-equivalent $\alpha$ roughness value
+ * into the exponents of the distribution. This is done in a way, such that
+ * the different distributions all produce a similar appearance for the
+ * same value of $\alpha$.
  *
  * \renderings{
  *     \rendering{Ground glass (GGX, $\alpha$=0.304,
@@ -191,26 +199,21 @@ public:
 		m_eta = intIOR / extIOR;
 		m_invEta = 1 / m_eta;
 
-		m_distribution = MicrofacetDistribution(
-			props.getString("distribution", "beckmann")
-		);
+		MicrofacetDistribution distr(props);
+		m_type = distr.getType();
+		m_sampleVisible = distr.getSampleVisible();
 
-		Float alpha = props.getFloat("alpha", 0.1f),
-			  alphaU = props.getFloat("alphaU", alpha),
-			  alphaV = props.getFloat("alphaV", alpha);
-
-		m_alphaU = new ConstantFloatTexture(alphaU);
-		if (alphaU == alphaV)
+		m_alphaU = new ConstantFloatTexture(distr.getAlphaU());
+		if (distr.getAlphaU() == distr.getAlphaV())
 			m_alphaV = m_alphaU;
 		else
-			m_alphaV = new ConstantFloatTexture(alphaV);
+			m_alphaV = new ConstantFloatTexture(distr.getAlphaV());
 	}
 
 	RoughDielectric(Stream *stream, InstanceManager *manager)
 	 : BSDF(stream, manager) {
-		m_distribution = MicrofacetDistribution(
-			(MicrofacetDistribution::EType) stream->readUInt()
-		);
+		m_type = (MicrofacetDistribution::EType) stream->readUInt();
+		m_sampleVisible = stream->readBool();
 		m_alphaU = static_cast<Texture *>(manager->getInstance(stream));
 		m_alphaV = static_cast<Texture *>(manager->getInstance(stream));
 		m_specularReflectance = static_cast<Texture *>(manager->getInstance(stream));
@@ -221,17 +224,22 @@ public:
 		configure();
 	}
 
+	void serialize(Stream *stream, InstanceManager *manager) const {
+		BSDF::serialize(stream, manager);
+
+		stream->writeUInt((uint32_t) m_type);
+		stream->writeBool(m_sampleVisible);
+		manager->serialize(stream, m_alphaU.get());
+		manager->serialize(stream, m_alphaV.get());
+		manager->serialize(stream, m_specularReflectance.get());
+		manager->serialize(stream, m_specularTransmittance.get());
+		stream->writeFloat(m_eta);
+	}
+
 	void configure() {
 		unsigned int extraFlags = 0;
-		if (m_alphaU != m_alphaV) {
+		if (m_alphaU != m_alphaV)
 			extraFlags |= EAnisotropic;
-			if (m_distribution.getType() !=
-				MicrofacetDistribution::EAshikhminShirley)
-				Log(EError, "Different roughness values along the tangent and "
-						"bitangent directions are only supported when using the "
-						"anisotropic Ashikhmin-Shirley microfacet distribution "
-						"(named \"as\")");
-		}
 
 		if (!m_alphaU->isConstant() || !m_alphaV->isConstant())
 			extraFlags |= ESpatiallyVarying;
@@ -293,14 +301,17 @@ public:
 		   same hemisphere as the macrosurface normal */
 		H *= math::signum(Frame::cosTheta(H));
 
-		/* Evaluate the roughness */
-		Float alphaU = m_distribution.transformRoughness(
-					m_alphaU->eval(bRec.its).average()),
-			  alphaV = m_distribution.transformRoughness(
-					m_alphaV->eval(bRec.its).average());
+		/* Construct the microfacet distribution matching the
+		   roughness values at the current surface position. */
+		MicrofacetDistribution distr(
+			m_type,
+			m_alphaU->eval(bRec.its).average(),
+			m_alphaV->eval(bRec.its).average(),
+			m_sampleVisible
+		);
 
-		/* Evaluate the microsurface normal distribution */
-		const Float D = m_distribution.eval(H, alphaU, alphaV);
+		/* Evaluate the microfacet normal distribution */
+		const Float D = distr.eval(H);
 		if (D == 0)
 			return Spectrum(0.0f);
 
@@ -308,7 +319,7 @@ public:
 		const Float F = fresnelDielectricExt(dot(bRec.wi, H), m_eta);
 
 		/* Smith's shadow-masking function */
-		const Float G = m_distribution.G(bRec.wi, bRec.wo, H, alphaU, alphaV);
+		const Float G = distr.G(bRec.wi, bRec.wo, H);
 
 		if (reflect) {
 			/* Calculate the total amount of reflection */
@@ -383,21 +394,24 @@ public:
 		   same hemisphere as the macrosurface normal */
 		H *= math::signum(Frame::cosTheta(H));
 
-		/* Evaluate the roughness */
-		Float alphaU = m_alphaU->eval(bRec.its).average(),
-			  alphaV = m_alphaV->eval(bRec.its).average();
+		/* Construct the microfacet distribution matching the
+		   roughness values at the current surface position. */
+		MicrofacetDistribution sampleDistr(
+			m_type,
+			m_alphaU->eval(bRec.its).average(),
+			m_alphaV->eval(bRec.its).average(),
+			m_sampleVisible
+		);
 
-#if ENLARGE_LOBE_TRICK == 1
-		Float factor = (1.2f - 0.2f * std::sqrt(
-			std::abs(Frame::cosTheta(bRec.wi))));
-		alphaU *= factor; alphaV *= factor;
-#endif
+		/* Trick by Walter et al.: slightly scale the roughness values to
+		   reduce importance sampling weights. Not needed for the
+		   Heitz and D'Eon sampling technique. */
+		if (!m_sampleVisible)
+			sampleDistr.scaleAlpha(1.2f - 0.2f * std::sqrt(
+				std::abs(Frame::cosTheta(bRec.wi))));
 
-		alphaU = m_distribution.transformRoughness(alphaU);
-		alphaV = m_distribution.transformRoughness(alphaV);
-
-		/* Evaluate the microsurface normal sampling density */
-		Float prob = m_distribution.pdf(H, alphaU, alphaV);
+		/* Evaluate the microfacet model sampling density function */
+		Float prob = sampleDistr.pdf(math::signum(Frame::cosTheta(bRec.wi)) * bRec.wi, H);
 
 		if (hasTransmission && hasReflection) {
 			Float F = fresnelDielectricExt(dot(bRec.wi, H), m_eta);
@@ -419,44 +433,42 @@ public:
 		if (!hasReflection && !hasTransmission)
 			return Spectrum(0.0f);
 
-		/* Evaluate the roughness */
-		Float alphaU = m_alphaU->eval(bRec.its).average(),
-		      alphaV = m_alphaV->eval(bRec.its).average(),
-		      sampleAlphaU = alphaU,
-		      sampleAlphaV = alphaV;
+		/* Construct the microfacet distribution matching the
+		   roughness values at the current surface position. */
+		MicrofacetDistribution distr(
+			m_type,
+			m_alphaU->eval(bRec.its).average(),
+			m_alphaV->eval(bRec.its).average(),
+			m_sampleVisible
+		);
 
-#if ENLARGE_LOBE_TRICK == 1
-		Float factor = (1.2f - 0.2f * std::sqrt(
-			std::abs(Frame::cosTheta(bRec.wi))));
-		sampleAlphaU *= factor; sampleAlphaV *= factor;
-#endif
+		/* Trick by Walter et al.: slightly scale the roughness values to
+		   reduce importance sampling weights. Not needed for the
+		   Heitz and D'Eon sampling technique. */
+		MicrofacetDistribution sampleDistr(distr);
+		if (!m_sampleVisible)
+			sampleDistr.scaleAlpha(1.2f - 0.2f * std::sqrt(
+				std::abs(Frame::cosTheta(bRec.wi))));
 
-		alphaU = m_distribution.transformRoughness(alphaU);
-		alphaV = m_distribution.transformRoughness(alphaV);
-		sampleAlphaU = m_distribution.transformRoughness(sampleAlphaU);
-		sampleAlphaV = m_distribution.transformRoughness(sampleAlphaV);
-
-		/* Sample M, the microsurface normal */
+		/* Sample M, the microfacet normal */
 		Float microfacetPDF;
-		const Normal m = m_distribution.sample(sample,
-				sampleAlphaU, sampleAlphaV, microfacetPDF);
-
+		const Normal m = sampleDistr.sample(math::signum(Frame::cosTheta(bRec.wi)) * bRec.wi, sample, microfacetPDF);
 		if (microfacetPDF == 0)
 			return Spectrum(0.0f);
 
-		Float cosThetaT, numerator = 1.0f;
+		Float cosThetaT;
 		Float F = fresnelDielectricExt(dot(bRec.wi, m), cosThetaT, m_eta);
+		Spectrum weight(1.0f);
 
 		if (hasReflection && hasTransmission) {
 			if (bRec.sampler->next1D() > F)
 				sampleReflection = false;
 		} else {
-			numerator = hasReflection ? F : (1-F);
+			weight = Spectrum(hasReflection ? F : (1-F));
 		}
 
-		Spectrum result;
 		if (sampleReflection) {
-			/* Perfect specular reflection based on the microsurface normal */
+			/* Perfect specular reflection based on the microfacet normal */
 			bRec.wo = reflect(bRec.wi, m);
 			bRec.eta = 1.0f;
 			bRec.sampledComponent = 0;
@@ -466,12 +478,12 @@ public:
 			if (Frame::cosTheta(bRec.wi) * Frame::cosTheta(bRec.wo) <= 0)
 				return Spectrum(0.0f);
 
-			result = m_specularReflectance->eval(bRec.its);
+			weight *= m_specularReflectance->eval(bRec.its);
 		} else {
 			if (cosThetaT == 0)
 				return Spectrum(0.0f);
 
-			/* Perfect specular transmission based on the microsurface normal */
+			/* Perfect specular transmission based on the microfacet normal */
 			bRec.wo = refract(bRec.wi, m, m_eta, cosThetaT);
 			bRec.eta = cosThetaT < 0 ? m_eta : m_invEta;
 			bRec.sampledComponent = 1;
@@ -486,17 +498,16 @@ public:
 			Float factor = (bRec.mode == ERadiance)
 				? (cosThetaT < 0 ? m_invEta : m_eta) : 1.0f;
 
-			result = m_specularTransmittance->eval(bRec.its) * (factor * factor);
+			weight *= m_specularTransmittance->eval(bRec.its) * (factor * factor);
 		}
 
-		numerator *= m_distribution.eval(m, alphaU, alphaV)
-			* m_distribution.G(bRec.wi, bRec.wo, m, alphaU, alphaV)
-			* dot(bRec.wi, m);
+		if (m_sampleVisible)
+			weight *= distr.smithG1(bRec.wo, m);
+		else
+			weight *= std::abs(distr.eval(m) * distr.G(bRec.wi, bRec.wo, m)
+				* dot(bRec.wi, m) / (microfacetPDF * Frame::cosTheta(bRec.wi)));
 
-		Float denominator = microfacetPDF
-			* Frame::cosTheta(bRec.wi);
-
-		return result * std::abs(numerator / denominator);
+		return weight;
 	}
 
 	Spectrum sample(BSDFSamplingRecord &bRec, Float &pdf, const Point2 &_sample) const {
@@ -511,35 +522,33 @@ public:
 		if (!hasReflection && !hasTransmission)
 			return Spectrum(0.0f);
 
-		/* Evaluate the roughness */
-		Float alphaU = m_alphaU->eval(bRec.its).average(),
-		      alphaV = m_alphaV->eval(bRec.its).average(),
-		      sampleAlphaU = alphaU,
-		      sampleAlphaV = alphaV;
+		/* Construct the microfacet distribution matching the
+		   roughness values at the current surface position. */
+		MicrofacetDistribution distr(
+			m_type,
+			m_alphaU->eval(bRec.its).average(),
+			m_alphaV->eval(bRec.its).average(),
+			m_sampleVisible
+		);
 
-#if ENLARGE_LOBE_TRICK == 1
-		Float factor = (1.2f - 0.2f * std::sqrt(
-			std::abs(Frame::cosTheta(bRec.wi))));
-		sampleAlphaU *= factor; sampleAlphaV *= factor;
-#endif
+		/* Trick by Walter et al.: slightly scale the roughness values to
+		   reduce importance sampling weights. Not needed for the
+		   Heitz and D'Eon sampling technique. */
+		MicrofacetDistribution sampleDistr(distr);
+		if (!m_sampleVisible)
+			sampleDistr.scaleAlpha(1.2f - 0.2f * std::sqrt(
+				std::abs(Frame::cosTheta(bRec.wi))));
 
-		alphaU = m_distribution.transformRoughness(alphaU);
-		alphaV = m_distribution.transformRoughness(alphaV);
-		sampleAlphaU = m_distribution.transformRoughness(sampleAlphaU);
-		sampleAlphaV = m_distribution.transformRoughness(sampleAlphaV);
-
-		/* Sample M, the microsurface normal */
+		/* Sample M, the microfacet normal */
 		Float microfacetPDF;
-		const Normal m = m_distribution.sample(sample,
-				sampleAlphaU, sampleAlphaV, microfacetPDF);
-
+		const Normal m = sampleDistr.sample(math::signum(Frame::cosTheta(bRec.wi)) * bRec.wi, sample, microfacetPDF);
 		if (microfacetPDF == 0)
 			return Spectrum(0.0f);
-
 		pdf = microfacetPDF;
 
-		Float cosThetaT, numerator = 1.0f;
+		Float cosThetaT;
 		Float F = fresnelDielectricExt(dot(bRec.wi, m), cosThetaT, m_eta);
+		Spectrum weight(1.0f);
 
 		if (hasReflection && hasTransmission) {
 			if (bRec.sampler->next1D() > F) {
@@ -549,14 +558,12 @@ public:
 				pdf *= F;
 			}
 		} else {
-			numerator = hasReflection ? F : (1-F);
+			weight *= hasReflection ? F : (1-F);
 		}
 
-		Spectrum result;
 		Float dwh_dwo;
-
 		if (sampleReflection) {
-			/* Perfect specular reflection based on the microsurface normal */
+			/* Perfect specular reflection based on the microfacet normal */
 			bRec.wo = reflect(bRec.wi, m);
 			bRec.eta = 1.0f;
 			bRec.sampledComponent = 0;
@@ -566,7 +573,7 @@ public:
 			if (Frame::cosTheta(bRec.wi) * Frame::cosTheta(bRec.wo) <= 0)
 				return Spectrum(0.0f);
 
-			result = m_specularReflectance->eval(bRec.its);
+			weight *= m_specularReflectance->eval(bRec.its);
 
 			/* Jacobian of the half-direction mapping */
 			dwh_dwo = 1.0f / (4.0f * dot(bRec.wo, m));
@@ -574,7 +581,7 @@ public:
 			if (cosThetaT == 0)
 				return Spectrum(0.0f);
 
-			/* Perfect specular transmission based on the microsurface normal */
+			/* Perfect specular transmission based on the microfacet normal */
 			bRec.wo = refract(bRec.wi, m, m_eta, cosThetaT);
 			bRec.eta = cosThetaT < 0 ? m_eta : m_invEta;
 			bRec.sampledComponent = 1;
@@ -589,22 +596,22 @@ public:
 			Float factor = (bRec.mode == ERadiance)
 				? (cosThetaT < 0 ? m_invEta : m_eta) : 1.0f;
 
-			result = m_specularTransmittance->eval(bRec.its) * (factor * factor);
+			weight *= m_specularTransmittance->eval(bRec.its) * (factor * factor);
 
 			/* Jacobian of the half-direction mapping */
 			Float sqrtDenom = dot(bRec.wi, m) + bRec.eta * dot(bRec.wo, m);
 			dwh_dwo = (bRec.eta*bRec.eta * dot(bRec.wo, m)) / (sqrtDenom*sqrtDenom);
 		}
 
-		numerator *= m_distribution.eval(m, alphaU, alphaV)
-			* m_distribution.G(bRec.wi, bRec.wo, m, alphaU, alphaV)
-			* dot(bRec.wi, m);
-
-		Float denominator = microfacetPDF * Frame::cosTheta(bRec.wi);
+		if (m_sampleVisible)
+			weight *= distr.smithG1(bRec.wo, m);
+		else
+			weight *= std::abs(distr.eval(m) * distr.G(bRec.wi, bRec.wo, m)
+				* dot(bRec.wi, m) / (microfacetPDF * Frame::cosTheta(bRec.wi)));
 
 		pdf *= std::abs(dwh_dwo);
 
-		return result * std::abs(numerator / denominator);
+		return weight;
 	}
 
 	void addChild(const std::string &name, ConfigurableObject *child) {
@@ -626,17 +633,6 @@ public:
 		}
 	}
 
-	void serialize(Stream *stream, InstanceManager *manager) const {
-		BSDF::serialize(stream, manager);
-
-		stream->writeUInt((uint32_t) m_distribution.getType());
-		manager->serialize(stream, m_alphaU.get());
-		manager->serialize(stream, m_alphaV.get());
-		manager->serialize(stream, m_specularReflectance.get());
-		manager->serialize(stream, m_specularTransmittance.get());
-		stream->writeFloat(m_eta);
-	}
-
 	Float getEta() const {
 		return m_eta;
 	}
@@ -650,7 +646,8 @@ public:
 		std::ostringstream oss;
 		oss << "RoughDielectric[" << endl
 			<< "  id = \"" << getID() << "\"," << endl
-			<< "  distribution = " << m_distribution.toString() << "," << endl
+			<< "  distribution = " << MicrofacetDistribution::distributionName(m_type) << "," << endl
+			<< "  sampleVisible = " << m_sampleVisible << "," << endl
 			<< "  eta = " << m_eta << "," << endl
 			<< "  alphaU = " << indent(m_alphaU->toString()) << "," << endl
 			<< "  alphaV = " << indent(m_alphaV->toString()) << "," << endl
@@ -664,11 +661,12 @@ public:
 
 	MTS_DECLARE_CLASS()
 private:
-	MicrofacetDistribution m_distribution;
+	MicrofacetDistribution::EType m_type;
 	ref<Texture> m_specularTransmittance;
 	ref<Texture> m_specularReflectance;
 	ref<Texture> m_alphaU, m_alphaV;
 	Float m_eta, m_invEta;
+	bool m_sampleVisible;
 };
 
 /* Fake glass shader -- it is really hopeless to visualize
diff --git a/src/bsdfs/roughplastic.cpp b/src/bsdfs/roughplastic.cpp
index 9747339b..1c18c00a 100644
--- a/src/bsdfs/roughplastic.cpp
+++ b/src/bsdfs/roughplastic.cpp
@@ -32,18 +32,20 @@ MTS_NAMESPACE_BEGIN
  *     \parameter{distribution}{\String}{
  *          Specifies the type of microfacet normal distribution
  *          used to model the surface roughness.
+ *          \vspace{-1mm}
  *       \begin{enumerate}[(i)]
  *           \item \code{beckmann}: Physically-based distribution derived from
- *               Gaussian random surfaces. This is the default.
- *           \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.
- *           \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.
- *              \vspace{-3mm}
+ *               Gaussian random surfaces. This is the default.\vspace{-1.5mm}
+ *           \item \code{ggx}: The GGX \cite{Walter07Microfacet} distribution (also known as
+ *               Trowbridge-Reitz \cite{Trowbridge19975Average} distribution)
+ *               was designed to better approximate the long tails observed in measurements
+ *               of ground surfaces, which are not modeled by the Beckmann distribution.
+ *           \vspace{-1.5mm}
+ *           \item \code{phong}: Classical Phong distribution.
+ *              In most cases, the \code{ggx} and \code{beckmann} distributions
+ *              should be preferred, since they provide better importance sampling
+ *              and accurate shadowing/masking computations.
+ *              \vspace{-4mm}
  *       \end{enumerate}
  *     }
  *     \parameter{alpha}{\Float\Or\Texture}{
@@ -52,10 +54,16 @@ MTS_NAMESPACE_BEGIN
  *         \emph{root mean square} (RMS) slope of the microfacets.
  *         \default{0.1}.
  *     }
+ *
  *     \parameter{intIOR}{\Float\Or\String}{Interior index of refraction specified
  *      numerically or using a known material name. \default{\texttt{polypropylene} / 1.49}}
  *     \parameter{extIOR}{\Float\Or\String}{Exterior index of refraction specified
  *      numerically or using a known material name. \default{\texttt{air} / 1.000277}}
+ *     \parameter{sampleVisible}{\Boolean}{
+ *	       Enables an improved importance sampling technique. Refer to
+ *	       pages \pageref{plg:roughconductor} and \pageref{sec:visiblenormal-sampling}
+ *	       for details. \default{\code{true}}
+ *     }
  *     \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}}
@@ -75,12 +83,12 @@ MTS_NAMESPACE_BEGIN
  * preferred over heuristic models like \pluginref{phong} and \pluginref{ward}
  * when possible.
  *
- * Microfacet theory describes rough surfaces as an arrangement of unresolved and
- * ideally specular facets, whose normal directions are given by a specially
- * chosen \emph{microfacet distribution}.
- * By accounting for shadowing and masking effects between these facets, it is
- * possible to reproduce the important off-specular reflections peaks observed
- * in real-world measurements of such materials.
+ * Microfacet theory describes rough surfaces as an arrangement of
+ * unresolved and ideally specular facets, whose normal directions are given by
+ * a specially chosen \emph{microfacet distribution}. By accounting for shadowing
+ * and masking effects between these facets, it is possible to reproduce the important
+ * off-specular reflections peaks observed in real-world measurements of such
+ * materials.
  *
  * \renderings{
  *     \rendering{Beckmann, $\alpha=0.1$}{bsdf_roughplastic_beckmann}
@@ -112,12 +120,10 @@ MTS_NAMESPACE_BEGIN
  * of this parameter given in the \pluginref{plastic} plugin section
  * on \pluginpage{plastic}.
  *
- *
- * To get an intuition about the effect of the surface roughness
- * parameter $\alpha$, consider the following approximate classification:
- * a value of $\alpha=0.001-0.01$ corresponds to a material
- * with slight imperfections on an
- * otherwise smooth surface finish, $\alpha=0.1$ is relatively rough,
+ * To get an intuition about the effect of the surface roughness parameter
+ * $\alpha$, consider the following approximate classification: a value of
+ * $\alpha=0.001-0.01$ corresponds to a material with slight imperfections
+ * on an otherwise smooth surface finish, $\alpha=0.1$ is relatively rough,
  * and $\alpha=0.3-0.7$ is \emph{extremely} rough (e.g. an etched or ground
  * finish). Values significantly above that are probably not too realistic.
  *
@@ -133,7 +139,6 @@ MTS_NAMESPACE_BEGIN
  *        \pluginref{roughplastic} and support for nonlinear color shifts.
  *     }
  * }
- * \newpage
  * \renderings{
  *     \rendering{Wood material with smooth horizontal stripes}{bsdf_roughplastic_roughtex1}
  *     \rendering{A material with imperfections at a much smaller scale than what
@@ -181,11 +186,11 @@ MTS_NAMESPACE_BEGIN
  * values of this function over a large range of parameter values. At runtime,
  * the relevant parts are extracted using tricubic interpolation.
  *
- * When rendering with the Phong microfacet distributions, a conversion
- * is used to turn the specified $\alpha$ roughness value into the Phong
- * exponent. This is done in a way, such that the different distributions
- * all produce a similar appearance for the same value of $\alpha$.
- *
+ * When rendering with the Phong microfacet distribution, a conversion is
+ * used to turn the specified Beckmann-equivalent $\alpha$ roughness value
+ * into the exponent parameter of this distribution. This is done in a way,
+ * such that the same value $\alpha$ will produce a similar appearance across
+ * different microfacet distributions.
  */
 class RoughPlastic : public BSDF {
 public:
@@ -207,27 +212,25 @@ public:
 
 		m_eta = intIOR / extIOR;
 
-		m_distribution = MicrofacetDistribution(
-			props.getString("distribution", "beckmann")
-		);
+		m_nonlinear = props.getBoolean("nonlinear", false);
 
-		if (m_distribution.isAnisotropic())
+		MicrofacetDistribution distr(props);
+		m_type = distr.getType();
+		m_sampleVisible = distr.getSampleVisible();
+
+		if (distr.isAnisotropic())
 			Log(EError, "The 'roughplastic' plugin currently does not support "
 				"anisotropic microfacet distributions!");
 
-		m_nonlinear = props.getBoolean("nonlinear", false);
-
-		m_alpha = new ConstantFloatTexture(
-			props.getFloat("alpha", 0.1f));
+		m_alpha = new ConstantFloatTexture(distr.getAlpha());
 
 		m_specularSamplingWeight = 0.0f;
 	}
 
 	RoughPlastic(Stream *stream, InstanceManager *manager)
 	 : BSDF(stream, manager) {
-		m_distribution = MicrofacetDistribution(
-			(MicrofacetDistribution::EType) stream->readUInt()
-		);
+		m_type = (MicrofacetDistribution::EType) stream->readUInt();
+		m_sampleVisible = stream->readBool();
 		m_specularReflectance = static_cast<Texture *>(manager->getInstance(stream));
 		m_diffuseReflectance = static_cast<Texture *>(manager->getInstance(stream));
 		m_alpha = static_cast<Texture *>(manager->getInstance(stream));
@@ -240,7 +243,8 @@ public:
 	void serialize(Stream *stream, InstanceManager *manager) const {
 		BSDF::serialize(stream, manager);
 
-		stream->writeUInt((uint32_t) m_distribution.getType());
+		stream->writeUInt((uint32_t) m_type);
+		stream->writeBool(m_sampleVisible);
 		manager->serialize(stream, m_specularReflectance.get());
 		manager->serialize(stream, m_diffuseReflectance.get());
 		manager->serialize(stream, m_alpha.get());
@@ -277,8 +281,7 @@ public:
 		if (!m_externalRoughTransmittance.get()) {
 			/* Load precomputed data used to compute the rough
 			   transmittance through the dielectric interface */
-			m_externalRoughTransmittance = new RoughTransmittance(
-				m_distribution.getType());
+			m_externalRoughTransmittance = new RoughTransmittance(m_type);
 
 			m_externalRoughTransmittance->checkEta(m_eta);
 			m_externalRoughTransmittance->checkAlpha(m_alpha->getMinimum().average());
@@ -332,23 +335,27 @@ public:
 			(!hasSpecular && !hasDiffuse))
 			return Spectrum(0.0f);
 
-		/* Evaluate the roughness texture */
-		Float alpha = m_alpha->eval(bRec.its).average();
-		Float alphaT = m_distribution.transformRoughness(alpha);
+		/* Construct the microfacet distribution matching the
+		   roughness values at the current surface position. */
+		MicrofacetDistribution distr(
+			m_type,
+			m_alpha->eval(bRec.its).average(),
+			m_sampleVisible
+		);
 
 		Spectrum result(0.0f);
 		if (hasSpecular) {
 			/* Calculate the reflection half-vector */
 			const Vector H = normalize(bRec.wo+bRec.wi);
 
-			/* Evaluate the microsurface normal distribution */
-			const Float D = m_distribution.eval(H, alphaT);
+			/* Evaluate the microfacet normal distribution */
+			const Float D = distr.eval(H);
 
 			/* Fresnel term */
 			const Float F = fresnelDielectricExt(dot(bRec.wi, H), m_eta);
 
 			/* Smith's shadow-masking function */
-			const Float G = m_distribution.G(bRec.wi, bRec.wo, H, alphaT);
+			const Float G = distr.G(bRec.wi, bRec.wo, H);
 
 			/* Calculate the specular reflection component */
 			Float value = F * D * G /
@@ -359,9 +366,9 @@ public:
 
 		if (hasDiffuse) {
 			Spectrum diff = m_diffuseReflectance->eval(bRec.its);
-			Float T12 = m_externalRoughTransmittance->eval(Frame::cosTheta(bRec.wi), alpha);
-			Float T21 = m_externalRoughTransmittance->eval(Frame::cosTheta(bRec.wo), alpha);
-			Float Fdr = 1-m_internalRoughTransmittance->evalDiffuse(alpha);
+			Float T12 = m_externalRoughTransmittance->eval(Frame::cosTheta(bRec.wi), distr.getAlpha());
+			Float T21 = m_externalRoughTransmittance->eval(Frame::cosTheta(bRec.wo), distr.getAlpha());
+			Float Fdr = 1-m_internalRoughTransmittance->evalDiffuse(distr.getAlpha());
 
 			if (m_nonlinear)
 				diff /= Spectrum(1.0f) - diff * Fdr;
@@ -386,9 +393,13 @@ public:
 			(!hasSpecular && !hasDiffuse))
 			return 0.0f;
 
-		/* Evaluate the roughness texture */
-		Float alpha = m_alpha->eval(bRec.its).average();
-		Float alphaT = m_distribution.transformRoughness(alpha);
+		/* Construct the microfacet distribution matching the
+		   roughness values at the current surface position. */
+		MicrofacetDistribution distr(
+			m_type,
+			m_alpha->eval(bRec.its).average(),
+			m_sampleVisible
+		);
 
 		/* Calculate the reflection half-vector */
 		const Vector H = normalize(bRec.wo+bRec.wi);
@@ -396,7 +407,7 @@ public:
 		Float probDiffuse, probSpecular;
 		if (hasSpecular && hasDiffuse) {
 			/* Find the probability of sampling the specular component */
-			probSpecular = 1-m_externalRoughTransmittance->eval(Frame::cosTheta(bRec.wi), alpha);
+			probSpecular = 1-m_externalRoughTransmittance->eval(Frame::cosTheta(bRec.wi), distr.getAlpha());
 
 			/* Reallocate samples */
 			probSpecular = (probSpecular*m_specularSamplingWeight) /
@@ -413,8 +424,8 @@ public:
 			/* Jacobian of the half-direction mapping */
 			const Float dwh_dwo = 1.0f / (4.0f * dot(bRec.wo, H));
 
-			/* Evaluate the microsurface normal distribution */
-			const Float prob = m_distribution.pdf(H, alphaT);
+			/* Evaluate the microfacet model sampling density function */
+			const Float prob = distr.pdf(bRec.wi, H);
 
 			result = prob * dwh_dwo * probSpecular;
 		}
@@ -437,14 +448,18 @@ public:
 		bool choseSpecular = hasSpecular;
 		Point2 sample(_sample);
 
-		/* Evaluate the roughness texture */
-		Float alpha = m_alpha->eval(bRec.its).average();
-		Float alphaT = m_distribution.transformRoughness(alpha);
+		/* Construct the microfacet distribution matching the
+		   roughness values at the current surface position. */
+		MicrofacetDistribution distr(
+			m_type,
+			m_alpha->eval(bRec.its).average(),
+			m_sampleVisible
+		);
 
 		Float probSpecular;
 		if (hasSpecular && hasDiffuse) {
 			/* Find the probability of sampling the specular component */
-			probSpecular = 1 - m_externalRoughTransmittance->eval(Frame::cosTheta(bRec.wi), alpha);
+			probSpecular = 1 - m_externalRoughTransmittance->eval(Frame::cosTheta(bRec.wi), distr.getAlpha());
 
 			/* Reallocate samples */
 			probSpecular = (probSpecular*m_specularSamplingWeight) /
@@ -460,8 +475,8 @@ public:
 		}
 
 		if (choseSpecular) {
-			/* Perfect specular reflection based on the microsurface normal */
-			Normal m = m_distribution.sample(sample, alphaT);
+			/* Perfect specular reflection based on the microfacet normal */
+			Normal m = distr.sample(bRec.wi, sample);
 			bRec.wo = reflect(bRec.wi, m);
 			bRec.sampledComponent = 0;
 			bRec.sampledType = EGlossyReflection;
@@ -518,7 +533,8 @@ public:
 		std::ostringstream oss;
 		oss << "RoughPlastic[" << endl
 			<< "  id = \"" << getID() << "\"," << endl
-			<< "  distribution = " << m_distribution.toString() << "," << endl
+			<< "  distribution = " << MicrofacetDistribution::distributionName(m_type) << "," << endl
+			<< "  sampleVisible = " << m_sampleVisible << "," << endl
 			<< "  alpha = " << indent(m_alpha->toString()) << "," << endl
 			<< "  specularReflectance = " << indent(m_specularReflectance->toString()) << "," << endl
 			<< "  diffuseReflectance = " << indent(m_diffuseReflectance->toString()) << "," << endl
@@ -534,7 +550,7 @@ public:
 
 	MTS_DECLARE_CLASS()
 private:
-	MicrofacetDistribution m_distribution;
+	MicrofacetDistribution::EType m_type;
 	ref<RoughTransmittance> m_externalRoughTransmittance;
 	ref<RoughTransmittance> m_internalRoughTransmittance;
 	ref<Texture> m_diffuseReflectance;
@@ -543,6 +559,7 @@ private:
 	Float m_eta, m_invEta2;
 	Float m_specularSamplingWeight;
 	bool m_nonlinear;
+	bool m_sampleVisible;
 };
 
 /**
diff --git a/src/tests/test_microfacet.cpp b/src/tests/test_microfacet.cpp
new file mode 100644
index 00000000..434c9d66
--- /dev/null
+++ b/src/tests/test_microfacet.cpp
@@ -0,0 +1,176 @@
+/*
+    This file is part of Mitsuba, a physically based rendering system.
+
+    Copyright (c) 2007-2014 by Wenzel Jakob and others.
+
+    Mitsuba is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License Version 3
+    as published by the Free Software Foundation.
+
+    Mitsuba is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program. If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <mitsuba/core/plugin.h>
+#include <mitsuba/core/statistics.h>
+#include <mitsuba/core/chisquare.h>
+#include <mitsuba/core/fresolver.h>
+#include <mitsuba/render/testcase.h>
+#include <boost/bind.hpp>
+#include "../bsdfs/microfacet.h"
+
+/* Statistical significance level of the test. Set to
+   1/4 percent by default -- we want there to be strong
+   evidence of an implementaiton error before failing
+   a test case */
+#define SIGNIFICANCE_LEVEL 0.0025f
+
+/* Relative bound on what is still accepted as roundoff
+   error -- be quite tolerant */
+#if defined(SINGLE_PRECISION)
+	#define ERROR_REQ 1e-2f
+#else
+	#define ERROR_REQ 1e-5
+#endif
+
+MTS_NAMESPACE_BEGIN
+
+class TestChiSquare : public TestCase {
+public:
+	MTS_BEGIN_TESTCASE()
+	MTS_DECLARE_TEST(test01_Microfacet)
+	MTS_DECLARE_TEST(test02_MicrofacetVisible)
+	MTS_END_TESTCASE()
+
+	class MicrofacetAdapter {
+	public:
+		MicrofacetAdapter(Sampler *sampler, const MicrofacetDistribution &distr, const Vector &wi = Vector(0.0f)) : m_sampler(sampler), m_distr(distr), m_wi(wi) { }
+
+		boost::tuple<Vector, Float, EMeasure> generateSample() {
+			Float pdf;
+
+			if (m_wi.lengthSquared() == 0) {
+				Normal m = m_distr.sampleAll(m_sampler->next2D(), pdf);
+				Float pdf_ref = m_distr.pdfAll(m);
+
+				SAssert(std::isfinite(pdf) && pdf > 0);
+				SAssert(std::isfinite(pdf_ref) && pdf_ref > 0);
+				SAssert(std::isfinite(m.x) && std::isfinite(m.y) && std::isfinite(m.z));
+				SAssert(std::abs(m.length() - 1) < 1e-4f);
+				SAssert(std::abs((pdf-pdf_ref)/pdf_ref) < 1e-4f);
+				return boost::make_tuple(m, 1.0f, ESolidAngle);
+			} else {
+				Normal m = m_distr.sampleVisible(m_wi, m_sampler->next2D());
+				SAssert(std::isfinite(m.x) && std::isfinite(m.y) && std::isfinite(m.z));
+				SAssert(std::abs(m.length() - 1) < 1e-4f);
+				return boost::make_tuple(m, 1.0f, ESolidAngle);
+			}
+		}
+
+		Float pdf(const Vector &d, EMeasure measure) const {
+			if (measure != ESolidAngle)
+				return 0.0f;
+
+			Float pdf = m_wi.lengthSquared() == 0 ? m_distr.pdfAll(d)
+				: m_distr.pdfVisible(m_wi, d);
+			SAssert(std::isfinite(pdf) && pdf >= 0);
+
+			return pdf;
+		}
+
+	private:
+		ref<Sampler> m_sampler;
+		MicrofacetDistribution m_distr;
+		Vector m_wi;
+	};
+
+	void test01_Microfacet() {
+		int thetaBins = 20;
+		std::vector<MicrofacetDistribution> distrs;
+
+		distrs.push_back(MicrofacetDistribution(MicrofacetDistribution::EBeckmann, 0.5f));
+		distrs.push_back(MicrofacetDistribution(MicrofacetDistribution::EBeckmann, 0.5f, 0.3f));
+		distrs.push_back(MicrofacetDistribution(MicrofacetDistribution::EGGX, 0.5f));
+		distrs.push_back(MicrofacetDistribution(MicrofacetDistribution::EGGX, 0.5f, 0.3f));
+		distrs.push_back(MicrofacetDistribution(MicrofacetDistribution::EPhong, 0.5f));
+		distrs.push_back(MicrofacetDistribution(MicrofacetDistribution::EPhong, 0.5f, 0.3f));
+
+
+		ref<Sampler> sampler = static_cast<Sampler *> (PluginManager::getInstance()->
+				createObject(MTS_CLASS(Sampler), Properties("independent")));
+		ref<ChiSquare> chiSqr = new ChiSquare(thetaBins, 2*thetaBins, (int) distrs.size());
+		chiSqr->setLogLevel(EDebug);
+
+		for (size_t i=0; i<distrs.size(); ++i) {
+			Log(EInfo, "Testing %s", distrs[i].toString().c_str());
+			// Initialize the tables used by the chi-square test
+			MicrofacetAdapter adapter(sampler, distrs[i]);
+			chiSqr->fill(
+				boost::bind(&MicrofacetAdapter::generateSample, &adapter),
+				boost::bind(&MicrofacetAdapter::pdf, &adapter, _1, _2)
+			);
+
+			// (the following assumes that the distribution has 1 parameter, e.g. exponent value)
+			ChiSquare::ETestResult result = chiSqr->runTest(SIGNIFICANCE_LEVEL);
+			if (result == ChiSquare::EReject) {
+				std::string filename = formatString("failure_%i.m", (int) i);
+				chiSqr->dumpTables(filename);
+				failAndContinue(formatString("Uh oh, the chi-square test indicates a potential "
+					"issue. Dumped the contingency tables to '%s' for user analysis",
+					filename.c_str()));
+			} else {
+				//chiSqr->dumpTables(formatString("success_%i.m", (int) i));
+				succeed();
+			}
+		}
+	}
+
+	void test02_MicrofacetVisible() {
+		int thetaBins = 10;
+		std::vector<std::pair<MicrofacetDistribution, Vector> > distrs;
+
+		ref<Sampler> sampler = static_cast<Sampler *> (PluginManager::getInstance()->
+				createObject(MTS_CLASS(Sampler), Properties("independent")));
+		for (int i=0; i<10; ++i) {
+			Vector wi = warp::squareToUniformHemisphere(sampler->next2D());
+			distrs.push_back(std::make_pair(MicrofacetDistribution(MicrofacetDistribution::EBeckmann, 0.3f), wi));
+			distrs.push_back(std::make_pair(MicrofacetDistribution(MicrofacetDistribution::EBeckmann, 0.5f, 0.3f), wi));
+			distrs.push_back(std::make_pair(MicrofacetDistribution(MicrofacetDistribution::EGGX, 0.1f), wi));
+			distrs.push_back(std::make_pair(MicrofacetDistribution(MicrofacetDistribution::EGGX, 0.2f, 0.3f), wi));
+		}
+
+		ref<ChiSquare> chiSqr = new ChiSquare(thetaBins, 2*thetaBins, (int) distrs.size());
+		chiSqr->setLogLevel(EDebug);
+
+		for (size_t i=0; i<distrs.size(); ++i) {
+			Log(EInfo, "Testing %s (wi=%s)", distrs[i].first.toString().c_str(), distrs[i].second.toString().c_str());
+			// Initialize the tables used by the chi-square test
+			MicrofacetAdapter adapter(sampler, distrs[i].first, distrs[i].second);
+			chiSqr->fill(
+				boost::bind(&MicrofacetAdapter::generateSample, &adapter),
+				boost::bind(&MicrofacetAdapter::pdf, &adapter, _1, _2)
+			);
+
+			// (the following assumes that the distribution has 1 parameter, e.g. exponent value)
+			ChiSquare::ETestResult result = chiSqr->runTest(SIGNIFICANCE_LEVEL);
+			if (result == ChiSquare::EReject) {
+				std::string filename = formatString("failure_%i.m", (int) i);
+				chiSqr->dumpTables(filename);
+				failAndContinue(formatString("Uh oh, the chi-square test indicates a potential "
+					"issue. Dumped the contingency tables to '%s' for user analysis",
+					filename.c_str()));
+			} else {
+				//chiSqr->dumpTables(formatString("success_%i.m", (int) i));
+				succeed();
+			}
+		}
+	}
+};
+
+MTS_EXPORT_TESTCASE(TestChiSquare, "Chi-square test for microfacet sampling")
+MTS_NAMESPACE_END