AO bench

http://lucille.atso-net.jp/aobench/

噂のAO benchをやってみようと、JavaScript版をXtalに移植してみました。


Luaのソースを落としてきて実行してみたところ、
私の実行環境では、Luaの実行時間259秒に対し、Xtalでは498秒かかりました。
およそ2倍遅いということですね。
挽回して同程度にはしていきたいところです。


一時オブジェクトを大量に作っているので、オブジェクトを使いまわすように書き直せば、かなり速くなると思われましたが、Luaのソースでも一時オブジェクトを作っていたので、比較のためにそれはしませんでした。

Squirrel版も書いて比較したいところ。


以下Xtalソース。
最新のリビジョンでコンパイル、実行できます。

IMAGE_WIDTH: 256;
IMAGE_HEIGHT: 256;
NSUBSAMPLES: 2;
NAO_SAMPLES: 8;

class vec{
	public _x;
	public _y;
	public _z;
	
	initialize(_x, _y, _z){
	}
	
	op_add(a){
		return vec(_x+a.x, _y+a.y, _z+a.z);
	}
	
	op_sub(a){
		return vec(_x-a.x, _y-a.y, _z-a.z);
	}
	
	length(){
		return math::sqrt(_x * _x + _y * _y + _z * _z);
	}
	
	normalize(){
		len: length();
		if(math::abs(len) > 1.0e-17) {
		   return vec(_x/len, _y/len, _z/len);
		}
		return this;	
	}
	
	cross(a){
		return vec(_y*a.z-_z*a.y, _z*a.x-_x*a.z, _x*a.y-_y*a.x);
	}

	dot(a){
		return _x * a.x + _y * a.y + _z * a.z;
	}
}

class Sphere{
	public _center;
	public _radius;
	
	initialize(_center, _radius){
	}
	
	intersect(ray, isect){
		rs:  ray.org - _center;
		B: rs.dot(ray.dir);
		C: rs.dot(rs) - (_radius * _radius);
		D: B * B - C;

		if(D > 0.0){
			t: -B - math::sqrt(D);

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

class Plane{
	public _p;
	public _n;

	initialize(_p, _n){
	}

	intersect(ray, isect){
		d: -_p.dot(_n);
		v: ray.dir.dot(_n);

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

		t: -(ray.org.dot(_n) + d) / v;

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

class Ray{
	public _org;
	public _dir;

	initialize(_org, _dir){
	}
}

class Isect{
	public _t: 1000000.0;
	public _hit: false;
	public _p: vec(0.0, 0.0, 0.0);
	public _n: vec(0.0, 0.0, 0.0);
}

fun clamp(f){
	i: f * 255.5;
	if(i > 255.0) i = 255.0;
	if(i < 0.0) i = 0.0;
	return (i+0.5).to_i; //math::round(i);
}

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

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

spheres: null;
plane: null;

fun init_scene(){
	spheres = [];
	spheres.resize(3);
	spheres[0] = Sphere(vec(-2.0, 0.0, -3.5), 0.5);
	spheres[1] = Sphere(vec(-0.5, 0.0, -3.0), 0.5);
	spheres[2] = Sphere(vec(1.0, 0.0, -2.2), 0.5);
	plane = Plane(vec(0.0, -0.5, 0.0), vec(0.0, 1.0, 0.0));
}

fun ambient_occlusion(isect){
	basis: [];
	basis.resize(3);
	orthoBasis(basis,  isect.n);

	ntheta: NAO_SAMPLES;
	nphi: NAO_SAMPLES;
	eps: 0.0001;
	occlusion: 0.0;

	p: 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++) {
			r: math::random();
			phi: 2.0 * math::PI * math::random();

			x: math::cos(phi) * math::sqrt(1.0 - r);
			y: math::sin(phi) * math::sqrt(1.0 - r);
			z: math::sqrt(r);

			rx: x * basis[0].x + y * basis[1].x + z * basis[2].x;
			ry: x * basis[0].y + y * basis[1].y + z * basis[2].y;
			rz: x * basis[0].z + y * basis[1].z + z * basis[2].z;

			raydir: vec(rx, ry, rz);
			ray: Ray(p, raydir);

			occIsect: 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;
			}
		}
	}

	occlusion = (ntheta * nphi - occlusion) / (ntheta * nphi);
	return vec(occlusion, occlusion, occlusion);
}

fun render(out, w, h, nsubsamples){
	cnt: 0;

	for(y: 0; y < h; y++) {
		for(x: 0; x < w; x++) {
			rad: vec(0.0, 0.0, 0.0);

			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: vec(px, py, -1.0).normalize;

					ray: Ray(vec(0.0, 0.0, 0.0), eye);

					isect: 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);
			
			out.put_u8(clamp(r));
			out.put_u8(clamp(g));
			out.put_u8(clamp(b));
		}
		
		%f!%s\n!(y).p;
		gc();
	}
}

init_scene();

open("out.ppm", "w"){
	it.put_s("P6\n");
	it.put_s(%f(%d %d\n)(IMAGE_WIDTH, IMAGE_HEIGHT));
	it.put_s("255\n");
	render(it, IMAGE_WIDTH, IMAGE_HEIGHT, NSUBSAMPLES);
}