do sphere & cylinder intersections in double precision

include/mitsuba/core/util.h

@@ -89,6 +89,12 @@ extern MTS_EXPORT_CORE Float lanczosSinc(Float t, Float tau = 2);
  */
 extern MTS_EXPORT_CORE bool solveQuadratic(Float a, Float b, Float c, Float &x0, Float &x1);
 
+/**
+ * \brief Similar to solveQuadratic(), but always uses double precision independent
+ * of the chosen compile-time precision.
+ */
+extern MTS_EXPORT_CORE bool solveQuadraticDouble(double a, double b, double c, double &x0, double &x1);
+
 /**
  * \brief Solve a 2x2 linear equation system using basic linear algebra
  */

src/libcore/util.cpp

@@ -374,6 +374,46 @@ bool solveQuadratic(Float a, Float b, Float c, Float &x0, Float &x1) {
 	return true;
 }
 
+bool solveQuadraticDouble(double a, double b, double c, double &x0, double &x1) {
+	/* Linear case */
+	if (a == 0) {
+		if (b != 0) {
+			x0 = x1 = -c / b;
+			return true;
+		}
+		return false;
+	}
+
+	double discrim = b*b - 4.0f*a*c;
+
+	/* Leave if there is no solution */
+	if (discrim < 0)
+		return false;
+
+	double temp, sqrtDiscrim = std::sqrt(discrim);
+
+	/* Numerically stable version of (-b (+/-) sqrtDiscrim) / (2 * a)
+	 *
+	 * Based on the observation that one solution is always
+	 * accurate while the other is not. Finds the solution of
+	 * greater magnitude which does not suffer from loss of
+	 * precision and then uses the identity x1 * x2 = c / a
+	 */
+	if (b < 0)
+		temp = -0.5f * (b - sqrtDiscrim);
+	else
+		temp = -0.5f * (b + sqrtDiscrim);
+
+	x0 = temp / a;
+	x1 = c / temp;
+
+	/* Return the results so that x0 < x1 */
+	if (x0 > x1)
+		std::swap(x0, x1);
+
+	return true;
+}
+
 bool solveLinearSystem2x2(const Float a[2][2], const Float b[2], Float x[2]) {
 	Float det = a[0][0] * a[1][1] - a[0][1] * a[1][0];
 

src/shapes/cylinder.cpp

@@ -69,29 +69,35 @@ public:
 	}
 
 	bool rayIntersect(const Ray &_ray, Float start, Float end, Float &t) const {
-		Float nearT, farT; Ray ray;
+		double nearT, farT; Ray ray;
 
 		/* Transform into the local coordinate system and normalize */
 		m_worldToObject(_ray, ray);
+		
+		const double
+			ox = ray.o.x,
+			oy = ray.o.y,
+			dx = ray.d.x, 
+			dy = ray.d.y;
 
-		const Float A = ray.d.x*ray.d.x + ray.d.y*ray.d.y;
-		const Float B = 2 * (ray.d.x*ray.o.x + ray.d.y*ray.o.y);
-		const Float C = ray.o.x*ray.o.x + ray.o.y*ray.o.y - m_radius*m_radius;
+		const double A = dx*dx + dy*dy;
+		const double B = 2 * (dx*ox + dy*oy);
+		const double C = ox*ox + oy*oy - m_radius*m_radius;
 
-		if (!solveQuadratic(A, B, C, nearT, farT))
+		if (!solveQuadraticDouble(A, B, C, nearT, farT))
 			return false;
 
 		if (nearT > end || farT < start)
 			return false;
 
-		const Float zPosNear = ray.o.z + ray.d.z * nearT;
-		const Float zPosFar = ray.o.z + ray.d.z * farT;
+		const double zPosNear = ray.o.z + ray.d.z * nearT;
+		const double zPosFar = ray.o.z + ray.d.z * farT;
 		if (zPosNear >= 0 && zPosNear <= m_length && nearT >= start) {
-			t = nearT;
+			t = (Float) nearT;
 		} else if (zPosFar >= 0 && zPosFar <= m_length) {
 			if (farT > end)
 				return false;
-			t = farT;
+			t = (Float) farT;
 		} else {
 			return false;
 		}
@@ -124,20 +130,6 @@ public:
 		its.hasUVPartials = false;
 		its.shape = this;
 
-		/* Intersection refinement step */
-		Vector2 localDir(normalize(Vector2(local.x, local.y)));
-		Vector rel = its.p - m_objectToWorld(Point(m_radius  * localDir.x, 
-			m_radius * localDir.y, local.z));
-		Float correction = -dot(rel, its.geoFrame.n)/dot(ray.d, its.geoFrame.n);
-
-		its.t += correction;
-		if (its.t < ray.mint || its.t > ray.maxt) {
-			its.t = std::numeric_limits<Float>::infinity(); 
-			return false;
-		}
-
-		its.p += ray.d * correction;
-
 		return true;
 	}
 

src/shapes/sphere.cpp

@@ -38,14 +38,16 @@ public:
 
 	bool rayIntersect(const Ray &ray, Float start, Float end, Float &t) const {
 		/* Transform into the local coordinate system and normalize */
-		Float nearT, farT;
-		const Float ox = ray.o.x - m_center.x, oy = ray.o.y - m_center.y, 
-			oz = ray.o.z - m_center.z;
-		const Float A = ray.d.x*ray.d.x + ray.d.y*ray.d.y + ray.d.z*ray.d.z;
-		const Float B = 2 * (ray.d.x*ox + ray.d.y*oy + ray.d.z*oz);
-		const Float C = ox*ox + oy*oy + oz*oz - m_radius * m_radius;
+		double nearT, farT;
+		const double ox = (double) ray.o.x - (double) m_center.x, 
+			oy = (double) ray.o.y - (double) m_center.y, 
+			oz = (double) ray.o.z - (double) m_center.z;
+		const double dx = ray.d.x, dy = ray.d.y, dz = ray.d.z;
+		const double A = dx*dx + dy*dy + dz*dz;
+		const double B = 2 * (dx*ox + dy*oy + dz*oz);
+		const double C = ox*ox + oy*oy + oz*oz - m_radius * m_radius;
 
-		if (!solveQuadratic(A, B, C, nearT, farT))
+		if (!solveQuadraticDouble(A, B, C, nearT, farT))
 			return false;
 
 		if (nearT > end || farT < start)
@@ -53,9 +55,9 @@ public:
 		if (nearT < start) {
 			if (farT > end)
 				return false;
-			t = farT;		
+			t = (Float) farT;		
 		} else {
-			t = nearT;
+			t = (Float) nearT;
 		}
 		return true;
 	}
@@ -99,19 +101,6 @@ public:
 			coordinateSystem(its.geoFrame.n, its.geoFrame.s, its.geoFrame.t);
 		}
 
-		/* Intersection refinement step */
-		Vector rel = its.p - (m_center + its.geoFrame.n * m_radius);
-		Float correction = -dot(rel, its.geoFrame.n)/dot(ray.d, its.geoFrame.n);
-
-		/* Move outside of the object */
-		its.t += correction;
-
-		if (its.t < ray.mint || its.t > ray.maxt) {
-			its.t = std::numeric_limits<Float>::infinity(); 
-			return false;
-		}
-
-		its.p += ray.d * correction;
  		its.shFrame = its.geoFrame;
  		its.wi = its.toLocal(-ray.d);
  		its.hasUVPartials = false;