renderer.cpp

Go to the documentation of this file.

/*
    toxic - A Global Illumination Renderer
    Copyright (C) 2003-2004 Francois Beaune
    Contact: http://toxicengine.sourceforge.net/

    This file is part of toxic.

    toxic is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    toxic is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with toxic; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include "renderer.h"   // include first
#include "common/math/customdistribution.h"
#include "common/misc/minmax.h"
#include "context.h"
#include "hit.h"
#include "ibdf.h"
#include "icamera.h"
#include "iedf.h"
#include "iobject.h"
#include "isurfaceshader.h"
#include "photonmap.h"
#include "scene.h"
#include "settings.h"
#include "statistics.h"
#include "surfacebasis.h"
#include "surfacesamplerfactory.h"
#include "utilities.h"

#include "common/misc/pentiumrtsc.h"

#include <cmath>
#include <iostream>

using namespace sheep;
using namespace toxic;

Renderer::Renderer() :
    m_scene(0),
    m_camera(0),
    m_global_pm(0),
    m_caustics_pm(0),
    m_framebuffer(0),
    m_pixel_ss(0),
    m_arealight_ss(0),
    m_primary_fg_hss(0),
    m_secondary_fg_hss(0)
#ifdef ENABLE_WHITTED_SAMPLING_CACHE
    , m_prev_pixels(0)
#endif  // ENABLE_WHITTED_SAMPLING_CACHE
{
}

Renderer::~Renderer() {
    cleanup();
}

void Renderer::SetScene(const Scene *scene) {
    assert(scene);
    m_scene = scene;
}

void Renderer::SetCamera(const ICamera *camera) {
    assert(camera);
    m_camera = camera;
}

void Renderer::SetGlobalPhotonMap(const PhotonMap *global_pm) {
    assert(global_pm);
    m_global_pm = global_pm;
}

void Renderer::SetCausticsPhotonMap(const PhotonMap *caustics_pm) {
    assert(caustics_pm);
    m_caustics_pm = caustics_pm;
}

void Renderer::SetFramebuffer(Framebuffer *framebuffer) {
    assert(framebuffer);
    m_framebuffer = framebuffer;

    ResetRenderArea();
}

void Renderer::SetRenderArea(int x0, int y0, int x1, int y1) {
    assert(m_framebuffer);

    const int xmax = m_framebuffer->GetWidth() - 1;
    const int ymax = m_framebuffer->GetHeight() - 1;

    if(x0 < 0) x0 = 0;
    if(x0 > xmax) x0 = xmax;
    if(y0 < 0) y0 = 0;
    if(y0 > ymax) y0 = ymax;

    if(x1 < 0) x1 = 0;
    if(x1 > xmax) x1 = xmax;
    if(y1 < 0) y1 = 0;
    if(y1 > ymax) y1 = ymax;

    assert(x1 >= x0);
    assert(y1 >= y0);

    m_area_x0 = x0;
    m_area_y0 = y0;
    m_area_x1 = x1;
    m_area_y1 = y1;
}

void Renderer::ResetRenderArea() {
    assert(m_framebuffer);

    m_area_x0 = 0;
    m_area_y0 = 0;
    m_area_x1 = m_framebuffer->GetWidth() - 1;
    m_area_y1 = m_framebuffer->GetHeight() - 1;
}

void Renderer::Restart(const Context &context) {
    cleanup();

    m_x = m_area_x0;
    m_y = m_area_y0;
    m_completed_pixels = 0;

    // Precompute some values.

    m_area_pixels =
        (m_area_x1 - m_area_x0 + 1) *
        (m_area_y1 - m_area_y0 + 1);

    assert(m_area_pixels > 0);
    m_inv_area_pixels = 1.0 / m_area_pixels;

    // Initialize pixel surface sampler.
    m_pixel_ss = 0;
#ifdef ENABLE_WHITTED_SAMPLING_CACHE
    m_prev_pixels = 0;
#endif  // ENABLE_WHITTED_SAMPLING_CACHE
    switch(context.m_settings->m_rendering.m_pixel_sampling.m_algorithm) {
    // Supersampling.
    case Settings::Rendering::PixelSampling::SUPERSAMPLING:
        m_pixel_ss = SurfaceSamplerFactory::Create(
            context.m_settings->m_rendering.m_pixel_sampling.m_supersampling
        );
        break;
    // Whitted adaptive sampling.
    case Settings::Rendering::PixelSampling::WHITTED_ADAPTIVE_SAMPLING:
#ifdef ENABLE_WHITTED_SAMPLING_CACHE
        m_prev_pixels = new Color3[m_framebuffer->GetWidth() + 2];
#endif  // ENABLE_WHITTED_SAMPLING_CACHE
        break;
    default:
        assert(!"Internal failure.");
    }

    if(context.m_settings->m_rendering.m_components.m_direct_lighting.m_enabled) {
        // Initialize area light surface sampler.
        m_arealight_ss = SurfaceSamplerFactory::Create(
            context.m_settings->m_rendering.m_components.m_direct_lighting.m_arealight_sampling
        );

        m_irradiance_samples.reserve(m_arealight_ss->GetSampleCount());
    } else {
        m_arealight_ss = 0;
    }

    if(context.m_settings->m_rendering.m_components.m_indirect_lighting.m_enabled) {
        // Initialize primary final gathering hemisphere sampler.
        m_primary_fg_hss = SurfaceSamplerFactory::Create(
            context.m_settings->m_rendering.m_components.m_indirect_lighting.m_primary_fg
        );

        // Initialize secondary final gathering hemisphere sampler.
        if(context.m_settings->m_rendering.m_components.m_indirect_lighting.m_secondary_fg.m_enabled) {
            m_secondary_fg_hss = SurfaceSamplerFactory::Create(
                context.m_settings->m_rendering.m_components.m_indirect_lighting.m_secondary_fg
            );
        } else {
            m_secondary_fg_hss = 0;
        }
    } else {
        m_primary_fg_hss = 0;
        m_secondary_fg_hss = 0;
    }
}

void Renderer::RenderNextPixel(const Context &context) {
    assert(m_scene);
    assert(m_camera);
    assert(m_framebuffer);

    if(!IsRenderComplete()) {
        render_pixel(context);

        ++context.m_statistics->m_total_pixels;

        ++m_completed_pixels;

        if(m_x < m_area_x1)
            ++m_x;
        else {
            m_x = m_area_x0;
            ++m_y;
        }
    }
}

void Renderer::Annotate(const Context &context) {
    m_scene->GetRootObject()->Annotate(context, m_camera, m_framebuffer);
}

void Renderer::cleanup() {
    delete m_pixel_ss;
    delete m_arealight_ss;
    delete m_primary_fg_hss;
    delete m_secondary_fg_hss;

#ifdef ENABLE_WHITTED_SAMPLING_CACHE
    delete [] m_prev_pixels;
#endif  // ENABLE_WHITTED_SAMPLING_CACHE
}

void Renderer::render_pixel(const Context &context) {
/*  static PentiumRTSC timer;

    timer.Start();*/

    switch(context.m_settings->m_rendering.m_pixel_sampling.m_algorithm) {
    // Supersampling.
    case Settings::Rendering::PixelSampling::SUPERSAMPLING:
        render_pixel_supersampling(context);
        break;
    // Whitted adaptive sampling.
    case Settings::Rendering::PixelSampling::WHITTED_ADAPTIVE_SAMPLING:
        render_pixel_whitted(context);
        break;
    default:
        assert(!"Internal failure.");
    }

/*  timer.Stop();

    m_framebuffer->SetPixel(
        m_x, m_y,
        Color3(static_cast<Real>(timer.GetResult()))
    );*/
}

void Renderer::render_pixel_supersampling(const Context &context) {
    const Point2 pixel_center = m_framebuffer->ConvertToNDC(m_x, m_y);
    const Real pixel_width = m_framebuffer->GetPixelWidth();
    const Real pixel_height = m_framebuffer->GetPixelHeight();

    m_pixel_ss->GenerateNewSamples(context);

    for(ISurfaceSampler::SampleVector::iterator i = m_pixel_ss->m_samples.begin(),
        e = m_pixel_ss->m_samples.end(); i != e; ++i)
    {
        i->m_x = pixel_center.m_x + (i->m_x - 0.5) * pixel_width;
        i->m_y = pixel_center.m_y + (i->m_y - 0.5) * pixel_height;
    }

    Color3 radiance(0.0);

    for(ISurfaceSampler::SampleVector::const_iterator i = m_pixel_ss->m_samples.begin(),
        e = m_pixel_ss->m_samples.end(); i != e; ++i)
    {
        // Trace a ray through the pixel to compute the radiance received from this direction.
        const Ray ray = m_camera->ComputeRay(context, *i);
        radiance += trace(context, ray, 0);
    }

    radiance *= m_pixel_ss->GetInvSampleCount();

    // Store the radiance into the framebuffer.
    m_framebuffer->SetPixel(m_x, m_y, radiance);
}

void Renderer::render_pixel_whitted(const Context &context) {
    const Point2 pixel_center = m_framebuffer->ConvertToNDC(m_x, m_y);

    const Real half_pix_w = 0.5 * m_framebuffer->GetPixelWidth();
    const Real half_pix_h = 0.5 * m_framebuffer->GetPixelHeight();

    const Point2 v0(pixel_center.m_x - half_pix_w, pixel_center.m_y + half_pix_h);
    const Point2 v1(pixel_center.m_x + half_pix_w, pixel_center.m_y + half_pix_h);
    const Point2 v2(pixel_center.m_x + half_pix_w, pixel_center.m_y - half_pix_h);
    const Point2 v3(pixel_center.m_x - half_pix_w, pixel_center.m_y - half_pix_h);

#ifdef ENABLE_WHITTED_SAMPLING_CACHE
    Color3 r0, r1, r2, r3;

    if(m_y == m_area_y0) {
        if(m_x == m_area_x0) {
            // First line, first column.

            r0 = trace(context, m_camera->ComputeRay(context, v0), 0);
            r1 = trace(context, m_camera->ComputeRay(context, v1), 0);
            r2 = trace(context, m_camera->ComputeRay(context, v2), 0);
            r3 = trace(context, m_camera->ComputeRay(context, v3), 0);

            m_prev_pixels[0] = r1;
            m_prev_pixels[m_x + 1] = r3;
            m_prev_pixels[m_x + 2] = r2;
        } else {
            // First line, any column.

            r0 = m_prev_pixels[0];
            r1 = trace(context, m_camera->ComputeRay(context, v1), 0);
            r2 = trace(context, m_camera->ComputeRay(context, v2), 0);
            r3 = m_prev_pixels[m_x + 1];

            m_prev_pixels[0] = r1;
            m_prev_pixels[m_x + 2] = r2;
        }
    } else {
        if(m_x == m_area_x0) {
            // Any line, first column.

            r0 = m_prev_pixels[m_x + 1];
            r1 = m_prev_pixels[m_x + 2];
            r2 = trace(context, m_camera->ComputeRay(context, v2), 0);
            r3 = trace(context, m_camera->ComputeRay(context, v3), 0);

            m_prev_pixels[0] = r1;
            m_prev_pixels[m_x + 1] = r3;
            m_prev_pixels[m_x + 2] = r2;
        } else {
            // Any line, any column (the most common case).

            r0 = m_prev_pixels[0];
            r1 = m_prev_pixels[m_x + 2];
            r2 = trace(context, m_camera->ComputeRay(context, v2), 0);
            r3 = m_prev_pixels[m_x + 1];

            m_prev_pixels[0] = r1;
            m_prev_pixels[m_x + 2] = r2;
        }
    }
#else
    const Color3 r0 = trace(context, m_camera->ComputeRay(context, v0), 0);
    const Color3 r1 = trace(context, m_camera->ComputeRay(context, v1), 0);
    const Color3 r2 = trace(context, m_camera->ComputeRay(context, v2), 0);
    const Color3 r3 = trace(context, m_camera->ComputeRay(context, v3), 0);
#endif  // ENABLE_WHITTED_SAMPLING_CACHE

    const Color3 radiance = whitted_recursive_sampling(
        context,
        pixel_center,
        v0, v1, v2, v3,
        r0, r1, r2, r3,
        0   // recursion level
    );

    // Store the radiance into the framebuffer.
    m_framebuffer->SetPixel(m_x, m_y, radiance);
}

namespace {
    inline
    Real contrast(const Color3 &c0, const Color3 &c1) {
        return fabs(c0.Luminance() - c1.Luminance());
    }
}

Color3 Renderer::whitted_recursive_sampling(const Context &context,
                                            const Point2 &center,
                                            const Point2 &v0,
                                            const Point2 &v1,
                                            const Point2 &v2,
                                            const Point2 &v3,
                                            const Color3 &r0,
                                            const Color3 &r1,
                                            const Color3 &r2,
                                            const Color3 &r3,
                                            int recursion_level)
{
    assert(recursion_level >= 0);
    assert(recursion_level <=
        context.m_settings->m_rendering.m_pixel_sampling.m_whitted_adaptive_sampling.m_max_depth);

    const Real contrast_threshold =
        context.m_settings->m_rendering.m_pixel_sampling.m_whitted_adaptive_sampling.m_contrast_threshold;
    const int max_depth =
        context.m_settings->m_rendering.m_pixel_sampling.m_whitted_adaptive_sampling.m_max_depth;

    if(recursion_level < max_depth &&
       (contrast(r0, r1) > contrast_threshold ||
        contrast(r1, r2) > contrast_threshold ||
        contrast(r2, r3) > contrast_threshold ||
        contrast(r3, r0) > contrast_threshold ||
        contrast(r0, r2) > contrast_threshold ||
        contrast(r1, r3) > contrast_threshold))
    {
        const Point2 v4(center.m_x, v0.m_y);
        const Point2 v5(v1.m_x, center.m_y);
        const Point2 v6(center.m_x, v2.m_y);
        const Point2 v7(v3.m_x, center.m_y);

        const Color3 rc = trace(context, m_camera->ComputeRay(context, center), 0);
        const Color3 r4 = trace(context, m_camera->ComputeRay(context, v4), 0);
        const Color3 r5 = trace(context, m_camera->ComputeRay(context, v5), 0);
        const Color3 r6 = trace(context, m_camera->ComputeRay(context, v6), 0);
        const Color3 r7 = trace(context, m_camera->ComputeRay(context, v7), 0);

        const Real half_subpix_w = 0.5 * (center.m_x - v0.m_x);
        const Real half_subpix_h = 0.5 * (v0.m_y - center.m_y);

        assert(half_subpix_w > 0.0);
        assert(half_subpix_h > 0.0);

        const int new_recursion_level = recursion_level + 1;

        const Color3 spr0 = whitted_recursive_sampling(
            context,
            Point2(v0.m_x + half_subpix_w, v0.m_y - half_subpix_h),
            v0, v4, center, v7,
            r0, r4, rc, r7,
            new_recursion_level
        );

        const Color3 spr1 = whitted_recursive_sampling(
            context,
            Point2(v1.m_x - half_subpix_w, v1.m_y - half_subpix_h),
            v4, v1, v5, center,
            r4, r1, r5, rc,
            new_recursion_level
        );

        const Color3 spr2 = whitted_recursive_sampling(
            context,
            Point2(v2.m_x - half_subpix_w, v2.m_y + half_subpix_h),
            center, v5, v2, v6,
            rc, r5, r2, r6,
            new_recursion_level
        );

        const Color3 spr3 = whitted_recursive_sampling(
            context,
            Point2(v3.m_x + half_subpix_w, v3.m_y + half_subpix_h),
            v7, center, v6, v3,
            r7, rc, r6, r3,
            new_recursion_level
        );

        return (1.0 / 4.0) * (spr0 + spr1 + spr2 + spr3);
    } else return (1.0 / 4.0) * (r0 + r1 + r2 + r3);
}

Color3 Renderer::trace(const Context &context, const Ray &ray, int recursion_level) {
    assert(ray.m_direction.IsUnitLength());

#ifdef _DEBUG
    if(context.m_settings->m_rendering.m_components.m_specular_reflections.m_enabled) {
        assert(recursion_level >= 0);
        assert(recursion_level <=
            context.m_settings->m_rendering.m_components.m_specular_reflections.m_max_depth);
    } else assert(recursion_level == 0);
#endif  // _DEBUG

    switch(recursion_level) {
    case 0:
        ++context.m_statistics->m_primary_rays;
        break;
    default:
        ++context.m_statistics->m_secondary_rays;
    }

    Hit hit;

    if(!m_scene->Trace(context, ray, &hit)) {
        // The ray did not hit any object: return the background color.
        return m_scene->GetBackgroundColor();
    }

    const Vector3 outgoing = -ray.m_direction;

    // Extract intersection data.
    Point3 point;
    Vector3 geometric_normal, shading_normal;
    hit.ExtractIntersection(ray, &point, &geometric_normal, &shading_normal);
    const ISurfaceShader *surface_shader = hit.m_object->GetSurfaceShader();

    ShadingData shadingdata;
    surface_shader->Shade(hit, &shadingdata);

// False colors rendering.
//srand(int(hit.m_object));
//return Color3(Real(rand() % 256) / 255.0, Real(rand() % 256) / 255.0, Real(rand() % 256) / 255.0);

    Color3 radiance(0.0);

    // Add direct lighting.
    if(context.m_settings->m_rendering.m_components.m_direct_lighting.m_enabled) {
        //!\todo Add an option to enable or disable direct incoming light.
        if(surface_shader->GetEDF()) {
            //!\todo This is a hack.
            radiance +=
                surface_shader->GetRadiantExitance() / surface_shader->GetRadiantExitance().Average();
//          radiance +=
//              surface_shader->GetEDF()->GetEmittedRadiance(context, normal, outgoing) *
//              surface_shader->GetRadiantExitance();
        }

        radiance +=
            compute_direct_illumination(
                context,
                point,
                geometric_normal,
                shading_normal,
                outgoing,
                surface_shader,
                shadingdata
            );
    }

    // Add indirect lighting.
    if(context.m_settings->m_rendering.m_components.m_indirect_lighting.m_enabled) {
        if(context.m_settings->m_rendering.m_components.m_indirect_lighting.m_primary_fg.m_enabled) {
            radiance +=
                final_gathering(
                    context,
                    point,
                    geometric_normal,
                    shading_normal,
                    outgoing,
                    shadingdata
                );
        } else {
            radiance +=
                compute_indirect_illumination(
                    context,
                    point,
                    geometric_normal,
                    shading_normal,
                    outgoing,
                    shadingdata
                );
        }
    }

    // Add specular reflections.
    if(context.m_settings->m_rendering.m_components.m_specular_reflections.m_enabled) {
        if(recursion_level < context.m_settings->m_rendering.m_components.m_specular_reflections.m_max_depth) {
            radiance +=
                compute_specular_reflections(
                    context,
                    point,
                    geometric_normal,
                    shading_normal,
                    outgoing,
                    shadingdata,
                    recursion_level
                );
        }
    }

    // Add caustics.
    if(context.m_settings->m_rendering.m_components.m_caustics.m_enabled) {
        radiance +=
            compute_caustics(
                context,
                point,
                geometric_normal,
                shading_normal,
                outgoing,
                shadingdata
            );
    }

    return radiance;
}

Color3 Renderer::compute_direct_illumination(const Context &context,
                                             const Point3 &point,
                                             const Vector3 &geometric_normal,
                                             const Vector3 &shading_normal,
                                             const Vector3 &outgoing,
                                             const ISurfaceShader *surface_shader,
                                             const ShadingData &shadingdata)
{
    assert(geometric_normal.IsUnitLength());
    assert(shading_normal.IsUnitLength());
    assert(outgoing.IsUnitLength());
    assert(surface_shader);

    if(!shadingdata.m_bdf->IsDiffuse() || shadingdata.m_reflectance.IsBlack())
        return Color3(0.0);

    // Build an orthonormal basis around the normal vector.
    //!\todo Move to caller.
    const SurfaceBasis surfacebasis(shading_normal);
    const Vector3 local_outgoing = surfacebasis.TransformToLocal(outgoing);

    // Start with ambient illumination.
    Color3 total_radiance = m_scene->GetAmbientIllumination();

    //!\todo Choose lights based on their occlusion-free contribution to the
    //! point being illuminated.
    for(int i = 0; i < m_scene->GetLightCount(); ++i) {
        const ILight *light = m_scene->GetLight(i);

        //!\todo Avoid generating samples for point lights.
        m_irradiance_samples.clear();
        m_arealight_ss->GenerateNewSamples(context);

        // Compute irradiance due to this light source.
        light->ComputeIrradiance(
            context,
            m_scene,
            point,
            geometric_normal,
            shading_normal,
            m_arealight_ss->m_samples,
            &m_irradiance_samples
        );

        assert(m_irradiance_samples.size() <= m_arealight_ss->m_samples.size());

        Color3 radiance(0.0);

        for(ILight::IrradianceSampleVector::const_iterator i = m_irradiance_samples.begin(),
            e = m_irradiance_samples.end(); i != e; ++i)
        {
            const ILight::IrradianceSample &sample = *i;

            if(sample.m_irradiance.IsBlack())
                continue;

            // Evaluate the BSDF. 'bsdf_value' is expressed in sr^-1.
            const Real bsdf_value =
                shadingdata.m_bdf->Evaluate(
                    context,
                    surfacebasis.TransformToLocal(sample.m_direction),
                    local_outgoing
                );

            // Compute and accumulate the reflected radiance due to this light source.
            radiance += bsdf_value * sample.m_irradiance;
        }

        if(m_irradiance_samples.size() > 1)
            radiance /= m_irradiance_samples.size();    //!\todo Precompute the division.

        total_radiance += radiance;
    }

    return shadingdata.m_reflectance * total_radiance;
}

Color3 Renderer::compute_specular_reflections(const Context &context,
                                              const Point3 &point,
                                              const Vector3 &geometric_normal,
                                              const Vector3 &shading_normal,
                                              const Vector3 &outgoing,
                                              const ShadingData &shadingdata,
                                              int recursion_level)
{
    assert(geometric_normal.IsUnitLength());
    assert(shading_normal.IsUnitLength());
    assert(outgoing.IsUnitLength());
    assert(recursion_level >= 0);
    assert(recursion_level <
        context.m_settings->m_rendering.m_components.m_specular_reflections.m_max_depth);

    if(!shadingdata.m_bdf->IsSpecular() || shadingdata.m_reflectance.IsBlack())
        return Color3(0.0);

    // Build an orthonormal basis around the normal vector.
    //!\todo Move to caller.
    const SurfaceBasis surfacebasis(shading_normal);
    const Vector3 local_outgoing = surfacebasis.TransformToLocal(outgoing);

    Vector3 scattereddir;

    const Real bsdf_value =
        shadingdata.m_bdf->EvaluateSpecular(
            context,
            local_outgoing,
            &scattereddir
        );

    if(bsdf_value == 0.0)
        return Color3(0.0);

    // Reflected rays will be casted from a point that is slighty shifted along
    // the geometric normal to the surface at the intersection point. This
    // is necessary to avoid false intersections.
    const Point3 shifted_point = point + 1.0e-8 * geometric_normal;

    const Color3 radiance =
        trace(
            context,
            Ray(Ray::REFLECTED_RAY, shifted_point, surfacebasis.TransformToWorld(scattereddir)),
            recursion_level + 1
        );

    //!\todo Warning: units inconsistency.
    return bsdf_value * shadingdata.m_reflectance * radiance;
}

Color3 Renderer::final_gathering(const Context &context,
                                 const Point3 &point,
                                 const Vector3 &geometric_normal,
                                 const Vector3 &shading_normal,
                                 const Vector3 &outgoing,
                                 const ShadingData &shadingdata,
                                 bool is_secondary /*= false*/) const
{
    assert(geometric_normal.IsUnitLength());
    assert(shading_normal.IsUnitLength());
    assert(outgoing.IsUnitLength());

//  if(is_secondary)
//      return Color3(100.0, 0.0, 100.0);

    // Build an orthonormal basis around the normal vector.
    //!\todo Move to caller.
    const SurfaceBasis surfacebasis(shading_normal);

    Ray reflected_ray(is_secondary ? Ray::SECONDARY_FG_RAY : Ray::PRIMARY_FG_RAY);
    reflected_ray.m_origin = point + 1.0e-8 * geometric_normal;

    ISurfaceSampler::SampleVector::const_iterator i, e;

    if(is_secondary) {
        m_secondary_fg_hss->GenerateNewSamples(context);

        i = m_secondary_fg_hss->m_samples.begin();
        e = m_secondary_fg_hss->m_samples.end();

        context.m_statistics->m_secondary_final_gathering_rays +=
            m_secondary_fg_hss->m_samples.size();
    } else {
        m_primary_fg_hss->GenerateNewSamples(context);

        i = m_primary_fg_hss->m_samples.begin();
        e = m_primary_fg_hss->m_samples.end();

        context.m_statistics->m_primary_final_gathering_rays +=
            m_primary_fg_hss->m_samples.size();
    }

    Color3 radiance(0.0);
    int samples = 0;

    for(; i != e; ++i) {
        // Compute the direction of the reflected ray in local space.
        const Vector3 direction = CosineHemisphereSampling(i->m_x, i->m_y);

        // Compute the direction of the reflected ray in world space.
        reflected_ray.NewId();
        reflected_ray.m_direction = surfacebasis.TransformToWorld(direction);
        assert(reflected_ray.m_direction.IsUnitLength());

        // Trace the ray.
        Hit hit;
        if(m_scene->Trace(context, reflected_ray, &hit))
            ++samples;
        else continue;  // next cell

        // Extract intersection data.
        Point3 point2;
        Vector3 geometric_normal2, shading_normal2;
        hit.ExtractIntersection(reflected_ray, &point2, &geometric_normal2, &shading_normal2);
        const ISurfaceShader *surface_shader2 = hit.m_object->GetSurfaceShader();

        ShadingData shadingdata2;
        surface_shader2->Shade(hit, &shadingdata2);

        if( !is_secondary &&
            context.m_settings->m_rendering.m_components.m_indirect_lighting.m_secondary_fg.m_enabled &&
            hit.m_abscissa <
            context.m_settings->m_rendering.m_components.m_indirect_lighting.m_secondary_fg.m_distance_threshold)
        {
            // Distance along this final gathering ray is below the threshold:
            // we must perform a secondary, simpler final gathering.
            radiance +=
                final_gathering(
                    context,
                    point2,
                    geometric_normal2,
                    shading_normal2,
                    -reflected_ray.m_direction,
                    shadingdata2,
                    true    // secondary final gathering
                );
        } else {
            radiance +=
                compute_indirect_illumination(
                    context,
                    point2,
                    geometric_normal2,
                    shading_normal2,
                    -reflected_ray.m_direction,
                    shadingdata2
                );
        }
    }

    if(samples > 1)
        radiance /= samples;

    return shadingdata.m_reflectance * radiance;
}

Color3 Renderer::compute_caustics(const Context &context,
                                  const Point3 &point,
                                  const Vector3 &geometric_normal,
                                  const Vector3 &shading_normal,
                                  const Vector3 &outgoing,
                                  const ShadingData &shadingdata) const
{
    assert(geometric_normal.IsUnitLength());
    assert(shading_normal.IsUnitLength());
    assert(outgoing.IsUnitLength());

    if(!shadingdata.m_bdf->IsDiffuse())
        return Color3(0.0);

    return m_caustics_pm->ComputeRadiance(
        context,
        point,
        geometric_normal,
        outgoing,
        shadingdata,
        context.m_settings->m_rendering.m_components.m_caustics.m_radiance_est.m_max_distance,
        context.m_settings->m_rendering.m_components.m_caustics.m_radiance_est.m_max_photons
    );
}

Color3 Renderer::compute_indirect_illumination(const Context &context,
                                               const Point3 &point,
                                               const Vector3 &geometric_normal,
                                               const Vector3 &shading_normal,
                                               const Vector3 &outgoing,
                                               const ShadingData &shadingdata) const
{
    assert(geometric_normal.IsUnitLength());
    assert(shading_normal.IsUnitLength());
    assert(outgoing.IsUnitLength());

    if(context.m_settings->m_rendering.m_components.m_indirect_lighting.m_radiance_precomp.m_enabled) {
        return
            m_global_pm->GetPrecomputedRadiance(
                point,
                geometric_normal,
                shadingdata,
                context.m_settings->m_rendering.m_components.m_indirect_lighting.m_radiance_precomp.m_max_search_dist
            );
    } else {
        //!\todo Let the user choose between exact and approximate radiance computation.
        return
            m_global_pm->ComputeRadiance(
                context,
                point,
                geometric_normal,
                outgoing,
                shadingdata,
                context.m_settings->m_rendering.m_components.m_indirect_lighting.m_radiance_est.m_max_distance,
                context.m_settings->m_rendering.m_components.m_indirect_lighting.m_radiance_est.m_max_photons
            );
    }
}

---