photonmap.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 "photonmap.h"  // include first
#include "common/misc/binarystream.h"
#include "common/misc/minmax.h"
#include "context.h"
#include "ibdf.h"
#include "isurfaceshader.h"
#include "surfacebasis.h"

#include <algorithm>    // std::swap<>(), std::min<>()
#include <cmath>
#include <cstdio>

using namespace sheep;
using namespace std;
using namespace toxic;

//! Signature string identifying photon map files.
//! Increment version number each time the file structure change.
const string PhotonMap::m_photon_map_sig = "toxic photon map file version 1";

//! Buffer size for photon map I/O, expressed in bytes.
const int PHOTONMAP_IO_BUFFER_SIZE = 1024 * 1024;

PhotonMap::PhotonMap() :
    m_stored_photons(0),
    m_half_stored_photons(0),
    m_prev_scale(1)
{
    m_photons.push_back(Photon());

    m_bbox_min[0] = m_bbox_min[1] = m_bbox_min[2] = 1.0e8f;
    m_bbox_max[0] = m_bbox_max[1] = m_bbox_max[2] = -1.0e8f;

    // Allocate the buffers used by IrradianceEstimate(), FastIrradianceEstimate()
    // and BuildPhotonDirectionHistogram() methods.
    m_dist2_buf = new float32[BUF_SIZE];
    m_index_buf = new const Photon *[BUF_SIZE];
}

PhotonMap::~PhotonMap() {
    delete [] m_index_buf;
    delete [] m_dist2_buf;
}

int PhotonMap::GetSizeInMemory() const {
    return
        sizeof(PhotonMap) +
        sizeof(Photon) * static_cast<int>(m_photons.capacity()) +
        sizeof(float32) * BUF_SIZE +
        sizeof(Photon *) * BUF_SIZE;
}

void PhotonMap::StorePhoton(const Point3 &position,
                            const Color3 &power,
                            const Vector3 &direction
#ifdef ENABLE_RADIANCES_PRECOMPUTATION
                            , const Vector3 &normal
#endif  // ENABLE_RADIANCES_PRECOMPUTATION
                            )
{
    Photon photon;

    photon.SetPosition(position);
    photon.SetPower(power);
    photon.SetIncidentDirection(direction);

#ifdef ENABLE_RADIANCES_PRECOMPUTATION
    photon.SetSurfaceNormal(normal);
#endif  // ENABLE_RADIANCES_PRECOMPUTATION

    m_photons.push_back(photon);

    ++m_stored_photons;

    // Update the bounding box.
    if(position.m_x < m_bbox_min[0]) m_bbox_min[0] = static_cast<float32>(position.m_x);
    if(position.m_x > m_bbox_max[0]) m_bbox_max[0] = static_cast<float32>(position.m_x);
    if(position.m_y < m_bbox_min[1]) m_bbox_min[1] = static_cast<float32>(position.m_y);
    if(position.m_y > m_bbox_max[1]) m_bbox_max[1] = static_cast<float32>(position.m_y);
    if(position.m_z < m_bbox_min[2]) m_bbox_min[2] = static_cast<float32>(position.m_z);
    if(position.m_z > m_bbox_max[2]) m_bbox_max[2] = static_cast<float32>(position.m_z);
}

void PhotonMap::ScalePhotonPower(Real scale) {
    for(int i = m_prev_scale; i <= m_stored_photons; ++i)
        m_photons[i].ScalePower(scale);

    m_prev_scale = m_stored_photons + 1;
}

void PhotonMap::Balance() {
    if(m_stored_photons > 1) {
        // Allocate two temporary arrays for the balancing procedure.
        Photon **pa1 = new Photon *[m_stored_photons + 1];
        Photon **pa2 = new Photon *[m_stored_photons + 1];

        for(int i = 0; i <= m_stored_photons; ++i)
            pa2[i] = &m_photons[i];

        balance_segment(pa1, pa2, 1, 1, m_stored_photons);

        delete [] pa2;

        // Reorganize balanced kd-tree (make a heap).
        int j = 1, foo = 1;
        Photon foo_photon = m_photons[j];

        for(int i = 1; i <= m_stored_photons; ++i) {
            //!\todo Warning: this is not portable.
            const int d = pa1[j] - &m_photons[0];

            pa1[j] = 0;

            if(d != foo)
                m_photons[j] = m_photons[d];
            else {
                m_photons[j] = foo_photon;

                if(i < m_stored_photons) {
                    for(; foo <= m_stored_photons; ++foo)
                        if(pa1[foo] != 0)
                            break;
                    foo_photon = m_photons[foo];
                    j = foo;
                }

                continue;
            }

            j = d;
        }

        delete [] pa1;
    }

    m_half_stored_photons = m_stored_photons / 2 - 1;
}

Color3 PhotonMap::ComputeRadiance(const Context &context,
                                  const Point3 &point,
                                  const Vector3 &normal,
                                  const Vector3 &outgoing,
                                  const ShadingData &shadingdata,
                                  Real max_dist,
                                  int max_photons) const
{
    assert(max_photons + 1 <= BUF_SIZE);

    nearest_photons np;
    np.m_dist2 = m_dist2_buf;
    np.m_index = m_index_buf;
    np.m_position[0] = static_cast<float32>(point.m_x);
    np.m_position[1] = static_cast<float32>(point.m_y);
    np.m_position[2] = static_cast<float32>(point.m_z);
    np.m_max = max_photons;
    np.m_found = 0;
    np.m_got_heap = 0;
    np.m_dist2[0] = static_cast<float32>(max_dist * max_dist);

    // Locate the nearest photons.
    locate_photons(&np, 1);

    // If less than 8 photons return.
    if(np.m_found < 8)
        return Color3(0.0);

    const SurfaceBasis surfacebasis(normal);
    const Vector3 local_outgoing = surfacebasis.TransformToLocal(outgoing);

    Color3 radiant_intensity(0.0);

    for(int i = 1; i <= np.m_found; ++i) {
        const Photon *p = np.m_index[i];

#ifdef CHECK_FOR_THIN_SURFACES
        if(p->GetIncidentDirection() * normal <= 0.0)
            continue;
#endif  // CHECK_FOR_THIN_SURFACES

#ifdef ENABLE_RADIANCES_PRECOMPUTATION
        // In some cases, this test significantly improves the quality of
        // the indirect lighting in corners and along edges.
        //!\todo Let THRESHOLD be a user setting.
        const Real THRESHOLD = 0.9;
        if(p->GetSurfaceNormal() * normal < THRESHOLD)
            continue;
#endif  // ENABLE_RADIANCES_PRECOMPUTATION

        // Evaluate the local BSDF at this photon's location.
        const Real bsdf_value =
            shadingdata.m_bdf->Evaluate(
                context,
                surfacebasis.TransformToLocal(p->GetIncidentDirection()),
                local_outgoing);

        // Accumulate radiant intensity.
        //!\todo Missing cosine factor?
        radiant_intensity += bsdf_value * p->GetPower();
    }

    // Compute the area of the disc defined by the intersection of the nearest
    // photons bounding sphere and the surface (locally approximated by a plane)
    // containing the point at which we evaluate the irradiance.
    assert(fnz(np.m_dist2[0]));
    const Real inv_disc_area = (1.0 / PI) / np.m_dist2[0];

    // Compute radiance from radiant intensity.
    const Color3 radiance = radiant_intensity * inv_disc_area;

    return shadingdata.m_reflectance * radiance;
}

Color3 PhotonMap::ApproximateRadiance(const Point3 &point,
                                      const Vector3 &normal,
                                      const ShadingData &shadingdata,
                                      Real max_dist,
                                      int max_photons) const
{
    assert(max_photons + 1 <= BUF_SIZE);

    nearest_photons np;
    np.m_dist2 = m_dist2_buf;
    np.m_index = m_index_buf;
    np.m_position[0] = static_cast<float32>(point.m_x);
    np.m_position[1] = static_cast<float32>(point.m_y);
    np.m_position[2] = static_cast<float32>(point.m_z);
    np.m_max = max_photons;
    np.m_found = 0;
    np.m_got_heap = 0;
    np.m_dist2[0] = static_cast<float32>(max_dist * max_dist);

    // Locate the nearest photons.
    locate_photons(&np, 1);

    // If less than 8 photons return.
    if(np.m_found < 8)
        return Color3(0.0);

    Color3 radiant_intensity(0.0);

    for(int i = 1; i <= np.m_found; ++i) {
        const Photon *p = np.m_index[i];

#ifdef CHECK_FOR_THIN_SURFACES
        if(p->GetIncidentDirection() * normal <= 0.0)
            continue;
#endif  // CHECK_FOR_THIN_SURFACES

#ifdef ENABLE_RADIANCES_PRECOMPUTATION
        // In some cases, this test significantly improves the quality of
        // the indirect lighting in corners and along edges.
        //!\todo Let THRESHOLD be a user setting.
        const Real THRESHOLD = 0.9;
        if(p->GetSurfaceNormal() * normal < THRESHOLD)
            continue;
#endif  // ENABLE_RADIANCES_PRECOMPUTATION

        // Approximate the local BSDF at this photon's location by the Lambertian BRDF.
        //!\todo Move out of the loop.
        const Real bsdf_value = (1.0 / PI);

        // Accumulate radiant intensity.
        radiant_intensity += bsdf_value * p->GetPower();
    }

    // Compute the area of the disc defined by the intersection of the nearest
    // photons bounding sphere and the surface (locally approximated by a plane)
    // containing the point at which we evaluate the irradiance.
    assert(fnz(np.m_dist2[0]));
    const Real inv_disc_area = (1.0 / PI) / np.m_dist2[0];

    // Compute radiance from radiant intensity.
    const Color3 radiance = radiant_intensity * inv_disc_area;

    return shadingdata.m_reflectance * radiance;
}

#ifdef ENABLE_RADIANCES_PRECOMPUTATION

void PhotonMap::PrecomputeRadiances(int spacing,
                                    Real max_dist,
                                    int max_photons,
                                    ProgressMonitor *progmon /*= 0*/)
{
    assert(spacing >= 1);
    assert(max_photons > 0);

    if(progmon && m_stored_photons > 0)
        progmon->StartJob(0.0, m_stored_photons);

    ShadingData fakeshadingdata;
    fakeshadingdata.m_reflectance = Color3(1.0);

    for(int i = 1; i <= m_stored_photons; i += spacing) {
        if(progmon)
            progmon->SetJobProgress(i);

        Photon *photon = &m_photons[i];

        const Color3 radiance =
            ApproximateRadiance(
                photon->GetPosition(),
                photon->GetSurfaceNormal(),
                fakeshadingdata,
                max_dist,
                max_photons);

        photon->SetPrecomputedRadiance(radiance);
    }

    if(progmon)
        progmon->Done();
}

Color3 PhotonMap::GetPrecomputedRadiance(const Point3 &point,
                                         const Vector3 &normal,
                                         const ShadingData &shadingdata,
                                         Real max_dist) const
{
    nearest_photons np;
    np.m_dist2 = m_dist2_buf;
    np.m_index = m_index_buf;
    np.m_position[0] = static_cast<float32>(point.m_x);
    np.m_position[1] = static_cast<float32>(point.m_y);
    np.m_position[2] = static_cast<float32>(point.m_z);
    np.m_dist2[0] = static_cast<float32>(max_dist * max_dist);

    np.m_index[0] = 0;  // no photon found yet

    // Locate the nearest photon.
    locate_nearest_photon_with_irradiance(normal, &np, 1);

    const Photon *nearest = np.m_index[0];

    if(!nearest)
        return Color3(0.0);

    return shadingdata.m_reflectance * nearest->GetPrecomputedRadiance();
}

#endif  // ENABLE_RADIANCES_PRECOMPUTATION

void PhotonMap::BuildPhotonDirectionHistogram(const Point3 &position,
                                              const SurfaceBasis &surfacebasis,
                                              Real max_dist,
                                              int max_photons,
                                              Real *histogram,
                                              int m, int n,
                                              Real *total_power) const
{
    assert(max_photons + 1 <= BUF_SIZE);
    assert(histogram);
    assert(total_power);

    for(int i = m * n - 1; i >= 0; --i)
        histogram[i] = 0.0;

    *total_power = 0.0;

    nearest_photons np;
    np.m_dist2 = m_dist2_buf;
    np.m_index = m_index_buf;
    np.m_position[0] = static_cast<float32>(position.m_x);
    np.m_position[1] = static_cast<float32>(position.m_y);
    np.m_position[2] = static_cast<float32>(position.m_z);
    np.m_max = max_photons;
    np.m_found = 0;
    np.m_got_heap = 0;
    np.m_dist2[0] = static_cast<float32>(max_dist * max_dist);

    // Locate the nearest photons.
    locate_photons(&np, 1);

    for(int i = 1; i <= np.m_found; ++i) {
        const Photon *p = np.m_index[i];

#ifdef CHECK_FOR_THIN_SURFACES
        if(p->GetIncidentDirection() * surfacebasis.GetNormal() <= 0.0)
            continue;
#endif  // CHECK_FOR_THIN_SURFACES

#ifdef ENABLE_RADIANCES_PRECOMPUTATION
        // In some cases, this test significantly improves the quality of
        // the indirect lighting in corners and along edges.
        const Real THRESHOLD = 0.0; // 0.9?
        if(p->GetSurfaceNormal() * surfacebasis.GetNormal() < THRESHOLD)
            continue;
#endif  // ENABLE_RADIANCES_PRECOMPUTATION

        // Get photon incident direction back.
        const Vector3 photon_dir = p->GetIncidentDirection();

        // Compute photon direction in surface space.
        const Vector3 local_photon_dir = surfacebasis.TransformToLocal(photon_dir);
        assert(local_photon_dir.IsUnitLength());
        assert(local_photon_dir.m_y > 0.0);

        // Compute spherical coordinates of the photon direction.
        Real phi, theta;
        VectorToSphericalCoords(local_photon_dir, &phi, &theta);
        assert(phi >= 0.0 && phi <= 2.0 * PI);
        assert(theta >= 0.0 && theta <= PI);

        const Real cos_theta = cos(theta);
        const Real u1 = phi / (2.0 * PI);
        const Real v1 = cos_theta * cos_theta;

        assert(u1 >= 0.0 && u1 <= 1.0);
        assert(v1 >= 0.0 && v1 <= 1.0);

        const int hi = min(m - 1, static_cast<int>(u1 * m));
        const int hj = min(n - 1, static_cast<int>(v1 * n));

        assert(hi >= 0 && hi < m);
        assert(hj >= 0 && hj < n);

        const Real photon_power = /*1.0*/p->GetPower().Average();

        histogram[hj * m + hi] += photon_power;
        *total_power += photon_power;
    }

    const Real frac = *total_power / (m * n * 10.0);

    for(int i = m * n - 1; i >= 0; --i) {
        if(histogram[i] == 0.0) {
            histogram[i] = frac;
            *total_power += frac;
        }
    }
}

bool PhotonMap::WriteToFile(const string &filename,
                            ProgressMonitor *progmon /*= 0*/) const
{
    try {
        BinaryStream os(
            filename.c_str(),
            BinaryStream::OutputStream,
            BinaryStream::BigEndian // by convention, photon map files are written using big endian convention
        );

        if(!os.IsOpen())
            return false;

        // Set buffer size.
        os.SetBufferSize(PHOTONMAP_IO_BUFFER_SIZE);

        // Write photon map file signature.
        os.Write(m_photon_map_sig.c_str(), m_photon_map_sig.size());

        // Write the number of stored photons.
        os.Write(m_stored_photons);

        if(progmon && m_stored_photons > 0)
            progmon->StartJob(0.0, m_stored_photons);

        // Write each photon individually.
        for(int i = 1; i <= m_stored_photons; ++i) {
            m_photons[i].WriteToStream(os);

            if(progmon && ((i & 16383) == 0))
                progmon->SetJobProgress(i);
        }

        if(progmon)
            progmon->Done();

        // Write the bounding box of the photon array.
        os.Write(m_bbox_min[0]);
        os.Write(m_bbox_min[1]);
        os.Write(m_bbox_min[2]);
        os.Write(m_bbox_max[0]);
        os.Write(m_bbox_max[1]);
        os.Write(m_bbox_max[2]);

        // Write remaining fields.
        os.Write(m_half_stored_photons);

        // The stream is automatically closed.
        return true;
    }
    catch(BinaryStream::BaseException) {
        // I/O error.
        return false;
    }
}

PhotonMap *PhotonMap::CreateFromFile(const string &filename,
                                     ProgressMonitor *progmon /*= 0*/)
{
    PhotonMap *photon_map = 0;

    try {
        BinaryStream is(
            filename.c_str(),
            BinaryStream::InputStream,
            BinaryStream::BigEndian // by convention, photon map files are written using big endian convention
        );

        if(!is.IsOpen())
            return 0;

        // Set buffer size.
        is.SetBufferSize(PHOTONMAP_IO_BUFFER_SIZE);

        // Read photon map file signature and make sure it is valid.
        char *buf = new char[m_photon_map_sig.size() + 1];
        is.Read(buf, m_photon_map_sig.size());
        buf[m_photon_map_sig.size()] = '\0';
        if(string(buf) != m_photon_map_sig)
            return 0;
        delete [] buf;

        // Create the photon map.
        photon_map = new PhotonMap();

        // Read the number of stored photons.
        is.Read(&photon_map->m_stored_photons);

        if(photon_map->m_stored_photons > 0)
            photon_map->m_photons.reserve(photon_map->m_stored_photons);

        if(progmon && photon_map->m_stored_photons > 0)
            progmon->StartJob(1.0, photon_map->m_stored_photons);

        // Read each photon individually.
        for(int i = 1; i <= photon_map->m_stored_photons; ++i) {
            Photon photon;
            photon.ReadFromStream(is);
            photon_map->m_photons.push_back(photon);

            if(progmon && ((i & 32767) == 0))
                progmon->SetJobProgress(i);
        }

        if(progmon)
            progmon->Done();

        // Read the bounding box of the photon array.
        is.Read(&photon_map->m_bbox_min[0]);
        is.Read(&photon_map->m_bbox_min[1]);
        is.Read(&photon_map->m_bbox_min[2]);
        is.Read(&photon_map->m_bbox_max[0]);
        is.Read(&photon_map->m_bbox_max[1]);
        is.Read(&photon_map->m_bbox_max[2]);

        // Read remaining fields.
        is.Read(&photon_map->m_half_stored_photons);

        // The stream is automatically closed.
        return photon_map;
    }
    catch(BinaryStream::BaseException) {
        // I/O error.
        delete photon_map;
        return 0;
    }
}

void PhotonMap::median_split(Photon **p,
                             int start, int end,
                             int median,
                             int axis) const
{
    assert(p);

    int left = start;
    int right = end;

    while(right > left) {
        const float32 v = p[right]->m_position[axis];
        int i = left - 1;
        int j = right;

        while(true) {
            while(p[++i]->m_position[axis] < v) ;
            while(p[--j]->m_position[axis] > v && j > left) ;

            if(i >= j)
                break;

            swap(p[i], p[j]);
        }

        swap(p[i], p[right]);

        if(i >= median)
            right = i - 1;
        if(i <= median)
            left = i + 1;
    }
}

void PhotonMap::balance_segment(Photon **pbal,
                                Photon **porg,
                                int index,
                                int start, int end)
{
    assert(pbal);
    assert(porg);

    // Compute new median.

    int median = 1;

    while(4 * median <= end - start + 1)
        median += median;

    if(3 * median <= end - start + 1) {
        median += median;
        median += start - 1;
    } else median = end - median + 1;

    // Find axis to split along.

    int axis = 2;

    if( (m_bbox_max[0] - m_bbox_min[0] > m_bbox_max[1] - m_bbox_min[1]) &&
        (m_bbox_max[0] - m_bbox_min[0] > m_bbox_max[2] - m_bbox_min[2]))
        axis = 0;
    else if(m_bbox_max[1] - m_bbox_min[1] > m_bbox_max[2] - m_bbox_min[2])
        axis = 1;

    // Partition photon block around the median.

    median_split(porg, start, end, median, axis);

    pbal[index] = porg[median];
    pbal[index]->SetSplittingPlaneAxis(axis);

    // Recursively balance the left and right block.

    if(median > start) {
        // Balance left segment.
        if(start < median - 1) {
            const float32 tmp = m_bbox_max[axis];
            m_bbox_max[axis] = pbal[index]->m_position[axis];
            balance_segment(pbal, porg, 2 * index, start, median - 1);
            m_bbox_max[axis] = tmp;
        } else pbal[2 * index] = porg[start];
    }

    if(median < end) {
        // Balance right segment.
        if(median + 1 < end) {
            const float32 tmp = m_bbox_min[axis];
            m_bbox_min[axis] = pbal[index]->m_position[axis];
            balance_segment(pbal, porg, 2 * index + 1, median + 1, end);
            m_bbox_min[axis] = tmp;
        } else pbal[2 * index + 1] = porg[end];
    }
}

void PhotonMap::locate_photons(nearest_photons *np, int index) const {
    const Photon *p = &m_photons[index];

    if(index < m_half_stored_photons) {
        // We're not in a leaf.

        const int plane = p->GetSplittingPlaneAxis();
        const float32 dist1 = np->m_position[plane] - p->m_position[plane];

        if(dist1 > 0.0f) {  // if dist1 is positive, search right plane first
            locate_photons(np, 2 * index + 1);
            // Also search left plane if updated bounding sphere crosses the splitting plane.
            if(dist1 * dist1 < np->m_dist2[0])
                locate_photons(np, 2 * index);
        } else {    // dist1 is negative, search left plane first
            locate_photons(np, 2 * index);
            // Also search right plane if updated bounding sphere crosses the splitting plane.
            if(dist1 * dist1 < np->m_dist2[0])
                locate_photons(np, 2 * index + 1);
        }
    }

    // Compute squared distance between current photon and np->m_position.
    float32 dist1;
    dist1 = p->m_position[0] - np->m_position[0];
    float32 dist2 = dist1 * dist1;
    dist1 = p->m_position[1] - np->m_position[1];
    dist2 += dist1 * dist1;
    dist1 = p->m_position[2] - np->m_position[2];
    dist2 += dist1 * dist1;

    if(dist2 < np->m_dist2[0]) {    // np->m_dist2[0] is the farthest photon
        // We found a photon. Insert it in the candidate list.

        if(np->m_found < np->m_max) {
            // Heap is not full; use array.
            ++np->m_found;
            np->m_dist2[np->m_found] = dist2;
            np->m_index[np->m_found] = p;
        } else {
            int j, parent;

            if(np->m_got_heap == 0) {   // do we need to build the heap?
                // Build heap.
                
                float32 dst2;
                const Photon *phot;
                const int half_found = np->m_found >> 1;

                for(int k = half_found; k >= 1; --k) {
                    parent = k;
                    phot = np->m_index[k];
                    dst2 = np->m_dist2[k];

                    while(parent <= half_found) {
                        j = parent + parent;

                        if(j < np->m_found && np->m_dist2[j] < np->m_dist2[j + 1])
                            ++j;

                        if(dst2 >= np->m_dist2[j])
                            break;

                        np->m_dist2[parent] = np->m_dist2[j];
                        np->m_index[parent] = np->m_index[j];

                        parent = j;
                    }

                    np->m_dist2[parent] = dst2;
                    np->m_index[parent] = phot;
                }

                np->m_got_heap = 1;
            }

            // Insert new photon into max heap. Delete largest element,
            // insert new, and reorder the heap.

            parent = 1;
            j = 2;  // left child of the root

            while(j <= np->m_found) {
                if(j < np->m_found && np->m_dist2[j] < np->m_dist2[j + 1])
                    ++j;

                if(dist2 > np->m_dist2[j])
                    break;

                np->m_dist2[parent] = np->m_dist2[j];
                np->m_index[parent] = np->m_index[j];
                parent = j;

                j += j;
            }

            //!\todo The 'pane' global illumination renderer (http://www.csit.fsu.edu/~beason/pane)
            //! does a check here. Check whether or not this is relevant.
            np->m_index[parent] = p;
            np->m_dist2[parent] = dist2;

            np->m_dist2[0] = np->m_dist2[1];
        }
    }
}

#ifdef ENABLE_RADIANCES_PRECOMPUTATION

void PhotonMap::locate_nearest_photon_with_irradiance(const Vector3 &normal,
                                                      nearest_photons *np,
                                                      int index) const
{
    const Photon *p = &m_photons[index];

    if(index < m_half_stored_photons) {
        // We're not in a leaf.

        const int plane = p->GetSplittingPlaneAxis();
        const float32 dist1 = np->m_position[plane] - p->m_position[plane];

        if(dist1 > 0.0f) {  // if dist1 is positive, search right plane first
            locate_nearest_photon_with_irradiance(normal, np, 2 * index + 1);
            // Also search left plane if updated bounding sphere crosses the splitting plane.
            if(dist1 * dist1 < np->m_dist2[0])
                locate_nearest_photon_with_irradiance(normal, np, 2 * index);
        } else {    // dist1 is negative, search left plane first
            locate_nearest_photon_with_irradiance(normal, np, 2 * index);
            // Also search right plane if updated bounding sphere crosses the splitting plane.
            if(dist1 * dist1 < np->m_dist2[0])
                locate_nearest_photon_with_irradiance(normal, np, 2 * index + 1);
        }
    }

    if(p->HasPrecomputedRadiance()) {
        const Real THRESHOLD = 0.9;

        if(normal * p->GetSurfaceNormal() >= THRESHOLD) {
            // Compute squared distance between current photon and np->m_position.
            float32 dist1;
            dist1 = p->m_position[0] - np->m_position[0];
            float32 dist2 = dist1 * dist1;
            dist1 = p->m_position[1] - np->m_position[1];
            dist2 += dist1 * dist1;
            dist1 = p->m_position[2] - np->m_position[2];
            dist2 += dist1 * dist1;

            if(dist2 < np->m_dist2[0]) {
                // We found a photon closer than the closest one until now.
                np->m_dist2[0] = dist2;
                np->m_index[0] = p;
            }
        }
    }
}

#endif  // ENABLE_RADIANCES_PRECOMPUTATION

---