--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/scene.cc	Thu Oct 25 16:40:22 2007 +0200
@@ -0,0 +1,167 @@
+/*
+ * C++ RayTracer
+ * file: scene.cc
+ *
+ * Radek Brich, 2006
+ */
+
+#include <math.h>
+#include "scene.h"
+
+bool Sphere::intersect(const Ray &ray, float &dist)
+{
+	Vector3 V = ((Ray)ray).a - 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;
+		if (t1 > 0)
+		{
+			if (t1 < dist) {
+				dist = t1;
+				return true;
+			}
+		} else {
+			float t2 = Vd + Det;
+			if (t2 > 0)
+			{
+				// ray from inside of the sphere
+				dist = t2;
+				return true;
+			}
+		}
+	}
+	return false;
+}
+
+bool Sphere::intersect_all(const Ray &ray, float dist, vector<float> &allts)
+{
+	//allts = new vector<float>();
+
+	Vector3 V = ((Ray)ray).a - 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 Plane::intersect(const Ray &ray, float &dist)
+{
+	float dir = dot(N, ray.dir);
+	if (dir != 0)
+	{
+		float newdist = -(dot(N, ray.a) + d) / dir;
+		if (newdist > 0 && newdist < dist) {
+			dist = newdist;
+			return true;
+		}
+	}
+	return false;
+}
+
+// this initialization and following intersection methods implements
+// Fast Triangle Intersection algorithm from
+// http://www.mpi-inf.mpg.de/~wald/PhD/
+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;
+
+	N = cross(c, b);
+
+	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;
+	}
+
+	int u = (k + 1) % 3;
+	int v = (k + 2) % 3;
+
+	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]);
+	bnu = b[u] * reci;
+	bnv = -b[v] * reci;
+
+	// second line equation
+	cnu = c[v] * reci;
+	cnv = -c[u] * reci;
+
+	// finalize normal
+	N.normalize();
+}
+
+// see comment for previous method
+bool Triangle::intersect(const Ray &ray, float &dist)
+{
+	Vector3 O = ray.a;
+	Vector3 D = ray.dir;
+
+	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 (!(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;
+
+	if (beta < 0)
+		return false;
+
+	float gamma = hu * cnu + hv * cnv;
+	if (gamma < 0)
+		return false;
+
+	if ((beta + gamma) > 1)
+		return false;
+
+	dist = t;
+	return true;
+}