Implementation of the paper 'Accurate computation of single scattering in participating media with refractive boundaries' contributed by Nicolas Holzschuch

data/schema/scene.xsd

@@ -133,7 +133,7 @@
 				<xsd:choice minOccurs="0" maxOccurs="unbounded">
 					<xsd:group ref="objectGroup"/>
 					<xsd:element name="bsdf" type="bsdf"/>
-					<xsd:element name="subsurface" type="object"/>
+                    <xsd:element name="subsurface" type="subsurface"/>
 					<xsd:element name="ref" type="reference"/>
 					<xsd:element name="sensor" type="sensor"/>
 					<xsd:element name="emitter" type="emitter"/>
@@ -146,6 +146,19 @@
 		</xsd:complexContent>
 	</xsd:complexType>
 
+    <!-- SINGLESCATTER Element -->
+    <xsd:complexType name="subsurface">
+        <xsd:complexContent>
+            <xsd:extension base="objectBase">
+                <xsd:choice minOccurs="0" maxOccurs="unbounded">
+                    <xsd:group ref="objectGroup"/>
+                    <xsd:element name="phase" type="phase"/>
+                    <xsd:element name="bsdf" type="bsdf"/>
+                </xsd:choice>
+            </xsd:extension>
+        </xsd:complexContent>
+    </xsd:complexType>
+
 	<!-- MEDIUM Element -->
 	<xsd:complexType name="medium">
 		<xsd:complexContent>

include/mitsuba/core/matrix.h

@@ -607,6 +607,13 @@ public:
 			m[2][0] * v.x + m[2][1] * v.y + m[2][2] * v.z);
 	}
 
+	/// multiply the matrix by transpose of v, returns the transpose of the line vector.
+	inline Vector preMult(const Vector &v) const {
+		return Vector(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0],
+					  v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1],
+					  v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2]);
+	}
+
 	/// Scalar multiplication (creates a temporary)
 	inline Matrix3x3 operator*(Float value) const {
 		Matrix3x3 result;
@@ -616,6 +623,7 @@ public:
 		return result;
 	}
 
+
 	/// Assignment operator
 	inline Matrix3x3 &operator=(const Matrix<3, 3, Float> &mat) {
 		for (int i=0; i<3; ++i)

include/mitsuba/core/triangle.h

@@ -58,6 +58,21 @@ struct MTS_EXPORT_CORE Triangle {
 	 */
 	AABB getClippedAABB(const Point *positions, const AABB &aabb) const;
 
+	// Returns the bounding sphere of the triangle
+	inline BSphere getBSphere(const Point *positions) const {
+		Vector a = (positions[idx[1]] - positions[idx[0]]);
+		Vector b = (positions[idx[2]] - positions[idx[0]]);
+		Float a2 = dot(a, a);
+		Float b2 = dot(b, b);
+		Float da = std::sqrt(a2);
+		Float db = std::sqrt(b2);
+		Vector axb = cross(a, b);
+		Float axb2 = dot(axb, axb);
+		Float daxb = std::sqrt(axb2);
+		return BSphere(positions[idx[0]] + cross(a2 * b - b2 * a, axb) / (2 * axb2),
+					   da * db * (a - b).length() / (2 * daxb));
+	}
+
 	/// Uniformly sample a point on the triangle and return its normal and UV coordinates
 	Point sample(const Point *positions, const Normal *normals,
 			const Point2 *texCoords, Normal &n, Point2 &uv,

include/mitsuba/render/gkdtree.h

@@ -762,6 +762,11 @@ public:
 		m_traversalCost = traversalCost;
 	}
 
+	/**
+	 * \brief Returns the underlying kd-tree index buffer
+	 */
+	inline IndexType *getIndices() const { return m_indices; }
+
 	/**
 	 * \brief Return the traversal cost used by the tree construction heuristic
 	 */

include/mitsuba/render/skdtree.h

@@ -71,6 +71,8 @@ class MTS_EXPORT_RENDER ShapeKDTree : public SAHKDTree3D<ShapeKDTree> {
 	friend class SAHKDTree3D<ShapeKDTree>;
 	friend class Instance;
 	friend class AnimatedInstance;
+	friend class SingleScatter;
+
 public:
 	// =============================================================
 	//! @{ \name Initialization and tree construction

src/subsurface/SConscript

@@ -3,4 +3,5 @@ Import('env', 'plugins')
 plugins += env.SharedLibrary('dipole', 
 	['dipole.cpp', 'irrproc.cpp', 'irrtree.cpp', 'bluenoise.cpp'])
 
+plugins += env.SharedLibrary('singlescatter', ['singlescatter.cpp'])
 Export('plugins')