Numerical bug in intersecting the equation of ray and torus when the camera is far from the torus

Question

I am trying to ray trace a torus without triangulating the torus and just by intersecting the ray and torus analytic equation. I did that with the following code:

void circularTorusIntersectFunc(const CircularTorus* circularToruses, RTCRay& ray, size_t item)
{
  const CircularTorus& torus = circularToruses[item];

  Vec3fa O = ray.org /*- sphere.p*/;
  Vec3fa Dir = ray.dir;
  O.w = 1.0f;
  Dir.w = 0.0f;
  O = torus.inv_transform.mult(O);
  Dir = torus.inv_transform.mult(Dir);

  // r1: cross section of torus
  // r2: the ring's radius
  //  _____                     ____
  // / r1  \------->r2<--------/    \
  // \_____/                   \____/

  float r2 = sqr(torus.r1);
  float R2 = sqr(torus.r2);

  double a4 = sqr(dot(Dir, Dir));
  double a3 = 4 * dot(Dir, Dir) * dot(O, Dir);
  double a2 = 4 * sqr(dot(O, Dir)) + 2 * dot(Dir, Dir) * (dot(O, O) - r2 - R2) + 4 * R2 * sqr(Dir.z);
  double a1 = 4 * dot(O, Dir) * (dot(O, O) - r2 - R2) + 8 * R2 * O.z * Dir.z;
  double a0 = sqr(dot(O, O) - r2 - R2) + 4 * R2 * sqr(O.z) - 4 * R2 * r2;

  a3 /= a4; a2 /= a4; a1 /= a4; a0 /= a4;

  double roots[4];
  int n_real_roots;
  n_real_roots = SolveP4(roots, a3, a2, a1, a0);

  if (n_real_roots == 0) return;

  Vec3fa intersect_point;
  for (int i = 0; i < n_real_roots; i++)
  {
    float root = static_cast<float>(roots[i]);
    intersect_point = root * Dir + O;

    if ((ray.tnear <= root) && (root <= ray.tfar)) {

      ray.u = 0.0f;
      ray.v = 0.0f;
      ray.tfar = root;
      ray.geomID = torus.geomID;
      ray.primID = item;
      Vec3fa normal(
        4.0 * intersect_point.x * (sqr(intersect_point.x) + sqr(intersect_point.y) + sqr(intersect_point.z) - r2 - R2),
        4.0 * intersect_point.y * (sqr(intersect_point.x) + sqr(intersect_point.y) + sqr(intersect_point.z) - r2 - R2),
        4.0 * intersect_point.z * (sqr(intersect_point.x) + sqr(intersect_point.y) + sqr(intersect_point.z) - r2 - R2) + 8 * R2*intersect_point.z,
        0.0f
        );

      ray.Ng = normalize(torus.transform.mult(normal));
    }
  }
}

The code to solve the equation for SolveP4 function is taken from Solution of cubic and quatric functions.

The problem is when we are looking at the torus closely, it works pretty nice as follows:

But when I zoom out the camera, so camera is looking at the torus far from it, it suddenly gets so noisy and it is shape is not well identified. I tried to use more than 1 samples per pixels but still I have the same problem. It is as follows:

It seems I am facing a numerical problem but I dont know how to solve it. Anyone can help me with that?

Also, it is good to mention that I am raytracing the torus with Intel's Embree Lib.

Update (Single Color):

Answer 1

I think a lot of the problem is using single precision float rather than double precision.

Define two functions

double dsqr(double x) { return x*x; }

double ddot(const Vec3fa &a,Vec3fa &b) {
  double x1 = a.x, y1 = a.y, z1 = a.z;
  double x2 = b.x, y2 = b.y, z2 = b.z;
  return x1*x2 + y1*y2 + z1*z2;
}

to find the square and the dot product but using double precision. Change the calculations of r2 R2 a4 a3 a2 a1 and a0 to use these

double r2 = dsqr(torus.r1);
double R2 = dsqr(torus.r2);

double a4 = dsqr(ddot(Dir, Dir));
double a3 = 4 * ddot(Dir, Dir) * ddot(O, Dir);
double a2 = 4 * dsqr(ddot(O, Dir)) + 2 * ddot(Dir, Dir) * (ddot(O, O) - r2 - R2)
    + 4 * R2 * dsqr(Dir.z);
double a1 = 4 * ddot(O, Dir) * (ddot(O, O) - r2 - R2) + 8 * R2 * O.z * Dir.z;
double a0 = dsqr(ddot(O, O) - r2 - R2) + 4 * R2 * dsqr(O.z) - 4 * R2 * r2;

all the remaining code is the same. In my test this made a fuzzy looking image look perfectly sharp.

Answer 2

When I was writing my raytracer (BTW I was using a great book called "Ray Tracing from the Ground Up") I had some problems with Toruses too. Back then I used algorithms from Graphics Gems github repo to compute ray torus intersection points. The solution was simply to use smaller toruses e.g. when my torus had outer radius greater than 100.0 and the ray started at (0,0,0) my raytracer experienced plenty of numerical errors. Using smaller torus radii like 1.0 solved my issues.

The source of these numerical errors lies in building coefficient for torus polynomial, with torus of size 100.0 some coefficient that are generated during that computation can exceed 1e20. With double precision that guarantees about 15 significant digits, that created substantial loss of precision.

Answer 3

PovRay give interesting and efficient solution for this. Just move origin of ray very close to torus and coefficient will have good values for polynomial solver. How close: origin should be on sphere with radius mayor+2*minor. ... and keep mayor radius to one as suggested by @csharpfolk