Low-discrepancy samplers

[tuomlarsen]

Hello,

I'm trying to understand how low discrepancy sequences. I think I got the basic idea but I have really difficulty to make it work inside a renderer. Please, can someone explain a few things to me?

1. Should I replace all random() functions with the one of low discrepancy or is it enough to replace just the one in cosine-weighted hemisphere?

2. Should I make only one LD generator for the whole image or should I somehow reser it on every row or pixel? I mean, in case of Halton sequences I saw that sometimes the base changes (according to current row), how does this work? And in case of stratified sampler, should I generate the same numbers for every row? The same stratas for every row?

3. I saw that in case of Sobol sampler there are many "dimensions" - if I only wanted the most simple thing ever, should I use just 2 dimensions for diffuse hemisphere sampling?

I tried to read pbrt and the sources of many renderers but I still cannot make it work on my own. For example, I wanted to add Halton sampler to smallpt but when I do the whole thing goes havoc, please see the attachments (I could not attach the source code so I'm posting it below).

Thanks for any pointer or help or line of code, this really is a bit difficult for me to understand.

PS: I got the code for Halton sequence from "Analysis of the Quasi-Monte Carlo Integration of the Rendering Equation" by Laszlo Szirmay-Kalos and Werner Purgathofer.

Code: Select all

#include <math.h>   // smallpt, a Path Tracer by Kevin Beason, 2008
#include <stdlib.h> // Make : g++ -O3 -fopenmp smallpt.cpp -o smallpt
#include <stdio.h>  //        Remove "-fopenmp" for g++ version < 4.2
class Halton {
    double value, inv_base;
public:
    void Number( long i, int base ) {
        double f = inv_base = 1.0/base;
        value = 0.0;
        while ( i > 0 ) {
            value += f * (double)(i % base);
            i /= base;
            f *= inv_base;
        }
    }
    void Next( ) {
        double r = 1.0 - value - 0.0000000001;
        if (inv_base < r) value += inv_base;
        else {
            double h = inv_base, hh;
            do {
                hh = h;
                h *= inv_base;
            } while ( h >= r );
            value += hh + h - 1.0;
        }
    }
    double get( ) { return value; }
    operator double( ) { return value; }
};
struct Vec {        // Usage: time ./smallpt 5000 && xv image.ppm
  double x, y, z;                  // position, also color (r,g,b)
  Vec(double x_=0, double y_=0, double 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*(double 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& norm(){ return *this = *this * (1/sqrt(x*x+y*y+z*z)); }
  double dot(const Vec &b) const { return x*b.x+y*b.y+z*b.z; } // cross:
  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 o, d; Ray(Vec o_, Vec d_) : o(o_), d(d_) {} };
enum Refl_t { DIFF, SPEC, REFR };  // material types, used in radiance()
struct Sphere {
  double rad;       // radius
  Vec p, e, c;      // position, emission, color
  Refl_t refl;      // reflection type (DIFFuse, SPECular, REFRactive)
  Sphere(double rad_, Vec p_, Vec e_, Vec c_, Refl_t refl_):
    rad(rad_), p(p_), e(e_), c(c_), refl(refl_) {}
  double intersect(const Ray &r) const { // returns distance, 0 if nohit
    Vec op = p-r.o; // Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0
    double t, eps=1e-4, b=op.dot(r.d), det=b*b-op.dot(op)+rad*rad;
    if (det<0) return 0; else det=sqrt(det);
    return (t=b-det)>eps ? t : ((t=b+det)>eps ? t : 0);
  }
};
Sphere spheres[] = {//Scene: radius, position, emission, color, material
  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 double clamp(double x){ return x<0 ? 0 : x>1 ? 1 : x; }
inline int toInt(double x){ return int(pow(clamp(x),1/2.2)*255+.5); }
inline bool intersect(const Ray &r, double &t, int &id){
  double n=sizeof(spheres)/sizeof(Sphere), d, inf=t=1e20;
  for(int i=int(n);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){
Vec radiance(const Ray &r, int depth, unsigned short *Xi, Halton h1, Halton h2){
  double t;                               // distance to intersection
  int id=0;                               // id of intersected object
  if (!intersect(r, t, id)) return Vec(); // if miss, return black
  const Sphere &obj = spheres[id];        // the hit object
  Vec x=r.o+r.d*t, n=(x-obj.p).norm(), nl=n.dot(r.d)<0?n:n*-1, f=obj.c;
  double 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.e; //R.R.
  if (obj.refl == DIFF){                  // Ideal DIFFUSE reflection
//    double r1=2*M_PI*erand48(Xi), r2=erand48(Xi), r2s=sqrt(r2);
    h1.Next();
    h2.Next();
    double r1=2*M_PI*(double)h1, r2=(double)h2, r2s=sqrt(r2);
    Vec w=nl, u=((fabs(w.x)>.1?Vec(0,1):Vec(1))%w).norm(), v=w%u;
    Vec d = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2)).norm();
//    return obj.e + f.mult(radiance(Ray(x,d),depth,Xi));
    return obj.e + f.mult(radiance(Ray(x,d),depth,Xi,h1,h2));
  } else if (obj.refl == SPEC)            // Ideal SPECULAR reflection
//    return obj.e + f.mult(radiance(Ray(x,r.d-n*2*n.dot(r.d)),depth,Xi));
    return obj.e + f.mult(radiance(Ray(x,r.d-n*2*n.dot(r.d)),depth,Xi,h1,h2));
  Ray reflRay(x, r.d-n*2*n.dot(r.d));     // Ideal dielectric REFRACTION
  bool into = n.dot(nl)>0;                // Ray from outside going in?
  double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t;
  if ((cos2t=1-nnt*nnt*(1-ddn*ddn))<0)    // Total internal reflection
//    return obj.e + f.mult(radiance(reflRay,depth,Xi));
    return obj.e + f.mult(radiance(reflRay,depth,Xi,h1,h2));
  Vec tdir = (r.d*nnt - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).norm();
  double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c = 1-(into?-ddn:tdir.dot(n));
  double Re=R0+(1-R0)*c*c*c*c*c,Tr=1-Re,P=.25+.5*Re,RP=Re/P,TP=Tr/(1-P);
  return obj.e + 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);
    radiance(reflRay,depth,Xi,h1,h2)*RP:radiance(Ray(x,tdir),depth,Xi,h1,h2)*TP) :
    radiance(reflRay,depth,Xi,h1,h2)*Re+radiance(Ray(x,tdir),depth,Xi,h1,h2)*Tr);
}
int main(int argc, char *argv[]){
  Halton h1, h2;
  h1.Number(0, 2);
  h2.Number(0, 3);
  int w=256, h=256, samps = argc==2 ? atoi(argv[1])/4 : 1; // # samples
  Ray cam(Vec(50,52,295.6), Vec(0,-0.042612,-1).norm()); // cam pos, dir
  Vec cx=Vec(w*.5135/h), cy=(cx%cam.d).norm()*.5135, r, *c=new Vec[w*h];
//#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%%",samps*4,100.*y/(h-1));
    for (unsigned short x=0, Xi[3]={0,0,y*y*y}; x<w; x++)   // Loop cols
      for (int sy=0, i=(h-y-1)*w+x; sy<2; sy++)     // 2x2 subpixel rows
        for (int sx=0; sx<2; sx++, r=Vec()){        // 2x2 subpixel cols
          for (int s=0; s<samps; s++){
            double r1=2*erand48(Xi), dx=r1<1 ? sqrt(r1)-1: 1-sqrt(2-r1);
            double r2=2*erand48(Xi), dy=r2<1 ? sqrt(r2)-1: 1-sqrt(2-r2);
            Vec d = cx*( ( (sx+.5 + dx)/2 + x)/w - .5) +
                    cy*( ( (sy+.5 + dy)/2 + y)/h - .5) + cam.d;
//            r = r + radiance(Ray(cam.o+d*140,d.norm()),0,Xi)*(1./samps);
            r = r + radiance(Ray(cam.o+d*140,d.norm()),0,Xi,h1,h2)*(1./samps);
          } // Camera rays are pushed ^^^^^ forward to start in interior
          c[i] = c[i] + Vec(clamp(r.x),clamp(r.y),clamp(r.z))*.25;
        }
  }
  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));
}


[jbikker]

I am working on this right now, so let me share a few bits. First of all, as suggested by "Megakernels Considered Harmful - Wavefront Path Tracing on GPUs" (Laine, Karras & Aila, 2013) I consulted "Constructing Sobol Sequences with Better Two-dimensional Projections" (Joe & Kuo, 2008), as well as some code on Kuo's page (http://web.maths.unsw.edu.au/~fkuo/sobol). I reduced this to the following easy-to-use class:

Code: Select all

class SobolGenerator
{
public:
   SobolGenerator()
   {
      memset( X, 0, sizeof( X ) );
   }
   void Generate( unsigned N )
   {
      unsigned int s[32] = { 0,  1,  2,  3,  3,  4,  4,  5,  5,  5,  5,  5,  5,  6,  6,  6,  6,  6,  6,  7,  7,  7,  7,  7,  7,  7,  7,  7,  7,  7,  7,  7 };
      unsigned int a[32] = { 0,  0,  1,  1,  2,  1,  4,  2, 13,  7, 14, 11,  4,  1, 16, 13, 22, 19, 25,  1, 32,  4,  8,  7, 56, 14, 28, 19, 50, 21, 42, 31 };
      unsigned int m[32][8] = { { 0, 1, 0, 0, 0, 0, 0, 0 }, { 0, 1, 1, 0, 0, 0, 0, 0 }, { 0, 1, 1, 1, 0, 0, 0, 0 }, { 0, 1, 3, 1, 0, 0, 0, 0 },
                          { 0, 1, 1, 7, 13, 0, 0, 0 }, { 0, 1, 1, 3, 7, 0, 0, 0 }, { 0, 1, 3, 1, 7, 21, 0, 0 }, { 0, 1, 3, 1, 3, 9, 0, 0 },
                          { 0, 1, 1, 5, 9, 13, 0, 0 }, { 0, 1, 1, 3, 9, 13, 0, 0 }, { 0, 1, 1, 5, 3, 7, 0, 0 }, { 0, 1, 1, 5, 7, 11, 0, 0 },
                          { 0, 1, 3, 3, 13, 15, 43, 0 }, { 0, 1, 3, 5, 11, 25, 45, 0 }, { 0, 1, 1, 3, 11, 3, 45, 0 }, { 0, 1, 3, 5, 1, 9, 21, 0 },
                          { 0, 1, 1, 3, 13, 9, 9, 0 }, { 0, 1, 3, 5, 7, 17, 53, 0 }, { 0, 1, 3, 1, 7, 11, 51, 115 }, { 0, 1, 1, 7, 9, 25, 35, 11 },
                          { 0, 1, 3, 1, 1, 31, 5, 1 }, { 0, 1, 1, 5, 9, 11, 1, 121 }, { 0, 1, 1, 1, 15, 11, 59, 21 }, { 0, 1, 3, 1, 3, 17, 49, 51 },
                          { 0, 1, 1, 5, 7, 15, 25, 13 }, { 0, 1, 3, 7, 7, 1, 45, 7 }, { 0, 1, 1, 5, 7, 21, 21, 37 }, { 0, 1, 3, 1, 13, 9, 49, 23 },
                          { 0, 1, 1, 1, 3, 3, 35, 123 }, { 0, 1, 1, 7, 5, 15, 47, 117 }, { 0, 1, 1, 3, 13, 9, 23, 33 } };
      unsigned int V[33], C = 1, value = N;
      const float rpf = 1.0f / powf( 2.0f, 32 );
      while (value & 1) value >>= 1, C++;
      for (unsigned int i = 1; i <= 32; i++ ) V[i] = 1 << (32 - i); // all m's = 1
      X[0] = X[0] ^ V[C], points[0] = X[0] * rpf;
      for ( unsigned int j = 1; j <= 31; X[j] = X[j] ^ V[C], points[j] = X[j] * rpf, j++ )
         if (32 <= s[j]) for( unsigned int i = 1; i <= 32; i++ ) V[i] = m[j][i] << (32 - i); else
         {
            for( unsigned int i = 1; i <= s[j]; i++ ) V[i] = m[j][i] << (32 - i);
            for( unsigned int i = s[j] + 1; i <= 32; i++ )
            {
               V[i] = V[i-s[j]] ^ (V[i - s[j]] >> s[j]);
               for( unsigned int k = 1; k <= s[j] - 1; k++ ) V[i] ^= (((a[j] >> (s[j] - 1 - k)) & 1) * V[i - k]);
            }
         }
   }
   float points[32];
   unsigned int X[32];
};

which is used as follows:
1) Instantiate a SobolGenerator (e.g., SobolGenerator sobol; );
2) For each path call 'sobol.Generate( frameNumber );';
3) After that, use the points on 'sobol.point[]' to replace your random number generator calls.
In your case, you would only be using sobol.point[0] and sobol.point[1].
The above code results in a very nice set of pseudo-random numbers; I tried many pairs out of the point array to plot random points in 2D, and although there are some apparent patterns, it's pretty decent.
If you use this in a path tracer, you should use the same set of 32 values for each pixel. To prevent each pixel using the same values, you xor the retrieved points with a fixed random number based on pixel index.
The above is assuming you do one path per pixel at a time, if you do multiple paths per pixel, step 2 becomes something like 'sobol.Generate( frameNumber * spp + currentSample )'.

I hope this makes sense,
- Jacco.

EDIT: Note that I am assuming here that 32 dimensions is sufficient for you. In my case it typically is; if you need more than 32 random numbers either see the link I provided for more data, or switch to a random number generator after dimension 32, as proposed by Lain et al. Regarding that paper: I recommend reading section 4 (starting at "For generating low-discrepancy quasirandom numbers needed in the samplers...") for additional insight.

[ingenious]

To prevent each pixel using the same values, you xor the retrieved points with a fixed random number based on pixel index.

So you XOR the bits of the float number returned by the Sobol generator with an integer computed as a hash of the pixel index? Are you sure this preserves the good stratification properties of the sequence?

What I do is a Cranley-Patterson rotation. That is, instead of XOR-ing the Sobol number, I add an offset to it, which again is computed from the pixel index. The final number is (SobolNumber + offset) % 1. I also found that it's good to use a different offset for every dimension. So what I do is, I initialize a simple congruential random generator seeding it with the pixel index, and then for every new dimension I take the next number from that generator as an offset. This has worked very well for me, and is still very fast.

The best way to go is to actually consider the image plane as a whole and just shoot paths through it using the Sobol sequence. Then you don't need any offsets or XORs because every pixel gets a different path. The one disadvantage of this approach is that the pixels are traversed in an guaranteed incoherent order For the Halton sequence, Gruenschloss et al. showed how to traverse the pixels in arbitrary order (see the most recent Quasi MC Siggraph course notes), but I'm not sure this has been done for the Sobol sequence. Does anybody know?

[tuomlarsen]

Jacco, thanks many times for the code you posted!

But I must be really stupid because even that you posted step-by-step instructions I cannot make it work.

1) Instantiate a SobolGenerator (e.g., SobolGenerator sobol; );

Ok, easy.

2) For each path call 'sobol.Generate( frameNumber );

smallpt uses several samples and it divides a pixel into 2x2 grid. So that would be:

Code: Select all

sobol.Generate(samps*(sy*2+sx)+s);

right?

3) After that, use the points on 'sobol.point[]' to replace your random number generator calls.

I guess just like:

Code: Select all

//double r1=2*M_PI*erand48(Xi), r2=erand48(Xi), r2s=sqrt(r2);
double r1=2*M_PI*sobol.points[0], r2=sobol.points[1], r2s=sqrt(r2);


After this I get the sobol1.png attached.

What the XORing goes, I'm not sure I followed because when I do:

Code: Select all

int index = (y*2+sy)*w+(x*2+sx);
...
double r1=2*M_PI*(float)((int)sobol.points[0]^index)
double r2=(float)((int)sobol.points[1]^index)


I get almost everything black (sobol2.png)

Also, why does the SobolGenerator uses so many dimensions (32)? I understand that there might be some case where one needs so high dimensions (like in statistics) but for path-tracing it is way less, no? I mean, I try to use only two (for hemisphere sampling), another two could be used for antialiasing, one for camera time, ... and whatnot but definitely way less than 32. Would it be faster if I knew in advance that I only want, say, two dimensions?

On related note, should I Sobol-/Halton-/...-sample only the hemisphere or also for the camera ray jittering? What I'm after is renderings which have less perceivable noise...

spectral, could you please tell me how to calculate the "offset"? Is it something like pixel position (x+y*width) divided by pixel count (width*height)?

Second try:

Code: Select all

#include <math.h>   // smallpt, a Path Tracer by Kevin Beason, 2008
#include <stdlib.h> // Make : g++ -O3 -fopenmp smallpt.cpp -o smallpt
#include <stdio.h>  //        Remove "-fopenmp" for g++ version < 4.2
#include <string.h>
class Halton {
    double value, inv_base;
public:
   void Number( long i, int base ) {
      double f = inv_base = 1.0/base;
      value = 0.0;
      while ( i > 0 ) {
         value += f * (double)(i % base);
         i /= base;
         f *= inv_base;
      }
   }
   void Next( ) {
      double r = 1.0 - value - 0.0000000001;
      if (inv_base < r) value += inv_base;
      else {
         double h = inv_base, hh;
         do {
            hh = h;
            h *= inv_base;
         } while ( h >= r );
         value += hh + h - 1.0;
      }
   }
   double get( ) { return value; }
   operator double( ) { return value; }
};

class SobolGenerator
{
public:
   SobolGenerator()
   {
      memset( X, 0, sizeof( X ) );
   }
   void Generate( unsigned N )
   {
      unsigned int s[32] = { 0,  1,  2,  3,  3,  4,  4,  5,  5,  5,  5,  5,  5,  6,  6,  6,  6,  6,  6,  7,  7,  7,  7,  7,  7,  7,  7,  7,  7,  7,  7,  7 };
      unsigned int a[32] = { 0,  0,  1,  1,  2,  1,  4,  2, 13,  7, 14, 11,  4,  1, 16, 13, 22, 19, 25,  1, 32,  4,  8,  7, 56, 14, 28, 19, 50, 21, 42, 31 };
      unsigned int m[32][8] = { { 0, 1, 0, 0, 0, 0, 0, 0 }, { 0, 1, 1, 0, 0, 0, 0, 0 }, { 0, 1, 1, 1, 0, 0, 0, 0 }, { 0, 1, 3, 1, 0, 0, 0, 0 },
                          { 0, 1, 1, 7, 13, 0, 0, 0 }, { 0, 1, 1, 3, 7, 0, 0, 0 }, { 0, 1, 3, 1, 7, 21, 0, 0 }, { 0, 1, 3, 1, 3, 9, 0, 0 },
                          { 0, 1, 1, 5, 9, 13, 0, 0 }, { 0, 1, 1, 3, 9, 13, 0, 0 }, { 0, 1, 1, 5, 3, 7, 0, 0 }, { 0, 1, 1, 5, 7, 11, 0, 0 },
                          { 0, 1, 3, 3, 13, 15, 43, 0 }, { 0, 1, 3, 5, 11, 25, 45, 0 }, { 0, 1, 1, 3, 11, 3, 45, 0 }, { 0, 1, 3, 5, 1, 9, 21, 0 },
                          { 0, 1, 1, 3, 13, 9, 9, 0 }, { 0, 1, 3, 5, 7, 17, 53, 0 }, { 0, 1, 3, 1, 7, 11, 51, 115 }, { 0, 1, 1, 7, 9, 25, 35, 11 },
                          { 0, 1, 3, 1, 1, 31, 5, 1 }, { 0, 1, 1, 5, 9, 11, 1, 121 }, { 0, 1, 1, 1, 15, 11, 59, 21 }, { 0, 1, 3, 1, 3, 17, 49, 51 },
                          { 0, 1, 1, 5, 7, 15, 25, 13 }, { 0, 1, 3, 7, 7, 1, 45, 7 }, { 0, 1, 1, 5, 7, 21, 21, 37 }, { 0, 1, 3, 1, 13, 9, 49, 23 },
                          { 0, 1, 1, 1, 3, 3, 35, 123 }, { 0, 1, 1, 7, 5, 15, 47, 117 }, { 0, 1, 1, 3, 13, 9, 23, 33 } };
      unsigned int V[33], C = 1, value = N;
      const float rpf = 1.0f / powf( 2.0f, 32 );
      while (value & 1) value >>= 1, C++;
      for (unsigned int i = 1; i <= 32; i++ ) V[i] = 1 << (32 - i); // all m's = 1
      X[0] = X[0] ^ V[C], points[0] = X[0] * rpf;
      for ( unsigned int j = 1; j <= 31; X[j] = X[j] ^ V[C], points[j] = X[j] * rpf, j++ )
         if (32 <= s[j]) for( unsigned int i = 1; i <= 32; i++ ) V[i] = m[j][i] << (32 - i); else
         {
            for( unsigned int i = 1; i <= s[j]; i++ ) V[i] = m[j][i] << (32 - i);
            for( unsigned int i = s[j] + 1; i <= 32; i++ )
            {
               V[i] = V[i-s[j]] ^ (V[i - s[j]] >> s[j]);
               for( unsigned int k = 1; k <= s[j] - 1; k++ ) V[i] ^= (((a[j] >> (s[j] - 1 - k)) & 1) * V[i - k]);
            }
         }
   }
   float points[32];
   unsigned int X[32];
};


struct Vec {        // Usage: time ./smallpt 5000 && xv image.ppm
  double x, y, z;                  // position, also color (r,g,b)
  Vec(double x_=0, double y_=0, double 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*(double 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& norm(){ return *this = *this * (1/sqrt(x*x+y*y+z*z)); }
  double dot(const Vec &b) const { return x*b.x+y*b.y+z*b.z; } // cross:
  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 o, d; Ray(Vec o_, Vec d_) : o(o_), d(d_) {} };
enum Refl_t { DIFF, SPEC, REFR };  // material types, used in radiance()
struct Sphere {
  double rad;       // radius
  Vec p, e, c;      // position, emission, color
  Refl_t refl;      // reflection type (DIFFuse, SPECular, REFRactive)
  Sphere(double rad_, Vec p_, Vec e_, Vec c_, Refl_t refl_):
    rad(rad_), p(p_), e(e_), c(c_), refl(refl_) {}
  double intersect(const Ray &r) const { // returns distance, 0 if nohit
    Vec op = p-r.o; // Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0
    double t, eps=1e-4, b=op.dot(r.d), det=b*b-op.dot(op)+rad*rad;
    if (det<0) return 0; else det=sqrt(det);
    return (t=b-det)>eps ? t : ((t=b+det)>eps ? t : 0);
  }
};
Sphere spheres[] = {//Scene: radius, position, emission, color, material
  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 double clamp(double x){ return x<0 ? 0 : x>1 ? 1 : x; }
inline int toInt(double x){ return int(pow(clamp(x),1/2.2)*255+.5); }
inline bool intersect(const Ray &r, double &t, int &id){
  double n=sizeof(spheres)/sizeof(Sphere), d, inf=t=1e20;
  for(int i=int(n);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){
//Vec radiance(const Ray &r, int depth, unsigned short *Xi, Halton h1, Halton h2){
Vec radiance(const Ray &r, int depth, unsigned short *Xi, SobolGenerator sobol, int index){
  double t;                               // distance to intersection
  int id=0;                               // id of intersected object
  if (!intersect(r, t, id)) return Vec(); // if miss, return black
  const Sphere &obj = spheres[id];        // the hit object
  Vec x=r.o+r.d*t, n=(x-obj.p).norm(), nl=n.dot(r.d)<0?n:n*-1, f=obj.c;
  double 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.e; //R.R.
  if (obj.refl == DIFF){                  // Ideal DIFFUSE reflection
//    double r1=2*M_PI*erand48(Xi), r2=erand48(Xi), r2s=sqrt(r2);
//   h1.Next();
//   h2.Next();
//    double r1=2*M_PI*(double)h1, r2=(double)h2, r2s=sqrt(r2);
//   double r1=2*M_PI*sobol.points[0], r2=sobol.points[1], r2s=sqrt(r2);
   double r1=2*M_PI*(float)((int)sobol.points[0]^index), r2=(float)((int)sobol.points[1]^index), r2s=sqrt(r2);
    Vec w=nl, u=((fabs(w.x)>.1?Vec(0,1):Vec(1))%w).norm(), v=w%u;
    Vec d = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2)).norm();
//    return obj.e + f.mult(radiance(Ray(x,d),depth,Xi));
//    return obj.e + f.mult(radiance(Ray(x,d),depth,Xi,h1,h2));
    return obj.e + f.mult(radiance(Ray(x,d),depth,Xi,sobol,index));
  } else if (obj.refl == SPEC)            // Ideal SPECULAR reflection
//    return obj.e + f.mult(radiance(Ray(x,r.d-n*2*n.dot(r.d)),depth,Xi));
//    return obj.e + f.mult(radiance(Ray(x,r.d-n*2*n.dot(r.d)),depth,Xi,h1,h2));
    return obj.e + f.mult(radiance(Ray(x,r.d-n*2*n.dot(r.d)),depth,Xi,sobol,index));
  Ray reflRay(x, r.d-n*2*n.dot(r.d));     // Ideal dielectric REFRACTION
  bool into = n.dot(nl)>0;                // Ray from outside going in?
  double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t;
  if ((cos2t=1-nnt*nnt*(1-ddn*ddn))<0)    // Total internal reflection
//    return obj.e + f.mult(radiance(reflRay,depth,Xi));
//    return obj.e + f.mult(radiance(reflRay,depth,Xi,h1,h2));
    return obj.e + f.mult(radiance(reflRay,depth,Xi,sobol,index));
  Vec tdir = (r.d*nnt - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).norm();
  double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c = 1-(into?-ddn:tdir.dot(n));
  double Re=R0+(1-R0)*c*c*c*c*c,Tr=1-Re,P=.25+.5*Re,RP=Re/P,TP=Tr/(1-P);
  return obj.e + 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);
//    radiance(reflRay,depth,Xi,h1,h2)*RP:radiance(Ray(x,tdir),depth,Xi,h1,h2)*TP) :
//    radiance(reflRay,depth,Xi,h1,h2)*Re+radiance(Ray(x,tdir),depth,Xi,h1,h2)*Tr);
    radiance(reflRay,depth,Xi,sobol,index)*RP:radiance(Ray(x,tdir),depth,Xi,sobol,index)*TP) :
    radiance(reflRay,depth,Xi,sobol,index)*Re+radiance(Ray(x,tdir),depth,Xi,sobol,index)*Tr);
}
int main(int argc, char *argv[]){
  //Halton h1, h2;
  //h1.Number(0, 2);
  //h2.Number(0, 3);
  SobolGenerator sobol;
  int w=256, h=256, samps = argc==2 ? atoi(argv[1])/4 : 1; // # samples
  Ray cam(Vec(50,52,295.6), Vec(0,-0.042612,-1).norm()); // cam pos, dir
  Vec cx=Vec(w*.5135/h), cy=(cx%cam.d).norm()*.5135, r, *c=new Vec[w*h];
//#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%%",samps*4,100.*y/(h-1));
    for (unsigned short x=0, Xi[3]={0,0,y*y*y}; x<w; x++)   // Loop cols
      for (int sy=0, i=(h-y-1)*w+x; sy<2; sy++)     // 2x2 subpixel rows
        for (int sx=0; sx<2; sx++, r=Vec()){        // 2x2 subpixel cols
          for (int s=0; s<samps; s++){
            double r1=2*erand48(Xi), dx=r1<1 ? sqrt(r1)-1: 1-sqrt(2-r1);
            double r2=2*erand48(Xi), dy=r2<1 ? sqrt(r2)-1: 1-sqrt(2-r2);
            Vec d = cx*( ( (sx+.5 + dx)/2 + x)/w - .5) +
                    cy*( ( (sy+.5 + dy)/2 + y)/h - .5) + cam.d;
//            r = r + radiance(Ray(cam.o+d*140,d.norm()),0,Xi)*(1./samps);
//            r = r + radiance(Ray(cam.o+d*140,d.norm()),0,Xi,h1,h2)*(1./samps);
         sobol.Generate(samps*(sy*2+sx)+s);
         int index = (y*2+sy)*w+(x*2+sx);
         //for (int i=0; i<32; i++){
         //   sobol.points[i] = (float)((int)sobol.points[i] ^ index);
         //}
            r = r + radiance(Ray(cam.o+d*140,d.norm()),0,Xi,sobol,index)*(1./samps);
          } // Camera rays are pushed ^^^^^ forward to start in interior
          c[i] = c[i] + Vec(clamp(r.x),clamp(r.y),clamp(r.z))*.25;
        }
  }
  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));
}


[tuomlarsen]

In other words, I more or less "get" that it is somehow possible to generate sequence of numbers which are nicely apart of each other. The problem I have, even if I generate such numbers and visually verify that they do approximately what I want (by generating many points on a plane), I still cannot make it work with a renderer. Because every time I replace proper PRNG with that sequence I either get everything black or weird artifacts...

[toxie]

The best way to go is to actually consider the image plane as a whole and just shoot paths through it using the Sobol sequence. Then you don't need any offsets or XORs because every pixel gets a different path. The one disadvantage of this approach is that the pixels are traversed in an guaranteed incoherent order For the Halton sequence, Gruenschloss et al. showed how to traverse the pixels in arbitrary order (see the most recent Quasi MC Siggraph course notes), but I'm not sure this has been done for the Sobol sequence. Does anybody know?

Yup, there's also source released by him, see http://gruenschloss.org/ (Enumerating Quasi-Monte Carlo Point Sequences in Elementary Intervals)

[friedlinguini]

Also, why does the SobolGenerator uses so many dimensions (32)? I understand that there might be some case where one needs so high dimensions (like in statistics) but for path-tracing it is way less, no? I mean, I try to use only two (for hemisphere sampling), another two could be used for antialiasing, one for camera time, ... and whatnot but definitely way less than 32. Would it be faster if I knew in advance that I only want, say, two dimensions?

I suspect this is be the root of your problem. Every random number you use within a sample should correspond to a different dimension. In your current code, for diffuse reflection you evaluate a direction based on sobol.points[0] and sobol.points[1]. For a given path, you will wind up with the same reflected direction (in tangent space) for every bounce.

What you probably want is something like...

Code: Select all

int dimension;
...
double r1=2*M_PI*(float)((int)sobol.points[dimension++]^index)
double r2=(float)((int)sobol.points[dimension++]^index)
...
for (int s=0; s<samps; s++){
    dimension = 0;

[tuomlarsen]

Thank you, friedlinguini! That tiny-bit clears things up - but what about, say, stratified sampling? What corresponds to sampling dimension in stratified sampling? I mean, with Halton/Sobol one uses ever higher primes with every light bounce - but in case of stratified sampling..what do the stratas correspond to actually? Is it possible to use stratified sampling with Russian roulette (as in, one does not know the path lengths in advance)?

Btw, the code I posted earlier is definitely wrong:

Code: Select all

double r2=(float)((int)sobol.points[dimension++]^index)

This still does not work perfectly but it's a bit better:

Code: Select all

double r2=(double)((unsigned int)(sobol.points[dimension++]*w*h*4.0)^index)/(w*h*4.0);

[friedlinguini]

Thank you, friedlinguini! That tiny-bit clears things up - but what about, say, stratified sampling? What corresponds to sampling dimension in stratified sampling? I mean, with Halton/Sobol one uses ever higher primes with every light bounce - but in case of stratified sampling..what do the stratas correspond to actually? Is it possible to use stratified sampling with Russian roulette (as in, one does not know the path lengths in advance)?

What about stratified sampling? You can replace a stratified sample value with a number from a low-discrepancy sequence. E.g.,

Naive:

Code: Select all

for (int i = 0; i < M * N; ++i) {
    r0 = random();
    r1 = random();
    ...

Stratified:

Code: Select all

for (int i = 0; i < M; ++i) {
    for (int j = 0; j < N; ++j) {
        r0 = (i + random()) / M;
        r1 = (j + random()) / N;
    ...

Sobol:

Code: Select all

SobolGenerator sobol;
for (int i = 0; i < M * N; ++i) {
    sobol.Generate(frameNumber * M * N + i);
    dimension = 0;
    r0 = sobol.point[dimension++];
    r1 = sobol.point[dimension++];
    ...

Samples taken from low-discrepancy sequences should just naturally interact nicely in multiple dimensions in much the same way that stratified samples do. I'm not sure I understand the question otherwise.

Yes, you can use it with Russian roulette as well, though the benefits diminish after the first couple of bounces. The benefits of low-discrepancy sampling aren't particularly visible in tertiary rays. That's one reason (in addition to the storage costs involved) that the Sobol generator falls back to ordinary random number generation after 32 dimensions. Halton sequences could theoretically be used for an unbounded number of dimensions, but if you don't have enough samples and don't take care to scramble the numbers, the correlation between dimensions outweighs the theoretical low discrepancy (http://marcoagd.usuarios.rdc.puc-rio.br ... _20x21.gif).

[spectral]

In other words, I more or less "get" that it is somehow possible to generate sequence of numbers which are nicely apart of each other. The problem I have, even if I generate such numbers and visually verify that they do approximately what I want (by generating many points on a plane), I still cannot make it work with a renderer. Because every time I replace proper PRNG with that sequence I either get everything black or weird artifacts...

So, notice that:
1) Use uniform RNG can "hide" some problems that you have in your renderer...
2) First, try to play with QMC sequence and check the result you got. By example, you can use the http://www.processing.or tool to play with.
I join a simple example of processing code (Java like) that play with QMC sequence. (It is an easy way to test it...)
3) Notice that QMC sequences are only interesting for the first 1000 iterations... (approx.) but it can help for better interactivity or better image approximation for the first MLT phase
4) Take a look at the following paper, seems interesting : http://ompf2.com/viewtopic.php?f=6&p=3849#p3849
5) I definitively agree with Toxie ... check : http://gruenschloss.org/

Hope it helps