Incremental 2D Sobol based point set

[Roquen]

Quick hack in a moment of work avoidance. Hammersley point set like 2D LDS using Sobol for one coordinate, computed incrementally.  Mimimally verified, written in Java with notes for C/C++.  The scrabbling isn't very good, but the computation is constant time instead of log(n)

Code:

  @SuppressWarnings("boxing")
  // n = number of points to dump out, s = scrambling value
  public static void quickHack(int n, int s)
  {
    // in C/C++ use unsigned
    int   si = s;
    float R = 1.f/n;
    float x,y;
    
    // hammersley style point set, sobel via graycode for second dimension
    for(int i=0; i<n; i++) {
      
      // the point S(i) = (x,y)
      x = i*R;
      /// with unsigned support lose the shift and multiply by 0x1p-32
      y = (si>>>8)*0x1p-24f;
      
      System.out.printf("%f %f\n",x, y);
      
      // update sobol value. GCC/VC have an intrinsic for trailing zeros
      // otherwise there's hacker's delight like methods
      si = si ^ 0x80000000 >>> Integer.numberOfTrailingZeros(~i);
    } 
    System.out.println();
  }

[Roquen]

Some info I skipped on the first time.  The scrambling is questionable and I haven't run any stats test to check the impact and my Z₂ knowledge isn't good enough at the theoretic level to have a feel for it.  The advantage of the above method is that it's about as fast as you can get.  It's also going to have a low star-discrepancy.  The pixar paper I mentioned in the other thread is here: Correlated Multi-Jittered Sampling  (As an aside the ad-hoc "squaring" domain mapping of the m-set that I tossed out here is mentioned in Figure 10)

The big downside of the above (and all Hammersley like sets) is that you have to decide in advance how many points you're going to use for a given computation (pixel, area, whatever) and cannot progressively add more.  Also to translating the above to GLSL requires 'findLSB'.

[Roquen]

In yet another attempt at work avoidance I dug up a progressive version...aka one where you can keep adding points until you're happy.  This time in C.  One other gotcha I forgot to mention is that neither this version nor the last deal properly with the index wrapping back to 0.  This also has a seek method which can be easily adapted to the point set version above.

Code:

#include <stdint.h>

#ifdef __GNUC__
#define ctz(x) __builtin_ctz(x)
#else
// fill-in
#endif

// sequence state data (move this and function prototypes to a header)
typedef struct
{
  uint32_t i,d0,d1;
} Sobol2D_t;


// direction table
static uint32_t D[33];

// one-time initialization
void Sobol2D_Init()
{
  uint32_t c = D[0] = 1 << 31;
  
  for(int i=1; i<32; i++) {
    c = c ^ (c >> 1);
    D[i] = c;
  }
}

// reset the sequence
void Sobol2D_Reset(Sobol2D_t* data)
{
  data->i = data->d0 = data->d1 = 0;
}

// next value in the sequence
void Sobol2D_Next(float* v, Sobol2D_t* state)
{
  uint32_t c = ctz(~state->i);

  v[0] = state->d0 * 0x1p-32f;
  v[1] = state->d1 * 0x1p-32f;
  
  state->d0 ^= 0x80000000U >> c;
  state->d1 ^= D[c];
  state->i  += 1;    
}

void Sobol2D_Seek(Sobol2D_t* state, uint32_t num)
{
  uint32_t i = state->i;
  uint32_t n = i + num;
  uint32_t a = i ^ (i >> 1);
  uint32_t b = n ^ (n >> 1);
  uint32_t d = a ^ b;
  uint32_t c = 0;
  
  while(d != 0) {
    if ((d & 1) != 0) {
      d0 ^= 0x80000000 >> c;
      d1 ^= D[c];
    }
    d >>= 1;
    c  += 1;
  }
  state->i = n;
}

#if 1
#include <stdio.h>

void main(int argc, char** argv)
{
  Sobol2D_t d;
  float r[2];

  Sobol2D_Init();
  Sobol2D_Reset(&d);
  
  for(int i=0; i<22; i++) {
    Sobol2D_Next(r, &d);
    printf("%f %f\n", r[0],r[1]);
  }
}
#endif

[Roquen]

But hey!  This is the programming sub-fourm!  Seriously this is in response to Alef's request for a Sobol reference and knightly's testing of buddhabrot rendering using a Halton sequence.

[eiffie]

I'm thankful for the code! It is true tho that a few examples would spark more interest. Sometimes people can only see the use of things when they SEE it.

[Roquen]

I know...I'm rather hoping someone else will use it and provide a demonstration. 

[Alef]

Nice, epecialy progressive.
You 'll see something only if you 'll put it in engine in place of random numbers. Throught in wikipedia there is some pictures showing smoother distribution if compared with random sequences.

[Roquen]

Tossed a full implementation (java..for C/C++ porting see notes above) of both flavors here: https://github.com/roquendm/JGO-Grabbag/tree/master/src/roquen/math/lds