POV Torus bug fix

The Persistence of Vision Raytracer (3.5) has a problem with the torus when the ray is cast from a long distance. The intersection with the torus is calculated by solving roots of a fourth order equation. When the origin of the ray is far away, the standard double-precision floating-point cannot represent the accuracy needed.

The result is that the torus can have black speckles or simply disappear!

The problem has been reported in POV newsgroups:
http://news.povray.org/povray.general/3330/
http://news.povray.org/povray.general/2610/
http://news.povray.org/povray.binaries.images/12239/

The problem also shows when rendering LDraw models, typically when using the L3P Camera Angle option -ca1 to have a rendering with almost orthogonal camera. This is achieved by narrowing the camera angle to 1 degree and moving the camera far away to maintain the model in the view.

Here is one of the models in the 1995 set 4555 Freight Loading Station, a racing car (l3p -b -f -ca1 -q3 m4555d.dat):
 

Original rendering

Rendering with bug fix applied

The stud logo is too dark, black speckles are showing.

Here is a more obvious example, the Plate 1 x 1 Round" (l3p -b15 -c2 -ca1 -q3 -cg45,45 4073):
 

Original rendering

Rendering with bug fix applied

Here is a close up of some of the tori:
 

Original rendering

Rendering with bug fix applied

Actually the problem was brought to my attention by Jeroen de Haan while he was beta testing the coming L3P. It features primitive substitution for the 1-4ccyli.dat "Cylinder Tube 0.25". The quarter torus was missing in some of his renderings. We tried an "isosurface { function { f_torus(...) ... }" and it did work on his specific rendering. However, while testing the isosurface further I saw it had the same problem. So I decided to strike at the root of the problem. I downloaded the POV source code and began debugging hinted by some of the replies in the postings mentioned above.

The bug fix

A sphere enclosing the torus has radius R+r. However, to avoid any problems I use radius R+r+r which safely encloses the torus. Then, if the origin of the ray P is outside that sphere, I move the origin close to the sphere. There's no need to calculate the exact intersection of the ray and the enclosing sphere. Close to the sphere is fine. This is for speed considerations.

As you can see in the images above the bug fix solves the problem. I have also tested the bug fix with positive results on the samples in the postings mentioned above.

Function intersect_torus in file torus.cpp (the changes are marked in red)

static int intersect_torus(RAY *Ray, TORUS *Torus, DBL *Depth)
{
  int i, n;
  DBL len, R2, Py2, Dy2, PDy2, k1, k2;
  DBL y1, y2, r1, r2;
  DBL c[5];
  DBL r[4];
  VECTOR P, D;
  DBL DistanceP; // Distance from P to torus center (origo).
  DBL BoundingSphereRadius; // Sphere fully (amply) enclosing torus.
  DBL Closer; // P is moved Closer*D closer to torus.

  Increase_Counter(stats[Ray_Torus_Tests]);

  /* Transform the ray into the torus space. */

  MInvTransPoint(P, Ray->Initial, Torus->Trans);

  MInvTransDirection(D, Ray->Direction, Torus->Trans);

  VLength(len, D);

  VInverseScaleEq(D, len);

  i = 0;

  y1 = -Torus->r;
  y2 =  Torus->r;
  r1 = Sqr(Torus->R - Torus->r);

  if ( Torus->R < Torus->r ) {
    r1 = 0;
  }
  
  r2 = Sqr(Torus->R + Torus->r);

#ifdef TORUS_EXTRA_STATS
  Increase_Counter(stats[Torus_Bound_Tests]);
#endif

  if (Test_Thick_Cylinder(P, D, y1, y2, r1, r2))
  {
#ifdef TORUS_EXTRA_STATS
    Increase_Counter(stats[Torus_Bound_Tests_Succeeded]);
#endif

    // Move P close to bounding sphere to have more precise root calculation.
    // Bounding sphere radius is R + r, we add r once more to ensure
    // that P is safely outside sphere.
    BoundingSphereRadius = Torus->R + Torus->r + Torus->r;
    DistanceP = VSumSqr(P); // Distance is currently squared.
    Closer = 0.0;
    if (DistanceP > Sqr(BoundingSphereRadius))
    {
      DistanceP = sqrt(DistanceP); // Now real distance.
      Closer = DistanceP - BoundingSphereRadius;
      VAddScaledEq(P, Closer, D);
    }

    R2   = Sqr(Torus->R);
    r2   = Sqr(Torus->r);

    Py2  = P[Y] * P[Y];
    Dy2  = D[Y] * D[Y];
    PDy2 = P[Y] * D[Y];

    k1   = P[X] * P[X] + P[Z] * P[Z] + Py2 - R2 - r2;
    k2   = P[X] * D[X] + P[Z] * D[Z] + PDy2;

    c[0] = 1.0;

    c[1] = 4.0 * k2;

    c[2] = 2.0 * (k1 + 2.0 * (k2 * k2 + R2 * Dy2));

    c[3] = 4.0 * (k2 * k1 + 2.0 * R2 * PDy2);

    c[4] = k1 * k1 + 4.0 * R2 * (Py2 - r2);

    n = Solve_Polynomial(4, c, r, Test_Flag(Torus, STURM_FLAG), ROOT_TOLERANCE);

    while(n--)
    {
      Depth[i++] = r[n] / len;
      Depth[i++] = (r[n] + Closer) / len;
    }
  }

  if (i)
  {
    Increase_Counter(stats[Ray_Torus_Tests_Succeeded]);
  }

  return(i);
}

Download

Lars C. Hassing's Homepage

Last updated September 25, 2004