sky luminaire by Tom Kazimiers, composite light sources, rayleigh scattering in media

data/tests/test_bsdf.xml

@@ -2,6 +2,19 @@
 	 to be tested for consistency. This is done
 	 using the testcase 'test_chisquare' -->
 <scene version="0.3.0">
+	<bsdf type="roughdiffuse">
+		<float name="alpha" value="1"/>
+	</bsdf>
+
+	<!-- Test the rough dielectric model with the anisotropic
+		 Ashikhmin-Shirley microfacet distribution -->
+	<bsdf type="roughdielectric">
+		<string name="distribution" value="ggx"/>
+		<float name="alpha" value="3"/>
+		<float name="intIOR" value="1.5"/>
+		<float name="extIOR" value="1.0"/>
+	</bsdf>
+
 	<!-- Test the coating model with a transmissive
 		 + reflective material -->
 	<bsdf type="coating">

data/tests/test_phase.xml

@@ -4,7 +4,10 @@
 <scene version="0.3.0">
 	<!-- Test the isotropic phase function -->
 	<phase type="isotropic"/>
-	
+
+	<!-- Test the Rayleigh phase function -->
+	<phase type="rayleigh"/>
+
 	<!-- Test the Henyey-Greenstein phase function
 		 configured as highly forward-scattering -->
 	<phase type="hg">

doc/acknowledgements.tex

@@ -13,6 +13,15 @@ and measurements of atomic scattering factors made by the Center For
 X-Ray Optics (CXRO) at Berkeley and the Lawrence Livermore National 
 Laboratory (LLNL).
 
+The following people have contributed code or bugfixes:
+\begin{itemize}
+\item Milo\^{s} Ha\^{s}an
+\item Tom Kazimiers
+\item Marios Papas
+\item Edgar Vel\'{a}zquez-Armend\'{a}riz
+\item Jirka Vorba
+\end{itemize}
+
 Mitsuba makes heavy use of the following amazing libraries and tools: 
 \begin{itemize}
 \item Qt 4 by Nokia

doc/gendoc.py

@@ -59,6 +59,8 @@ f = open('plugins_generated.tex', 'w')
 f.write('\section{Plugin reference}\n')
 f.write('\input{section_bsdf}\n')
 os.path.walk('../src/bsdfs', traverse, f)
+f.write('\input{section_phase}\n')
+os.path.walk('../src/phase', traverse, f)
 f.write('\input{section_subsurface}\n')
 os.path.walk('../src/subsurface', traverse, f)
 f.write('\input{section_media}\n')

doc/main.bib

@@ -154,3 +154,44 @@
 	address = {New York, NY, USA},
 } 
 
+@article{Jakob2010Radiative,
+	title={A radiative transfer framework for rendering materials with anisotropic structure},
+	author={Jakob, W. and Arbree, A. and Moon, J.T. and Bala, K. and Marschner, S.},
+	journal={ACM Transactions on Graphics (TOG), Proceedings of SIGGRAPH 2010},
+	volume={29},
+	number={4},
+	pages={53},
+	year={2010},
+	publisher={ACM}
+}
+
+@article{Zhao2011Building,
+	title={{Building Volumetric Appearance Models of Fabric using Micro CT Imaging}},
+	author={Zhao, Shuang and Jakob, Wenzel and Marschner, Steve and Bala, Kavita},
+	journal={ACM Transactions on Graphics (TOG), Proceedings of SIGGRAPH 2011},
+	volume={30},
+	number={4},
+	pages={53},
+	year={2011},
+	publisher={ACM}
+}
+
+@article{Henyey1941Diffuse,
+	title={Diffuse radiation in the galaxy},
+	author={Henyey, L.G. and Greenstein, J.L.},
+	journal={The Astrophysical Journal},
+	volume={93},
+	pages={70--83},
+	year={1941}
+}
+
+@inproceedings{Kajiya1989Rendering,
+	title={Rendering fur with three dimensional textures},
+	author={Kajiya, J.T. and Kay, T.L.},
+	booktitle={ACM Siggraph Computer Graphics},
+	volume={23},
+	number={3},
+	pages={271--280},
+	year={1989},
+	organization={ACM}
+}

doc/section_phase.tex

@@ -0,0 +1,18 @@
+\newpage
+\subsection{Phase functions}
+\label{sec:phase}
+This section contains a description of all implemented medium scattering models, which 
+are also known as \emph{phase functions}. These are very similar in principle to surface 
+scattering models (or \emph{BSDF}s), and essentially describe where light travels after 
+hitting a particle within the medium.
+
+The most commonly used models for smoke, fog, and other homogeneous media
+are isotropic scattering (\pluginref{isotropic}) and the Henyey-Greenstein 
+phase function (\pluginref{hg}). Mitsuba also supports \emph{anisotropic}
+media, where the behavior of the medium changes depending on the direction 
+of light propagation (e.g. in volumetric representations of fabric). These 
+are the Kajiya-Kay (\pluginref{kkay}) and Micro-flake (\pluginref{microflake}) 
+models.
+
+Finally, there is also a phase function for simulating scattering in
+planetary atmospheres (\pluginref{rayleigh}).

include/mitsuba/core/properties.h

@@ -30,14 +30,30 @@ MTS_NAMESPACE_BEGIN
  */
 class MTS_EXPORT_CORE Properties {
 public:
+	/// Supported types of properties
 	enum EPropertyType {
+		/// Boolean value (true/false)
 		EBoolean = 0,
+		/// 64-bit signed integer
 		EInteger,
+		/// Floating point value
 		EFloat,
+		/// 3D point
 		EPoint,
+		/// 4x4 transform for homogeneous coordinates
 		ETransform,
+		/// Discretized color spectrum
 		ESpectrum,
-		EString
+		/// Arbitrary-length string
+		EString,
+		/// Arbitrary data (pointer+size)
+		EData
+	};
+
+	/// Simple pointer-size pair for passing arbitrary data (e.g. between plugins)
+	struct Data {
+		uint8_t *ptr;
+		size_t size;
 	};
 
 	/// Construct an empty property container
@@ -88,6 +104,13 @@ public:
 	/// Get a single precision floating point value (with default)
 	Float getFloat(const std::string &name, Float defVal) const;
 
+	/// Set an arbitrary data value
+	void setData(const std::string &name, Data value, bool warnDuplicates = true);
+	/// Get an arbitrary data value
+	Data getData(const std::string &name) const;
+	/// Get an arbitrary data value (with default)
+	Data getData(const std::string &name, Data defVal) const;
+
 	/// Set a linear transformation
 	void setTransform(const std::string &name, const Transform &value, bool warnDuplicates = true);
 	/// Get a linear transformation
@@ -146,6 +169,7 @@ private:
 			bool v_boolean;
 			int64_t v_long;
 			Float v_float;
+			Data v_data;
 		};
 		// not allowed in union (constructor)
 		Point v_point; 

include/mitsuba/core/spectrum.h

@@ -110,6 +110,31 @@ private:
 	Float m_temperature;
 };
 
+/**
+ * \brief Spectral power distribution used to model Rayleigh scattering
+ *
+ * This distribution captures the wavelength dependence of Rayleigh
+ * scattering.
+ *
+ * \ingroup libcore
+ */
+class MTS_EXPORT_CORE RayleighSpectrum : public ContinuousSpectrum {
+public:
+	/// Compute a new rayleigh spectrum
+	inline RayleighSpectrum() { }
+
+	virtual ~RayleighSpectrum() { }
+
+	/** \brief Return the value of the spectral power distribution
+	 * at the given wavelength.
+	 *
+	 * The units are Watts per unit surface area (m^-2) 
+	 * per unit wavelength (nm^-1) per steradian (sr^-1)
+	 */
+	virtual Float eval(Float lambda) const;
+};
+
+
 /**
  * \brief This spectral power distribution is defined as the
  * product of two other continuous spectra.

include/mitsuba/render/bsdf.h

@@ -336,7 +336,7 @@ public:
 	virtual Spectrum sample(BSDFQueryRecord &bRec, const Point2 &sample) const = 0;
 
 	/**
-	 * \c Sample the BSDF and explicitly provide the probability density
+	 * \brief Sample the BSDF and explicitly provide the probability density
 	 * of the sampled direction. 
 	 *
 	 * If a component mask or a specific component index is given, the 

include/mitsuba/render/luminaire.h

@@ -342,6 +342,18 @@ public:
 	//! @{ \name Miscellaneous
 	// =============================================================
 
+	/// Is this a compound luminaire consisting of several sub-objects?
+	virtual bool isCompound() const;
+
+	/**
+	 * \brief Return a sub-element of a compound luminaire. 
+	 *
+	 * When expanding luminaires, the scene will repeatedly call this
+	 * function with increasing indices. Returning \a NULL indicates
+	 * that no more are available.
+	 */
+	virtual Luminaire *getElement(int i);
+
 	/// Serialize this luminaire to a binary data stream
 	virtual void serialize(Stream *stream, InstanceManager *manager) const;
 

include/mitsuba/render/shape.h

@@ -195,7 +195,6 @@ public:
 	 */
 	virtual Shape *getElement(int i);
 
-
 	/**
 	 * \brief Return the shape's surface area
 	 *

src/bsdfs/bump.cpp

@@ -59,11 +59,11 @@ MTS_NAMESPACE_BEGIN
  *     <!-- The bump map is applied to a rough metal BRDF -->
  *     <bsdf type="roughconductor"/>
  *
- *     <texture name="scale">
+ *     <texture type="scale">
  *         <!-- The scale of the displacement gets multiplied by 10x -->
  *         <float name="scale" value="10"/>
  *
- *         <texture name="bitmap">
+ *         <texture type="bitmap">
  *             <string name="filename" value="bumpmap.png"/>
  *         </texture>
  *     </texture>

src/bsdfs/mask.cpp

@@ -45,6 +45,7 @@ MTS_NAMESPACE_BEGIN
  * \begin{xml}[caption=Material configuration for a transparent leaf,
  *    label=lst:mask-leaf]
  * <bsdf type="mask">
+ *     <!-- Base material: a two-sided textured diffuse BSDF -->
  *     <bsdf type="twosided">
  *         <bsdf type="diffuse">
  *             <texture name="reflectance" type="bitmap">
@@ -52,6 +53,8 @@ MTS_NAMESPACE_BEGIN
  *             </texture>
  *         </bsdf>
  *     </bsdf>
+ *
+ *     <!-- Fetch the opacity mask from a bitmap -->
  *     <texture name="opacity" type="bitmap">
  *         <string name="filename" value="leaf_opacity.jpg"/>
  *         <float name="gamma" value="1"/>

src/libcore/properties.cpp

@@ -160,6 +160,35 @@ size_t Properties::getSize(const std::string &name, size_t defVal) const {
 	return (size_t) (*it).second.v_long;
 }
 
+void Properties::setData(const std::string &name, Data value, bool warnDuplicates) {
+	if (hasProperty(name) && warnDuplicates)
+		SLog(EWarn, "Property \"%s\" has already been specified!", name.c_str());
+	m_elements[name].type = EData;
+	m_elements[name].v_data = value;
+	m_elements[name].queried = false;
+}
+
+Properties::Data Properties::getData(const std::string &name) const {
+	if (!hasProperty(name))
+		SLog(EError, "Property \"%s\" missing", name.c_str());
+	std::map<std::string, Element>::const_iterator it = m_elements.find(name);
+	if ((*it).second.type != EData)
+		SLog(EError, "The property \"%s\" has the wrong type (expected <data>). The detailed "
+				"listing was:\n%s", name.c_str(), toString().c_str());
+	(*it).second.queried = true;
+	return (*it).second.v_data;
+}
+
+Properties::Data Properties::getData(const std::string &name, Data defVal) const {
+	if (!hasProperty(name))
+		return defVal;
+	std::map<std::string, Element>::const_iterator it = m_elements.find(name);
+	if ((*it).second.type != EData)
+		SLog(EError, "The property \"%s\" has the wrong type (expected <data>). The detailed "
+				"listing was:\n%s", name.c_str(), toString().c_str());
+	(*it).second.queried = true;
+	return (*it).second.v_data;
+}
 
 void Properties::setFloat(const std::string &name, Float value, bool warnDuplicates) {
 	if (hasProperty(name) && warnDuplicates)
@@ -351,6 +380,9 @@ std::string Properties::toString() const {
 			case EString:
 				oss << "\"" << el.v_string << "\"";
 				break;
+			case EData:
+				oss << el.v_data.ptr << " (size=" << el.v_data.size << ")";
+				break;
 			default:
 				SLog(EError, "Encountered an unknown property type!");
 		}

src/libcore/spectrum.cpp

@@ -444,6 +444,12 @@ Float ProductSpectrum::eval(Float lambda) const {
 	return m_spec1.eval(lambda) * m_spec2.eval(lambda);
 }
 
+Float RayleighSpectrum::eval(Float lambda) {
+	Float lambdaSqr = lambda*lambda;
+
+	return 1.0f / (lambdaSqr*lambdaSqr);
+}
+
 Float ContinuousSpectrum::average(Float lambdaMin, Float lambdaMax) const {
 	GaussLobattoIntegrator integrator(10000, Epsilon, Epsilon, false, false);
 

src/librender/luminaire.cpp

@@ -95,6 +95,14 @@ Spectrum Luminaire::Le(const ShapeSamplingRecord &sRec, const Vector &d) const {
 	return Spectrum(0.0f);
 }
 
+bool Luminaire::isCompound() const {
+	return false;
+}
+
+Luminaire *Luminaire::getElement(int i) {
+	return NULL;
+}
+
 std::string EmissionRecord::toString() const {
 	std::ostringstream oss;
 	oss << "EmissionRecord[" << std::endl

src/librender/scene.cpp

@@ -147,44 +147,44 @@ Scene::Scene(Stream *stream, InstanceManager *manager)
 	m_aabb = AABB(stream);
 	m_bsphere = BSphere(stream);
 	m_backgroundLuminaire = static_cast<Luminaire *>(manager->getInstance(stream));
-	int count = stream->readInt();
+	size_t count = stream->readSize();
-	for (int i=0; i<count; ++i) {
+	for (size_t i=0; i<count; ++i) {
 		Shape *shape = static_cast<Shape *>(manager->getInstance(stream));
 		shape->incRef();
 		m_shapes.push_back(shape);
 	}
-	count = stream->readInt();
+	count = stream->readSize();
-	for (int i=0; i<count; ++i) {
+	for (size_t i=0; i<count; ++i) {
 		TriMesh *trimesh = static_cast<TriMesh *>(manager->getInstance(stream));
 		trimesh->incRef();
 		m_meshes.push_back(trimesh);
 	}
-	count = stream->readInt();
+	count = stream->readSize();
-	for (int i=0; i<count; ++i) {
+	for (size_t i=0; i<count; ++i) {
 		Luminaire *luminaire = static_cast<Luminaire *>(manager->getInstance(stream));
 		luminaire->incRef();
 		m_luminaires.push_back(luminaire);
 	}
-	count = stream->readInt();
+	count = stream->readSize();
-	for (int i=0; i<count; ++i) {
+	for (size_t i=0; i<count; ++i) {
 		Medium *medium = static_cast<Medium *>(manager->getInstance(stream));
 		medium->incRef();
 		m_media.insert(medium);
 	}
-	count = stream->readInt();
+	count = stream->readSize();
-	for (int i=0; i<count; ++i) {
+	for (size_t i=0; i<count; ++i) {
 		Subsurface *ssIntegrator = static_cast<Subsurface *>(manager->getInstance(stream));
 		ssIntegrator->incRef();
 		m_ssIntegrators.push_back(ssIntegrator);
 	}
-	count = stream->readInt();
+	count = stream->readSize();
-	for (int i=0; i<count; ++i) {
+	for (size_t i=0; i<count; ++i) {
 		ConfigurableObject *obj = static_cast<ConfigurableObject *>(manager->getInstance(stream));
 		obj->incRef();
 		m_objects.push_back(obj);
 	}
-	count = stream->readInt();
+	count = stream->readSize();
-	for (int i=0; i<count; ++i) {
+	for (size_t i=0; i<count; ++i) {
 		NetworkedObject *obj = static_cast<NetworkedObject *>(manager->getInstance(stream));
 		m_netObjects.push_back(obj); // Do not increase the ref. count
 	}
@@ -577,6 +577,18 @@ void Scene::addChild(const std::string &name, ConfigurableObject *child) {
 		}
 	} else if (cClass->derivesFrom(MTS_CLASS(Luminaire))) {
 		Luminaire *luminaire = static_cast<Luminaire *>(child);
+
+		if (luminaire->isCompound()) {
+			int index = 0;
+			do {
+				ref<Luminaire> element = luminaire->getElement(index++);
+				if (element == NULL)
+					break;
+				addChild(name, element);
+			} while (true);
+			return;
+		}
+
 		luminaire->incRef();
 		m_luminaires.push_back(luminaire);
 		if (luminaire->isBackgroundLuminaire()) {
@@ -700,26 +712,26 @@ void Scene::serialize(Stream *stream, InstanceManager *manager) const {
 	m_aabb.serialize(stream);
 	m_bsphere.serialize(stream);
 	manager->serialize(stream, m_backgroundLuminaire.get());
-	stream->writeUInt((uint32_t) m_shapes.size());
+	stream->writeSize(m_shapes.size());
 	for (size_t i=0; i<m_shapes.size(); ++i) 
 		manager->serialize(stream, m_shapes[i]);
-	stream->writeUInt((uint32_t) m_meshes.size());
+	stream->writeSize(m_meshes.size());
 	for (size_t i=0; i<m_meshes.size(); ++i) 
 		manager->serialize(stream, m_meshes[i]);
-	stream->writeUInt((uint32_t) m_luminaires.size());
+	stream->writeSize(m_luminaires.size());
 	for (size_t i=0; i<m_luminaires.size(); ++i) 
 		manager->serialize(stream, m_luminaires[i]);
-	stream->writeUInt((uint32_t) m_media.size());
+	stream->writeSize(m_media.size());
 	for (std::set<Medium *>::const_iterator it = m_media.begin();
 			it != m_media.end(); ++it)
 		manager->serialize(stream, *it);
-	stream->writeUInt((uint32_t) m_ssIntegrators.size());
+	stream->writeSize(m_ssIntegrators.size());
 	for (size_t i=0; i<m_ssIntegrators.size(); ++i) 
 		manager->serialize(stream, m_ssIntegrators[i]);
-	stream->writeUInt((uint32_t) m_objects.size());
+	stream->writeSize(m_objects.size());
 	for (size_t i=0; i<m_objects.size(); ++i) 
 		manager->serialize(stream, m_objects[i]);
-	stream->writeUInt((uint32_t) m_netObjects.size());
+	stream->writeSize(m_netObjects.size());
 	for (size_t i=0; i<m_netObjects.size(); ++i) 
 		manager->serialize(stream, m_netObjects[i]);
 }

src/luminaires/SConscript

@@ -7,5 +7,6 @@ plugins += env.SharedLibrary('spot', ['spot.cpp'])
 plugins += env.SharedLibrary('point', ['point.cpp'])
 plugins += env.SharedLibrary('collimated', ['collimated.cpp'])
 plugins += env.SharedLibrary('directional', ['directional.cpp'])
+plugins += env.SharedLibrary('sky', ['sky.cpp'])
 
 Export('plugins')

src/luminaires/sky.cpp

@@ -0,0 +1,560 @@
+/*
+	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/core/util.h>
+#include <mitsuba/core/bitmap.h>
+#include <mitsuba/core/fstream.h>
+
+#define SAMPLE_UNIFORMLY 1
+
+/* Define this to automatically set the viewing angles phi to zero if
+ * its theta is zero. This means that there is only one possible
+ * viewing zenith angle. Some implementitions do this, but it is
+ * unclear why one should use it. */
+// #define ENFORCE_SINGLE_ZENITH_ANGLE
+
+MTS_NAMESPACE_BEGIN
+
+/*
+ * A sun and skylight luminaire. In its local coordinate system, the sun
+ * is "above" on positive Y direction. So when configuring, keep in mind
+ * that south = x, east = y and up = z. All times in decimal form (6.25
+ * = 6:15 AM) and all angles in radians.
+ *
+ * The model behind it is described by Preetham et al. (2002).
+ */
+class SkyLuminaire : public Luminaire {
+public:
+	/**
+	 * Creates a new Sky luminaire. The defaults values origate
+	 * from the sample code of Preetham et al. and create an over
+	 * head sun on a clear day.
+	 */
+	SkyLuminaire(const Properties &props)
+			: Luminaire(props) {
+		/* Transformation from the luminaire's local coordinates to
+		 *  world coordiantes */
+		m_luminaireToWorld = Transform::rotate(Vector(1, 0, 0),-90)
+			* props.getTransform("toWorld", Transform());
+		m_worldToLuminaire = m_luminaireToWorld.inverse();
+
+		m_skyScale = props.getFloat("skyScale", Float(1.0));
+		m_turbidity = props.getFloat("turbidity", Float(2.0));
+		m_aConst = props.getFloat("aConst", 1.0);
+		m_bConst = props.getFloat("bConst", 1.0);
+		m_cConst = props.getFloat("cConst", 1.0);
+		m_dConst = props.getFloat("dConst", 1.0);
+		m_eConst = props.getFloat("eConst", 1.0);
+		m_clipBelowHorizon = props.getBoolean("clipBelowHorizon", true);
+		m_exposure = props.getFloat("exposure", 1.0/15.0);
+
+		/* Do some input checks for sun position information */
+		bool hasSunDir = props.hasProperty("sunDirection");
+		bool hasLatitude = props.hasProperty("latitude");
+		bool hasLongitude = props.hasProperty("longitude");
+		bool hasStdMrd = props.hasProperty("standardMeridian");
+		bool hasJulDay = props.hasProperty("julianDay");
+		bool hasTimeOfDay = props.hasProperty("timeOfDay");
+
+		bool hasSomeLocInfo = hasLatitude || hasLongitude || hasStdMrd || hasJulDay || hasTimeOfDay;
+		bool hasAllLocInfo = hasLatitude && hasLongitude && hasStdMrd && hasJulDay && hasTimeOfDay;
+
+		if (!(hasSunDir || hasSomeLocInfo)) {
+			/* If no sun position has been specified, use a default one. */
+			hasSomeLocInfo = hasAllLocInfo = true;
+		} else if (hasSunDir && hasSomeLocInfo) {
+			/* We don't allow the use of both position formats, raise error. */
+			throw std::runtime_error("Please decide for either positioning the sun by a direction vector or by some location and time information.");
+		} else if (hasSomeLocInfo && !hasAllLocInfo) {
+			/* The sun positioning by location and time information is missing
+			* some parameters in the input, raise error. */
+			throw std::runtime_error("Please give all required parameters for specifing the sun's position by location and time information. At least one is missing.");
+		} /* else we have complete positioning information */
+
+		/* The direction should always relate to "up" beeing in positive Z */
+		Vector sunDir = props.getVector("sunDirection", Vector(0.0, 0.0, 1.0));
+		Float lat = props.getFloat("latitude", 51.050891);
+		Float lon  = props.getFloat("longitude", 13.694458);
+		Float stdMrd = props.getFloat("standardMeridian", 0);
+		Float julianDay = props.getFloat("julianDay", 200);
+		Float timeOfDay = props.getFloat("timeOfDay", 15.00);
+		
+		/* configure position of sun */
+		if (hasSunDir)
+			configureSunPosition(sunDir);
+		else
+			configureSunPosition(lat, lon, stdMrd, julianDay, timeOfDay);
+
+		configure();
+	}
+
+	SkyLuminaire(Stream *stream, InstanceManager *manager) 
+		    : Luminaire(stream, manager) {
+		m_skyScale = stream->readFloat();
+		m_turbidity = stream->readFloat();
+		m_thetaS = stream->readFloat();
+		m_phiS = stream->readFloat();
+		m_aConst = stream->readFloat();
+		m_bConst = stream->readFloat();
+		m_cConst = stream->readFloat();
+		m_dConst = stream->readFloat();
+		m_eConst = stream->readFloat();
+		m_clipBelowHorizon = stream->readBool();
+
+		configure();
+	}
+
+	void serialize(Stream *stream, InstanceManager *manager) const {
+		Luminaire::serialize(stream, manager);
+		stream->writeFloat(m_skyScale);
+		stream->writeFloat(m_turbidity);
+		stream->writeFloat(m_thetaS);
+		stream->writeFloat(m_phiS);
+		stream->writeFloat(m_aConst);
+		stream->writeFloat(m_bConst);
+		stream->writeFloat(m_cConst);
+		stream->writeFloat(m_dConst);
+		stream->writeFloat(m_eConst);
+		stream->writeBool(m_clipBelowHorizon);
+	}
+
+	/**
+	 * Precalculates some inernal varibles. It needs a valid sun position.
+	 * Hence it needs a configureSunPos... method called before itself.
+	 */
+	void configure() {
+		Float theta2 = m_thetaS * m_thetaS;
+		Float theta3 = theta2 * m_thetaS;
+
+		Float turb2 = m_turbidity * m_turbidity;
+
+		/* calculate zenith chromaticity */
+		m_zenithX =
+		(+0.00165*theta3 - 0.00374*theta2 + 0.00208*m_thetaS + 0)       * turb2 +
+		(-0.02902*theta3 + 0.06377*theta2 - 0.03202*m_thetaS + 0.00394) * m_turbidity  +
+		(+0.11693*theta3 - 0.21196*theta2 + 0.06052*m_thetaS + 0.25885);
+
+		m_zenithY =
+		(+0.00275*theta3 - 0.00610*theta2 + 0.00316*m_thetaS  + 0)       * turb2 +
+		(-0.04214*theta3 + 0.08970*theta2 - 0.04153*m_thetaS  + 0.00515) * m_turbidity  +
+		(+0.15346*theta3 - 0.26756*theta2 + 0.06669*m_thetaS  + 0.26688);
+
+		/* calculate zenith luminance */
+		Float chi = (4.0/9.0 - m_turbidity / 120.0) * (M_PI - 2 * m_thetaS);
+
+		m_zenithL = (4.0453 * m_turbidity - 4.9710) * tan(chi)
+			- 0.2155 * m_turbidity + 2.4192;
+
+		m_perezL[0] =  ( 0.17872 * m_turbidity  - 1.46303) * m_aConst;
+		m_perezL[1] =  (-0.35540 * m_turbidity  + 0.42749) * m_bConst;
+		m_perezL[2] =  (-0.02266 * m_turbidity  + 5.32505) * m_cConst;
+		m_perezL[3] =  ( 0.12064 * m_turbidity  - 2.57705) * m_dConst;
+		m_perezL[4] =  (-0.06696 * m_turbidity  + 0.37027) * m_eConst;
+
+		m_perezX[0] =  (-0.01925 * m_turbidity - 0.25922) * m_aConst;
+		m_perezX[1] =  (-0.06651 * m_turbidity + 0.00081) * m_bConst;
+		m_perezX[2] =  (-0.00041 * m_turbidity + 0.21247) * m_cConst;
+		m_perezX[3] =  (-0.06409 * m_turbidity - 0.89887) * m_dConst;
+		m_perezX[4] =  (-0.00325 * m_turbidity + 0.04517) * m_eConst;
+
+		m_perezY[0] =  (-0.01669 * m_turbidity - 0.26078) * m_aConst;
+		m_perezY[1] =  (-0.09495 * m_turbidity + 0.00921) * m_bConst;
+		m_perezY[2] =  (-0.00792 * m_turbidity + 0.21023) * m_cConst;
+		m_perezY[3] =  (-0.04405 * m_turbidity - 1.65369) * m_dConst;
+		m_perezY[4] =  (-0.01092 * m_turbidity + 0.05291) * m_eConst;
+
+
+
+		int thetaBins = 512, phiBins = thetaBins*2;
+		ref<Bitmap> bitmap = new Bitmap(phiBins, thetaBins, 128);
+		Point2 factor(M_PI / thetaBins, (2*M_PI) / phiBins);
+		for (int i=0; i<thetaBins; ++i) {
+			Float theta = (i+.5f)*factor.x;
+			for (int j=0; j<phiBins; ++j) {
+				Float phi = (j+.5f)*factor.x;
+				Spectrum s = Le(sphericalDirection(theta, phi));
+				Float r, g, b;
+				s.toLinearRGB(r, g, b);
+				bitmap->getFloatData()[(j+i*phiBins)*4 + 0] = r;
+				bitmap->getFloatData()[(j+i*phiBins)*4 + 1] = g;
+				bitmap->getFloatData()[(j+i*phiBins)*4 + 2] = b;
+				bitmap->getFloatData()[(j+i*phiBins)*4 + 3] = 1;
+			}
+		}
+		ref<FileStream> fs = new FileStream("out.exr", FileStream::ETruncReadWrite);
+		bitmap->save(Bitmap::EEXR, fs);
+	}
+
+	/**
+	 * Configures the position of the sun. This calculation is based on
+	 * your position on the world and time of day.
+	 * From IES Lighting Handbook pg 361.
+	 */
+	void configureSunPosition(const Float lat, const Float lon,
+			const int stdMrd, const int julDay, const Float timeOfDay) {
+		const Float solarTime = timeOfDay
+			+ (0.170 * sin(4.0 * M_PI * (julDay - 80.0) / 373.0)
+			- 0.129 * sin(2.0 * M_PI * (julDay - 8.0) / 355.0))
+			+ (stdMrd - lon) / 15.0;
+
+		const Float solarDeclination = (0.4093 * sin(2 * M_PI * (julDay - 81) / 368));
+
+		const Float solarAltitude = asin(sin(degToRad(lat))
+			* sin(solarDeclination) - cos(degToRad(lat))
+			* cos(solarDeclination) * cos(M_PI * solarTime / 12.0));
+
+		const Float opp = -cos(solarDeclination) * sin(M_PI * solarTime / 12.0);
+		const Float adj = -(cos(degToRad(lat)) * sin(solarDeclination)
+			+ sin(degToRad(lat)) * cos(solarDeclination)
+			* cos(M_PI * solarTime / 12.0));
+
+		const Float solarAzimuth = atan2(opp,adj);
+
+		m_phiS = -solarAzimuth;
+		m_thetaS = M_PI / 2.0 - solarAltitude;
+	}
+
+	/**
+	 * Configures the position of the sun by using a vector that points
+	 * to the sun. It is expected, that +x = south, +y = east, +z = up.
+	 */
+	void configureSunPosition(const Vector& sunDir) {
+		Vector wh = normalize(sunDir);
+		Point2 sunPos = toSphericalCoordinates(wh);
+		m_thetaS = sunPos.x;
+		m_phiS = sunPos.y;
+	}
+
+	void preprocess(const Scene *scene) {
+		/* Get the scene's bounding sphere and slightly enlarge it */
+		m_bsphere = scene->getBSphere();
+		m_bsphere.radius *= 1.01f;
+		m_surfaceArea = m_bsphere.radius * m_bsphere.radius * M_PI;
+		m_invSurfaceArea = 1/m_surfaceArea;
+	}
+
+	Spectrum getPower() const {
+		/* TODO */
+		return m_average * (M_PI * 4 * M_PI
+		* m_bsphere.radius * m_bsphere.radius);
+	}
+
+	inline Spectrum Le(const Vector &direction) const {
+		/* Compute sky light radiance for direction */
+		Vector d = normalize(m_worldToLuminaire(direction));
+
+		if (m_clipBelowHorizon) {
+			/* if sun is below horizon, return black */
+			if (d.z < 0.0f)
+				return Spectrum(0.0f);
+		}
+
+		/* make all "zero values" the same */
+		if (d.z < 0.001f)
+			d = normalize(Vector(d.x, d.y,  0.001f));
+
+		const Point2 dSpherical = toSphericalCoordinates(d);
+		const Float theta = dSpherical.x;
+
+#ifndef ENFORCE_SINGLE_ZENITH_ANGLE
+		const Float phi = dSpherical.y;
+#else
+		Float phi;
+		if (fabs(theta) < 1e-5)
+			phi = 0.0f;
+		else
+			phi = dSpherical.y;
+#endif
+
+		Spectrum L;
+		getSkySpectralRadiance(theta, phi, L);
+		L *= m_skyScale;
+		
+		return L;
+	}
+
+	inline Spectrum Le(const Ray &ray) const {
+		return Le(normalize(ray.d));
+	}
+
+	Spectrum Le(const LuminaireSamplingRecord &lRec) const {
+		return Le(-lRec.d);
+	}
+
+	inline void sample(const Point &p, LuminaireSamplingRecord &lRec,
+			const Point2 &sample) const {
+		lRec.d = sampleDirection(sample, lRec.pdf, lRec.value);
+		lRec.sRec.p = p - lRec.d * (2 * m_bsphere.radius);
+	}
+
+	void sample(const Intersection &its, LuminaireSamplingRecord &lRec,
+		const Point2 &sample) const {
+		SkyLuminaire::sample(its.p, lRec, sample);
+	}
+
+	inline Float pdf(const Point &p, const LuminaireSamplingRecord &lRec, bool delta) const {
+#if defined(SAMPLE_UNIFORMLY)
+		return 1.0f / (4 * M_PI);
+#endif
+	}
+
+	Float pdf(const Intersection &its, const LuminaireSamplingRecord &lRec, bool delta) const {
+		return SkyLuminaire::pdf(its.p, lRec, delta);
+	}
+
+	/**
+	 * This is the tricky bit - we want to sample a ray that
+	 * has uniform density over the set of all rays passing
+	 * through the scene.
+	 * For more detail, see "Using low-discrepancy sequences and 
+	 * the Crofton formula to compute surface areas of geometric models"
+	 * by Li, X. and Wang, W. and Martin, R.R. and Bowyer, A. 
+	 * (Computer-Aided Design vol 35, #9, pp. 771--782)
+	 */
+	void sampleEmission(EmissionRecord &eRec, 
+		const Point2 &sample1, const Point2 &sample2) const {
+		Assert(eRec.type == EmissionRecord::ENormal);
+		/* Chord model - generate the ray passing through two uniformly
+		   distributed points on a sphere containing the scene */
+		Vector d = squareToSphere(sample1);
+		eRec.sRec.p = m_bsphere.center + d * m_bsphere.radius;
+		eRec.sRec.n = Normal(-d);
+		Point p2 = m_bsphere.center + squareToSphere(sample2) * m_bsphere.radius;
+		eRec.d = p2 - eRec.sRec.p;
+		Float length = eRec.d.length();
+
+		if (length == 0) {
+			eRec.value = Spectrum(0.0f);
+			eRec.pdfArea = eRec.pdfDir = 1.0f;
+			return;
+		}
+
+		eRec.d /= length;
+		eRec.pdfArea = 1.0f / (4 * M_PI * m_bsphere.radius * m_bsphere.radius);
+		eRec.pdfDir = INV_PI * dot(eRec.sRec.n, eRec.d);
+		eRec.value = Le(-eRec.d);
+	}
+
+	void sampleEmissionArea(EmissionRecord &eRec, const Point2 &sample) const {
+		if (eRec.type == EmissionRecord::ENormal) {
+			Vector d = squareToSphere(sample);
+			eRec.sRec.p = m_bsphere.center + d * m_bsphere.radius;
+			eRec.sRec.n = Normal(-d);
+			eRec.pdfArea = 1.0f / (4 * M_PI * m_bsphere.radius * m_bsphere.radius);
+			eRec.value = Spectrum(M_PI);
+		} else {
+			/* Preview mode, which is more suitable for VPL-based rendering: approximate 
+			   the infinitely far-away source with set of diffuse point sources */
+			const Float radius = m_bsphere.radius * 1.5f;
+			Vector d = squareToSphere(sample);
+			eRec.sRec.p = m_bsphere.center + d * radius;
+			eRec.sRec.n = Normal(-d);
+			eRec.pdfArea = 1.0f / (4 * M_PI * radius * radius);
+			eRec.value = Le(d) * M_PI;
+		}
+	}
+
+	Spectrum sampleEmissionDirection(EmissionRecord &eRec, const Point2 &sample) const {
+		Float radius = m_bsphere.radius;
+		if (eRec.type == EmissionRecord::EPreview) 
+			radius *= 1.5f;
+		Point p2 = m_bsphere.center + squareToSphere(sample) * radius;
+		eRec.d = p2 - eRec.sRec.p;
+		Float length = eRec.d.length();
+
+		if (length == 0.0f) {
+			eRec.pdfDir = 1.0f;
+			return Spectrum(0.0f);
+		}
+
+		eRec.d /= length;
+		eRec.pdfDir = INV_PI * dot(eRec.sRec.n, eRec.d);
+		if (eRec.type == EmissionRecord::ENormal)
+			return Le(-eRec.d) * INV_PI;
+		else
+			return Spectrum(INV_PI);
+	}
+    
+	Spectrum fDirection(const EmissionRecord &eRec) const {
+		if (eRec.type == EmissionRecord::ENormal)
+			return Le(-eRec.d) * INV_PI;
+		else
+			return Spectrum(INV_PI);
+	}
+
+	Spectrum fArea(const EmissionRecord &eRec) const {
+		Assert(eRec.type == EmissionRecord::ENormal);
+		return Spectrum(M_PI);
+	}
+
+	Spectrum f(const EmissionRecord &eRec) const {
+		if (eRec.type == EmissionRecord::ENormal)
+			return Le(-eRec.d) * INV_PI;
+		else
+			return Spectrum(INV_PI);
+	}
+
+	void pdfEmission(EmissionRecord &eRec, bool delta) const {
+		Assert(eRec.type == EmissionRecord::ENormal);
+		Float dp = dot(eRec.sRec.n, eRec.d);
+		if (dp > 0)
+			eRec.pdfDir = delta ? 0.0f : INV_PI * dp;
+		else
+			eRec.pdfDir = 0;
+		eRec.pdfArea = delta ? 0.0f : m_invSurfaceArea;
+	}
+
+	std::string toString() const {
+		std::ostringstream oss;
+		oss << "SkyLuminaire[" << std::endl
+			<< "  sky sscale = " << m_skyScale << "," << std::endl
+			<< "  power = " << getPower().toString() << "," << std::endl
+			<< "  sun pos = theta: " << m_thetaS << ", phi: "<< m_phiS << ","  << std::endl 
+			<< "  turbidity = " << m_turbidity << "," << std::endl 
+			<< "]";
+		return oss.str();
+	}
+
+	bool isBackgroundLuminaire() const {
+		return true;
+	}
+
+	Vector sampleDirection(Point2 sample, Float &pdf, Spectrum &value) const {
+#if defined(SAMPLE_UNIFORMLY)
+		pdf = 1.0f / (4*M_PI);
+		Vector d = squareToSphere(sample);
+		value = Le(-d);
+		return d;
+#endif
+	}
+
+private:
+	/**
+	 * Calculates the anlgle between two spherical cooridnates. All
+	 * angles in radians, theta angles measured from up/zenith direction.
+	 */
+	inline Float getAngleBetween(const Float thetav, const Float phiv,
+			const Float theta, const Float phi) const {
+		const Float cospsi = sin(thetav) * sin(theta) * cos(phi - phiv)
+			+ cos(thetav) * cos(theta);
+
+		if (cospsi > 1.0f)
+			return 0.0f;
+		if (cospsi < -1.0f)
+			return M_PI;
+		return acos(cospsi);
+	}
+
+	/**
+	 * Calculates the distribution of the sky radiance with two
+	 * Perez functions:
+	 *
+	 *        Perez(Theta, Gamma)
+	 *   d = ---------------------
+	 *        Perez(0, ThetaSun)
+	 *
+	 * From IES Lighting Handbook pg 361
+	 */
+	inline Float getDistribution(const Float *lam, const Float theta,
+			const Float gamma) const {
+		const Float cosGamma = cos(gamma);
+		const Float num = ( (1 + lam[0] * exp(lam[1] / cos(theta)))
+			* (1 + lam[2] * exp(lam[3] * gamma)
+			+ lam[4] * cosGamma * cosGamma));
+
+		const Float cosTheta = cos(m_thetaS);
+		const Float den = ( (1 + lam[0] * exp(lam[1] /* / cos 0 */))
+			* (1 + lam[2] * exp(lam[3] * m_thetaS)
+			+ lam[4] * cosTheta * cosTheta));
+
+		return (num / den);
+	}
+
+	/**
+	 * Calculates the spectral radiance of the sky in the specified directiono.
+	 */
+	void getSkySpectralRadiance(const Float theta, const Float phi, Spectrum &dstSpect) const {
+		/* add bottom half of hemisphere with horizon colour */
+		const Float theta_fin = std::min(theta, (M_PI * 0.5f) - 0.001f);
+		/* get angle between sun (zenith is 0, 0) and point (theta, phi) */
+		const Float gamma = getAngleBetween(theta, phi, m_thetaS, m_phiS);
+		/* Compute xyY values by calculating the distribution for the point
+		 * point of interest and multiplying it with the the components
+		 * zenith value. */
+		const Float x = m_zenithX * getDistribution(m_perezX, theta_fin, gamma);
+		const Float y = m_zenithY * getDistribution(m_perezY, theta_fin, gamma);
+		const Float Y = m_zenithL * getDistribution(m_perezL, theta_fin, gamma);
+
+		/* Apply an exponential exposure function */
+		// Y = 1.0 - exp(-m_exposure * Y);
+		/* Convert xyY to XYZ */
+		const Float yFrac = Y / y;
+		const Float X = yFrac * x;
+		/* It seems the following is necassary to stay always above zero */
+		const Float z = std::max(0.0f, 1.0f - x - y);
+		const Float Z = yFrac * z;
+
+		/* Create spectrum from XYZ values */
+		dstSpect.fromXYZ(X, Y, Z);
+		/* The produced spectrum contains out-of-gamut colors.
+		 * It is common to clamp resulting values to zero. */
+		dstSpect.clampNegative();
+	}
+
+	MTS_DECLARE_CLASS()
+protected:
+	Spectrum m_average;
+	Float m_surfaceArea;
+	Float m_invSurfaceArea;
+	BSphere m_bsphere;
+	Float m_exposure;
+	Float m_skyScale;
+	/* The turbidity of the sky ranges normally from 0 to 30+.
+	 * For clear skies values in range [2,6] are useful. */
+	Float m_turbidity;
+	Float m_thetaS, m_phiS;
+	Float m_zenithL;
+	Float m_zenithX;
+	Float m_zenithY;
+	/* The distribution coefficints are called A, B, C, D and E by
+	 * Preedham. Since they exist for x, y and Y (here called L)
+	 * we save the precalculated version of it
+	 *
+	 * distribution coefficients for luminance distribution function */
+	Float m_perezL[5];
+	/* distribution coefficients for x distribution function */
+	Float m_perezX[5];
+	/* distribution coefficients for y distribution function */
+	Float m_perezY[5];
+	/* distribution function tuning coefficients a, b, c, d and e.
+	 * They can be used to scale the luminance, x and y distribution
+	 * coefficints. */
+	Float m_aConst, m_bConst, m_cConst, m_dConst, m_eConst;
+
+	/* To disable clipping of incoming sky light below the horizon, set this
+	 * to false. If set to true black is returned for queries. Be awere that
+	 * a huge amount of additional radiance is coming in (i.e. it gets a
+	 * lot brighter). */
+	bool m_clipBelowHorizon;
+};
+
+MTS_IMPLEMENT_CLASS_S(SkyLuminaire, false, Luminaire)
+MTS_EXPORT_PLUGIN(SkyLuminaire, "Sky luminaire");
+MTS_NAMESPACE_END
+

src/phase/SConscript

@@ -2,6 +2,7 @@ Import('env', 'plugins')
 
 plugins += env.SharedLibrary('isotropic', ['isotropic.cpp'])
 plugins += env.SharedLibrary('hg', ['hg.cpp'])
+plugins += env.SharedLibrary('rayleigh', ['rayleigh.cpp'])
 plugins += env.SharedLibrary('kkay', ['kkay.cpp'])
 plugins += env.SharedLibrary('microflake', ['microflake.cpp'])
 

src/phase/hg.cpp

@@ -23,9 +23,23 @@
 
 MTS_NAMESPACE_BEGIN
 
-/**
+/*!\plugin{hg}{Henyey-Greenstein phase function}
- * Phase function by Henyey and Greenstein (1941). Parameterizable
+ * \order{2}
- * from backward- through isotropic- to forward scattering.
+ * \parameters{
+ *     \parameter{g}{\Float}{
+ *       This parameter must be somewhere in the range $-1$ to $1$
+ *       (but not equal to $-1$ or $1$). It denotes the \emph{mean cosine} 
+ *       of scattering interactions. A value greater than zero indicates that
+ *       medium interactions predominantly scatter incident light into a similar
+ *       direction (i.e. the medium is \emph{forward-scattering}), whereas
+ *       values smaller than zero cause the medium to be
+ *       scatter more light in the opposite direction.
+ *     }
+ * }
+ * This plugin implements the phase function model proposed by 
+ * Henyey and Greenstein \cite{Henyey1941Diffuse}. It is 
+ * parameterizable from backward- ($g<0$) through 
+ * isotropic- ($g=0$) to forward ($g>0$) scattering.
  */
 class HGPhaseFunction : public PhaseFunction {
 public:
@@ -35,6 +49,8 @@ public:
 		   lie in [-1, 1] where >0 is forward scattering and <0 is backward
 		   scattering. */
 		m_g = props.getFloat("g", 0.8f);
+		if (m_g >= 1 || m_g <= -1)
+			Log(EError, "Asymmetry parameter must be in the interval (-1, 1)!");
 	}
 
 	HGPhaseFunction(Stream *stream, InstanceManager *manager) 

src/phase/isotropic.cpp

@@ -21,8 +21,12 @@
 
 MTS_NAMESPACE_BEGIN
 
-/**
+/*!\plugin{isotropic}{Isotropic phase function}
- * Isotropic phase function (i.e. f=1/(4*pi))
+ * \order{1}
+ *
+ * This phase function simulates completely uniform scattering,
+ * where all directionality is lost after a single scattering 
+ * interaction. It does not have any parameters.
  */
 class IsotropicPhaseFunction : public PhaseFunction {
 public:

src/phase/kkay.cpp

@@ -19,14 +19,15 @@
 #include <mitsuba/render/phase.h>
 #include <mitsuba/render/medium.h>
 #include <mitsuba/render/sampler.h>
-#include <mitsuba/core/random.h>
 #include <mitsuba/core/properties.h>
 #include <mitsuba/core/frame.h>
 
 MTS_NAMESPACE_BEGIN
 
-/**
+/*!\plugin{kkay}{Kajiya-Kay phase function}
- * \brief Kajiya-Kay phase function
+ * This plugin implements the Kajiya-Kay \cite{Kajiya1989Rendering}
+ * phase function for volumetric rendering of fibers, e.g. 
+ * hair or cloth.
  *
  * The function is normalized so that it has no energy loss when 
  * \a ks=1 and illumination arrives perpendicularly to the surface.
@@ -118,41 +119,6 @@ public:
 		return "KajiyaKayPhaseFunction[]";
 	}
 
-#if 0
-	/// For testing purposes
-	Float integrateOverOutgoing(const Vector &wi) {
-		MediumSamplingRecord mRec;
-		mRec.orientation = Vector(0,0,1);
-		int res = 100;
-		/* Nested composite Simpson's rule */
-		Float hExt = M_PI / res,
-		      hInt = (2*M_PI)/(res*2);
-
-		Float result = 0;
-
-		for (int i=0; i<=res; ++i) {
-			Float theta = hExt*i;
-			Float weightExt = (i & 1) ? 4.0f : 2.0f;
-			if (i == 0 || i == res)
-				weightExt = 1;
-
-			for (int j=0; j<=res*2; ++j) {
-				Float phi = hInt*j;
-				Float weightInt = (j & 1) ? 4.0f : 2.0f;
-				if (j == 0 || j == 2*res)
-					weightInt = 1;
-
-				Float value = f(mRec, wi, sphericalDirection(theta, phi))[0]*std::sin(theta)
-					* weightExt*weightInt;
-
-				result += value;
-			}
-		}
-		
-		return hExt*hInt/9;
-	}
-#endif
-
 	MTS_DECLARE_CLASS()
 private:
 	Float m_ks, m_kd, m_exponent, m_normalization;

src/phase/microflake.cpp

@@ -24,6 +24,8 @@
 /// Generate a few statistics related to the implementation?
 // #define MICROFLAKE_STATISTICS 1
 
+/// The following file implements the micro-flake distribution
+/// for rough fibers
 #include "microflake_fiber.h"
 
 MTS_NAMESPACE_BEGIN
@@ -33,24 +35,32 @@ static StatsCounter avgSampleIterations("Micro-flake model",
 		"Average rejection sampling iterations", EAverage);
 #endif
 
-/**
+/*!\plugin{microflake}{Micro-flake phase function}
- * Implements the anisotropic micro-flake phase function described in
+ * \parameters{
+ *     \parameter{stddev}{\Float}{
+ *       Standard deviation of the micro-flake normals. This
+ *       specifies the roughness of the fibers in the medium. 
+ *     }
+ * }
+ * This plugin implements the anisotropic micro-flake phase function 
+ * described in
+ * ``A radiative transfer framework for rendering materials with
+ * anisotropic structure'' by Wenzel Jakob, Adam Arbree, 
+ * Jonathan T. Moon, Kavita Bala, and Steve Marschner 
+ * \cite{Jakob2010Radiative}.
  *
- * "A radiative transfer framework for rendering materials with
+ * The implementation in this plugin is specific to rough fibers
- * anisotropic structure" by Wenzel Jakob, Adam Arbree, 
+ * and uses a Gaussian-type flake distribution. It is much faster
- * Jonathan T. Moon, Kavita Bala, and Steve Marschner,
+ * than the spherical harmonics approach proposed in the original 
- * ACM SIGGRAPH 2010
+ * paper. This distribution, as well as the implemented sampling 
+ * method, are described in the paper
+ * ``Building Volumetric Appearance Models of Fabric using 
+ * Micro CT Imaging'' by Shuang Zhao, Wenzel Jakob, Steve Marschner,
+ * and Kavita Bala \cite{Zhao2011Building}.
  *
- * The optimized implementation here works without the use of 
+ * Note: this phase function must be used with a medium that specifies
- * spherical harmonics and is specific to rough fibers and a 
+ * the local fiber orientation at different points in space. Please
- * Gaussian-type distribution. This distribution, as well as the
+ * refer to \pluginref{heterogeneous} for details.
- * implemented sampling method are described in the paper
- * 
- * "Building Volumetric Appearance Models of Fabric using 
- * Micro CT Imaging" by Shuang Zhao, Wenzel Jakob, Steve Marschner,
- * and Kavita Bala, ACM SIGGRAPH 2011
- *
- * \author Wenzel Jakob
  */
 class MicroflakePhaseFunction : public PhaseFunction {
 public:

src/phase/rayleigh.cpp

@@ -0,0 +1,99 @@
+/*
+    This file is part of Mitsuba, a physically based rendering system.
+
+    Copyright (c) 2007-2011 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/phase.h>
+#include <mitsuba/render/sampler.h>
+#include <mitsuba/core/properties.h>
+#include <mitsuba/core/frame.h>
+
+MTS_NAMESPACE_BEGIN
+
+/*!\plugin{rayleigh}{Rayleigh phase function}
+ * \order{3}
+ *
+ * Scattering by particles that are much smaller than the wavelength
+ * of light (e.g. in aerosols) is well-approxiamted by the Rayleigh
+ * phase function. This plugin implements an unpolarized version of this 
+ * scattering model (i.e the effects of polarization are ignored).
+ * This plugin is useful for simulating scattering in planetary 
+ * atmospheres.
+ *
+ * This model has no parameters.
+ */
+class RayleighPhaseFunction : public PhaseFunction {
+public:
+	RayleighPhaseFunction(const Properties &props) 
+		: PhaseFunction(props) { }
+
+	RayleighPhaseFunction(Stream *stream, InstanceManager *manager) 
+		: PhaseFunction(stream, manager) { }
+
+	virtual ~RayleighPhaseFunction() { }
+
+	void serialize(Stream *stream, InstanceManager *manager) const {
+		PhaseFunction::serialize(stream, manager);
+	}
+
+	inline Float sample(PhaseFunctionQueryRecord &pRec,
+			Sampler *sampler) const {
+		Point2 sample(sampler->next2D());
+
+		Float z = 2 * (2*sample.x - 1),
+			  tmp = std::sqrt(z*z+1),
+			  A = std::pow(z+tmp, 1.0f/3.0f),
+			  B = std::pow(z-tmp, 1.0f/3.0f),
+			  cosTheta = A + B,
+			  sinTheta = std::sqrt(std::max((Float) 0, 1.0f-cosTheta*cosTheta)),
+			  phi = 2*M_PI*sample.y,
+			  cosPhi = std::cos(phi),
+			  sinPhi = std::sin(phi);
+
+		Vector dir(
+			sinTheta * cosPhi,
+			sinTheta * sinPhi,
+			cosTheta);
+
+		pRec.wo = Frame(-pRec.wi).toWorld(dir);
+		return 1.0f;
+	}
+
+	Float sample(PhaseFunctionQueryRecord &pRec,
+			Float &pdf, Sampler *sampler) const {
+		RayleighPhaseFunction::sample(pRec, sampler);
+		pdf = RayleighPhaseFunction::eval(pRec);
+		return pdf;
+	}
+
+
+	Float eval(const PhaseFunctionQueryRecord &pRec) const {
+		Float mu = dot(pRec.wi, pRec.wo);
+		return (3.0f/(16.0f*M_PI)) * (1+mu*mu);
+	}
+
+	std::string toString() const {
+		std::ostringstream oss;
+		oss << "RayleighPhaseFunction[]";
+		return oss.str();
+	}
+
+	MTS_DECLARE_CLASS()
+};
+
+MTS_IMPLEMENT_CLASS_S(RayleighPhaseFunction, false, PhaseFunction)
+MTS_EXPORT_PLUGIN(RayleighPhaseFunction, "Rayleigh phase function");
+MTS_NAMESPACE_END