Edit

thodg/smallpt/smallpt.cpp

Branch :

  • smallpt.cpp 
  • /* smallpt, a Path Tracer by Kevin Beason, 2008
 * Make : g++ -O3 -fopenmp smallpt.cpp -o smallpt
 *        Remove "-fopenmp" for g++ version < 4.2
 * Usage: time ./smallpt 5000 && xv image.ppm
 */
#include <math.h>
#include <stdlib.h>
#include <stdio.h>

typedef double F;
inline F Fsqrt (F x) { return sqrt(x); }
inline F Fabs  (F x) { return fabs(x); }

struct Vec {
  F x; // position, also color (r,g,b)
  F y;
  F z;

  Vec (F x_ = 0, F y_ = 0, F z_ = 0)
  {
    x = x_;
    y = y_;
    z = z_;
  }

  Vec operator + (const Vec &b) const
  {
    return Vec(x + b.x, y + b.y, z + b.z);
  }

  Vec operator - (const Vec &b) const
  {
    return Vec(x - b.x, y - b.y, z - b.z);
  }

  Vec operator * (F b) const
  {
    return Vec(x * b, y * b, z * b);
  }

  Vec mult (const Vec &b) const
  {
    return Vec(x * b.x, y * b.y, z * b.z);
  }

  Vec & normalize () {
    return *this = *this * (1 / Fsqrt(x * x + y * y + z * z));
  }

  F dot (const Vec &b) const
  {
    return x * b.x + y * b.y + z * b.z;
  }

  // cross product
  Vec operator % (Vec &b)
  {
    return Vec(y * b.z - z * b.y,
               z * b.x - x * b.z,
               x * b.y - y * b.x);
  }
};

struct Ray {
  Vec origin;
  Vec direction;

  Ray (Vec origin_, Vec direction_)
    : origin(origin_), direction(direction_)
  {}
};

// material types, used in radiance()
enum ReflectionType {
  DIFFUSE,
  SPECULAR,
  REFRACTION
};

struct Sphere {
  F radius;            // radius
  Vec position;
  Vec emission;
  Vec color;
  ReflectionType reflection_type;
  
  Sphere (F radius_, Vec position_, Vec emission_, Vec color_,
          ReflectionType reflection_type_)
    : radius(radius_), position(position_), emission(emission_),
      color(color_), reflection_type(reflection_type_)
  {}

  // returns distance, 0 if nohit
  F intersect (const Ray &r) const
  {
    /* Solve :
     * t^2 * d . d + 2 * t * (o - p) . d + (o - p) . (o - p) - R^2 = 0
     */
    Vec op = position - r.origin;
    F t;
    F eps = 1e-4;
    F b = op.dot(r.direction);
    F det = b * b - op.dot(op) + radius * radius;
    if (det < 0)
      return 0;
    else
      det = Fsqrt(det);
    t = b - det;
    if (t > eps)
      return t;
    t = b + det;
    if (t > eps)
      return t;
    return 0;
  }
};

#define SPHERES_COUNT       9
#define SPHERES_ROT_RADIUS 20.0
#define SPHERES_ROT_SPEED   0.25

struct Scene {
  Sphere spheres[SPHERES_COUNT] = {
    Sphere(16.5,                          // Mirror ball
           Vec(50.0, 16.5, 50.0),
           Vec(),
           Vec(1, 1, 1) * 0.999,
           SPECULAR),
    Sphere(16.5,                          // Glass
           Vec(50, 16.5, 50.0),
           Vec(),
           Vec(1, 1, 1) * 0.999,
           REFRACTION),
    Sphere(1e5,                           // Left
           Vec(1e5 + 1, 40.8, 81.6),
           Vec(),
           Vec(1, 1, 1) * 0.7,
           SPECULAR),
    Sphere(1e5,                           // Right
           Vec(-1e5 + 99, 40.8, 81.6),
           Vec(),
           Vec(1, 1, 1) * 0.7,
           SPECULAR),
    Sphere(1e5,                           // Back
           Vec(50, 40.8, 1e5),
           Vec(),
           Vec(1, 1, 1) * 0.7,
           SPECULAR),
    Sphere(1e5,                           // Front
           Vec(50, 40.8, -1e5 + 170),
           Vec(),
           Vec(1, 1, 1) * 0.7,
           SPECULAR),
    Sphere(1e5,                           // Bottom
           Vec(50, 1e5, 81.6),
           Vec(),
           Vec(1, 0.9, 0.7) * 0.9,
           SPECULAR),
    Sphere(1e5,                           // Top
           Vec(50, -1e5 + 81.6, 81.6),
           Vec(),
           Vec(1, 0.9, 0.9) * 0.8,
           DIFFUSE),
    Sphere(600,                           // Light
           Vec(50, 681.6 - 0.27, 81.6),
           Vec(12, 12, 12),
           Vec(),
           DIFFUSE)
  };
  
  void update_time (F time)
  {
    F theta = time * SPHERES_ROT_SPEED * M_PI * 2.0;
    spheres[0].position = Vec(50.0, 16.5, 50.0) +
      Vec(cos(theta) * SPHERES_ROT_RADIUS, 0.0,
          sin(theta) * SPHERES_ROT_RADIUS);
    theta += M_PI;
    spheres[1].position = Vec(50.0, 16.5, 50.0) +
      Vec(cos(theta) * SPHERES_ROT_RADIUS, 0.0,
          sin(theta) * SPHERES_ROT_RADIUS);
  }
};

Scene g_scene;

inline F clamp (F x)
{
  return x < 0 ? 0 : x > 1 ? 1 : x;
}

inline int toInt (F x)
{
  return int(pow(clamp(x), 1 / 2.2) * 255 + .5);
}

inline bool intersect (const Ray &r, F &t, int &id)
{
  int i = SPHERES_COUNT;
  F d;
  F inf = 1e20;
  t = inf;
  while (i--) {
    if ((d = g_scene.spheres[i].intersect(r)) && d < t) {
      t = d;
      id = i;
    }
  }
  return t < inf;
}

Vec radiance (const Ray &r, int depth, unsigned short *Xi)
{
  F t;                                  // Distance to intersection
  int id = 0;                           // Id of intersected object
  if (! intersect(r, t, id))            // If miss, return black
    return Vec();
  const Sphere &obj = g_scene.spheres[id];      // The hit object
  Vec x = r.origin + r.direction * t;
  Vec n = (x - obj.position).normalize();
  Vec nl = n.dot(r.direction) < 0 ? n : n * -1;
  Vec f = obj.color;
  F p = f.x > f.y && f.x > f.z ? f.x : f.y > f.z ? f.y : f.z; // max refl
  if (++depth > 5) {
    if (erand48(Xi) < p)
      f = f * (1 / p);
    else
      return obj.emission; //R.R.
  }
  if (obj.reflection_type == DIFFUSE) {
    F r1 = 2 * M_PI * erand48(Xi);
    F r2 = erand48(Xi);
    F r2s = Fsqrt(r2);
    Vec w = nl;
    Vec u = ((Fabs(w.x) > .1 ? Vec(0,1) : Vec(1)) % w).normalize();
    Vec v= w % u;
    Vec d = (u * cos(r1) * r2s +
             v * sin(r1) * r2s +
             w * Fsqrt(1 - r2)).normalize();
    return obj.emission + f.mult(radiance(Ray(x, d), depth, Xi));
  }
  else if (obj.reflection_type == SPECULAR) {
    return obj.emission +
      f.mult(radiance(Ray(x, r.direction - n * 2 * n.dot(r.direction)),
                      depth, Xi));
  }
  else if (obj.reflection_type == REFRACTION) {
    Ray reflRay(x, r.direction - n * 2 * n.dot(r.direction));
    bool into = n.dot(nl) > 0;            // Ray from outside going in?
    F nc = 1;
    F nt = 1.5;
    F nnt = into ? nc / nt : nt / nc;
    F ddn = r.direction.dot(nl);
    F cos2t;
    cos2t = 1 - nnt * nnt * (1 - ddn * ddn);
    if (cos2t < 0)                        // Total internal reflection
      return obj.emission + f.mult(radiance(reflRay, depth, Xi));
    Vec tdir = (r.direction * nnt - n *
                ((into ? 1 : -1) *
                 (ddn * nnt + Fsqrt(cos2t)))).normalize();
    F a = nt - nc;
    F b = nt + nc;
    F R0 = a * a / (b * b);
    F c = 1 - (into ? -ddn : tdir.dot(n));
    F Re = R0 + (1 - R0) * c * c * c * c * c;
    F Tr = 1 - Re;
    F P = .25 + .5 * Re;
    F RP = Re / P;
    F TP = Tr / (1 - P);
    return obj.emission +
      f.mult(depth > 2 ?
             (erand48(Xi) < P ?           // Russian roulette
              radiance(reflRay, depth, Xi) * RP :
              radiance(Ray(x, tdir), depth, Xi) * TP) :
             radiance(reflRay, depth, Xi) * Re +
             radiance(Ray(x, tdir), depth, Xi) * Tr);
  }
  return Vec();
}

int main (int argc, char *argv[])
{
  int w = 800;
  int h = 600;
  int samples = argc >= 2 ? atoi(argv[1]) / 4 : 1;
  F   time = argc >= 3 ? atof(argv[2]) : 0.0;
  Ray camera(Vec(50,52,295.6), Vec(0,-0.042612,-1).normalize());
  Vec cx = Vec(w * .5135 / h);
  Vec cy = (cx % camera.direction).normalize() * .5135;
  Vec r;
  Vec *c = new Vec[w * h];
  g_scene.update_time(time);
#pragma omp parallel for schedule(dynamic, 1) private(r)    // OpenMP
  for (int y = 0; y < h; y++) {             // Loop over image rows
    fprintf(stderr, "\rRendering (%d spp) %5.2f%%",
            samples * 4,
            100.0 * y / (h - 1));
    for (unsigned short x = 0,
           Xi[3] = {0, 0, static_cast<unsigned short>(y * y * y)};
           x < w; x++) {                // Loop over image columns
      int i = (h - y - 1) * w + x;
      for (int sy = 0; sy < 2; sy++) {          // 2x2 subpixel rows
        for (int sx = 0; sx < 2; sx++) {        // 2x2 subpixel cols
          for (int s = 0; s < samples; s++){
            F r1 = 2 * erand48(Xi);
            F dx = r1 < 1 ? Fsqrt(r1) - 1 : 1 - Fsqrt(2 - r1);
            F r2 = 2 * erand48(Xi);
            F dy = r2 < 1 ? Fsqrt(r2)-1 : 1 - Fsqrt(2 - r2);
            Vec d = cx * (((sx + 0.5 + dx) / 2 + x) / w - 0.5) +
              cy * (((sy + 0.5 + dy) / 2 + y) / h - 0.5) +
              camera.direction;
            // Camera rays are pushed forward to start in interior
            r = r +
              radiance(Ray(camera.origin + d * 140,
                           d.normalize()), 0, Xi) *
              (1.0 / samples);
          }
          c[i] = c[i] + Vec(clamp(r.x), clamp(r.y), clamp(r.z)) * 0.25;
          r = Vec();
        }
      }
    }
  }
  char a[1024] = {0};
  snprintf(a, sizeof(a) - 1, "image_%09.9f.ppm", time);
  FILE *f = fopen(a, "wb");
  fprintf(f, "P3\n%d %d\n%d\n", w, h, 255);
  for (int i = 0; i < w * h; i++)
    fprintf(f, "%d %d %d ", toInt(c[i].x), toInt(c[i].y),
            toInt(c[i].z));
  delete[] c;
  delete g_scene;
}