# HG changeset patch
# User Radek Brich <radek.brich@devl.cz>
# Date 1197035799 -3600
# Node ID b8232edee786092bb1be0c342d0b7b523583b192
# Parent  d0d76e8a5203c2147a55ad7d19172c30c77f2d3e
tuned ray-triangle intersection, now there are three algorithms to choose from:
  Plucker, Barycentric and Barycentric with preprocessing (Wald)
methods in Vector and Shape (and derivates) made const

diff -r d0d76e8a5203 -r b8232edee786 ccdemos/Makefile
--- a/ccdemos/Makefile	Wed Dec 05 18:54:23 2007 +0100
+++ b/ccdemos/Makefile	Fri Dec 07 14:56:39 2007 +0100
@@ -36,4 +36,4 @@
 	$(CXX) -c -o $@ $*.cc $(CCFLAGS) $(DEFS)
 
 clean:
-	rm -f spheres_shadow realtime *.o
+	rm -f spheres_shadow realtime realtime_dragon *.o
diff -r d0d76e8a5203 -r b8232edee786 demos/dragon.py
--- a/demos/dragon.py	Wed Dec 05 18:54:23 2007 +0100
+++ b/demos/dragon.py	Fri Dec 07 14:56:39 2007 +0100
@@ -46,8 +46,11 @@
 mat = Material(colour=(0.9, 0.9, 0.9))
 LoadStanfordPlyFile(rt, mat, "../models/dragon/dragon_vrip_res4.ply", 29.0)
 
-light = Light(position=(-5.0, 2.0, 8.0), colour=(0.9, 0.3, 0.6))
-rt.addlight(light)
+light1 = Light(position=(-5.0, 2.0, 8.0), colour=(0.9, 0.3, 0.2))
+rt.addlight(light1)
+
+light2 = Light(position=(3.0, 0.0, 9.0), colour=(0.0, 1.0, 0.2))
+rt.addlight(light2)
 
 imagesize = (800, 600)
 data = rt.render(imagesize)
diff -r d0d76e8a5203 -r b8232edee786 include/scene.h
--- a/include/scene.h	Wed Dec 05 18:54:23 2007 +0100
+++ b/include/scene.h	Fri Dec 07 14:56:39 2007 +0100
@@ -11,8 +11,18 @@
 #include <vector>
 
 #include "noise.h"
+#include "vector.h"
 
-#include "vector.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;
 
@@ -140,15 +150,15 @@
 	virtual ~Shape() {};
 
 	// first intersection point
-	virtual bool intersect(const Ray &ray, Float &dist) = 0;
+	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) = 0;
+	virtual bool intersect_all(const Ray &ray, Float dist, vector<Float> &allts) const = 0;
 
 	// normal at point P
-	virtual Vector3 normal(Vector3 &P) = 0;
+	virtual Vector3 normal(Vector3 &P) const = 0;
 
-	virtual BBox get_bbox() = 0;
+	virtual BBox get_bbox() const = 0;
 };
 
 class ShapeList: public vector<Shape*>
@@ -166,10 +176,10 @@
 	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();
+	bool intersect(const Ray &ray, Float &dist) const;
+	bool intersect_all(const Ray &ray, Float dist, vector<Float> &allts) const;
+	Vector3 normal(Vector3 &P) const { return (P - center) * inv_radius; };
+	BBox get_bbox() const;
 };
 
 class Box: public Shape
@@ -184,26 +194,34 @@
 				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); };
+	bool intersect(const Ray &ray, Float &dist) const;
+	bool intersect_all(const Ray &ray, Float dist, vector<Float> &allts) const { return false; };
+	Vector3 normal(Vector3 &P) const;
+	BBox get_bbox() const { return BBox(L, H); };
 };
 
 class Triangle: public Shape
 {
+#ifdef TRI_BARI_PRE
+	Float nu, nv, nd;
 	int k; // dominant axis
-	Float nu, nv, nd;
 	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
 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();
+	bool intersect(const Ray &ray, Float &dist) const;
+	bool intersect_all(const Ray &ray, Float dist, vector<Float> &allts) const {return false;};
+	Vector3 normal(Vector3 &) const { return N; };
+	BBox get_bbox() const;
 };
 
 #endif
diff -r d0d76e8a5203 -r b8232edee786 include/vector.h
--- a/include/vector.h	Wed Dec 05 18:54:23 2007 +0100
+++ b/include/vector.h	Fri Dec 07 14:56:39 2007 +0100
@@ -32,9 +32,9 @@
 	Vector3(Float ax, Float ay, Float az): x(ax), y(ay), z(az) {};
 
 	// index operator
-	Float &operator[](int index)	{ return cell[index]; };
+	const Float &operator[](int index) const { return cell[index]; };
 
-	bool operator==(Vector3 &v) { return x==v.x && y==v.y && z==v.z; };
+	bool operator==(Vector3 &v) const { return x==v.x && y==v.y && z==v.z; };
 
 	// normalize
 	Vector3 normalize()
@@ -47,18 +47,18 @@
 	}
 
 	// get normalized copy
-	Vector3 unit()
+	Vector3 unit() const
 	{
 		Vector3 u(*this);
 		return u.normalize();;
 	}
 
 	// square magnitude, magnitude
-	Float mag2()	{ return x * x + y * y + z * z; }
-	Float mag()	{ return sqrtf(mag2()); }
+	Float mag2() const	{ return x * x + y * y + z * z; }
+	Float mag() const	{ return sqrtf(mag2()); }
 
 	// negative
-	Vector3 operator-()	{ return Vector3(-x, -y, -z); }
+	Vector3 operator-() const { return Vector3(-x, -y, -z); }
 
 	// accumulate
 	Vector3 operator+=(const Vector3 &v)
diff -r d0d76e8a5203 -r b8232edee786 src/kdtree.cc
--- a/src/kdtree.cc	Wed Dec 05 18:54:23 2007 +0100
+++ b/src/kdtree.cc	Fri Dec 07 14:56:39 2007 +0100
@@ -357,9 +357,9 @@
 
 			/* setup the new exit point and push it onto stack */
 			exPt = new StackElem(farchild, t, Vector3());
-			exPt->pb[axis] = splitVal;
-			exPt->pb[(axis+1)%3] = ray.o.cell[(axis+1)%3] + t * ray.dir.cell[(axis+1)%3];
-			exPt->pb[(axis+2)%3] = ray.o.cell[(axis+2)%3] + t * ray.dir.cell[(axis+2)%3];
+			exPt->pb.cell[axis] = splitVal;
+			exPt->pb.cell[(axis+1)%3] = ray.o.cell[(axis+1)%3] + t * ray.dir.cell[(axis+1)%3];
+			exPt->pb.cell[(axis+2)%3] = ray.o.cell[(axis+2)%3] + t * ray.dir.cell[(axis+2)%3];
 			st.push_back(exPt);
 		} /* while */
 
diff -r d0d76e8a5203 -r b8232edee786 src/raytracer.cc
--- a/src/raytracer.cc	Wed Dec 05 18:54:23 2007 +0100
+++ b/src/raytracer.cc	Fri Dec 07 14:56:39 2007 +0100
@@ -100,6 +100,9 @@
 		Colour acc = Colour();
 		Vector3 P = ray.o + ray.dir * nearest_distance; // point of intersection
 		Vector3 normal = nearest_shape->normal(P);
+		// make shapes double sided
+		if (dot(normal, ray.dir) > 0.0)
+			normal = - normal;
 		acc = PhongShader_ambient(*nearest_shape->material, P);
 
 		vector<Light*>::iterator light;
diff -r d0d76e8a5203 -r b8232edee786 src/scene.cc
--- 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;