hslt: more consistent shading frame computation in various shapes

include/mitsuba/core/util.h

@@ -329,13 +329,42 @@ extern MTS_EXPORT_CORE bool solveLinearSystem2x2(const Float a[2][2], const Floa
  */
 extern MTS_EXPORT_CORE void coordinateSystem(const Vector &a, Vector &b, Vector &c);
 
+
 /**
- * \brief Derivatives of a frame formed by coordinateSystem
+ * \brief Given a smoothly varying shading normal and a tangent of a shape parameterization,
- * \param n Source tangent frame that was created with coordinateSystem
+ * compute a smoothly varying orthonormal frame
- * \param ds derivative of the frame with respect to s. ds.n should already contain normal derivative along s
+ *
- * \param dt derivative of the frame with respect to t. dt.n should already contain normal derivative along t
+ * \param n
+ *    A shading normal at a surface position
+ * \param dpdu
+ *    Position derivative of the underlying parameterization with respect to the 'u' coordinate
+ * \param frame
+ *    Used to return the computed frame
+ *
+ * Mitsuba uses this function to compute the field \ref Intersection::shFrame
  */
-extern MTS_EXPORT_CORE void coordinateSystemDerivatives(const Frame &frame, Frame &ds, Frame &dt);
+extern MTS_EXPORT_CORE void computeShadingFrame(const Vector &n, const Vector &dpdu, Frame &frame);
+
+/**
+ * \brief Compute the spatial derivative of \ref computeShadingFrame
+ *
+ * This is used by Manifold Exploration and Half Vector Light Transport
+ *
+ * \param n
+ *    A shading normal at a surface position
+ * \param dpdu
+ *    Position derivative of the underlying parameterization with respect to the 'u' coordinate
+ * \param dndu
+ *    Derivative of the shading normal along the 'u' coordinate
+ * \param dndv
+ *    Derivative of the shading normal along the 'v' coordinate
+ * \param du
+ *    Used to return the 'u' derivative of the frame
+ * \param dv
+ *    Used to return the 'v' derivative of the frame
+ */
+extern MTS_EXPORT_CORE void computeShadingFrameDerivative(const Vector &n, const Vector &dpdu, 
+		const Vector &dndu, const Vector &dndv, Frame &du, Frame &dv);
 
 /**
  * \brief Generate (optionally jittered) stratified 1D samples

include/mitsuba/render/skdtree.h

@@ -376,6 +376,7 @@ protected:
 				its.dpdu = side1;
 				its.dpdv = side2;
 			}
+
 			if (EXPECT_TAKEN(vertexNormals)) {
 				const Normal
 					&n0 = vertexNormals[idx0],
@@ -384,23 +385,13 @@ protected:
 
 				its.shFrame.n = normalize(n0 * b.x + n1 * b.y + n2 * b.z);
 
-				if (EXPECT_TAKEN(!vertexTangents)) {
-					coordinateSystem(its.shFrame.n, its.shFrame.s, its.shFrame.t);
-				} else {
-					/* Align shFrame.s with dpdu, use Gram-Schmidt to orthogonalize */
-					its.shFrame.s = normalize(its.dpdu - its.shFrame.n
-						* dot(its.shFrame.n, its.dpdu));
-					its.shFrame.t = cross(its.shFrame.n, its.shFrame.s);
-				}
-
 				/* Ensure that the geometric & shading normals face the same direction */
 				if (dot(faceNormal, its.shFrame.n) < 0)
 					faceNormal = -faceNormal;
-
-				its.geoFrame = Frame(faceNormal);
 			} else {
-				its.shFrame = its.geoFrame = Frame(faceNormal);
+				its.shFrame.n = faceNormal;
 			}
+			its.geoFrame = Frame(faceNormal);
 
 			if (EXPECT_TAKEN(vertexTexcoords)) {
 				const Point2 &t0 = vertexTexcoords[idx0];
@@ -420,7 +411,6 @@ protected:
 					result[2], Spectrum::EReflectance);
 			}
 
-			its.wi = its.toLocal(-ray.d);
 			its.shape = trimesh;
 			its.hasUVPartials = false;
 			its.primIndex = cache->primIndex;
@@ -430,6 +420,9 @@ protected:
 			shape->fillIntersectionRecord(ray,
 				reinterpret_cast<const uint8_t*>(temp) + 2*sizeof(IndexType), its);
 		}
+
+		computeShadingFrame(its.shFrame.n, its.dpdu, its.shFrame);
+		its.wi = its.toLocal(-ray.d);
 	}
 
 	/// Plain shadow ray query (used by the 'instance' plugin)

src/libcore/util.cpp

@@ -600,27 +600,29 @@ void coordinateSystem(const Vector &a, Vector &b, Vector &c) {
 	b = cross(c, a);
 }
 
-void coordinateSystemDerivatives(const Frame &frame, Frame &ds, Frame &dt) {
+void computeShadingFrame(const Vector &n, const Vector &dpdu, Frame &frame) {
-
+	frame.n = n;
-	const Vector n = frame.n;
+	frame.s = normalize(dpdu - frame.n
-	const Vector s = frame.s;
+		* dot(frame.n, dpdu));
-
+	frame.t = cross(frame.n, frame.s);
-	if(std::abs(n.x) > std::abs(n.y)) {
-		const Float invLen = 1 / std::sqrt(n.x * n.x + n.z * n.z);
-		ds.s = Vector(ds.n.z * invLen, 0, -ds.n.x * invLen);
-		ds.s -= s * dot(ds.s, s);
-		dt.s = Vector(dt.n.z * invLen, 0, -dt.n.x * invLen);
-		dt.s -= s * dot(dt.s, s);
-	} else {
-		const Float invLen = 1 / std::sqrt(n.y * n.y + n.z * n.z);
-		ds.s = Vector(0, ds.n.z * invLen, -ds.n.y * invLen);
-		ds.s -= s * dot(ds.s, s);
-		dt.s = Vector(0, dt.n.z * invLen, -dt.n.y * invLen);
-		dt.s -= s * dot(dt.s, s);
 }
 
-	dt.t = cross(s, ds.n) + cross(dt.s, n);
+void computeShadingFrameDerivative(const Vector &n, const Vector &dpdu, const Vector &dndu, const Vector &dndv, Frame &du, Frame &dv) {
-	ds.t = cross(s, dt.n) + cross(ds.s, n);
+	Vector s = dpdu - n * dot(n, dpdu);
+	Float invLen_s = 1.0f / s.length();
+	s *= invLen_s;
+
+	du.s = invLen_s * (-dndu * dot(n, dpdu) - n * dot(dndu, dpdu));
+	dv.s = invLen_s * (-dndv * dot(n, dpdu) - n * dot(dndv, dpdu));
+
+	du.s -= s * dot(du.s, s);
+	dv.s -= s * dot(dv.s, s);
+
+	du.t = cross(dndu, s) + cross(n, du.s);
+	dv.t = cross(dndv, s) + cross(n, dv.s);
+
+	du.n = dndu;
+	dv.n = dndv;
 }
 
 Point2 toSphericalCoordinates(const Vector &v) {

src/librender/bsdf.cpp

@@ -65,31 +65,14 @@ Float BSDF::getEta() const {
 
 Frame BSDF::getFrame(const Intersection &its) const {
 	Frame result;
-
+	computeShadingFrame(its.shFrame.n, its.dpdu, result);
-	result.n = its.shFrame.n;
-	result.s = normalize(its.dpdu - result.n
-		* dot(result.n, its.dpdu));
-	result.t = cross(result.n, result.s);
-
 	return result;
 }
 
 void BSDF::getFrameDerivative(const Intersection &its, Frame &du, Frame &dv) const {
-	(its.instance ? its.instance : its.shape)->getNormalDerivative(its, du.n, dv.n, true);
+	Vector dndu, dndv;
-
+	(its.instance ? its.instance : its.shape)->getNormalDerivative(its, dndu, dndv, true);
-	Vector n = its.shFrame.n;
+	computeShadingFrameDerivative(its.shFrame.n, its.dpdu, dndu, dndv, du, dv);
-	Vector s = its.dpdu - n * dot(n, its.dpdu);
-	Float invLen_s = 1.0f / s.length();
-	s *= invLen_s;
-
-	du.s = invLen_s * (-du.n * dot(n, its.dpdu) - n * dot(du.n, its.dpdu));
-	dv.s = invLen_s * (-dv.n * dot(n, its.dpdu) - n * dot(dv.n, its.dpdu));
-
-	du.s -= s * dot(du.s, s);
-	dv.s -= s * dot(dv.s, s);
-
-	du.t = cross(du.n, s) + cross(n, du.s);
-	dv.t = cross(dv.n, s) + cross(n, dv.s);
 }
 
 Float BSDF::getRoughness(const Intersection &its, int component) const {

src/shapes/cylinder.cpp

@@ -223,8 +223,7 @@ public:
 
 		if (m_flipNormals)
 			its.geoFrame.n *= -1;
-		its.shFrame = its.geoFrame;
+		its.shFrame.n = its.geoFrame.n;
-		its.wi = its.toLocal(-ray.d);
 		its.hasUVPartials = false;
 		its.instance = NULL;
 		its.time = ray.time;

src/shapes/deformable.cpp

@@ -477,23 +477,13 @@ public:
 
 			its.shFrame.n = normalize(n0 * b.x + n1 * b.y + n2 * b.z);
 
-			if (EXPECT_TAKEN(!vertexTangents0)) {
-				coordinateSystem(its.shFrame.n, its.shFrame.s, its.shFrame.t);
-			} else {
-				/* Align shFrame.s with dpdu, use Gram-Schmidt to orthogonalize */
-				its.shFrame.s = normalize(its.dpdu - its.shFrame.n
-					* dot(its.shFrame.n, its.dpdu));
-				its.shFrame.t = cross(its.shFrame.n, its.shFrame.s);
-			}
-
 			/* Ensure that the geometric & shading normals face the same direction */
 			if (dot(faceNormal, its.shFrame.n) < 0)
 				faceNormal = -faceNormal;
-
-			its.geoFrame = Frame(faceNormal);
 		} else {
-			its.shFrame = its.geoFrame = Frame(faceNormal);
+			its.shFrame.n = faceNormal;
 		}
+		its.geoFrame = Frame(faceNormal);
 
 		if (EXPECT_TAKEN(vertexTexcoords0)) {
 			Point2
@@ -515,7 +505,6 @@ public:
 				result[2], Spectrum::EReflectance);
 		}
 
-		its.wi = its.toLocal(-ray.d);
 		its.shape = m_kdtree->getMesh(0, cache->shapeIndex);
 		its.hasUVPartials = false;
 		its.primIndex = cache->primIndex;

src/shapes/disk.cpp

@@ -189,12 +189,9 @@ public:
 			its.dpdv = trafo(Vector(0, 1, 0));
 		}
 
-		its.shFrame = its.geoFrame = Frame(
+		its.shFrame.n = normalize(trafo(Normal(0, 0, 1)));
-			normalize(its.dpdu), normalize(its.dpdv),
-			normalize(trafo(Normal(0, 0, 1))));
 		its.uv = Point2(r, phi * INV_TWOPI);
 		its.p = ray(its.t);
-		its.wi = its.toLocal(-ray.d);
  		its.hasUVPartials = false;
 		its.instance = NULL;
 		its.time = ray.time;

src/shapes/hair.cpp

@@ -845,8 +845,7 @@ void HairShape::fillIntersectionRecord(const Ray &ray,
 	const Vector local = its.geoFrame.toLocal(relHitPoint);
 	its.p += its.geoFrame.n * (m_kdtree->getRadius() - std::sqrt(local.y*local.y+local.z*local.z));
 
-	its.shFrame = its.geoFrame;
+	its.shFrame.n = its.geoFrame.n;
-	its.wi = its.toLocal(-ray.d);
 	its.hasUVPartials = false;
 	its.instance = this;
 	its.time = ray.time;

src/shapes/heightfield.cpp

@@ -402,11 +402,8 @@ public:
 			its.shFrame.n = normalize(m_objectToWorld(Normal(
 				(1 - temp.p.x) * ((1-temp.p.y) * n00 + temp.p.y * n01)
 				   + temp.p.x  * ((1-temp.p.y) * n10 + temp.p.y * n11))));
-
-			its.shFrame.s = normalize(its.geoFrame.s - dot(its.geoFrame.s, its.shFrame.n) * its.shFrame.n);
-			its.shFrame.t = cross(its.shFrame.n, its.shFrame.s);
 		} else {
-			its.shFrame = its.geoFrame;
+			its.shFrame.n = its.geoFrame.n;
 		}
 
 		if (m_flipNormals) {
@@ -415,7 +412,6 @@ public:
 		}
 
 		its.shape = this;
- 		its.wi = its.toLocal(-ray.d);
  		its.hasUVPartials = false;
 		its.instance = NULL;
 		its.time = ray.time;

src/shapes/instance.cpp

@@ -133,15 +133,11 @@ void Instance::adjustTime(Intersection &its, Float time) const {
 	Transform trafo = m_transform->eval(its.time).inverse();
 	trafo = m_transform->eval(time) * trafo;
 
-	Vector s = trafo(its.shFrame.s);
-	its.shFrame.n = normalize(trafo(its.shFrame.n));
-	its.shFrame.s = normalize(s - its.shFrame.n
-		* dot(its.shFrame.n, s));
-	its.shFrame.t = cross(its.shFrame.n, its.shFrame.s);
-	its.geoFrame = Frame(normalize(trafo(its.geoFrame.n)));
 	its.dpdu = trafo(its.dpdu);
 	its.dpdv = trafo(its.dpdv);
+	its.geoFrame = Frame(normalize(trafo(its.geoFrame.n)));
 	its.p = trafo(its.p);
+	computeShadingFrame(normalize(trafo(its.shFrame.n)), its.dpdu, its.shFrame);
 	its.wi = normalize(trafo(its.wi));
 	its.instance = this;
 	its.time = time;
@@ -155,16 +151,11 @@ void Instance::fillIntersectionRecord(const Ray &_ray,
 	trafo.inverse()(_ray, ray);
 	kdtree->fillIntersectionRecord<false>(ray, temp, its);
 
-	Vector s = trafo(its.shFrame.s);
 	its.shFrame.n = normalize(trafo(its.shFrame.n));
-	its.shFrame.s = normalize(s - its.shFrame.n
-		* dot(its.shFrame.n, s));
-	its.shFrame.t = cross(its.shFrame.n, its.shFrame.s);
 	its.geoFrame = Frame(normalize(trafo(its.geoFrame.n)));
 	its.dpdu = trafo(its.dpdu);
 	its.dpdv = trafo(its.dpdv);
 	its.p = trafo(its.p);
-	its.wi = normalize(its.shFrame.toLocal(-_ray.d));
 	its.instance = this;
 }
 

src/shapes/rectangle.cpp

@@ -155,13 +155,13 @@ public:
 	void fillIntersectionRecord(const Ray &ray,
 			const void *temp, Intersection &its) const {
 		const Float *data = static_cast<const Float *>(temp);
-		its.shFrame = its.geoFrame = m_frame;
+		its.geoFrame = m_frame;
+		its.shFrame = its.geoFrame.n;
 		its.shape = this;
 		its.dpdu = m_dpdu;
 		its.dpdv = m_dpdv;
 		its.uv = Point2(0.5f * (data[0]+1), 0.5f * (data[1]+1));
 		its.p = ray(its.t);
-		its.wi = its.toLocal(-ray.d);
  		its.hasUVPartials = false;
 		its.instance = NULL;
 		its.time = ray.time;

src/shapes/sphere.cpp

@@ -248,8 +248,7 @@ public:
 		if (m_flipNormals)
 			its.geoFrame.n *= -1;
 
- 		its.shFrame = its.geoFrame;
+ 		its.shFrame.n = its.geoFrame.n;
- 		its.wi = its.toLocal(-ray.d);
  		its.hasUVPartials = false;
 		its.instance = NULL;
 		its.time = ray.time;