Eikonal 3D implementation

float neighborMin(float3 uv, float3 delta)
{
    float a = unpack(sdfSample(uv + delta));
    float b = unpack(sdfSample(uv - delta));
    return min(abs(a), abs(b));
}

float eikonal1d(float h, float v, float g)
{
    return min(h, v) + g;
}

float eikonal2d(float h, float v, float g)
{
    if (abs(h - v) >= g)
    {
        return eikonal1d(h, v, g);
    }
    else
    {
        float hv = h + v;
        float d = hv * hv - 2.0 * (h * h + v * v - g * g);
        return 0.5 * (hv + sqrt(d));
    }
}

bool eikonal3d(float3 u, out float e, float g)
{
    float xyz = u.x + u.y + u.z;
    float n = 3.0;
    float d = (xyz * xyz) - n * ((u.x * u.x + u.y * u.y + u.z * u.z) - (g * g));
    if (d < 0.0)
    {
        return false;
    }
    else
    {
        e = (1.0 / n) * (xyz + sqrt(d));
        return true;
    }
}

float eikonal(float3 u, float g)
{
    // 3d
    float e = 0.0;
    if (eikonal3d(u, e, g))
    {
        return e;
    }
    else
    {
        // 2d
        float e1 = eikonal2d(u.x, u.y, g);
        float e2 = eikonal2d(u.x, u.z, g);
        float e3 = eikonal2d(u.y, u.z, g);
        return min(min(e1, e2), e3);
    }
}

[numthreads(8, 8, 8)]
void SDFFix(uint3 id : SV_DispatchThreadID)
{
    float3 kernel = float3(float(id.x), float(id.y), float(id.z));
    float3 sdfSize = float3(float(_sdfTexSizeX), float(_sdfTexSizeY), float(_sdfTexSizeZ));
    float3 uv = (kernel + float3(0.5, 0.5, 0.5)) / sdfSize;
    float3 worldPos = (uv - float3(0.5, 0.5, 0.5)) * _sdfRadius * 2.0;
    float current = unpack(sdfSample(uv));

    float3 u;
    u.x = neighborMin(uv, float3(1.0 / _sdfTexSizeX, 0.0, 0.0));
    u.y = neighborMin(uv, float3(0.0, 1.0 / _sdfTexSizeY, 0.0));
    u.z = neighborMin(uv, float3(0.0, 0.0, 1.0 / _sdfTexSizeZ));

    float g = (_sdfRadius * 2.0) / _sdfTexSizeX;
    float e = eikonal(u, g);

    if (current > g) current = e;

    _sdfTexOut[id.xyz] = pack(current);
}