Welcome to Fractal Forums

Fractal Software => Programming => Topic started by: Roquen on October 11, 2013, 11:04:09 AM




Title: Incremental 2D Sobol based point set
Post by: Roquen on October 11, 2013, 11:04:09 AM
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();
  }


Title: Re: Incremental 2D Sobol based point set
Post by: Roquen on October 14, 2013, 11:25:26 AM
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 Z2 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 (http://graphics.pixar.com/library/MultiJitteredSampling/)  (As an aside the ad-hoc "squaring" domain mapping of the m-set that I tossed out here (http://www.fractalforums.com/fragmentarium/fractal-plus-change-of-domain-mashup/msg65054/#msg65054) 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'.


Title: Re: Incremental 2D Sobol based point set
Post by: Roquen on October 15, 2013, 11:09:56 PM
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


Title: Re: Incremental 2D Sobol based point set
Post by: Roquen on October 16, 2013, 10:50:51 AM
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.


Title: Re: Incremental 2D Sobol based point set
Post by: eiffie on October 17, 2013, 06:57:06 PM
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.


Title: Re: Incremental 2D Sobol based point set
Post by: Roquen on October 18, 2013, 09:55:05 AM
I know...I'm rather hoping someone else will use it and provide a demonstration.  ;)


Title: Re: Incremental 2D Sobol based point set
Post by: Alef on October 26, 2013, 03:52:03 PM
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.



Title: Re: Incremental 2D Sobol based point set
Post by: Roquen on March 28, 2014, 10:15:39 PM
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