Численная ошибка в пересечении уравнения луча и тора, когда камера находится далеко от тора

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 функция берется из решения кубической и катрической функций.

Проблема в том, что когда мы внимательно смотрим на тор, он работает довольно хорошо:

Но когда я уменьшаю камеру, чтобы камера смотрела на тору вдали от нее, она внезапно становится такой шумной, и ее форма плохо определена. Я пытался использовать более 1 выборки на пиксели, но все же у меня та же проблема. Это выглядит следующим образом:

Кажется, я столкнулся с численной проблемой, но я не знаю, как ее решить. Кто-нибудь может мне помочь с этим?

Также стоит упомянуть, что я отслеживаю тору с помощью Intel 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

PovRay дает для этого интересное и эффективное решение. Просто переместите начало луча очень близко к тору, и коэффициент будет иметь хорошие значения для полиномиального решателя. Насколько близко: начало координат должно быть на сфере с радиусом mayor+2*minor. ... и оставьте радиус мэра равным одному, как это предлагает @csharpfolk

Answer 3

Когда я писал свой raytracer (кстати, я использовал замечательную книгу под названием "Трассировка лучей с нуля"), у меня тоже были некоторые проблемы с Торусом. Тогда я использовал алгоритмы github repo от Graphics Gems для вычисления точек пересечения луча и тора. Решение состояло в том, чтобы просто использовать меньшие торы, например, когда мой тор имел внешний радиус, превышающий 100.0 и луч начался в (0,0,0) мой raytracer столкнулся с множеством ошибок. Используя меньшие радиусы тора, как 1.0 решил мои проблемы.

Источник этих числовых ошибок лежит в построении коэффициента для многочлена тора с тором размера 100.0 некоторый коэффициент, который генерируется во время этого вычисления, может превышать 1e20, С double точность, которая гарантирует около 15 значащих цифр, что привело к значительной потере точности.