world.C

/*
 * This file is part of the "Archon" framework.
 * (http://files3d.sourceforge.net)
 *
 * Copyright © 2002 by Kristian Spangsege and Brian Kristiansen.
 *
 * Permission to use, copy, modify, and distribute this software and
 * its documentation under the terms of the GNU General Public License is
 * hereby granted. No representations are made about the suitability of
 * this software for any purpose. It is provided "as is" without express
 * or implied warranty. See the GNU General Public License
 * (http://www.gnu.org/copyleft/gpl.html) for more details.
 *
 * The characters in this file are ISO8859-1 encoded.
 *
 * The documentation in this file is in "Doxygen" style
 * (http://www.doxygen.org).
 */

#include <float.h>
#include <vector>

#include <archon/util/random.H>
#include <archon/util/color.H>
#include <archon/util/image.H>

#include "light.H"
#include "object.H"
#include "world.H"

#if ! defined M_PI
#define M_E           2.7182818284590452353602874713526625  /* e */
#define M_LOG2E       1.4426950408889634073599246810018922  /* log_2 e */
#define M_LOG10E      0.4342944819032518276511289189166051  /* log_10 e */
#define M_LN2         0.6931471805599453094172321214581766  /* log_e 2 */
#define M_LN10        2.3025850929940456840179914546843642  /* log_e 10 */
#define M_PI          3.1415926535897932384626433832795029  /* pi */
#define M_PI_2        1.5707963267948966192313216916397514  /* pi/2 */
#define M_PI_4        0.7853981633974483096156608458198757  /* pi/4 */
#define M_1_PI        0.3183098861837906715377675267450287  /* 1/pi */
#define M_2_PI        0.6366197723675813430755350534900574  /* 2/pi */
#define M_2_SQRTPI    1.1283791670955125738961589031215452  /* 2/sqrt(pi) */
#define M_SQRT2       1.4142135623730950488016887242096981  /* sqrt(2) */
#define M_SQRT1_2     0.7071067811865475244008443621048490  /* 1/sqrt(2) */
#endif

namespace Archon
{
  using namespace Utilities;
  using namespace Math;

  namespace Raytracer
  {
    World::World(const Vector4 &backgroundColor,
                 double significanceThreshold,
                 int maxDepth,
                 bool enableTransmission,
                 unsigned photonsInEstimate,
                 bool showOnlyPhotonMap):
      globalAmbientIntencity(0.4),
      backgroundColor(backgroundColor),
      significanceThreshold(significanceThreshold),
      maxDepth(maxDepth),
      enableTransmission(enableTransmission),
      photonsInEstimate(photonsInEstimate),
      showOnlyPhotonMap(showOnlyPhotonMap) {}

    World::World(double globalAmbientIntencity,
                 const Vector4 &backgroundColor,
                 double significanceThreshold,
                 int maxDepth,
                 bool enableTransmission,
                 unsigned photonsInEstimate,
                 bool showOnlyPhotonMap):
      globalAmbientIntencity(globalAmbientIntencity),
      backgroundColor(backgroundColor),
      significanceThreshold(significanceThreshold),
      maxDepth(maxDepth),
      enableTransmission(enableTransmission),
      photonsInEstimate(photonsInEstimate),
      showOnlyPhotonMap(showOnlyPhotonMap)
    {
      textureList.push_back(Image("texture/leafpolar.ppm"));     // 0
      textureList.push_back(Image("texture/moon1polar.ppm"));    // 1
      textureList.push_back(Image("texture/kristian.ppm"));      // 2
      textureList.push_back(Image("texture/granite.ppm"));       // 3
      textureList.push_back(Image("texture/skin.ppm"));          // 4
      textureList.push_back(Image("texture/spotty.ppm"));        // 5
      textureList.push_back(Image("texture/bkristia.ppm"));      // 6
      textureList.push_back(Image("texture/chess-marble.ppm"));  // 7
      textureList.push_back(Image("texture/glass.ppm"));         // 8

      surfaceList.push_back(new StandardSurface(Color::seagreen,  0.05, 0.1, 0.05, 0.9, 20, 20, 1.15));    // 0
      surfaceList.push_back(new StandardSurface(Color::crimson,   0.1, 0.2, 0.9, 0.1, 60, 20, 1.2));    // 1
      surfaceList.push_back(new StandardSurface(Color::steelblue, 0.1, 0.2, 0.9, 0.1, 60, 20, 1.2));    // 2
      surfaceList.push_back(new StandardSurface(Color::darkolivegreen));                                // 3
      surfaceList.push_back(new StandardSurface(Color::gold));                                          // 4
      surfaceList.push_back(new Texture(textureList[3], 0.1, 0.8, 0.0, 0.0, 20, 20));                   // 5
      surfaceList.push_back(new Texture(textureList[4], 0.9, 0.9, 0.1, 0.1, 20, 20, 1, 3, 3));          // 6
      surfaceList.push_back(new Texture(textureList[5], 0.1, 0.5, 0.5, 0.1, 20, 20, 1, 2, 2));          // 7
      surfaceList.push_back(new Texture(textureList[2], 0.1, 0.9, 0.1, 0.1, 20, 20));                   // 8
      surfaceList.push_back(new Texture(textureList[6], 0.1, 0.2, 0.8, 0.8, 20, 20));                   // 9
      surfaceList.push_back(new Texture(textureList[6], 0.1, 0.9, 0.1, 0.1, 20, 20));                   // 10
      surfaceList.push_back(new StandardSurface(Vector4(1, 1, 1), 0.0, 0.1, 0.2, 0.6, 20, 20, 1.33));  // 11
      surfaceList.push_back(new Texture(textureList[8], 0.9, 0, 0, 0, 20, 20));                         // 12
      surfaceList.push_back(new Texture(textureList[7], 0.1, 0.0, 0.5, 0.5, 20, 20, 1, 3, 3));          // 13
      surfaceList.push_back(new Texture(textureList[8], 0.1, 0.18, 0.05, 0.9, 60, 20, 1.2));            // 14
      surfaceList.push_back(new StandardSurface(Vector4(1, 1, 1), 0, 0.8, 0.1, 0, 20, 20, 1.33));  // 15


      //    objectList.push_back(new Sphere(surfaceList[7], Vector3(-1.5,  1.5, -3), Vector3::unitY, -Vector3::unitZ, 1));
      //    objectList.push_back(new Sphere(surfaceList[6], Vector3( 0,    0.5, -2), Vector3::unitY, -Vector3::unitZ, 1));
      //    objectList.push_back(new Sphere(surfaceList[5], Vector3( 1.5, -0.5, -1), Vector3::unitY, -Vector3::unitZ, 1));
      //    objectList.push_back(new Sphere(surfaceList[0], Vector3(-3, 1.6, 1.7), Vector3::unitY, -Vector3::unitZ, 1.2));
      //    objectList.push_back(new Sphere(surfaceList[2], Vector3(0.5, 0, 0), Vector3::unitY, -Vector3::unitZ, 1));
      //objectList.push_back(new Sphere(surfaceList[11], Vector3(0, 4, 0), Vector3::unitY, -Vector3::unitZ, 2));

      //    objectList.push_back(new Torus(surfaceList[1], Vector3(-1, 3, -2), Vector3::unitY, Vector3::unitX, 2, 0.8));
      objectList.push_back(new Torus(surfaceList[11], CoordSystem3(Vector3(0, 4, 0), Matrix3x3::identity), 2, 1));

      // box
      /*
        {
        Vector3 v0(-1, -1, -1);
        Vector3 v1(-1, -1,  1);
        Vector3 v2( 1, -1,  1);
        Vector3 v3( 1, -1, -1);
        Vector3 v4(-1,  1, -1);
        Vector3 v5(-1,  1,  1);
        Vector3 v6( 1,  1,  1);
        Vector3 v7( 1,  1, -1);

        CoordSystem3 modelSystem(Vector3::zero, Matrix3x3::identity);
        modelSystem.origin += Vector3(0.7, 2, 3.4);
        //     modelSystem.origin += normalize(Vector3(0, 1, 5) - modelSystem.origin)*1;
        modelSystem.basis.rotateRect(Vector3(0, 1, 0), M_PI/4);
        modelSystem.basis.rotateRect(Vector3(1, 0, 0), M_PI/4);
        //modelSystem.basis *= 0.5;

        vertexList.push_back(new Vector3(modelSystem * v0));
        vertexList.push_back(new Vector3(modelSystem * v1));
        vertexList.push_back(new Vector3(modelSystem * v2));
        vertexList.push_back(new Vector3(modelSystem * v3));
        vertexList.push_back(new Vector3(modelSystem * v4));
        vertexList.push_back(new Vector3(modelSystem * v5));
        vertexList.push_back(new Vector3(modelSystem * v6));
        vertexList.push_back(new Vector3(modelSystem * v7));

        std::vector<Vector3 *> left;
        left.push_back(vertexList[0]);
        left.push_back(vertexList[1]);
        left.push_back(vertexList[5]);
        left.push_back(vertexList[4]);

        std::vector<Vector3 *> right;
        right.push_back(vertexList[2]);
        right.push_back(vertexList[3]);
        right.push_back(vertexList[7]);
        right.push_back(vertexList[6]);

        std::vector<Vector3 *> down;
        down.push_back(vertexList[3]);
        down.push_back(vertexList[2]);
        down.push_back(vertexList[1]);
        down.push_back(vertexList[0]);

        std::vector<Vector3 *> up;
        up.push_back(vertexList[4]);
        up.push_back(vertexList[5]);
        up.push_back(vertexList[6]);
        up.push_back(vertexList[7]);

        std::vector<Vector3 *> far;
        far.push_back(vertexList[0]);
        far.push_back(vertexList[4]);
        far.push_back(vertexList[7]);
        far.push_back(vertexList[3]);

        std::vector<Vector3 *> near;
        near.push_back(vertexList[1]);
        near.push_back(vertexList[2]);
        near.push_back(vertexList[6]);
        near.push_back(vertexList[5]);


        std::vector<Polygon *> box1;
        box1.push_back(new Polygon(surfaceList[14],  left,  Vector2::zero, Vector2::unitX));
        box1.push_back(new Polygon(surfaceList[14],  right, Vector2::zero, Vector2::unitX));
        box1.push_back(new Polygon(surfaceList[14],  down,  Vector2::zero, Vector2::unitX));
        box1.push_back(new Polygon(surfaceList[14],  up,    Vector2::zero, Vector2::unitX));
        box1.push_back(new Polygon(surfaceList[14],  far,   Vector2::zero, Vector2::unitX));
        box1.push_back(new Polygon(surfaceList[14],  near,  Vector2::zero, Vector2::unitX));

        objectList.push_back(new ConvexPolyhedron(box1));
        }
      */
      // floor
      {
        Vector3 v0(-1, 0, -1);
        Vector3 v1(-1, 0,  1);
        Vector3 v2( 1, 0,  1);
        Vector3 v3( 1, 0, -1);

        CoordSystem3 modelSystem(Vector3(0, 0, 0), Matrix3x3::identity);
        modelSystem.origin += Vector3(0, 0, 0);
        modelSystem.basis *= 5;

        const int n = vertexList.size();

        vertexList.push_back(new Vector3(modelSystem(v1)));
        vertexList.push_back(new Vector3(modelSystem(v2)));
        vertexList.push_back(new Vector3(modelSystem(v3)));
        vertexList.push_back(new Vector3(modelSystem(v0)));

        std::vector<Vector3 *> floor;
        floor.push_back(vertexList[n + 0]);
        floor.push_back(vertexList[n + 1]);
        floor.push_back(vertexList[n + 2]);
        floor.push_back(vertexList[n + 3]);

        objectList.push_back(new Polygon(surfaceList[15], floor, Vector2(0, 0), Vector2(1, 0)));
      }

      // right wall
      {
        Vector3 v0(-1, 0, -1);
        Vector3 v1(-1, 0,  1);
        Vector3 v2( 1, 0,  1);
        Vector3 v3( 1, 0, -1);

        CoordSystem3 modelSystem(Vector3(0, 0, 0), Matrix3x3::identity);
        modelSystem.origin += Vector3(5, 5, 0);
        modelSystem.basis.roll(M_PI/2);
        modelSystem.basis *= 5;

        const int n = vertexList.size();

        vertexList.push_back(new Vector3(modelSystem(v1)));
        vertexList.push_back(new Vector3(modelSystem(v2)));
        vertexList.push_back(new Vector3(modelSystem(v3)));
        vertexList.push_back(new Vector3(modelSystem(v0)));

        std::vector<Vector3 *> polygon;
        polygon.push_back(vertexList[n + 0]);
        polygon.push_back(vertexList[n + 1]);
        polygon.push_back(vertexList[n + 2]);
        polygon.push_back(vertexList[n + 3]);

        objectList.push_back(new Polygon(surfaceList[15], polygon, Vector2(0, 0), Vector2(1, 0)));
      }


      // background
      /*
        {
        Vector3 v0(-1, 0, -1);
        Vector3 v1(-1, 0,  1);
        Vector3 v2( 1, 0,  1);
        Vector3 v3( 1, 0, -1);

        CoordSystem3 modelSystem(Vector3(0, 1, 5), Matrix3x3::identity);
        modelSystem.basis.pitchRect(M_PI/2-M_PI/10);
        modelSystem.origin += modelSystem.basis.z * -10;
        modelSystem.origin += modelSystem.basis.y * -40;
        modelSystem.basis *= 25;

        const int n = vertexList.size();

        vertexList.push_back(new Vector3(modelSystem.transform(v1)));
        vertexList.push_back(new Vector3(modelSystem.transform(v2)));
        vertexList.push_back(new Vector3(modelSystem.transform(v3)));
        vertexList.push_back(new Vector3(modelSystem.transform(v0)));

        std::vector<Vector3 *> background;
        background.push_back(vertexList[n + 0]);
        background.push_back(vertexList[n + 1]);
        background.push_back(vertexList[n + 2]);
        background.push_back(vertexList[n + 3]);

        objectList.push_back(new Polygon(surfaceList[12], background, Vector2::zero, Vector2::unitX));
        }
      */

      //    lightList.push_back(new Light(Vector3(10, 0, 0), Color::darkmagenta));
      //    lightList.push_back(new Light(Vector3(10, 10, 0), Color::darkred * 1.2));
      //    lightList.push_back(new Light(Vector3(0, 0, 10), Color::goldenrod * 1.2));

      // lightList.push_back(new Light(Vector3(-7, 10, 5), Color::white));
      // lightList.push_back(new Light(Vector3( 0,  5, 5), Color::white));
      // lightList.push_back(new Light(Vector3( 7, 10, 5), Color::white));
      lightList.push_back(new Light(Vector3(0, 8, -2), Color::red));
      lightList.push_back(new Light(Vector3(-2, 8, 0), Color::lime));
      lightList.push_back(new Light(Vector3(0, 8, 2), Color::blue));
    }



    World::~World()
    {
      for(std::vector<Light   *>::iterator i =   lightList.begin();
          i !=   lightList.end(); ++i) delete *i;
      for(std::vector<Object  *>::iterator i =  objectList.begin();
          i !=  objectList.end(); ++i) delete *i;
      for(std::vector<Vector3 *>::iterator i =  vertexList.begin();
          i !=  vertexList.end(); ++i) delete *i;
      for(std::vector<Surface *>::iterator i = surfaceList.begin();
          i != surfaceList.end(); ++i) delete *i;
    }


    Vector4 World::trace(const Line3 &ray,
                       const Object *originObject,
                       const Object::Part *originPart,
                       int depth, double significance,
                       bool transmission)
    {
      /*
       * Determine the object which is hit first by 'ray' taking the
       * origin point and direction of 'ray' into account. At the same
       * time, determine the "part" of that object which was hit, and
       * the distance from the ray origin to the point of intersection.
       *
       * Note: We pass to the intersection methods a reference to the
       * object and object-part wherefrom the ray originates. This
       * enables the individual type of object to discriminate reliably
       * between the originating intersection point and the intersection
       * point of interrest. It also enables some types of objects to
       * reduce the cost of the intersection operation. Especially
       * convex objects need not check for external intersection against
       * themselves, like concave objects do.
       */
      const Object *hitObject = 0;
      const Object::Part *hitPart = 0;
      double hitDist = 0;
      if(transmission)
      {
        const Solid *s = dynamic_cast<const Solid *>(originObject);
        if(!s) throw InternalErrorException("Transmission initiated "
                                            "on non-solid object");
        hitObject = originObject;
        s->interiorIntersect(ray, originPart, hitDist, hitPart);
      }
      else
      {
        // Search for the nearest intersection
        const Object::Part *p;
        double d;
        for(std::vector<Object *>::iterator i=objectList.begin();
            i != objectList.end(); ++i)
        {
          if((*i)->intersect(ray, originObject, originPart, d, &p) &&
             (!hitObject || d < hitDist))
          {
            hitDist = d;
            hitPart = p;
            hitObject = *i;
          }
        }
      }

      /*
       * If the ray hits nothing, we think of the light as coming from
       * the "background" and therefore default to a fixed background
       * color.
       */
      if(!hitObject) return backgroundColor;


      // Determine the intersection point
      Vector3 intersection = ray.direction;
      intersection *= hitDist;
      intersection += ray.point;

      // Fetch surface, normal and possibly texture coordinates from the intersected object-part
      const Surface *surface = hitPart->getSurface();
      Vector3 normal;
      Vector2 texturePoint;
      if(dynamic_cast<const Texture *>(surface))
        hitPart->map(intersection, normal, texturePoint);
      else
        hitPart->map(intersection, normal);
      Vector4 surfaceColor;
      surface->getColor(texturePoint, surfaceColor);

      std::vector<const KDTree::Element *> nearest;
      double radius;
      causticPhotonMap.findNearest(intersection, photonsInEstimate, 30, radius, nearest);

      radius/=60;

      Vector4 color = Color::black;

      if(nearest.size() == photonsInEstimate)
      {
        Vector3 caustic(0, 0, 0);
        for(unsigned i=0; i<photonsInEstimate; ++i) caustic += ((Photon *)nearest[i])->power;
        caustic *= surface->getDiffuseReflection() / (M_PI * M_PI * radius * radius);
        color[0] = caustic[0] * surfaceColor[0];
        color[1] = caustic[1] * surfaceColor[1];
        color[2] = caustic[2] * surfaceColor[2];
      }

      if(showOnlyPhotonMap) return nearest.size() < photonsInEstimate ? Color::black : Color::white * (1 - (atan(radius*10) / (M_PI/2)));

      // Add up color contributions from visible light sources.
      for(std::vector<Light *>::iterator i=lightList.begin();
          i != lightList.end(); ++i)
      {
        /*
         * When testing for eclipsed light sources we actually generate
         * a ray from the intersection point towards the light source,
         * and test for intersection between that line and all the
         * objects in the world. We can terminate the check imediately,
         * though, as soon as we encounter an intersection, because this
         * guarantees eclipse.
         */
        Vector3 lightDirection = (*i)->getPosition();
        lightDirection -= intersection;
        const double lightDist = lightDirection.length();
        lightDirection.normalize();
        const Line3 rayToLight(intersection, lightDirection);

        /*
         * Since in practice it always is so, that the intersected object
         * shadows at least half of the world space from the point of
         * view of the intersection point, it generally pays to check
         * for eclipse by this object before checking other objects.
         *
         * Note that this check will be particularly simple for convex
         * objects since in this case, if the angle between the surface
         * normal and the light direction is greater than 90 degrees,
         * then the light source is eclipsed.
         */
        double d;
        if(dot(normal, lightDirection) < 0) continue;

        // Check for eclipse by other objects.
        bool eclipsed = false;
        for(std::vector<Object *>::iterator j=objectList.begin();
            j != objectList.end(); ++j)
        {
          if((*j)->intersect(rayToLight, hitObject, hitPart, d) &&
             d < lightDist)
          {
            eclipsed = true;
            break;
          }
        }
        if(eclipsed) continue;

        // Add color contribution to pixel
        Vector4 c = surface->shade(**i, normal,
                                 transmission ? ray.direction : -ray.direction,
                                 lightDirection, surfaceColor);
        color[0] += c[0];
        color[1] += c[1];
        color[2] += c[2];
      }

      // Add diffuse reflection of ambient light
      const double ambientReflection =
        surface->getAmbientReflection() * globalAmbientIntencity;
      color[0] += surfaceColor[0] * ambientReflection;
      color[1] += surfaceColor[1] * ambientReflection;
      color[2] += surfaceColor[2] * ambientReflection;

      /*
       * Follow specular reflection
       */
      const double specularReflection =
        surface->getSpecularReflection();
      const double reflectionSignificance =
        significance * specularReflection;
      if(depth < maxDepth && reflectionSignificance > significanceThreshold)
      {
        Vector3 reflectedDirection = normal;
        reflectedDirection *= -2 * dot(ray.direction, normal);
        reflectedDirection += ray.direction;
        reflectedDirection.normalize(); // In theory not needed, but for the sake of stability

        const Vector4 c =
          trace(Line3(intersection, reflectedDirection), hitObject, hitPart,
                depth + 1, reflectionSignificance, transmission);

        color[0] += c[0] * surfaceColor[0] * specularReflection;
        color[1] += c[1] * surfaceColor[1] * specularReflection;
        color[2] += c[2] * surfaceColor[2] * specularReflection;
      }

      /*
       * Note: We currently cannot allow thin objects to be transparent
       * due to lack of a consistent way to deal with such.
       */
      if(enableTransmission && dynamic_cast<const Solid *>(hitObject))
      {
        /*
         * Follow specular refraction into the surface
         * The determination of the direction vector of the refraction
         * is based on: http://www.acm.org/crossroads/xrds3-4/raytracing.html
         */
        const double specularRefraction =
          surface->getSpecularRefraction();
        const double refractionSignificance =
          significance * specularRefraction;
        if(depth < maxDepth && refractionSignificance > significanceThreshold)
        {
          const double cosine = dot(ray.direction, transmission ? -normal : normal);
          const double n = transmission ? 1 / surface->getRefractiveIndex() : surface->getRefractiveIndex();
          const double a = 1 - (1 - cosine*cosine) / (n*n);

          if(a > 0)
          {
            Vector3 v = transmission ? -normal : normal;
            v *= (cosine/n + sqrt(a));
            Vector3 refractedDirection = ray.direction;
            refractedDirection /= n;
            refractedDirection -= v;
            refractedDirection.normalize();

            const Vector4 c =
              trace(Line3(intersection, refractedDirection),
                    hitObject, hitPart,
                    depth + 1, refractionSignificance, !transmission);

            color[0] += c[0] * surfaceColor[0] * specularRefraction;
            color[1] += c[1] * surfaceColor[1] * specularRefraction;
            color[2] += c[2] * surfaceColor[2] * specularRefraction;
          }
        }
      }

      return color;
    }

    void World::photonTrace(unsigned totalPhotons)
    {
      double totalEmission = 0;
      for(std::vector<Light *>::iterator i=lightList.begin();
          i != lightList.end(); ++i)
      {
        const Vector4 &c = (*i)->getColor();
        totalEmission+=sqrt(c[0]*c[0] + c[1]*c[1] + c[1]*c[1]);
      }
    
      for(std::vector<Light *>::iterator i=lightList.begin();
          i != lightList.end(); ++i)
      {
        const Vector4 &c = (*i)->getColor();
        const double emission = sqrt(c[0]*c[0] + c[1]*c[1] + c[2]*c[2]);
        const unsigned numberPhotons = unsigned(emission/totalEmission*totalPhotons);

        Vector3 power;
        power[0] = c[0];
        power[1] = c[1];
        power[2] = c[2];
        power *= 240.0 / numberPhotons;

        for(unsigned j=0; j<numberPhotons; j++)
        {
          Vector3 v;
          drawUnitVector3(v[0], v[1], v[2]);
          Line3 ray((*i)->getPosition(), v);
          photonTrace(ray, photons, 0, 0, 0, false, true, power);
        }  
      }

      /*
        photons.push_back(Photon(Vector3(0, 0, 0), Vector3(1, 1, 1)));
        photons.push_back(Photon(Vector3(1, 0, 0), Vector3(1, 1, 1)));
        photons.push_back(Photon(Vector3(2, 0, 0), Vector3(1, 1, 1)));
        photons.push_back(Photon(Vector3(3, 0, 0), Vector3(1, 1, 1)));
        photons.push_back(Photon(Vector3(4, 0, 0), Vector3(1, 1, 1)));
      */

      std::cout << "Number of photons in map: " << photons.size() << "\n";

      for(std::vector<Photon>::iterator i=photons.begin();
          i != photons.end(); ++i) causticPhotonMap.insert(&*i);

      causticPhotonMap.balanceTree();
    }

    void World::photonTrace(const Line3 &ray, std::vector<Photon> &photons,
                            const Object *originObject,
                            const Object::Part *originPart,
                            int depth, bool transmission, 
                            bool caustic, const Vector3 &power)
    {
      const Object *hitObject = 0;
      const Object::Part *hitPart = 0;
      double hitDist = 0;

      if(transmission)
      {
        const Solid *s = dynamic_cast<const Solid *>(originObject);
        if(!s) throw InternalErrorException("Transmission initiated "
                                            "on non-solid object");
        hitObject = originObject;
        s->interiorIntersect(ray, originPart, hitDist, hitPart);
      }
      else
      {
        const Object::Part *p;
        double d;
        for(std::vector<Object *>::iterator i=objectList.begin();
            i != objectList.end(); ++i)
        {
          if((*i)->intersect(ray, originObject, originPart, d, &p) &&
             (!hitObject || d < hitDist))
          {
            hitDist = d;
            hitPart = p;
            hitObject = *i;
          }
        }
      }

      if(!hitObject) return;

      // Determine the intersection point
      Vector3 intersection = ray.direction;
      intersection *= hitDist;
      intersection += ray.point;


      // Fetch surface normal from the intersected object-part
      const Surface *surface = hitPart->getSurface();
      Vector3 normal;
      hitPart->map(intersection, normal);


      double r = Utilities::drawFloat();
      if((r-=surface->getDiffuseReflection())<0)
      {
        if(caustic && depth>=1)
          photons.push_back(Photon(intersection, power));

        if(caustic) return;
        // else random direction
      }
      else if((r-=surface->getSpecularReflection())<0)
      {
        Vector3 reflectedDirection = normal;
        reflectedDirection *= -2 * dot(ray.direction, normal);
        reflectedDirection += ray.direction;

        photonTrace(Line3(intersection, reflectedDirection), photons, hitObject, hitPart,
                    depth + 1, transmission, caustic, power);
      }
      else if((r-=surface->getSpecularRefraction())<0)
      {
        const double cosine = dot(ray.direction, transmission ? -normal : normal);
        const double n = transmission ? 1 / surface->getRefractiveIndex() : surface->getRefractiveIndex();
        const double a = 1 - (1 - cosine*cosine) / (n*n);

        if(a > 0)
        {
          Vector3 v = transmission ? -normal : normal;
          v *= (cosine/n + sqrt(a));
          Vector3 refractedDirection = ray.direction;
          refractedDirection /= n;
          refractedDirection -= v;

          photonTrace(Line3(intersection, refractedDirection), photons,
                      hitObject, hitPart, depth + 1, !transmission, 
                      caustic, power);
        }
        else
        {
          Vector3 reflectedDirection = normal;
          reflectedDirection *= -2 * dot(ray.direction, normal);
          reflectedDirection += ray.direction;

          photonTrace(Line3(intersection, reflectedDirection), photons, hitObject, hitPart,
                      depth + 1, transmission, caustic, power);
        }
      }
      else
      {
        // Absorbed.
        return;
      }
    }
  }
}

---