gistfile1.java

class vec
{
  float x, y, z; 

  vec(float x, float y, float z) {
    this.x = x;
    this.y = y;
    this.z = z;
  }

  vec(vec v) {
    this.x = v.x;
    this.y = v.y;
    this.z = v.z;
  }

  vec add(vec b) {
    return new vec(this.x + b.x, this.y + b.y, this.z + b.z);
  }

  vec sub(vec b) {
    return new vec(this.x - b.x, this.y - b.y, this.z - b.z);
  }

  vec cross(vec b) {
    float u = this.y * b.z - b.y * this.z;
    float v = this.z * b.x - b.z * this.x;
    float w = this.x * b.y - b.x * this.y;

    return new vec(u, v, w);
  }

  void normalize() {
    float d = this.len();
    float invlen;

    if (abs(d) > 1.0e-6) {
      invlen = 1.0 / d;
      this.x *= invlen;
      this.y *= invlen;
      this.z *= invlen;
    }  
  }

  float len() {
    float d = x * x + y * y + z * z;
    return sqrt(d);
  }

  float dot(vec b) {
    return this.x * b.x + this.y * b.y + this.z * b.z;
  }

  vec neg() {
    return new vec(-x, -y, -z);
  }

}

class Ray
{
  vec org;
  vec dir;

  Ray(vec o, vec d) {
    org = o;
    dir = d;
  }

}

class Intersection
{
  float    t;
  vec      p;     // hit point
  vec      n;     // normal
  boolean  hit;

  Intersection() {
    hit = false;
    t   = 1.0e+30;
    n   = new vec(0.0, 0.0, 0.0);
  }

}

class Sphere
{
  vec   center;
  float radius; 

  Sphere(vec center, float radius) {
    this.center = center;
    this.radius = radius;
  }

  void intersect(Intersection isect, Ray ray)
  {

    // rs = ray.org - sphere.center
    vec rs  = ray.org.sub(this.center);
    float B = rs.dot(ray.dir);
    float C = rs.dot(rs) - (this.radius * this.radius);
    float D = B * B - C;

    if (D > 0.0) {
      float t = -B - sqrt(D);
      if ( (t > 0.0) && (t < isect.t) ) {
        isect.t = t;
        isect.hit = true;

        // calculate normal.
        vec p = new vec(ray.org.x + ray.dir.x * t,
                        ray.org.y + ray.dir.y * t,
                        ray.org.z + ray.dir.z * t);
        vec n = p.sub(center);
        n.normalize();
        isect.n = n;

        isect.p = p;
      }
    }

  }  
}

class Plane
{
  vec p;  // point on the plane
  vec n;  // normal to the plane
  
  Plane(vec p, vec n) {
    this.p = p;
    this.n = n;
  }
  
  void intersect(Intersection isect, Ray ray) {
    // d = -(p . n)
    // t = -(ray.org . n + d) / (ray.dir . n)
    float d = -p.dot(n);
    float v = ray.dir.dot(n);

    if (abs(v) < 1.0e-6) return; // the plane is parallel to the ray.

    float t = -(ray.org.dot(n) + d) / v;
    
    if ( (t > 0) && (t < isect.t) ) {
      isect.hit = true;
      isect.t   = t;
      isect.n   = n;

      vec p = new vec(ray.org.x + t * ray.dir.x,
                      ray.org.y + t * ray.dir.y,
                      ray.org.z + t * ray.dir.z);

      isect.p   = p;
    }
  }
  
}

class Triangle
{
  vec v0, v1, v2;

  Triangle(vec v0, vec v1, vec v2) {
    this.v0 = v0;
    this.v1 = v1;
    this.v2 = v2;
  }

  void intersect(Intersection isect, Ray ray) {

    vec e1 = v1.sub(v0);
    vec e2 = v2.sub(v0);

    vec p  = ray.dir.cross(e2);

    float det = e1.dot(p);
    float invdet;

    if (abs(det) > 1.0e-6) {
      invdet = 1.0 / det;
    } 
    else {
      return;
    }

    vec s = ray.org.sub(v0);
    vec q = s.cross(e1);

    float u = s.dot(p) * invdet;
    float v = q.dot(ray.dir) * invdet;
    float t = e2.dot(q) * invdet;

    if ( (u < 0.0) || (u > 1.0)     ) return;
    if ( (v < 0.0) || ((u+v) > 1.0) ) return;
    if ( (t < 0.0) || (t > isect.t) ) return;

    isect.hit = true;
    isect.t   = t;

  } 

}

PFont      font;
boolean    rendererd = false;
int        renderY   = 0;
int        renderLine[];
float      startTime, endTime, elapsedTime;

// render settings
final int  IMAGE_WIDTH   = 256;
final int  IMAGE_HEIGHT  = 256;
final int  NSUBSAMPLES   = 2;
final int  NAO_SAMPLES   = 8;     // # of sample rays for ambient occlusion

// scene data
Sphere  spheres[];
Plane   plane;

void setup() {

  size(IMAGE_WIDTH, IMAGE_HEIGHT);
  font = loadFont("Monospaced-16.vlw");
  textFont(font,16);

  renderLine = new int[width];

  spheres    = new Sphere[3];
  spheres[0] = new Sphere(new vec(-2.0, 0.0, -3.5), 0.5);
  spheres[1] = new Sphere(new vec(-0.5, 0.0, -3.0), 0.5);
  spheres[2] = new Sphere(new vec(1.0, 0.0, -2.2), 0.5);
  plane  = new Plane(new vec(0.0, -0.5, 0.0), new vec(0.0, 1.0, 0.0));

  startTime = millis();
}

void draw() {
  int i;

  if (rendererd) {
    String s = str(elapsedTime / 1000.0) + "sec";
    text(s, 10, 20);
    return;
  }

  render(renderLine, width, height, renderY, NSUBSAMPLES);
  
  loadPixels();
  
  for (i = 0; i < width; i++) {
    pixels[renderY * width + i] = renderLine[i];
  }
  
  renderY = renderY + 1;
  if (renderY >= height) {
    rendererd   = true;
    endTime     = millis();
    elapsedTime = endTime - startTime;

  }
  
  updatePixels();
}

char clamp(float f)
{
  int i = (int)(f * 255.5);

  if (i < 0) i = 0;
  if (i > 255) i = 255;

  return (char)i;  

}

void orthoBasis(vec basis[], vec n)
{
  int i;

  basis[2] = new vec(n);
  basis[1] = new vec(0.0, 0.0, 0.0);

  if ((n.x < 0.6) && (n.x > -0.6)) {
    basis[1].x = 1.0;
  } else if ((n.y < 0.6) && (n.y > -0.6)) {
    basis[1].y = 1.0;
  } else if ((n.z < 0.6) && (n.z > -0.6)) {
    basis[1].z = 1.0;
  } else {
    basis[1].x = 1.0;
  }

  basis[0] = basis[1].cross(basis[2]);
  basis[0].normalize();

  basis[1] = basis[2].cross(basis[0]);
  basis[1].normalize();

}

vec ambientOcclusion(Intersection isect)
{
  int   i, j;
  int   ntheta = NAO_SAMPLES;
  int   nphi   = NAO_SAMPLES;
  float eps = 0.0001;

  // Slightly move ray org towards ray dir to avoid numerical probrem.
  vec p = new vec(isect.p.x + eps * isect.n.x,
                  isect.p.y + eps * isect.n.y,
                  isect.p.z + eps * isect.n.z);

  // Calculate orthogonal basis.
  vec basis[]; basis = new vec[3];
  orthoBasis(basis, isect.n);

  float occlusion = 0.0;

  for (j = 0; j < ntheta; j++) {
    for (i = 0; i < nphi; i++) {

      // Pick a random ray direction with importance sampling.
      // p = cos(theta) / PI
      float r     = random(1.0);
      float phi   = 2.0 * PI * random(1.0);

      float x = cos(phi) * sqrt(1.0 - r);
      float y = sin(phi) * sqrt(1.0 - r);
      float z = sqrt(r);

      // local -> global
      float rx = x * basis[0].x + y * basis[1].x + z * basis[2].x;
      float ry = x * basis[0].y + y * basis[1].y + z * basis[2].y;
      float rz = x * basis[0].z + y * basis[1].z + z * basis[2].z;

      vec raydir = new vec(rx, ry, rz);

      Ray ray = new Ray(p, raydir);

      Intersection occIsect = new Intersection();
      spheres[0].intersect(occIsect, ray);
      spheres[1].intersect(occIsect, ray);
      spheres[2].intersect(occIsect, ray);
      plane.intersect(occIsect, ray);

      if (occIsect.hit) occlusion += 1.0;

    }
  }

  // [0.0, 1.0]
  occlusion = (ntheta * nphi - occlusion) / (float)(ntheta * nphi);
  
  return new vec(occlusion, occlusion, occlusion);
}

void render(int[] img, int width, int height, int y, int nsubsamples)
{
  int i, j;
  int u, v;

  float[] fimg = new float[3 * img.length];

  // clear floating point image
  for (i = 0; i < img.length; i++) {
    fimg[i] = 0.0;
  }


  for (i = 0; i < width; i++ ) {

    // subsampling
    for (v = 0; v < nsubsamples; v++) {
      for (u = 0; u < nsubsamples; u++) {
        float px = (i + (u / (float)nsubsamples) - (width / 2.0))/(width / 2.0);
        float py = (y + (v / (float)nsubsamples) - (height / 2.0))/(height / 2.0);
        py = -py;   // flip Y

        float t = 10000.0;                
        vec eye = new vec(px, py, -1.0);
        eye.normalize();

        Ray ray = new Ray(new vec(0.0, 0.0, 0.0), new vec(eye));

        Intersection isect = new Intersection();
        spheres[0].intersect(isect, ray);
        spheres[1].intersect(isect, ray);
        spheres[2].intersect(isect, ray);
        plane.intersect(isect, ray);
          
        if (isect.hit) {

          t = isect.t;

          vec col = ambientOcclusion(isect);

          fimg[3 * i + 0] += col.x;
          fimg[3 * i + 1] += col.y;
          fimg[3 * i + 2] += col.z;
        }        

      }
    }  

    fimg[3 * i + 0] /= (nsubsamples * nsubsamples);
    fimg[3 * i + 1] /= (nsubsamples * nsubsamples);
    fimg[3 * i + 2] /= (nsubsamples * nsubsamples);

    float r = fimg[3 * i + 0];
    float g = fimg[3 * i + 1];
    float b = fimg[3 * i + 2];

    img[i] = color(clamp(r), clamp(g), clamp(b));
    //img[i] = color(255, 255, 255);

  }

}