immutable.gen.Vec3d.scala

package jp.peppermint.vec.immutable.gen

object Real {
  trait DoubleIsReal extends Real[Double] {
    def sqrt(x: Double): Double = math.sqrt(x)
  }
  implicit object DoubleIsReal extends DoubleIsReal
}

trait Real[T] {
  //def fromDouble(x: Double): T
  
  def sqrt(x: T): T
  
  class Ops(lhs: T) {
    def sqrt(): T = Real.this.sqrt(lhs)
  }
  implicit def mkRealOps(lhs: T): Ops = new Ops(lhs)
}


case class Vec3[T](x: T, y: T, z: T)(implicit num: Fractional[T], real: Real[T]) {
  import num._
  import real._

  def +(v: Vec3[T]) = Vec3[T](x + v.x, y + v.y, z + v.z)
  def -(v: Vec3[T]) = Vec3[T](x - v.x, y - v.y, z - v.z)
  def *(s: T)=Vec3[T](x * s, y * s, z * s)
  def *(v: Vec3[T]) = x * v.x + y * v.y + z * v.z
  def **(v: Vec3[T]) = Vec3[T](y * v.z - z * v.y,
			       z * v.x - x * v.z,
			       x * v.y - y * v.x)


  def unary_- = Vec3[T](-x, -y, -z)

  def normalize: Vec3[T] = {
    val length:T = (this * this).sqrt

    if (length.abs > zero) {
      return Vec3[T](x / length, y / length, z / length)
    } else {
      return this
    }
  }
}

class Ao {

  val WIDTH       = 256
  val HEIGHT      = 256
  val NSUBSAMPLES = 2
  val NAO_SAMPLES = 8

  class Isect {
    var t: Double = 0
    var p: Vec3[Double] = null
    var n: Vec3[Double] = null
    var hit: Int = 0
  }
  
  class Sphere {
    var center: Vec3[Double] = null
    var radius: Double = 0
  }

  class Plane {
    var p: Vec3[Double] = null
    var n: Vec3[Double] = null
  }

  class Ray {
    var org: Vec3[Double] = null
    var dir: Vec3[Double] = null
  }

  var spheres: Array[Sphere] = null
  var plane: Plane = null


  def ray_sphere_intersect(isect: Isect, ray: Ray, sphere: Sphere) {
    val rs = ray.org - sphere.center

    val B = rs * ray.dir
    val C = rs * rs - sphere.radius * sphere.radius
    val D = B * B - C

    if (D > 0.0) {
      val t = -B - math.sqrt(D);
        
      if ((t > 0.0) && (t < isect.t)) {
        isect.t = t;
        isect.hit = 1;
            
        isect.p = ray.org + ray.dir * t
        isect.n = (isect.p - sphere.center).normalize
      }
    }
  }

  def ray_plane_intersect(isect: Isect, ray: Ray, plane: Plane) {
    val d = - plane.p * plane.n
    val v = ray.dir * plane.n

    if (v.abs < 1.0e-17) return;

    val t = -(ray.org * plane.n + d) / v;

    if ((t > 0.0) && (t < isect.t)) {
        isect.t = t;
        isect.hit = 1;
        
        isect.p = ray.org + ray.dir * t;

        isect.n = plane.n;
    }
  }

  def orthoBasis(basis: Array[Vec3[Double]], n: Vec3[Double]) {
    basis(2) = n;

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

    basis(0) = (basis(1) ** basis(2)).normalize

    basis(1) = (basis(2) ** basis(0)).normalize
  }


  def ambient_occlusion(col: Array[Double], isect: Isect) {
    val ntheta = NAO_SAMPLES
    val nphi   = NAO_SAMPLES
    val eps = 0.0001

    val p = isect.p + isect.n * eps

    val basis = new Array[Vec3[Double]](3)
    orthoBasis(basis, isect.n);

    var occlusion = 0.0;

    for (j <- 0 until ntheta) {
      for (i <- 0 until nphi) {
        val theta = math.sqrt(math.random)
        val phi   = 2.0 * math.Pi * math.random

        val x = math.cos(phi) * theta;
        val y = math.sin(phi) * theta;
        val z = math.sqrt(1.0 - theta * theta);

        // local -> global
        val ray= new Ray
	
        ray.org = p;
        ray.dir = basis(0) * x + basis(1) * y + basis(2) * z

        val occIsect = new Isect
        occIsect.t   = 1.0e+17
        occIsect.hit = 0

        ray_sphere_intersect(occIsect, ray, spheres(0))
        ray_sphere_intersect(occIsect, ray, spheres(1))
        ray_sphere_intersect(occIsect, ray, spheres(2))
        ray_plane_intersect (occIsect, ray, plane); 

        if (occIsect.hit > 0)
	  occlusion += 1.0;
            
      }
    }

    occlusion = (ntheta * nphi - occlusion) / (ntheta * nphi)

    col(0) = occlusion;
    col(1) = occlusion;
    col(2) = occlusion;
  }

  def clamp(f: Double): Byte = {
    val i = (f * 255.5).toInt

    if (i < 0)
      return 0.toByte
    if (i > 255)
      return 255.toByte

    return i.toByte
  }


  def render(img: Array[Byte], w: Int, h: Int, nsubsamples: Int) {
    val fimg = new Array[Double](w * h * 3);

    for (y <- 0 until h) {
      for (x <- 0 until w) {
            
        for (v <- 0 until nsubsamples) {
          for (u <- 0 until nsubsamples) {
            val px = (x + (u / nsubsamples.toDouble) - (w / 2.0)) / (w / 2.0);
            val py = -(y + (v / nsubsamples.toDouble) - (h / 2.0)) / (h / 2.0);

            val ray = new Ray

            ray.org = Vec3[Double](0.0, 0.0, 0.0)
            ray.dir = Vec3[Double](px, py, -1.0).normalize

            val isect = new Isect
            isect.t   = 1.0e+17;
            isect.hit = 0;

            ray_sphere_intersect(isect, ray, spheres(0));
            ray_sphere_intersect(isect, ray, spheres(1));
            ray_sphere_intersect(isect, ray, spheres(2));
            ray_plane_intersect (isect, ray, plane);

            if (isect.hit > 0) {
              val col =new Array[Double](3);
              ambient_occlusion(col, isect);

              fimg(3 * (y * w + x) + 0) += col(0);
              fimg(3 * (y * w + x) + 1) += col(1);
              fimg(3 * (y * w + x) + 2) += col(2);
            }
	    
          }
        }
	
        fimg(3 * (y * w + x) + 0) /= nsubsamples * nsubsamples;
        fimg(3 * (y * w + x) + 1) /= nsubsamples * nsubsamples;
        fimg(3 * (y * w + x) + 2) /= nsubsamples * nsubsamples;
        
        img(3 * (y * w + x) + 0) = clamp(fimg(3 *(y * w + x) + 0));
        img(3 * (y * w + x) + 1) = clamp(fimg(3 *(y * w + x) + 1));
        img(3 * (y * w + x) + 2) = clamp(fimg(3 *(y * w + x) + 2));
      }
    }
    
  }

  def init_scene() {
    spheres = new Array[Sphere](3)
    spheres(0) = new Sphere
    spheres(0).center = Vec3[Double](-2.0, 0.0, -3.5)
    spheres(0).radius = 0.5;
    
    spheres(1) = new Sphere
    spheres(1).center = Vec3[Double](-0.5, 0.0, -3.0)
    spheres(1).radius = 0.5;
    
    spheres(2) = new Sphere
    spheres(2).center = Vec3[Double](1.0, 0.0, -2.2)
    spheres(2).radius = 0.5;

    
    plane = new Plane
    plane.p = Vec3[Double](0.0, -0.5, 0.0)

    plane.n = Vec3[Double](0.0, 1.0, 0.0)
  }

  def saveppm(fname: String, w: Int, h: Int, img: Array[Byte]) {
    var out = new java.io.DataOutputStream(new java.io.FileOutputStream(fname))
    out.writeBytes("P6\n")
    out.writeBytes("%d %d\n".format(w, h))
    out.writeBytes("255\n")
    out.write(img, 0, w * h * 3);
  }

  def exec(save: Boolean) {
    val img = new Array[Byte](WIDTH * HEIGHT * 3);

    init_scene();

    render(img, WIDTH, HEIGHT, NSUBSAMPLES);

    if (save)
      saveppm("ao.ppm", WIDTH, HEIGHT, img); 
  }
}

object App {

  def benchmark(times: Int, mult: Int = 1)(f: =>Any) = 
    new scala.testing.Benchmark{ multiplier = mult; def run() = f }.runBenchmark(times)

  implicit def benchmark2report(list:List[Long]) =
    new { def report() = {
      println(list);
      println((list.sum - list.max - list.min).toDouble / (list.size-2))
    }}


  def main(argv: Array[String]) {
    val ao = new Ao
    ao.exec(true)

    benchmark(5){ (new Ao).exec(false) }.report
  }
}

immutable.Vec3d.scala

package jp.peppermint.vec.immutable

case class Vec3d(x:Double, y:Double, z: Double) {
  def +(v: Vec3d) = Vec3d(x + v.x, y + v.y, z + v.z)
  def -(v: Vec3d) = Vec3d(x - v.x, y - v.y, z - v.z)
  def *(s: Double) = Vec3d(x * s, y * s, z * s)
  def *(v: Vec3d) = x * v.x + y * v.y + z * v.z
  def **(v: Vec3d) = Vec3d(y * v.z - z * v.y,
			   z * v.x - x * v.z,
			   x * v.y - y * v.x)

  def unary_- = Vec3d(-x, -y, -z)

  def normalize: Vec3d = {
    val length = math.sqrt((this * this).toDouble)

    if (math.abs(length) > 1.0e-17) {
      return Vec3d(x / length, y / length, z / length)
    } else {
      return this
    }
  }
}

class Ao {

  val WIDTH       = 256
  val HEIGHT      = 256
  val NSUBSAMPLES = 2
  val NAO_SAMPLES = 8

  class Isect {
    var t: Double = 0
    var p: Vec3d = null
    var n: Vec3d = null
    var hit: Int = 0
  }
  
  class Sphere {
    var center: Vec3d = null
    var radius: Double = 0
  }

  class Plane {
    var p: Vec3d = null
    var n: Vec3d = null
  }

  class Ray {
    var org: Vec3d = null
    var dir: Vec3d = null
  }

  var spheres: Array[Sphere] = null
  var plane: Plane = null


  def ray_sphere_intersect(isect: Isect, ray: Ray, sphere: Sphere) {
    val rs = ray.org - sphere.center

    val B = rs * ray.dir
    val C = rs * rs - sphere.radius * sphere.radius
    val D = B * B - C

    if (D > 0.0) {
      val t = -B - math.sqrt(D);
        
      if ((t > 0.0) && (t < isect.t)) {
        isect.t = t;
        isect.hit = 1;
            
        isect.p = ray.org + ray.dir * t
        isect.n = (isect.p - sphere.center).normalize
      }
    }
  }

  def ray_plane_intersect(isect: Isect, ray: Ray, plane: Plane) {
    val d = - plane.p * plane.n
    val v = ray.dir * plane.n

    if (math.abs(v) < 1.0e-17) return;

    val t = -(ray.org * plane.n + d) / v;

    if ((t > 0.0) && (t < isect.t)) {
        isect.t = t;
        isect.hit = 1;
        
        isect.p = ray.org + ray.dir * t;

        isect.n = plane.n;
    }
  }

  def orthoBasis(basis: Array[Vec3d], n: Vec3d) {
    basis(2) = n;

    basis(1) =
      if ((n.x < 0.6) && (n.x > -0.6)) {
	Vec3d(1.0, 0, 0)
      } else if ((n.y < 0.6) && (n.y > -0.6)) {
	Vec3d(0, 1.0, 0)
      } else if ((n.z < 0.6) && (n.z > -0.6)) {
        Vec3d(0, 0, 1.0)
      } else {
	Vec3d(1.0, 0, 0)
      }

    basis(0) = (basis(1) ** basis(2)).normalize

    basis(1) = (basis(2) ** basis(0)).normalize
  }


  def ambient_occlusion(col: Array[Double], isect: Isect) {
    val ntheta = NAO_SAMPLES
    val nphi   = NAO_SAMPLES
    val eps = 0.0001

    val p = isect.p + isect.n * eps

    val basis = new Array[Vec3d](3)
    orthoBasis(basis, isect.n);

    var occlusion = 0.0;

    for (j <- 0 until ntheta) {
      for (i <- 0 until nphi) {
        val theta = math.sqrt(math.random)
        val phi   = 2.0 * math.Pi * math.random

        val x = math.cos(phi) * theta;
        val y = math.sin(phi) * theta;
        val z = math.sqrt(1.0 - theta * theta);

        // local -> global
        val ray= new Ray
	
        ray.org = p;
        ray.dir = basis(0) * x + basis(1) * y + basis(2) * z

        val occIsect = new Isect
        occIsect.t   = 1.0e+17
        occIsect.hit = 0

        ray_sphere_intersect(occIsect, ray, spheres(0))
        ray_sphere_intersect(occIsect, ray, spheres(1))
        ray_sphere_intersect(occIsect, ray, spheres(2))
        ray_plane_intersect (occIsect, ray, plane); 

        if (occIsect.hit > 0)
	  occlusion += 1.0;
            
      }
    }

    occlusion = (ntheta * nphi - occlusion) / (ntheta * nphi)

    col(0) = occlusion;
    col(1) = occlusion;
    col(2) = occlusion;
  }

  def clamp(f: Double): Byte = {
    val i = (f * 255.5).toInt

    if (i < 0)
      return 0.toByte
    if (i > 255)
      return 255.toByte

    return i.toByte
  }


  def render(img: Array[Byte], w: Int, h: Int, nsubsamples: Int) {
    val fimg = new Array[Double](w * h * 3);

    for (y <- 0 until h) {
      for (x <- 0 until w) {
            
        for (v <- 0 until nsubsamples) {
          for (u <- 0 until nsubsamples) {
            val px = (x + (u / nsubsamples.toDouble) - (w / 2.0)) / (w / 2.0);
            val py = -(y + (v / nsubsamples.toDouble) - (h / 2.0)) / (h / 2.0);

            val ray = new Ray

            ray.org = Vec3d(0.0, 0.0, 0.0)
            ray.dir = Vec3d(px, py, -1.0).normalize

            val isect = new Isect
            isect.t   = 1.0e+17;
            isect.hit = 0;

            ray_sphere_intersect(isect, ray, spheres(0));
            ray_sphere_intersect(isect, ray, spheres(1));
            ray_sphere_intersect(isect, ray, spheres(2));
            ray_plane_intersect (isect, ray, plane);

            if (isect.hit > 0) {
              val col =new Array[Double](3);
              ambient_occlusion(col, isect);

              fimg(3 * (y * w + x) + 0) += col(0);
              fimg(3 * (y * w + x) + 1) += col(1);
              fimg(3 * (y * w + x) + 2) += col(2);
            }
	    
          }
        }
	
        fimg(3 * (y * w + x) + 0) /= nsubsamples * nsubsamples;
        fimg(3 * (y * w + x) + 1) /= nsubsamples * nsubsamples;
        fimg(3 * (y * w + x) + 2) /= nsubsamples * nsubsamples;
        
        img(3 * (y * w + x) + 0) = clamp(fimg(3 *(y * w + x) + 0));
        img(3 * (y * w + x) + 1) = clamp(fimg(3 *(y * w + x) + 1));
        img(3 * (y * w + x) + 2) = clamp(fimg(3 *(y * w + x) + 2));
      }
    }
    
  }

  def init_scene() {
    spheres = new Array[Sphere](3)
    spheres(0) = new Sphere
    spheres(0).center = Vec3d(-2.0, 0.0, -3.5)
    spheres(0).radius = 0.5;
    
    spheres(1) = new Sphere
    spheres(1).center = Vec3d(-0.5, 0.0, -3.0)
    spheres(1).radius = 0.5;
    
    spheres(2) = new Sphere
    spheres(2).center = Vec3d(1.0, 0.0, -2.2)
    spheres(2).radius = 0.5;

    
    plane = new Plane
    plane.p = Vec3d(0.0, -0.5, 0.0)

    plane.n = Vec3d(0.0, 1.0, 0.0)
  }

  def saveppm(fname: String, w: Int, h: Int, img: Array[Byte]) {
    var out = new java.io.DataOutputStream(new java.io.FileOutputStream(fname))
    out.writeBytes("P6\n")
    out.writeBytes("%d %d\n".format(w, h))
    out.writeBytes("255\n")
    out.write(img, 0, w * h * 3);
  }

  def exec(save: Boolean) {
    val img = new Array[Byte](WIDTH * HEIGHT * 3);

    init_scene();

    render(img, WIDTH, HEIGHT, NSUBSAMPLES);

    if (save)
      saveppm("ao.ppm", WIDTH, HEIGHT, img); 
  }
}

object App {

  def benchmark(times: Int, mult: Int = 1)(f: =>Any) = 
    new scala.testing.Benchmark{ multiplier = mult; def run() = f }.runBenchmark(times)

  implicit def benchmark2report(list:List[Long]) =
    new { def report() = {
      println(list);
      println((list.sum - list.max - list.min).toDouble / (list.size-2))
    }}

  def main(argv: Array[String]) {
    val ao = new Ao
    ao.exec(true)

    benchmark(5){ (new Ao).exec(false) }.report
  }
}

mutable.Vec3d.scala

package jp.peppermint.vec.mutable

class Vec3d(var x:Double = 0, var y:Double = 0, var z: Double = 0) {
  def this(v: Vec3d) =
    this(v.x, v.y, v.z)

  def +=(v: Vec3d) {
    x += v.x
    y += v.y
    z += v.z
  }
  def -=(v: Vec3d) {
    x -= v.x
    y -= v.y
    z -= v.z
  }
  def *=(s: Double) {
    x *= s
    y *= s
    z *= s
  }
  def *(v: Vec3d) = x * v.x + y * v.y + z * v.z
  def **=(v: Vec3d) {
    val xx = y * v.z - z * v.y
    val yy = z * v.x - x * v.z
    val zz = x * v.y - y * v.x
    x = xx
    y = yy
    z = zz
  }

  def set(v: Vec3d) {
    x = v.x
    y = v.y
    z = v.z
  }

  def set(xx: Double, yy: Double, zz: Double) {
    x = xx
    y = yy
    z = zz
  }

  def normalize {
    val length = math.sqrt((this * this).toDouble)

    if (math.abs(length) > 1.0e-17) {
      x /= length
      y /= length
      z /= length
    }
  }
}

class Ao {

  val WIDTH       = 256
  val HEIGHT      = 256
  val NSUBSAMPLES = 2
  val NAO_SAMPLES = 8

  class Isect {
    var t: Double = 0
    val p: Vec3d = new Vec3d
    val n: Vec3d = new Vec3d
    var hit: Int = 0
  }
  
  class Sphere {
    var center: Vec3d = null
    var radius: Double = 0
  }

  class Plane {
    var p: Vec3d = null
    var n: Vec3d = null
  }

  class Ray {
    val org: Vec3d = new Vec3d
    val dir: Vec3d = new Vec3d
  }

  var spheres: Array[Sphere] = null
  var plane: Plane = null


  def ray_sphere_intersect(isect: Isect, ray: Ray, sphere: Sphere) {
    val rs = new Vec3d(ray.org)
    rs -= sphere.center

    val B = rs * ray.dir
    val C = rs * rs - sphere.radius * sphere.radius
    val D = B * B - C

    if (D > 0.0) {
      val t = -B - math.sqrt(D);
        
      if ((t > 0.0) && (t < isect.t)) {
        isect.t = t;
        isect.hit = 1;

        //isect.p = ray.org + ray.dir * t
        isect.p.set(ray.dir)
	isect.p *= t
	isect.p += ray.org
        //isect.n = (isect.p - sphere.center).normalize
	isect.n set isect.p
	isect.n -= sphere.center
	isect.n.normalize
      }
    }
  }

  def ray_plane_intersect(isect: Isect, ray: Ray, plane: Plane) {
    val d = - (plane.p * plane.n)
    val v = ray.dir * plane.n

    if (math.abs(v) < 1.0e-17) return;

    val t = -(ray.org * plane.n + d) / v;

    if ((t > 0.0) && (t < isect.t)) {
        isect.t = t;
        isect.hit = 1;
        
	//isect.p = ray.org + ray.dir * t;
        isect.p set ray.dir
	isect.p *= t
	isect.p += ray.org

        isect.n set plane.n;
    }
  }

  def orthoBasis(basis: Array[Vec3d], n: Vec3d) {
    basis(2) set n;

    if ((n.x < 0.6) && (n.x > -0.6)) {
      basis(1).set(1.0, 0, 0)
    } else if ((n.y < 0.6) && (n.y > -0.6)) {
      basis(1).set(0, 1.0, 0)
    } else if ((n.z < 0.6) && (n.z > -0.6)) {
      basis(1).set(0, 0, 1.0)
    } else {
      basis(1).set(1.0, 0, 0)
    }

    //basis(0) = (basis(1) ** basis(2)).normalize
    basis(0) set basis(1)
    basis(0) **= basis(2)
    basis(0).normalize

    //basis(1) = (basis(2) ** basis(0)).normalize
    basis(1) set basis(2)
    basis(1) **= basis(0)
    basis(1).normalize
  }

  def ambient_occlusion(col: Array[Double], isect: Isect) {
    val ntheta = NAO_SAMPLES
    val nphi   = NAO_SAMPLES
    val eps = 0.0001

    //val p = isect.p + isect.n * eps
    val p = new Vec3d(isect.n)
    p *= eps
    p += isect.p

    val basis = Array(new Vec3d, new Vec3d, new Vec3d)
    orthoBasis(basis, isect.n);

    var occlusion = 0.0;

    val ray = new Ray
    val occIsect = new Isect

    for (j <- 0 until ntheta) {
      for (i <- 0 until nphi) {
        val theta = math.sqrt(math.random)
        val phi   = 2.0 * math.Pi * math.random

        val x = math.cos(phi) * theta;
        val y = math.sin(phi) * theta;
        val z = math.sqrt(1.0 - theta * theta);

        // local -> global
        val rx = x * basis(0).x + y * basis(1).x + z * basis(2).x;
        val ry = x * basis(0).y + y * basis(1).y + z * basis(2).y;
        val rz = x * basis(0).z + y * basis(1).z + z * basis(2).z;

        ray.org set p;
        ray.dir set(rx, ry, rz)

        occIsect.t   = 1.0e+17
        occIsect.hit = 0

        ray_sphere_intersect(occIsect, ray, spheres(0))
        ray_sphere_intersect(occIsect, ray, spheres(1))
        ray_sphere_intersect(occIsect, ray, spheres(2))
        ray_plane_intersect (occIsect, ray, plane); 

        if (occIsect.hit > 0)
	  occlusion += 1.0;
      }
    }

    occlusion = (ntheta * nphi - occlusion) / (ntheta * nphi)

    col(0) = occlusion;
    col(1) = occlusion;
    col(2) = occlusion;
  }

  def clamp(f: Double): Byte = {
    val i = (f * 255.5).toInt

    if (i < 0)
      return 0.toByte
    if (i > 255)
      return 255.toByte

    return i.toByte
  }


  def render(img: Array[Byte], w: Int, h: Int, nsubsamples: Int) {
    val fimg = new Array[Double](w * h * 3);

    val ray = new Ray
    val isect = new Isect
    val col =new Array[Double](3);

    for (y <- 0 until h) {
      for (x <- 0 until w) {
            
        for (v <- 0 until nsubsamples) {
          for (u <- 0 until nsubsamples) {
            val px = (x + (u / nsubsamples.toDouble) - (w / 2.0)) / (w / 2.0);
            val py = -(y + (v / nsubsamples.toDouble) - (h / 2.0)) / (h / 2.0);

            ray.org set(0.0, 0.0, 0.0)
            ray.dir set(px, py, -1.0)
	    ray.dir.normalize

            isect.t   = 1.0e+17;
            isect.hit = 0;

            ray_sphere_intersect(isect, ray, spheres(0));
            ray_sphere_intersect(isect, ray, spheres(1));
            ray_sphere_intersect(isect, ray, spheres(2));
            ray_plane_intersect (isect, ray, plane);

            if (isect.hit > 0) {
              col(0) = 0
              col(1) = 0
              col(2) = 0
              ambient_occlusion(col, isect);

              fimg(3 * (y * w + x) + 0) += col(0);
              fimg(3 * (y * w + x) + 1) += col(1);
              fimg(3 * (y * w + x) + 2) += col(2);
            }
	    
          }
        }
	
        fimg(3 * (y * w + x) + 0) /= nsubsamples * nsubsamples;
        fimg(3 * (y * w + x) + 1) /= nsubsamples * nsubsamples;
        fimg(3 * (y * w + x) + 2) /= nsubsamples * nsubsamples;
        
        img(3 * (y * w + x) + 0) = clamp(fimg(3 *(y * w + x) + 0));
        img(3 * (y * w + x) + 1) = clamp(fimg(3 *(y * w + x) + 1));
        img(3 * (y * w + x) + 2) = clamp(fimg(3 *(y * w + x) + 2));
      }
    }
    
  }

  def init_scene() {
    spheres = new Array[Sphere](3)
    spheres(0) = new Sphere
    spheres(0).center = new Vec3d(-2.0, 0.0, -3.5)
    spheres(0).radius = 0.5;
    
    spheres(1) = new Sphere
    spheres(1).center = new Vec3d(-0.5, 0.0, -3.0)
    spheres(1).radius = 0.5;
    
    spheres(2) = new Sphere
    spheres(2).center = new Vec3d(1.0, 0.0, -2.2)
    spheres(2).radius = 0.5;

    
    plane = new Plane
    plane.p = new Vec3d(0.0, -0.5, 0.0)
    plane.n = new Vec3d(0.0, 1.0, 0.0)
  }

  def saveppm(fname: String, w: Int, h: Int, img: Array[Byte]) {
    var out = new java.io.DataOutputStream(new java.io.FileOutputStream(fname))
    out.writeBytes("P6\n")
    out.writeBytes("%d %d\n".format(w, h))
    out.writeBytes("255\n")
    out.write(img, 0, w * h * 3);
  }

  def exec(save: Boolean) {
    val img = new Array[Byte](WIDTH * HEIGHT * 3);

    init_scene();

    render(img, WIDTH, HEIGHT, NSUBSAMPLES);

    if (save)
      saveppm("ao.ppm", WIDTH, HEIGHT, img); 
  }
}

object App {

  def benchmark(times: Int, mult: Int = 1)(f: =>Any) = 
    new scala.testing.Benchmark{ multiplier = mult; def run() = f }.runBenchmark(times)

  implicit def benchmark2report(list:List[Long]) =
    new { def report() = {
      println(list);
      println((list.sum - list.max - list.min).toDouble / (list.size-2))
    }}

  def main(argv: Array[String]) {
    val ao = new Ao
    ao.exec(true)

    benchmark(5){ (new Ao).exec(false) }.report
  }
}