diff --git a/include/mitsuba/core/chisquare.h b/include/mitsuba/core/chisquare.h
new file mode 100644
index 00000000..b506abff
--- /dev/null
+++ b/include/mitsuba/core/chisquare.h
@@ -0,0 +1,278 @@
+/*
+    This file is part of Mitsuba, a physically based rendering system.
+
+    Copyright (c) 2007-2010 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/>.
+*/
+
+#if !defined(__CHI_SQUARE_TEST_H)
+#define __CHI_SQUARE_TEST_H
+
+#include <mitsuba/core/fresolver.h>
+#include <mitsuba/core/quad.h>
+#include <mitsuba/core/timer.h>
+#include <mitsuba/core/random.h>
+#include <boost/bind.hpp>
+#include <boost/math/distributions/chi_squared.hpp>
+
+MTS_NAMESPACE_BEGIN
+
+/**
+ * \brief Chi-square goodness-of-fit test on the sphere
+ *
+ * This class performs a chi-square goodness-of-fit test of the null hypothesis
+ * that a specified sampling procedure produces samples that are distributed
+ * according to a supplied density function. This is very useful to verify BRDF
+ * and phase function sampling codes for their correctness.
+ *
+ * This implementation works by generating a large batch of samples, which are
+ * then accumulated into rectangular bins in spherical coordinates. To obtain
+ * reference bin counts, the provided density function is numerically
+ * integrated over the area of each bin. Comparing the actual and reference
+ * bin counts yields the desired test statistic.
+ *
+ * Given a BSDF/phase function with the following interface
+ * <code>
+ * class MyDistribution {
+ *     Vector sample(const Point2 &uniformSample) const;
+ *     Float pdf(const Vector &direction) const;
+ * };
+ * </code>
+ *
+ * the code in this class might be used as follows
+ * 
+ * <code>
+ * MyDistribution myDistrInstance;
+ * ChiSquareTest chiSqr;
+ *
+ * // Initialize the tables used by the chi-square test
+ * chiSqr.fill(
+ *    boost::bind(&MyDistribution::sample, myDistrInstance, _1),
+ *    boost::bind(&MyDistribution::pdf, myDistrInstance, _1)
+ * );
+ *
+ * // Optional: dump the tables to a MATLAB file for external analysis
+ * chiSqr.dumpTables("debug.m");
+ *
+ * // (the folowing assumes that the distribution has 1 parameter, e.g. exponent value)
+ * if (!chiSqr.runTest(1))
+ *    Log(EError, "Uh oh -- test failed, the implementation is probably incorrect!");
+ * </code>
+ */
+class ChiSquareTest {
+public:
+	/**
+	 * \brief Create a new Chi-square test instance with the given
+	 * resolution and sample count
+	 *
+	 * \param thetaBins
+	 *    Number of bins wrt. latitude. The default is 10
+	 *
+	 * \param phiBins
+	 *    Number of bins wrt. azimuth. The default is to use
+	 *    twice the number of \c thetaBins
+	 *
+	 * \param sampleCount
+	 *    Number of samples to be used when computing the bin
+	 *    values. The default is \c thetaBins*phiBins*1000
+	 */
+	ChiSquareTest(int thetaBins = 10, int phiBins = 0, size_t sampleCount = 0)
+	    	: m_thetaBins(thetaBins), m_phiBins(phiBins), m_sampleCount(sampleCount) {
+		if (m_phiBins == 0)
+			m_phiBins = 2*m_thetaBins;
+		if (m_sampleCount == 0)
+			m_sampleCount = m_thetaBins * m_phiBins * 1000;
+		m_table = new uint32_t[m_thetaBins*m_phiBins];
+		m_refTable = new Float[m_thetaBins*m_phiBins];
+	}
+
+	/// Release all memory
+	~ChiSquareTest() {
+		delete[] m_table;
+		delete[] m_refTable;
+	}
+
+	/**
+	 * \brief Fill the actual and reference bin counts
+	 *
+	 * Please see the class documentation for a description
+	 * on how to invoke this function
+	 */
+	void fill(
+		const boost::function<Vector (const Point2 &)> &sampleFn,
+		const boost::function<Float (const Vector &)> &pdfFn) {
+		ref<Random> random = new Random();
+		memset(m_table, 0, m_thetaBins*m_phiBins*sizeof(uint32_t));
+
+		SLog(EInfo, "Accumulating " SIZE_T_FMT " samples into a %ix%i"
+				" contingency table", m_sampleCount, m_thetaBins, m_phiBins);
+		Point2 factor(m_thetaBins / M_PI, m_phiBins / (2*M_PI));
+
+		ref<Timer> timer = new Timer();
+		for (size_t i=0; i<m_sampleCount; ++i) {
+			Point2 sample(random->nextFloat(), random->nextFloat());
+			Point2 sphCoords = toSphericalCoordinates(sampleFn(sample));
+
+			int thetaBin = std::min(std::max(0,
+				floorToInt(sphCoords.x * factor.x)), m_thetaBins-1);
+			int phiBin = std::min(std::max(0,
+				floorToInt(sphCoords.y * factor.y)), m_phiBins-1);
+
+			++m_table[thetaBin * m_phiBins + phiBin];
+		}
+		SLog(EInfo, "Done, took %i ms.", timer->getMilliseconds());
+		factor = Point2(M_PI / m_thetaBins, (2*M_PI) / m_phiBins);
+
+		SLog(EInfo, "Integrating reference contingency table");
+		timer->reset();
+		Float min[2], max[2];
+		size_t idx = 0;
+
+		NDIntegrator integrator(1, 2, 50000, 0, 1e-6f);
+		Float maxError = 0, integral = 0;
+		for (int i=0; i<m_thetaBins; ++i) {
+			min[0] = i * factor.x;
+			max[0] = (i+1) * factor.x;
+			for (int j=0; j<m_phiBins; ++j) {
+				min[1] = j * factor.y;
+				max[1] = (j+1) * factor.y;
+				Float result, error;
+				size_t evals;
+
+				integrator.integrateVectorized(
+					boost::bind(&ChiSquareTest::integrand, pdfFn, _1, _2, _3),
+					min, max, &result, &error, evals
+				);
+
+				integral += result;
+				m_refTable[idx++] = result * m_sampleCount;
+				maxError = std::max(maxError, error);
+			}
+		}
+		SLog(EInfo, "Done, took %i ms (max error = %f, integral=%f).", 
+				timer->getMilliseconds(), maxError, integral);
+	}
+
+	/**
+	 * \brief Dump the bin counts to a file using MATLAB format
+	 */
+	void dumpTables(const fs::path &filename) {
+		fs::ofstream out(filename);
+		out << "tbl_counts = [ ";
+		for (int i=0; i<m_thetaBins; ++i) {
+			for (int j=0; j<m_phiBins; ++j) {
+				out << m_table[i*m_phiBins+j];
+				if (j+1 < m_phiBins)
+					out << ", ";
+			}
+			if (i+1 < m_thetaBins)
+				out << "; ";
+		}
+		out << " ];" << endl
+			<< "tbl_ref = [ ";
+		for (int i=0; i<m_thetaBins; ++i) {
+			for (int j=0; j<m_phiBins; ++j) {
+				out << m_refTable[i*m_phiBins+j];
+				if (j+1 < m_phiBins)
+					out << ", ";
+			}
+			if (i+1 < m_thetaBins)
+				out << "; ";
+		}
+		out << " ];" << endl;
+		out.close();
+	}
+
+	/**
+	 * \brief Perform the actual chi-square test
+	 *
+	 * \param distParams
+	 *     Number of parameters of the distribution in question.
+	 *     Anything such as lobe width, indices of refraction, etc.
+	 *     should be counted.
+	 *
+	 * \param pvalThresh
+	 *     The implementation will reject the null hypothesis
+	 *     when the computed p-value lies below this parameter
+	 *     (default: 0.01f)
+	 *
+	 * \return \c false if the null hypothesis was rejected
+	 *     and \c true otherwise.
+	 */
+	bool runTest(int distParams, Float pvalThresh = 0.01f) {
+		/* Compute the chi-square statistic */
+		Float chsq = 0.0f;
+		Float pooledCounts = 0, pooledRef = 0;
+		int pooledCells = 0;
+		int df = 0;
+
+		for (int i=0; i<m_thetaBins*m_phiBins; ++i) {
+			if (m_refTable[i] < 5) {
+				pooledCounts += m_table[i];
+				pooledRef += m_refTable[i];
+				++pooledCells;
+			} else {
+				Float diff = m_table[i]-m_refTable[i];
+				chsq += (diff*diff) / m_refTable[i];
+				++df;
+			}
+		}
+
+		if (pooledCells > 0) {
+			SLog(EInfo, "Pooled %i cells with an expected "
+				"number of < 5 entries!", pooledCells);
+			if (pooledRef < 5) {
+				SLog(EWarn, "Even after pooling, the expected "
+					"number of entries is < 5 (%f), expect badness!",
+					pooledRef);
+			}
+			Float diff = pooledCounts - pooledRef;
+			chsq += (diff*diff) / pooledRef;
+			++df;
+		}
+
+		df -= distParams + 1;
+		SLog(EInfo, "Chi-square statistic = %e (df=%i)", chsq, df);
+		boost::math::chi_squared chSqDist(df);
+		/* Probability of obtaining a test statistic at least
+		   as extreme as the one observed under the assumption
+		   that the distributions match */
+		Float pval = 1 - boost::math::cdf(chSqDist, chsq);
+		SLog(EInfo, "P-value = %e", pval);
+
+		if (pval < pvalThresh) {
+			SLog(EWarn, "Rejecting the null hypothesis");
+			return false;
+		}
+		return true;
+
+	}
+protected:
+	static void integrand(
+		const boost::function<Float (const Vector &)> &pdfFn,
+			size_t nPts, const Float *in, Float *out) {
+		#pragma omp parallel for
+		for (size_t i=0; i<nPts; ++i)
+			out[i] = pdfFn(sphericalDirection(in[2*i], in[2*i+1])) * std::sin(in[2*i]);
+	}
+private:
+	int m_thetaBins, m_phiBins;
+	size_t m_sampleCount;
+	uint32_t *m_table;
+	Float *m_refTable;
+};
+
+MTS_NAMESPACE_END
+
+#endif /* __CHI_SQUARE_TEST_H */
diff --git a/include/mitsuba/core/stl.h b/include/mitsuba/core/stl.h
index b381f1d1..cbaaeb5a 100644
--- a/include/mitsuba/core/stl.h
+++ b/include/mitsuba/core/stl.h
@@ -88,76 +88,8 @@ namespace std {
 };
 using std::select2nd;
 using std::compose1;
-
 #endif	
 
-/* Forward declarations */
-MTS_NAMESPACE_BEGIN
-extern MTS_EXPORT_CORE void * __restrict allocAligned(size_t size);
-extern MTS_EXPORT_CORE void freeAligned(void *ptr);
-
-/**
- * \brief Aligned memory allocator for use with SSE2-based code
- * \headerfile mitsuba/core/stl.h mitsuba/mitsuba.h
- * 
- * Basic implementaiton, which forwards all calls to \ref allocAligned.
- */
-template <typename T> class aligned_allocator {
-public:
-	typedef size_t    size_type;
-	typedef ptrdiff_t difference_type;
-	typedef T*        pointer;
-	typedef const T*  const_pointer;
-	typedef T&        reference;
-	typedef const T&  const_reference;
-	typedef T         value_type;
-
-	/// \cond
-	template <class U> struct rebind {
-		typedef aligned_allocator<U> other;
-	};
-	/// \endcond
-
-	pointer address (reference value) const {
-		return &value;
-	}
-
-	const_pointer address (const_reference value) const {
-		return &value;
-	}
-	
-	aligned_allocator() throw() { }	
-
-	aligned_allocator(const aligned_allocator&) throw() { }
-
-	template <class U> aligned_allocator (const aligned_allocator<U>&) throw()  { }
-
-	~aligned_allocator() throw() { }
-
-	size_type max_size () const throw() {
-		return INT_MAX;
-	}
-
-
-	T* __restrict allocate (size_type num, const_pointer *hint = 0) {
-		return (T *) mitsuba::allocAligned(num*sizeof(T));
-	}
-
-	void construct (pointer p, const T& value) {
-		*p=value;
-	}
-
-	void destroy (pointer p) {
-		p->~T();
-	};
-	
-	void deallocate (pointer p, size_type num) {
-		freeAligned(p);
-	}
-};
-
-MTS_NAMESPACE_END
-
 #if defined(WIN32)
 inline bool mts_isnan(float f) {
 	int classification = ::_fpclass(f);
diff --git a/src/libcore/util.cpp b/src/libcore/util.cpp
index 9b1e3b3f..088b7256 100644
--- a/src/libcore/util.cpp
+++ b/src/libcore/util.cpp
@@ -616,10 +616,13 @@ Point2 squareToTriangle(const Point2 &sample) {
 }
 
 Point2 toSphericalCoordinates(const Vector &v) {
-	return Point2(
+	Point2 result(
 		std::acos(v.z),
 		std::atan2(v.y, v.x)
 	);
+	if (result.y < 0)
+		result.y += 2*M_PI;
+	return result;
 }
 
 Point2 squareToDiskConcentric(const Point2 &sample) {
diff --git a/src/medium/heterogeneous-flake.cpp b/src/medium/heterogeneous-flake.cpp
deleted file mode 100644
index e7c72c06..00000000
--- a/src/medium/heterogeneous-flake.cpp
+++ /dev/null
@@ -1,736 +0,0 @@
-/*
-    This file is part of Mitsuba, a physically based rendering system.
-
-    Copyright (c) 2007-2010 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/render/scene.h>
-#include <mitsuba/render/volume.h>
-#include <mitsuba/core/fstream.h>
-#include <mitsuba/core/brent.h>
-#include <mitsuba/core/shvector4d.h>
-#include <mitsuba/core/statistics.h>
-#include <fstream>
-
-/*
- * When inverting an integral of the form f(t):=\int_0^t [...] dt
- * using composite Simpson's rule quadrature, the following constant
- * specified how many times the step size will be reduced when we
- * have almost reached the desired value.
- */
-#define NUM_REFINES 8 
-
-/**
- * Allow to stop integrating densities when the resulting segment
- * has a throughput of less than 'Epsilon'
- */
-#define HET_EARLY_EXIT 1
-
-/// Header identification record for phase function caches
-#define PHASE_FUNCTION_HEADER 0x40044004
-
-#define USE_REJECTION_SAMPLING 1
-#define USE_COSTHETA_DISTR 1
-
-MTS_NAMESPACE_BEGIN
-
-static StatsCounter totalSamples("Micro-flake model", "Sample generations");
-static StatsCounter avgSampleIterations("Micro-flake model", "Average rejection sampling iterations", EAverage);
-static StatsCounter brentSolves("Micro-flake model", "Brent solver calls");
-static StatsCounter avgBrentFunEvals("Micro-flake model", "Average Brent solver function evaluations", EAverage);
-
-struct AbsCos {
-	Float operator()(const Vector &w) const { return std::abs(w.z); }
-};
-
-#define USE_COSTHETA_DISTR 1
-
-enum EMode {
-	ESinP = 0,
-	ECosP
-};
-
-struct FlakeDistr {
-	Float exp;
-	EMode mode;
-
-	FlakeDistr() { }
-	FlakeDistr(EMode mode, Float exp) : exp(exp), mode(mode) { }
-
-	inline Float operator()(const Vector &w) const {
-		switch (mode) {
-			case ESinP:
-				return std::pow(1-std::pow(w.z, 2), exp/2.0f);
-			case ECosP:
-				return std::pow(std::abs(w.z), exp);
-			default:
-				SLog(EError, "Invalid flake distribution!");
-				return 0.0f;
-		}
-	}
-};
-
-template <typename Distr> struct FlakePhaseFunctor {
-	Distr d;
-
-	FlakePhaseFunctor(Distr d) : d(d) { }
-
-	inline Float operator()(const Vector &wi, const Vector &wo) const {
-		Vector H = normalize(wi + wo);
-		return d(H) + d(-H);
-	}
-};
-
-class HeterogeneousFlakeMedium;
-
-
-class FlakePhaseFunction : public PhaseFunction {
-public:
-	FlakePhaseFunction(const SHVector &sigmaT, int samplingRecursions, const SHVector4D *phaseExpansion, 
-		Float exponent, Float normalization, EMode mode, HeterogeneousFlakeMedium *medium) : PhaseFunction(Properties()), 
-		m_sigmaT(sigmaT), m_samplingRecursions(samplingRecursions), m_distr(mode, exponent), m_phaseExpansion(phaseExpansion), 
-		m_exponent(exponent), m_normalization(normalization), m_mode(mode), m_medium(medium) {
-		initialize();
-	}
-
-	FlakePhaseFunction(Stream *stream, InstanceManager *manager) :
-			PhaseFunction(stream, manager) {
-		/* Not implemented */
-	}
-
-	void serialize(Stream *stream, InstanceManager *manager) const {
-		PhaseFunction::serialize(stream, manager);
-		/* Not implemented */
-	}
-
-	void initialize() {
-		m_shSampler = new SHSampler(m_phaseExpansion->getBands(), m_samplingRecursions);
-		Log(EInfo, "Constructing a SH sampler: %s", m_shSampler->toString().c_str());
-	}
-
-	virtual ~FlakePhaseFunction() { }
-
-	/* Implementations are below */
-	Spectrum f(const PhaseFunctionQueryRecord &pRec) const;
-
-	Spectrum sample(PhaseFunctionQueryRecord &pRec,
-		Float &pdf, Sampler *sampler) const;
-	
-	Spectrum sample(PhaseFunctionQueryRecord &pRec,
-		Sampler *sampler) const;
-
-	Float pdf(const PhaseFunctionQueryRecord &pRec) const;
-
-	std::string toString() const {
-		return "FlakePhaseFunction[]";
-	}
-
-	MTS_DECLARE_CLASS()
-private:
-	SHVector m_sigmaT;
-	ref<SHSampler> m_shSampler;
-	int m_samplingRecursions;
-	FlakeDistr m_distr;
-	const SHVector4D *m_phaseExpansion;
-	mutable PrimitiveThreadLocal<SHVector> m_temp;
-	Float m_exponent;
-	Float m_normalization;
-	EMode m_mode;
-	const HeterogeneousFlakeMedium *m_medium; 
-};
-
-class HeterogeneousFlakeMedium : public Medium {
-public:
-	HeterogeneousFlakeMedium(const Properties &props) 
-		: Medium(props) {
-		m_stepSize = props.getFloat("stepSize", 0);
-
-		/* Flake model configuration */
-		int bands = props.getInteger("bands", 7);
-		int phaseBands = props.getInteger("phaseBands", bands);
-		Assert(bands > 0 && phaseBands > 0);
-
-		m_samplingRecursions = props.getInteger("samplingRecursions", 10);
-		std::string mode = props.getString("mode", "sinp");
-
-		if (mode == "sinp") {
-			m_mode = ESinP;
-		} else if (mode == "cosp") {
-			m_mode = ECosP;
-		} else {
-			Log(EError, "Unknown mode: %s!", mode.c_str());
-		}	
-
-		const int numSamples = 80;
-		ref<FileStream> stream;
-
-		SHVector absCos = SHVector(bands);
-		absCos.project(AbsCos(), numSamples);
-		Float absCosError = absCos.l2Error(AbsCos(), numSamples);
-		
-		m_exponent = props.getFloat("exponent", 4);
-		FlakeDistr distrFunc(m_mode, m_exponent);
-
-		m_D = SHVector(bands);
-		m_D.project(distrFunc, numSamples);
-		m_normalization = 1 / (m_D(0, 0) * 2 * std::sqrt(M_PI));
-		Float dError = m_D.l2Error(distrFunc, numSamples);
-		m_D.normalize();
-
-		Log(EInfo, "=======  Flake medium properties =======");
-		Log(EInfo, "       distribution type = %s", mode.c_str());
-		Log(EInfo, "   distribution exponent = %f", m_exponent);
-		Log(EInfo, "    SH bands (D, sigmaT) = %i", bands);
-		Log(EInfo, "    SH bands (phase fct) = %i", phaseBands);
-		Log(EInfo, "  |cos| projection error = %f", absCosError);
-		Log(EInfo, "     D  projection error = %f", dError);
-		Log(EInfo, "     sampling recursions = %i", m_samplingRecursions);
-
-		bool computePhaseProjection = true;
-		if (fs::exists("flake-phase.dat")) {
-			/* Avoid recomputing this every time */
-			stream = new FileStream("flake-phase.dat", FileStream::EReadOnly);
-			unsigned int header = stream->readUInt();
-			int nBands = stream->readInt();
-			int mode = stream->readInt();
-			Float exponent = stream->readFloat();
-			if (header == PHASE_FUNCTION_HEADER && nBands == phaseBands 
-				&& exponent == m_exponent && mode == m_mode) {
-				m_phaseExpansion = SHVector4D(stream);
-				computePhaseProjection = false;
-			}
-			stream->close();
-		}
-
-		if (computePhaseProjection) {
-			int res = 50;
-			FlakePhaseFunctor<FlakeDistr> phaseFunctor(distrFunc);
-			m_phaseExpansion = SHVector4D(res, 2*res, phaseBands);
-			Log(EDebug, "Phase function parameterization: %s", m_phaseExpansion.toString().c_str());
-			m_phaseExpansion.project(phaseFunctor, numSamples);
-			m_phaseExpansion.normalize();
-			stream = new FileStream("flake-phase.dat", FileStream::ETruncReadWrite);
-			stream->writeInt(PHASE_FUNCTION_HEADER);
-			stream->writeInt(phaseBands);
-			stream->writeInt(m_mode);
-			stream->writeFloat(m_exponent);
-			m_phaseExpansion.serialize(stream);
-			stream->close();
-		}
-
-		for (int i=0; i<phaseBands; ++i)
-			Log(EInfo, "Energy in phase function band %i: %f", i, m_phaseExpansion.energy(i));
-
-		/* Compute sigma_t except for the a*rho factor */
-		m_sigmaT = SHVector(m_D);
-		m_sigmaT.convolve(absCos);
-		Assert(m_sigmaT.isAzimuthallyInvariant());
-		m_kdTree = new ShapeKDTree();
-	}
-	
-	/* Unserialize from a binary data stream */
-	HeterogeneousFlakeMedium(Stream *stream, InstanceManager *manager) 
-		: Medium(stream, manager) {
-		m_kdTree = new ShapeKDTree();
-		size_t shapeCount = stream->readUInt();
-		for (size_t i=0; i<shapeCount; ++i) 
-			addChild("", static_cast<Shape *>(manager->getInstance(stream)));
-		m_densities = static_cast<VolumeDataSource *>(manager->getInstance(stream));
-		m_albedo = static_cast<VolumeDataSource *>(manager->getInstance(stream));
-		m_orientations = static_cast<VolumeDataSource *>(manager->getInstance(stream));
-		m_stepSize = stream->readFloat();
-		m_exponent = stream->readFloat();
-		m_normalization = stream->readFloat();
-		m_mode = (EMode) stream->readInt();
-		m_D = SHVector(stream);
-		m_sigmaT = SHVector(stream);
-		m_phaseExpansion = SHVector4D(stream);
-		m_samplingRecursions = stream->readInt();
-		configure();
-	}
-
-	virtual ~HeterogeneousFlakeMedium() {
-		for (size_t i=0; i<m_shapes.size(); ++i)
-			m_shapes[i]->decRef();
-	}
-
-	/* Serialize the volume to a binary data stream */
-	void serialize(Stream *stream, InstanceManager *manager) const {
-		Medium::serialize(stream, manager);
-		stream->writeUInt((uint32_t) m_shapes.size());
-		for (size_t i=0; i<m_shapes.size(); ++i)
-			manager->serialize(stream, m_shapes[i]);
-		manager->serialize(stream, m_densities.get());
-		manager->serialize(stream, m_albedo.get());
-		manager->serialize(stream, m_orientations.get());
-		stream->writeFloat(m_stepSize);
-		stream->writeFloat(m_exponent);
-		stream->writeFloat(m_normalization);
-		stream->writeInt(m_mode);
-		m_D.serialize(stream);
-		m_sigmaT.serialize(stream);
-		m_phaseExpansion.serialize(stream);
-		stream->writeInt(m_samplingRecursions);
-	}
-
-	void configure() {
-		Medium::configure();
-		if (m_densities.get() == NULL)
-			Log(EError, "No densities specified!");
-		if (m_albedo.get() == NULL)
-			Log(EError, "No albedo specified!");
-		if (m_orientations.get() == NULL)
-			Log(EError, "No orientations specified!");
-
-		if (m_stepSize == 0) {
-			m_stepSize = std::min(std::min(
-				m_densities->getStepSize(),
-				m_albedo->getStepSize()),
-				m_orientations->getStepSize());
-
-			if (m_stepSize == std::numeric_limits<Float>::infinity()) 
-				Log(EError, "Unable to infer step size, please specify!");
-		}
-
-		m_phaseFunction = new FlakePhaseFunction(m_sigmaT, m_samplingRecursions,
-			&m_phaseExpansion, m_exponent, m_normalization, m_mode, this);
-
-		if (m_shapes.size() > 0) {
-			m_kdTree->build();
-			m_aabb = m_kdTree->getAABB();
-		} else {
-			m_aabb = m_densities->getAABB();
-		}
-	}
-
-	void setParent(ConfigurableObject *parent) {
-		if (parent->getClass()->derivesFrom(MTS_CLASS(Shape)))
-			Log(EError, "Medium cannot be a parent of a shape");
-	}
-
-
-	void addChild(const std::string &name, ConfigurableObject *child) {
-		if (child->getClass()->derivesFrom(MTS_CLASS(Shape))) {
-			Shape *shape = static_cast<Shape *>(child);
-			if (shape->isCompound()) {
-				int ctr = 0;
-				while (true) {
-					ref<Shape> childShape = shape->getElement(ctr++);
-					if (!childShape)
-						break;
-					addChild("", childShape);
-				}
-			} else {
-				m_kdTree->addShape(shape);
-				shape->incRef();
-				m_shapes.push_back(shape);
-			}
-		} else if (child->getClass()->derivesFrom(MTS_CLASS(VolumeDataSource))) {
-			VolumeDataSource *volume = static_cast<VolumeDataSource *>(child);
-
-			if (name == "albedo") {
-				Assert(volume->supportsSpectrumLookups());
-				m_albedo = volume;
-			} else if (name == "density") {
-				Assert(volume->supportsFloatLookups());
-				m_densities = volume;
-			} else if (name == "orientations") {
-				Assert(volume->supportsVectorLookups());
-				m_orientations = volume;
-			}
-		} else {
-			Medium::addChild(name, child);
-		}
-	}
-
-	/**
-	 * \brief Intersect a ray with the stencil and 
-	 * determine whether it started on the inside
-	 */
-	inline bool rayIntersect(const Ray &ray, Float &t, bool &inside) const {
-		if (m_shapes.size() != 0) {
-			Intersection its;
-			if (!m_kdTree->rayIntersect(ray, its)) 
-				return false;
-			inside = dot(its.geoFrame.n, ray.d) > 0;
-			t = its.t;
-		} else {
-			Float mint, maxt;
-			if (!m_aabb.rayIntersect(ray, mint, maxt))
-				return false;
-			if (mint <= ray.mint && maxt <= ray.mint)
-				return false;
-			inside = mint <= 0 && maxt > 0;
-			t = (mint <= 0) ? maxt : mint;
-		}
-		return true;
-	}
-
-	inline Vector lookupOrientation(const Point &p) const {
-		Vector orientation = m_orientations->lookupVector(p);
-		Float lengthSqr = orientation.lengthSquared();
-		if (lengthSqr != 0)
-			return orientation / std::sqrt(lengthSqr);
-		else
-			return Vector(0.0f);
-	}
-
-	inline Float sigmaT(const Point &p, const Vector &d) const {
-		/**
-		 * Double the densities so that a uniform flake distribution
-		 * reproduces isotropic scattering
-		 */
-		Float density = 2 * m_densities->lookupFloat(p) * m_densityMultiplier;
-		if (density == 0)
-			return 0.0f;
-		Vector orientation = lookupOrientation(p);
-		if (orientation == Vector(0.0f))
-			return 0.0f;
-
-		Vector localD = Frame(orientation).toLocal(d);
-
-		return density * m_sigmaT.evalAzimuthallyInvariant(localD);
-	}
-
-	/// Integrate densities using the composite Simpson's rule
-	inline Float integrateDensity(Ray ray, Float length) const {
-		if (length == 0)
-			return 0.0f;
-		int nParts = (int) std::ceil(length/m_stepSize);
-		nParts += nParts % 2;
-		const Float stepSize = length/nParts;
-		const Vector increment = ray.d * stepSize;
-
-		Float m = 4;
-		Float accumulatedDensity = 
-			sigmaT(ray.o, ray.d) + sigmaT(ray(length), ray.d);
-#ifdef HET_EARLY_EXIT
-		const Float stopAfterDensity = -std::log(Epsilon);
-		const Float stopValue = stopAfterDensity*3/stepSize;
-#endif
-
-		ray.o += increment;
-		for (int i=1; i<nParts; ++i) {
-			Float value = sigmaT(ray.o, ray.d);
-			accumulatedDensity += value*m;
-			Point tmp = ray.o;
-			ray.o += increment;
-
-			if (ray.o == tmp) {
-				Log(EWarn, "integrateDensity(): failed to make forward progress -- "
-						"round-off error issues? The step size was %.18f", stepSize);
-				break;
-			}
-
-#ifdef HET_EARLY_EXIT
-			if (accumulatedDensity > stopValue) // Stop early
-				return std::numeric_limits<Float>::infinity();
-#endif
-			m = 6 - m;
-		}
-		accumulatedDensity *= stepSize/3;
-
-		return accumulatedDensity;
-	}
-
-	/**
-	 * \brief Integrate the density covered by the specified ray
-	 * segment, while clipping the volume to the underlying
-	 * stencil volume.
-	 */
-	inline Float integrateDensity(const Ray &r) const {
-		Ray ray(r(r.mint), r.d,
-			0, std::numeric_limits<Float>::infinity(), 0.0f);
-		Float integral = 0, remaining = r.maxt - r.mint;
-		int iterations = 0;
-		bool inside = false;
-		Float t = 0;
-
-		while (remaining > 0 && rayIntersect(ray, t, inside)) {
-			if (inside) 
-				integral += integrateDensity(ray, std::min(remaining, t));
-
-			remaining -= t;
-			ray.o = ray(t);
-			ray.mint = Epsilon;
-
-			if (++iterations > 100) {
-				/// Just a precaution..
-				Log(EWarn, "integrateDensity(): round-off error issues?");
-				break;
-			}
-		}
-		return integral;
-	}
-	
-	/**
-	 * Attempts to solve the following 1D integral equation for 't'
-	 * \int_0^t density(ray.o + x * ray.d) * dx + accumulatedDensity == desiredDensity.
-	 * When no solution can be found in [0, maxDist] the function returns
-	 * false. For convenience, the function returns the current values of sigmaT 
-	 * and the albedo, as well as the 3D position 'ray.o+t*ray.d' upon
-	 * success.
-	 */
-	bool invertDensityIntegral(Ray ray, Float maxDist, 
-			Float &accumulatedDensity, Float desiredDensity,
-			Float &currentSigmaT, Spectrum &currentAlbedo, 
-			Point &currentPoint) const {
-		if (maxDist == 0)
-			return std::numeric_limits<Float>::infinity();
-		int nParts = (int) std::ceil(maxDist/m_stepSize);
-		Float stepSize = maxDist/nParts;
-		Vector fullIncrement = ray.d * stepSize,
-			   halfIncrement = fullIncrement * .5f;
-		Float t = 0;
-		int numRefines = 0;
-		bool success = false;
-		Float node1 = sigmaT(ray.o, ray.d);
-
-		while (t <= maxDist) {
-			Float node2 = sigmaT(ray.o + halfIncrement, ray.d);
-			Float node3 = sigmaT(ray.o + fullIncrement, ray.d);
-			const Float newDensity = accumulatedDensity + 
-				(node1+node2*4+node3) * (stepSize/6);
-
-			if (newDensity > desiredDensity) {
-				if (t+stepSize/2 <= maxDist) {
-					/* Record the last "good" scattering event */
-					success = true;
-					if (EXPECT_TAKEN(node2 != 0)) {
-						currentPoint = ray.o + halfIncrement;
-						currentSigmaT = node2;
-					} else if (node3 != 0 && t+stepSize <= maxDist) {
-						currentPoint = ray.o + fullIncrement;
-						currentSigmaT = node3;
-					} else if (node1 != 0) {
-						currentPoint = ray.o;
-						currentSigmaT = node1;
-					}
-				}
-
-				if (++numRefines > NUM_REFINES)
-					break;
-
-				stepSize *= .5f;
-				fullIncrement = halfIncrement;
-				halfIncrement *= .5f;
-				continue;
-			}
-
-			if (ray.o+fullIncrement == ray.o) /* Unable to make forward progress - stop */
-				break;
-
-			accumulatedDensity = newDensity;
-			node1 = node3;
-			t += stepSize;
-			ray.o += fullIncrement;
-		}
-
-		if (success) 
-			currentAlbedo = m_albedo->lookupSpectrum(currentPoint);
-
-		return success;
-	}
-
-	/// Evaluate the transmittance along a ray segment
-	Spectrum tau(const Ray &ray) const {
-		return Spectrum(std::exp(-integrateDensity(ray)));
-	}
-
-	bool sampleDistance(const Ray &r, MediumSamplingRecord &mRec,
-			Sampler *sampler) const {
-		Ray ray(r(r.mint), r.d, 0.0f);
-		Float distSurf       = r.maxt - r.mint,
-			  desiredDensity     = -std::log(1-sampler->next1D()),
-			  accumulatedDensity = 0.0f,
-			  currentSigmaT  = 0.0f;
-		Spectrum currentAlbedo(0.0f);
-		int iterations = 0;
-		bool inside = false, success = false;
-		Float t = 0;
-		Float sigmaTOrigin = sigmaT(ray.o, ray.d);
-
-		while (rayIntersect(ray, t, inside)) {
-			if (inside) {
-				Point currentPoint(0.0f);
-				success = invertDensityIntegral(ray, std::min(t, distSurf), 
-						accumulatedDensity, desiredDensity, currentSigmaT, currentAlbedo,
-						currentPoint);
-				if (success) {
-					/* A medium interaction occurred */
-					mRec.p = currentPoint;
-					mRec.t = (mRec.p-r.o).length();
-					success = true;
-					break;
-				}
-			}
-
-			distSurf -= t;
-
-			if (distSurf < 0) {
-				/* A surface interaction occurred */
-				break;
-			}
-			ray.o = ray(t);
-			ray.mint = Epsilon;
-
-			if (++iterations > 100) {
-				/// Just a precaution..
-				Log(EWarn, "sampleDistance(): round-off error issues?");
-				break;
-			}
-		}
-		Float expVal = std::max(Epsilon, std::exp(-accumulatedDensity));
-
-		if (success)
-			mRec.orientation = lookupOrientation(mRec.p);
-
-		mRec.pdfFailure = expVal;
-		mRec.pdfSuccess = expVal * std::max((Float) Epsilon, currentSigmaT);
-		mRec.pdfSuccessRev = expVal * std::max((Float) Epsilon, sigmaTOrigin);
-		mRec.transmittance = Spectrum(expVal);
-
-		if (!success)
-			return false;
-
-		mRec.sigmaS = currentAlbedo * currentSigmaT;
-		mRec.sigmaA = Spectrum(currentSigmaT) - mRec.sigmaS;
-		mRec.albedo = currentAlbedo.max();
-		mRec.medium = this;
-		return true;
-	}
-
-	void pdfDistance(const Ray &ray, Float t, MediumSamplingRecord &mRec) const {
-		Float expVal = std::exp(-integrateDensity(Ray(ray, 0, t)));
-
-		mRec.transmittance = Spectrum(expVal);
-		mRec.pdfFailure = expVal;
-		mRec.pdfSuccess = expVal * std::max((Float) Epsilon, 
-			sigmaT(ray(t), ray.d));
-		mRec.pdfSuccessRev = expVal * std::max((Float) Epsilon, 
-			sigmaT(ray(ray.mint), ray.d));
-	}
-
-	std::string toString() const {
-		std::ostringstream oss;
-		oss << "HeterogeneousFlakeMedium[" << endl
-			<< "  albedo=" << indent(m_albedo.toString()) << "," << endl
-			<< "  orientations=" << indent(m_orientations.toString()) << "," << endl
-			<< "  densities=" << indent(m_densities.toString()) << "," << endl
-			<< "  stepSize=" << m_stepSize << "," << endl
-			<< "  densityMultiplier=" << m_densityMultiplier << endl
-			<< "]";
-		return oss.str();
-	}
-	MTS_DECLARE_CLASS()
-private:
-	ref<VolumeDataSource> m_densities;
-	ref<VolumeDataSource> m_albedo;
-	ref<VolumeDataSource> m_orientations;
-	ref<ShapeKDTree> m_kdTree;
-	std::vector<Shape *> m_shapes;
-	Float m_stepSize;
-	/* Flake model-related */
-	SHVector m_D, m_sigmaT;
-	SHVector4D m_phaseExpansion;
-	int m_samplingRecursions;
-	EMode m_mode;
-	Float m_exponent, m_normalization;
-};
-
-Spectrum FlakePhaseFunction::f(const PhaseFunctionQueryRecord &pRec) const {
-	Frame frame(pRec.mRec.orientation);
-	Vector wi = frame.toLocal(pRec.wi);
-	Vector wo = frame.toLocal(pRec.wo);
-#if defined(FLAKE_EVAL_APPROXIMATION) 
-	/* Evaluate spherical harmonics representation - might be inaccurate */
-	SHVector temp(m_phaseExpansion->getBands());
-	m_phaseExpansion->lookup(wi, temp);
-	Float result = std::max((Float) 0, temp.eval(wo));
-	return Spectrum(result);
-#else
-	/* Evaluate the real phase function */
-	Float sigmaT = m_sigmaT.evalAzimuthallyInvariant(wi);
-	Vector H = wi + wo;
-	Float length = H.length();
-
-	if (length == 0)
-		return Spectrum(0.0f);
-	else
-		H /= length;
-
-	Float result = m_distr(H) / (2*sigmaT);
-	result *= m_normalization;
-	return Spectrum(result);
-#endif
-}
-
-Spectrum FlakePhaseFunction::sample(PhaseFunctionQueryRecord &pRec,
-		Float &_pdf, Sampler *sampler) const {
-#if defined(FLAKE_SAMPLE_UNIFORM)
-	/* Uniform sampling */
-	pRec.wo = squareToSphere(sampler->next2D());
-	_pdf = 1/(4*M_PI);
-	return f(pRec);
-#else
-	Point2 sample(sampler->next2D());
-	Frame frame(pRec.mRec.orientation);
-	Vector wi = frame.toLocal(pRec.wi);
-	/* Importance sampling using the interpolated phase function */
-	SHVector &temp = m_temp.get();
-	if (EXPECT_NOT_TAKEN(temp.getBands() == 0))
-		temp = SHVector(m_phaseExpansion->getBands());
-	m_phaseExpansion->lookup(wi, temp);
-	_pdf = m_shSampler->warp(temp, sample);
-
-	Vector wo = sphericalDirection(sample.x, sample.y);
-	pRec.wo = frame.toWorld(wo);
-	#if defined(FLAKE_EVAL_APPROXIMATION) 
-		return Spectrum(std::max((Float) 0, temp.eval(sample.x, sample.y)));
-	#else
-		return f(pRec);
-	#endif
-#endif
-}
-
-Spectrum FlakePhaseFunction::sample(PhaseFunctionQueryRecord &pRec,
-		Sampler *sampler) const {
-	Float pdf;
-	Spectrum result = FlakePhaseFunction::sample(pRec, pdf, sampler);
-	if (!result.isZero())
-		return result / pdf;
-	else
-		return Spectrum(0.0f);
-}
-
-Float FlakePhaseFunction::pdf(const PhaseFunctionQueryRecord &pRec) const {
-#if defined(FLAKE_SAMPLE_UNIFORM)
-	return 1/(4*M_PI);
-#else
-	Frame frame(pRec.mRec.orientation);
-	Vector wi = frame.toLocal(pRec.wi);
-	Vector wo = frame.toLocal(pRec.wo);
-	SHVector temp(m_phaseExpansion->getBands());
-	m_phaseExpansion->lookup(wi, temp);
-
-	return std::max((Float) 0, temp.eval(wo));
-#endif
-}
-
-MTS_IMPLEMENT_CLASS_S(HeterogeneousFlakeMedium, false, Medium)
-MTS_IMPLEMENT_CLASS_S(FlakePhaseFunction, false, PhaseFunction)
-MTS_EXPORT_PLUGIN(HeterogeneousFlakeMedium, "Heterogeneous micro-flake medium");
-MTS_NAMESPACE_END
-
diff --git a/src/medium/heterogeneous.cpp b/src/medium/heterogeneous.cpp
index db1e3e3b..e93a2b5b 100644
--- a/src/medium/heterogeneous.cpp
+++ b/src/medium/heterogeneous.cpp
@@ -356,7 +356,7 @@ public:
 
 					x -= fx/dfx;
 
-					if (x <= a || x >= b || dfx == 0) 
+					if (EXPECT_NOT_TAKEN(x <= a || x >= b || dfx == 0)) 
 						x = 0.5f * (b + a);
 
 					/* Integrated version of the above Lagrange polynomial */
diff --git a/src/phase/microflake.cpp b/src/phase/microflake.cpp
index 6c3aeecf..0b428b80 100644
--- a/src/phase/microflake.cpp
+++ b/src/phase/microflake.cpp
@@ -16,16 +16,24 @@
     along with this program. If not, see <http://www.gnu.org/licenses/>.
 */
 
+#include <mitsuba/core/chisquare.h>
 #include <mitsuba/render/phase.h>
 #include <mitsuba/render/sampler.h>
-#include "microflake.h"
+#include "microflake_fiber.h"
 
 MTS_NAMESPACE_BEGIN
 
 class MicroflakePhaseFunction : public PhaseFunction {
 public:
-	MicroflakePhaseFunction(const Properties &props) 
-		: PhaseFunction(props) {
+	MicroflakePhaseFunction(const Properties &props) : PhaseFunction(props) {
+		m_fiberDistr = GaussianFiberDistribution(props.getFloat("stddev"));
+		ChiSquareTest test(7);
+		test.fill(
+			boost::bind(&GaussianFiberDistribution::sample, m_fiberDistr, _1),
+			boost::bind(&GaussianFiberDistribution::pdf, m_fiberDistr, _1)
+		);
+		test.dumpTables("test.m");
+		test.runTest(1);
 	}
 
 	MicroflakePhaseFunction(Stream *stream, InstanceManager *manager) 
@@ -60,10 +68,16 @@ public:
 	}
 
 	std::string toString() const {
-		return "MicroflakePhaseFunction[]";
+		std::ostringstream oss;
+		oss << "MicroflakePhaseFunction[" << endl
+			<< "   fiberDistr = " << indent(m_fiberDistr.toString()) << endl
+			<< "]";
+		return oss.str();
 	}
 
 	MTS_DECLARE_CLASS()
+private:
+	GaussianFiberDistribution m_fiberDistr;
 };
 
 
diff --git a/src/phase/microflake.h b/src/phase/microflake_fiber.h
similarity index 88%
rename from src/phase/microflake.h
rename to src/phase/microflake_fiber.h
index 05d8a033..dfd3f8ff 100644
--- a/src/phase/microflake.h
+++ b/src/phase/microflake_fiber.h
@@ -16,20 +16,30 @@
     along with this program. If not, see <http://www.gnu.org/licenses/>.
 */
 
-#if !defined(__MICROFLAKE_SIGMA_T_H)
-#define __MICROFLAKE_SIGMA_T_H
+#if !defined(__MICROFLAKE_FIBER_DIST_H)
+#define __MICROFLAKE_FIBER_DIST_H
 
 #include <mitsuba/core/statistics.h>
 #include <mitsuba/core/brent.h>
-#include <boost/math/special_functions/erf.hpp>
 #include <boost/bind.hpp>
+#include <boost/math/special_functions/erf.hpp>
 
 MTS_NAMESPACE_BEGIN
 
-#define SIGMA_T_RANGE_MIN 4e-08
-#define SIGMA_T_RANGE_MAX 4
-#define SIGMA_T_ELEMENTS  100
-#define SIGMA_T_COEFFS    10
+FINLINE Float mts_erf(Float arg) {
+	#if defined(__GNUC__) && defined(SINGLE_PRECISION)
+		return erff(arg);
+	#elif defined(__GNUC__) && defined(DOUBLE_PRECISION)
+		return erf(arg);
+	#else
+		return boost::math::erf<Float>(arg);
+	#endif
+}
+
+#define FIBERDIST_STDDEV_MIN        4e-08
+#define FIBERDIST_STDDEV_MAX        4
+#define FIBERDIST_SIGMA_T_ELEMENTS  100
+#define FIBERDIST_SIGMA_T_COEFFS    10
 
 /**
  * Each row of this table table stores an expansion of \sigma_t in terms of
@@ -37,11 +47,11 @@ MTS_NAMESPACE_BEGIN
  * fiber axis and the incident direction. The standard deviation of the
  * underlying distribution increases from row to row using a nonlinear 
  * mapping to capture the `interesting' parts more efficiently (see
- * \ref sigmaT_fiberGaussian())
+ * \ref sigmaT_fiberSigmaT())
  *
  * The implementation here has an average absolute error of 3.17534e-05.
  */
-const double fiberGaussianCoeffs[SIGMA_T_ELEMENTS][SIGMA_T_COEFFS] = {
+const double fiberSigmaTCoeffs[FIBERDIST_SIGMA_T_ELEMENTS][FIBERDIST_SIGMA_T_COEFFS] = {
 	{ 2.5037347588631811e-08,   6.3661860705430207e-01,   1.7827351092236654e-05,  -1.3015084365974872e-04,   5.2848718115683369e-04,  -1.2810825925271274e-03,   1.8970200528656278e-03,  -1.6818874867112754e-03,   8.1991128831759852e-04,  -1.6898472595983094e-04 },
 	{ 4.0059757909494120e-07,   6.3660112034736627e-01,   2.8523763553955916e-04,  -2.0824136415171779e-03,   8.4557955473769653e-03,  -2.0497323170502568e-02,   3.0352323203715059e-02,  -2.6910202150020268e-02,   1.3118581833396092e-02,  -2.7037558688078889e-03 },
 	{ 2.0280252492606374e-06,   6.3652534713591225e-01,   1.4440155345667449e-03,  -1.0542219099960448e-02,   4.2807465147063795e-02,  -1.0376769900017280e-01,   1.5365863698661997e-01,  -1.3623289910560743e-01,   6.6412820861387445e-02,  -1.3687764171763206e-02 },
@@ -144,27 +154,27 @@ const double fiberGaussianCoeffs[SIGMA_T_ELEMENTS][SIGMA_T_COEFFS] = {
 	{ 4.9740122861160618e-01,   4.9788891152147130e-06,   3.7427454517739989e-03,   1.2662117029549336e-03,  -3.7139444921194809e-03,   2.9874687143092160e-03,   6.6568975580594270e-03,  -1.6932579956119298e-02,   1.3976501563774946e-02,  -4.0922961870819563e-03 }
 };
 
-double sigmaT_fiberGaussian(double stddev, double sinTheta) {
+double sigmaT_fiberDist(double stddev, double sinTheta) {
 	/* Undo the nonlinear mapping */
-	const double factor = SIGMA_T_ELEMENTS / std::pow(SIGMA_T_RANGE_MAX, 0.25);
+	const double factor = FIBERDIST_SIGMA_T_ELEMENTS / std::pow(FIBERDIST_STDDEV_MAX, 0.25);
 	double pos = std::pow(stddev, 0.25) * factor - 1;
 
 	/* Clamp to the supported parameter range */
-	pos = std::min(std::max(0.0, pos), (double) (SIGMA_T_ELEMENTS-1));
+	pos = std::min(std::max(0.0, pos), (double) (FIBERDIST_SIGMA_T_ELEMENTS-1));
 
 	/* Linearly interpolate coefficients using \c alpha */
 	int idx0 = (int) std::floor(pos), idx1 = (int) std::ceil(pos);
 	double alpha = pos - idx0, base = 1, result = 0;
 
 	SAssert(
-		idx0 >= 0 && idx0 < SIGMA_T_ELEMENTS &&
-		idx1 >= 0 && idx1 < SIGMA_T_ELEMENTS
+		idx0 >= 0 && idx0 < FIBERDIST_SIGMA_T_ELEMENTS &&
+		idx1 >= 0 && idx1 < FIBERDIST_SIGMA_T_ELEMENTS
 	);
 
-	/* Compute the expansion */
-	for (int i=0; i<SIGMA_T_COEFFS; ++i) {
-		result += base * ((1-alpha) * fiberGaussianCoeffs[idx0][i]
-				+ alpha * fiberGaussianCoeffs[idx1][i]);
+	/* Evaluate the expansion */
+	for (int i=0; i<FIBERDIST_SIGMA_T_COEFFS; ++i) {
+		result += base * ((1-alpha) * fiberSigmaTCoeffs[idx0][i]
+				+ alpha * fiberSigmaTCoeffs[idx1][i]);
 		base *= sinTheta;
 	}
 
@@ -176,56 +186,106 @@ static StatsCounter brentSolves("Micro-flake model",
 static StatsCounter avgBrentFunEvals("Micro-flake model",
 		"Average Brent solver function evaluations", EAverage);
 
-struct GaussianFiberDistribution {
+class GaussianFiberDistribution {
 public:
+	inline GaussianFiberDistribution() {}
+
 	inline GaussianFiberDistribution(Float stddev) : m_stddev(stddev) {
-		m_normalization = 1/(std::sqrt(2*M_PI) * m_stddev * 
-				boost::math::erf<Float>(1/(SQRT_TWO * m_stddev)));
-		m_c1 = 1.0f/boost::math::erf<Float>(1/(SQRT_TWO * m_stddev));
+		m_normalization = 1/(std::pow(2*M_PI, 3.0f/2.0f) * m_stddev * 
+				mts_erf(1/(SQRT_TWO * m_stddev)));
+		m_c1 = 1.0f/mts_erf(1/(SQRT_TWO * m_stddev));
+
+		if (stddev < FIBERDIST_STDDEV_MIN || stddev > FIBERDIST_STDDEV_MAX)
+			SLog(EError, "Standard deviation parameter is out of range (must "
+				"be in [%f, %f])!", FIBERDIST_STDDEV_MIN, FIBERDIST_STDDEV_MAX);
+
+		/* Determine expansion coefficients of sigma_t for a fixed stddev */
+		Float pos = std::pow(stddev / FIBERDIST_STDDEV_MAX, 0.25) * 
+			FIBERDIST_SIGMA_T_ELEMENTS - 1;
+
+		pos = std::min(std::max((Float) 0, pos), (Float) (FIBERDIST_SIGMA_T_ELEMENTS-1));
+
+		int idx0 = (int) std::floor(pos), idx1 = (int) std::ceil(pos);
+		Float alpha = pos - idx0;
+
+		for (int i=0; i<FIBERDIST_SIGMA_T_COEFFS; ++i)
+			m_coeffs[i] = (Float) (((1-alpha) * fiberSigmaTCoeffs[idx0][i]
+				+ alpha * fiberSigmaTCoeffs[idx1][i]));
 	}
 
 	/// Evaluate \sigma_t as a function of \cos\theta
-	Float sigmaT(Float cosTheta) const {
-		return sigmaT_fiberGaussian(m_stddev, 
-			std::sqrt(std::max((Float) 0, 1-cosTheta*cosTheta)));
+	inline Float sigmaT(Float cosTheta) const {
+		Float sinTheta = std::sqrt(std::max(
+				(Float) 0, 1-cosTheta*cosTheta)),
+			  base = 1.0f, result = 0.0f;
+
+		/* Evaluate the expansion */
+		for (int i=0; i<FIBERDIST_SIGMA_T_COEFFS; ++i) {
+			result += base * m_coeffs[i];
+			base *= sinTheta;
+		}
+
+		return result;
 	}
 
-	/// Evaluate the longitudinal density as a function of \cos\theta
-	Float pdfCosTheta(Float cosTheta) const {
-		return std::exp(-cosTheta*cosTheta/(2*m_stddev*m_stddev)) * m_normalization;
+	/// Evaluate the density as a function of \cos\theta
+	inline Float pdfCosTheta(Float cosTheta) const {
+		return std::exp(-cosTheta*cosTheta 
+			/ (2*m_stddev*m_stddev)) * m_normalization;
 	}
 
-	/// Evaluate the longitudinal CDF as a function of \cos\theta
-	inline Float cdf(Float theta) const {
-		return 0.5f * (1.0f - 
-			boost::math::erf<Float>(std::cos(theta)/(SQRT_TWO * m_stddev))*m_c1);
+	/// Evaluate the density as a function of direction
+	inline Float pdf(const Vector &d) const {
+		return std::exp(-d.z*d.z
+			/ (2*m_stddev*m_stddev)) * m_normalization;
 	}
 
 	/**
-	 * Apply the inversion method to sample \theta given 
-	 * a uniformly distributed \xi on [0, 1]
+	 * \brief Apply the inversion method to sample \cos\theta given 
+	 * a uniformly distributed r.v. \xi on [0, 1]
 	 */
-	Float sample(Float xi) const {
+	Vector sample(const Point2 &sample) const {
 		BrentSolver brentSolver(100, 1e-6f);
 		BrentSolver::Result result = brentSolver.solve(
-			boost::bind(&GaussianFiberDistribution::cdfFunctor, this, xi, _1), 0, M_PI);
+			boost::bind(&GaussianFiberDistribution::cdfFunctor, 
+				this, sample.x, _1), -1, 1);
 		SAssert(result.success);
 		avgBrentFunEvals.incrementBase();
 		++brentSolves;
-		return result.x;
-	}
-private:
-	Float cdfFunctor(Float theta, Float xi) const {
-		++avgBrentFunEvals;
-		return cdf(theta)-xi;
+
+		Float cosTheta = result.x, 
+			  sinTheta = std::sqrt(std::max((Float) 0, 1-cosTheta*cosTheta)),
+			  phi = 2 * M_PI * sample.y,
+			  sinPhi = std::sin(phi), cosPhi = std::cos(phi);
+
+		return Vector(sinTheta * cosPhi, sinTheta * sinPhi, cosTheta);
 	}
 
+	std::string toString() const {
+		std::ostringstream oss;
+		oss << "GaussianFiberDistribution[stddev="
+			<< m_stddev << "]";
+		return oss.str();
+	}
+protected:
+	/// Evaluate the longitudinal CDF as a function of \cos\theta
+	inline Float cdf(Float cosTheta) const {
+		return 0.5f * (1.0f - 
+			mts_erf(cosTheta / (SQRT_TWO * m_stddev)) * m_c1);
+	}
+
+	Float cdfFunctor(Float xi, Float cosTheta) const {
+		++avgBrentFunEvals;
+		return cdf(cosTheta)-xi;
+	}
+
+protected:
 	Float m_stddev;
 	Float m_normalization;
 	Float m_c1;
+	Float m_coeffs[FIBERDIST_SIGMA_T_COEFFS];
 };
 
-
 MTS_NAMESPACE_END
 
-#endif /* __MICROFLAKE_SIGMA_T_H */
+#endif /* __MICROFLAKE_FIBER_DIST_H */