Мне нужно повысить точность для функции в одном из моихАтмосферное рассеяние фрагмента GLSL который вычисляет пересечение между одиночным лучом и выровненным по оси эллипсоидом.

Это основная функция моего шейдера атмосферного рассеяния. Старый оригинальный шейдер был включенfloats и для нормального рендеринга было нормально, но после добавления зума я обнаружил, что при относительно небольших расстояниях точность теряется. На поплавках полезные расстояния для Земли составляли только до 0,005 а.е. (астрономическая единица). Поэтому я попытался перенести критическую функцию наdouble и это помогает, так что теперь полезное расстояние составляет около 1,0 AU (с небольшими артефактами)

Этоdouble версия функции внутри Fragment Shader (в старом стиле используются устаревшие вещи !!!)

#extension GL_ARB_gpu_shader_fp64 : enable
double abs(double x) { if (x<0.0) x=-x; return x; }
// compute length of ray(p0,dp) to intersection with ellipsoid((0,0,0),r) -> view_depth_l0,1
// where r.x is elipsoid rx^-2, r.y = ry^-2 and r.z=rz^-2
float view_depth_l0=-1.0,view_depth_l1=-1.0;
bool _view_depth(vec3 _p0,vec3 _dp,vec3 _r)
    {
    double a,b,c,d,l0,l1;
    dvec3 p0,dp,r;
    p0=dvec3(_p0);
    dp=dvec3(_dp);
    r =dvec3(_r );
    view_depth_l0=-1.0;
    view_depth_l1=-1.0;
    a=(dp.x*dp.x*r.x)
     +(dp.y*dp.y*r.y)
     +(dp.z*dp.z*r.z); a*=2.0;
    b=(p0.x*dp.x*r.x)
     +(p0.y*dp.y*r.y)
     +(p0.z*dp.z*r.z); b*=2.0;
    c=(p0.x*p0.x*r.x)
     +(p0.y*p0.y*r.y)
     +(p0.z*p0.z*r.z)-1.0;
    d=((b*b)-(2.0*a*c));
    if (d<0.0) return false;
    d=sqrt(d);
    l0=(-b+d)/a;
    l1=(-b-d)/a;
    if (abs(l0)>abs(l1)) { a=l0; l0=l1; l1=a; }
    if (l0<0.0)          { a=l0; l0=l1; l1=a; }
    if (l0<0.0) return false;
    view_depth_l0=float(l0);
    view_depth_l1=float(l1);
    return true;
    }

вход луча и радиусов ^ -2 эллипсоида

выходное расстояние от р0 до пересечений

Точность входных и выходных переменных float (достаточно)

Вот так это выглядит после портирования на Double

Итак, вопрос: Q1. Как повысить точность этой функции?

целевая точность дляview_depth_l0,view_depth_l1 является+/- 20 m для расстояния|p0|=100 AU

Это было бы идеально, сейчас кажется, что +/- 5 км на расстоянии 10 а.е., что плохо. Даже 10-кратное точное вычисление будет огромным шагом вперед, есть идеи?

[edit1] l0, l1 range

Я ошибся с плавающей точкой преобразованияview_depth_l0,view_depth_l1 является причиной артефактов. После смещения на относительное расстояние точность значительно улучшилась. Я только добавил это:

    // relative shift to preserve accuracy
    const double m0=1000000000.0; // >= max view depth !!!
    if (l0>m0){ a=floor(l0/m0)*m0; a-=m0; if (l1>l0) l1-=a; l0-=a; }

до этого:

    view_depth_l0=float(l0);
    view_depth_l1=float(l1);
    return true;
    }

Остальная часть шейдерной ручкиl0,l1 В качестве относительных значений в любом случае получим следующее:

для расстояний до 10,0 а.е. это нормально (артефакты заметны только в очень больших масштабах), новые артефакты возникают скорее всего в других местах, поэтому мне придется продолжить исследования, когда у меня будет время и желание.

[edit2] смещение p0 вдоль dp ближе к (0,0,0)

Реализация требует относительно дорогих функций нормализации и длины, и результат без смещения диапазона (edit1) был немного лучше, чем необработанная функция, но улучшение не слишком большое. С изменением диапазона (edit1) результат тот же, что и раньше, поэтому это не так. Мой вывод заключается в том, что все оставшиеся артефакты не вызваны самой функцией просмотра.

Я попробую портировать шейдер на#version 400 + fp64 в целом, чтобы проверить, не слишком ли округлены входные данные с плавающей точкой

[edit3] фактический исходный код

#extension GL_ARB_gpu_shader_fp64 : enable
double abs(double x) { if (x<0.0) x=-x; return x; }
// compute length of ray(p0,dp) to intersection with ellipsoid((0,0,0),r) -> view_depth_l0,1
// where r.x is elipsoid rx^-2, r.y = ry^-2 and r.z=rz^-2
float view_depth_ll= 0.0, // shift to boost accuracy
      view_depth_l0=-1.0, // view_depth_ll+view_depth_l0 first hit
      view_depth_l1=-1.0; // view_depth_ll+view_depth_l1 second hit
const double view_depth_max=100000000.0; // > max view depth
bool _view_depth(vec3 _p0,vec3 _dp,vec3 _r)
    {
    dvec3 p0,dp,r;
    double a,b,c,d,l0,l1;
    view_depth_ll= 0.0;
    view_depth_l0=-1.0;
    view_depth_l1=-1.0;
    // conversion to double
    p0=dvec3(_p0);
    dp=dvec3(_dp);
    r =dvec3(_r );
    // quadratic equation a.l.l+b.l+c=0; l0,l1=?;
    a=(dp.x*dp.x*r.x)
     +(dp.y*dp.y*r.y)
     +(dp.z*dp.z*r.z);
    b=(p0.x*dp.x*r.x)
     +(p0.y*dp.y*r.y)
     +(p0.z*dp.z*r.z); b*=2.0;
    c=(p0.x*p0.x*r.x)
     +(p0.y*p0.y*r.y)
     +(p0.z*p0.z*r.z)-1.0;
    // discriminant d=sqrt(b.b-4.a.c)
    d=((b*b)-(4.0*a*c));
    if (d<0.0) return false;
    d=sqrt(d);
    // standard solution l0,l1=(-b +/- d)/2.a
    a*=2.0;
    l0=(-b+d)/a;
    l1=(-b-d)/a;
    // alternative solution q=-0.5*(b+sign(b).d) l0=q/a; l1=c/q; (should be more accurate sometimes)
//  if (b<0.0) d=-d; d=-0.5*(b+d);
//  l0=d/a;
//  l1=c/d;
    // sort l0,l1 asc
    if (abs(l0)>abs(l1)) { a=l0; l0=l1; l1=a; }
    if (l0<0.0)          { a=l0; l0=l1; l1=a; }
    if (l0<0.0) return false;
    // relative shift to preserve accuracy after conversion back float
    if (l0>view_depth_max){ a=floor(l0/view_depth_max)*view_depth_max; a-=view_depth_max; view_depth_ll=float(a); if (l1>l0) l1-=a; l0-=a; }
    // conversion back float
    view_depth_l0=float(l0);
    view_depth_l1=float(l1);
    return true;
    }

Перенос остальных шейдеров на удвоение не имеет никакого эффекта. Единственное, что может улучшить этоdouble входные данные (входdouble но GL преобразовать его вfloat), но мой нынешний GLSL HW не позволяет64 bit интерполяторы

Q2. Есть ли способ пройтиdouble интерполяторы из вершины в фрагментный шейдер?

Что-то вродеvarying dvec4 pixel_pos; в старом стиле GLSL илиout smooth dvec4 pixel_pos; в профиле ядра