From 995810faaf6e2e2def1910213edcc0b4130ab38f Mon Sep 17 00:00:00 2001 From: Wenzel Jakob Date: Thu, 2 Sep 2010 22:32:33 +0200 Subject: [PATCH] do sphere & cylinder intersections in double precision --- include/mitsuba/core/util.h | 6 ++++++ src/libcore/util.cpp | 40 +++++++++++++++++++++++++++++++++++++ src/shapes/cylinder.cpp | 38 ++++++++++++++--------------------- src/shapes/sphere.cpp | 33 ++++++++++-------------------- 4 files changed, 72 insertions(+), 45 deletions(-) diff --git a/include/mitsuba/core/util.h b/include/mitsuba/core/util.h index 1fdbbc5c..19d8b91a 100644 --- a/include/mitsuba/core/util.h +++ b/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 */ diff --git a/src/libcore/util.cpp b/src/libcore/util.cpp index bc092816..76630128 100644 --- a/src/libcore/util.cpp +++ b/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]; diff --git a/src/shapes/cylinder.cpp b/src/shapes/cylinder.cpp index 91d19846..1072e8b2 100644 --- a/src/shapes/cylinder.cpp +++ b/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::infinity(); - return false; - } - - its.p += ray.d * correction; - return true; } diff --git a/src/shapes/sphere.cpp b/src/shapes/sphere.cpp index a06c0ff5..636abad9 100644 --- a/src/shapes/sphere.cpp +++ b/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::infinity(); - return false; - } - - its.p += ray.d * correction; its.shFrame = its.geoFrame; its.wi = its.toLocal(-ray.d); its.hasUVPartials = false;