diff --git a/include/mitsuba/core/util.h b/include/mitsuba/core/util.h
index 46089548..c31e6df7 100644
--- a/include/mitsuba/core/util.h
+++ b/include/mitsuba/core/util.h
@@ -358,11 +358,14 @@ extern MTS_EXPORT_CORE bool solveQuadraticDouble(double a, double b,
  *      Position of the last knot
  * \param size
  *      Denotes the size of the \c data array
+ * \param extrapolate
+ *      Extrapolate data values when \c x is out of range? (default: \c false)
  * \return
- *      The interpolated value or zero when \c x lies outside of [\c min, \c max]
+ *      The interpolated value or zero when <tt>extrapolate=false</tt>tt>
+ *      and \c x lies outside of [\c min, \c max]
  */
 extern MTS_EXPORT_CORE Float interpCubic1D(Float x, const Float *data,
-		Float min, Float max, size_t size);
+		Float min, Float max, size_t size, bool extrapolate = false);
 
 /**
  * \brief Evaluate a cubic spline interpolant of an \a irregularly sampled 1D function
@@ -382,11 +385,14 @@ extern MTS_EXPORT_CORE Float interpCubic1D(Float x, const Float *data,
  *      the entries of \c nodes.
  * \param size
  *      Denotes the size of the \c data array
+ * \param extrapolate
+ *      Extrapolate data values when \c x is out of range? (default: \c false)
  * \return
- *      The interpolated value or zero when \c x lies outside of \a [\c min, \c max]
+ *      The interpolated value or zero when <tt>extrapolate=false</tt>tt>
+ *      and \c x lies outside of \a [\c min, \c max]
  */
 extern MTS_EXPORT Float interpCubic1DIrregular(Float x, const Float *nodes,
-		const Float *data, size_t size);
+		const Float *data, size_t size, bool extrapolate = false);
 
 /**
  * \brief Evaluate a cubic spline interpolant of a regularly sampled 2D function
@@ -407,11 +413,14 @@ extern MTS_EXPORT Float interpCubic1DIrregular(Float x, const Float *nodes,
  *      Position of the last knot on each dimension
  * \param size
  *      Denotes the size of the \c data array (along each dimension)
+ * \param extrapolate
+ *      Extrapolate data values when \c p is out of range? (default: \c false)
  * \return
- *      The interpolated value or zero when \c p lies outside of the knot range
+ *      The interpolated value or zero when <tt>extrapolate=false</tt>tt> and
+ *      \c p lies outside of the knot range
  */
 extern MTS_EXPORT_CORE Float interpCubic2D(const Point2 &p, const Float *data,
-		const Point2 &min, const Point2 &max, const Size2 &size);
+		const Point2 &min, const Point2 &max, const Size2 &size, bool extrapolate = false);
 
 /**
  * \brief Evaluate a cubic spline interpolant of an \a irregularly sampled 2D function
@@ -435,11 +444,14 @@ extern MTS_EXPORT_CORE Float interpCubic2D(const Point2 &p, const Float *data,
  *      Consecutive entries of this array correspond to increments in the 'x' coordinate.
  * \param size
  *      Denotes the size of the \c data array (along each dimension)
+ * \param extrapolate
+ *      Extrapolate data values when \c p is out of range? (default: \c false)
  * \return
- *      The interpolated value or zero when \c p lies outside of the knot range
+ *      The interpolated value or zero when <tt>extrapolate=false</tt>tt> and
+ *      \c p lies outside of the knot range
  */
 extern MTS_EXPORT_CORE Float interpCubic2DIrregular(const Point2 &p, const Float **nodes,
-		const Float *data, const Size2 &size);
+		const Float *data, const Size2 &size, bool extrapolate = false);
 
 /**
  * \brief Evaluate a cubic spline interpolant of a regularly sampled 3D function
@@ -461,11 +473,14 @@ extern MTS_EXPORT_CORE Float interpCubic2DIrregular(const Point2 &p, const Float
  *      Position of the last knot on each dimension
  * \param size
  *      Denotes the size of the \c data array (along each dimension)
+ * \param extrapolate
+ *      Extrapolate data values when \c p is out of range? (default: \c false)
  * \return
- *      The interpolated value or zero when \c p lies outside of the knot range
+ *      The interpolated value or zero when <tt>extrapolate=false</tt>tt> and
+ *      \c p lies outside of the knot range
  */
 extern MTS_EXPORT_CORE Float interpCubic3D(const Point3 &p, const Float *data,
-		const Point3 &min, const Point3 &max, const Size3 &size);
+		const Point3 &min, const Point3 &max, const Size3 &size, bool extrapolate = false);
 
 /**
  * \brief Evaluate a cubic spline interpolant of an \a irregularly sampled 3D function
@@ -490,11 +505,14 @@ extern MTS_EXPORT_CORE Float interpCubic3D(const Point3 &p, const Float *data,
  *      then 'y', and finally 'z' increments.
  * \param size
  *      Denotes the size of the \c data array (along each dimension)
+ * \param extrapolate
+ *      Extrapolate data values when \c p is out of range? (default: \c false)
  * \return
- *      The interpolated value or zero when \c p lies outside of the knot range
+ *      The interpolated value or zero when <tt>extrapolate=false</tt>tt> and
+ *      \c p lies outside of the knot range
  */
 extern MTS_EXPORT_CORE Float interpCubic3DIrregular(const Point3 &p, const Float **nodes,
-		const Float *data, const Size3 &size);
+		const Float *data, const Size3 &size, bool extrapolate = false);
 
 //// Convert radians to degrees
 inline Float radToDeg(Float value) { return value * (180.0f / M_PI); }
diff --git a/include/mitsuba/render/trimesh.h b/include/mitsuba/render/trimesh.h
index c4d62365..a68535cd 100644
--- a/include/mitsuba/render/trimesh.h
+++ b/include/mitsuba/render/trimesh.h
@@ -206,8 +206,15 @@ public:
 	 */
 	void computeUVTangents();
 
-	/// Generate surface normals
-	void computeNormals();
+	/**
+	 * \brief Generate smooth vertex normals?
+	 *
+	 * \param force
+	 *   When this parameter is set to true, the function
+	 *   generates normals <em>even</em> when there are
+	 *   already existing ones.
+	 */
+	void computeNormals(bool force = false);
 
 	/**
 	 * \brief Rebuild the mesh so that adjacent faces
diff --git a/src/libcore/util.cpp b/src/libcore/util.cpp
index 3f2bf951..fcadfb84 100644
--- a/src/libcore/util.cpp
+++ b/src/libcore/util.cpp
@@ -469,9 +469,9 @@ bool solveLinearSystem2x2(const Float a[2][2], const Float b[2], Float x[2]) {
 	return true;
 }
 
-Float interpCubic1D(Float x, const Float *data, Float min, Float max, size_t size) {
+Float interpCubic1D(Float x, const Float *data, Float min, Float max, size_t size, bool extrapolate) {
 	/* Give up when given an out-of-range or NaN argument */
-	if (!(x >= min && x <= max))
+	if (!(x >= min && x <= max) && !extrapolate)
 		return 0.0f;
 
 	/* Transform 'x' so that knots lie at integer positions */
@@ -508,9 +508,9 @@ Float interpCubic1D(Float x, const Float *data, Float min, Float max, size_t siz
 		(   t3 - t2)       * d1;
 }
 
-Float interpCubic1DIrregular(Float x, const Float *nodes, const Float *data, size_t size) {
+Float interpCubic1DIrregular(Float x, const Float *nodes, const Float *data, size_t size, bool extrapolate) {
 	/* Give up when given an out-of-range or NaN argument */
-	if (!(x >= nodes[0] && x <= nodes[size-1]))
+	if (!(x >= nodes[0] && x <= nodes[size-1]) && !extrapolate)
 		return 0.0f;
 
 	size_t k = (size_t) std::max((ptrdiff_t) 0, std::min((ptrdiff_t) size - 2,
@@ -545,7 +545,7 @@ Float interpCubic1DIrregular(Float x, const Float *nodes, const Float *data, siz
 
 
 Float interpCubic2D(const Point2 &p, const Float *data,
-		const Point2 &min, const Point2 &max, const Size2 &size) {
+		const Point2 &min, const Point2 &max, const Size2 &size, bool extrapolate) {
 	Float knotWeights[2][4];
 	Size2 knot;
 
@@ -553,7 +553,7 @@ Float interpCubic2D(const Point2 &p, const Float *data,
 	for (int dim=0; dim<2; ++dim) {
 		Float *weights = knotWeights[dim];
 		/* Give up when given an out-of-range or NaN argument */
-		if (!(p[dim] >= min[dim] && p[dim] <= max[dim]))
+		if (!(p[dim] >= min[dim] && p[dim] <= max[dim]) && !extrapolate)
 			return 0.0f;
 
 		/* Transform 'p' so that knots lie at integer positions */
@@ -615,7 +615,7 @@ Float interpCubic2D(const Point2 &p, const Float *data,
 }
 
 Float interpCubic2DIrregular(const Point2 &p, const Float **nodes_,
-			const Float *data, const Size2 &size) {
+			const Float *data, const Size2 &size, bool extrapolate) {
 	Float knotWeights[2][4];
 	Size2 knot;
 
@@ -625,7 +625,7 @@ Float interpCubic2DIrregular(const Point2 &p, const Float **nodes_,
 		Float *weights = knotWeights[dim];
 
 		/* Give up when given an out-of-range or NaN argument */
-		if (!(p[dim] >= nodes[0] && p[dim] <= nodes[size[dim]-1]))
+		if (!(p[dim] >= nodes[0] && p[dim] <= nodes[size[dim]-1]) && !extrapolate)
 			return 0.0f;
 
 		/* Find the index of the left knot in the queried subinterval, be
@@ -689,7 +689,7 @@ Float interpCubic2DIrregular(const Point2 &p, const Float **nodes_,
 }
 
 Float interpCubic3D(const Point3 &p, const Float *data,
-		const Point3 &min, const Point3 &max, const Size3 &size) {
+		const Point3 &min, const Point3 &max, const Size3 &size, bool extrapolate) {
 	Float knotWeights[3][4];
 	Size3 knot;
 
@@ -697,7 +697,7 @@ Float interpCubic3D(const Point3 &p, const Float *data,
 	for (int dim=0; dim<3; ++dim) {
 		Float *weights = knotWeights[dim];
 		/* Give up when given an out-of-range or NaN argument */
-		if (!(p[dim] >= min[dim] && p[dim] <= max[dim]))
+		if (!(p[dim] >= min[dim] && p[dim] <= max[dim]) && !extrapolate)
 			return 0.0f;
 
 		/* Transform 'p' so that knots lie at integer positions */
@@ -763,7 +763,7 @@ Float interpCubic3D(const Point3 &p, const Float *data,
 }
 
 Float interpCubic3DIrregular(const Point3 &p, const Float **nodes_,
-			const Float *data, const Size3 &size) {
+			const Float *data, const Size3 &size, bool extrapolate) {
 	Float knotWeights[3][4];
 	Size3 knot;
 
@@ -773,7 +773,7 @@ Float interpCubic3DIrregular(const Point3 &p, const Float **nodes_,
 		Float *weights = knotWeights[dim];
 
 		/* Give up when given an out-of-range or NaN argument */
-		if (!(p[dim] >= nodes[0] && p[dim] <= nodes[size[dim]-1]))
+		if (!(p[dim] >= nodes[0] && p[dim] <= nodes[size[dim]-1]) && !extrapolate)
 			return 0.0f;
 
 		/* Find the index of the left knot in the queried subinterval, be
diff --git a/src/librender/trimesh.cpp b/src/librender/trimesh.cpp
index 56b747c4..98e7f757 100644
--- a/src/librender/trimesh.cpp
+++ b/src/librender/trimesh.cpp
@@ -542,7 +542,7 @@ void TriMesh::rebuildTopology(Float maxAngle) {
 	configure();
 }
 
-void TriMesh::computeNormals() {
+void TriMesh::computeNormals(bool force) {
 	int invalidNormals = 0;
 	if (m_faceNormals) {
 		if (m_normals) {
@@ -558,7 +558,7 @@ void TriMesh::computeNormals() {
 			}
 		}
 	} else {
-		if (m_normals) {
+		if (m_normals && !force) {
 			if (m_flipNormals) {
 				for (size_t i=0; i<m_vertexCount; i++)
 					m_normals[i] *= -1;
@@ -566,7 +566,8 @@ void TriMesh::computeNormals() {
 				/* Do nothing */
 			}
 		} else {
-			m_normals = new Normal[m_vertexCount];
+			if (!m_normals)
+				m_normals = new Normal[m_vertexCount];
 			memset(m_normals, 0, sizeof(Normal)*m_vertexCount);
 
 			/* Well-behaved vertex normal computation based on