# HG changeset patch
# User Radek Brich <radek.brich@devl.cz>
# Date 1196379891 -3600
# Node ID 79b516a3803d92a01b52850b2cd124a5ae766be5
# Parent  f22952603f2959a47311e59818108c57a6f21f59
naive color driven sub-sampling
slightly optimized KdTree::nearest_intersection
fixed bug in Box::intersect
fixed stripes on spheres in spheres_ao.py (was caused by AO distance)
new KdTree property: max_depth
minor changes in realtime.py

diff -r f22952603f29 -r 79b516a3803d TODO
--- a/TODO	Thu Nov 29 18:30:16 2007 +0100
+++ b/TODO	Fri Nov 30 00:44:51 2007 +0100
@@ -5,7 +5,7 @@
  * update Python binding, Camera, new classes
  * more complex demos
  * check/update triangle routines
-
+ * namespace
 
 New Classes?
 ============
diff -r f22952603f29 -r 79b516a3803d ccdemos/realtime.cc
--- a/ccdemos/realtime.cc	Thu Nov 29 18:30:16 2007 +0100
+++ b/ccdemos/realtime.cc	Fri Nov 30 00:44:51 2007 +0100
@@ -2,8 +2,8 @@
 
 #include "raytracer.h"
 
-int w = 512;
-int h = 200;
+int w = 480;
+int h = 288;
 float *render_buffer;
 
 Raytracer rt;
@@ -13,11 +13,9 @@
 {
 	rt.render(w, h, render_buffer);
 
-	if ( SDL_MUSTLOCK(screen) ) {
-		if ( SDL_LockSurface(screen) < 0 ) {
+	if (SDL_MUSTLOCK(screen))
+		if (SDL_LockSurface(screen) < 0)
 			return;
-		}
-	}
 
 	Uint32 *bufp = (Uint32 *)screen->pixels;
 	unsigned char c[3];
@@ -34,12 +32,13 @@
 		bufp++;
 	}
 
-	if ( SDL_MUSTLOCK(screen) ) {
+	if (SDL_MUSTLOCK(screen))
 		SDL_UnlockSurface(screen);
-	}
 
-	SDL_UpdateRect(screen, 0, 0, w, h);
-	SDL_Flip(screen);
+	if (screen->flags & SDL_DOUBLEBUF)
+		SDL_Flip(screen);
+	else
+		SDL_UpdateRect(screen, 0, 0, w, h);
 }
 
 int main()
@@ -64,25 +63,29 @@
 	render_buffer = (float *) malloc(w*h*3*sizeof(float));
 
 	rt.setThreads(2);
+	rt.setMaxDepth(1);
 
 	KdTree top;
 	rt.setTop(&top);
 
 	Light light1(Vector3(2.0, -5.0, -5.0), Colour(0.7, 0.3, 0.6));
+	light1.castShadows(false);
 	rt.addlight(&light1);
 
 	Light light2(Vector3(-2.0, 10.0, 2.0), Colour(0.4, 0.6, 0.3));
+	light2.castShadows(false);
 	rt.addlight(&light2);
 
 	Material mat_sph(Colour(1.0, 1.0, 1.0));
-	for (int y=0; y<20; y++)
-		for (int x=0; x<20; x++)
-			rt.addshape(new Sphere(Vector3(x-10, (float)random()/RAND_MAX*5.0, y-10), 0.45, &mat_sph));
+	for (int y=0; y<10; y++)
+		for (int x=0; x<10; x++)
+			rt.addshape(new Sphere(Vector3(x*2-10, (float)random()/RAND_MAX*5.0, y*2-10), 0.45, &mat_sph));
 
 	rt.setCamera(&cam);
 	cam.setEye(Vector3(0,0,10));
 
 	/* build kd-tree */
+	top.setMaxDepth(100);
 	top.optimize();
 
 	/* loop... */
diff -r f22952603f29 -r 79b516a3803d demos/spheres_ao.py
--- a/demos/spheres_ao.py	Thu Nov 29 18:30:16 2007 +0100
+++ b/demos/spheres_ao.py	Fri Nov 30 00:44:51 2007 +0100
@@ -7,7 +7,7 @@
 import Image
 
 rt = Raytracer()
-rt.ambientocclusion(samples=60, distance=6.0, angle=0.5)
+rt.ambientocclusion(samples=100, distance=16.0, angle=0.5)
 
 light1 = Light(position=(0.0, 5.0, -5.0), colour=(0.7, 0.3, 0.6))
 light1.castshadows(False)
diff -r f22952603f29 -r 79b516a3803d src/kdtree.cc
--- a/src/kdtree.cc	Thu Nov 29 18:30:16 2007 +0100
+++ b/src/kdtree.cc	Fri Nov 30 00:44:51 2007 +0100
@@ -50,11 +50,14 @@
 };
 
 // stack element for kd-tree traversal
-struct StackElem
+class StackElem
 {
+public:
 	KdNode* node; /* pointer to far child */
 	float t; /* the entry/exit signed distance */
 	Vector3 pb; /* the coordinates of entry/exit point */
+	StackElem(KdNode *anode, const float &at, const Vector3 &apb):
+		node(anode), t(at), pb(apb) {};
 };
 
 // ----------------------------------------
@@ -96,9 +99,9 @@
 		delete[] children;
 }
 
-void KdNode::subdivide(BBox bbox, int depth)
+void KdNode::subdivide(BBox bbox, int maxdepth)
 {
-	if (depth >= 20 || shapes->size() <= 4)
+	if (maxdepth <= 0 || shapes->size() <= 2)
 	{
 		setLeaf();
 		return;
@@ -295,8 +298,8 @@
 	lbb.H.cell[axis] = split;
 	rbb.L.cell[axis] = split;
 
-	children[0].subdivide(lbb, depth+1);
-	children[1].subdivide(rbb, depth+1);
+	children[0].subdivide(lbb, maxdepth-1);
+	children[1].subdivide(rbb, maxdepth-1);
 }
 
 void KdTree::build()
@@ -306,7 +309,7 @@
 	ShapeList::iterator shape;
 	for (shape = shapes.begin(); shape != shapes.end(); shape++)
 		root->addShape(*shape);
-	root->subdivide(bbox, 0);
+	root->subdivide(bbox, max_depth);
 	built = true;
 }
 
@@ -328,29 +331,23 @@
 	KdNode *farchild, *node;
 	node = root; /* start from the kd-tree root node */
 
-	stack<StackElem*> st;
-
-	StackElem *enPt = new StackElem();
-	enPt->t = a; /* set the signed distance */
-	enPt->node = NULL;
+	/* std vector is much faster than stack */
+	vector<StackElem*> st;
 
-	/* distinguish between internal and external origin */
-	if (a >= 0.0) /* a ray with external origin */
-		enPt->pb = ray.o + ray.dir * a;
-	else /* a ray with internal origin */
-		enPt->pb = ray.o;
+	StackElem *enPt = new StackElem(NULL, a,
+		/* distinguish between internal and external origin of a ray*/
+		a >= 0.0 ?
+			ray.o + ray.dir * a : /* external */
+			ray.o);               /* internal */
 
 	/* setup initial exit point in the stack */
-	StackElem *exPt = new StackElem();
-	exPt->t = b;
-	exPt->pb = ray.o + ray.dir * b;
-	exPt->node = NULL; /* set termination flag */
-	st.push(exPt);
+	StackElem *exPt = new StackElem(NULL, b, ray.o + ray.dir * b);
+	st.push_back(exPt);
 
 	/* loop, traverse through the whole kd-tree, until an object is intersected or ray leaves the scene */
 	while (node)
 	{
-		exPt = st.top();
+		exPt = st.back();
 		/* loop until a leaf is found */
 		while (!node->isLeaf())
 		{
@@ -393,13 +390,11 @@
 			t = (splitVal - ray.o.cell[axis]) / ray.dir.cell[axis];
 
 			/* setup the new exit point and push it onto stack */
-			exPt = new StackElem();
-			exPt->t = t;
-			exPt->node = farchild;
+			exPt = new StackElem(farchild, t, Vector3());
 			exPt->pb[axis] = splitVal;
 			exPt->pb[(axis+1)%3] = ray.o.cell[(axis+1)%3] + t * ray.dir.cell[(axis+1)%3];
 			exPt->pb[(axis+2)%3] = ray.o.cell[(axis+2)%3] + t * ray.dir.cell[(axis+2)%3];
-			st.push(exPt);
+			st.push_back(exPt);
 		} /* while */
 
 		/* current node is the leaf . . . empty or full */
@@ -419,8 +414,8 @@
 		{
 			while (!st.empty())
 			{
-				delete st.top();
-				st.pop();
+				delete st.back();
+				st.pop_back();
 			}
 			nearest_distance = dist;
 			return nearest_shape;
@@ -430,7 +425,7 @@
 
 		/* retrieve the pointer to the next node, it is possible that ray traversal terminates */
 		node = enPt->node;
-		st.pop();
+		st.pop_back();
 	} /* while */
 
 	delete enPt;
diff -r f22952603f29 -r 79b516a3803d src/kdtree.h
--- a/src/kdtree.h	Thu Nov 29 18:30:16 2007 +0100
+++ b/src/kdtree.h	Fri Nov 30 00:44:51 2007 +0100
@@ -62,8 +62,9 @@
 {
 	KdNode *root;
 	bool built;
+	int max_depth;
 public:
-	KdTree() : Container(), root(NULL), built(false) {};
+	KdTree() : Container(), root(NULL), built(false), max_depth(32) {};
 	~KdTree() { if (root) delete root; };
 	void addShape(Shape* aShape) { Container::addShape(aShape); built = false; };
 	Shape *nearest_intersection(const Shape *origin_shape, const Ray &ray,
@@ -72,6 +73,7 @@
 	void build();
 	void save(ostream &str, KdNode *node = NULL);
 	void load(istream &str, KdNode *node = NULL);
+	void setMaxDepth(int md) { max_depth = md; };
 };
 
 #endif
diff -r f22952603f29 -r 79b516a3803d src/raytracer.cc
--- a/src/raytracer.cc	Thu Nov 29 18:30:16 2007 +0100
+++ b/src/raytracer.cc	Fri Nov 30 00:44:51 2007 +0100
@@ -109,7 +109,7 @@
 			float L_dot_N = dot(L, normal);
 			if (L_dot_N > 0) {
 				// test if this light is occluded (sharp shadows)
-				if ((*light)->shadows) {
+				if ((*light)->cast_shadows) {
 					Ray shadow_ray = Ray(P, L);
 					float dist = FLT_MAX;
 					if (top->nearest_intersection(nearest_shape, shadow_ray, dist))
@@ -124,9 +124,8 @@
 		}
 
 		// reflection
-		int trace_max_depth = 4;
 		Vector3 newdir = ray.dir - 2.0 * dot(ray.dir, normal) * normal;
-		if (depth < trace_max_depth && nearest_shape->material->reflection > 0.01) {
+		if (depth < max_depth && nearest_shape->material->reflection > 0.01) {
 			Ray newray = Ray(P, newdir);
 			Colour refl_col = raytrace(newray, depth + 1, nearest_shape);
 			acc += nearest_shape->material->reflection * refl_col;
@@ -160,8 +159,11 @@
 static void *renderrow(void *data)
 {
 	RenderrowData *d = (RenderrowData*) data;
+	int subsample = d->rt->getSubsample();
+	float subsample2 = 1.0/(subsample*subsample);
+	int ww = d->w*3;
 	Vector3 dir = d->dfix;
-	for (int x = 0; x < d->w; x++) {
+	for (int x = 0; x < d->w; x += subsample) {
 		// generate a ray from eye passing through this pixel
 #if OVERSAMPLING
 		// 5x oversampling
@@ -182,11 +184,117 @@
 		dir.normalize();
 		Ray ray(d->eye, dir);
 		Colour c = d->rt->raytrace(ray, 0, NULL);
+		if (subsample > 1)
+		{
+			Colour ic;
+			// top-left is 'c'
+			// top-right
+			Vector3 tmpdir = dir + (subsample-1)*d->dx;
+			tmpdir.normalize();
+			Ray ray(d->eye, tmpdir);
+			Colour c2 = d->rt->raytrace(ray, 0, NULL);
+			// bottom right
+			tmpdir += (subsample-1)*d->dy;
+			tmpdir.normalize();
+			ray.dir = tmpdir;
+			Colour c4 = d->rt->raytrace(ray, 0, NULL);
+			// bottom left
+			tmpdir = dir + (subsample-1)*d->dy;
+			tmpdir.normalize();
+			ray.dir = tmpdir;
+			Colour c3 = d->rt->raytrace(ray, 0, NULL);
+			// are the colors similar?
+			float m = (c-c2).mag2();
+			m = max(m, (c2-c3).mag2());
+			m = max(m, (c3-c4).mag2());
+			m = max(m, (c4-c).mag2());
+			if (m < 0.001)
+			{
+				// interpolate
+				float *i = d->iter;
+				for (int x = 0; x < subsample; x++)
+				{
+					for (int y = 0; y < subsample; y++)
+					{
+						ic = c*(subsample-x)*(subsample-y)*subsample2
+							+ c2*(x)*(subsample-y)*subsample2
+							+ c3*(subsample-x)*(y)*subsample2
+							+ c4*(x)*(y)*subsample2;
+						*(i + ww*y) = ic.r;
+						*(i + ww*y + 1) = ic.g;
+						*(i + ww*y + 2) = ic.b;
+					}
+					i += 3;
+				}
+				d->iter = i;
+			}
+			else
+			{
+				// render
+				Vector3 tmpdir = dir;
+				// first column
+				*(d->iter) = c.r;
+				*(d->iter + 1) = c.g;
+				*(d->iter + 2) = c.b;
+				for (int y = 1; y < subsample-1; y++)
+				{
+					Vector3 tmp2dir = tmpdir + y*d->dy;
+					tmp2dir.normalize();
+					ray.dir = tmp2dir;
+					ic = d->rt->raytrace(ray, 0, NULL);
+					*(d->iter + ww*y) = ic.r;
+					*(d->iter + ww*y + 1) = ic.g;
+					*(d->iter + ww*y + 2) = ic.b;
+				}
+				*(d->iter + ww*(subsample-1)) = c3.r;
+				*(d->iter + ww*(subsample-1) + 1) = c3.g;
+				*(d->iter + ww*(subsample-1) + 2) = c3.b;
+				d->iter += 3;
+				tmpdir += d->dx;
+				// middle
+				for (int x = 1; x < subsample-1; x++)
+				{
+					for (int y = 0; y < subsample; y++)
+					{
+						Vector3 tmp2dir = tmpdir + y*d->dy;
+						tmp2dir.normalize();
+						ray.dir = tmp2dir;
+						ic = d->rt->raytrace(ray, 0, NULL);
+						*(d->iter + ww*y) = ic.r;
+						*(d->iter + ww*y + 1) = ic.g;
+						*(d->iter + ww*y + 2) = ic.b;
+					}
+					d->iter += 3;
+					tmpdir += d->dx;
+				}
+				// last column
+				*(d->iter) = c2.r;
+				*(d->iter + 1) = c2.g;
+				*(d->iter + 2) = c2.b;
+				for (int y = 1; y < subsample-1; y++)
+				{
+					Vector3 tmp2dir = tmpdir + y*d->dy;
+					tmp2dir.normalize();
+					ray.dir = tmp2dir;
+					ic = d->rt->raytrace(ray, 0, NULL);
+					*(d->iter + ww*y) = ic.r;
+					*(d->iter + ww*y + 1) = ic.g;
+					*(d->iter + ww*y + 2) = ic.b;
+				}
+				*(d->iter + ww*(subsample-1)) = c4.r;
+				*(d->iter + ww*(subsample-1) + 1) = c4.g;
+				*(d->iter + ww*(subsample-1) + 2) = c4.b;
+				d->iter += 3;
+			}
+		}
 #endif
-		*(d->iter++) = c.r;
-		*d->iter++ = c.g;
-		*d->iter++ = c.b;
-		dir += d->dx;
+		if (subsample <= 1)
+		{
+			*d->iter++ = c.r;
+			*d->iter++ = c.g;
+			*d->iter++ = c.b;
+		}
+		dir += d->dx*subsample;
 	}
 #ifdef PTHREADS
 	pthread_exit((void *)d);
@@ -217,7 +325,7 @@
 
 	/* for each pixel... */
 	infomsg("* rendering row   0 (  0%% done)");
-	for (int y = 0; y < h; y++)
+	for (int y = 0; y < h; y += subsample)
 	{
 		d = (RenderrowData*) malloc(sizeof(RenderrowData));
 		d->rt = this;
@@ -225,9 +333,7 @@
 		d->eye = camera->eye;
 		d->dfix = dfix;
 		d->dx = dx;
-#if OVERSAMPLING
 		d->dy = dy;
-#endif
 		d->iter = buffer + y*3*w;
 #ifdef PTHREADS
 		/* create new thread and increase 't' */
@@ -248,7 +354,7 @@
 		renderrow((void *)d);
 		free(d);
 #endif
-		dfix += dy;
+		dfix += dy*subsample;
 		infomsg("\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b%4d (%2d%% done)", y, y*100/(h-1));
 	}
 	infomsg("\n");
diff -r f22952603f29 -r 79b516a3803d src/raytracer.h
--- a/src/raytracer.h	Thu Nov 29 18:30:16 2007 +0100
+++ b/src/raytracer.h	Fri Nov 30 00:44:51 2007 +0100
@@ -20,10 +20,7 @@
 struct RenderrowData {
 	Raytracer *rt;
 	int w;
-	Vector3 eye, dfix, dx;
-#if OVERSAMPLING
-	Vector3 dy;
-#endif
+	Vector3 eye, dfix, dx, dy;
 	float *iter;
 };
 
@@ -36,11 +33,13 @@
 	int ao_samples;
 	float ao_distance, ao_angle;
 	int num_threads;
+	int subsample;
+	int max_depth;
 
 	Vector3 SphereDistribute(int i, int n, float extent, Vector3 &normal);
 public:
 	Raytracer(): top(NULL), camera(NULL), lights(), bg_colour(0.0, 0.0, 0.0),
-		ao_samples(0), num_threads(4) {};
+		ao_samples(0), num_threads(4), subsample(8), max_depth(4) {};
 	void render(int w, int h, float *buffer);
 	Colour raytrace(Ray &ray, int depth, Shape *origin_shape);
 	void addshape(Shape *shape) { top->addShape(shape); };
@@ -48,6 +47,8 @@
 	void setCamera(Camera *cam) { camera = cam; };
 	void setTop(Container *atop) { top = atop; };
 	Container *getTop() { return top; };
+	int getSubsample() { return subsample; };
+	void setMaxDepth(int newdepth) { max_depth = newdepth; };
 
 	void ambientocclusion(int samples, float distance, float angle);
 	void setThreads(int num) { num_threads = num; };
diff -r f22952603f29 -r 79b516a3803d src/raytracermodule.cc
--- a/src/raytracermodule.cc	Thu Nov 29 18:30:16 2007 +0100
+++ b/src/raytracermodule.cc	Fri Nov 30 00:44:51 2007 +0100
@@ -88,7 +88,7 @@
 	if (!PyArg_ParseTuple(args, "i", &shadows))
 		return NULL;
 
-	((LightObject *)self)->light->castshadows(shadows);
+	((LightObject *)self)->light->castShadows(shadows);
 
 	Py_INCREF(Py_None);
 	return Py_None;
diff -r f22952603f29 -r 79b516a3803d src/scene.cc
--- a/src/scene.cc	Thu Nov 29 18:30:16 2007 +0100
+++ b/src/scene.cc	Fri Nov 30 00:44:51 2007 +0100
@@ -12,6 +12,8 @@
 
 void Camera::rotate(const Quaternion &q)
 {
+	/*
+	//non-optimized
 	Quaternion res;
 	res = q * Quaternion(u) * conjugate(q);
 	u = res.toVector();
@@ -19,10 +21,9 @@
 	v = res.toVector();
 	res = q * Quaternion(p) * conjugate(q);
 	p = res.toVector();
-	p.normalize();
-	u.normalize();
-	v.normalize();
-/*  // optimized version
+	*/
+
+	// optimized
 	float t2 =   q.a*q.b;
 	float t3 =   q.a*q.c;
 	float t4 =   q.a*q.d;
@@ -44,7 +45,10 @@
 	x = 2*( (t8 + t10)*v.x + (t6 -  t4)*v.y + (t3 + t7)*v.z ) + v.x;
 	y = 2*( (t4 +  t6)*v.x + (t5 + t10)*v.y + (t9 - t2)*v.z ) + v.y;
 	z = 2*( (t7 -  t3)*v.x + (t2 +  t9)*v.y + (t5 + t8)*v.z ) + v.z;
-	v = Vector3(x,y,z);*/
+	v = Vector3(x,y,z);
+	p.normalize();
+	u.normalize();
+	v.normalize();
 }
 
 void Camera::move(const float fw, const float left, const float up)
@@ -92,28 +96,15 @@
 
 bool Sphere::intersect(const Ray &ray, float &dist)
 {
-	Vector3 V = ((Ray)ray).o - center;
-
-	float Vd = - dot(V, ray.dir);
-	float Det = Vd * Vd - (dot(V,V) - sqr_radius);
-
+	Vector3 V = ray.o - center;
+	register float d = -dot(V, ray.dir);
+	register float Det = d * d - (dot(V,V) - sqr_radius);
 	if (Det > 0) {
-		Det = sqrtf(Det);
-		float t1 = Vd - Det;
-		if (t1 > 0)
+		d -= sqrtf(Det);
+		if (d > 0 && d < dist)
 		{
-			if (t1 < dist) {
-				dist = t1;
-				return true;
-			}
-		} else {
-			float t2 = Vd + Det;
-			if (t2 > 0)
-			{
-				// ray from inside of the sphere
-				dist = t2;
-				return true;
-			}
+			dist = d;
+			return true;
 		}
 	}
 	return false;
@@ -163,8 +154,15 @@
 
 bool Box::intersect(const Ray &ray, float &dist)
 {
-	float b;
-	return get_bbox().intersect(ray, dist, b);
+	float a,b;
+	bool res = get_bbox().intersect(ray, a, b);
+	if (res && a < dist)
+	{
+		dist = a;
+		return true;
+	}
+	else
+		return false;
 }
 
 Vector3 Box::normal(Vector3 &P)
diff -r f22952603f29 -r 79b516a3803d src/scene.h
--- a/src/scene.h	Thu Nov 29 18:30:16 2007 +0100
+++ b/src/scene.h	Fri Nov 30 00:44:51 2007 +0100
@@ -92,11 +92,11 @@
 public:
 	Vector3 pos;
 	Colour colour;
-	int shadows;
+	bool cast_shadows;
 
 	Light(const Vector3 &position, const Colour &acolour):
-		pos(position), colour(acolour), shadows(1) {};
-	void castshadows(bool ashadows) { shadows = ashadows; };
+		pos(position), colour(acolour), cast_shadows(true) {};
+	void castShadows(bool cast) { cast_shadows = cast; };
 };
 
 class Texture
diff -r f22952603f29 -r 79b516a3803d src/vector.h
--- a/src/vector.h	Thu Nov 29 18:30:16 2007 +0100
+++ b/src/vector.h	Fri Nov 30 00:44:51 2007 +0100
@@ -34,6 +34,8 @@
 	// index operator
 	float &operator[](int index)	{ return cell[index]; };
 
+	bool operator==(Vector3 &v) { return x==v.x && y==v.y && z==v.z; };
+
 	// normalize
 	Vector3 normalize()
 	{