Archive for the ‘adaptive multisampling’ tag

Adaptive Multisampling for Area Lights

without comments

I recently added a first pass at area lights in the ray tracer. They are rather limited at the moment as the only point lights are supported and their shape is fixed as a sphere of user-controllable radius r. The sampling is only an approximation and non-uniform, but the results are good enough for now:

Basic Lighting Equation and Shadow Term

The basic lighting equation we’re using is something like this:

$$ I = I_A + \sum S\cdot (I_D + I_S) $$

where $S$ is our shadow term. Prior to adding area lights, the $S$ term was a boolean value that determined whether a light should or shouldn’t contribute light (i.e. energy) to the surface (i.e. “add color” to it). This is basic Computer Graphics 101 Shadows – nothing special.

With area lights, $S$ now represents a “coverage” value from 0.0 to 1.0 representing an approximation of what percentage of the light surface is visible at the point being illuminated. That value is then used to scale the energy of the light passed into the rest of the equation.

The Pseudo-Code

Here’s the basic algorithm being used:

  • If shadows are not enabled, shadow term = 1.0
  • else…
  • If the light radius is zero (i.e. a true point light, not an area light), return the classic boolean 0/1 value for the shadow term
  • else…
  • Create a disc (circular region in 3-space) around the light, oriented toward the intersection point
  • Sample at the disc center and at N points around the circumference
  • If all those samples yield the same value (i.e. 0 or 1), then return that value
  • else…
  • Do many more samples to random locations on the disc representing the light surface and use the mean (i.e. average) value as the shadow term

The Code

A goal of LxEngine is to keep the code base as self-explanatory as possible, so hopefully the code largely speaks for itself:

Adaptive Shadow Term Sampling

float _shadowTerm (const point_light_f& light, const intersection3f& intersection)
{
    if (environment.bShadowsEnabled)
    {
        const float radius   = light.area.radius;
        const float baseTerm = _shadowTermSingle(light.position, intersection);
 
        //
        // Check if the light is an area light, if so take multiple samples
        //
        if (radius > 0.0f)
        {
            //
            // Compute a 3-space disc about the light oriented toward the interesction point
            // 
            const vector3f L    (normalize(light.position - intersection.positionWc));
            const disc3f   disc (light.position, L, radius);
 
            auto sampler = [this, &intersection](const glgeom::point3f& pt) -> float
            {
                return _shadowTermSingle (pt, intersection);
            };
 
            //
            // Take several samples along the circumference of the disc to get an 
            // some sort of guess at the variance.  If the light is not completely
            // visble or completely obscured, then generate far more samples using
            // a random distribution on the disc to come up with a estimate as to
            // what percentage of the disc is visible from the intersection point.
            //
            const size_t kInitial = 6;
            const size_t kFull = 512;
            const float  kEpsilon = 1e-4f;
 
            const float term = sample_disc_circumference<float>(disc, kInitial, sampler);
            float value = glm::mix(baseTerm, term, 1.0f / float(kInitial + 1)); 
            if (value > 1.0f - kEpsilon || value < kEpsilon )
                return value;
            else
                return sample_disc_random<float>(disc, kFull, lx0::random_unit, sampler);
        }
        else
            return baseTerm;
    }
    else
        return 1.0f;
}

Code for Sampling the Circumference

Why do we sample the circumference? The assumption, which certainly isn’t true in the most general case, is that it’s most likely that if the light is partially obscured, one of the boundary points on the light surface will have a different shadow term than some other point on the boundary (or the surface center). Again, that’s not mathematically correct, but we’re assuming it’s accurate enough of the times for the kind of data sets we’re dealing with…

Sampling Along the Circumference of a Disc

GLGeom provides the functions we need to easily generate a set of samples along the circumference:

template <typename T>
T 
sample_disc_circumference (
    const glgeom::disc3t<T>&                     disc,
    size_t                                       samples, 
    std::function<T (const glgeom::point3t<T>&)> sampleFunc)
{
    auto offsets = perpendicular_circular_set(disc.normal, samples);
 
    auto sum = T(0);
    for (auto it = offsets.begin(); it != offsets.end(); ++it)
    {
        sum += sampleFunc(disc.origin + disc.radius * (*it));
    }
    return sum / T(samples);
}

…which in turn uses a function to generate a set of vector orthogonal to a base vector…

Generating a Set of Equally-Spaced Vectors Perpendicular to a Base Vector

template <typename T>
std::vector<vector3t<T>> 
perpendicular_circular_set (const vector3t<T>& w, int N)
{
    typedef vector3t<T> vector3;
    typedef T           scalar;
 
    vector3 u,v;
    perpendicular_axes_smooth(w, u, v);
 
    std::vector<vector3t<T>> results;
    results.reserve(N);
 
    scalar step = glgeom::two_pi().value / N;
    for (int i = 0; i < N; ++i)
    {
        scalar ang = (glgeom::pi().value * i) / N;
        scalar x = cos(ang);
        scalar y = sin(ang);
 
        vector3 p = x * u + y * v;
        results.push_back(p);
    }
    return results;
}

…which in turn uses a function to generate an arbitrary, but consistent and “continuous”, perpendicular vector from the base…

Generating an Continuously-Defined, Arbitrary Basis About a Vector

template <typename T>
vector3t<T>
perpendicular_axis_smooth (const vector3t<T>& w)
{
    vector3t<T> sum;
 
    auto q = abs(w);
    sum += (T(1) - q.x) * cross_with_x(w);
    sum += (T(1) - q.y) * cross_with_y(w);
    sum += (T(1) - q.z) * cross_with_z(w);
 
    return normalize(sum);
}
 
template <typename T>
void
perpendicular_axes_smooth (const vector3t<T>& w, vector3t<T>& u, vector3t<T>& v)
{
    u = perpendicular_axis_smooth(w);
    v = normalize(cross( normalize(w), u ));
}

…and lastly we have the case where we want to randomly sample from the disc…

Randomly Sampling from a Disc

One point worth noting: this is not a uniform sampling from the disc. A uniform sampling would mean that given an infinite number of samples for any given area of the disc, the same number of samples would fall in that area as any other same-sized area within the disc.

Assuming our randomFunc below does return uniform values ranged from $[0,1)$, the below function clearly is not uniform across the disc as the area of the disc varies with $r^2$ and the radius value has a uniform, linear distribution.

Uniform sampling from the disc is being saved for another day. One thing at a time.

template <typename T>
T
sample_disc_random (
    const glgeom::disc3t<T>&                     disc, 
    size_t                                       samples, 
    std::function<T ()>                          randomFunc,
    std::function<T (const glgeom::point3t<T>&)> sampleFunc)
{
    // Create a basis from the normal direction
    glgeom::vector3t<T> u, v;
    perpendicular_axes_smooth(disc.normal, u, v);
 
    //
    // Sample from within the disc
    //
    T sum = T(0);
    for (size_t i = 0; i < samples; ++i)
    {
        // Generate a random point within the disc, then transform to 3-space
        glm::detail::tvec2<T> offsetDisc (randomFunc(), randomFunc());
        offsetDisc = (2 * randomFunc() - 1) * disc.radius * glm::normalize(offsetDisc);
 
        const glgeom::vector3t<T> offsetWs = u * offsetDisc.x + v * offsetDisc.y;
 
        sum += sampleFunc(disc.origin + offsetWs);
    }
    return sum / T(samples);
}