kmx git

/* 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 Refl_t { DIFF, SPEC, REFR };

struct Sphere {
  F radius;            // radius
  Vec position;
  Vec emission;
  Vec color;
  Refl_t reflection_type;
  
  Sphere (F radius_, Vec position_, Vec emission_, Vec color_,
          Refl_t reflection_type_)
    : radius(radius_), p(p_), e(e_), c(c_), refl(refl_)
  {}

  // 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 = p - 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);
    return (t = b-det) > eps ? t : ((t = b + det) > eps ? t : 0);
  }
};

//Scene: radius, position, emission, color, material
Sphere g_spheres[] = {
  Sphere(1e5, Vec( 1e5+1,40.8,81.6), Vec(),Vec(.75,.25,.25),DIFF),//Left
  Sphere(1e5, Vec(-1e5+99,40.8,81.6),Vec(),Vec(.25,.25,.75),DIFF),//Rght
  Sphere(1e5, Vec(50,40.8, 1e5),     Vec(),Vec(.75,.75,.75),DIFF),//Back
  Sphere(1e5, Vec(50,40.8,-1e5+170), Vec(),Vec(),           DIFF),//Frnt
  Sphere(1e5, Vec(50, 1e5, 81.6),    Vec(),Vec(.75,.75,.75),DIFF),//Botm
  Sphere(1e5, Vec(50,-1e5+81.6,81.6),Vec(),Vec(.75,.75,.75),DIFF),//Top
  Sphere(16.5,Vec(27,16.5,47),       Vec(),Vec(1,1,1)*.999, SPEC),//Mirr
  Sphere(16.5,Vec(73,16.5,78),       Vec(),Vec(1,1,1)*.999, REFR),//Glas
  Sphere(600, Vec(50,681.6-.27,81.6),Vec(12,12,12),  Vec(), DIFF) //Lite
};

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 = sizeof(g_spheres) / sizeof(Sphere);
  F d;
  F inf = 1e20;
  t = inf;
  while (i--) {
    if ((d = 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 = 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.c;
  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 == DIFF) {    // Ideal diffuse reflection
    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.refl == SPEC)            // Ideal specular reflection
    return obj.emission +
      f.mult(radiance(Ray(x, r.direction - n * 2 * n.dot(r.direction)),
                      depth, Xi));
  else {                                // Ideal dielectric 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);
}

int main (int argc, char *argv[])
{
  int w = 1024;
  int h = 768;
  int samples = argc == 2 ? atoi(argv[1]) / 4 : 1;
  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];
#pragma omp parallel for schedule(dynamic, 1) private(r)    // OpenMP
  unsigned short x;
  unsigned short Xi[3];
  int y = 0;
  while (y < h) {                       // Loop over image rows
    fprintf(stderr, "\rRendering (%d spp) %5.2f%%",
            samples * 4,
            100.0 * y / (h - 1));
    x = 0;
    *Xi = (unsigned short [])
      {0, 0, static_cast<unsigned short>(y * y * y)};
    while (x < w) {                     // Loop over image columns
      int sy = 0;
      int i = (h - y - 1) * w + x;
      while (sy < 2) {                  // 2x2 subpixel rows
        int sx = 0;
        while (sx < 2) {                // 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();
          sx++;
        }
        sy++;
      }
      x++;
    }
    y++;
  }
  FILE *f = fopen("image.ppm", "w");    // Write image to PPM file.
  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));
}
thodg/smallpt/smallpt.cpp

Commit

smallpt.cpp
