aobench.scm

(define srand (js-call (js-eval "require") "/Users/yhara/proj/aobench/biwas/node_modules/srand/index.js"))
(js-invoke srand "seed" 12345)
(define (random-real) 
  (js-invoke srand "random"))

(define IMGOUT #t)

(define dumpval print)
(define (dumpvec v)
  (receive (x y z) (vec-xyz v)
    (format #t "~s ~s ~s\n" x y z)))

(define IMAGE_WIDTH 32)
;  8 : 30sec
; 16 : 2min
; 32 : 9min19sec (vs 0.056sec)
; 64 : 32min?
(define IMAGE_HEIGHT IMAGE_WIDTH)
(define NSUBSAMPLES 2)
(define NAO_SAMPLES 8)
(define PI (js-eval "Math.PI"))

(define-record-type vec
  (fields (mutable x) (mutable y) (mutable z)))

(define (vec-xyz v)
  (values (vec-x v) (vec-y v) (vec-z v)))

(define (vadd a b)
  (ake-vec
    (+ (vec-x a) (vec-x b))
    (+ (vec-y a) (vec-y b))
    (+ (vec-z a) (vec-z b))))

(define (vsub a b)
  (make-vec
    (- (vec-x a) (vec-x b))
    (- (vec-y a) (vec-y b))
    (- (vec-z a) (vec-z b))))

(define (vcross a b)
  (make-vec
    (- (* (vec-y a) (vec-z b)) (* (vec-z a) (vec-y b)))
    (- (* (vec-z a) (vec-x b)) (* (vec-x a) (vec-z b)))
    (- (* (vec-x a) (vec-y b)) (* (vec-y a) (vec-x b)))))

(define (vdot a b)
  (+ (* (vec-x a) (vec-x b))
     (* (vec-y a) (vec-y b))
     (* (vec-z a) (vec-z b))))

(define (vlength a)
  (sqrt (+ (* (vec-x a) (vec-x a))
           (* (vec-y a) (vec-y a))
           (* (vec-z a) (vec-z a)))))

(define (vnormalize a)
  (let ((len (vlength a))
        (v (make-vec (vec-x a) (vec-y a) (vec-z a))))
    (when (> (abs len) 1.0e-17)
      (vec-x-set! v (/ (vec-x a) len))
      (vec-y-set! v (/ (vec-y a) len))
      (vec-z-set! v (/ (vec-z a) len)))
    v))

(define-record-type sphere
  (fields center radius))

(define (sphere-intersect! @ ray isect)
  (let ((org (ray-org ray))
        (dir (ray-dir ray)))
    (let* ((rs (vsub org (sphere-center @)))
           (B  (vdot rs dir))
           (C  (- (vdot rs rs) (* (sphere-radius @) (sphere-radius @))))
           (D  (- (* B B) C)))
      (when (> D 0.0)
        (let1 t (- (- B) (sqrt D))
          (when (< 0.0 t (isect-t isect))
            (isect-t-set!   isect t)
            (isect-hit-set! isect #t)
            (isect-p-set!   isect (make-vec (+ (vec-x org) (* (vec-x dir) t))
                                            (+ (vec-y org) (* (vec-y dir) t))
                                            (+ (vec-z org) (* (vec-z dir) t))))
            (let1 n (vsub (isect-p isect) (sphere-center @))
              (isect-n-set! isect (vnormalize n)))))))))

(define-record-type plane
  (fields p n))

(define (plane-intersect! @ ray isect)
  (let ((org (ray-org ray))
        (dir (ray-dir ray)))
    (let ((d (- (vdot (plane-p @) (plane-n @))))
          (v (vdot dir (plane-n @))))
      (unless (< (abs v) 1.0e-17)
        (let1 t (/ (- (+ (vdot org (plane-n @)) d)) v)
          (when (< 0.0 t (isect-t isect))
            (isect-hit-set! isect #t)
            (isect-t-set!   isect t)
            (isect-n-set!   isect (plane-n @))
            (isect-p-set!   isect (make-vec (+ (vec-x org) (* (vec-x dir) t))
                                            (+ (vec-y org) (* (vec-y dir) t))
                                            (+ (vec-z org) (* (vec-z dir) t))))))))))

(define-record-type (ray _make-ray ray?)
  (fields org dir))

(define (make-ray org dir)
  (_make-ray org dir))

(define-record-type isect
  (fields (mutable t) (mutable hit) (mutable p) (mutable n))
  (protocol (lambda (p)
              (lambda ()
                (p 1000000.0 #f (make-vec 0.0 0.0 0.0) (make-vec 0.0 0.0 0.0))))))

(define (clamp f)
  (let1 i (* f 255.5)
    (when (> i 255.0) (set! i 255.0))
    (when (< i 0.0)   (set! i 0.0))
    (round i)))

(define (orthoBasis! basis n)
  (vector-set! basis 2 (make-vec (vec-x n) (vec-y n) (vec-z n)))
  (vector-set! basis 1 (make-vec 0.0 0.0 0.0))
  (let1 basis1 (vector-ref basis 1)
    (cond ((< -0.6 (vec-x n) 0.6)
           (vec-x-set! basis1 1.0))
          ((< -0.6 (vec-y n) 0.6)
           (vec-y-set! basis1 1.0))
          ((< -0.6 (vec-z n) 0.6)
           (vec-z-set! basis1 1.0))
          (else
           (vec-x-set! basis1 1.0))))
  (vector-set! basis 0 (vnormalize
                         (vcross (vector-ref basis 1) (vector-ref basis 2))))
  (vector-set! basis 1 (vnormalize
                         (vcross (vector-ref basis 2) (vector-ref basis 0)))))

(define *spheres* #f)
(define *plane* #f)

(define (init-scene!)
  (set! *spheres* (vector (make-sphere (make-vec -2.0 0.0 -3.5) 0.5)
                          (make-sphere (make-vec -0.5 0.0 -3.0) 0.5)
                          (make-sphere (make-vec  1.0 0.0 -2.2) 0.5)))
  (set! *plane* (make-plane (make-vec 0.0 -0.5 0.0)
                            (make-vec 0.0 1.0 0.0))))

(define (ambient-occlusion isect)
  (let ((basis (vector #f #f #f))
        (ntheta NAO_SAMPLES)
        (nphi NAO_SAMPLES)
        (eps 0.0001)
        (occlusion 0.0))
    (orthoBasis! basis (isect-n isect))
    (let1 p (make-vec (+ (vec-x (isect-p isect)) (* eps (vec-x (isect-n isect))))
                      (+ (vec-y (isect-p isect)) (* eps (vec-y (isect-n isect))))
                      (+ (vec-z (isect-p isect)) (* eps (vec-z (isect-n isect)))))
      (dotimes (j nphi)
        (dotimes (i ntheta)
          (let ((r (random-real))
                (phi (* 2.0 PI (random-real))))
            (let ((x (* (cos phi) (sqrt (- 1.0 r))))
                  (y (* (sin phi) (sqrt (- 1.0 r))))
                  (z (sqrt r)))
              (let-values (((b0x b0y b0z) (vec-xyz (vector-ref basis 0)))
                           ((b1x b1y b1z) (vec-xyz (vector-ref basis 1)))
                           ((b2x b2y b2z) (vec-xyz (vector-ref basis 2))))
                (let ((rx (+ (* x b0x) (* y b1x) (* z b2x)))
                      (ry (+ (* x b0y) (* y b1y) (* z b2y)))
                      (rz (+ (* x b0z) (* y b1z) (* z b2z))))
                  (let* ((raydir (make-vec rx ry rz))
                         (ray (make-ray p raydir))
                         (occ-isect (make-isect)))
                    (sphere-intersect! (vector-ref *spheres* 0) ray occ-isect)
                    (sphere-intersect! (vector-ref *spheres* 1) ray occ-isect)
                    (sphere-intersect! (vector-ref *spheres* 2) ray occ-isect)
                    (plane-intersect! *plane* ray occ-isect)
                    (when (isect-hit occ-isect)
                      (set! occlusion (+ occlusion 1.0)))))))))))

    (let1 occ (/ (- (* ntheta nphi) occlusion) (* ntheta nphi))
      (make-vec occ occ occ))))

(define (render! buf w h nsubsamples)
  (let1 cnt 0
    (dotimes (y h)
      (dotimes (x w)
        (let1 rad (make-vec 0.0 0.0 0.0)
          (dotimes (v nsubsamples)
            (dotimes (u nsubsamples)
              (inc! cnt)
              (let ((px    (/ (+ x (- (/ u nsubsamples) (/ w 2.0))) (/ w 2.0)))
                    (py (- (/ (+ y (- (/ v nsubsamples) (/ h 2.0))) (/ h 2.0)))))
                (let* ((eye (vnormalize (make-vec px py -1.0)))
                       (ray (make-ray (make-vec 0.0 0.0 0.0) eye))
                       (isect (make-isect)))
                  (sphere-intersect! (vector-ref *spheres* 0) ray isect)
                  (sphere-intersect! (vector-ref *spheres* 1) ray isect)
                  (sphere-intersect! (vector-ref *spheres* 2) ray isect)
                  (plane-intersect! *plane* ray isect)

                  (when (isect-hit isect)
                    (let1 col (ambient-occlusion isect)
                      (vec-x-set! rad (+ (vec-x rad) (vec-x col)))
                      (vec-y-set! rad (+ (vec-y rad) (vec-y col)))
                      (vec-z-set! rad (+ (vec-z rad) (vec-z col)))))))))

          (let ((r (/ (vec-x rad) (* nsubsamples nsubsamples)))
                (g (/ (vec-y rad) (* nsubsamples nsubsamples)))
                (b (/ (vec-z rad) (* nsubsamples nsubsamples))))
            (vector-push! buf (clamp r) (clamp g) (clamp b))))))))

(define (dump-binary buf)
  (let1 buffer (js-new "Buffer" (vector-length buf))
    (dotimes (i (vector-length buf))
      (js-invoke buffer "writeUInt8" (vector-ref buf i) i))
    (js-invoke (js-eval "process.stdout") "write" buffer)))
 ; (vector-for-each (lambda (i) (display (integer->char i))) buf))

(define (main)
  (format #t "P6\n")
  (format #t "~a ~a\n" IMAGE_WIDTH IMAGE_HEIGHT)
  (format #t "255\n")
  (let ((buf (make-vector 0))
        (elapsed 0)
        (start (current-date)))
    (init-scene!)
    (render! buf IMAGE_WIDTH IMAGE_HEIGHT 1)
    (if IMGOUT
      (dump-binary buf)
      (vector-for-each print buf))))

(main)

aonode.js

var util = require('util'),
    srand = require('srand');
srand.seed(12345);

var IMGOUT = true

function dumpval(x) { util.puts(x); }
function dumpvec(v) { util.puts(""+v.x+" "+v.y+ " "+v.z); }

var IMAGE_WIDTH  = 32
var IMAGE_HEIGHT = IMAGE_WIDTH
var NSUBSAMPLES  = 2
var NAO_SAMPLES  = 8

Math.random = srand.random;

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

}

function vadd(a, b)
{
    return new vec(a.x + b.x, a.y + b.y, a.z + b.z);
}

function vsub(a, b)
{
    return new vec(a.x - b.x, a.y - b.y, a.z - b.z);
}

function vcross(a, b)
{
    return new vec(a.y * b.z - a.z * b.y,
                   a.z * b.x - a.x * b.z,
                   a.x * b.y - a.y * b.x);
}

function vdot(a, b)
{
    return (a.x * b.x + a.y * b.y + a.z * b.z);
}

function vlength(a)
{
    return Math.sqrt(a.x * a.x + a.y * a.y + a.z * a.z);
}

function vnormalize(a)
{
    var len = vlength(a);
    var v = new vec(a.x, a.y, a.z);

    if (Math.abs(len) > 1.0e-17) {
        v.x /= len;
        v.y /= len;
        v.z /= len;
    }

    return v;
}

function Sphere(center, radius)
{
    this.center = center;
    this.radius = radius;

    this.intersect = function (ray, isect) {

        // rs = ray.org - sphere.center
        var rs  = vsub(ray.org, this.center);
        var B = vdot(rs, ray.dir);
        var C = vdot(rs, rs) - (this.radius * this.radius);
        var D = B * B - C;


        if (D > 0.0) {
            var t = -B - Math.sqrt(D);

            if ( (t > 0.0) && (t < isect.t) ) {
                isect.t   = t;
                isect.hit = true;

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

                // calculate normal.
                var n = vsub(isect.p, this.center);
                isect.n = vnormalize(n);

            }
        }
    }
}

function Plane(p, n)
{
    this.p = p;
    this.n = n;

    this.intersect = function (ray, isect) {
        var d = -vdot(this.p, this.n);
        var v = vdot(ray.dir, this.n);

        if (Math.abs(v) < 1.0e-17) return;      // no hit

        var t = -(vdot(ray.org, this.n) + d) / v;

        if ( (t > 0.0) && (t < isect.t) ) {
            isect.hit = true;
            isect.t   = t;
            isect.n   = this.n;
            isect.p   = new vec( ray.org.x + t * ray.dir.x,
                                 ray.org.y + t * ray.dir.y,
                                 ray.org.z + t * ray.dir.z );
        }
    }
}

function Ray(org, dir)
{
    this.org = org;
    this.dir = dir;
}

function Isect()
{
    this.t   = 1000000.0;     // far away
    this.hit = false;
    this.p   = new vec(0.0, 0.0, 0.0)
    this.n   = new vec(0.0, 0.0, 0.0)
}

function clamp(f)
{
    
    i = f * 255.5;
    if (i > 255.0) i = 255.0;
    if (i < 0.0) i = 0.0;
    return Math.round(i)
}

function orthoBasis(basis, n)
{
    basis[2] = new vec(n.x, n.y, n.z)
    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] = vcross(basis[1], basis[2]);
    basis[0] = vnormalize(basis[0]);

    basis[1] = vcross(basis[2], basis[0]);
    basis[1] = vnormalize(basis[1]);
}

var spheres;
var plane;

function init_scene()
{
    spheres    = new Array(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));
}


function ambient_occlusion(isect)
{
    var basis = new Array(3);
    orthoBasis(basis,  isect.n);

    var ntheta    = NAO_SAMPLES;
    var nphi      = NAO_SAMPLES;
    var eps       = 0.0001;
    var occlusion = 0.0;

    var p = new vec(isect.p.x + eps * isect.n.x,
                    isect.p.y + eps * isect.n.y,
                    isect.p.z + eps * isect.n.z);

    for (j = 0; j < nphi; j++) {
        for (i = 0; i < ntheta; i++) {
            var r = Math.random();
            var phi = 2.0 * Math.PI * Math.random();

            var x = Math.cos(phi) * Math.sqrt(1.0 - r);
            var y = Math.sin(phi) * Math.sqrt(1.0 - r);
            var z = Math.sqrt(r);


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

            var raydir = new vec(rx, ry, rz);
            var ray = new Ray(p, raydir);

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

            if (occIsect.hit) occlusion += 1.0;
        }
    }
    
    // [0.0, 1.0]
    occlusion = (ntheta * nphi - occlusion) / (ntheta * nphi);

    return new vec(occlusion, occlusion, occlusion);
}


function render(buf, w, h, nsubsamples)
{
    cnt = 0;
    for (y = 0; y < h; y++) {
        for (x = 0; x < w; x++) {

            rad = new vec(0.0, 0.0, 0.0);

            // subsampling
            for (v = 0; v < nsubsamples; v++) {
                for (u = 0; u < nsubsamples; u++) {

                    cnt++;
                    px =  (x + (u / nsubsamples) - (w / 2.0)) / (w / 2.0);
                    py = -(y + (v / nsubsamples) - (h / 2.0)) / (h / 2.0);

                    eye = vnormalize(new vec(px, py, -1.0));

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

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


                    if (isect.hit) {
                            
                        col = ambient_occlusion(isect);

                        rad.x += col.x;
                        rad.y += col.y;
                        rad.z += col.z;
                    }
                }
            }

            r = rad.x / (nsubsamples * nsubsamples);
            g = rad.y / (nsubsamples * nsubsamples);
            b = rad.z / (nsubsamples * nsubsamples);

            // use fill rect
//            ctx.fillStyle = "rgb(" + clamp(r) + "," + clamp(g) + "," + clamp(b) + ")";
//            ctx.fillRect (x, y, 1, 1);
            buf.push(clamp(r), clamp(g), clamp(b));
            
        }
    }

}


var onNode = (typeof(process) == 'object');

function page_onload()
{

    if (onNode){
      // ao.ppm
      process.stdout.write("P6\n");
      process.stdout.write(IMAGE_WIDTH + " " + IMAGE_HEIGHT + "\n");
      process.stdout.write("255\n");
    }
    else {
      document.elapform.elapsec.value = "Running..."
    }

    var buf = [];

    var elapsed = 0;
    var start = new Date();
    init_scene();
    render(buf, IMAGE_WIDTH, IMAGE_HEIGHT, 1)

    elapsed = new Date() - start;

    if (onNode) {
      process.stderr.write("Time: " + (elapsed / 1000.0) + "\n");

      var buffer = new Buffer(buf.length);
      for(var i=0; i<buf.length; i++){
        buffer.writeUInt8(buf[i], i);
        if(!IMGOUT) util.puts(buf[i]);
      }
      if(IMGOUT) process.stdout.write(buffer);
    }
    else {
      document.elapform.elapsec.value = elapsed / 1000.0
    }


    //img = ctx.getImageData(10, 10, 50, 50)
    //document.write(img.data[41]);
    //ret = ctx.putImagedata(img, 10, 10);
    //print(ret);
}

if (onNode) {
  page_onload();
}