Обратный отбор проб с поверхности Ламберта

Question

В книге Физически обоснованный рендеринг , поверхность Ламберта отбирается следующим образом (см. http://www.pbr-book.org/3ed-2018/Light_Transport_I_Surface_Reflection/Sampling_Reflection_Functions.html#):

void Sample_f(Vector3f const& wo, Vector3f* wi, const Point2f& u)
{
    // Cosine-sample the hemisphere, flipping the direction if necessary
    *wi = CosineSampleHemisphere(u);
    if (wo.z < 0) wi->z *= -1;
}

inline Vector3f CosineSampleHemisphere(Point2f const& u)
{
    Point2f d = ConcentricSampleDisk(u);
    Float z = std::sqrt(std::max((Float)0, 1 - d.x * d.x - d.y * d.y));
    return Vector3f(d.x, d.y, z);
}

Point2f ConcentricSampleDisk(Point2f const& u)
{
    // Map uniform random numbers to $[-1,1]^2$
    Point2f uOffset = 2.f * u - Vector2f(1, 1);

    // Handle degeneracy at the origin
    if (uOffset.x == 0 && uOffset.y == 0) return Point2f(0, 0);

    // Apply concentric mapping to point
    Float theta, r;
    if (std::abs(uOffset.x) > std::abs(uOffset.y)) {
        r = uOffset.x;
        theta = PiOver4 * (uOffset.y / uOffset.x);
    } else {
        r = uOffset.y;
        theta = PiOver2 - PiOver4 * (uOffset.x / uOffset.y);
    }
    return r * Point2f(std::cos(theta), std::sin(theta));
}

Что я хочу сейчас сделать, так это, учитывая wo и wi, вычислить u так, чтобы вызов Sample_f(wo, &wi_other, u) дал wi_other == wi (по крайней мере приблизительно).

Пока это не трудно в принципе решить эту проблему, мое решение страдает от неточности с плавающей точкой. Если вы знакомы с трассировкой лучей: если луч, следующий за точно вычисленным направлением wi, достигает точки поверхности p, может оказаться, что приблизительно вычисленное направление wi_other близко не попадает на всю поверхность, на которой расположен p .

Это мое решение до сих пор:

Point2f invert_sample_f(pbrt::Vector3f wi, pbrt::Vector3f const& wo)
{
    if (wo.z < 0)
        wi.z *= -1;
    return cosine_sample_hemisphere_inverse(wi);
}

template<typename RealType = pbrt::Float>
pbrt::Point2<RealType> cosine_sample_hemisphere_inverse(pbrt::Vector3<RealType> const& w) {
    return concentric_map_inverse<RealType>({ w.x, w.y });
}

template<typename RealType = pbrt::Float>
pbrt::Point2<RealType> concentric_map_inverse(pbrt::Point2<RealType> u)
{
    u = cartesian_to_polar(u);
    auto const& r = u.x;
    auto& phi = u.y;

    if (r == 0)
        return { 0, 0 };

//  wrap ϕ -> [-π/4, 7π/4)
    if (phi >= 7 * pbrt::PiOver4)
        phi -= 2 * pbrt::Pi;

    if (-pbrt::PiOver4 < phi && phi < pbrt::PiOver4)
    {// sector 1
        u = { r, r * phi / pbrt::PiOver4 };
    }
    else if (pbrt::PiOver4 <= phi && phi <= 3 * pbrt::PiOver4)
    {// sector 2
        u = { r * (2 - phi / pbrt::PiOver4), r };
    }
    else if (3 * pbrt::PiOver4 < phi && phi < 5 * pbrt::PiOver4)
    {// sector 3
        u = { -r, r * (4 - phi / pbrt::PiOver4) };
    }
    else // 5 * pbrt::PiOver4 <= phi && phi <= -pbrt::PiOver4
    {// sector 4
        u = { r * (phi / pbrt::PiOver4 - 6), -r };
    }

    return (u + pbrt::Vector2<RealType>{ 1, 1 }) / 2;
}

template<typename RealType = pbrt::Float>
pbrt::Point2<RealType> cartesian_to_polar(pbrt::Point2<RealType> const& p)
{
    auto const &x = p.x,
        &y = p.y;

    RealType phi;
    if (x < 0)
        phi = pbrt::Pi + std::atan(y / x);
    else if (x > 0)
        phi = y < 0 ? 2 * pbrt::Pi + std::atan(y / x) : std::atan(y / x);
    else // x == 0
        phi = y < 0 ? 3 * pbrt::PiOver2 : pbrt::PiOver2;

    RealType const r = std::sqrt(x * x + y * y);
    return { r, phi };
}

Можем ли мы как-то уменьшить ошибку решения?