当相机远离环面时，射线和环面方程相交的数值错误

Question

我试图在不对环面进行三角剖分的情况下对环面进行光线追踪，而只是通过将光线和环面解析方程相交。我用以下代码做到了这一点：

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));
    }
  }
}

求解SolveP4函数方程的代码取自Solution of cubic and quatric functions。

问题是当我们仔细观察环面时，它工作得很好，如下所示：

但是当我缩小相机时，相机正在看着远离它的环面，它突然变得如此嘈杂并且形状无法很好地识别。我尝试每个像素使用超过 1 个样本，但我仍然遇到同样的问题。具体如下：

看来我遇到了一个数值问题，但我不知道如何解决。任何人都可以帮助我吗？

此外，值得一提的是，我正在使用英特尔的 Embree Lib 对环面进行光线追踪。

更新（单色）：

Answer 1

我认为很多问题是使用单精度浮点数而不是双精度。

定义两个函数

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;
}

求平方和点积，但使用双精度。更改 r2 R2 a4 a3 a2 a1 和 a0 的计算以使用这些

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;

其余代码全部相同。在我的测试中，这使得看起来模糊的图像看起来非常清晰。

Answer 2

当我写我的 raytracer 时（顺便说一句，我使用了一本名为 "Ray Tracing from the Ground Up" 的好书）我也遇到了 Toruses 的一些问题。那时我使用来自 Graphics Gems github repo 的算法来计算射线环相交点。解决方案是简单地使用较小的圆环，例如当我的环面的外半径大于 100.0 并且光线开始于 (0,0,0) 时，我的光线追踪器出现了大量的数值错误。 使用像 1.0 这样更小的环面半径解决了我的问题。

这些数值错误的根源在于环面多项式的构建系数，环面的大小为 100.0，在计算过程中生成的一些系数可能会超过 1e20。 double 保证大约 15 位有效数字的精度，造成了很大的精度损失。

Answer 3

PovRay 对此提供了有趣且有效的解决方案。只需将射线的原点移动到非常靠近环面的位置，多项式求解器的系数就会有很好的值。接近程度：原点应位于半径为 mayor+2*minor 的球体上。 ...并按照@csharpfolk

的建议将市长半径保持为 1

当相机远离环面时，射线和环面方程相交的数值错误

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

math

numerical

raytracing

numerical-methods

embree