James Bigler: CS6650 Homework 2

Using this function L(x,y) = 0.5 * (1 + sin(0.01 * x*x + 0.01 * y*y)) plotted over the domain of x:[0,255], y:[0,255] the artifacts of aliasing is immediately noticed. The image can be antialiased by using a filter. The distribution of points can have a big impact on the quality of the sampling. Using only 64 samples per pixels I examine the absolute error and average relative error of different distribution of points as compared to an image generated by a cubic bspline filter with 100,000 samples per pixel.

The sampled image is on the left and the difference image is on the right. Error is reported on the image space of [0,1].

Regular

Random

Jitter

Multi-jitter

Hammersley base 2

Hammersley base 5

Hammersley base 7

// this is sample code

// the function that evalulates L
static float L(float x, float y) {
  return (0.5 * (1 + sinf(0.01*x*x + 0.01*y*y))); 
}

// genfrand is a random number [0,1)

RGB RandomSample::render(const int x, const int y, Worker * wrk) {
  float val = 0;
  for (int index = 0; index < 64; index++) {
    float i = -0.5 + wrk->get_rng()->genfrand();
    float j = -0.5 + wrk->get_rng()->genfrand();
    val += L(x + i, y + j);
  }
  val /= 64;
  return RGB(val,val,val);
}

RGB RegularSample::render(const int x, const int y, Worker * /*wrk*/) {
  float val = 0;
  for (float i = -0.5; i < 0.5; i+=0.125)
    for (float j = -0.5; j < 0.5; j+=0.125) {
      val += L(x + i, y + j);
    }
  val /= 64;
  return RGB(val,val,val);
}

RGB JitterSample::render(const int x, const int y, Worker * wrk) {
  float val = 0;
  for (float i = -0.5; i < 0.5; i+=0.125)
    for (float j = -0.5; j < 0.5; j+=0.125) {
      // [-0.5,0.5] * 0.125 = [-0.0625, 0.0625]
      float i_jitter = (-0.5 + wrk->get_rng()->genfrand()) * 0.125;
      float j_jitter = (-0.5 + wrk->get_rng()->genfrand()) * 0.125;
      val += L(x + i + i_jitter, y + j + j_jitter);
    }
  val /= 64;
  return RGB(val,val,val);
}

void MultiJitterSample::fill_points(Point *p, int m, int n, MT_RNG *rng) {
  // code adapted from Graphics Gems IV pg.370-374 by Chui et.al.
  int i,j;
  float subcell_width = 1./(m*n);
  // Initialize points to the "canonical" multi-jittered pattern
  for (i = 0; i < m; i++)
    for (j = 0; j < n; j++) {
      p[i*n + j].x = (i*n + j + rng->genfrand()) * subcell_width;
      p[i*n + j].y = (j*m + i + rng->genfrand()) * subcell_width;
    }
  // Shuffle x coordinates within each column of cells
  for (i = 0; i < m; i++)
    for (j = 0; j < n; j++) {
      int k = (int) (rng->genfrand() * (n-1-j) + j);
      float t = p[i*n + j].x;
      p[i*n + j].x = p[i*n + k].x;
      p[i*n + k].x = t;
    }  
  for (i = 0; i < m; i++)
    for (j = 0; j < n; j++) {
      int k = (int) (rng->genfrand() * (m-1-j) + j);
      float t = p[j*n + i].y;
      p[j*n + i].y = p[k*n + i].y;
      p[k*n + i].y = t;
    }  
}

RGB MultiJitterSample::render(const int x, const int y, Worker * wrk) {
#define MULTI_M 8
#define MULTI_N 8
  Point p[MULTI_M * MULTI_N];
  fill_points(p,MULTI_M, MULTI_N, wrk->get_rng());
  float val = 0;
  for (int i = 0; i < MULTI_M; i++)
    for (int j = 0; j < MULTI_N; j++) {
      val += L(x + p[i*MULTI_N + j].x - 0.5, y + p[i*MULTI_N + j].y - 0.5);
    }
  val /= 64;
  return RGB(val,val,val);
}

/* Based on code taken from 
 * Copyright (c) Tien-Tsin Wong, 1997.  
 * All Rights Reserved.
 * 
 * 11 Feb 1997
 *
 * ip = 1.0/base;
 */
RGB HammersleySample::render(const int x, const int y, Worker * wrk) {
  float val = 0;

  float a, p;
  int k, kk;
  
  for (k=0; k< iters ; k++)
  {
    float u = 0;
    for (p=ip, kk=k ; kk ; p*=ip, kk/=base)  // kk = (int)(kk/base)
      if ((a = kk % base))
	u += a * p;

    float v = (k + 0.5) / iters;

    val += L(x + u - 0.5, y + v - 0.5);
  }
  
  val /= iters;
  return RGB(val,val,val);
}