--- a/.bzrignore	Fri Nov 30 00:44:51 2007 +0100
+++ b/.bzrignore	Mon Dec 03 01:49:23 2007 +0100
@@ -1,6 +1,8 @@
+demos/ModulePath
 demos/*.png
 demos/kdtree.obj
 demos/*.ply
 ccdemos/*.png
 ccdemos/spheres_shadow
 ccdemos/realtime
+bin/*
\ No newline at end of file

--- a/Makefile	Fri Nov 30 00:44:51 2007 +0100
+++ b/Makefile	Mon Dec 03 01:49:23 2007 +0100
@@ -1,39 +1,38 @@
-CCFLAGS=-I./src -Wall -Wno-write-strings -fno-strict-aliasing -DPTHREADS
-LDFLAGS=
+ROOT=$(shell pwd)
+include config.mk
+
+all: python-module demos ccdemos
+
+python-module: libs-float
+	$(MAKE) -C src python-module
+
+demos: python-module
+	$(MAKE) -C demos
 
-ifeq ($(OS), Windows_NT)
-  CCFLAGS+=-I"C:/Program Files/Python25/include"
-  LDFLAGS+=-L"C:\Program Files\Python25\libs" -lpython25 -lpthreadGC2
-  MODULENAME=raytracer.pyd
-else
-  CCFLAGS+=-pthread -fPIC `python-config --includes`
-  MODULENAME=raytracermodule.so
-endif
+ccdemos: libs-float libs-double
+	$(MAKE) -C ccdemos
+
+libs-float:
+	$(MAKE) -C src libs-float
 
-# optimizations
-#CCFLAGS+=-g -O0
-CCFLAGS+=-O3 -pipe -fomit-frame-pointer -ffast-math -msse3
+libs-double:
+	$(MAKE) -C src libs-double
+
+clean:
+	$(MAKE) -C src clean
+	$(MAKE) -C demos clean
+	$(MAKE) -C ccdemos clean
 
 
 # TARGETS
 #########
 
-all: python-module
-
-python-module: $(MODULENAME)
-
 tests: testvector testmatrix
 
-clean:
-	rm -f *.o $(MODULENAME)
-
 
 # RULES
 #######
 
-%.o: src/%.cc
-	$(CXX) -c -o $@ src/$*.cc $(CCFLAGS)
-
 test%: tests/%.cc
 	$(CXX) -o $@ tests/$*.cc $(CCFLAGS)
 	./$@
@@ -42,18 +41,6 @@
 # DEPENDENCIES
 ##############
 
-# C++ raytracer
-vector.o: src/vector.cc src/vector.h
-matrix.o: src/matrix.cc src/matrix.h src/vector.h
-noise.o: src/noise.cc src/noise.h
-scene.o: src/scene.cc src/scene.h src/vector.h src/noise.h src/common.h
-kdtree.o: src/kdtree.cc src/kdtree.h src/scene.h
-raytracer.o: src/raytracer.cc src/raytracer.h src/scene.h src/vector.h src/noise.h
-
-# python module
-raytracermodule.o: src/raytracermodule.cc src/raytracer.h src/scene.h src/vector.h
-$(MODULENAME): raytracermodule.o raytracer.o scene.o noise.o kdtree.o
-	$(CXX) $^ -shared -o $@ $(LDFLAGS)
 
 # library tests
 testvector: tests/vector.cc src/vector.h

--- a/ccdemos/Makefile	Fri Nov 30 00:44:51 2007 +0100
+++ b/ccdemos/Makefile	Mon Dec 03 01:49:23 2007 +0100
@@ -1,28 +1,33 @@
-CCFLAGS=-g -O3 -I../src
-LDFLAGS=-L.. -pthread
-RGBLIB_LDFLAGS=$(LDFLAGS) -lpng
-SDL_CCFLAGS=$(CCFLAGS) $(shell sdl-config --cflags)
-SDL_LDFLAGS=$(LDFLAGS) $(shell sdl-config --libs)
-PYRIT_OBJS=$(shell ls ../*.o | grep -v raytracermodule)
+ifndef $(ROOT)
+	ROOT=$(shell pwd)/..
+endif
+
+include $(ROOT)/config.mk
+
+
+### Targets ###
+all: realtime spheres_shadow
+
+realtime: realtime.o libs-double
+	$(CXX) -o $@ $(ROOT)/bin/libs-double/*.o $< $(LDFLAGS) $(SDL_LDFLAGS)
 
-all: spheres_shadow realtime
+spheres_shadow: spheres_shadow.o image.o libs-float
+	$(CXX) -o $@ $(ROOT)/bin/libs-float/*.o $< image.o $(LDFLAGS) -lpng
+
+libs-float:
+	$(MAKE) -C ../src libs-float
 
-%.o: %.c
+libs-double:
+	$(MAKE) -C ../src libs-double
+
+realtime.o: realtime.cc
+	$(CXX) -c -o $@ $(CCFLAGS) $(SDL_CCFLAGS) $< $(DEFS) -DPYRIT_DOUBLE
+
+image.o: image.c
 	$(CXX) -c -o $@ $*.c
 
-%.o: %.cc
-	$(CXX) -c -o $@ $*.cc $(CCFLAGS)
-
-%: %.o
-	(cd .. && make)
-	$(CXX) -o $@ $(PYRIT_OBJS) $^ $(RGBLIB_LDFLAGS)
-
-image.o: image.c
 spheres_shadow.o: spheres_shadow.cc
-spheres_shadow: spheres_shadow.o image.o
-
-realtime: realtime.cc
-	$(CXX) -o $@ $@.cc $(SDL_CCFLAGS) -L.. $(PYRIT_OBJS) $(SDL_LDFLAGS)
+	$(CXX) -c -o $@ $*.cc $(CCFLAGS) $(DEFS)
 
 clean:
 	rm -f spheres_shadow realtime *.o

--- a/ccdemos/realtime.cc	Fri Nov 30 00:44:51 2007 +0100
+++ b/ccdemos/realtime.cc	Mon Dec 03 01:49:23 2007 +0100
@@ -4,7 +4,7 @@
 
 int w = 480;
 int h = 288;
-float *render_buffer;
+Float *render_buffer;
 
 Raytracer rt;
 Camera cam;
@@ -19,7 +19,7 @@
 
 	Uint32 *bufp = (Uint32 *)screen->pixels;
 	unsigned char c[3];
-	for (float *fd = render_buffer; fd != render_buffer + w*h*3; fd += 3)
+	for (Float *fd = render_buffer; fd != render_buffer + w*h*3; fd += 3)
 	{
 		for (int i = 0; i < 3; i++)
 		{
@@ -60,10 +60,10 @@
 	}
 
 	/* initialize raytracer and prepare scene */
-	render_buffer = (float *) malloc(w*h*3*sizeof(float));
+	render_buffer = (Float *) malloc(w*h*3*sizeof(Float));
 
 	rt.setThreads(2);
-	rt.setMaxDepth(1);
+	rt.setMaxDepth(3);
 
 	KdTree top;
 	rt.setTop(&top);
@@ -79,7 +79,7 @@
 	Material mat_sph(Colour(1.0, 1.0, 1.0));
 	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.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));
@@ -91,7 +91,7 @@
 	/* loop... */
 	SDL_Event event;
 	bool quit = false;
-	float roty = 0.0, rotx = 0.0, move = 0.0;
+	Float roty = 0.0, rotx = 0.0, move = 0.0;
 	while (!quit)
 	{
 		while (SDL_PollEvent(&event))
@@ -100,7 +100,7 @@
 				case SDL_VIDEORESIZE:
 					w = event.resize.w;
 					h = event.resize.h;
-					render_buffer = (float *) realloc(render_buffer, w*h*3*sizeof(float));
+					render_buffer = (Float *) realloc(render_buffer, w*h*3*sizeof(Float));
 					screen = SDL_SetVideoMode(w, h, 32, SDL_HWSURFACE|SDL_DOUBLEBUF|SDL_RESIZABLE);
 					break;
 				case SDL_KEYDOWN:

--- a/ccdemos/spheres_shadow.cc	Fri Nov 30 00:44:51 2007 +0100
+++ b/ccdemos/spheres_shadow.cc	Mon Dec 03 01:49:23 2007 +0100
@@ -4,7 +4,10 @@
 int main()
 {
 	Raytracer rt;
-	rt.setThreads(1);
+	rt.setThreads(2);
+
+	KdTree top;
+	rt.setTop(&top);
 
 	Light light1(Vector3(0.0, 5.0, -5.0), Colour(0.7, 0.3, 0.6));
 	rt.addlight(&light1);
@@ -29,19 +32,21 @@
 	Sphere tinysphere(Vector3(-1.2, 0.0, -2.0), 0.5, &mat3);
 	rt.addshape(&tinysphere);
 
+	top.optimize();
+
 	Camera cam;
 	cam.setEye(Vector3(0,0,15));
 	rt.setCamera(&cam);
 
 	int w = 800;
 	int h = 600;
-	float *fdata = (float *) malloc(w*h*3*sizeof(float));
+	Float *fdata = (Float *) malloc(w*h*3*sizeof(Float));
 	rt.render(w, h, fdata);
 
 	struct image *img;
 	new_image(&img, w, h, 3);
 
-	float *fd = fdata;
+	Float *fd = fdata;
 	for (char *cd = img->data; cd != img->data + w*h*3; cd++, fd++) {
 		if (*fd > 1.0)
 			*cd = 255;

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/config.mk	Mon Dec 03 01:49:23 2007 +0100
@@ -0,0 +1,25 @@
+DEFS=-DPTHREADS
+
+CCFLAGS=-Wall -Wno-write-strings -fno-strict-aliasing -I$(ROOT)/include
+LDFLAGS=
+
+PY_CCFLAGS=$(shell python-config --includes)
+PY_LDFLAGS=$(shell python-config --libs)
+
+SDL_CCFLAGS=$(shell sdl-config --cflags)
+SDL_LDFLAGS=$(shell sdl-config --libs)
+
+ifeq ($(OS), Windows_NT)
+  LDFLAGS+=-lpthreadGC2
+  PY_CCFLAGS=-I"C:/Program Files/Python25/include"
+  PY_LDFLAGS=-L"C:\Program Files\Python25\libs" -lpython25
+  MODULENAME=raytracer.pyd
+else
+  CCFLAGS+=-pthread -fPIC
+  LDFLAGS+=-pthread
+  MODULENAME=raytracermodule.so
+endif
+
+# optimizations
+CCFLAGS+=-g -O0
+#CCFLAGS+=-O3 -pipe -fomit-frame-pointer -ffast-math -msse3

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/demos/Makefile	Mon Dec 03 01:49:23 2007 +0100
@@ -0,0 +1,16 @@
+ifndef $(ROOT)
+    ROOT=$(shell pwd)/..
+endif
+    
+
+all: ModulePath reqs models
+
+ModulePath:
+	echo "$(ROOT)/bin/python-module" > ModulePath
+
+reqs:
+	$(MAKE) -C .. libs-float python-module
+
+models: ;
+
+clean: ;

--- a/demos/spheres_ao.py	Fri Nov 30 00:44:51 2007 +0100
+++ b/demos/spheres_ao.py	Mon Dec 03 01:49:23 2007 +0100
@@ -1,7 +1,7 @@
 #!/usr/bin/python
 
 import sys
-sys.path.append("..")
+sys.path.append(open('ModulePath').read().strip())
 
 from raytracer import Raytracer, Material, Box, Sphere, Light
 import Image

--- a/demos/spheres_shadow.py	Fri Nov 30 00:44:51 2007 +0100
+++ b/demos/spheres_shadow.py	Mon Dec 03 01:49:23 2007 +0100
@@ -1,7 +1,7 @@
 #!/usr/bin/python
 
 import sys
-sys.path.append("..")
+sys.path.append(open('ModulePath').read().strip())
 
 from raytracer import Raytracer, Material, Box, Sphere, Light
 import Image

--- a/demos/triangles_monkey.py	Fri Nov 30 00:44:51 2007 +0100
+++ b/demos/triangles_monkey.py	Mon Dec 03 01:49:23 2007 +0100
@@ -1,7 +1,7 @@
 #!/usr/bin/python
 
 import sys
-sys.path.append("..")
+sys.path.append(open('ModulePath').read().strip())
 
 from raytracer import Raytracer, Light, Sphere, Triangle, Material
 import Image

--- a/demos/triangles_sphere.py	Fri Nov 30 00:44:51 2007 +0100
+++ b/demos/triangles_sphere.py	Mon Dec 03 01:49:23 2007 +0100
@@ -1,7 +1,7 @@
 #!/usr/bin/python
 
 import sys
-sys.path.append("..")
+sys.path.append(open('ModulePath').read().strip())
 
 from raytracer import Raytracer, Light, Sphere, Triangle, Material
 import Image

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/common.h	Mon Dec 03 01:49:23 2007 +0100
@@ -0,0 +1,28 @@
+#ifndef COMMON_H
+#define COMMON_H
+
+#include <stdio.h>
+#include <stdarg.h>
+#include <float.h>
+
+#ifdef PYRIT_DOUBLE
+# define Float double
+# define Eps DBL_EPSILON
+# define Inf DBL_MAX
+#else
+# define Float float
+# define Eps 1e-6
+# define Inf FLT_MAX
+#endif
+
+inline void infomsg(const char *format, ...)
+{
+#ifndef PYRIT_QUIET
+	va_list ap;
+	va_start(ap, format);
+	vprintf(format, ap);
+	va_end(ap);
+#endif
+}
+
+#endif

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/kdtree.h	Mon Dec 03 01:49:23 2007 +0100
@@ -0,0 +1,79 @@
+#ifndef KDTREE_H
+#define KDTREE_H
+
+#include <vector>
+#include <iostream>
+#include <fstream>
+
+#include "scene.h"
+
+using namespace std;
+
+class ShapeList: public vector<Shape*>
+{
+};
+
+class Container
+{
+protected:
+	BBox bbox;
+public:
+	ShapeList shapes;
+	Container(): bbox(), shapes() {};
+	virtual ~Container() {};
+	virtual void addShape(Shape* aShape);
+	//void addShapeNoExtend(Shape* aShape) { shapes.push_back(aShape); };
+	virtual Shape *nearest_intersection(const Shape *origin_shape, const Ray &ray,
+		Float &nearest_distance);
+	virtual void optimize() {};
+};
+
+class KdNode
+{
+	Float split;
+	short axis; /* 0,1,2 => x,y,z; 3 => leaf */
+public:
+	union {
+		KdNode *children;
+		ShapeList *shapes;
+	};
+
+	KdNode() : axis(3) { shapes = new ShapeList(); };
+	~KdNode();
+
+	void setAxis(short aAxis) { axis = aAxis; };
+	short getAxis() { return axis; };
+
+	void setSplit(Float aSplit) { split = aSplit; };
+	Float getSplit() { return split; };
+
+	void setLeaf() { axis = 3; };
+	bool isLeaf() { return axis == 3; };
+
+	KdNode *getLeftChild() { return children; };
+	KdNode *getRightChild() { return children+1; };
+
+	void addShape(Shape* aShape) { shapes->push_back(aShape); };
+
+	void subdivide(BBox bbox, int depth);
+};
+
+class KdTree: public Container
+{
+	KdNode *root;
+	bool built;
+	int max_depth;
+public:
+	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,
+		Float &nearest_distance);
+	void optimize() { build(); };
+	void build();
+	void save(ostream &str, KdNode *node = NULL);
+	void load(istream &str, KdNode *node = NULL);
+	void setMaxDepth(int md) { max_depth = md; };
+};
+
+#endif

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/matrix.h	Mon Dec 03 01:49:23 2007 +0100
@@ -0,0 +1,94 @@
+/*
+ * C++ RayTracer
+ * file: matrix.h
+ *
+ * Radek Brich, 2006
+ */
+
+/* not used at this time */
+
+#ifndef MATRIX_H
+#define MATRIX_H
+
+#include "vector.h"
+
+using namespace std;
+
+class Matrix
+{
+public:
+	Float data[4][4];
+
+	Matrix(): {};
+
+	// sum
+	friend Matrix operator+(const Matrix &a, const Matrix &b)
+	{
+		Matrix m = Matrix();
+		for (int i = 0; i < 4; i++)
+			for (int j = 0; j < 4; j++)
+				m.data[i][j] = a.data[i][j] + b.data[i][j];
+		return m;
+	}
+
+	// difference
+	friend Matrix operator-(const Matrix &a, const Matrix &b)
+	{
+		Matrix m = Matrix();
+		for (int i = 0; i < 4; i++)
+			for (int j = 0; j < 4; j++)
+				m.data[i][j] = a.data[i][j] - b.data[i][j];
+		return m;
+	}
+
+	// product
+	friend Matrix operator*(const Matrix &a, const Matrix &b)
+	{
+		Matrix m = Matrix();
+		for (int i = 0; i < 4; i++)
+			for (int j = 0; j < 4; j++)
+				m.data[i][j] =
+					a.data[i][0] * b.data[0][j] +
+					a.data[i][1] * b.data[1][j] +
+					a.data[i][2] * b.data[2][j] +
+					a.data[i][3] * b.data[3][j];
+		return m;
+	}
+
+	// negative
+	Matrix operator-()
+	{
+		Matrix m = Matrix();
+		for (int i = 0; i < 4; i++)
+			for (int j = 0; j < 4; j++)
+				m.data[i][j] = -data[i][j];
+		return m;
+	}
+
+	// product of matrix and scalar
+	Matrix operator*(const Float &f)
+	{
+		Matrix m = Matrix();
+		for (int i = 0; i < 4; i++)
+			for (int j = 0; j < 4; j++)
+				m.data[i][j] = data[i][j] * f;
+		return m;
+	}
+
+	friend Matrix operator*(const Float &f, Matrix &m) { return m * f; };
+
+	// product of matrix and vector
+	Vector3 operator*(const Vector3 &v)
+	{
+		Vector3 u = Vector3();
+		u.x = data[0][0] * v.x + data[0][1] * v.y + data[0][2] * v.z + data[0][3] * v.w;
+		u.y = data[1][0] * v.x + data[1][1] * v.y + data[1][2] * v.z + data[1][3] * v.w;
+		u.z = data[2][0] * v.x + data[2][1] * v.y + data[2][2] * v.z + data[2][3] * v.w;
+		u.w = data[3][0] * v.x + data[3][1] * v.y + data[3][2] * v.z + data[3][3] * v.w;
+		return u;
+	}
+
+	friend Matrix operator*(const Vector3 &v, Matrix &m) { return m * v; };
+};
+
+#endif

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/noise.h	Mon Dec 03 01:49:23 2007 +0100
@@ -0,0 +1,8 @@
+#ifndef NOISE_H
+#define NOISE_H
+
+#include "common.h"
+
+Float perlin(Float x, Float y, Float z);
+
+#endif

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/raytracer.h	Mon Dec 03 01:49:23 2007 +0100
@@ -0,0 +1,57 @@
+/*
+ * C++ RayTracer
+ * file: raytracer.h
+ *
+ * Radek Brich, 2006
+ */
+
+#ifndef RAYTRACER_H
+#define RAYTRACER_H
+
+#include <vector>
+
+#include "common.h"
+#include "kdtree.h"
+#include "scene.h"
+
+using namespace std;
+
+class Raytracer;
+struct RenderrowData {
+	Raytracer *rt;
+	int w;
+	Vector3 eye, dfix, dx, dy;
+	Float *iter;
+};
+
+class Raytracer
+{
+	Container *top;
+	Camera *camera;
+	vector<Light*> lights;
+	Colour bg_colour;
+	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), 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); };
+	void addlight(Light *light);
+	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; };
+};
+
+#endif

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/scene.h	Mon Dec 03 01:49:23 2007 +0100
@@ -0,0 +1,205 @@
+/*
+ * C++ RayTracer
+ * file: scene.h
+ *
+ * Radek Brich, 2006
+ */
+
+#ifndef SCENE_H
+#define SCENE_H
+
+#include <vector>
+
+#include "noise.h"
+
+#include "vector.h"
+
+using namespace std;
+
+class Ray
+{
+public:
+	Vector3 o, dir;
+	Ray(const Vector3 &ao, const Vector3 &adir):
+		o(ao), dir(adir) {};
+};
+
+class Quaternion
+{
+public:
+	Float a,b,c,d;
+	Quaternion(): a(0), b(0), c(0), d(0) {};
+	Quaternion(const Float aa, const Float ab, const Float ac, const Float ad):
+		a(aa), b(ab), c(ac), d(ad) {};
+	Quaternion(const Vector3& v): a(1), b(v.x), c(v.y), d(v.z) {};
+
+	Vector3 toVector() { return Vector3(b/a, c/a, d/a); };
+
+	Quaternion normalize()
+	{
+		Float f = 1.0f / sqrtf(a * a + b * b + c * c + d * d);
+		a *= f;
+		b *= f;
+		c *= f;
+		d *= f;
+		return *this;
+	};
+	friend Quaternion operator*(const Quaternion &q1, const Quaternion &q2)
+	{
+		return Quaternion(
+			q1.a*q2.a - q1.b*q2.b - q1.c*q2.c - q1.d*q2.d,
+			q1.a*q2.b + q1.b*q2.a + q1.c*q2.d - q1.d*q2.c,
+			q1.a*q2.c + q1.c*q2.a + q1.d*q2.b - q1.b*q2.d,
+			q1.a*q2.d + q1.d*q2.a + q1.b*q2.c - q1.c*q2.b);
+	};
+	friend Quaternion conjugate(const Quaternion &q)
+	{
+		return Quaternion(q.a, -q.b, -q.c, -q.d);
+	}
+};
+
+class Camera
+{
+public:
+	Vector3 eye, p, u, v;
+	Float f;
+
+	Camera(): eye(0,0,10), p(0,0,-1), u(-1,0,0), v(0,1,0), f(3.14/4.0) {};
+	Camera(const Vector3 &C, const Vector3 &ap, const Vector3 &au, const Vector3 &av):
+		eye(C), p(ap), u(au), v(av), f(3.14/4.0) {};
+	void setEye(const Vector3 &aeye) { eye = aeye; };
+	void setFocalLength(const Float af) { f = af; };
+	void rotate(const Quaternion &q);
+	void move(const Float fw, const Float left, const Float up);
+};
+
+/* axis-aligned bounding box */
+class BBox
+{
+public:
+	Vector3 L;
+	Vector3 H;
+	BBox(): L(), H() {};
+	BBox(const Vector3 aL, const Vector3 aH): L(aL), H(aH) {};
+	Float w() { return H.x-L.x; };
+	Float h() { return H.y-L.y; };
+	Float d() { return H.z-L.z; };
+	bool intersect(const Ray &ray, Float &a, Float &b);
+};
+
+class Light
+{
+public:
+	Vector3 pos;
+	Colour colour;
+	bool cast_shadows;
+
+	Light(const Vector3 &position, const Colour &acolour):
+		pos(position), colour(acolour), cast_shadows(true) {};
+	void castShadows(bool cast) { cast_shadows = cast; };
+};
+
+class Texture
+{
+public:
+	Colour colour;
+	Colour evaluate(Vector3 point)
+	{
+		Float sum = 0.0;
+		for (int i = 1; i < 5; i++)
+			sum += fabsf(perlin(point.x*i, point.y*i, point.z*i))/i;
+		Float value = sinf(point.x + sum)/2 + 0.5;
+		return Colour(value*colour.r, value*colour.g, value*colour.b);
+	};
+};
+
+class Material
+{
+public:
+	Float ambient, diffuse, specular, shininess; // Phong constants
+	Float reflection; // how much reflectife is the surface
+	Float refraction; // refraction index
+	Float transmitivity;
+	Texture texture;
+
+	Material(const Colour &acolour) {
+		texture.colour = acolour;
+		ambient = 0.1;
+		diffuse = 0.5;
+		specular = 0.1;
+		shininess = 0.5;
+		reflection = 0.5;
+	}
+};
+
+class Shape
+{
+public:
+	Material *material;
+	Shape() {};
+	virtual ~Shape() {};
+
+	// first intersection point
+	virtual bool intersect(const Ray &ray, Float &dist) = 0;
+
+	// all intersections (only for CSG)
+	virtual bool intersect_all(const Ray &ray, Float dist, vector<Float> &allts) = 0;
+
+	// normal at point P
+	virtual Vector3 normal(Vector3 &P) = 0;
+
+	virtual BBox get_bbox() = 0;
+};
+
+class Sphere: public Shape
+{
+	Float sqr_radius;
+	Float inv_radius;
+public:
+	Vector3 center;
+	Float radius;
+
+	Sphere(const Vector3 &acenter, const Float aradius, Material *amaterial):
+		sqr_radius(aradius*aradius), inv_radius(1.0f/aradius),
+		center(acenter), radius(aradius) { material = amaterial; }
+	bool intersect(const Ray &ray, Float &dist);
+	bool intersect_all(const Ray &ray, Float dist, vector<Float> &allts);
+	Vector3 normal(Vector3 &P) { return (P - center) * inv_radius; };
+	BBox get_bbox();
+};
+
+class Box: public Shape
+{
+	Vector3 L;
+	Vector3 H;
+public:
+	Box(const Vector3 &aL, const Vector3 &aH, Material *amaterial): L(aL), H(aH)
+	{
+		for (int i = 0; i < 3; i++)
+			if (L.cell[i] > H.cell[i])
+				swap(L.cell[i], H.cell[i]);
+		material = amaterial;
+	};
+	bool intersect(const Ray &ray, Float &dist);
+	bool intersect_all(const Ray &ray, Float dist, vector<Float> &allts) {return false;};
+	Vector3 normal(Vector3 &P);
+	BBox get_bbox() { return BBox(L, H); };
+};
+
+class Triangle: public Shape
+{
+	int k; // dominant axis
+	Float nu, nv, nd;
+	Float bnu, bnv;
+	Float cnu, cnv;
+public:
+	Vector3 A, B, C, N;
+
+	Triangle(const Vector3 &aA, const Vector3 &aB, const Vector3 &aC, Material *amaterial);
+	bool intersect(const Ray &ray, Float &dist);
+	bool intersect_all(const Ray &ray, Float dist, vector<Float> &allts) {return false;};
+	Vector3 normal(Vector3 &) { return N; };
+	BBox get_bbox();
+};
+
+#endif

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/vector.h	Mon Dec 03 01:49:23 2007 +0100
@@ -0,0 +1,136 @@
+/*
+ * C++ RayTracer
+ * file: vector.h
+ *
+ * Radek Brich, 2006
+ */
+
+#ifndef VECTOR_H
+#define VECTOR_H
+
+#include <math.h>
+#include <iostream>
+
+using namespace std;
+
+class Vector3
+{
+public:
+	// data
+	union {
+		struct {
+			Float x, y, z;
+		};
+		struct {
+			Float r, g, b;
+		};
+		Float cell[3];
+	};
+
+	// constructors
+	Vector3(): x(0.0f), y(0.0f), z(0.0f) {};
+	Vector3(Float ax, Float ay, Float az): x(ax), y(ay), z(az) {};
+
+	// 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()
+	{
+		Float f = 1.0f / mag();
+		x *= f;
+		y *= f;
+		z *= f;
+		return *this;
+	}
+
+	// get normalized copy
+	Vector3 unit()
+	{
+		Vector3 u(*this);
+		return u.normalize();;
+	}
+
+	// square magnitude, magnitude
+	Float mag2()	{ return x * x + y * y + z * z; }
+	Float mag()	{ return sqrtf(mag2()); }
+
+	// negative
+	Vector3 operator-()	{ return Vector3(-x, -y, -z); }
+
+	// accumulate
+	Vector3 operator+=(const Vector3 &v)
+	{
+		x += v.x;
+		y += v.y;
+		z += v.z;
+		return *this;
+	};
+
+	// sum
+	friend Vector3 operator+(const Vector3 &a, const Vector3 &b)
+	{
+		return Vector3(a.x + b.x, a.y + b.y, a.z + b.z);
+	};
+
+	// difference
+	friend Vector3 operator-(const Vector3 &a, const Vector3 &b)
+	{
+		return Vector3(a.x - b.x, a.y - b.y, a.z - b.z);
+	};
+
+	// dot product
+	friend Float dot(const Vector3 &a, const Vector3 &b)
+	{
+		return a.x * b.x + a.y * b.y + a.z * b.z;
+	};
+
+	// cross product
+	friend Vector3 cross(const Vector3 &a, const Vector3 &b)
+	{
+		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);
+	};
+
+	// product of vector and scalar
+	friend Vector3 operator*(const Vector3 &v, const Float &f)
+	{
+		return Vector3(f * v.x, f * v.y, f * v.z);
+	}
+
+	friend Vector3 operator*(const Float &f, const Vector3 &v)
+	{
+		return v * f;
+	};
+
+	// vector plus scalar
+	friend Vector3 operator+(const Vector3 &v, const Float &f)
+	{
+		return Vector3(v.x + f, v.y + f, v.z + f);
+	}
+
+	// vector minus scalar
+	friend Vector3 operator-(const Vector3 &v, const Float &f)
+	{
+		return Vector3(v.x - f, v.y - f, v.z - f);
+	}
+
+	// cell by cell product (only usable for colours)
+	friend Vector3 operator*(const Vector3 &a,  const Vector3 &b)
+	{
+		return Vector3(a.x * b.x, a.y * b.y, a.z * b.z);
+	};
+
+	// print
+	friend ostream & operator<<(ostream &st, const Vector3 &v)
+	{
+		return st << "(" << v.x << ", " << v.y  << ", " << v.z << ")";
+	}
+};
+
+typedef Vector3 Colour;
+
+#endif

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Makefile	Mon Dec 03 01:49:23 2007 +0100
@@ -0,0 +1,50 @@
+ifndef $(ROOT)
+	ROOT=$(shell pwd)/..
+endif
+
+include $(ROOT)/config.mk
+
+vpath %.cc $(ROOT)/src
+vpath %.h $(ROOT)/include
+LIBOBJS=raytracer.o scene.o noise.o kdtree.o
+CCFLAGS+=-I$(ROOT)/include
+
+### Targets ###
+all: libs-float libs-double python-module
+
+libs-float:
+	mkdir -p $(ROOT)/bin/$@
+	$(MAKE) -C $(ROOT)/bin/$@ -f $(ROOT)/src/Makefile libs ROOT="$(ROOT)"
+
+libs-double:
+	mkdir -p $(ROOT)/bin/$@
+	$(MAKE) -C $(ROOT)/bin/$@ -f $(ROOT)/src/Makefile libs ROOT="$(ROOT)" DEFS="$(DEFS) -DPYRIT_DOUBLE"
+
+libs: $(LIBOBJS)
+
+python-module: libs-float
+	mkdir -p $(ROOT)/bin/$@
+	$(MAKE) -C $(ROOT)/bin/$@ -f $(ROOT)/src/Makefile $(MODULENAME) ROOT="$(ROOT)"
+
+$(MODULENAME): raytracermodule.o
+	$(CXX) -shared -o $@ $< $(ROOT)/bin/libs-float/*.o $(LDFLAGS) $(PY_LDFLAGS)
+
+clean:
+	rm -rf $(ROOT)/bin/libs-*
+	rm -rf $(ROOT)/bin/python-module
+
+
+### Rules ###
+%.o: %.cc
+	$(CXX) -c -o $@ $(DEFS) $(CCFLAGS) $<
+
+
+### Dependencies ###
+matrix.o: matrix.cc matrix.h vector.h common.h
+noise.o: noise.cc noise.h common.h
+scene.o: scene.cc scene.h vector.h noise.h common.h
+kdtree.o: kdtree.cc kdtree.h scene.h common.h
+raytracer.o: raytracer.cc raytracer.h scene.h vector.h noise.h common.h
+
+raytracermodule.o: raytracermodule.cc raytracer.h scene.h vector.h common.h
+	$(CXX) -c -o $@ $(DEFS) $(CCFLAGS) $(PY_CCFLAGS) $<

--- a/src/common.h	Fri Nov 30 00:44:51 2007 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,21 +0,0 @@
-#ifndef COMMON_H
-#define COMMON_H
-
-#include <stdio.h>
-#include <stdarg.h>
-#include <float.h>
-
-#define Eps 1e-6
-#define Inf FLT_MAX
-
-inline void infomsg(const char *format, ...)
-{
-#ifndef PYRIT_QUITE
-	va_list ap;
-	va_start(ap, format);
-	vprintf(format, ap);
-	va_end(ap);
-#endif
-}
-
-#endif

--- a/src/kdtree.cc	Fri Nov 30 00:44:51 2007 +0100
+++ b/src/kdtree.cc	Mon Dec 03 01:49:23 2007 +0100
@@ -35,10 +35,10 @@
 class SplitPos
 {
 public:
-	float pos;
+	Float pos;
 	int lnum, rnum;
 	SplitPos(): pos(0.0), lnum(0), rnum(0) {};
-	SplitPos(float &aPos): pos(aPos), lnum(0), rnum(0) {};
+	SplitPos(Float &aPos): pos(aPos), lnum(0), rnum(0) {};
 	friend bool operator<(const SplitPos& a, const SplitPos& b)
 		{ return a.pos < b.pos; };
 	friend bool operator==(const SplitPos& a, const SplitPos& b)
@@ -54,9 +54,9 @@
 {
 public:
 	KdNode* node; /* pointer to far child */
-	float t; /* the entry/exit signed distance */
+	Float t; /* the entry/exit signed distance */
 	Vector3 pb; /* the coordinates of entry/exit point */
-	StackElem(KdNode *anode, const float &at, const Vector3 &apb):
+	StackElem(KdNode *anode, const Float &at, const Vector3 &apb):
 		node(anode), t(at), pb(apb) {};
 };
 
@@ -81,7 +81,7 @@
 };
 
 Shape *Container::nearest_intersection(const Shape *origin_shape, const Ray &ray,
-	float &nearest_distance)
+	Float &nearest_distance)
 {
 	Shape *nearest_shape = NULL;
 	ShapeList::iterator shape;
@@ -183,9 +183,9 @@
 	}
 
 	// choose best split pos
-	const float K = 1.4; // constant, K = cost of traversal / cost of ray-triangle intersection
-	float SAV = 2*(bbox.w()*bbox.h() + bbox.w()*bbox.d() + bbox.h()*bbox.d()); // surface area of node
-	float cost = SAV * (K + shapes->size()); // initial cost = non-split cost
+	const Float K = 1.4; // constant, K = cost of traversal / cost of ray-triangle intersection
+	Float SAV = 2*(bbox.w()*bbox.h() + bbox.w()*bbox.d() + bbox.h()*bbox.d()); // surface area of node
+	Float cost = SAV * (K + shapes->size()); // initial cost = non-split cost
 	SplitPos *splitpos = NULL;
 	bool leaf = true;
 	BBox lbb = bbox;
@@ -195,9 +195,9 @@
 		// calculate SAH cost of this split
 		lbb.H.cell[axis] = spl->pos;
 		rbb.L.cell[axis] = spl->pos;
-		float SAL = 2*(lbb.w()*lbb.h() + lbb.w()*lbb.d() + lbb.h()*lbb.d());
-		float SAR = 2*(rbb.w()*rbb.h() + rbb.w()*rbb.d() + rbb.h()*rbb.d());
-		float splitcost = K + SAL/SAV*(K+spl->lnum) + SAR/SAV*(K+spl->rnum);
+		Float SAL = 2*(lbb.w()*lbb.h() + lbb.w()*lbb.d() + lbb.h()*lbb.d());
+		Float SAR = 2*(rbb.w()*rbb.h() + rbb.w()*rbb.d() + rbb.h()*rbb.d());
+		Float splitcost = K + SAL/SAV*(K+spl->lnum) + SAR/SAV*(K+spl->rnum);
 
 		if (splitcost < cost)
 		{
@@ -237,8 +237,8 @@
 #endif
 
 	split = splitpos->pos;
-	float lnum = splitpos->lnum;
-	float rnum = splitpos->rnum;
+	Float lnum = splitpos->lnum;
+	Float rnum = splitpos->rnum;
 
 	// split this node
 	delete shapes;
@@ -315,10 +315,10 @@
 
 /* algorithm by Vlastimil Havran, Heuristic Ray Shooting Algorithms, appendix C */
 Shape *KdTree::nearest_intersection(const Shape *origin_shape, const Ray &ray,
-	float &nearest_distance)
+	Float &nearest_distance)
 {
-	float a, b; /* entry/exit signed distance */
-	float t;    /* signed distance to the splitting plane */
+	Float a, b; /* entry/exit signed distance */
+	Float t;    /* signed distance to the splitting plane */
 
 	/* if we have no tree, fall back to naive test */
 	if (!built)
@@ -352,7 +352,7 @@
 		while (!node->isLeaf())
 		{
 			/* retrieve position of splitting plane */
-			float splitVal = node->getSplit();
+			Float splitVal = node->getSplit();
 			short axis = node->getAxis();
 
 			if (enPt->pb[axis] <= splitVal)
@@ -401,7 +401,7 @@
 		/* "intersect ray with each object in the object list, discarding "
 		"those lying before stack[enPt].t or farther than stack[exPt].t" */
 		Shape *nearest_shape = NULL;
-		float dist = exPt->t;
+		Float dist = exPt->t;
 		ShapeList::iterator shape;
 		for (shape = node->shapes->begin(); shape != node->shapes->end(); shape++)
 			if (*shape != origin_shape && (*shape)->intersect(ray, dist)

--- a/src/kdtree.h	Fri Nov 30 00:44:51 2007 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,79 +0,0 @@
-#ifndef KDTREE_H
-#define KDTREE_H
-
-#include <vector>
-#include <iostream>
-#include <fstream>
-
-#include "scene.h"
-
-using namespace std;
-
-class ShapeList: public vector<Shape*>
-{
-};
-
-class Container
-{
-protected:
-	BBox bbox;
-public:
-	ShapeList shapes;
-	Container(): bbox(), shapes() {};
-	virtual ~Container() {};
-	virtual void addShape(Shape* aShape);
-	//void addShapeNoExtend(Shape* aShape) { shapes.push_back(aShape); };
-	virtual Shape *nearest_intersection(const Shape *origin_shape, const Ray &ray,
-		float &nearest_distance);
-	virtual void optimize() {};
-};
-
-class KdNode
-{
-	float split;
-	short axis; /* 0,1,2 => x,y,z; 3 => leaf */
-public:
-	union {
-		KdNode *children;
-		ShapeList *shapes;
-	};
-
-	KdNode() : axis(3) { shapes = new ShapeList(); };
-	~KdNode();
-
-	void setAxis(short aAxis) { axis = aAxis; };
-	short getAxis() { return axis; };
-
-	void setSplit(float aSplit) { split = aSplit; };
-	float getSplit() { return split; };
-
-	void setLeaf() { axis = 3; };
-	bool isLeaf() { return axis == 3; };
-
-	KdNode *getLeftChild() { return children; };
-	KdNode *getRightChild() { return children+1; };
-
-	void addShape(Shape* aShape) { shapes->push_back(aShape); };
-
-	void subdivide(BBox bbox, int depth);
-};
-
-class KdTree: public Container
-{
-	KdNode *root;
-	bool built;
-	int max_depth;
-public:
-	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,
-		float &nearest_distance);
-	void optimize() { build(); };
-	void build();
-	void save(ostream &str, KdNode *node = NULL);
-	void load(istream &str, KdNode *node = NULL);
-	void setMaxDepth(int md) { max_depth = md; };
-};
-
-#endif

--- a/src/matrix.h	Fri Nov 30 00:44:51 2007 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,94 +0,0 @@
-/*
- * C++ RayTracer
- * file: matrix.h
- *
- * Radek Brich, 2006
- */
-
-/* not used at this time */
-
-#ifndef MATRIX_H
-#define MATRIX_H
-
-#include "vector.h"
-
-using namespace std;
-
-class Matrix
-{
-public:
-	float data[4][4];
-
-	Matrix(): {};
-
-	// sum
-	friend Matrix operator+(const Matrix &a, const Matrix &b)
-	{
-		Matrix m = Matrix();
-		for (int i = 0; i < 4; i++)
-			for (int j = 0; j < 4; j++)
-				m.data[i][j] = a.data[i][j] + b.data[i][j];
-		return m;
-	}
-
-	// difference
-	friend Matrix operator-(const Matrix &a, const Matrix &b)
-	{
-		Matrix m = Matrix();
-		for (int i = 0; i < 4; i++)
-			for (int j = 0; j < 4; j++)
-				m.data[i][j] = a.data[i][j] - b.data[i][j];
-		return m;
-	}
-
-	// product
-	friend Matrix operator*(const Matrix &a, const Matrix &b)
-	{
-		Matrix m = Matrix();
-		for (int i = 0; i < 4; i++)
-			for (int j = 0; j < 4; j++)
-				m.data[i][j] =
-					a.data[i][0] * b.data[0][j] +
-					a.data[i][1] * b.data[1][j] +
-					a.data[i][2] * b.data[2][j] +
-					a.data[i][3] * b.data[3][j];
-		return m;
-	}
-
-	// negative
-	Matrix operator-()
-	{
-		Matrix m = Matrix();
-		for (int i = 0; i < 4; i++)
-			for (int j = 0; j < 4; j++)
-				m.data[i][j] = -data[i][j];
-		return m;
-	}
-
-	// product of matrix and scalar
-	Matrix operator*(const float &f)
-	{
-		Matrix m = Matrix();
-		for (int i = 0; i < 4; i++)
-			for (int j = 0; j < 4; j++)
-				m.data[i][j] = data[i][j] * f;
-		return m;
-	}
-
-	friend Matrix operator*(const float &f, Matrix &m) { return m * f; };
-
-	// product of matrix and vector
-	Vector3 operator*(const Vector3 &v)
-	{
-		Vector3 u = Vector3();
-		u.x = data[0][0] * v.x + data[0][1] * v.y + data[0][2] * v.z + data[0][3] * v.w;
-		u.y = data[1][0] * v.x + data[1][1] * v.y + data[1][2] * v.z + data[1][3] * v.w;
-		u.z = data[2][0] * v.x + data[2][1] * v.y + data[2][2] * v.z + data[2][3] * v.w;
-		u.w = data[3][0] * v.x + data[3][1] * v.y + data[3][2] * v.z + data[3][3] * v.w;
-		return u;
-	}
-
-	friend Matrix operator*(const Vector3 &v, Matrix &m) { return m * v; };
-};
-
-#endif

--- a/src/noise.cc	Fri Nov 30 00:44:51 2007 +0100
+++ b/src/noise.cc	Mon Dec 03 01:49:23 2007 +0100
@@ -8,20 +8,20 @@
 
 #include "noise.h"
 
-double fade(double t)
+Float fade(Float t)
 {
       return t * t * t * (t * (t * 6 - 15) + 10);
 }
 
-double lerp(double t, double a, double b)
+Float lerp(Float t, Float a, Float b)
 {
       return a + t * (b - a);
 }
 
-double grad(unsigned char hash, double x, double y, double z)
+Float grad(unsigned char hash, Float x, Float y, Float z)
 {
       unsigned char h = hash & 15;            // CONVERT LO 4 BITS OF HASH CODE
-      double u = h<8 ? x : y,                 // INTO 12 GRADIENT DIRECTIONS.
+      Float u = h<8 ? x : y,                 // INTO 12 GRADIENT DIRECTIONS.
              v = h<4 ? y : h==12||h==14 ? x : z;
       return ((h&1) == 0 ? u : -u) + ((h&2) == 0 ? v : -v);
 }
@@ -56,7 +56,7 @@
    };
 
 
-double perlin(double x, double y, double z)
+Float perlin(Float x, Float y, Float z)
 {
       int X = (int)floor(x) & 255,                  // FIND UNIT CUBE THAT
           Y = (int)floor(y) & 255,                  // CONTAINS POINT.
@@ -64,7 +64,7 @@
       x -= floor(x);                                // FIND RELATIVE X,Y,Z
       y -= floor(y);                                // OF POINT IN CUBE.
       z -= floor(z);
-      double u = fade(x),                                // COMPUTE FADE CURVES
+      Float u = fade(x),                                // COMPUTE FADE CURVES
              v = fade(y),                                // FOR EACH OF X,Y,Z.
              w = fade(z);
       int A = p[X  ]+Y, AA = p[A]+Z, AB = p[A+1]+Z,      // HASH COORDINATES OF

--- a/src/noise.h	Fri Nov 30 00:44:51 2007 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,6 +0,0 @@
-#ifndef NOISE_H
-#define NOISE_H
-
-double perlin(double x, double y, double z);
-
-#endif

--- a/src/raytracer.cc	Fri Nov 30 00:44:51 2007 +0100
+++ b/src/raytracer.cc	Mon Dec 03 01:49:23 2007 +0100
@@ -15,9 +15,9 @@
 
 // Hammersley spherical point distribution
 // http://www.cse.cuhk.edu.hk/~ttwong/papers/udpoint/udpoints.html
-Vector3 Raytracer::SphereDistribute(int i, int n, float extent, Vector3 &normal)
+Vector3 Raytracer::SphereDistribute(int i, int n, Float extent, Vector3 &normal)
 {
-	float p, t, st, phi, phirad;
+	Float p, t, st, phi, phirad;
 	int kk;
 
 	t = 0;
@@ -31,7 +31,7 @@
 
 	st = sqrt(1.0 - t*t);
 
-	float x, y, z, xx, yy, zz, q;
+	Float x, y, z, xx, yy, zz, q;
 	x = st * cos(phirad);
 	y = st * sin(phirad);
 	z = t;
@@ -75,8 +75,8 @@
 	Colour I = Colour();
 	Vector3 L = light.pos - P;
 	L.normalize();
-	float L_dot_N = dot(L, N);
-	float R_dot_V = dot(R, V);
+	Float L_dot_N = dot(L, N);
+	Float R_dot_V = dot(R, V);
 
 	Colour col = mat.texture.colour; //mat.texture.evaluate(P);
 
@@ -91,7 +91,7 @@
 
 Colour Raytracer::raytrace(Ray &ray, int depth, Shape *origin_shape)
 {
-	float nearest_distance = FLT_MAX; //Infinity
+	Float nearest_distance = Inf;
 	Shape *nearest_shape = top->nearest_intersection(origin_shape, ray, nearest_distance);
 
 	if (nearest_shape == NULL) {
@@ -106,12 +106,12 @@
 		for (light = lights.begin(); light != lights.end(); light++) {
 			Vector3 jo, L = (*light)->pos - P; // direction vector to light
 			L.normalize();
-			float L_dot_N = dot(L, normal);
+			Float L_dot_N = dot(L, normal);
 			if (L_dot_N > 0) {
 				// test if this light is occluded (sharp shadows)
 				if ((*light)->cast_shadows) {
 					Ray shadow_ray = Ray(P, L);
-					float dist = FLT_MAX;
+					Float dist = FLT_MAX;
 					if (top->nearest_intersection(nearest_shape, shadow_ray, dist))
 						continue;
 				}
@@ -137,18 +137,18 @@
 		// ambient occlusion
 		if (ao_samples)
 		{
-			float miss = 0;
+			Float miss = 0;
 			for (int i = 0; i < ao_samples; i++) {
 				Vector3 dir = SphereDistribute(i, ao_samples, ao_angle, normal);
 				Ray ao_ray = Ray(P, dir);
-				float dist = ao_distance;
+				Float dist = ao_distance;
 				Shape *shape_in_way = top->nearest_intersection(nearest_shape, ao_ray, dist);
 				if (shape_in_way == NULL)
 					miss += 1.0;
 				else
 					miss += dist / ao_distance;
 			}
-			float ao_intensity = miss / ao_samples;
+			Float ao_intensity = miss / ao_samples;
 			acc = acc * ao_intensity;
 		}
 
@@ -160,7 +160,7 @@
 {
 	RenderrowData *d = (RenderrowData*) data;
 	int subsample = d->rt->getSubsample();
-	float subsample2 = 1.0/(subsample*subsample);
+	Float subsample2 = 1.0/(subsample*subsample);
 	int ww = d->w*3;
 	Vector3 dir = d->dfix;
 	for (int x = 0; x < d->w; x += subsample) {
@@ -171,8 +171,8 @@
 
 		for (int i = 0; i < 5; i++)
 		{
-			float osax[] = {0.0, -0.4, +0.4, +0.4, -0.4};
-			float osay[] = {0.0, -0.4, -0.4, +0.4, +0.4};
+			Float osax[] = {0.0, -0.4, +0.4, +0.4, -0.4};
+			Float osay[] = {0.0, -0.4, -0.4, +0.4, +0.4};
 			Vector3 tmpdir = dir + osax[i]*d->dx + osay[i]*d->dy;
 			tmpdir.normalize();
 			Ray ray(d->eye, tmpdir);
@@ -204,14 +204,14 @@
 			ray.dir = tmpdir;
 			Colour c3 = d->rt->raytrace(ray, 0, NULL);
 			// are the colors similar?
-			float m = (c-c2).mag2();
+			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;
+				Float *i = d->iter;
 				for (int x = 0; x < subsample; x++)
 				{
 					for (int y = 0; y < subsample; y++)
@@ -302,14 +302,14 @@
 	return (void *)d;
 }
 
-void Raytracer::render(int w, int h, float *buffer)
+void Raytracer::render(int w, int h, Float *buffer)
 {
 	if (!camera || !top || !buffer)
 		return;
 
 	RenderrowData *d;
 
-	float S = 0.5/w;
+	Float S = 0.5/w;
 	Vector3 dfix = camera->u*(-w/2.0*S/camera->f)
 		+ camera->v*(h/2.0*S/camera->f) + camera->p;
 	Vector3 dx = camera->u * (S/camera->f);
@@ -372,7 +372,7 @@
 	lights.push_back(light);
 }
 
-void Raytracer::ambientocclusion(int samples, float distance, float angle)
+void Raytracer::ambientocclusion(int samples, Float distance, Float angle)
 {
 	ao_samples = samples;
 	ao_distance = distance;

--- a/src/raytracer.h	Fri Nov 30 00:44:51 2007 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,57 +0,0 @@
-/*
- * C++ RayTracer
- * file: raytracer.h
- *
- * Radek Brich, 2006
- */
-
-#ifndef RAYTRACER_H
-#define RAYTRACER_H
-
-#include <vector>
-
-#include "common.h"
-#include "kdtree.h"
-#include "scene.h"
-
-using namespace std;
-
-class Raytracer;
-struct RenderrowData {
-	Raytracer *rt;
-	int w;
-	Vector3 eye, dfix, dx, dy;
-	float *iter;
-};
-
-class Raytracer
-{
-	Container *top;
-	Camera *camera;
-	vector<Light*> lights;
-	Colour bg_colour;
-	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), 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); };
-	void addlight(Light *light);
-	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; };
-};
-
-#endif

--- a/src/raytracermodule.cc	Fri Nov 30 00:44:51 2007 +0100
+++ b/src/raytracermodule.cc	Mon Dec 03 01:49:23 2007 +0100
@@ -53,8 +53,8 @@
 	LightObject *v;
 	static char *kwdlist[] = {"position", "colour", NULL};
 	PyObject *TPos, *TCol = NULL;
-	float px, py, pz;
-	float cr = 1.0, cg = 1.0, cb = 1.0;
+	Float px, py, pz;
+	Float cr = 1.0, cg = 1.0, cb = 1.0;
 
 	if (!PyArg_ParseTupleAndKeywords(args, kwd, "O!|O!", kwdlist,
 		&PyTuple_Type, &TPos, &PyTuple_Type, &TCol))
@@ -133,7 +133,7 @@
 	MaterialObject *v;
 	static char *kwdlist[] = {"colour", NULL};
 	PyObject *TCol = NULL;
-	float cr=1.0, cg=1.0, cb=1.0;
+	Float cr=1.0, cg=1.0, cb=1.0;
 
 	if (!PyArg_ParseTupleAndKeywords(args, kwd, "|O!", kwdlist,
 		&PyTuple_Type, &TCol))
@@ -198,7 +198,7 @@
 	MaterialObject *material;
 	static char *kwdlist[] = {"centre", "radius", "material", NULL};
 	PyObject *TCentre = NULL;
-	float cx, cy, cz, radius;
+	Float cx, cy, cz, radius;
 
 	if (!PyArg_ParseTupleAndKeywords(args, kwd, "O!fO!", kwdlist,
 		&PyTuple_Type, &TCentre, &radius, &MaterialType, &material))
@@ -265,7 +265,7 @@
 	static char *kwdlist[] = {"L", "H", "material", NULL};
 	PyObject *TL = NULL;
 	PyObject *TH = NULL;
-	float lx, ly, lz, hx, hy, hz;
+	Float lx, ly, lz, hx, hy, hz;
 
 	if (!PyArg_ParseTupleAndKeywords(args, kwd, "O!O!O!", kwdlist,
 		&PyTuple_Type, &TL, &PyTuple_Type, &TH, &MaterialType, &material))
@@ -334,7 +334,7 @@
 	MaterialObject *material;
 	static char *kwdlist[] = {"A", "B", "C", "material", NULL};
 	PyObject *A = NULL, *B = NULL, *C = NULL;
-	float ax, ay, az, bx, by, bz, cx, cy, cz;
+	Float ax, ay, az, bx, by, bz, cx, cy, cz;
 
 	if (!PyArg_ParseTupleAndKeywords(args, kwd, "O!O!O!O!", kwdlist,
 		&PyTuple_Type, &A, &PyTuple_Type, &B, &PyTuple_Type, &C,
@@ -447,7 +447,7 @@
 {
 	int w = 0, h = 0;
 	char *chardata;
-	float *data;
+	Float *data;
 	PyObject *o;
 
 	if (!PyArg_ParseTuple(args, "(ii)", &w, &h))
@@ -455,7 +455,7 @@
 
 	printf("[pyrit] Raytracing...\n");
 	((RaytracerObject *)self)->raytracer->getTop()->optimize();
-	data = (float *) malloc(w*h*3*sizeof(float));
+	data = (Float *) malloc(w*h*3*sizeof(Float));
 	((RaytracerObject *)self)->raytracer->render(w, h, data);
 	if (!data) {
 		Py_INCREF(Py_None);
@@ -465,7 +465,7 @@
 	// convert data to char
 	printf("[pyrit] Converting image data (float to char)...\n");
 	chardata = (char *) malloc(w*h*3);
-	float *d = data;
+	Float *d = data;
 	for (char *c = chardata; c != chardata + w*h*3; c++, d++) {
 		if (*d > 1.0)
 			*c = 255;
@@ -511,7 +511,7 @@
 static PyObject* Raytracer_ambientocclusion(PyObject* self, PyObject* args, PyObject *kwd)
 {
 	int samples = 0;
-	float distance = 0.0, angle = 0.0;
+	Float distance = 0.0, angle = 0.0;
 	static char *kwdlist[] = {"samples", "distance", "angle", NULL};
 
 	if (!PyArg_ParseTupleAndKeywords(args, kwd, "iff", kwdlist,

--- a/src/scene.cc	Fri Nov 30 00:44:51 2007 +0100
+++ b/src/scene.cc	Mon Dec 03 01:49:23 2007 +0100
@@ -24,16 +24,16 @@
 	*/
 
 	// optimized
-	float t2 =   q.a*q.b;
-	float t3 =   q.a*q.c;
-	float t4 =   q.a*q.d;
-	float t5 =  -q.b*q.b;
-	float t6 =   q.b*q.c;
-	float t7 =   q.b*q.d;
-	float t8 =  -q.c*q.c;
-	float t9 =   q.c*q.d;
-	float t10 = -q.d*q.d;
-	float x,y,z;
+	Float t2 =   q.a*q.b;
+	Float t3 =   q.a*q.c;
+	Float t4 =   q.a*q.d;
+	Float t5 =  -q.b*q.b;
+	Float t6 =   q.b*q.c;
+	Float t7 =   q.b*q.d;
+	Float t8 =  -q.c*q.c;
+	Float t9 =   q.c*q.d;
+	Float t10 = -q.d*q.d;
+	Float x,y,z;
 	x = 2*( (t8 + t10)*p.x + (t6 -  t4)*p.y + (t3 + t7)*p.z ) + p.x;
 	y = 2*( (t4 +  t6)*p.x + (t5 + t10)*p.y + (t9 - t2)*p.z ) + p.y;
 	z = 2*( (t7 -  t3)*p.x + (t2 +  t9)*p.y + (t5 + t8)*p.z ) + p.z;
@@ -51,17 +51,17 @@
 	v.normalize();
 }
 
-void Camera::move(const float fw, const float left, const float up)
+void Camera::move(const Float fw, const Float left, const Float up)
 {
 	eye = eye + fw*p + left*u + up*v;
 }
 
 /* http://www.siggraph.org/education/materials/HyperGraph/raytrace/rtinter3.htm */
-bool BBox::intersect(const Ray &ray, float &a, float &b)
+bool BBox::intersect(const Ray &ray, Float &a, Float &b)
 {
-	float tnear = -FLT_MAX;
-	float tfar = FLT_MAX;
-	float t1, t2;
+	Float tnear = -FLT_MAX;
+	Float tfar = FLT_MAX;
+	Float t1, t2;
 
 	for (int i = 0; i < 3; i++)
 	{
@@ -94,11 +94,11 @@
 	return true;
 }
 
-bool Sphere::intersect(const Ray &ray, float &dist)
+bool Sphere::intersect(const Ray &ray, Float &dist)
 {
 	Vector3 V = ray.o - center;
-	register float d = -dot(V, ray.dir);
-	register float Det = d * d - (dot(V,V) - sqr_radius);
+	register Float d = -dot(V, ray.dir);
+	register Float Det = d * d - (dot(V,V) - sqr_radius);
 	if (Det > 0) {
 		d -= sqrtf(Det);
 		if (d > 0 && d < dist)
@@ -110,18 +110,18 @@
 	return false;
 }
 
-bool Sphere::intersect_all(const Ray &ray, float dist, vector<float> &allts)
+bool Sphere::intersect_all(const Ray &ray, Float dist, vector<Float> &allts)
 {
-	//allts = new vector<float>();
+	//allts = new vector<Float>();
 
 	Vector3 V = ((Ray)ray).o - center;
-	float Vd = - dot(V, ray.dir);
-	float Det = Vd * Vd - (dot(V,V) - sqr_radius);
+	Float Vd = - dot(V, ray.dir);
+	Float Det = Vd * Vd - (dot(V,V) - sqr_radius);
 
 	if (Det > 0) {
 		Det = sqrtf(Det);
-		float t1 = Vd - Det;
-		float t2 = Vd + Det;
+		Float t1 = Vd - Det;
+		Float t2 = Vd + Det;
 		if (t1 < 0)
 		{
 			if (t2 > 0)
@@ -152,9 +152,9 @@
 	return bbox;
 }
 
-bool Box::intersect(const Ray &ray, float &dist)
+bool Box::intersect(const Ray &ray, Float &dist)
 {
-	float a,b;
+	Float a,b;
 	bool res = get_bbox().intersect(ray, a, b);
 	if (res && a < dist)
 	{
@@ -216,13 +216,13 @@
 	int u = (k + 1) % 3;
 	int v = (k + 2) % 3;
 
-	float krec = 1.0f / N[k];
+	Float krec = 1.0f / N[k];
 	nu = N[u] * krec;
 	nv = N[v] * krec;
 	nd = dot(N, A) * krec;
 
 	// first line equation
-	float reci = 1.0f / (b[u] * c[v] - b[v] * c[u]);
+	Float reci = 1.0f / (b[u] * c[v] - b[v] * c[u]);
 	bnu = b[u] * reci;
 	bnv = -b[v] * reci;
 
@@ -235,7 +235,7 @@
 }
 
 // see comment for previous method
-bool Triangle::intersect(const Ray &ray, float &dist)
+bool Triangle::intersect(const Ray &ray, Float &dist)
 {
 	Vector3 O = ray.o;
 	Vector3 D = ray.dir;
@@ -243,20 +243,20 @@
 	const int modulo3[5] = {0,1,2,0,1};
 	const int ku = modulo3[k+1];
 	const int kv = modulo3[k+2];
-	const float lnd = 1.0f / (D[k] + nu * D[ku] + nv * D[kv]);
-	const float t = (nd - O[k] - nu * O[ku] - nv * O[kv]) * lnd;
+	const Float lnd = 1.0f / (D[k] + nu * D[ku] + nv * D[kv]);
+	const Float t = (nd - O[k] - nu * O[ku] - nv * O[kv]) * lnd;
 
 	if (!(t < dist && t > 0))
 		return false;
 
-	float hu = O[ku] + t * D[ku] - A[ku];
-	float hv = O[kv] + t * D[kv] - A[kv];
-	float beta = hv * bnu + hu * bnv;
+	Float hu = O[ku] + t * D[ku] - A[ku];
+	Float hv = O[kv] + t * D[kv] - A[kv];
+	Float beta = hv * bnu + hu * bnv;
 
 	if (beta < 0)
 		return false;
 
-	float gamma = hu * cnu + hv * cnv;
+	Float gamma = hu * cnu + hv * cnv;
 	if (gamma < 0)
 		return false;
 

--- a/src/scene.h	Fri Nov 30 00:44:51 2007 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,205 +0,0 @@
-/*
- * C++ RayTracer
- * file: scene.h
- *
- * Radek Brich, 2006
- */
-
-#ifndef SCENE_H
-#define SCENE_H
-
-#include <vector>
-
-#include "noise.h"
-
-#include "vector.h"
-
-using namespace std;
-
-class Ray
-{
-public:
-	Vector3 o, dir;
-	Ray(const Vector3 &ao, const Vector3 &adir):
-		o(ao), dir(adir) {};
-};
-
-class Quaternion
-{
-public:
-	float a,b,c,d;
-	Quaternion(): a(0), b(0), c(0), d(0) {};
-	Quaternion(const float aa, const float ab, const float ac, const float ad):
-		a(aa), b(ab), c(ac), d(ad) {};
-	Quaternion(const Vector3& v): a(1), b(v.x), c(v.y), d(v.z) {};
-
-	Vector3 toVector() { return Vector3(b/a, c/a, d/a); };
-
-	Quaternion normalize()
-	{
-		float f = 1.0f / sqrtf(a * a + b * b + c * c + d * d);
-		a *= f;
-		b *= f;
-		c *= f;
-		d *= f;
-		return *this;
-	};
-	friend Quaternion operator*(const Quaternion &q1, const Quaternion &q2)
-	{
-		return Quaternion(
-			q1.a*q2.a - q1.b*q2.b - q1.c*q2.c - q1.d*q2.d,
-			q1.a*q2.b + q1.b*q2.a + q1.c*q2.d - q1.d*q2.c,
-			q1.a*q2.c + q1.c*q2.a + q1.d*q2.b - q1.b*q2.d,
-			q1.a*q2.d + q1.d*q2.a + q1.b*q2.c - q1.c*q2.b);
-	};
-	friend Quaternion conjugate(const Quaternion &q)
-	{
-		return Quaternion(q.a, -q.b, -q.c, -q.d);
-	}
-};
-
-class Camera
-{
-public:
-	Vector3 eye, p, u, v;
-	float f;
-
-	Camera(): eye(0,0,10), p(0,0,-1), u(-1,0,0), v(0,1,0), f(3.14/4.0) {};
-	Camera(const Vector3 &C, const Vector3 &ap, const Vector3 &au, const Vector3 &av):
-		eye(C), p(ap), u(au), v(av), f(3.14/4.0) {};
-	void setEye(const Vector3 &aeye) { eye = aeye; };
-	void setFocalLength(const float af) { f = af; };
-	void rotate(const Quaternion &q);
-	void move(const float fw, const float left, const float up);
-};
-
-/* axis-aligned bounding box */
-class BBox
-{
-public:
-	Vector3 L;
-	Vector3 H;
-	BBox(): L(), H() {};
-	BBox(const Vector3 aL, const Vector3 aH): L(aL), H(aH) {};
-	float w() { return H.x-L.x; };
-	float h() { return H.y-L.y; };
-	float d() { return H.z-L.z; };
-	bool intersect(const Ray &ray, float &a, float &b);
-};
-
-class Light
-{
-public:
-	Vector3 pos;
-	Colour colour;
-	bool cast_shadows;
-
-	Light(const Vector3 &position, const Colour &acolour):
-		pos(position), colour(acolour), cast_shadows(true) {};
-	void castShadows(bool cast) { cast_shadows = cast; };
-};
-
-class Texture
-{
-public:
-	Colour colour;
-	Colour evaluate(Vector3 point)
-	{
-		float sum = 0.0;
-		for (int i = 1; i < 5; i++)
-			sum += fabsf(perlin(point.x*i, point.y*i, point.z*i))/i;
-		float value = sinf(point.x + sum)/2 + 0.5;
-		return Colour(value*colour.r, value*colour.g, value*colour.b);
-	};
-};
-
-class Material
-{
-public:
-	float ambient, diffuse, specular, shininess; // Phong constants
-	float reflection; // how much reflectife is the surface
-	float refraction; // refraction index
-	float transmitivity;
-	Texture texture;
-
-	Material(const Colour &acolour) {
-		texture.colour = acolour;
-		ambient = 0.1;
-		diffuse = 0.5;
-		specular = 0.1;
-		shininess = 0.5;
-		reflection = 0.5;
-	}
-};
-
-class Shape
-{
-public:
-	Material *material;
-	Shape() {};
-	virtual ~Shape() {};
-
-	// first intersection point
-	virtual bool intersect(const Ray &ray, float &dist) = 0;
-
-	// all intersections (only for CSG)
-	virtual bool intersect_all(const Ray &ray, float dist, vector<float> &allts) = 0;
-
-	// normal at point P
-	virtual Vector3 normal(Vector3 &P) = 0;
-
-	virtual BBox get_bbox() = 0;
-};
-
-class Sphere: public Shape
-{
-	float sqr_radius;
-	float inv_radius;
-public:
-	Vector3 center;
-	float radius;
-
-	Sphere(const Vector3 &acenter, const float aradius, Material *amaterial):
-		sqr_radius(aradius*aradius), inv_radius(1.0f/aradius),
-		center(acenter), radius(aradius) { material = amaterial; }
-	bool intersect(const Ray &ray, float &dist);
-	bool intersect_all(const Ray &ray, float dist, vector<float> &allts);
-	Vector3 normal(Vector3 &P) { return (P - center) * inv_radius; };
-	BBox get_bbox();
-};
-
-class Box: public Shape
-{
-	Vector3 L;
-	Vector3 H;
-public:
-	Box(const Vector3 &aL, const Vector3 &aH, Material *amaterial): L(aL), H(aH)
-	{
-		for (int i = 0; i < 3; i++)
-			if (L.cell[i] > H.cell[i])
-				swap(L.cell[i], H.cell[i]);
-		material = amaterial;
-	};
-	bool intersect(const Ray &ray, float &dist);
-	bool intersect_all(const Ray &ray, float dist, vector<float> &allts) {return false;};
-	Vector3 normal(Vector3 &P);
-	BBox get_bbox() { return BBox(L, H); };
-};
-
-class Triangle: public Shape
-{
-	int k; // dominant axis
-	float nu, nv, nd;
-	float bnu, bnv;
-	float cnu, cnv;
-public:
-	Vector3 A, B, C, N;
-
-	Triangle(const Vector3 &aA, const Vector3 &aB, const Vector3 &aC, Material *amaterial);
-	bool intersect(const Ray &ray, float &dist);
-	bool intersect_all(const Ray &ray, float dist, vector<float> &allts) {return false;};
-	Vector3 normal(Vector3 &) { return N; };
-	BBox get_bbox();
-};
-
-#endif

--- a/src/vector.h	Fri Nov 30 00:44:51 2007 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,136 +0,0 @@
-/*
- * C++ RayTracer
- * file: vector.h
- *
- * Radek Brich, 2006
- */
-
-#ifndef VECTOR_H
-#define VECTOR_H
-
-#include <math.h>
-#include <iostream>
-
-using namespace std;
-
-class Vector3
-{
-public:
-	// data
-	union {
-		struct {
-			float x, y, z;
-		};
-		struct {
-			float r, g, b;
-		};
-		float cell[3];
-	};
-
-	// constructors
-	Vector3(): x(0.0f), y(0.0f), z(0.0f) {};
-	Vector3(float ax, float ay, float az): x(ax), y(ay), z(az) {};
-
-	// 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()
-	{
-		float f = 1.0f / mag();
-		x *= f;
-		y *= f;
-		z *= f;
-		return *this;
-	}
-
-	// get normalized copy
-	Vector3 unit()
-	{
-		Vector3 u(*this);
-		return u.normalize();;
-	}
-
-	// square magnitude, magnitude
-	float mag2()	{ return x * x + y * y + z * z; }
-	float mag()	{ return sqrtf(mag2()); }
-
-	// negative
-	Vector3 operator-()	{ return Vector3(-x, -y, -z); }
-
-	// accumulate
-	Vector3 operator+=(const Vector3 &v)
-	{
-		x += v.x;
-		y += v.y;
-		z += v.z;
-		return *this;
-	};
-
-	// sum
-	friend Vector3 operator+(const Vector3 &a, const Vector3 &b)
-	{
-		return Vector3(a.x + b.x, a.y + b.y, a.z + b.z);
-	};
-
-	// difference
-	friend Vector3 operator-(const Vector3 &a, const Vector3 &b)
-	{
-		return Vector3(a.x - b.x, a.y - b.y, a.z - b.z);
-	};
-
-	// dot product
-	friend float dot(const Vector3 &a, const Vector3 &b)
-	{
-		return a.x * b.x + a.y * b.y + a.z * b.z;
-	};
-
-	// cross product
-	friend Vector3 cross(const Vector3 &a, const Vector3 &b)
-	{
-		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);
-	};
-
-	// product of vector and scalar
-	friend Vector3 operator*(const Vector3 &v, const float &f)
-	{
-		return Vector3(f * v.x, f * v.y, f * v.z);
-	}
-
-	friend Vector3 operator*(const float &f, const Vector3 &v)
-	{
-		return v * f;
-	};
-
-	// vector plus scalar
-	friend Vector3 operator+(const Vector3 &v, const float &f)
-	{
-		return Vector3(v.x + f, v.y + f, v.z + f);
-	}
-
-	// vector minus scalar
-	friend Vector3 operator-(const Vector3 &v, const float &f)
-	{
-		return Vector3(v.x - f, v.y - f, v.z - f);
-	}
-
-	// cell by cell product (only usable for colours)
-	friend Vector3 operator*(const Vector3 &a,  const Vector3 &b)
-	{
-		return Vector3(a.x * b.x, a.y * b.y, a.z * b.z);
-	};
-
-	// print
-	friend ostream & operator<<(ostream &st, const Vector3 &v)
-	{
-		return st << "(" << v.x << ", " << v.y  << ", " << v.z << ")";
-	}
-};
-
-typedef Vector3 Colour;
-
-#endif