--- a/src/scene.cc	Wed Dec 05 18:54:23 2007 +0100
+++ b/src/scene.cc	Fri Dec 07 14:56:39 2007 +0100
@@ -94,7 +94,7 @@
 	return true;
 }
 
-bool Sphere::intersect(const Ray &ray, Float &dist)
+bool Sphere::intersect(const Ray &ray, Float &dist) const
 {
 	Vector3 V = ray.o - center;
 	register Float d = -dot(V, ray.dir);
@@ -110,7 +110,7 @@
 	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) const
 {
 	//allts = new vector<Float>();
 
@@ -144,7 +144,7 @@
 	return false;
 }
 
-BBox Sphere::get_bbox()
+BBox Sphere::get_bbox() const
 {
 	BBox bbox = BBox();
 	bbox.L = center - radius;
@@ -152,7 +152,7 @@
 	return bbox;
 }
 
-bool Box::intersect(const Ray &ray, Float &dist)
+bool Box::intersect(const Ray &ray, Float &dist) const
 {
 	Float a,b;
 	bool res = get_bbox().intersect(ray, a, b);
@@ -165,7 +165,7 @@
 		return false;
 }
 
-Vector3 Box::normal(Vector3 &P)
+Vector3 Box::normal(Vector3 &P) const
 {
 	Vector3 N;
 	for (int i = 0; i < 3; i++)
@@ -186,18 +186,42 @@
 	return N;
 }
 
-// this initialization and following intersection methods implements
-// Fast Triangle Intersection algorithm from
-// http://www.mpi-inf.mpg.de/~wald/PhD/
+#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(const Vector3 &aA, const Vector3 &aB, const Vector3 &aC, Material *amaterial)
 	: A(aA), B(aB), C(aC)
 {
 	material = amaterial;
-	Vector3 c = B - A;
-	Vector3 b = C - A;
+	material->reflection = 0;
+
+	const Vector3 c = B - A;
+	const Vector3 b = C - A;
+
+	N = -cross(c, b);
+	N.normalize();
 
-	N = cross(c, b);
+#ifdef TRI_PLUCKER
+	Plucker(B,C,pla);
+	Plucker(C,A,plb);
+	Plucker(A,B,plc);
+#endif
 
+#if defined(TRI_BARI) || defined(TRI_BARI_PRE)
 	if (fabsf(N.x) > fabsf(N.y))
 	{
 		if (fabsf(N.x) > fabsf(N.z))
@@ -212,11 +236,12 @@
 		else
 			k = 2;
 	}
-
+#endif
+#ifdef TRI_BARI_PRE
 	int u = (k + 1) % 3;
 	int v = (k + 2) % 3;
 
-	Float krec = 1.0f / N[k];
+	Float krec = 1.0 / N[k];
 	nu = N[u] * krec;
 	nv = N[v] * krec;
 	nd = dot(N, A) * krec;
@@ -227,47 +252,85 @@
 	bnv = -b[v] * reci;
 
 	// second line equation
-	cnu = c[v] * reci;
-	cnv = -c[u] * reci;
-
-	// finalize normal
-	N.normalize();
+	cnu = -c[u] * reci;
+	cnv = c[v] * reci;
+#endif
 }
 
-// see comment for previous method
-bool Triangle::intersect(const Ray &ray, Float &dist)
+bool Triangle::intersect(const Ray &ray, Float &dist) const
 {
-	Vector3 O = ray.o;
-	Vector3 D = ray.dir;
+#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), N) / dot(ray.dir,N);
+	if(t <= Eps || t >= dist)
+		return false;
+	dist = t;
+	return true;
+#endif
 
-	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;
+#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 < 0 || t >= dist)
+	if (t >= dist || t < Eps)
 		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;
+	const Float hu = O[u] + t * D[u] - A[u];
+	const Float hv = O[v] + t * D[v] - A[v];
+	const Float beta = hv * bnu + hu * bnv;
 
-	if (beta < 0)
+	if (beta < 0.)
 		return false;
 
-	Float gamma = hu * cnu + hv * cnv;
-	if (gamma < 0)
-		return false;
-
-	if ((beta + gamma) > 1)
+	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
 }
 
-BBox Triangle::get_bbox()
+BBox Triangle::get_bbox() const
 {
 	BBox bbox = BBox();
 	bbox.L = A;