sphere.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 "sphere.h" // include first
#include "common/math/real.h"
#include "common/math/sampling.h"
#include "context.h"
#include "hit.h"
#include "iedf.h"
#include "isurfaceshader.h"
#include "ray.h"
#include "scene.h"
#include "settings.h"
#include "statistics.h"
#include "surfacebasis.h"
#include "utilities.h"

#include <cassert>
#include <cmath>

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

Sphere::Sphere(const Matrix4 &m,
               const ISurfaceShader *surface_shader,
               IntersectionMask intersection_mask /*= INTERSECT_ALL_RAYS*/) :
    IAreaLight(m, surface_shader, intersection_mask),
    m_center(TransformToWorld(Point3(0.0))),
    m_ex(TransformToWorld(Vector3(1.0, 0.0, 0.0))),
    m_ey(TransformToWorld(Vector3(0.0, 1.0, 0.0))),
    m_ez(TransformToWorld(Vector3(0.0, 0.0, 1.0))),
    m_bs(m_center)
{
    m_aabb = AABB3(Point3(-1.0), Point3(1.0)).Transform(m);

    // Enlarge the object bounding box to avoid numerical instabilities.
    m_aabb.Extend(OBJECT_AABB_EXTENSION);

    m_bs.Include(TransformToWorld(Point3(1.0, 0.0, 0.0)));
    m_bs.Include(TransformToWorld(Point3(-1.0, 0.0, 0.0)));
    m_bs.Include(TransformToWorld(Point3(0.0, 1.0, 0.0)));
    m_bs.Include(TransformToWorld(Point3(0.0, -1.0, 0.0)));
    m_bs.Include(TransformToWorld(Point3(0.0, 0.0, 1.0)));
    m_bs.Include(TransformToWorld(Point3(0.0, 0.0, -1.0)));

    // Compute the square of the radius of the bounding sphere
    // (required by the ComputeIrradiance() method).
    m_bs.ComputeSquareRadius();
}

Real Sphere::ComputeSurfaceArea() const {
    const Real r = m_ex.Norm();

    if(!(r == m_ey.Norm() && r == m_ez.Norm())) {
        assert(!"Ellipsoid surface area computation not implemented yet.");
        //!\todo Implement.
        //! Ellipsoid Surface Area: http://documents.wolfram.com/v4/MainBook/G.1.7.html
        //! Source Code: http://www.tcs.auckland.ac.nz/dgt/Source_Code.php?id=2
    }

    return 4.0 * PI * r * r;
}

void Sphere::EvaluateSurface(const Point2 &input,
                             Point3 *point,
                             Vector3 *geometric_normal) const
{
    assert(point);
    assert(geometric_normal);

    const Vector3 direction = UniformSphereSampling(input.m_x, input.m_y);

    *point =
        m_center +
        m_ex * direction.m_x +
        m_ey * direction.m_y +
        m_ez * direction.m_z;

    *geometric_normal = TransformNormalToWorld(direction);
}

bool Sphere::Intersect(const Context &context,
                       const Ray &ray,  //!< 'ray' is expressed in world space.
                       Hit *hit /*= 0*/) const
{
    // Test whether this object intersects with this type of ray.
    if(!(m_intersection_mask & ray.GetType()))
        return false;

    ++context.m_statistics->m_tested_intersections;

    const Ray r = TransformToLocal(ray);
    const Real inv_scale2 = 1.0 / r.m_direction.SquareNorm();

    const Real a = r.m_origin * r.m_direction;
    const Real b = r.m_origin * r.m_origin - 1.0;

    const Real delta = a * a * inv_scale2 - b;

    if(delta < 0.0)
        return false;   // the ray does not intersect the sphere

    const Real inv_scale = sqrt(inv_scale2);
    const Real c = a * inv_scale;

    const Real sqrt_delta = sqrt(delta);

    Real t = -c - sqrt_delta;

    if(t < 0.0) {
        t = -c + sqrt_delta;

        if(t < 0.0)
            return false;   // the intersection point is not on the ray
    }

    // The ray does intersect the sphere.
    ++context.m_statistics->m_intersections_found;

    if(hit) {
        t *= inv_scale;

        hit->m_object = this;
        hit->m_abscissa = t;    // world space
        //!\todo Not sure: should t in r.GetPointAt(t) be in object space or in world space?
        hit->m_geometric_normal = r.GetPointAt(t);          // object space
        hit->m_shading_normal = hit->m_geometric_normal;    // object space

        // Normals will be made unit-length later, when absolutely necessary.
    }

    return true;
}

void Sphere::ComputeIrradiance(const Context &context,
                               const Scene *scene,
                               const Point3 &point,
                               const Vector3 &geometric_normal,
                               const Vector3 &shading_normal,
                               const ISurfaceSampler::SampleVector &input,
                               IrradianceSampleVector *output) const
{
    assert(scene);
    assert(geometric_normal.IsUnitLength());
    assert(shading_normal.IsUnitLength());
    assert(output);

    const Point3 shifted_point = point + 1.0e-8 * geometric_normal;

    const Vector3 psi(m_bs.m_center - point);
    const Real psi_inv_norm2 = 1.0 / psi.SquareNorm();
    const Real cos_theta_max = sqrt(1.0 - m_bs.m_radius2 * psi_inv_norm2);
    const Real k = 1.0 - cos_theta_max;
    const Vector3 unit_psi(psi * sqrt(psi_inv_norm2));
    const Real sphere_solid_angle = 2.0 * PI * k;

    // Build an orthonormal basis around the Psi vector.
    const SurfaceBasis surfacebasis(unit_psi);

    for(ISurfaceSampler::SampleVector::const_iterator i = input.begin(), e = input.end(); i != e; ++i) {
        //!\todo Optimize (inline and simplify expressions).
        const Real phi = 2.0 * PI * i->m_x;
        const Real theta = acos(1.0 - i->m_y * k);
        const Vector3 direction = SphericalCoordsToVector(phi, theta);

        const Vector3 incoming = surfacebasis.TransformToWorld(direction);
        assert(incoming.IsUnitLength());

        Hit hit;

        // Make sure that this sample really is on the emitter. It is always the
        // case for perfectly spherical light sources, but it might not be the case
        // for ellipsoid light sources, because we sample the cone sustained by
        // the bounding sphere of the light source.
        if(!Intersect(context, Ray(Ray::VISIBILITY_RAY, point, incoming), &hit))
            continue;   // this sample is not on the emitter: next sample

        const Real distance_to_emitter = hit.m_abscissa;

        if(m_cast_shadows) {
            ++context.m_statistics->m_shadow_rays;

            const Ray ray(Ray::SHADOW_RAY, shifted_point, incoming);

            if(scene->Trace(context, ray, &hit)) {
                const bool light_source_visible =
                    hit.m_object == this ||
                    hit.m_abscissa > distance_to_emitter;

                if(!light_source_visible) {
                    output->push_back(IrradianceSample(Color3(0.0), incoming));
                    continue;   // this light is occluded from this direction: next sample
                }
            }
        }

        //!\todo Flip the normal? Call Hit::ExtractIntersection()?
        const Vector3 sample_geometric_normal =
            hit.m_object->TransformNormalToWorld(hit.m_geometric_normal);
        assert(sample_geometric_normal.IsUnitLength());

        //!\todo Not sure which normal to use.
        const Real cos_theta = shading_normal * incoming;

        if(cos_theta <= 0.0) {
            // The surface is not directed toward the light: next sample.
            output->push_back(IrradianceSample(Color3(0.0), incoming));
            continue;
        }

        const Real emitted_radiance = GetSurfaceShader()->GetEDF()->GetEmittedRadiance(
            context,
            sample_geometric_normal,
            -incoming);

        assert(emitted_radiance >= 0.0);

        if(emitted_radiance <= 0.0) {
            // The light surface element is not directed toward the point: next sample.
            output->push_back(IrradianceSample(Color3(0.0), incoming));
            continue;
        }

        const Color3 irradiance = (sphere_solid_angle * cos_theta * emitted_radiance) *
            m_surface_shader->GetRadiantExitance();

        output->push_back(IrradianceSample(irradiance, incoming));
    }
}

---