# HG changeset patch
# User John Tsiombikas <nuclear@member.fsf.org>
# Date 1284133620 -3600
# Node ID 6a30f27fa1e6a7425d5d4a25730f2883fce2a8bc
# Parent  54a96b738afe71c06ae62ca0eb69493ff4cebb79
separated the OpenGL visualization and added a CPU raytracing mode

diff -r 54a96b738afe -r 6a30f27fa1e6 rt.cl
--- a/rt.cl	Sun Sep 05 16:43:55 2010 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,395 +0,0 @@
-/* vim: set ft=opencl:ts=4:sw=4 */
-#include "common.h"
-
-struct RendInfo {
-	float4 ambient;
-	int xsz, ysz;
-	int num_faces, num_lights;
-	int max_iter;
-	int cast_shadows;
-};
-
-struct Vertex {
-	float4 pos;
-	float4 normal;
-	float4 tex;
-	float4 padding;
-};
-
-struct Face {
-	struct Vertex v[3];
-	float4 normal;
-	int matid;
-	int padding[3];
-};
-
-struct Material {
-	float4 kd, ks;
-	float kr, kt;
-	float spow;
-	float padding;
-};
-
-struct Light {
-	float4 pos, color;
-};
-
-struct Ray {
-	float4 origin, dir;
-};
-
-struct SurfPoint {
-	float t;
-	float4 pos, norm, dbg;
-	global const struct Face *obj;
-	struct Material mat;
-};
-
-struct Scene {
-	float4 ambient;
-	global const struct Face *faces;
-	int num_faces;
-	global const struct Light *lights;
-	int num_lights;
-	global const struct Material *matlib;
-	//global const struct KDNode *kdtree;
-	bool cast_shadows;
-};
-
-struct AABBox {
-	float4 min, max;
-};
-
-struct KDNode {
-	struct AABBox aabb;
-	int face_idx[MAX_NODE_FACES];
-	int num_faces;
-	int left, right;
-	int padding;
-};
-
-#define MIN_ENERGY	0.001
-#define EPSILON		1e-5
-
-float4 shade(struct Ray ray, struct Scene *scn, const struct SurfPoint *sp, read_only image2d_t kdimg);
-bool find_intersection(struct Ray ray, const struct Scene *scn, struct SurfPoint *sp, read_only image2d_t kdimg);
-bool intersect(struct Ray ray, global const struct Face *face, struct SurfPoint *sp);
-bool intersect_aabb(struct Ray ray, struct AABBox aabb);
-
-float4 reflect(float4 v, float4 n);
-float4 transform(float4 v, global const float *xform);
-void transform_ray(struct Ray *ray, global const float *xform, global const float *invtrans);
-float4 calc_bary(float4 pt, global const struct Face *face, float4 norm);
-float mean(float4 v);
-
-void read_kdnode(int idx, struct KDNode *node, read_only image2d_t kdimg);
-
-
-kernel void render(write_only image2d_t fb,
-		global const struct RendInfo *rinf,
-		global const struct Face *faces,
-		global const struct Material *matlib,
-		global const struct Light *lights,
-		global const struct Ray *primrays,
-		global const float *xform,
-		global const float *invtrans,
-		//global const struct KDNode *kdtree
-		read_only image2d_t kdtree_img)
-{
-	int idx = get_global_id(0);
-
-	struct Scene scn;
-	scn.ambient = rinf->ambient;
-	scn.faces = faces;
-	scn.num_faces = rinf->num_faces;
-	scn.lights = lights;
-	scn.num_lights = rinf->num_lights;
-	scn.matlib = matlib;
-	scn.cast_shadows = rinf->cast_shadows;
-
-	struct Ray ray = primrays[idx];
-	transform_ray(&ray, xform, invtrans);
-
-	float4 pixel = (float4)(0, 0, 0, 0);
-	float4 energy = (float4)(1.0, 1.0, 1.0, 0.0);
-	int iter = 0;
-
-	while(iter++ <= rinf->max_iter && mean(energy) > MIN_ENERGY) {
-		struct SurfPoint sp;
-		if(find_intersection(ray, &scn, &sp, kdtree_img)) {
-			pixel += shade(ray, &scn, &sp, kdtree_img) * energy;
-
-			float4 refl_col = sp.mat.ks * sp.mat.kr;
-
-			ray.origin = sp.pos;
-			ray.dir = reflect(-ray.dir, sp.norm);
-
-			energy *= refl_col;
-		} else {
-			energy = (float4)(0.0, 0.0, 0.0, 0.0);
-		}
-	}
-
-	int2 coord;
-	coord.x = idx % rinf->xsz;
-	coord.y = idx / rinf->xsz;
-
-	write_imagef(fb, coord, pixel);
-}
-
-float4 shade(struct Ray ray, struct Scene *scn, const struct SurfPoint *sp, read_only image2d_t kdimg)
-{
-	float4 norm = sp->norm;
-	//bool entering = true;
-
-	if(dot(ray.dir, norm) >= 0.0) {
-		norm = -norm;
-		//entering = false;
-	}
-
-	float4 dcol = scn->ambient * sp->mat.kd;
-	float4 scol = (float4)(0, 0, 0, 0);
-
-	for(int i=0; i<scn->num_lights; i++) {
-		float4 ldir = scn->lights[i].pos - sp->pos;
-
-		struct Ray shadowray;
-		shadowray.origin = sp->pos;
-		shadowray.dir = ldir;
-
-		if(!scn->cast_shadows || !find_intersection(shadowray, scn, 0, kdimg)) {
-			ldir = normalize(ldir);
-			float4 vdir = -ray.dir;
-			vdir.x = native_divide(vdir.x, RAY_MAG);
-			vdir.y = native_divide(vdir.y, RAY_MAG);
-			vdir.z = native_divide(vdir.z, RAY_MAG);
-			float4 vref = reflect(vdir, norm);
-
-			float diff = fmax(dot(ldir, norm), 0.0f);
-			dcol += sp->mat.kd /* scn->lights[i].color*/ * diff;
-
-			float spec = native_powr(fmax(dot(ldir, vref), 0.0f), sp->mat.spow);
-			scol += sp->mat.ks /* scn->lights[i].color*/ * spec;
-		}
-	}
-
-	return dcol + scol;
-}
-
-#define STACK_SIZE	MAX_TREE_DEPTH
-bool find_intersection(struct Ray ray, const struct Scene *scn, struct SurfPoint *spres, read_only image2d_t kdimg)
-{
-	struct SurfPoint sp0;
-	sp0.t = 1.0;
-	sp0.obj = 0;
-
-	int idxstack[STACK_SIZE];
-	int top = 0;			// points after the topmost element of the stack
-	idxstack[top++] = 0;	// root at tree[0]
-
-	while(top > 0) {
-		int idx = idxstack[--top];	// remove this index from the stack and process it
-
-		struct KDNode node;
-		read_kdnode(idx, &node, kdimg);
-
-		if(intersect_aabb(ray, node.aabb)) {
-			if(node.left == -1) {
-				// leaf node... check each face in turn and update the nearest intersection as needed
-				for(int i=0; i<node.num_faces; i++) {
-					struct SurfPoint spt;
-					int fidx = node.face_idx[i];
-
-					if(intersect(ray, scn->faces + fidx, &spt) && spt.t < sp0.t) {
-						sp0 = spt;
-					}
-				}
-			} else {
-				// internal node... recurse to the children
-				idxstack[top++] = node.left;
-				idxstack[top++] = node.right;
-			}
-		}
-	}
-
-	if(!sp0.obj) {
-		return false;
-	}
-
-	if(spres) {
-		*spres = sp0;
-		spres->mat = scn->matlib[sp0.obj->matid];
-	}
-	return true;
-}
-
-bool intersect(struct Ray ray, global const struct Face *face, struct SurfPoint *sp)
-{
-	float4 origin = ray.origin;
-	float4 dir = ray.dir;
-	float4 norm = face->normal;
-
-	float ndotdir = dot(dir, norm);
-
-	if(fabs(ndotdir) <= EPSILON) {
-		return false;
-	}
-
-	float4 pt = face->v[0].pos;
-	float4 vec = pt - origin;
-
-	float ndotvec = dot(norm, vec);
-	float t = native_divide(ndotvec, ndotdir);
-
-	if(t < EPSILON || t > 1.0) {
-		return false;
-	}
-	pt = origin + dir * t;
-
-
-	float4 bc = calc_bary(pt, face, norm);
-	float bc_sum = bc.x + bc.y + bc.z;
-
-	if(bc_sum < 1.0 - EPSILON || bc_sum > 1.0 + EPSILON) {
-		return false;
-		bc *= 1.2;
-	}
-
-	sp->t = t;
-	sp->pos = pt;
-	sp->norm = normalize(face->v[0].normal * bc.x + face->v[1].normal * bc.y + face->v[2].normal * bc.z);
-	sp->obj = face;
-	sp->dbg = bc;
-	return true;
-}
-
-bool intersect_aabb(struct Ray ray, struct AABBox aabb)
-{
-	if(ray.origin.x >= aabb.min.x && ray.origin.y >= aabb.min.y && ray.origin.z >= aabb.min.z &&
-			ray.origin.x < aabb.max.x && ray.origin.y < aabb.max.y && ray.origin.z < aabb.max.z) {
-		return true;
-	}
-
-	float4 bbox[2] = {
-		aabb.min.x, aabb.min.y, aabb.min.z, 0,
-		aabb.max.x, aabb.max.y, aabb.max.z, 0
-	};
-
-	int xsign = (int)(ray.dir.x < 0.0);
-	float invdirx = native_recip(ray.dir.x);
-	float tmin = (bbox[xsign].x - ray.origin.x) * invdirx;
-	float tmax = (bbox[1 - xsign].x - ray.origin.x) * invdirx;
-
-	int ysign = (int)(ray.dir.y < 0.0);
-	float invdiry = native_recip(ray.dir.y);
-	float tymin = (bbox[ysign].y - ray.origin.y) * invdiry;
-	float tymax = (bbox[1 - ysign].y - ray.origin.y) * invdiry;
-
-	if(tmin > tymax || tymin > tmax) {
-		return false;
-	}
-
-	if(tymin > tmin) tmin = tymin;
-	if(tymax < tmax) tmax = tymax;
-
-	int zsign = (int)(ray.dir.z < 0.0);
-	float invdirz = native_recip(ray.dir.z);
-	float tzmin = (bbox[zsign].z - ray.origin.z) * invdirz;
-	float tzmax = (bbox[1 - zsign].z - ray.origin.z) * invdirz;
-
-	if(tmin > tzmax || tzmin > tmax) {
-		return false;
-	}
-
-	return tmin < 1.0 && tmax > 0.0;
-}
-
-float4 reflect(float4 v, float4 n)
-{
-	return 2.0f * dot(v, n) * n - v;
-}
-
-float4 transform(float4 v, global const float *xform)
-{
-	float4 res;
-	res.x = v.x * xform[0] + v.y * xform[4] + v.z * xform[8] + xform[12];
-	res.y = v.x * xform[1] + v.y * xform[5] + v.z * xform[9] + xform[13];
-	res.z = v.x * xform[2] + v.y * xform[6] + v.z * xform[10] + xform[14];
-	res.w = 0.0;
-	return res;
-}
-
-void transform_ray(struct Ray *ray, global const float *xform, global const float *invtrans)
-{
-	ray->origin = transform(ray->origin, xform);
-	ray->dir = transform(ray->dir, invtrans);
-}
-
-float4 calc_bary(float4 pt, global const struct Face *face, float4 norm)
-{
-	float4 bc = (float4)(0, 0, 0, 0);
-
-	// calculate area of the whole triangle
-	float4 v1 = face->v[1].pos - face->v[0].pos;
-	float4 v2 = face->v[2].pos - face->v[0].pos;
-	float4 xv1v2 = cross(v1, v2);
-
-	float area = fabs(dot(xv1v2, norm)) * 0.5;
-	if(area < EPSILON) {
-		return bc;
-	}
-
-	float4 pv0 = face->v[0].pos - pt;
-	float4 pv1 = face->v[1].pos - pt;
-	float4 pv2 = face->v[2].pos - pt;
-
-	// calculate the area of each sub-triangle
-	float4 x12 = cross(pv1, pv2);
-	float4 x20 = cross(pv2, pv0);
-	float4 x01 = cross(pv0, pv1);
-
-	float a0 = fabs(dot(x12, norm)) * 0.5;
-	float a1 = fabs(dot(x20, norm)) * 0.5;
-	float a2 = fabs(dot(x01, norm)) * 0.5;
-
-	bc.x = native_divide(a0, area);
-	bc.y = native_divide(a1, area);
-	bc.z = native_divide(a2, area);
-	return bc;
-}
-
-float mean(float4 v)
-{
-	return native_divide(v.x + v.y + v.z, 3.0);
-}
-
-
-const sampler_t kdsampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_NONE | CLK_FILTER_NEAREST;
-
-// read a KD-tree node from a texture scanline
-void read_kdnode(int idx, struct KDNode *node, read_only image2d_t kdimg)
-{
-	int startx = KDIMG_NODE_WIDTH * (idx / KDIMG_MAX_HEIGHT);
-
-	int2 tc;
-	tc.x = startx;
-	tc.y = idx % KDIMG_MAX_HEIGHT;
-
-	node->aabb.min = read_imagef(kdimg, kdsampler, tc); tc.x++;
-	node->aabb.max = read_imagef(kdimg, kdsampler, tc);
-
-	tc.x = startx + 2 + MAX_NODE_FACES / 4;
-	float4 pix = read_imagef(kdimg, kdsampler, tc);
-	node->num_faces = (int)pix.x;
-	node->left = (int)pix.y;
-	node->right = (int)pix.z;
-
-	tc.x = startx + 2;
-	for(int i=0; i<node->num_faces; i+=4) {
-		float4 pix = read_imagef(kdimg, kdsampler, tc); tc.x++;
-		node->face_idx[i] = (int)pix.x;
-		node->face_idx[i + 1] = (int)pix.y;
-		node->face_idx[i + 2] = (int)pix.z;
-		node->face_idx[i + 3] = (int)pix.w;
-	}
-}
diff -r 54a96b738afe -r 6a30f27fa1e6 src/clray.cc
--- a/src/clray.cc	Sun Sep 05 16:43:55 2010 +0100
+++ b/src/clray.cc	Fri Sep 10 16:47:00 2010 +0100
@@ -33,14 +33,21 @@
 static float cam_theta, cam_phi = 25.0;
 static float cam_dist = 10.0;
 
-static bool dbg_glrender = false;
-static bool dbg_show_kdtree = false;
+static bool dbg_glrender;
+static bool dbg_nocl;
+static bool dbg_show_kdtree;
 static bool dbg_show_obj = true;
 bool dbg_frame_time = true;
 
 static Scene scn;
 static unsigned int tex;
 
+static Light lightlist[] = {
+	{{-8, 15, 18, 0}, {1, 1, 1, 1}}
+};
+
+
+
 int main(int argc, char **argv)
 {
 	glutInitWindowSize(800, 600);
@@ -81,6 +88,10 @@
 				dbg_glrender = true;
 				break;
 
+			case 'n':
+				dbg_nocl = true;
+				break;
+
 			default:
 				fprintf(stderr, "unrecognized option: %s\n", argv[i]);
 				return 1;
@@ -103,6 +114,11 @@
 		return false;
 	}
 
+	int num_lights = sizeof lightlist / sizeof *lightlist;
+	for(int i=0; i<num_lights; i++) {
+		scn.add_light(lightlist[i]);
+	}
+
 	glutInitDisplayMode(GLUT_RGB | GLUT_DEPTH | GLUT_DOUBLE);
 	glutCreateWindow("OpenCL Raytracer");
 
@@ -171,11 +187,6 @@
 
 		glGetFloatv(GL_MODELVIEW_MATRIX, mat.m);
 
-		inv_mat = mat;
-		inv_mat.invert();
-
-		/*inv_trans = inv_mat;
-		inv_trans.transpose();*/
 		inv_trans = mat;
 		inv_trans.m[3] = inv_trans.m[7] = inv_trans.m[11] = 0.0;
 		inv_trans.m[12] = inv_trans.m[13] = inv_trans.m[14] = 0.0;
@@ -185,14 +196,21 @@
 		glPopMatrix();
 
 		if(!dbg_glrender) {
-			if(!render()) {
-				exit(1);
+			if(dbg_nocl) {
+				dbg_render(mat.m, inv_trans.m);
+			} else {
+				if(!render()) {
+					exit(1);
+				}
 			}
 			need_update = false;
 		}
 	}
 
 	if(dbg_glrender) {
+		inv_mat = mat;
+		inv_mat.invert();
+
 		glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
 		glLoadMatrixf(inv_mat.m);
 		dbg_render_gl(&scn, dbg_show_kdtree, dbg_show_obj);
@@ -322,6 +340,13 @@
 		dbg_frame_time = !dbg_frame_time;
 		break;
 
+	case 'n':
+		dbg_nocl = !dbg_nocl;
+		printf("switching to %s rendering\n", dbg_nocl ? "debug CPU" : "OpenCL");
+		need_update = true;
+		glutPostRedisplay();
+		break;
+
 	default:
 		break;
 	}
diff -r 54a96b738afe -r 6a30f27fa1e6 src/common.h
--- a/src/common.h	Sun Sep 05 16:43:55 2010 +0100
+++ b/src/common.h	Fri Sep 10 16:47:00 2010 +0100
@@ -1,6 +1,12 @@
 #ifndef COMMON_H_
 #define COMMON_H_
 
+/* error threshold */
+#define EPSILON		1e-5
+
+/* minimum trace energy threshold */
+#define MIN_ENERGY	0.001
+
 /* primary ray magnitude */
 #define RAY_MAG		500.0
 
diff -r 54a96b738afe -r 6a30f27fa1e6 src/dbggl.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/dbggl.cc	Fri Sep 10 16:47:00 2010 +0100
@@ -0,0 +1,77 @@
+#include <string.h>
+#include <assert.h>
+#include "rt.h"
+#include "ogl.h"
+
+#define MIN(a, b)	((a) < (b) ? (a) : (b))
+static void dbg_set_gl_material(Material *mat)
+{
+	static Material def_mat = {{0.7, 0.7, 0.7, 1}, {0, 0, 0, 0}, 0, 0, 0};
+
+	if(!mat) mat = &def_mat;
+
+	glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat->kd);
+	glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, mat->ks);
+	glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, MIN(mat->spow, 128.0f));
+}
+
+void dbg_render_gl(Scene *scn, bool show_tree, bool show_obj)
+{
+	const RendInfo *rinf = get_render_info();
+	const Face *faces = scn->get_face_buffer();
+
+	glPushAttrib(GL_ENABLE_BIT | GL_TRANSFORM_BIT | GL_LIGHTING_BIT);
+
+	for(int i=0; i<scn->get_num_lights(); i++) {
+		float lpos[4];
+
+		memcpy(lpos, scn->lights[i].pos, sizeof lpos);
+		lpos[3] = 1.0;
+
+		glLightfv(GL_LIGHT0 + i, GL_POSITION, lpos);
+		glLightfv(GL_LIGHT0 + i, GL_DIFFUSE, scn->lights[i].color);
+		glEnable(GL_LIGHT0 + i);
+	}
+
+	glDisable(GL_TEXTURE_2D);
+	glEnable(GL_DEPTH_TEST);
+	glEnable(GL_LIGHTING);
+
+	glMatrixMode(GL_PROJECTION);
+	glPushMatrix();
+	glLoadIdentity();
+	gluPerspective(45.0, (float)rinf->xsz / (float)rinf->ysz, 0.5, 1000.0);
+
+	if(show_obj) {
+		Material *materials = scn->get_materials();
+
+		int num_faces = scn->get_num_faces();
+		int cur_mat = -1;
+
+		for(int i=0; i<num_faces; i++) {
+			if(faces[i].matid != cur_mat) {
+				if(cur_mat != -1) {
+					glEnd();
+				}
+				dbg_set_gl_material(materials ? materials + faces[i].matid : 0);
+				cur_mat = faces[i].matid;
+				glBegin(GL_TRIANGLES);
+			}
+
+			for(int j=0; j<3; j++) {
+				glNormal3fv(faces[i].v[j].normal);
+				glVertex3fv(faces[i].v[j].pos);
+			}
+		}
+		glEnd();
+	}
+
+	if(show_tree) {
+		scn->draw_kdtree();
+	}
+
+	glPopMatrix();
+	glPopAttrib();
+
+	assert(glGetError() == GL_NO_ERROR);
+}
diff -r 54a96b738afe -r 6a30f27fa1e6 src/dbgray.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/dbgray.cc	Fri Sep 10 16:47:00 2010 +0100
@@ -0,0 +1,342 @@
+#include <string.h>
+#include <assert.h>
+#include "rt.h"
+#include "ogl.h"
+#include "vector.h"
+#include "timer.h"
+
+struct SurfPoint {
+	float t;
+	Vector3 pos, norm;
+	const Face *face;
+};
+
+static void trace_ray(float *pixel, const Ray &ray, int iter, float energy = 1.0f);
+static void shade(float *pixel, const Ray &ray, const SurfPoint &sp, int iter, float energy = 1.0f);
+static bool find_intersection(const Ray &ray, const Scene *scn, const KDNode *kd, SurfPoint *spret);
+static bool ray_aabb_test(const Ray &ray, const AABBox &aabb);
+static bool ray_triangle_test(const Ray &ray, const Face *face, SurfPoint *sp);
+static Vector3 calc_bary(const Vector3 &pt, const Face *face, const Vector3 &norm);
+static void transform(float *res, const float *v, const float *xform);
+static void transform_ray(Ray *ray, const float *xform, const float *invtrans_xform);
+
+static int xsz, ysz;
+static float *fb;
+static unsigned int tex;
+static Scene *scn;
+static const Ray *prim_rays;
+
+
+bool init_dbg_renderer(int width, int height, Scene *scene, unsigned int texid)
+{
+	try {
+		fb = new float[3 * width * height];
+	}
+	catch(...) {
+		return false;
+	}
+
+	xsz = width;
+	ysz = height;
+	tex = texid;
+	scn = scene;
+	return true;
+}
+
+void destroy_dbg_renderer()
+{
+	delete [] fb;
+	delete [] prim_rays;
+}
+
+void dbg_set_primary_rays(const Ray *rays)
+{
+	prim_rays = rays;
+}
+
+void dbg_render(const float *xform, const float *invtrans_xform, int num_threads)
+{
+	unsigned long t0 = get_msec();
+
+	int iter = get_render_option_int(ROPT_ITER);
+
+	int offs = 0;
+	for(int i=0; i<ysz; i++) {
+		for(int j=0; j<xsz; j++) {
+			Ray ray = prim_rays[offs];
+			transform_ray(&ray, xform, invtrans_xform);
+
+			trace_ray(fb + offs * 3, ray, iter, 1.0);
+			offs++;
+		}
+	}
+
+	glPushAttrib(GL_TEXTURE_BIT);
+	glBindTexture(GL_TEXTURE_2D, tex);
+	glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, xsz, ysz, GL_RGB, GL_FLOAT, fb);
+	glPopAttrib();
+
+	printf("rendered in %lu msec\n", get_msec() - t0);
+}
+
+static void trace_ray(float *pixel, const Ray &ray, int iter, float energy)
+{
+	SurfPoint sp;
+
+	if(find_intersection(ray, scn, scn->kdtree, &sp)) {
+		shade(pixel, ray, sp, iter, energy);
+	} else {
+		pixel[0] = pixel[1] = pixel[2] = 0.05f;
+	}
+}
+
+#define MAX(a, b)	((a) > (b) ? (a) : (b))
+
+static void shade(float *pixel, const Ray &ray, const SurfPoint &sp, int iter, float energy)
+{
+	const Material *mat = scn->get_materials() + sp.face->matid;
+
+	bool cast_shadows = get_render_option_bool(ROPT_SHAD);
+	Vector3 raydir(ray.dir);
+	Vector3 norm = sp.norm;
+
+	if(dot(raydir, norm) >= 0.0) {
+		norm = -norm;
+	}
+
+	float dcol[3] = {0, 0, 0};
+	float scol[3] = {0, 0, 0};
+
+	for(int i=0; i<scn->get_num_lights(); i++) {
+		Vector3 lpos(scn->lights[i].pos);
+		Vector3 ldir = lpos - sp.pos;
+
+		Ray shadowray;
+		shadowray.origin[0] = sp.pos.x;
+		shadowray.origin[1] = sp.pos.y;
+		shadowray.origin[2] = sp.pos.z;
+		shadowray.dir[0] = ldir.x;
+		shadowray.dir[1] = ldir.y;
+		shadowray.dir[2] = ldir.z;
+
+		if(!cast_shadows || !find_intersection(shadowray, scn, scn->kdtree, 0)) {
+
+			ldir.normalize();
+
+			Vector3 vdir = -raydir / RAY_MAG;
+			Vector3 vref = reflect(vdir, norm);
+
+			float ndotl = dot(ldir, norm);
+			float diff = MAX(ndotl, 0.0f);
+
+			dcol[0] += mat->kd[0] * diff;
+			dcol[1] += mat->kd[1] * diff;
+			dcol[2] += mat->kd[2] * diff;
+
+			float ldotvr = dot(ldir, vref);
+			float spec = pow(MAX(ldotvr, 0.0f), mat->spow);
+
+			scol[0] += mat->ks[0] * spec;
+			scol[1] += mat->ks[1] * spec;
+			scol[2] += mat->ks[2] * spec;
+		}
+	}
+
+	float refl_color[3];
+	refl_color[0] = mat->ks[0] * mat->kr;
+	refl_color[1] = mat->ks[1] * mat->kr;
+	refl_color[2] = mat->ks[2] * mat->kr;
+
+	energy *= (refl_color[0] + refl_color[1] + refl_color[2]) / 3.0;
+	if(iter >= 0 && energy > MIN_ENERGY) {
+		Vector3 rdir = reflect(-raydir, norm);
+
+		Ray refl;
+		refl.origin[0] = sp.pos.x;
+		refl.origin[1] = sp.pos.y;
+		refl.origin[2] = sp.pos.z;
+		refl.dir[0] = rdir.x;
+		refl.dir[1] = rdir.y;
+		refl.dir[2] = rdir.z;
+
+		float rcol[3];
+		trace_ray(rcol, refl, iter - 1, energy);
+		scol[0] += rcol[0] * mat->ks[0] * mat->kr;
+		scol[1] += rcol[1] * mat->ks[1] * mat->kr;
+		scol[2] += rcol[2] * mat->ks[2] * mat->kr;
+	}
+
+	pixel[0] = dcol[0] + scol[0];
+	pixel[1] = dcol[1] + scol[1];
+	pixel[2] = dcol[2] + scol[2];
+}
+
+static bool find_intersection(const Ray &ray, const Scene *scn, const KDNode *kd, SurfPoint *spret)
+{
+	if(!ray_aabb_test(ray, kd->aabb)) {
+		return false;
+	}
+
+	SurfPoint sp, sptmp;
+	if(!spret) {
+		spret = &sptmp;
+	}
+
+	spret->t = RAY_MAG;
+	spret->face = 0;
+
+	if(kd->left) {
+		assert(kd->right);
+
+		bool found = find_intersection(ray, scn, kd->left, spret);
+		if(find_intersection(ray, scn, kd->right, &sp)) {
+			if(!found || sp.t < spret->t) {
+				*spret = sp;
+			}
+			found = true;
+		}
+		return found;
+	}
+
+	const Face *faces = scn->get_face_buffer();
+
+	for(size_t i=0; i<kd->face_idx.size(); i++) {
+		if(ray_triangle_test(ray, faces + kd->face_idx[i], &sp) && sp.t < spret->t) {
+			*spret = sp;
+		}
+	}
+	return spret->face != 0;
+}
+
+static bool ray_aabb_test(const Ray &ray, const AABBox &aabb)
+{
+	if(ray.origin[0] >= aabb.min[0] && ray.origin[1] >= aabb.min[1] && ray.origin[2] >= aabb.min[2] &&
+			ray.origin[0] < aabb.max[0] && ray.origin[1] < aabb.max[1] && ray.origin[2] < aabb.max[2]) {
+		return true;
+	}
+
+	float bbox[][3] = {
+		{aabb.min[0], aabb.min[1], aabb.min[2]},
+		{aabb.max[0], aabb.max[1], aabb.max[2]}
+	};
+
+	int xsign = (int)(ray.dir[0] < 0.0);
+	float invdirx = 1.0 / ray.dir[0];
+	float tmin = (bbox[xsign][0] - ray.origin[0]) * invdirx;
+	float tmax = (bbox[1 - xsign][0] - ray.origin[0]) * invdirx;
+
+	int ysign = (int)(ray.dir[1] < 0.0);
+	float invdiry = 1.0 / ray.dir[1];
+	float tymin = (bbox[ysign][1] - ray.origin[1]) * invdiry;
+	float tymax = (bbox[1 - ysign][1] - ray.origin[1]) * invdiry;
+
+	if(tmin > tymax || tymin > tmax) {
+		return false;
+	}
+
+	if(tymin > tmin) tmin = tymin;
+	if(tymax < tmax) tmax = tymax;
+
+	int zsign = (int)(ray.dir[2] < 0.0);
+	float invdirz = 1.0 / ray.dir[2];
+	float tzmin = (bbox[zsign][2] - ray.origin[2]) * invdirz;
+	float tzmax = (bbox[1 - zsign][2] - ray.origin[2]) * invdirz;
+
+	if(tmin > tzmax || tzmin > tmax) {
+		return false;
+	}
+
+	return tmin < 1.0 && tmax > 0.0;
+
+}
+
+static bool ray_triangle_test(const Ray &ray, const Face *face, SurfPoint *sp)
+{
+	Vector3 origin = ray.origin;
+	Vector3 dir = ray.dir;
+	Vector3 norm = face->normal;
+
+	float ndotdir = dot(dir, norm);
+
+	if(fabs(ndotdir) <= EPSILON) {
+		return false;
+	}
+
+	Vector3 pt = face->v[0].pos;
+	Vector3 vec = pt - origin;
+
+	float ndotvec = dot(norm, vec);
+	float t = ndotvec / ndotdir;
+
+	if(t < EPSILON || t > 1.0) {
+		return false;
+	}
+	pt = origin + dir * t;
+
+
+	Vector3 bc = calc_bary(pt, face, norm);
+	float bc_sum = bc.x + bc.y + bc.z;
+
+	if(bc_sum < 1.0 - EPSILON || bc_sum > 1.0 + EPSILON) {
+		return false;
+	}
+
+	Vector3 n0(face->v[0].normal);
+	Vector3 n1(face->v[1].normal);
+	Vector3 n2(face->v[2].normal);
+
+	sp->t = t;
+	sp->pos = pt;
+	sp->norm = n0 * bc.x + n1 * bc.y + n2 * bc.z;
+	sp->norm.normalize();
+	sp->face = face;
+	return true;
+}
+
+static Vector3 calc_bary(const Vector3 &pt, const Face *face, const Vector3 &norm)
+{
+	Vector3 bc(0.0f, 0.0f, 0.0f);
+
+	Vector3 v1 = Vector3(face->v[1].pos) - Vector3(face->v[0].pos);
+	Vector3 v2 = Vector3(face->v[2].pos) - Vector3(face->v[0].pos);
+	Vector3 xv1v2 = cross(v1, v2);
+
+	float area = fabs(dot(xv1v2, norm)) * 0.5;
+	if(area < EPSILON) {
+		return bc;
+	}
+
+	Vector3 pv0 = face->v[0].pos - pt;
+	Vector3 pv1 = face->v[1].pos - pt;
+	Vector3 pv2 = face->v[2].pos - pt;
+
+	// calculate the area of each sub-triangle
+	Vector3 x12 = cross(pv1, pv2);
+	Vector3 x20 = cross(pv2, pv0);
+	Vector3 x01 = cross(pv0, pv1);
+
+	float a0 = fabs(dot(x12, norm)) * 0.5;
+	float a1 = fabs(dot(x20, norm)) * 0.5;
+	float a2 = fabs(dot(x01, norm)) * 0.5;
+
+	bc.x = a0 / area;
+	bc.y = a1 / area;
+	bc.z = a2 / area;
+	return bc;
+
+}
+
+static void transform(float *res, const float *v, const float *xform)
+{
+	float tmp[3];
+	tmp[0] = v[0] * xform[0] + v[1] * xform[4] + v[2] * xform[8] + xform[12];
+	tmp[1] = v[0] * xform[1] + v[1] * xform[5] + v[2] * xform[9] + xform[13];
+	tmp[2] = v[0] * xform[2] + v[1] * xform[6] + v[2] * xform[10] + xform[14];
+	memcpy(res, tmp, sizeof tmp);
+}
+
+static void transform_ray(Ray *ray, const float *xform, const float *invtrans_xform)
+{
+	transform(ray->origin, ray->origin, xform);
+	transform(ray->dir, ray->dir, invtrans_xform);
+}
diff -r 54a96b738afe -r 6a30f27fa1e6 src/rt.cc
--- a/src/rt.cc	Sun Sep 05 16:43:55 2010 +0100
+++ b/src/rt.cc	Fri Sep 10 16:47:00 2010 +0100
@@ -24,22 +24,6 @@
 	NUM_KERNEL_ARGS
 };
 
-struct RendInfo {
-	float ambient[4];
-	int xsz, ysz;
-	int num_faces, num_lights;
-	int max_iter;
-	int cast_shadows;
-};
-
-struct Ray {
-	float origin[4], dir[4];
-};
-
-struct Light {
-	float pos[4], color[4];
-};
-
 static void update_render_info();
 static Ray get_primary_ray(int x, int y, int w, int h, float vfov_deg);
 static float *create_kdimage(const KDNodeGPU *kdtree, int num_nodes, int *xsz_ret, int *ysz_ret);
@@ -49,10 +33,6 @@
 static CLProgram *prog;
 static int global_size;
 
-static Light lightlist[] = {
-	{{-8, 15, 18, 0}, {1, 1, 1, 1}}
-};
-
 
 static RendInfo rinf;
 static int saved_iter_val;
@@ -72,7 +52,7 @@
 	rinf.xsz = xsz;
 	rinf.ysz = ysz;
 	rinf.num_faces = scn->get_num_faces();
-	rinf.num_lights = sizeof lightlist / sizeof *lightlist;
+	rinf.num_lights = scn->get_num_lights();
 	rinf.max_iter = saved_iter_val = 6;
 	rinf.cast_shadows = true;
 
@@ -84,10 +64,11 @@
 			prim_rays[i * xsz + j] = get_primary_ray(j, i, xsz, ysz, 45.0);
 		}
 	}
+	dbg_set_primary_rays(prim_rays);	// give them to the debug renderer
 
 	/* setup opencl */
 	prog = new CLProgram("render");
-	if(!prog->load("rt.cl")) {
+	if(!prog->load("src/rt.cl")) {
 		return false;
 	}
 
@@ -114,7 +95,7 @@
 	prog->set_arg_buffer(KARG_RENDER_INFO, ARG_RD, sizeof rinf, &rinf);
 	prog->set_arg_buffer(KARG_FACES, ARG_RD, rinf.num_faces * sizeof(Face), faces);
 	prog->set_arg_buffer(KARG_MATLIB, ARG_RD, scn->get_num_materials() * sizeof(Material), scn->get_materials());
-	prog->set_arg_buffer(KARG_LIGHTS, ARG_RD, sizeof lightlist, lightlist);
+	prog->set_arg_buffer(KARG_LIGHTS, ARG_RD, scn->get_num_lights() * sizeof(Light), scn->get_lights());
 	prog->set_arg_buffer(KARG_PRIM_RAYS, ARG_RD, xsz * ysz * sizeof *prim_rays, prim_rays);
 	prog->set_arg_buffer(KARG_XFORM, ARG_RD, 16 * sizeof(float));
 	prog->set_arg_buffer(KARG_INVTRANS_XFORM, ARG_RD, 16 * sizeof(float));
@@ -133,9 +114,12 @@
 		return false;
 	}
 
-	delete [] prim_rays;
+	//delete [] prim_rays; now dbg_renderer handles them
 
 	global_size = xsz * ysz;
+
+
+	init_dbg_renderer(xsz, ysz, scn, tex);
 	return true;
 }
 
@@ -143,7 +127,11 @@
 {
 	delete prog;
 
-	printf("rendertime mean: %ld msec\n", timing_sample_sum / num_timing_samples);
+	destroy_dbg_renderer();
+
+	if(num_timing_samples) {
+		printf("rendertime mean: %ld msec\n", timing_sample_sum / num_timing_samples);
+	}
 }
 
 bool render()
@@ -201,75 +189,6 @@
 	return true;
 }
 
-#define MIN(a, b)	((a) < (b) ? (a) : (b))
-static void dbg_set_gl_material(Material *mat)
-{
-	static Material def_mat = {{0.7, 0.7, 0.7, 1}, {0, 0, 0, 0}, 0, 0, 0};
-
-	if(!mat) mat = &def_mat;
-
-	glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat->kd);
-	glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, mat->ks);
-	glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, MIN(mat->spow, 128.0f));
-}
-
-void dbg_render_gl(Scene *scn, bool show_tree, bool show_obj)
-{
-	glPushAttrib(GL_ENABLE_BIT | GL_TRANSFORM_BIT | GL_LIGHTING_BIT);
-
-	for(int i=0; i<rinf.num_lights; i++) {
-		float lpos[4];
-
-		memcpy(lpos, lightlist[i].pos, sizeof lpos);
-		lpos[3] = 1.0;
-
-		glLightfv(GL_LIGHT0 + i, GL_POSITION, lpos);
-		glLightfv(GL_LIGHT0 + i, GL_DIFFUSE, lightlist[i].color);
-		glEnable(GL_LIGHT0 + i);
-	}
-
-	glDisable(GL_TEXTURE_2D);
-	glEnable(GL_DEPTH_TEST);
-	glEnable(GL_LIGHTING);
-
-	glMatrixMode(GL_PROJECTION);
-	glPushMatrix();
-	glLoadIdentity();
-	gluPerspective(45.0, (float)rinf.xsz / (float)rinf.ysz, 0.5, 1000.0);
-
-	if(show_obj) {
-		Material *materials = scn->get_materials();
-
-		int num_faces = scn->get_num_faces();
-		int cur_mat = -1;
-
-		for(int i=0; i<num_faces; i++) {
-			if(faces[i].matid != cur_mat) {
-				if(cur_mat != -1) {
-					glEnd();
-				}
-				dbg_set_gl_material(materials ? materials + faces[i].matid : 0);
-				cur_mat = faces[i].matid;
-				glBegin(GL_TRIANGLES);
-			}
-
-			for(int j=0; j<3; j++) {
-				glNormal3fv(faces[i].v[j].normal);
-				glVertex3fv(faces[i].v[j].pos);
-			}
-		}
-		glEnd();
-	}
-
-	if(show_tree) {
-		scn->draw_kdtree();
-	}
-
-	glPopMatrix();
-	glPopAttrib();
-
-	assert(glGetError() == GL_NO_ERROR);
-}
 
 void set_xform(float *matrix, float *invtrans)
 {
@@ -286,6 +205,12 @@
 	unmap_mem_buffer(mbuf_invtrans);
 }
 
+
+const RendInfo *get_render_info()
+{
+	return &rinf;
+}
+
 void set_render_option(int opt, bool val)
 {
 	switch(opt) {
@@ -403,6 +328,8 @@
 	return ray;
 }
 
+#define MIN(a, b)	((a) < (b) ? (a) : (b))
+
 static float *create_kdimage(const KDNodeGPU *kdtree, int num_nodes, int *xsz_ret, int *ysz_ret)
 {
 	int ysz = MIN(num_nodes, KDIMG_MAX_HEIGHT);
diff -r 54a96b738afe -r 6a30f27fa1e6 src/rt.cl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/rt.cl	Fri Sep 10 16:47:00 2010 +0100
@@ -0,0 +1,393 @@
+/* vim: set ft=opencl:ts=4:sw=4 */
+#include "common.h"
+
+struct RendInfo {
+	float4 ambient;
+	int xsz, ysz;
+	int num_faces, num_lights;
+	int max_iter;
+	int cast_shadows;
+};
+
+struct Vertex {
+	float4 pos;
+	float4 normal;
+	float4 tex;
+	float4 padding;
+};
+
+struct Face {
+	struct Vertex v[3];
+	float4 normal;
+	int matid;
+	int padding[3];
+};
+
+struct Material {
+	float4 kd, ks;
+	float kr, kt;
+	float spow;
+	float padding;
+};
+
+struct Light {
+	float4 pos, color;
+};
+
+struct Ray {
+	float4 origin, dir;
+};
+
+struct SurfPoint {
+	float t;
+	float4 pos, norm, dbg;
+	global const struct Face *obj;
+	struct Material mat;
+};
+
+struct Scene {
+	float4 ambient;
+	global const struct Face *faces;
+	int num_faces;
+	global const struct Light *lights;
+	int num_lights;
+	global const struct Material *matlib;
+	//global const struct KDNode *kdtree;
+	bool cast_shadows;
+};
+
+struct AABBox {
+	float4 min, max;
+};
+
+struct KDNode {
+	struct AABBox aabb;
+	int face_idx[MAX_NODE_FACES];
+	int num_faces;
+	int left, right;
+	int padding;
+};
+
+#define MIN_ENERGY	0.001
+
+float4 shade(struct Ray ray, struct Scene *scn, const struct SurfPoint *sp, read_only image2d_t kdimg);
+bool find_intersection(struct Ray ray, const struct Scene *scn, struct SurfPoint *sp, read_only image2d_t kdimg);
+bool intersect(struct Ray ray, global const struct Face *face, struct SurfPoint *sp);
+bool intersect_aabb(struct Ray ray, struct AABBox aabb);
+
+float4 reflect(float4 v, float4 n);
+float4 transform(float4 v, global const float *xform);
+void transform_ray(struct Ray *ray, global const float *xform, global const float *invtrans);
+float4 calc_bary(float4 pt, global const struct Face *face, float4 norm);
+float mean(float4 v);
+
+void read_kdnode(int idx, struct KDNode *node, read_only image2d_t kdimg);
+
+
+kernel void render(write_only image2d_t fb,
+		global const struct RendInfo *rinf,
+		global const struct Face *faces,
+		global const struct Material *matlib,
+		global const struct Light *lights,
+		global const struct Ray *primrays,
+		global const float *xform,
+		global const float *invtrans,
+		//global const struct KDNode *kdtree
+		read_only image2d_t kdtree_img)
+{
+	int idx = get_global_id(0);
+
+	struct Scene scn;
+	scn.ambient = rinf->ambient;
+	scn.faces = faces;
+	scn.num_faces = rinf->num_faces;
+	scn.lights = lights;
+	scn.num_lights = rinf->num_lights;
+	scn.matlib = matlib;
+	scn.cast_shadows = rinf->cast_shadows;
+
+	struct Ray ray = primrays[idx];
+	transform_ray(&ray, xform, invtrans);
+
+	float4 pixel = (float4)(0, 0, 0, 0);
+	float4 energy = (float4)(1.0, 1.0, 1.0, 0.0);
+	int iter = 0;
+
+	while(iter++ <= rinf->max_iter && mean(energy) > MIN_ENERGY) {
+		struct SurfPoint sp;
+		if(find_intersection(ray, &scn, &sp, kdtree_img)) {
+			pixel += shade(ray, &scn, &sp, kdtree_img) * energy;
+
+			float4 refl_col = sp.mat.ks * sp.mat.kr;
+
+			ray.origin = sp.pos;
+			ray.dir = reflect(-ray.dir, sp.norm);
+
+			energy *= refl_col;
+		} else {
+			energy = (float4)(0.0, 0.0, 0.0, 0.0);
+		}
+	}
+
+	int2 coord;
+	coord.x = idx % rinf->xsz;
+	coord.y = idx / rinf->xsz;
+
+	write_imagef(fb, coord, pixel);
+}
+
+float4 shade(struct Ray ray, struct Scene *scn, const struct SurfPoint *sp, read_only image2d_t kdimg)
+{
+	float4 norm = sp->norm;
+	//bool entering = true;
+
+	if(dot(ray.dir, norm) >= 0.0) {
+		norm = -norm;
+		//entering = false;
+	}
+
+	float4 dcol = scn->ambient * sp->mat.kd;
+	float4 scol = (float4)(0, 0, 0, 0);
+
+	for(int i=0; i<scn->num_lights; i++) {
+		float4 ldir = scn->lights[i].pos - sp->pos;
+
+		struct Ray shadowray;
+		shadowray.origin = sp->pos;
+		shadowray.dir = ldir;
+
+		if(!scn->cast_shadows || !find_intersection(shadowray, scn, 0, kdimg)) {
+			ldir = normalize(ldir);
+			float4 vdir = -ray.dir;
+			vdir.x = native_divide(vdir.x, RAY_MAG);
+			vdir.y = native_divide(vdir.y, RAY_MAG);
+			vdir.z = native_divide(vdir.z, RAY_MAG);
+			float4 vref = reflect(vdir, norm);
+
+			float diff = fmax(dot(ldir, norm), 0.0f);
+			dcol += sp->mat.kd /* scn->lights[i].color*/ * diff;
+
+			float spec = native_powr(fmax(dot(ldir, vref), 0.0f), sp->mat.spow);
+			scol += sp->mat.ks /* scn->lights[i].color*/ * spec;
+		}
+	}
+
+	return dcol + scol;
+}
+
+#define STACK_SIZE	MAX_TREE_DEPTH
+bool find_intersection(struct Ray ray, const struct Scene *scn, struct SurfPoint *spres, read_only image2d_t kdimg)
+{
+	struct SurfPoint sp0;
+	sp0.t = 1.0;
+	sp0.obj = 0;
+
+	int idxstack[STACK_SIZE];
+	int top = 0;			// points after the topmost element of the stack
+	idxstack[top++] = 0;	// root at tree[0]
+
+	while(top > 0) {
+		int idx = idxstack[--top];	// remove this index from the stack and process it
+
+		struct KDNode node;
+		read_kdnode(idx, &node, kdimg);
+
+		if(intersect_aabb(ray, node.aabb)) {
+			if(node.left == -1) {
+				// leaf node... check each face in turn and update the nearest intersection as needed
+				for(int i=0; i<node.num_faces; i++) {
+					struct SurfPoint spt;
+					int fidx = node.face_idx[i];
+
+					if(intersect(ray, scn->faces + fidx, &spt) && spt.t < sp0.t) {
+						sp0 = spt;
+					}
+				}
+			} else {
+				// internal node... recurse to the children
+				idxstack[top++] = node.left;
+				idxstack[top++] = node.right;
+			}
+		}
+	}
+
+	if(!sp0.obj) {
+		return false;
+	}
+
+	if(spres) {
+		*spres = sp0;
+		spres->mat = scn->matlib[sp0.obj->matid];
+	}
+	return true;
+}
+
+bool intersect(struct Ray ray, global const struct Face *face, struct SurfPoint *sp)
+{
+	float4 origin = ray.origin;
+	float4 dir = ray.dir;
+	float4 norm = face->normal;
+
+	float ndotdir = dot(dir, norm);
+
+	if(fabs(ndotdir) <= EPSILON) {
+		return false;
+	}
+
+	float4 pt = face->v[0].pos;
+	float4 vec = pt - origin;
+
+	float ndotvec = dot(norm, vec);
+	float t = native_divide(ndotvec, ndotdir);
+
+	if(t < EPSILON || t > 1.0) {
+		return false;
+	}
+	pt = origin + dir * t;
+
+
+	float4 bc = calc_bary(pt, face, norm);
+	float bc_sum = bc.x + bc.y + bc.z;
+
+	if(bc_sum < 1.0 - EPSILON || bc_sum > 1.0 + EPSILON) {
+		return false;
+	}
+
+	sp->t = t;
+	sp->pos = pt;
+	sp->norm = normalize(face->v[0].normal * bc.x + face->v[1].normal * bc.y + face->v[2].normal * bc.z);
+	sp->obj = face;
+	sp->dbg = bc;
+	return true;
+}
+
+bool intersect_aabb(struct Ray ray, struct AABBox aabb)
+{
+	if(ray.origin.x >= aabb.min.x && ray.origin.y >= aabb.min.y && ray.origin.z >= aabb.min.z &&
+			ray.origin.x < aabb.max.x && ray.origin.y < aabb.max.y && ray.origin.z < aabb.max.z) {
+		return true;
+	}
+
+	float4 bbox[2] = {
+		aabb.min.x, aabb.min.y, aabb.min.z, 0,
+		aabb.max.x, aabb.max.y, aabb.max.z, 0
+	};
+
+	int xsign = (int)(ray.dir.x < 0.0);
+	float invdirx = native_recip(ray.dir.x);
+	float tmin = (bbox[xsign].x - ray.origin.x) * invdirx;
+	float tmax = (bbox[1 - xsign].x - ray.origin.x) * invdirx;
+
+	int ysign = (int)(ray.dir.y < 0.0);
+	float invdiry = native_recip(ray.dir.y);
+	float tymin = (bbox[ysign].y - ray.origin.y) * invdiry;
+	float tymax = (bbox[1 - ysign].y - ray.origin.y) * invdiry;
+
+	if(tmin > tymax || tymin > tmax) {
+		return false;
+	}
+
+	if(tymin > tmin) tmin = tymin;
+	if(tymax < tmax) tmax = tymax;
+
+	int zsign = (int)(ray.dir.z < 0.0);
+	float invdirz = native_recip(ray.dir.z);
+	float tzmin = (bbox[zsign].z - ray.origin.z) * invdirz;
+	float tzmax = (bbox[1 - zsign].z - ray.origin.z) * invdirz;
+
+	if(tmin > tzmax || tzmin > tmax) {
+		return false;
+	}
+
+	return tmin < 1.0 && tmax > 0.0;
+}
+
+float4 reflect(float4 v, float4 n)
+{
+	return 2.0f * dot(v, n) * n - v;
+}
+
+float4 transform(float4 v, global const float *xform)
+{
+	float4 res;
+	res.x = v.x * xform[0] + v.y * xform[4] + v.z * xform[8] + xform[12];
+	res.y = v.x * xform[1] + v.y * xform[5] + v.z * xform[9] + xform[13];
+	res.z = v.x * xform[2] + v.y * xform[6] + v.z * xform[10] + xform[14];
+	res.w = 0.0;
+	return res;
+}
+
+void transform_ray(struct Ray *ray, global const float *xform, global const float *invtrans)
+{
+	ray->origin = transform(ray->origin, xform);
+	ray->dir = transform(ray->dir, invtrans);
+}
+
+float4 calc_bary(float4 pt, global const struct Face *face, float4 norm)
+{
+	float4 bc = (float4)(0, 0, 0, 0);
+
+	// calculate area of the whole triangle
+	float4 v1 = face->v[1].pos - face->v[0].pos;
+	float4 v2 = face->v[2].pos - face->v[0].pos;
+	float4 xv1v2 = cross(v1, v2);
+
+	float area = fabs(dot(xv1v2, norm)) * 0.5;
+	if(area < EPSILON) {
+		return bc;
+	}
+
+	float4 pv0 = face->v[0].pos - pt;
+	float4 pv1 = face->v[1].pos - pt;
+	float4 pv2 = face->v[2].pos - pt;
+
+	// calculate the area of each sub-triangle
+	float4 x12 = cross(pv1, pv2);
+	float4 x20 = cross(pv2, pv0);
+	float4 x01 = cross(pv0, pv1);
+
+	float a0 = fabs(dot(x12, norm)) * 0.5;
+	float a1 = fabs(dot(x20, norm)) * 0.5;
+	float a2 = fabs(dot(x01, norm)) * 0.5;
+
+	bc.x = native_divide(a0, area);
+	bc.y = native_divide(a1, area);
+	bc.z = native_divide(a2, area);
+	return bc;
+}
+
+float mean(float4 v)
+{
+	return native_divide(v.x + v.y + v.z, 3.0);
+}
+
+
+const sampler_t kdsampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_NONE | CLK_FILTER_NEAREST;
+
+// read a KD-tree node from a texture scanline
+void read_kdnode(int idx, struct KDNode *node, read_only image2d_t kdimg)
+{
+	int startx = KDIMG_NODE_WIDTH * (idx / KDIMG_MAX_HEIGHT);
+
+	int2 tc;
+	tc.x = startx;
+	tc.y = idx % KDIMG_MAX_HEIGHT;
+
+	node->aabb.min = read_imagef(kdimg, kdsampler, tc); tc.x++;
+	node->aabb.max = read_imagef(kdimg, kdsampler, tc);
+
+	tc.x = startx + 2 + MAX_NODE_FACES / 4;
+	float4 pix = read_imagef(kdimg, kdsampler, tc);
+	node->num_faces = (int)pix.x;
+	node->left = (int)pix.y;
+	node->right = (int)pix.z;
+
+	tc.x = startx + 2;
+	for(int i=0; i<node->num_faces; i+=4) {
+		float4 pix = read_imagef(kdimg, kdsampler, tc); tc.x++;
+		node->face_idx[i] = (int)pix.x;
+		node->face_idx[i + 1] = (int)pix.y;
+		node->face_idx[i + 2] = (int)pix.z;
+		node->face_idx[i + 3] = (int)pix.w;
+	}
+}
diff -r 54a96b738afe -r 6a30f27fa1e6 src/rt.h
--- a/src/rt.h	Sun Sep 05 16:43:55 2010 +0100
+++ b/src/rt.h	Fri Sep 10 16:47:00 2010 +0100
@@ -11,11 +11,26 @@
 	NUM_RENDER_OPTIONS
 };
 
+struct RendInfo {
+	float ambient[4];
+	int xsz, ysz;
+	int num_faces, num_lights;
+	int max_iter;
+	int cast_shadows;
+};
+
+struct Ray {
+	float origin[4], dir[4];
+};
+
+
 bool init_renderer(int xsz, int ysz, Scene *scn, unsigned int tex);
 void destroy_renderer();
 bool render();
 void set_xform(float *matrix, float *invtrans);
 
+const RendInfo *get_render_info();
+
 void set_render_option(int opt, bool val);
 void set_render_option(int opt, int val);
 void set_render_option(int opt, float val);
@@ -24,6 +39,13 @@
 int get_render_option_int(int opt);
 float get_render_option_float(int opt);
 
+// raytrace in the CPU
+bool init_dbg_renderer(int xsz, int ysz, Scene *scn, unsigned int texid);
+void destroy_dbg_renderer();
+void dbg_set_primary_rays(const Ray *rays);
+void dbg_render(const float *xform, const float *invtrans_xform, int num_threads = -1);
+
+// visualize the scene using OpenGL
 void dbg_render_gl(Scene *scn, bool show_tree = false, bool show_obj = true);
 
 #endif	/* RT_H_ */
diff -r 54a96b738afe -r 6a30f27fa1e6 src/scene.cc
--- a/src/scene.cc	Sun Sep 05 16:43:55 2010 +0100
+++ b/src/scene.cc	Fri Sep 10 16:47:00 2010 +0100
@@ -109,11 +109,27 @@
 	return true;
 }
 
+bool Scene::add_light(const Light &lt)
+{
+	try {
+		lights.push_back(lt);
+	}
+	catch(...) {
+		return false;
+	}
+	return true;
+}
+
 int Scene::get_num_meshes() const
 {
 	return (int)meshes.size();
 }
 
+int Scene::get_num_lights() const
+{
+	return (int)lights.size();
+}
+
 int Scene::get_num_faces() const
 {
 	if(num_faces >= 0) {
@@ -137,6 +153,38 @@
 	return kdtree_nodes(kdtree);
 }
 
+Mesh **Scene::get_meshes()
+{
+	if(meshes.empty()) {
+		return 0;
+	}
+	return &meshes[0];
+}
+
+const Mesh * const *Scene::get_meshes() const
+{
+	if(meshes.empty()) {
+		return 0;
+	}
+	return &meshes[0];
+}
+
+Light *Scene::get_lights()
+{
+	if(lights.empty()) {
+		return 0;
+	}
+	return &lights[0];
+}
+
+const Light *Scene::get_lights() const
+{
+	if(lights.empty()) {
+		return 0;
+	}
+	return &lights[0];
+}
+
 Material *Scene::get_materials()
 {
 	if(matlib.empty()) {
diff -r 54a96b738afe -r 6a30f27fa1e6 src/scene.h
--- a/src/scene.h	Sun Sep 05 16:43:55 2010 +0100
+++ b/src/scene.h	Fri Sep 10 16:47:00 2010 +0100
@@ -34,6 +34,10 @@
 	int matid;
 };
 
+struct Light {
+	float pos[4], color[4];
+};
+
 class AABBox {
 public:
 	float min[4], max[4];
@@ -70,6 +74,7 @@
 
 public:
 	std::vector<Mesh*> meshes;
+	std::vector<Light> lights;
 	std::vector<Material> matlib;
 	KDNode *kdtree;
 
@@ -77,11 +82,20 @@
 	~Scene();
 
 	bool add_mesh(Mesh *m);
+	bool add_light(const Light &lt);
+
 	int get_num_meshes() const;
+	int get_num_lights() const;
 	int get_num_faces() const;
 	int get_num_materials() const;
 	int get_num_kdnodes() const;
 
+	Mesh **get_meshes();
+	const Mesh * const *get_meshes() const;
+
+	Light *get_lights();
+	const Light *get_lights() const;
+
 	Material *get_materials();
 	const Material *get_materials() const;
 
diff -r 54a96b738afe -r 6a30f27fa1e6 src/vector.cc
--- a/src/vector.cc	Sun Sep 05 16:43:55 2010 +0100
+++ b/src/vector.cc	Fri Sep 10 16:47:00 2010 +0100
@@ -18,6 +18,13 @@
 	this->z = z;
 }
 
+Vector3::Vector3(const float *arr)
+{
+	x = arr[0];
+	y = arr[1];
+	z = arr[2];
+}
+
 void Vector3::normalize()
 {
 	float len = sqrt(x * x + y * y + z * z);
diff -r 54a96b738afe -r 6a30f27fa1e6 src/vector.h
--- a/src/vector.h	Sun Sep 05 16:43:55 2010 +0100
+++ b/src/vector.h	Fri Sep 10 16:47:00 2010 +0100
@@ -15,6 +15,7 @@
 
 	Vector3();
 	Vector3(float x, float y, float z);
+	Vector3(const float *arr);
 
 	void normalize();
 	inline float length();
@@ -28,10 +29,13 @@
 
 inline Vector3 operator -(const Vector3 &vec);
 inline Vector3 operator *(const Vector3 &vec, float s);
+inline Vector3 operator /(const Vector3 &vec, float s);
 
 inline float dot(const Vector3 &a, const Vector3 &b);
 inline Vector3 cross(const Vector3 &a, const Vector3 &b);
 
+inline Vector3 reflect(const Vector3 &v, const Vector3 &n);
+
 #include "vector.inl"
 
 #endif	/* VECTOR_H_ */
diff -r 54a96b738afe -r 6a30f27fa1e6 src/vector.inl
--- a/src/vector.inl	Sun Sep 05 16:43:55 2010 +0100
+++ b/src/vector.inl	Fri Sep 10 16:47:00 2010 +0100
@@ -41,6 +41,10 @@
 	return Vector3(vec.x * s, vec.y * s, vec.z * s);
 }
 
+inline Vector3 operator /(const Vector3 &vec, float s)
+{
+	return Vector3(vec.x / s, vec.y / s, vec.z / s);
+}
 
 inline float dot(const Vector3 &a, const Vector3 &b)
 {
@@ -51,3 +55,8 @@
 {
 	return Vector3(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x);
 }
+
+inline Vector3 reflect(const Vector3 &v, const Vector3 &n)
+{
+	return n * (2.0 * dot(v, n)) - v;
+}