bdpt and mtl integrators render sphere incorrectly due to precision limitations

[BrianSwift]

# When 32-bit float pbrt renders this scene of a sphere positioned away from origin
# with bdpt or mtl integrators a dark ring is produced near the limb.
# Using pbrt built with -DPBRT_FLOAT_AS_DOUBLE=1 eliminates the issue.
#
# Other integrators do not exhibit the issue.
#
# If position and size are decreased by facor of 10, issues does not appear.
# If position and size are increased by facor of 10, pbrt aborts with traceback from
#   transform.h:228] Check failed: wp != 0 (0 vs. 0)

# Environment: MacOS 10.14.6, Xcode 11.2.1
#   pbrt-v3 checkout on 12/5/2019
#   commit 416a080082a27e1f204e31ffe58e5e5cab2d179c (HEAD -> master, origin/master, origin/HEAD)

Film "image" "string filename" "s6.png"
     "integer xresolution" [400] "integer yresolution" [400]

Sampler "halton" "integer pixelsamples" 128

# These integrators produce expected results
#Integrator "directlighting"
#Integrator "path"
#Integrator "sppm"
#Integrator "whitted"
#Integrator "volpath"

# These integrators show problem
Integrator "bdpt"
#Integrator "mlt"

# 800000 scale FAILS to render correctly. Has dark ring near limb.
LookAt 796800. 0 0   # eye
       800000. 0 0  # look at point
       0 0 1    # up vector

Camera "perspective" "float fov" 4

WorldBegin

LightSource "distant"  "point to" [ 1 0  0 ]
   "blackbody L" [6500 1.5]

AttributeBegin
Material "matte" "rgb Kd" [ .9 .9 .9 ]
  Translate 800000. 0 0
  Shape "sphere" "float radius" 80.
AttributeEnd

WorldEnd

[mmp]

This is a tough case numerically--quite far away spheres and a narrow fov. As it turns out, the underlying precision issue is nicely described in this Ray Tracing Gems chapter.

With the fix described there (diff of a hacked in version of it below), results are much better (but not perfect, on my system at least.) However, my inclination is to leave the github pbrt-v3 code as is, in the interests of minimizing divergence from the book text. This is already fixed in the eventually forthcoming pbrt-v4.

Thanks!

diff --git a/src/shapes/sphere.cpp b/src/shapes/sphere.cpp
index 40b819a..51d5ed9 100755
--- a/src/shapes/sphere.cpp
+++ b/src/shapes/sphere.cpp
@@ -66,7 +66,32 @@ bool Sphere::Intersect(const Ray &r, Float *tHit, SurfaceInteraction *isect,

     // Solve quadratic equation for _t_ values
     EFloat t0, t1;
+#if 0
     if (!Quadratic(a, b, c, &t0, &t1)) return false;
+#else
+    // RT Gems
+    EFloat rootDiscrim;
+    {
+    EFloat f = b / (2 * a);  // (o . d) / LengthSquared(d)
+    EFloat fp[3] = { ox - f * dx, oy - f * dy, oz - f * dz };
+    EFloat sqrtf = sqrt(fp[0] * fp[0] + fp[1] * fp[1] + fp[2] * fp[2]);
+    EFloat discrim = 4 * a * ((EFloat(radius)) - sqrtf) *
+        ((EFloat(radius)) + sqrtf);
+
+    if (discrim.LowerBound() < 0) return false;
+    rootDiscrim = sqrt(discrim);
+    }
+
+    // Compute quadratic _t_ values
+    EFloat q;
+    if ((Float)b < 0)
+        q = -.5 * (b - rootDiscrim);
+    else
+        q = -.5 * (b + rootDiscrim);
+    t0 = q / a;
+    t1 = c / q;
+    if (t0.LowerBound() > t1.LowerBound()) std::swap(t0, t1);
+#endif

     // Check quadric shape _t0_ and _t1_ for nearest intersection
     if (t0.UpperBound() > ray.tMax || t1.LowerBound() <= 0) return false;
@@ -174,7 +199,32 @@ bool Sphere::IntersectP(const Ray &r, bool testAlphaTexture) const {

     // Solve quadratic equation for _t_ values
     EFloat t0, t1;
+#if 0
     if (!Quadratic(a, b, c, &t0, &t1)) return false;
+#else
+    // RT Gems
+    EFloat rootDiscrim;
+    {
+    EFloat f = b / (2 * a);  // (o . d) / LengthSquared(d)
+    EFloat fp[3] = { ox - f * dx, oy - f * dy, oz - f * dz };
+    EFloat sqrtf = sqrt(fp[0] * fp[0] + fp[1] * fp[1] + fp[2] * fp[2]);
+    EFloat discrim = 4 * a * ((EFloat(radius)) - sqrtf) *
+        ((EFloat(radius)) + sqrtf);
+
+    if (discrim.LowerBound() < 0) return false;
+    rootDiscrim = sqrt(discrim);
+    }
+
+    // Compute quadratic _t_ values
+    EFloat q;
+    if ((Float)b < 0)
+        q = -.5 * (b - rootDiscrim);
+    else
+        q = -.5 * (b + rootDiscrim);
+    t0 = q / a;
+    t1 = c / q;
+    if (t0.LowerBound() > t1.LowerBound()) std::swap(t0, t1);
+#endif

     // Check quadric shape _t0_ and _t1_ for nearest intersection
     if (t0.UpperBound() > ray.tMax || t1.LowerBound() <= 0) return false;

[BrianSwift]

Tried the above patch with a 32-bit float pbrt, and only difference appears to be a black pixel in the middle of image.