--- a/TODO	Tue Apr 22 13:33:12 2008 +0200
+++ b/TODO	Wed Apr 23 10:38:33 2008 +0200
@@ -7,10 +7,6 @@
  * textures (3D procedural, pixmaps)
  * generalization: Camera "shader" (ray generator), surface shader and maybe light & background shaders
  * namespace
- * kd-tree:
-   - optimize structures
-   - optimize construction: use box-shape intersection instead of bounding boxes of shapes
-   - save/load
  * Python binding for all classes
  * stochastic oversampling
  * absorbtion of refracted rays in dense materials (can be computed using shape distance and some 'absorbance' constant)

--- a/include/common.h	Tue Apr 22 13:33:12 2008 +0200
+++ b/include/common.h	Wed Apr 23 10:38:33 2008 +0200
@@ -31,6 +31,9 @@
 #include <stdarg.h>
 #include <float.h>
 #include <pthread.h>
+#include <string>
+
+using namespace std;
 
 #ifdef PYRIT_DOUBLE
 # define Float double
@@ -107,4 +110,19 @@
 	}
 }
 
+inline void trim(string& s)
+{
+	const char *ws = " \n";
+	string::size_type pos = s.find_last_not_of(ws);
+	if (pos != string::npos)
+	{
+		s.erase(pos + 1);
+		pos = s.find_first_not_of(ws);
+		if (pos != string::npos)
+			s.erase(0, pos);
+	}
+	else
+		s.erase(s.begin(), s.end());
+}
+
 #endif

--- a/include/container.h	Tue Apr 22 13:33:12 2008 +0200
+++ b/include/container.h	Wed Apr 23 10:38:33 2008 +0200
@@ -29,7 +29,7 @@
 
 #include <vector>
 
-#include "scene.h"
+#include "shapes.h"
 
 using namespace std;
 
@@ -49,6 +49,10 @@
 	virtual Shape *nearest_intersection(const Shape *origin_shape, const Ray &ray,
 		Float &nearest_distance);
 	virtual void optimize() {};
+
+	ShapeList & getShapes() { return shapes; };
+
+	virtual ostream & dump(ostream &st);
 };
 
 #endif

--- a/include/kdtree.h	Tue Apr 22 13:33:12 2008 +0200
+++ b/include/kdtree.h	Wed Apr 23 10:38:33 2008 +0200
@@ -53,19 +53,19 @@
 	~KdNode();
 
 	void setLeaf() { flags |= 3; };
-	bool isLeaf() { return (flags & 3) == 3; };
+	const bool isLeaf() const { return (flags & 3) == 3; };
 
 	void setAxis(int aAxis) { flags &= ~3; flags |= aAxis; };
-	short getAxis() { return flags & 3; };
+	const int getAxis() const { return flags & 3; };
 
 	void setSplit(Float aSplit) { split = aSplit; };
-	Float& getSplit() { return split; };
+	const Float& getSplit() const { return split; };
 
 	void setChildren(KdNode *node) { children = node; assert((flags & 3) == 0); };
-	KdNode* getLeftChild() { return (KdNode*)((off_t)children & ~3); };
-	KdNode* getRightChild() { return (KdNode*)((off_t)children & ~3) + 1; };
+	KdNode* getLeftChild() const { return (KdNode*)((off_t)children & ~3); };
+	KdNode* getRightChild() const { return (KdNode*)((off_t)children & ~3) + 1; };
 
-	ShapeList* getShapes() { return (ShapeList*)((off_t)shapes & ~3); };
+	ShapeList* getShapes() const { return (ShapeList*)((off_t)shapes & ~3); };
 	void addShape(Shape* aShape) { getShapes()->push_back(aShape); };
 };
 
@@ -79,6 +79,7 @@
 	int max_depth;
 
 	void recursive_build(KdNode *node, BBox bbox, int maxdepth);
+	void recursive_load(istream &st, KdNode *node);
 public:
 	KdTree() : Container(), root(NULL), built(false), max_depth(32) {};
 	~KdTree() { if (root) delete root; };
@@ -87,9 +88,12 @@
 		Float &nearest_distance);
 	void optimize() { build(); };
 	void build();
-	void save(ostream &str, KdNode *node = NULL);
-	void load(istream &str, KdNode *node = NULL);
+	const bool isBuilt() const { return built; };
+	KdNode *getRootNode() const { return root; };
 	void setMaxDepth(int md) { max_depth = md; };
+
+	ostream & dump(ostream &st);
+	istream & load(istream &st);
 };
 
 #endif

--- a/include/raytracer.h	Tue Apr 22 13:33:12 2008 +0200
+++ b/include/raytracer.h	Wed Apr 23 10:38:33 2008 +0200
@@ -34,8 +34,6 @@
 #include "container.h"
 #include "scene.h"
 
-using namespace std;
-
 class Raytracer;
 
 struct RenderrowData {

--- a/include/sampler.h	Tue Apr 22 13:33:12 2008 +0200
+++ b/include/sampler.h	Wed Apr 23 10:38:33 2008 +0200
@@ -75,12 +75,12 @@
 class DefaultSampler: public Sampler
 {
 	int phase;
-	int subsample;
+	int subsample;  // 0,1 = no, 1+ = size of sampling grid
 	int oversample; // 0 = no, 1 = 5x, 2 = 9x, 3 = 16x
 	int sx,sy,osa_samp; // current sample properties
 public:
 	DefaultSampler(Float *abuffer, int &aw, int &ah):
-		Sampler(abuffer, aw, ah), phase(-1), subsample(8), oversample(0) {};
+		Sampler(abuffer, aw, ah), phase(-1), subsample(4), oversample(0) {};
 	void init();
 	int initSampleSet();
 	bool nextSample(Sample *s);

--- a/include/scene.h	Tue Apr 22 13:33:12 2008 +0200
+++ b/include/scene.h	Wed Apr 23 10:38:33 2008 +0200
@@ -30,24 +30,15 @@
 #include <vector>
 #include <typeinfo>
 
+#include "common.h"
 #include "sampler.h"
 #include "noise.h"
 #include "vector.h"
 #include "quaternion.h"
 
-/*
-triangle intersection alghoritm
-chooses are:
-TRI_PLUCKER
-TRI_BARI
-TRI_BARI_PRE
-*/
-#if !defined(TRI_PLUCKER) && !defined(TRI_BARI) && !defined(TRI_BARI_PRE)
-#	define TRI_BARI_PRE
-#endif
-
-using namespace std;
-
+/**
+ * ray
+ */
 class Ray
 {
 public:
@@ -164,167 +155,4 @@
 	void setSmooth(bool sm) { smooth = sm; };
 };
 
-/**
- * shape
- */
-class Shape
-{
-public:
-	Material *material;
-	Shape() {};
-	virtual ~Shape() {};
-
-	// first intersection point
-	virtual bool intersect(const Ray &ray, Float &dist) const = 0;
-
-	// all intersections (only for CSG)
-	virtual bool intersect_all(const Ray &ray, Float dist, vector<Float> &allts) const = 0;
-
-	// intersection with AABB
-	virtual bool intersect_bbox(const BBox &bbox) const = 0;
-
-	// normal at point P
-	virtual const Vector3 normal(const Vector3 &P) const = 0;
-
-	virtual BBox get_bbox() const = 0;
-};
-
-/**
- * list of shapes
- */
-class ShapeList: public vector<Shape*>
-{
-};
-
-/**
- * sphere shape
- */
-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) const;
-	bool intersect_all(const Ray &ray, Float dist, vector<Float> &allts) const;
-	bool intersect_bbox(const BBox &bbox) const;
-	const Vector3 normal(const Vector3 &P) const { return (P - center) * inv_radius; };
-	BBox get_bbox() const;
-};
-
-/**
- * box shape
- */
-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) const;
-	bool intersect_all(const Ray &ray, Float dist, vector<Float> &allts) const { return false; };
-	bool intersect_bbox(const BBox &bbox) const;
-	const Vector3 normal(const Vector3 &P) const;
-	BBox get_bbox() const { return BBox(L, H); };
-};
-
-/**
- * triangle vertex
- */
-class Vertex
-{
-public:
-	Vector3 P;
-	Vertex(const Vector3 &aP): P(aP) {};
-};
-
-/**
- * triangle vertex with normal
- */
-class NormalVertex: public Vertex
-{
-public:
-	Vector3 N;
-	NormalVertex(const NormalVertex *v): Vertex(v->P), N(v->N) {};
-	NormalVertex(const Vector3 &aP): Vertex(aP) {};
-	NormalVertex(const Vector3 &aP, const Vector3 &aN): Vertex(aP), N(aN) {};
-	const Vector3 &getNormal() { return N; };
-	void setNormal(const Vector3 &aN) { N = aN; };
-};
-
-/**
- * triangle shape
- */
-class Triangle: public Shape
-{
-#ifdef TRI_BARI_PRE
-	Float nu, nv, nd;
-	int k; // dominant axis
-	Float bnu, bnv;
-	Float cnu, cnv;
 #endif
-#ifdef TRI_BARI
-	int k; // dominant axis
-#endif
-#ifdef TRI_PLUCKER
-	Float pla[6], plb[6], plc[6];
-#endif
-	Vector3 N;
-	const Vector3 smooth_normal(const Vector3 &P) const
-	{
-#ifdef TRI_BARI_PRE
-		const Vector3 &NA = static_cast<NormalVertex*>(A)->N;
-		const Vector3 &NB = static_cast<NormalVertex*>(B)->N;
-		const Vector3 &NC = static_cast<NormalVertex*>(C)->N;
-		static const int modulo3[5] = {0,1,2,0,1};
-		register const int ku = modulo3[k+1];
-		register const int kv = modulo3[k+2];
-		const Float pu = P[ku] - A->P[ku];
-		const Float pv = P[kv] - A->P[kv];
-		const Float u = pv * bnu + pu * bnv;
-		const Float v = pu * cnv + pv * cnu;
-		Vector3 n = NA + u * (NB - NA) + v * (NC - NA);
-		n.normalize();
-		return n;
-#else
-		return N; // not implemented for other algorithms
-#endif
-	};
-public:
-	Vertex *A, *B, *C;
-
-	Triangle() {};
-	Triangle(Vertex *aA, Vertex *aB, Vertex *aC, Material *amaterial);
-	bool intersect(const Ray &ray, Float &dist) const;
-	bool intersect_all(const Ray &ray, Float dist, vector<Float> &allts) const {return false;};
-	bool intersect_bbox(const BBox &bbox) const;
-	const Vector3 normal(const Vector3 &P) const { return (material->smooth ? smooth_normal(P) : N); };
-	const Vector3 getNormal() const { return N; };
-	BBox get_bbox() const;
-};
-
-template <class T> class Array
-{
-	T *array;
-public:
-	Array(int n) { array = new T[n]; };
-	~Array() { delete[] array; };
-	const T &operator[](int i) const { return array[i]; };
-};
-
-typedef Array<Vertex> VertexArray;
-typedef Array<NormalVertex> NormalVertexArray;
-typedef Array<Triangle> TriangleArray;
-
-#endif

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/serialize.h	Wed Apr 23 10:38:33 2008 +0200
@@ -0,0 +1,57 @@
+/*
+ * serialize.h: object serialization functions
+ *
+ * This file is part of Pyrit Ray Tracer.
+ *
+ * Copyright 2008  Radek Brich
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ * THE SOFTWARE.
+ */
+
+#ifndef SERIALIZE_H
+#define SERIALIZE_H
+
+#include <map>
+
+#include "shapes.h"
+#include "container.h"
+#include "kdtree.h"
+
+class Indexer
+{
+	map <void *, int> indexmap;
+	int index;
+public:
+	Indexer(): indexmap(), index(0) {};
+	void reset() { indexmap.clear(); index = 0; };
+	bool get(void *o, int &retidx);
+	const int &operator[](void *o) { return indexmap[o]; };
+};
+
+extern Indexer vertex_index, shape_index;
+
+void resetSerializer();
+Shape *loadShape(istream &st);
+
+ostream & operator<<(ostream &st, Shape &o);
+ostream & operator<<(ostream &st, Vertex &o);
+ostream & operator<<(ostream &st, Container &o);
+istream & operator>>(istream &st, Vector3 &v);
+
+#endif

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/shapes.h	Wed Apr 23 10:38:33 2008 +0200
@@ -0,0 +1,217 @@
+/*
+ * shapes.h: shape classes
+ *
+ * This file is part of Pyrit Ray Tracer.
+ *
+ * Copyright 2006, 2007, 2008  Radek Brich
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ * THE SOFTWARE.
+ */
+
+#ifndef SHAPES_H
+#define SHAPES_H
+
+#include "common.h"
+#include "scene.h"
+
+/*
+triangle intersection alghoritm
+options are:
+TRI_PLUCKER
+TRI_BARI
+TRI_BARI_PRE
+*/
+#if !defined(TRI_PLUCKER) && !defined(TRI_BARI) && !defined(TRI_BARI_PRE)
+#	define TRI_BARI_PRE
+#endif
+
+/**
+ * shape
+ */
+class Shape
+{
+public:
+	Material *material;
+	Shape() {};
+	virtual ~Shape() {};
+
+	// first intersection point
+	virtual bool intersect(const Ray &ray, Float &dist) const = 0;
+
+	// all intersections (only for CSG)
+	virtual bool intersect_all(const Ray &ray, Float dist, vector<Float> &allts) const = 0;
+
+	// intersection with AABB
+	virtual bool intersect_bbox(const BBox &bbox) const = 0;
+
+	// normal at point P
+	virtual const Vector3 normal(const Vector3 &P) const = 0;
+
+	virtual BBox get_bbox() const = 0;
+
+	virtual ostream & dump(ostream &st) const = 0;
+};
+
+/**
+ * list of shapes
+ */
+typedef vector<Shape*> ShapeList;
+
+/**
+ * sphere shape
+ */
+class Sphere: public Shape
+{
+	Vector3 center;
+	Float radius;
+
+	Float sqr_radius;
+	Float inv_radius;
+public:
+	Sphere(const Vector3 &acenter, const Float aradius, Material *amaterial):
+		center(acenter), radius(aradius),
+		sqr_radius(aradius*aradius), inv_radius(1.0f/aradius)
+		{ material = amaterial; }
+	bool intersect(const Ray &ray, Float &dist) const;
+	bool intersect_all(const Ray &ray, Float dist, vector<Float> &allts) const;
+	bool intersect_bbox(const BBox &bbox) const;
+	const Vector3 normal(const Vector3 &P) const { return (P - center) * inv_radius; };
+	BBox get_bbox() const;
+	const Vector3 getCenter() const { return center; };
+	const Float getRadius() const { return radius; };
+	ostream & dump(ostream &st) const;
+};
+
+/**
+ * box shape
+ */
+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) const;
+	bool intersect_all(const Ray &ray, Float dist, vector<Float> &allts) const { return false; };
+	bool intersect_bbox(const BBox &bbox) const;
+	const Vector3 normal(const Vector3 &P) const;
+	BBox get_bbox() const { return BBox(L, H); };
+	const Vector3 getL() const { return L; };
+	const Vector3 getH() const { return H; };
+	ostream & dump(ostream &st) const;
+};
+
+/**
+ * triangle vertex
+ */
+class Vertex
+{
+public:
+	Vector3 P;
+	Vertex(const Vector3 &aP): P(aP) {};
+	virtual ostream & dump(ostream &st) const;
+};
+
+/**
+ * triangle vertex with normal
+ */
+class NormalVertex: public Vertex
+{
+public:
+	Vector3 N;
+	NormalVertex(const NormalVertex *v): Vertex(v->P), N(v->N) {};
+	NormalVertex(const Vector3 &aP): Vertex(aP) {};
+	NormalVertex(const Vector3 &aP, const Vector3 &aN): Vertex(aP), N(aN) {};
+	const Vector3 &getNormal() { return N; };
+	void setNormal(const Vector3 &aN) { N = aN; };
+	ostream & dump(ostream &st) const;
+};
+
+/**
+ * triangle shape
+ */
+class Triangle: public Shape
+{
+#ifdef TRI_BARI_PRE
+	Float nu, nv, nd;
+	int k; // dominant axis
+	Float bnu, bnv;
+	Float cnu, cnv;
+#endif
+#ifdef TRI_BARI
+	int k; // dominant axis
+#endif
+#ifdef TRI_PLUCKER
+	Float pla[6], plb[6], plc[6];
+#endif
+	Vector3 N;
+	const Vector3 smooth_normal(const Vector3 &P) const
+	{
+#ifdef TRI_BARI_PRE
+		const Vector3 &NA = static_cast<NormalVertex*>(A)->N;
+		const Vector3 &NB = static_cast<NormalVertex*>(B)->N;
+		const Vector3 &NC = static_cast<NormalVertex*>(C)->N;
+		static const int modulo3[5] = {0,1,2,0,1};
+		register const int ku = modulo3[k+1];
+		register const int kv = modulo3[k+2];
+		const Float pu = P[ku] - A->P[ku];
+		const Float pv = P[kv] - A->P[kv];
+		const Float u = pv * bnu + pu * bnv;
+		const Float v = pu * cnv + pv * cnu;
+		Vector3 n = NA + u * (NB - NA) + v * (NC - NA);
+		n.normalize();
+		return n;
+#else
+		return N; // not implemented for other algorithms
+#endif
+	};
+public:
+	Vertex *A, *B, *C;
+
+	Triangle() {};
+	Triangle(Vertex *aA, Vertex *aB, Vertex *aC, Material *amaterial);
+	bool intersect(const Ray &ray, Float &dist) const;
+	bool intersect_all(const Ray &ray, Float dist, vector<Float> &allts) const {return false;};
+	bool intersect_bbox(const BBox &bbox) const;
+	const Vector3 normal(const Vector3 &P) const { return (material->smooth ? smooth_normal(P) : N); };
+	const Vector3 getNormal() const { return N; };
+	BBox get_bbox() const;
+	ostream & dump(ostream &st) const;
+};
+
+template <class T> class Array
+{
+	T *array;
+public:
+	Array(int n) { array = new T[n]; };
+	~Array() { delete[] array; };
+	const T &operator[](int i) const { return array[i]; };
+};
+
+typedef Array<Vertex> VertexArray;
+typedef Array<NormalVertex> NormalVertexArray;
+typedef Array<Triangle> TriangleArray;
+
+#endif

--- a/include/vector.h	Tue Apr 22 13:33:12 2008 +0200
+++ b/include/vector.h	Wed Apr 23 10:38:33 2008 +0200
@@ -30,6 +30,8 @@
 #include <math.h>
 #include <iostream>
 
+#include "common.h"
+
 using namespace std;
 
 /**
@@ -166,10 +168,24 @@
 		return Vector3(a.x * b.x, a.y * b.y, a.z * b.z);
 	};
 
-	// print
+	// write
 	friend ostream & operator<<(ostream &st, const Vector3 &v)
 	{
-		return st << "(" << v.x << ", " << v.y  << ", " << v.z << ")";
+		return st << "(" << v.x << "," << v.y  << "," << v.z << ")";
+	};
+
+	// read
+	friend istream & operator>>(istream &st, Vector3 &v)
+	{
+		char s[10];
+		st.getline(s, 10, '(');
+		st >> v.x;
+		st.getline(s, 10, ',');
+		st >> v.y;
+		st.getline(s, 10, ',');
+		st >> v.z;
+		st.getline(s, 10, ')');
+		return st;
 	};
 };
 

--- a/src/SConscript	Tue Apr 22 13:33:12 2008 +0200
+++ b/src/SConscript	Wed Apr 23 10:38:33 2008 +0200
@@ -16,8 +16,9 @@
 	pyenv.ParseConfig('python-config --includes --libs')
 
 sources = [
-	'raytracer.cc', 'scene.cc', 'sampler.cc',
-	'container.cc', 'kdtree.cc', 'octree.cc', 'noise.cc']
+	'raytracer.cc', 'scene.cc', 'shapes.cc', 'sampler.cc',
+	'container.cc', 'kdtree.cc', 'octree.cc', 'noise.cc',
+	'serialize.cc']
 
 objs = []
 shared_objs = []

--- a/src/container.cc	Tue Apr 22 13:33:12 2008 +0200
+++ b/src/container.cc	Wed Apr 23 10:38:33 2008 +0200
@@ -26,11 +26,11 @@
 
 #include "common.h"
 #include "container.h"
+#include "serialize.h"
 
 void Container::addShape(Shape* aShape)
 {
-	const Float e = 10*Eps;
-	shapes.push_back(aShape);
+	const Float e = Eps;
 	if (shapes.size() == 0) {
 		/* initialize bounding box */
 		bbox = aShape->get_bbox();
@@ -46,6 +46,7 @@
 		if (shapebb.H.y + e > bbox.H.y)  bbox.H.y = shapebb.H.y + e;
 		if (shapebb.H.z + e > bbox.H.z)  bbox.H.z = shapebb.H.z + e;
 	}
+	shapes.push_back(aShape);
 };
 
 Shape *Container::nearest_intersection(const Shape *origin_shape, const Ray &ray,
@@ -58,3 +59,16 @@
 			nearest_shape = *shape;
 	return nearest_shape;
 }
+
+ostream & Container::dump(ostream &st)
+{
+	st << "(container," << shapes.size();
+	ShapeList::iterator shape;
+	for (shape = shapes.begin(); shape != shapes.end(); shape++)
+	{
+		int idx;
+		if (!shape_index.get(*shape, idx))
+			st << "," << **shape;
+	}
+	return st << ")";
+}

--- a/src/kdtree.cc	Tue Apr 22 13:33:12 2008 +0200
+++ b/src/kdtree.cc	Wed Apr 23 10:38:33 2008 +0200
@@ -26,9 +26,11 @@
 
 #include <algorithm>
 #include <stack>
+#include <string>
+#include <sstream>
 
-#include "common.h"
 #include "kdtree.h"
+#include "serialize.h"
 
 class ShapeBound
 {
@@ -265,7 +267,7 @@
 			}
 			else
 			{ /* (stack[entry].pb[axis] > splitVal) */
-				if (splitVal < stack[exit].pb[axis])
+				if (stack[exit].pb[axis] > splitVal)
 				{
 					/* case P1, P2, P3, and N5 */
 					node = node->getRightChild();
@@ -305,7 +307,7 @@
 		ShapeList::iterator shape;
 		for (shape = node->getShapes()->begin(); shape != node->getShapes()->end(); shape++)
 			if (*shape != origin_shape && (*shape)->intersect(ray, dist)
-			&& dist >= stack[entry].t)
+			&& dist >= stack[entry].t - Eps)
 			{
 				nearest_shape = *shape;
 				nearest_distance = dist;
@@ -326,27 +328,118 @@
 	return NULL;
 }
 
-// this should save whole kd-tree with triangles distributed into leaves
-void KdTree::save(ostream &str, KdNode *node)
+ostream & operator<<(ostream &st, KdNode &node)
+{
+	if (node.isLeaf())
+	{
+		st << "(l," << node.getShapes()->size();
+		ShapeList::iterator shape;
+		for (shape = node.getShapes()->begin(); shape != node.getShapes()->end(); shape++)
+			st << "," << shape_index[*shape];
+		st << ")";
+	}
+	else
+	{
+		st << "(s," << (char)('x'+node.getAxis()) << "," << node.getSplit() << ",";
+		st << *node.getLeftChild() << ",";
+		st << *node.getRightChild() << ")";
+	}
+	return st;
+}
+
+ostream & KdTree::dump(ostream &st)
 {
 	if (!built)
-		return;
-	if (node == NULL)
-		node = root;
-	if (node->isLeaf())
-		str << "(leaf: " << node->getShapes()->size() << " shapes)";
-	else
+		return Container::dump(st);
+
+	st << "(kdtree," << shapes.size();
+	ShapeList::iterator shape;
+	for (shape = shapes.begin(); shape != shapes.end(); shape++)
+	{
+		int idx;
+		if (!shape_index.get(*shape, idx))
+			st << "," << endl << **shape;
+	}
+	return st << "," << endl << *getRootNode() << ")";
+}
+
+void KdTree::recursive_load(istream &st, KdNode *node)
+{
+	string s;
+	istringstream is;
+
+	getline(st, s, ',');
+	trim(s);
+
+	if (s.compare("(s") == 0)
 	{
-		str << "(split " << (char)('x'+node->getAxis()) << " at " << node->getSplit() << "; L=";
-		save(str, node->getLeftChild());
-		str << "; R=";
-		save(str, node->getRightChild());
-		str << ";)";
+		// split
+		int axis;
+		Float split;
+
+		delete node->getShapes();
+		node->setChildren(new KdNode[2]);
+
+		getline(st, s, ',');
+		axis = s.c_str()[0]-'x';
+		node->setAxis(axis);
+
+		st >> split;
+		getline(st, s, ',');
+		node->setSplit(split);
+
+		recursive_load(st, node->getLeftChild());
+		getline(st, s, ',');
+		recursive_load(st, node->getRightChild());
+		getline(st, s, ')');
+	}
+
+	if (s.compare("(l") == 0)
+	{
+		// leaf
+		int count, idx;
+
+		node->setLeaf();
+
+		st >> count;
+		for (int i = 0; i < count; i++)
+		{
+			getline(st, s, ',');
+			st >> idx;
+			node->addShape(shapes[idx]);
+		}
+		getline(st, s, ')');
 	}
 }
 
-// load kd-tree from file/stream
-void KdTree::load(istream &str, KdNode *node)
+istream & KdTree::load(istream &st)
 {
+	string s;
+	istringstream is;
 
+	getline(st, s, ',');
+	if (s.compare("(kdtree") != 0)
+		return st;
+
+	shapes.clear();
+	if (root) delete root;
+	root = new KdNode();
+
+	getline(st, s, ',');
+	int shape_count;
+	is.str(s);
+	is >> shape_count;
+
+	Shape *shape;
+	for (int i = 0; i < shape_count; i++)
+	{
+		shape = loadShape(st);
+		Container::addShape(shape);
+		getline(st, s, ',');
+	}
+
+	recursive_load(st, root);
+
+	built = true;
+	return st;
 }

--- a/src/raytracermodule.cc	Tue Apr 22 13:33:12 2008 +0200
+++ b/src/raytracermodule.cc	Wed Apr 23 10:38:33 2008 +0200
@@ -28,7 +28,6 @@
 
 #include <vector>
 #include "raytracer.h"
-#include "scene.h"
 #include "octree.h"
 
 //=========================== Light Source Object ===========================

--- a/src/scene.cc	Tue Apr 22 13:33:12 2008 +0200
+++ b/src/scene.cc	Wed Apr 23 10:38:33 2008 +0200
@@ -24,9 +24,6 @@
  * THE SOFTWARE.
  */
 
-#include <math.h>
-
-#include "common.h"
 #include "scene.h"
 
 void Camera::rotate(const Quaternion &q)
@@ -110,474 +107,3 @@
 	b = tfar;
 	return true;
 }
-
-bool Sphere::intersect(const Ray &ray, Float &dist) const
-{
-	Vector3 V = ray.o - center;
-	register Float d = -dot(V, ray.dir);
-	register Float Det = d * d - (dot(V,V) - sqr_radius);
-	register Float t1,t2;
-	if (Det > 0) {
-		Det = sqrtf(Det);
-		t1 = d - Det;
-		t2 = d + Det;
-		if (t1 > 0)
-		{
-			if (t1 < dist)
-			{
-				dist = t1;
-				return true;
-			}
-		}
-		else if (t2 > 0 && t2 < dist)
-		{
-			dist = t2;
-			return true;
-		}
-	}
-	return false;
-}
-
-/* if there should be CSG sometimes, this may be needed... */
-bool Sphere::intersect_all(const Ray &ray, Float dist, vector<Float> &allts) const
-{
-	//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);
-
-	if (Det > 0) {
-		Det = sqrtf(Det);
-		Float t1 = Vd - Det;
-		Float t2 = Vd + Det;
-		if (t1 < 0)
-		{
-			if (t2 > 0)
-			{
-				// ray from inside of the sphere
-				allts.push_back(0.0);
-				allts.push_back(t2);
-				return true;
-			}
-			else
-				return false;
-		}
-		else
-		{
-			allts.push_back(t1);
-			allts.push_back(t2);
-			return true;
-		}
-	}
-	return false;
-}
-
-bool Sphere::intersect_bbox(const BBox &bbox) const
-{
-	register float dmin = 0;
-	for (int i = 0; i < 3; i++)
-	{
-		if (center[i] < bbox.L[i])
-			dmin += (center[i] - bbox.L[i])*(center[i] - bbox.L[i]);
-		else
-		if (center[i] > bbox.H[i])
-			dmin += (center[i] - bbox.H[i])*(center[i] - bbox.H[i]);
-	}
-	if (dmin <= sqr_radius)
-		return true;
-	return false;
-};
-
-BBox Sphere::get_bbox() const
-{
-	return BBox(center - radius, center + radius);
-}
-
-bool Box::intersect(const Ray &ray, Float &dist) const
-{
-	register Float tnear = -Inf;
-	register Float tfar = Inf;
-	register Float t1, t2;
-
-	for (int i = 0; i < 3; i++)
-	{
-		if (ray.dir[i] == 0) {
-			/* ray is parallel to these planes */
-			if (ray.o[i] < L[i] || ray.o[i] > H[i])
-				return false;
-		}
-		else
-		{
-			/* compute the intersection distance of the planes */
-			t1 = (L[i] - ray.o[i]) / ray.dir[i];
-			t2 = (H[i] - ray.o[i]) / ray.dir[i];
-
-			if (t1 > t2)
-				swap(t1, t2);
-
-			if (t1 > tnear)
-				tnear = t1; /* want largest Tnear */
-			if (t2 < tfar)
-				tfar = t2; /* want smallest Tfar */
-			if (tnear > tfar || tfar < 0)
-				return false; /* box missed; box is behind ray */
-		}
-	}
-
-	if (tnear < dist)
-	{
-		dist = tnear;
-		return true;
-	}
-	return false;
-}
-
-bool Box::intersect_bbox(const BBox &bbox) const
-{
-	return (
-	H.x > bbox.L.x && L.x < bbox.H.x &&
-	H.y > bbox.L.y && L.y < bbox.H.y &&
-	H.z > bbox.L.z && L.z < bbox.H.z);
-}
-
-const Vector3 Box::normal(const Vector3 &P) const
-{
-	register Vector3 l = P - L;
-	register Vector3 h = H - P;
-
-	if (l.x < h.x)
-		h.x = -1;
-	else
-	{
-		l.x = h.x;
-		h.x = +1;
-	}
-
-	if (l.y < h.y)
-		h.y = -1;
-	else
-	{
-		l.y = h.y;
-		h.y = +1;
-	}
-
-	if (l.z < h.z)
-		h.z = -1;
-	else
-	{
-		l.z = h.z;
-		h.z = +1;
-	}
-
-	if (l.x > l.y)
-	{
-		h.x = 0;
-		if (l.y > l.z)
-			h.y = 0;
-		else
-			h.z = 0;
-	}
-	else
-	{
-		h.y = 0;
-		if (l.x > l.z)
-			h.x = 0;
-		else
-			h.z = 0;
-	}
-	return h;
-}
-
-#ifdef TRI_PLUCKER
-inline void Plucker(const Vector3 &p, const Vector3 &q, Float* pl)
-{
-    pl[0] = p.x*q.y - q.x*p.y;
-    pl[1] = p.x*q.z - q.x*p.z;
-    pl[2] = p.x - q.x;
-    pl[3] = p.y*q.z - q.y*p.z;
-    pl[4] = p.z - q.z;
-    pl[5] = q.y - p.y;
-}
-
-inline Float Side(const Float* pla, const Float* plb)
-{
-    return pla[0]*plb[4] + pla[1]*plb[5] + pla[2]*plb[3] + pla[4]*plb[0] + pla[5]*plb[1] + pla[3]*plb[2];
-}
-#endif
-
-Triangle::Triangle(Vertex *aA, Vertex *aB, Vertex *aC, Material *amaterial)
-	: A(aA), B(aB), C(aC)
-{
-	material = amaterial;
-
-	const Vector3 c = B->P - A->P;
-	const Vector3 b = C->P - A->P;
-
-	N = cross(c, b);
-	N.normalize();
-
-#ifdef TRI_PLUCKER
-	Plucker(B->P,C->P,pla);
-	Plucker(C->P,A->P,plb);
-	Plucker(A->P,B->P,plc);
-#endif
-
-#if defined(TRI_BARI) || defined(TRI_BARI_PRE)
-	if (fabsf(N.x) > fabsf(N.y))
-	{
-		if (fabsf(N.x) > fabsf(N.z))
-			k = 0;
-		else
-			k = 2;
-	}
-	else
-	{
-		if (fabsf(N.y) > fabsf(N.z))
-			k = 1;
-		else
-			k = 2;
-	}
-#endif
-#ifdef TRI_BARI_PRE
-	int u = (k + 1) % 3;
-	int v = (k + 2) % 3;
-
-	Float krec = 1.0 / N[k];
-	nu = N[u] * krec;
-	nv = N[v] * krec;
-	nd = dot(N, A->P) * krec;
-
-	// first line equation
-	Float reci = 1.0f / (b[u] * c[v] - b[v] * c[u]);
-	bnu = b[u] * reci;
-	bnv = -b[v] * reci;
-
-	// second line equation
-	cnu = -c[u] * reci;
-	cnv = c[v] * reci;
-#endif
-}
-
-bool Triangle::intersect(const Ray &ray, Float &dist) const
-{
-#ifdef TRI_PLUCKER
-	Float plr[6];
-	Plucker(ray.o, ray.o+ray.dir, plr);
-	const bool side0 = Side(plr, pla) >= 0.0;
-	const bool side1 = Side(plr, plb) >= 0.0;
-	if (side0 != side1)
-		return false;
-	const bool side2 = Side(plr, plc) >= 0.0;
-	if (side0 != side2)
-		return false;
-	const Float t = - dot( (ray.o - A->P), N) / dot(ray.dir,N);
-	if(t <= Eps || t >= dist)
-		return false;
-	dist = t;
-	return true;
-#endif
-
-#if defined(TRI_BARI) || defined(TRI_BARI_PRE)
-	static const int modulo3[5] = {0,1,2,0,1};
-	const Vector3 &O = ray.o;
-	const Vector3 &D = ray.dir;
-	register const int u = modulo3[k+1];
-	register const int v = modulo3[k+2];
-#endif
-#ifdef TRI_BARI_PRE
-	const Float t = (nd - O[k] - nu * O[u] - nv * O[v]) / (D[k] + nu * D[u] + nv * D[v]);
-
-	if (t >= dist || t < Eps)
-		return false;
-
-	const Float hu = O[u] + t * D[u] - A->P[u];
-	const Float hv = O[v] + t * D[v] - A->P[v];
-	const Float beta = hv * bnu + hu * bnv;
-
-	if (beta < 0.)
-		return false;
-
-	const Float gamma = hu * cnv + hv * cnu;
-	if (gamma < 0. || beta + gamma > 1.)
-		return false;
-
-	dist = t;
-	return true;
-#endif
-
-#ifdef TRI_BARI
-	// original barycentric coordinates based intesection
-	// not optimized, just for reference
-	const Vector3 c = B - A;
-	const Vector3 b = C - A;
-	// distance test
-	const Float t = - dot( (O-A), N) / dot(D,N);
-	if (t < Eps || t > dist)
-		return false;
-
-	// calc hitpoint
-	const Float Hu = O[u] + t * D[u] - A[u];
-	const Float Hv = O[v] + t * D[v] - A[v];
-	const Float beta = (b[u] * Hv - b[v] * Hu) / (b[u] * c[v] - b[v] * c[u]);
-	if (beta < 0)
-		return false;
-	const Float gamma = (c[v] * Hu - c[u] * Hv) / (b[u] * c[v] - b[v] * c[u]);
-	if (gamma < 0)
-		return false;
-	if (beta+gamma > 1)
-		return false;
-	dist = t;
-	return true;
-#endif
-}
-
-bool Triangle::intersect_bbox(const BBox &bbox) const
-{
-	const Vector3 boxcenter = (bbox.L+bbox.H)*0.5;
-	const Vector3 boxhalfsize = (bbox.H-bbox.L)*0.5;
-	const Vector3 v0 = A->P - boxcenter;
-	const Vector3 v1 = B->P - boxcenter;
-	const Vector3 v2 = C->P - boxcenter;
-	const Vector3 e0 = v1-v0;
-	const Vector3 e1 = v2-v1;
-	const Vector3 e2 = v0-v2;
-
-	Float fex = fabsf(e0.x);
-	Float fey = fabsf(e0.y);
-	Float fez = fabsf(e0.z);
-
-	Float p0,p1,p2,min,max,rad;
-
-	p0 = e0.z*v0.y - e0.y*v0.z;
-	p2 = e0.z*v2.y - e0.y*v2.z;
-	if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
-	rad = fez * boxhalfsize.y + fey * boxhalfsize.z;
-	if(min>rad || max<-rad) return false;
-
-	p0 = -e0.z*v0.x + e0.x*v0.z;
-	p2 = -e0.z*v2.x + e0.x*v2.z;
-	if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
-	rad = fez * boxhalfsize.x + fex * boxhalfsize.z;
-	if(min>rad || max<-rad) return false;
-
-	p1 = e0.y*v1.x - e0.x*v1.y;
-	p2 = e0.y*v2.x - e0.x*v2.y;
-	if(p2<p1) {min=p2; max=p1;} else {min=p1; max=p2;}
-	rad = fey * boxhalfsize.x + fex * boxhalfsize.y;
-	if(min>rad || max<-rad) return false;
-
-	fex = fabsf(e1.x);
-	fey = fabsf(e1.y);
-	fez = fabsf(e1.z);
-
-	p0 = e1.z*v0.y - e1.y*v0.z;
-	p2 = e1.z*v2.y - e1.y*v2.z;
-	if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
-	rad = fez * boxhalfsize.y + fey * boxhalfsize.z;
-	if(min>rad || max<-rad) return false;
-
-	p0 = -e1.z*v0.x + e1.x*v0.z;
-	p2 = -e1.z*v2.x + e1.x*v2.z;
-	if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
-	rad = fez * boxhalfsize.x + fex * boxhalfsize.z;
-	if(min>rad || max<-rad) return false;
-
-	p0 = e1.y*v0.x - e1.x*v0.y;
-	p1 = e1.y*v1.x - e1.x*v1.y;
-	if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;}
-	rad = fey * boxhalfsize.x + fex * boxhalfsize.y;
-	if(min>rad || max<-rad) return false;
-
-	fex = fabsf(e2.x);
-	fey = fabsf(e2.y);
-	fez = fabsf(e2.z);
-
-	p0 = e2.z*v0.y - e2.y*v0.z;
-	p1 = e2.z*v1.y - e2.y*v1.z;
-	if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;}
-	rad = fez * boxhalfsize.y + fey * boxhalfsize.z;
-	if(min>rad || max<-rad) return false;
-
-	p0 = -e2.z*v0.x + e2.x*v0.z;
-	p1 = -e2.z*v1.x + e2.x*v1.z;
-	if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;}
-	rad = fez * boxhalfsize.x + fex * boxhalfsize.z;
-	if(min>rad || max<-rad) return false;
-
-	p1 = e2.y*v1.x - e2.x*v1.y;
-	p2 = e2.y*v2.x - e2.x*v2.y;
-	if(p2<p1) {min=p2; max=p1;} else {min=p1; max=p2;}
-	rad = fey * boxhalfsize.x + fex * boxhalfsize.y;
-	if(min>rad || max<-rad) return false;
-
-	/* test overlap in the {x,y,z}-directions */
-	/* test in X-direction */
-	min = v0.x;
-	if (v1.x < min) min = v1.x;
-	if (v2.x < min) min = v2.x;
-	max = v0.x;
-	if (v1.x > max) max = v1.x;
-	if (v2.x > max) max = v2.x;
-	if(min>boxhalfsize.x || max<-boxhalfsize.x) return false;
-
-	/* test in Y-direction */
-	min = v0.y;
-	if (v1.y < min) min = v1.y;
-	if (v2.y < min) min = v2.y;
-	max = v0.y;
-	if (v1.y > max) max = v1.y;
-	if (v2.y > max) max = v2.y;
-	if(min>boxhalfsize.y || max<-boxhalfsize.y) return false;
-
-	/* test in Z-direction */
-	min = v0.z;
-	if (v1.z < min) min = v1.z;
-	if (v2.z < min) min = v2.z;
-	max = v0.z;
-	if (v1.z > max) max = v1.z;
-	if (v2.z > max) max = v2.z;
-	if(min>boxhalfsize.z || max<-boxhalfsize.z) return false;
-
-	/*  test if the box intersects the plane of the triangle */
-	Vector3 vmin,vmax;
-	Float v;
-	for(int q=0;q<3;q++)
-	{
-		v=v0[q];
-		if(N[q]>0.0f)
-		{
-			vmin.cell[q]=-boxhalfsize[q] - v;
-			vmax.cell[q]= boxhalfsize[q] - v;
-		}
-		else
-		{
-			vmin.cell[q]= boxhalfsize[q] - v;
-			vmax.cell[q]=-boxhalfsize[q] - v;
-		}
-	}
-	if(dot(N,vmin)>0.0f) return false;
-	if(dot(N,vmax)>=0.0f) return true;
-
-	return false;
-}
-
-BBox Triangle::get_bbox() const
-{
-	BBox bbox = BBox();
-	bbox.L = A->P;
-	if (B->P.x < bbox.L.x)  bbox.L.x = B->P.x;
-	if (C->P.x < bbox.L.x)  bbox.L.x = C->P.x;
-	if (B->P.y < bbox.L.y)  bbox.L.y = B->P.y;
-	if (C->P.y < bbox.L.y)  bbox.L.y = C->P.y;
-	if (B->P.z < bbox.L.z)  bbox.L.z = B->P.z;
-	if (C->P.z < bbox.L.z)  bbox.L.z = C->P.z;
-	bbox.H = A->P;
-	if (B->P.x > bbox.H.x)  bbox.H.x = B->P.x;
-	if (C->P.x > bbox.H.x)  bbox.H.x = C->P.x;
-	if (B->P.y > bbox.H.y)  bbox.H.y = B->P.y;
-	if (C->P.y > bbox.H.y)  bbox.H.y = C->P.y;
-	if (B->P.z > bbox.H.z)  bbox.H.z = B->P.z;
-	if (C->P.z > bbox.H.z)  bbox.H.z = C->P.z;
-	return bbox;
-};

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/serialize.cc	Wed Apr 23 10:38:33 2008 +0200
@@ -0,0 +1,97 @@
+/*
+ * serialize.cc: object serialization functions
+ *
+ * This file is part of Pyrit Ray Tracer.
+ *
+ * Copyright 2008  Radek Brich
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ * THE SOFTWARE.
+ */
+
+#include "serialize.h"
+#include <string>
+#include <sstream>
+
+Indexer vertex_index, shape_index;
+
+void resetSerializer()
+{
+	vertex_index.reset();
+	shape_index.reset();
+}
+
+bool Indexer::get(void *o, int &retidx)
+{
+	map <void *, int>::iterator i;
+	i = indexmap.find(o);
+	if (i == indexmap.end())
+	{
+		retidx = index++;
+		indexmap[o] = retidx;
+		return false;
+	}
+	else
+	{
+		retidx = i->second;
+		return true;
+	}
+}
+
+Shape *loadShape(istream &st)
+{
+	string s;
+	istringstream is;
+	getline(st, s, ',');
+	trim(s);
+	if (s.compare("(box") == 0)
+	{
+		Vector3 L,H;
+		st >> L;
+		getline(st, s, ',');
+		st >> H;
+		getline(st, s, ')');
+		return new Box(L, H, new Material(Colour(1,1,1)));
+	}
+	if (s.compare("(sphere") == 0)
+	{
+		Vector3 center;
+		Float radius;
+		st >> center;
+		getline(st, s, ',');
+		st >> radius;
+		getline(st, s, ')');
+		return new Sphere(center, radius, new Material(Colour(1,1,1)));
+	}
+	return NULL;
+}
+
+ostream & operator<<(ostream &st, Shape &o)
+{
+	return o.dump(st);
+}
+
+ostream & operator<<(ostream &st, Vertex &o)
+{
+	return o.dump(st);
+}
+
+ostream & operator<<(ostream &st, Container &o)
+{
+	return o.dump(st);
+}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/shapes.cc	Wed Apr 23 10:38:33 2008 +0200
@@ -0,0 +1,531 @@
+/*
+ * shapes.cc: shape classes
+ *
+ * This file is part of Pyrit Ray Tracer.
+ *
+ * Copyright 2006, 2007, 2008  Radek Brich
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ * THE SOFTWARE.
+ */
+
+#include "shapes.h"
+#include "serialize.h"
+
+bool Sphere::intersect(const Ray &ray, Float &dist) const
+{
+	Vector3 V = ray.o - center;
+	register Float d = -dot(V, ray.dir);
+	register Float Det = d * d - (dot(V,V) - sqr_radius);
+	register Float t1,t2;
+	if (Det > 0) {
+		Det = sqrtf(Det);
+		t1 = d - Det;
+		t2 = d + Det;
+		if (t1 > 0)
+		{
+			if (t1 < dist)
+			{
+				dist = t1;
+				return true;
+			}
+		}
+		else if (t2 > 0 && t2 < dist)
+		{
+			dist = t2;
+			return true;
+		}
+	}
+	return false;
+}
+
+/* if there should be CSG sometimes, this may be needed... */
+bool Sphere::intersect_all(const Ray &ray, Float dist, vector<Float> &allts) const
+{
+	//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);
+
+	if (Det > 0) {
+		Det = sqrtf(Det);
+		Float t1 = Vd - Det;
+		Float t2 = Vd + Det;
+		if (t1 < 0)
+		{
+			if (t2 > 0)
+			{
+				// ray from inside of the sphere
+				allts.push_back(0.0);
+				allts.push_back(t2);
+				return true;
+			}
+			else
+				return false;
+		}
+		else
+		{
+			allts.push_back(t1);
+			allts.push_back(t2);
+			return true;
+		}
+	}
+	return false;
+}
+
+bool Sphere::intersect_bbox(const BBox &bbox) const
+{
+	register float dmin = 0;
+	for (int i = 0; i < 3; i++)
+	{
+		if (center[i] < bbox.L[i])
+			dmin += (center[i] - bbox.L[i])*(center[i] - bbox.L[i]);
+		else
+		if (center[i] > bbox.H[i])
+			dmin += (center[i] - bbox.H[i])*(center[i] - bbox.H[i]);
+	}
+	if (dmin <= sqr_radius)
+		return true;
+	return false;
+};
+
+BBox Sphere::get_bbox() const
+{
+	return BBox(center - radius, center + radius);
+}
+
+bool Box::intersect(const Ray &ray, Float &dist) const
+{
+	register Float tnear = -Inf;
+	register Float tfar = Inf;
+	register Float t1, t2;
+
+	for (int i = 0; i < 3; i++)
+	{
+		if (ray.dir[i] == 0) {
+			/* ray is parallel to these planes */
+			if (ray.o[i] < L[i] || ray.o[i] > H[i])
+				return false;
+		}
+		else
+		{
+			/* compute the intersection distance of the planes */
+			t1 = (L[i] - ray.o[i]) / ray.dir[i];
+			t2 = (H[i] - ray.o[i]) / ray.dir[i];
+
+			if (t1 > t2)
+				swap(t1, t2);
+
+			if (t1 > tnear)
+				tnear = t1; /* want largest Tnear */
+			if (t2 < tfar)
+				tfar = t2; /* want smallest Tfar */
+			if (tnear > tfar || tfar < 0)
+				return false; /* box missed; box is behind ray */
+		}
+	}
+
+	if (tnear < dist)
+	{
+		dist = tnear;
+		return true;
+	}
+	return false;
+}
+
+bool Box::intersect_bbox(const BBox &bbox) const
+{
+	return (
+	H.x > bbox.L.x && L.x < bbox.H.x &&
+	H.y > bbox.L.y && L.y < bbox.H.y &&
+	H.z > bbox.L.z && L.z < bbox.H.z);
+}
+
+const Vector3 Box::normal(const Vector3 &P) const
+{
+	register Vector3 l = P - L;
+	register Vector3 h = H - P;
+
+	if (l.x < h.x)
+		h.x = -1;
+	else
+	{
+		l.x = h.x;
+		h.x = +1;
+	}
+
+	if (l.y < h.y)
+		h.y = -1;
+	else
+	{
+		l.y = h.y;
+		h.y = +1;
+	}
+
+	if (l.z < h.z)
+		h.z = -1;
+	else
+	{
+		l.z = h.z;
+		h.z = +1;
+	}
+
+	if (l.x > l.y)
+	{
+		h.x = 0;
+		if (l.y > l.z)
+			h.y = 0;
+		else
+			h.z = 0;
+	}
+	else
+	{
+		h.y = 0;
+		if (l.x > l.z)
+			h.x = 0;
+		else
+			h.z = 0;
+	}
+	return h;
+}
+
+#ifdef TRI_PLUCKER
+inline void Plucker(const Vector3 &p, const Vector3 &q, Float* pl)
+{
+    pl[0] = p.x*q.y - q.x*p.y;
+    pl[1] = p.x*q.z - q.x*p.z;
+    pl[2] = p.x - q.x;
+    pl[3] = p.y*q.z - q.y*p.z;
+    pl[4] = p.z - q.z;
+    pl[5] = q.y - p.y;
+}
+
+inline Float Side(const Float* pla, const Float* plb)
+{
+    return pla[0]*plb[4] + pla[1]*plb[5] + pla[2]*plb[3] + pla[4]*plb[0] + pla[5]*plb[1] + pla[3]*plb[2];
+}
+#endif
+
+Triangle::Triangle(Vertex *aA, Vertex *aB, Vertex *aC, Material *amaterial)
+	: A(aA), B(aB), C(aC)
+{
+	material = amaterial;
+
+	const Vector3 c = B->P - A->P;
+	const Vector3 b = C->P - A->P;
+
+	N = cross(c, b);
+	N.normalize();
+
+#ifdef TRI_PLUCKER
+	Plucker(B->P,C->P,pla);
+	Plucker(C->P,A->P,plb);
+	Plucker(A->P,B->P,plc);
+#endif
+
+#if defined(TRI_BARI) || defined(TRI_BARI_PRE)
+	if (fabsf(N.x) > fabsf(N.y))
+	{
+		if (fabsf(N.x) > fabsf(N.z))
+			k = 0;
+		else
+			k = 2;
+	}
+	else
+	{
+		if (fabsf(N.y) > fabsf(N.z))
+			k = 1;
+		else
+			k = 2;
+	}
+#endif
+#ifdef TRI_BARI_PRE
+	int u = (k + 1) % 3;
+	int v = (k + 2) % 3;
+
+	Float krec = 1.0 / N[k];
+	nu = N[u] * krec;
+	nv = N[v] * krec;
+	nd = dot(N, A->P) * krec;
+
+	// first line equation
+	Float reci = 1.0f / (b[u] * c[v] - b[v] * c[u]);
+	bnu = b[u] * reci;
+	bnv = -b[v] * reci;
+
+	// second line equation
+	cnu = -c[u] * reci;
+	cnv = c[v] * reci;
+#endif
+}
+
+bool Triangle::intersect(const Ray &ray, Float &dist) const
+{
+#ifdef TRI_PLUCKER
+	Float plr[6];
+	Plucker(ray.o, ray.o+ray.dir, plr);
+	const bool side0 = Side(plr, pla) >= 0.0;
+	const bool side1 = Side(plr, plb) >= 0.0;
+	if (side0 != side1)
+		return false;
+	const bool side2 = Side(plr, plc) >= 0.0;
+	if (side0 != side2)
+		return false;
+	const Float t = - dot( (ray.o - A->P), N) / dot(ray.dir,N);
+	if(t <= Eps || t >= dist)
+		return false;
+	dist = t;
+	return true;
+#endif
+
+#if defined(TRI_BARI) || defined(TRI_BARI_PRE)
+	static const int modulo3[5] = {0,1,2,0,1};
+	const Vector3 &O = ray.o;
+	const Vector3 &D = ray.dir;
+	register const int u = modulo3[k+1];
+	register const int v = modulo3[k+2];
+#endif
+#ifdef TRI_BARI_PRE
+	const Float t = (nd - O[k] - nu * O[u] - nv * O[v]) / (D[k] + nu * D[u] + nv * D[v]);
+
+	if (t >= dist || t < Eps)
+		return false;
+
+	const Float hu = O[u] + t * D[u] - A->P[u];
+	const Float hv = O[v] + t * D[v] - A->P[v];
+	const Float beta = hv * bnu + hu * bnv;
+
+	if (beta < 0.)
+		return false;
+
+	const Float gamma = hu * cnv + hv * cnu;
+	if (gamma < 0. || beta + gamma > 1.)
+		return false;
+
+	dist = t;
+	return true;
+#endif
+
+#ifdef TRI_BARI
+	// original barycentric coordinates based intesection
+	// not optimized, just for reference
+	const Vector3 c = B - A;
+	const Vector3 b = C - A;
+	// distance test
+	const Float t = - dot( (O-A), N) / dot(D,N);
+	if (t < Eps || t > dist)
+		return false;
+
+	// calc hitpoint
+	const Float Hu = O[u] + t * D[u] - A[u];
+	const Float Hv = O[v] + t * D[v] - A[v];
+	const Float beta = (b[u] * Hv - b[v] * Hu) / (b[u] * c[v] - b[v] * c[u]);
+	if (beta < 0)
+		return false;
+	const Float gamma = (c[v] * Hu - c[u] * Hv) / (b[u] * c[v] - b[v] * c[u]);
+	if (gamma < 0)
+		return false;
+	if (beta+gamma > 1)
+		return false;
+	dist = t;
+	return true;
+#endif
+}
+
+bool Triangle::intersect_bbox(const BBox &bbox) const
+{
+	const Vector3 boxcenter = (bbox.L+bbox.H)*0.5;
+	const Vector3 boxhalfsize = (bbox.H-bbox.L)*0.5;
+	const Vector3 v0 = A->P - boxcenter;
+	const Vector3 v1 = B->P - boxcenter;
+	const Vector3 v2 = C->P - boxcenter;
+	const Vector3 e0 = v1-v0;
+	const Vector3 e1 = v2-v1;
+	const Vector3 e2 = v0-v2;
+
+	Float fex = fabsf(e0.x);
+	Float fey = fabsf(e0.y);
+	Float fez = fabsf(e0.z);
+
+	Float p0,p1,p2,min,max,rad;
+
+	p0 = e0.z*v0.y - e0.y*v0.z;
+	p2 = e0.z*v2.y - e0.y*v2.z;
+	if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
+	rad = fez * boxhalfsize.y + fey * boxhalfsize.z;
+	if(min>rad || max<-rad) return false;
+
+	p0 = -e0.z*v0.x + e0.x*v0.z;
+	p2 = -e0.z*v2.x + e0.x*v2.z;
+	if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
+	rad = fez * boxhalfsize.x + fex * boxhalfsize.z;
+	if(min>rad || max<-rad) return false;
+
+	p1 = e0.y*v1.x - e0.x*v1.y;
+	p2 = e0.y*v2.x - e0.x*v2.y;
+	if(p2<p1) {min=p2; max=p1;} else {min=p1; max=p2;}
+	rad = fey * boxhalfsize.x + fex * boxhalfsize.y;
+	if(min>rad || max<-rad) return false;
+
+	fex = fabsf(e1.x);
+	fey = fabsf(e1.y);
+	fez = fabsf(e1.z);
+
+	p0 = e1.z*v0.y - e1.y*v0.z;
+	p2 = e1.z*v2.y - e1.y*v2.z;
+	if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
+	rad = fez * boxhalfsize.y + fey * boxhalfsize.z;
+	if(min>rad || max<-rad) return false;
+
+	p0 = -e1.z*v0.x + e1.x*v0.z;
+	p2 = -e1.z*v2.x + e1.x*v2.z;
+	if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
+	rad = fez * boxhalfsize.x + fex * boxhalfsize.z;
+	if(min>rad || max<-rad) return false;
+
+	p0 = e1.y*v0.x - e1.x*v0.y;
+	p1 = e1.y*v1.x - e1.x*v1.y;
+	if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;}
+	rad = fey * boxhalfsize.x + fex * boxhalfsize.y;
+	if(min>rad || max<-rad) return false;
+
+	fex = fabsf(e2.x);
+	fey = fabsf(e2.y);
+	fez = fabsf(e2.z);
+
+	p0 = e2.z*v0.y - e2.y*v0.z;
+	p1 = e2.z*v1.y - e2.y*v1.z;
+	if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;}
+	rad = fez * boxhalfsize.y + fey * boxhalfsize.z;
+	if(min>rad || max<-rad) return false;
+
+	p0 = -e2.z*v0.x + e2.x*v0.z;
+	p1 = -e2.z*v1.x + e2.x*v1.z;
+	if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;}
+	rad = fez * boxhalfsize.x + fex * boxhalfsize.z;
+	if(min>rad || max<-rad) return false;
+
+	p1 = e2.y*v1.x - e2.x*v1.y;
+	p2 = e2.y*v2.x - e2.x*v2.y;
+	if(p2<p1) {min=p2; max=p1;} else {min=p1; max=p2;}
+	rad = fey * boxhalfsize.x + fex * boxhalfsize.y;
+	if(min>rad || max<-rad) return false;
+
+	/* test overlap in the {x,y,z}-directions */
+	/* test in X-direction */
+	min = v0.x;
+	if (v1.x < min) min = v1.x;
+	if (v2.x < min) min = v2.x;
+	max = v0.x;
+	if (v1.x > max) max = v1.x;
+	if (v2.x > max) max = v2.x;
+	if(min>boxhalfsize.x || max<-boxhalfsize.x) return false;
+
+	/* test in Y-direction */
+	min = v0.y;
+	if (v1.y < min) min = v1.y;
+	if (v2.y < min) min = v2.y;
+	max = v0.y;
+	if (v1.y > max) max = v1.y;
+	if (v2.y > max) max = v2.y;
+	if(min>boxhalfsize.y || max<-boxhalfsize.y) return false;
+
+	/* test in Z-direction */
+	min = v0.z;
+	if (v1.z < min) min = v1.z;
+	if (v2.z < min) min = v2.z;
+	max = v0.z;
+	if (v1.z > max) max = v1.z;
+	if (v2.z > max) max = v2.z;
+	if(min>boxhalfsize.z || max<-boxhalfsize.z) return false;
+
+	/*  test if the box intersects the plane of the triangle */
+	Vector3 vmin,vmax;
+	Float v;
+	for(int q=0;q<3;q++)
+	{
+		v=v0[q];
+		if(N[q]>0.0f)
+		{
+			vmin.cell[q]=-boxhalfsize[q] - v;
+			vmax.cell[q]= boxhalfsize[q] - v;
+		}
+		else
+		{
+			vmin.cell[q]= boxhalfsize[q] - v;
+			vmax.cell[q]=-boxhalfsize[q] - v;
+		}
+	}
+	if(dot(N,vmin)>0.0f) return false;
+	if(dot(N,vmax)>=0.0f) return true;
+
+	return false;
+}
+
+BBox Triangle::get_bbox() const
+{
+	BBox bbox = BBox();
+	bbox.L = A->P;
+	if (B->P.x < bbox.L.x)  bbox.L.x = B->P.x;
+	if (C->P.x < bbox.L.x)  bbox.L.x = C->P.x;
+	if (B->P.y < bbox.L.y)  bbox.L.y = B->P.y;
+	if (C->P.y < bbox.L.y)  bbox.L.y = C->P.y;
+	if (B->P.z < bbox.L.z)  bbox.L.z = B->P.z;
+	if (C->P.z < bbox.L.z)  bbox.L.z = C->P.z;
+	bbox.H = A->P;
+	if (B->P.x > bbox.H.x)  bbox.H.x = B->P.x;
+	if (C->P.x > bbox.H.x)  bbox.H.x = C->P.x;
+	if (B->P.y > bbox.H.y)  bbox.H.y = B->P.y;
+	if (C->P.y > bbox.H.y)  bbox.H.y = C->P.y;
+	if (B->P.z > bbox.H.z)  bbox.H.z = B->P.z;
+	if (C->P.z > bbox.H.z)  bbox.H.z = C->P.z;
+	return bbox;
+}
+
+ostream & Sphere::dump(ostream &st) const
+{
+	return st << "(sphere," << center << "," << radius << ")";
+}
+
+ostream & Box::dump(ostream &st) const
+{
+	return st << "(box," << L << "," << H << ")";
+}
+
+ostream & Vertex::dump(ostream &st) const
+{
+	return st << "(v," << P << ")";
+}
+
+ostream & NormalVertex::dump(ostream &st) const
+{
+	return st << "(vn," << P << "," << N << ")";
+}
+
+ostream & Triangle::dump(ostream &st) const
+{
+	int idxA, idxB, idxC;
+	if (!vertex_index.get(A, idxA))
+		st << *A << "," << endl;
+	if (!vertex_index.get(B, idxB))
+		st << *B << "," << endl;
+	if (!vertex_index.get(C, idxC))
+		st << *C << "," << endl;
+	return st << "(t," << idxA << "," << idxB << "," << idxC << ")";
+}