xref: /dports/graphics/appleseed/appleseed-2.1.0-beta/src/appleseed/foundation/math/sampling/mappings.h

//
// This source file is part of appleseed.
// Visit https://appleseedhq.net/ for additional information and resources.
//
// This software is released under the MIT license.
//
// Copyright (c) 2010-2013 Francois Beaune, Jupiter Jazz Limited
// Copyright (c) 2014-2018 Francois Beaune, The appleseedhq Organization
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
// THE SOFTWARE.
//

#pragma once

// appleseed.foundation headers.
#include "foundation/math/scalar.h"
#include "foundation/math/sphericaltriangle.h"
#include "foundation/math/vector.h"

// Standard headers.
#include <cassert>
#include <cmath>
#include <cstddef>

namespace foundation
{

//
// References:
//
//   http://www.cs.kuleuven.ac.be/~phil/GI/TotalCompendium.pdf
//   http://jgt.akpeters.com/papers/ShirleyChiu97/
//   http://psgraphics.blogspot.com/2011/01/improved-code-for-concentric-map.html
//


//
// Sphere sampling functions.
//

// Map a uniform sample in [0,1)^2 to a direction over the unit sphere
// with a uniform probability density p(theta) = 1/(4*Pi).
template <typename T>
Vector<T, 3> sample_sphere_uniform(const Vector<T, 2>& s);


//
// Hemisphere sampling functions.
//

// Map a uniform sample in [0,1)^2 to a direction over the unit hemisphere
// with a uniform probability density p(theta) = 1/(2*Pi).
template <typename T>
Vector<T, 3> sample_hemisphere_uniform(const Vector<T, 2>& s);

// Map a uniform sample in [0,1)^2 to a direction over the unit hemisphere
// with a cosine-weighted probability density p(theta) = cos(theta)/Pi.
template <typename T>
Vector<T, 3> sample_hemisphere_cosine(const Vector<T, 2>& s);

// Map a uniform sample in [0,1)^2 to a direction over the unit hemisphere
// with a cosine lobe probability density p(theta) = (n+1)/(2*Pi)*cos(theta)^n.
template <typename T>
Vector<T, 3> sample_hemisphere_cosine_power(const Vector<T, 2>& s, const T n);
template <typename T>
T sample_hemisphere_cosine_power_pdf(const T cos_theta, const T n);


//
// Disk sampling functions.
//

// Map a uniform sample in [0,1)^2 to a point on the surface of the unit disk
// with a uniform probability density p(x) = 1/Pi.
template <typename T>
Vector<T, 2> sample_disk_uniform(const Vector<T, 2>& s);

// An alternate implementation of foundation::sample_disk_uniform() based on
// a polar parameterization of the unit disk. This implementation is about
// twice as fast as foundation::sample_disk_uniform() but has high distortion
// around the center of the disk. In practice it is fine for random or even
// QMC sampling.
template <typename T>
Vector<T, 2> sample_disk_uniform_alt(const Vector<T, 2>& s);

// Map a uniform random sample in [0,1) to a radius on the unit disk such
// that samples follow a 2D Gaussian distribution with probability density
//   pdf(x, y) = 4*a/Pi*exp(-a*(x^2+y^2))
//   pdf(r, theta) = r*pdf(x, y) = a/Pi*r*exp(-a*r^2)
template <typename T>
T sample_disk_gaussian(
    const T s,
    const T a);
template <typename T>
T sample_disk_gaussian_polar_pdf(
    const T r,
    const T a);
template <typename T>
T sample_disk_gaussian_area_pdf(
    const T r,
    const T a);


//
// Other sampling functions.
//

// Map a uniform sample in [0,1)^2 to a direction in the cone with its apex at
// the origin and extending toward Y+.
template <typename T>
Vector<T, 3> sample_cone_uniform(const Vector<T, 2>& s, const T cos_theta_max);
template <typename T>
T sample_cone_uniform_pdf(const T cos_theta_max);

// Map a uniform sample in [0,1) to a point on the unit circle with a uniform
// probability density p(x) = 1/(2*Pi).
template <typename T>
Vector<T, 2> sample_circle_uniform(const T s);

// Map a uniform sample in [0,1)^2 to a point on the surface of a triangle
// with a uniform probability density p(x) = 1/A. Return the barycentric
// coordinates of the point inside the triangle.
template <typename T>
Vector<T, 3> sample_triangle_uniform(const Vector<T, 2>& s);

// An alternate implementation of sample_triangle_uniform() due to Eric Heitz.
// See https://drive.google.com/file/d/1J-183vt4BrN9wmqItECIjjLIKwm29qSg/view
// for details.
template <typename T>
Vector<T, 3> sample_triangle_uniform_heitz(const Vector<T, 2>& s);

// Build a regular N-sided polygon to use with sample_regular_polygon_uniform().
template <typename T>
void build_regular_polygon(
    const size_t        vertex_count,
    const T             tilt_angle,
    Vector<T, 2>        vertices[]);

// Map a uniform sample in [0,1)^3 to a point on the surface of a regular
// N-sided polygon with a uniform probability density p(x) = 1/A.
template <typename T>
Vector<T, 2> sample_regular_polygon_uniform(
    const Vector<T, 3>& s,
    const size_t        vertex_count,
    const Vector<T, 2>  vertices[]);

// Map a uniform sample in [0,1)^2 to a point on the surface of a regular
// N-sided polygon with a uniform probability density p(x) = 1/A.
template <typename T>
Vector<T, 2> sample_regular_polygon_uniform(
    const Vector<T, 2>& s,
    const size_t        vertex_count,
    const Vector<T, 2>  vertices[]);

// Map a uniform sample in [0,1)^2 to a point on the surface of a spherical
// triangle on the unit sphere.
template <typename T>
Vector<T, 3> sample_spherical_triangle_uniform(
    const Vector<T, 3>& v0,
    const Vector<T, 3>& v1,
    const Vector<T, 3>& v2,
    const Vector<T, 2>& eta);


//
// Exponential sampling functions.
//
// Reference:
//
//   Physically Based Rendering, first edition, page 641.
//

// Map a uniform random sample in [0,1) to an exponential distribution
// of the form exp(-a*x).
template <typename T>
T sample_exponential_distribution(
    const T s,
    const T a);
template <typename T>
T exponential_distribution_pdf(
    const T x,
    const T a);

// Map a uniform random sample in [0,1) to an exponential distribution on segment [l, r].
// of the form exp(-a*(x - l)) - exp(-a*(x - r)).
template <typename T>
T sample_exponential_distribution_on_segment(
    const T s,
    const T a,
    const T l,
    const T r);
template <typename T>
T exponential_distribution_on_segment_pdf(
    const T x,
    const T a,
    const T l,
    const T r);


//
// Equiangular sampling along a ray with respect to a point.
//
// Reference:
//
//   Importance Sampling Techniques for Path Tracing in Participating Media
//   Christopher Kulla and Marcos Fajardo
//   https://www.solidangle.com/research/egsr2012_volume.pdf
//

template <typename T>
T sample_equiangular_distribution(
    const T u,                      // uniform random sample in [0,1)
    const T theta_a,                // start angle
    const T theta_b,                // end angle
    const T distance);              // distance from point to the ray
template <typename T>
T equiangular_distribution_pdf(
    const T t,
    const T theta_a,                // start angle
    const T theta_b,                // end angle
    const T distance);              // distance from point to the ray


//
// Reciprocal distribution:
// p(X) ~ 1 / X, where X belongs to [l, r).
//

template <typename T>
T sample_rcp_distribution(
    const T s,
    const T l,
    const T r);
template <typename T>
T rcp_distribution_pdf(
    const T x,
    const T l,
    const T r);


//
// Implementation.
//

template <typename T>
inline Vector<T, 3> sample_sphere_uniform(const Vector<T, 2>& s)
{
    assert(s[0] >= T(0.0) && s[0] < T(1.0));
    assert(s[1] >= T(0.0) && s[1] < T(1.0));

    const T phi = TwoPi<T>() * s[0];
    const T cos_theta = T(1.0) - T(2.0) * s[1];
    const T sin_theta = std::sqrt(T(1.0) - cos_theta * cos_theta);

    return
        Vector<T, 3>::make_unit_vector(
            cos_theta,
            sin_theta,
            std::cos(phi),
            std::sin(phi));
}

template <typename T>
inline Vector<T, 3> sample_hemisphere_uniform(const Vector<T, 2>& s)
{
    assert(s[0] >= T(0.0) && s[0] < T(1.0));
    assert(s[1] >= T(0.0) && s[1] < T(1.0));

    const T phi = TwoPi<T>() * s[0];
    const T cos_theta = T(1.0) - s[1];
    const T sin_theta = std::sqrt(T(1.0) - cos_theta * cos_theta);

    return
        Vector<T, 3>::make_unit_vector(
            cos_theta,
            sin_theta,
            std::cos(phi),
            std::sin(phi));
}

template <typename T>
inline Vector<T, 3> sample_hemisphere_cosine(const Vector<T, 2>& s)
{
    assert(s[0] >= T(0.0) && s[0] < T(1.0));
    assert(s[1] >= T(0.0) && s[1] < T(1.0));

    const T phi = TwoPi<T>() * s[0];
    const T cos_theta = std::sqrt(T(1.0) - s[1]);
    const T sin_theta = std::sqrt(s[1]);

    return
        Vector<T, 3>::make_unit_vector(
            cos_theta,
            sin_theta,
            std::cos(phi),
            std::sin(phi));
}

template <typename T>
inline Vector<T, 3> sample_hemisphere_cosine_power(const Vector<T, 2>& s, const T n)
{
    assert(s[0] >= T(0.0) && s[0] < T(1.0));
    assert(s[1] >= T(0.0) && s[1] < T(1.0));
    assert(n >= T(0.0));

    const T phi = TwoPi<T>() * s[0];
    const T cos_theta = std::pow(T(1.0) - s[1], T(1.0) / (n + T(1.0)));
    const T sin_theta = std::sqrt(T(1.0) - cos_theta * cos_theta);

    return
        Vector<T, 3>::make_unit_vector(
            cos_theta,
            sin_theta,
            std::cos(phi),
            std::sin(phi));
}

template <typename T>
inline T sample_hemisphere_cosine_power_pdf(const T cos_theta, const T n)
{
    return (n + T(1.0)) * RcpTwoPi<T>() * std::pow(cos_theta, n);
}

template <typename T>
inline Vector<T, 2> sample_disk_uniform(const Vector<T, 2>& s)
{
    assert(s[0] >= T(0.0) && s[0] < T(1.0));
    assert(s[1] >= T(0.0) && s[1] < T(1.0));

    const T a = T(2.0) * s[0] - T(1.0);
    const T b = T(2.0) * s[1] - T(1.0);

    T phi, r;

    if (a * a > b * b)
    {
        r = a;
        phi = PiOverFour<T>() * (b / a);
    }
    else if (b != T(0.0))
    {
        r = b;
        phi = HalfPi<T>() - PiOverFour<T>() * (a / b);
    }
    else return Vector<T, 2>(0.0);

    return Vector<T, 2>(r * std::cos(phi), r * std::sin(phi));
}

template <typename T>
inline Vector<T, 2> sample_disk_uniform_alt(const Vector<T, 2>& s)
{
    assert(s[0] >= T(0.0) && s[0] < T(1.0));
    assert(s[1] >= T(0.0) && s[1] < T(1.0));

    const T phi = TwoPi<T>() * s[0];
    const T r = std::sqrt(T(1.0) - s[1]);

    return r * Vector<T, 2>(std::cos(phi), std::sin(phi));
}

template <typename T>
inline T sample_disk_gaussian(
    const T s,
    const T a)
{
    assert(s >= T(0.0));
    assert(s < T(1.0));

    return std::sqrt(-std::log(T(1.0) - s) / a);
}

template <typename T>
inline T sample_disk_gaussian_polar_pdf(
    const T r,
    const T a)
{
    assert(r >= T(0.0));

    return a * RcpPi<T>() * r * std::exp(-a * r * r);
}

template <typename T>
inline T sample_disk_gaussian_area_pdf(
    const T r,
    const T a)
{
    assert(r >= T(0.0));

    return a * FourOverPi<T>() * std::exp(-a * r * r);
}

template <typename T>
inline Vector<T, 3> sample_cone_uniform(const Vector<T, 2>& s, const T cos_theta_max)
{
    assert(s[0] >= T(0.0) && s[0] < T(1.0));
    assert(s[1] >= T(0.0) && s[1] < T(1.0));

    const T phi = TwoPi<T>() * s[0];
    const T cos_theta = lerp(T(1.0), cos_theta_max, s[1]);
    const T sin_theta = std::sqrt(T(1.0) - cos_theta * cos_theta);

    return
        Vector<T, 3>::make_unit_vector(
            cos_theta,
            sin_theta,
            std::cos(phi),
            std::sin(phi));
}

template <typename T>
inline T sample_cone_uniform_pdf(const T cos_theta_max)
{
    return T(1.0) / (TwoPi<T>() * (T(1.0) - cos_theta_max));
}

template <typename T>
inline Vector<T, 2> sample_circle_uniform(const T s)
{
    assert(s >= T(0.0) && s < T(1.0));

    const T phi = s * TwoPi<T>();

    return Vector<T, 2>(std::cos(phi), std::sin(phi));
}

template <typename T>
inline Vector<T, 3> sample_triangle_uniform(const Vector<T, 2>& s)
{
    assert(s[0] >= T(0.0) && s[0] < T(1.0));
    assert(s[1] >= T(0.0) && s[1] < T(1.0));

    const T sqrt_s0 = std::sqrt(s[0]);

    Vector<T, 3> b;
    b[0] = T(1.0) - sqrt_s0;
    b[1] = (T(1.0) - s[1]) * sqrt_s0;
    b[2] = T(1.0) - b[0] - b[1];

    assert(feq(b[0] + b[1] + b[2], T(1.0)));

    return b;
}

template <typename T>
inline Vector<T, 3> sample_triangle_uniform_heitz(const Vector<T, 2>& s)
{
    assert(s[0] >= T(0.0) && s[0] < T(1.0));
    assert(s[1] >= T(0.0) && s[1] < T(1.0));

    T t0 = T(0.5) * s[0];
    T t1 = T(0.5) * s[1];

    const T offset = t1 - t0;

    if (offset > T(0.0))
        t1 += offset;
    else t0 -= offset;

    Vector<T, 3> b;
    b[0] = t0;
    b[1] = t1;
    b[2] = T(1.0) - b[0] - b[1];

    assert(feq(b[0] + b[1] + b[2], T(1.0)));

    return b;
}

template <typename T>
void build_regular_polygon(
    const size_t        vertex_count,
    const T             tilt_angle,
    Vector<T, 2>        vertices[])
{
    assert(vertex_count >= 3);

    for (size_t i = 0; i < vertex_count; ++i)
    {
        const T a = (i * TwoPi<T>()) / vertex_count + tilt_angle;
        vertices[i] = Vector<T, 2>(std::cos(a), std::sin(a));
    }
}

template <typename T>
inline Vector<T, 2> sample_regular_polygon_uniform(
    const Vector<T, 3>& s,
    const size_t        vertex_count,
    const Vector<T, 2>  vertices[])
{
    assert(s[0] >= T(0.0) && s[0] < T(1.0));
    assert(s[1] >= T(0.0) && s[1] < T(1.0));
    assert(s[2] >= T(0.0) && s[2] < T(1.0));
    assert(vertex_count >= 3);

    const size_t v0_index = truncate<size_t>(s[0] * vertex_count);

    size_t v1_index = v0_index + 1;
    if (v1_index == vertex_count)
        v1_index = 0;

    const Vector<T, 2>& v0 = vertices[v0_index];
    const Vector<T, 2>& v1 = vertices[v1_index];

    T t0 = T(0.5) * s[1];
    T t1 = T(0.5) * s[2];

    const T offset = t1 - t0;

    if (offset > T(0.0))
        t1 += offset;
    else t0 -= offset;

    return t0 * v0 + t1 * v1;
}

template <typename T>
inline Vector<T, 2> sample_regular_polygon_uniform(
    const Vector<T, 2>& s,
    const size_t        vertex_count,
    const Vector<T, 2>  vertices[])
{
    assert(s[0] >= T(0.0) && s[0] < T(1.0));
    assert(s[1] >= T(0.0) && s[1] < T(1.0));
    assert(vertex_count >= 3);

    const size_t v0_index = truncate<size_t>(s[0] * vertex_count);
    const T rescaled_s0 = s[0] * vertex_count - v0_index;
    assert(rescaled_s0 >= T(0.0) && rescaled_s0 < T(1.0));

    size_t v1_index = v0_index + 1;
    if (v1_index == vertex_count)
        v1_index = 0;

    const Vector<T, 2>& v0 = vertices[v0_index];
    const Vector<T, 2>& v1 = vertices[v1_index];

    T t0 = T(0.5) * rescaled_s0;
    T t1 = T(0.5) * s[1];

    const T offset = t1 - t0;

    if (offset > T(0.0))
        t1 += offset;
    else t0 -= offset;

    return t0 * v0 + t1 * v1;
}

template <typename T>
Vector<T, 3> sample_spherical_triangle_uniform(
    const Vector<T, 3>& A,
    const Vector<T, 3>& B,
    const Vector<T, 3>& C,
    const Vector<T, 2>& eta)
{
    // Compute the arc lengths of the sides of the spherical triangle.
    T a, b, c;
    compute_spherical_triangle_edge_lengths(A, B, C, a, b, c);

    // Compute the interior angles of the spherical triangle.
    T alpha, beta, gamma;
    compute_spherical_triangle_interior_angles(a, b, c, alpha, beta, gamma);
    const T cos_alpha = std::cos(alpha);
    const T sin_alpha = std::sin(alpha);

    // Compute the area of the spherical triangle.
    const T area = compute_spherical_triangle_area(alpha, beta, gamma);

    // Use one random variable to select the new area.
    const T area_hat = eta[0] * area;

    // Save the sine and cosine of the angle phi.
    const T s = std::sin(area_hat - alpha);
    const T t = std::cos(area_hat - alpha);

    // Compute the pair (u, v) that determines beta_hat.
    const T u = t - cos_alpha;
    const T v = s + sin_alpha * std::cos(c);

    // Let q be the cosine of the new edge length b_hat.
    const T q_num = (v * t - u * s) * cos_alpha - v;
    const T q_den = (v * s + u * t) * sin_alpha;
    const T q = clamp(q_num / q_den, T(-1.0), T(1.0));

    // Compute the third vertex of the sub-triangle.
    const Vector<T, 3> C_hat = q * A + std::sqrt(T(1.0) - q * q) * normalize(C - dot(C, A) * A);
    const T dot_C_hat_B = dot(C_hat, B);

    // Use the other random variable to select cos(theta).
    const T z = T(1.0) - eta[1] * (T(1.0) - dot_C_hat_B);

    // Construct the corresponding point on the sphere.
    return z * B + std::sqrt(T(1.0) - z * z) * normalize(C_hat - dot_C_hat_B * B);
}

template <typename T>
inline T sample_exponential_distribution(
    const T s,
    const T a)
{
    assert(s >= T(0.0));
    assert(s < T(1.0));

    return -std::log(T(1.0) - s) / a;
}

template <typename T>
inline T exponential_distribution_pdf(
    const T x,
    const T a)
{
    assert(x >= T(0.0));

    return a * std::exp(-a * x);
}

template <typename T>
inline T sample_exponential_distribution_on_segment(
    const T s,
    const T a,
    const T l,
    const T r)
{
    assert(s >= T(0.0));
    assert(s < T(1.0));
    assert(a > T(0.0));

    const T tail = T(1.0) - std::exp(-(r - l) * a);
    return clamp(l - std::log(T(1.0) - s * tail) / a, l, r);
}

template <typename T>
inline T exponential_distribution_on_segment_pdf(
    const T x,
    const T a,
    const T l,
    const T r)
{
    assert(x >= l);
    assert(x <= r);

    const T left = std::exp(a * (x - l));
    const T right = std::exp(a * (x - r));
    return a / (left - right);
}

template <typename T>
inline T sample_equiangular_distribution(
    const T s,
    const T theta_a,
    const T theta_b,
    const T distance)
{
    assert(s >= T(0.0));
    assert(s < T(1.0));

    return distance * std::tan(lerp(theta_a, theta_b, s));
}

template <typename T>
inline T equiangular_distribution_pdf(
    const T t,
    const T theta_a,
    const T theta_b,
    const T distance)
{
    return distance / ((theta_b - theta_a) * (square(distance) + square(t)));
}

template <typename T>
inline T sample_rcp_distribution(
    const T s,
    const T l,
    const T r)
{
    assert(l > T(0.0));
    assert(r > l);

    const T d = r / l;
    return
        d < T(1.01)
            ? l + s * (r - l)
            : clamp(l * std::pow(d, s), l, r);
}

template <typename T>
inline T rcp_distribution_pdf(
    const T x,
    const T l,
    const T r)
{
    assert(l > T(0.0));
    assert(x >= l);
    assert(x < r);

    const T d = r / l;
    return
        d < T(1.01)
            ? rcp(r - l)
            : rcp(x * std::log(d));
}

}   // namespace foundation