lots of bugfixes

include/mitsuba/render/gkdtree.h

@@ -512,7 +512,7 @@ public:
 		m_stopPrims = 4;
 		m_maxBadRefines = 3;
 		m_exactPrimThreshold = 65536;
-		m_maxDepth = 1;
+		m_maxDepth = 0;
 		m_retract = true;
 		m_parallelBuild = true;
 	}
@@ -899,6 +899,13 @@ public:
 		expPrimitivesIntersected /= rootSA;
 		sahCost /= rootSA;
 	
+		/* Slightly enlarge the bounding box 
+		   (necessary e.g. when the scene is planar) */
+		m_aabb.min -= (m_aabb.max-m_aabb.min) * Epsilon
+			+ Vector(Epsilon, Epsilon, Epsilon);
+		m_aabb.max += (m_aabb.max-m_aabb.min) * Epsilon
+			+ Vector(Epsilon, Epsilon, Epsilon);
+
 		Log(EDebug, "Detailed kd-tree statistics:");
 		Log(EDebug, "   Inner nodes               : %i", ctx.innerNodeCount);
 		Log(EDebug, "   Leaf nodes                : %i", ctx.leafNodeCount);
@@ -923,42 +930,12 @@ public:
 	/**
 	 * \brief Intersect a ray against all primitives stored in the kd-tree
 	 */
-	bool rayIntersect(const Ray &ray, Intersection &its) {
+	bool rayIntersect(const Ray &ray, Intersection &its);
-		uint32_t temp[MTS_KD_INTERSECTION_TEMP];
-		its.t = std::numeric_limits<Float>::infinity(); 
-		Float mint, maxt;
-
-		if (m_aabb.rayIntersect(ray, mint, maxt)) {
-			if (ray.mint > mint) mint = ray.mint;
-			if (ray.maxt < maxt) maxt = ray.maxt;
-
-			if (EXPECT_TAKEN(maxt > mint)) {
-				if (rayIntersect<false>(ray, mint, maxt, its.t, temp)) {
-					cast()->fillIntersectionDetails(ray, its.t, temp, its);
-					return true;
-				}
-			}
-		}
-		return false;
-	}
 
 	/**
 	 * \brief Test a ray for intersection against all primitives stored in the kd-tree
 	 */
-	bool rayIntersect(const Ray &ray) {
+	bool rayIntersect(const Ray &ray);
-		uint32_t temp[MTS_KD_INTERSECTION_TEMP];
-		Float mint, maxt, t;
-
-		if (m_aabb.rayIntersect(ray, mint, maxt)) {
-			if (ray.mint > mint) mint = ray.mint;
-			if (ray.maxt < maxt) maxt = ray.maxt;
-
-			if (EXPECT_TAKEN(maxt > mint)) 
-				if (rayIntersect<true>(ray, mint, maxt, t, temp)) 
-					return true;
-		}
-		return false;
-	}
 
 protected:
 	/// Primitive classification during tree-construction
@@ -1302,7 +1279,7 @@ protected:
 	 * \brief Hashed mailbox implementation
 	 */
 	struct HashedMailbox {
-		HashedMailbox() {
+		inline HashedMailbox() {
 			memset(entries, 0xFF, sizeof(index_type)*MTS_KD_MAILBOX_SIZE);
 		}
 
@@ -2455,7 +2432,7 @@ protected:
 	};
 
 	/// Ray traversal stack entry for incoherent ray tracing
-	struct KDStackEntry {
+	struct KDStackEntryHavran {
 		/* Pointer to the far child */
 		const KDNode * __restrict node;
 		/* Distance traveled along the ray (entry or exit) */
@@ -2467,16 +2444,17 @@ protected:
 	};
 
 	/**
-	 * \brief Internal kd-tree traversal implementation
+	 * \brief Internal kd-tree traversal implementation (Havran variant)
 	 */
-	template<bool shadowRay> FINLINE bool rayIntersect(const Ray &ray, 
+	template<bool shadowRay> FINLINE bool rayIntersectHavran(const Ray &ray, 
 			Float mint, Float maxt, Float &t, void *temp) {
-		static const int prevAxisTable[] = { 2, 0, 1 };
+		KDStackEntryHavran stack[MTS_KD_MAXDEPTH];
-		static const int nextAxisTable[] = { 1, 2, 0 };
-
-		KDStackEntry stack[MTS_KD_MAXDEPTH*2];
 		const KDNode * __restrict farChild,
 					 * __restrict currNode = m_nodes;
+		#if 0
+		static const int prevAxisTable[] = { 2, 0, 1 };
+		static const int nextAxisTable[] = { 1, 2, 0 };
+		#endif
 
 		#if defined(MTS_KD_MAILBOX_ENABLED)
 		HashedMailbox mailbox;
@@ -2540,14 +2518,18 @@ protected:
 				stack[exPt].prev = tmp;
 				stack[exPt].t = t;
 				stack[exPt].node = farChild;
+
 				#if 1
+				/* Faster than the original code with the 
+				   prevAxis & nextAxis table */
+				stack[exPt].p = ray(t);
+				#else
 				const int nextAxis = nextAxisTable[axis];
 				const int prevAxis = prevAxisTable[axis];
 				stack[exPt].p[nextAxis] = ray.o[nextAxis] + t*ray.d[nextAxis];
 				stack[exPt].p[prevAxis] = ray.o[prevAxis] + t*ray.d[prevAxis];
-#else
-				stack[exPt].p = ray(t);
 				#endif
+
 				stack[exPt].p[axis] = splitVal;
 			}
 
@@ -2560,6 +2542,8 @@ protected:
 
 			/* Floating-point arithmetic.. - use both absolute and relative 
 			   epsilons when looking for intersections in the subinterval */
+//			const Float searchStart = std::max(mint, stack[enPt].t * m_eps - eps);
+//			Float searchEnd = std::min(maxt, stack[exPt].t * p_eps + eps);
 			const Float searchStart = std::max(mint, stack[enPt].t * m_eps - eps);
 			Float searchEnd = std::min(maxt, stack[exPt].t * p_eps + eps);
 
@@ -2570,7 +2554,7 @@ protected:
 				const index_type primIdx = m_indices[entry];
 
 				#if defined(MTS_KD_MAILBOX_ENABLED)
-				if (mailbox.get(primIdx))
+				if (mailbox.contains(primIdx)) 
 					continue;
 				#endif
 
@@ -2590,10 +2574,8 @@ protected:
 				#endif
 			}
 
-			if (foundIntersection) {
+			if (foundIntersection) 
-				t = searchEnd;
 				return true;
-			}
 
 			/* Pop from the stack and advance to the next node on the interval */
 			enPt = exPt;
@@ -2604,6 +2586,33 @@ protected:
 		return false;
 	}
 
+	struct KDStackEntry {
+		const KDNode * __restrict node;
+		Float mint, maxt;
+	};
+
+	/**
+	 * \brief Internal kd-tree traversal implementation (Havran variant)
+	 */
+	template<bool shadowRay> FINLINE bool rayIntersect(const Ray &ray, 
+			Float mint, Float maxt, Float &t, void *temp) {
+		KDStackEntry stack[MTS_KD_MAXDEPTH];
+		int stackPos = 0;
+		KDNode *node = m_nodes;
+
+		while (node != NULL) {
+			if (EXPECT_TAKEN(!node->isLeaf())) {
+				float split = (Float) node->getSplit();
+				int axis = node->getAxis();
+				float t = (split - ray.o[axis]) * ray.dRcp[axis];
+
+				const KDNode * __restrict first, * __restrict second;
+
+			}
+		}
+		return false;
+	}
+
 private:
 	KDNode *m_nodes;
 	index_type *m_indices;
@@ -2622,6 +2631,44 @@ private:
 	BuildInterface m_interface;
 };
 
+
+template <typename Derived> bool GenericKDTree<Derived>::rayIntersect(
+		const Ray &ray, Intersection &its) {
+
+	uint32_t temp[MTS_KD_INTERSECTION_TEMP];
+	its.t = std::numeric_limits<Float>::infinity(); 
+	Float mint, maxt;
+
+	if (m_aabb.rayIntersect(ray, mint, maxt)) {
+		if (ray.mint > mint) mint = ray.mint;
+		if (ray.maxt < maxt) maxt = ray.maxt;
+
+		if (EXPECT_TAKEN(maxt > mint)) {
+			if (rayIntersectHavran<false>(ray, mint, maxt, its.t, temp)) {
+				cast()->fillIntersectionDetails(ray, its.t, temp, its);
+				return true;
+			}
+		}
+	}
+	return false;
+}
+
+template <typename Derived> bool GenericKDTree<Derived>::rayIntersect(
+		const Ray &ray) {
+	uint32_t temp[MTS_KD_INTERSECTION_TEMP];
+	Float mint, maxt, t;
+
+	if (m_aabb.rayIntersect(ray, mint, maxt)) {
+		if (ray.mint > mint) mint = ray.mint;
+		if (ray.maxt < maxt) maxt = ray.maxt;
+
+		if (EXPECT_TAKEN(maxt > mint)) 
+			if (rayIntersectHavran<true>(ray, mint, maxt, t, temp)) 
+				return true;
+	}
+	return false;
+}
+
 MTS_NAMESPACE_END
 
 #endif /* __KDTREE_GENERIC_H */

include/mitsuba/render/kdtree.h

@@ -30,10 +30,10 @@
  * whole following block:
  */
 #ifdef MTS_SSE
-#define MTS_USE_TRIACCEL4 1
+//#define MTS_USE_TRIACCEL4 1
-#define MTS_USE_TRIACCEL 1
+//#define MTS_USE_TRIACCEL 1
 #else
-#define MTS_USE_TRIACCEL 1
+//#define MTS_USE_TRIACCEL 1
 #endif
 
 #if defined(MTS_HAS_COHERENT_RT)

src/libcore/triangle.cpp

@@ -35,7 +35,7 @@ bool Triangle::rayIntersect(const Vertex *buffer, const Ray &ray,
 	/* if determinant is near zero, ray lies in plane of triangle */
 	Float det = dot(edge1, pvec);
 
-	if (det > -Epsilon && det < Epsilon)
+	if (det > -1e-8f && det < 1e-8f)
 		return false;
 	Float inv_det = 1.0f / det;
 

src/librender/kdtree_traversal.cpp

@@ -186,7 +186,7 @@ bool KDTree::rayIntersect(const Ray &ray, Intersection &its, Float mint, Float m
 			   Moeller and Trumbore (without having done any pre-computation) */
 			const KDTriangle &kdTri = m_triangles[m_indices[entry]];
 			if (EXPECT_TAKEN(kdTri.index != KNoTriangleFlag)) {
-				const TriMesh *mesh = m_meshes[kdTri.shapeIndex];
+				const TriMesh *mesh = (const TriMesh *) m_shapes[kdTri.shapeIndex];
 				const Triangle &tri = mesh->getTriangles()[kdTri.index];
 				if (tri.rayIntersect(mesh->getVertexBuffer(), ray, u, v, t)) {
 					if (t >= searchStart && t < searchEnd) {
@@ -202,7 +202,7 @@ bool KDTree::rayIntersect(const Ray &ray, Intersection &its, Float mint, Float m
 					}
 				}
 			} else {
-				if (m_shapes[kdTri.shapeIndex].rayIntersect(ray, searchStart, searchEnd, t)) {
+				if (m_shapes[kdTri.shapeIndex]->rayIntersect(ray, searchStart, searchEnd, t)) {
 					if (shadowRay)
 						return true;
 					foundIntersection = true;

src/tests/test_kd.cpp

@@ -107,7 +107,6 @@ public:
 				Float mint, Float maxt, Float &t, void *tmp) {
 			Float tempT, tempU, tempV;
 			if (m_triangles[idx].rayIntersect(m_vertexBuffer, ray, tempU, tempV, tempT)) {
-				cout << "Got one!" << endl;
 				if (tempT >= mint && tempT <= maxt) {
 					index_type *indexPtr = reinterpret_cast<index_type *>(tmp);
 					Float *floatPtr = reinterpret_cast<Float *>(indexPtr + 1);
@@ -129,17 +128,17 @@ public:
 				Float t, const void *tmp, Intersection &its) const {
 			its.p = ray(t);
 
-			const index_type *indexPtr = reinterpret_cast<const index_type *>(tmp);
+			//const index_type *indexPtr = reinterpret_cast<const index_type *>(tmp);
-			const Float *floatPtr = reinterpret_cast<const Float *>(indexPtr + 1);
+			//const Float *floatPtr = reinterpret_cast<const Float *>(indexPtr + 1);
 
-			const Triangle &tri = m_triangles[*indexPtr];
+			//const Triangle &tri = m_triangles[*indexPtr];
-			const Vertex &v0 = m_vertexBuffer[tri.idx[0]];
+			//const Vertex &v0 = m_vertexBuffer[tri.idx[0]];
-			const Vertex &v1 = m_vertexBuffer[tri.idx[1]];
+			//const Vertex &v1 = m_vertexBuffer[tri.idx[1]];
-			const Vertex &v2 = m_vertexBuffer[tri.idx[2]];
+			//const Vertex &v2 = m_vertexBuffer[tri.idx[2]];
-
-			const Float u = *floatPtr++, v = *floatPtr++;
-			const Vector b(1 - u - v, u, v);
 
+			//const Float u = *floatPtr++, v = *floatPtr++;
+			//const Vector b(1 - u - v, u, v);
+/*
 			its.uv = v0.uv * b.x + v1.uv * b.y + v2.uv * b.z;
 			its.dpdu = v0.dpdu * b.x + v1.dpdu * b.y + v2.dpdu * b.z;
 			its.dpdv = v0.dpdv * b.x + v1.dpdv * b.y + v2.dpdv * b.z;
@@ -150,6 +149,7 @@ public:
 				* dot(its.shFrame.n, its.dpdu));
 			its.shFrame.t = cross(its.shFrame.n, its.shFrame.s);
 			its.wi = its.toLocal(-ray.d);
+*/
 			its.hasUVPartials = false;
 		}
 
@@ -172,18 +172,19 @@ public:
 		tree.build();
 		BSphere bsphere(mesh->getBSphere());
 
+		/*
 		ref<KDTree> oldTree = new KDTree();
 		oldTree->addShape(mesh);
 		oldTree->build();
+		*/
 
 		for (int j=0; j<3; ++j) {
 			ref<Random> random = new Random();
 			ref<Timer> timer = new Timer();
-			size_t nRays = 100, nIntersections = 0, nIntersectionsBF = 0, nIntersectionsOld = 0;
+			size_t nRays = 10000000, nIntersections = 0;
 
 			Log(EInfo, "Bounding sphere: %s", bsphere.toString().c_str());
 			Log(EInfo, "Shooting " SIZE_T_FMT " rays ..", nRays);
-			uint32_t tmp[8];
 
 			for (size_t i=0; i<nRays; ++i) {
 				Point2 sample1(random->nextFloat(), random->nextFloat()),
@@ -193,25 +194,37 @@ public:
 				Ray r(p1, normalize(p2-p1));
 				Intersection its;
 
-				for (uint32_t j=0; j<tree.getPrimitiveCount(); ++j)
-					if (tree.intersect(r, j, r.mint, r.maxt, its.t, tmp) == TriKDTree::EYes)
-						nIntersectionsBF++;
-
 				if (tree.rayIntersect(r, its))
 					nIntersections++;
 
-				if (oldTree->rayIntersect(r, its))
-					nIntersectionsOld++;
 			}
 
-			Log(EInfo, "KD:  Found " SIZE_T_FMT " intersections in %i ms",
+			Log(EInfo, "New: Found " SIZE_T_FMT " intersections in %i ms",
 				nIntersections, timer->getMilliseconds());
-			Log(EInfo, "BF:  Found " SIZE_T_FMT " intersections in %i ms",
+			Log(EInfo, "New: %.3f MRays/s", 
-				nIntersectionsBF, timer->getMilliseconds());
-			Log(EInfo, "Old: Found " SIZE_T_FMT " intersections in %i ms",
-				nIntersectionsOld, timer->getMilliseconds());
-			Log(EInfo, "%.3f MRays/s", 
 				nRays / (timer->getMilliseconds() * (Float) 1000));
+/*
+			random = new Random();
+			timer->reset();
+			nIntersections=0;
+
+			for (size_t i=0; i<nRays; ++i) {
+				Point2 sample1(random->nextFloat(), random->nextFloat()),
+					sample2(random->nextFloat(), random->nextFloat());
+				Point p1 = bsphere.center + squareToSphere(sample1) * bsphere.radius;
+				Point p2 = bsphere.center + squareToSphere(sample2) * bsphere.radius;
+				Ray r(p1, normalize(p2-p1));
+				Intersection its;
+
+				if (oldTree->rayIntersect(r, its))
+					nIntersections++;
+			}
+
+			Log(EInfo, "Old: Found " SIZE_T_FMT " intersections in %i ms",
+				nIntersections, timer->getMilliseconds());
+			Log(EInfo, "Old: %.3f MRays/s", 
+				nRays / (timer->getMilliseconds() * (Float) 1000));
+*/
 		}
 	}
 };