Nicer testcase framework, eigendecomposition support

SConstruct

@@ -431,7 +431,6 @@ if hasQt:
 	qtEnv = mainEnv.Clone()
 	qtEnv.Append(CPPPATH=['src/qtgui'])
 	qtEnv.EnableQt4Modules(['QtGui', 'QtCore', 'QtOpenGL', 'QtXml', 'QtNetwork'])
-	print(qtEnv.Dump())
 	if sys.platform == 'win32':
 		index = qtEnv['CXXFLAGS'].index('_CONSOLE')
 		del qtEnv['CXXFLAGS'][index-1]
@@ -615,9 +614,11 @@ plugins += env.SharedLibrary('plugins/vpl', ['src/integrators/vpl/vpl.cpp'])
 #])
 
 # Testcases
+testEnv = env.Clone()
+testEnv.Append(CPPDEFINES = [['MTS_TESTCASE', '1']])
 for plugin in glob.glob('src/tests/test_*.cpp'):
 	name = os.path.basename(plugin)
-	plugins += env.SharedLibrary('plugins/' + name[0:len(name)-4], plugin)
+	plugins += testEnv.SharedLibrary('plugins/' + name[0:len(name)-4], plugin)
 
 installTargets = []
 

include/mitsuba/core/transform.h

@@ -53,6 +53,9 @@ public:
 	/// Return the determinant of the upper left 3x3 sub-matrix
 	Float det3x3() const;
 
+	/// Perform a symmetric 4x4 eigendecomposition into Q and D.
+	void symmEigenDecomp(Matrix4x4 *Q, Vector4 &d);
+
 	/// Return the inverse of this matrix
 	ref<Matrix4x4> inverse() const;
 

include/mitsuba/render/testcase.h

@@ -23,6 +23,16 @@
 
 MTS_NAMESPACE_BEGIN
 
+#if defined(MTS_TESTCASE)
+/**
+ * When a testcase is being compiled, define the following preprocessor macros for convenience
+ */
+#define assertEquals(expected, actual) assertEqualsImpl(expected, actual, Epsilon, __FILE__, __LINE__)
+#define assertEqualsEpsilon(expected, actual, epsilon) assertEqualsImpl(expected, actual, epsilon, __FILE__, __LINE__)
+#define assertTrue(expr) assertTrueImpl(expr, #expr, __FILE__, __LINE__)
+#define assertFalse(expr) assertFalseImpl(expr, #expr, __FILE__, __LINE__)
+#endif
+
 /** \brief Base class of all testcases. Implementations of this
  * interface can be executed using the 'mtsutil' command. The execution
  * order is as follows: after initializaiton using init(), any tests
@@ -54,6 +64,36 @@ public:
 protected:
 	/// Virtual destructor
 	virtual ~TestCase() { }
+
+	/// Asserts that the two integer values are equal
+	void assertEqualsImpl(int expected, int actual, const char *file, int line);
+
+	/// Asserts that the two floating point values are equal
+	void assertEqualsImpl(Float expected, Float actual, Float epsilon, const char *file, int line);
+
+	/// Asserts that the two 2D vectors are equal
+	void assertEqualsImpl(const Vector2 &expected, const Vector2 &actual, Float epsilon, const char *file, int line);
+
+	/// Asserts that the two 3D vectors are equal
+	void assertEqualsImpl(const Vector &expected, const Vector &actual, Float epsilon, const char *file, int line);
+
+	/// Asserts that the two 4D vectors are equal
+	void assertEqualsImpl(const Vector4 &expected, const Vector4 &actual, Float epsilon, const char *file, int line);
+
+	/// Asserts that the two 2D points are equal
+	void assertEqualsImpl(const Point2 &expected, const Point2 &actual, Float epsilon, const char *file, int line);
+
+	/// Asserts that the two 3D points are equal
+	void assertEqualsImpl(const Point &expected, const Point &actual, Float epsilon, const char *file, int line);
+
+	/// Asserts that the two 4x4 matrices are equal
+	void assertEqualsImpl(const Matrix4x4 *expected, const Matrix4x4 *actual, Float epsilon, const char *file, int line);
+
+	/// Asserts that a condition is true
+	void assertTrueImpl(bool condition, const char *expr, const char *file, int line);
+
+	/// Asserts that a condition is false
+	void assertFalseImpl(bool condition, const char *expr, const char *file, int line);
 protected:
 	int m_executed, m_succeeded;
 };

src/libcore/transform.cpp

@@ -135,7 +135,6 @@ Transform Transform::perspective(Float fov, Float clipNear, Float clipFar) {
 
 	/* Perform a scale so that the field of view is mapped
 	 * to the interval [-1, 1] */
-
 	Float cot = 1.0f / std::tan(degToRad(fov / 2.0f));
 
 	return Transform::scale(Vector(cot, cot, 1.0f)) * Transform(trafo);
@@ -350,7 +349,249 @@ void Matrix4x4::serialize(Stream *stream) const {
 	stream->writeFloatArray(&m[0][0], 16);
 }
 
+static void tred2(Float V[4][4], Float d[4], Float e[4]);
+static void tql2(Float V[4][4], Float d[4], Float e[4]);
 
+void Matrix4x4::symmEigenDecomp(Matrix4x4 *Q, Vector4 &d) {
+	Float e[4];
+	*Q = *this;
+	tred2(Q->m, (Float *) &d, e);
+	tql2(Q->m, (Float *) &d, e);
+}
+
+/* Eigendecomposition code for symmetric 3x3 matrices
+   Copied from the public domain Java Matrix library JAMA. */
+
+inline Float hypot2(Float x, Float y){
+	return std::sqrt(x * x + y * y);
+}
+
+// Symmetric Householder reduction to tridiagonal form.
+static void tred2(Float V[4][4], Float d[4], Float e[4]) {
+	//  This is derived from the Algol procedures tred2 by
+	//  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
+	//  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
+	//  Fortran subroutine in EISPACK.
+
+	const int n = 4;
+	for (int j = 0; j < n; j++)
+		d[j] = V[n - 1][j];
+	
+	// Householder reduction to tridiagonal form.
+	for (int i = n - 1; i > 0; i--) {
+		// Scale to avoid under/overflow.
+		Float scale = 0.0f, h = 0.0f;
+
+		for (int k = 0; k < i; k++)
+			scale = scale + std::abs(d[k]);
+
+		if (scale == 0.0f) {
+			e[i] = d[i - 1];
+			for (int j = 0; j < i; j++) {
+				d[j] = V[i - 1][j];
+				V[i][j] = 0.0f;
+				V[j][i] = 0.0f;
+			}
+		} else {
+			// Generate Householder vector.
+
+			for (int k = 0; k < i; k++) {
+				d[k] /= scale;
+				h += d[k] * d[k];
+			}
+			Float f = d[i - 1], 
+				  g = std::sqrt(h);
+
+			if (f > 0)
+				g = -g;
+			e[i] = scale * g;
+			h = h - f * g;
+			d[i - 1] = f - g;
+
+			for (int j = 0; j < i; j++)
+				e[j] = 0.0f;
+
+			 // Apply similarity transformation to remaining columns.
+			for (int j = 0; j < i; j++) {
+				f = d[j];
+				V[j][i] = f;
+				g = e[j] + V[j][j] * f;
+				for (int k = j + 1; k <= i - 1; k++) {
+					g += V[k][j] * d[k];
+					e[k] += V[k][j] * f;
+				}
+				e[j] = g;
+			}
+
+			f = 0.0f;
+			for (int j = 0; j < i; j++) {
+				e[j] /= h;
+				f += e[j] * d[j];
+			}
+			Float hh = f / (h + h);
+
+			for (int j = 0; j < i; j++)
+				e[j] -= hh * d[j];
+			for (int j = 0; j < i; j++) {
+				f = d[j];
+				g = e[j];
+				for (int k = j; k <= i - 1; k++)
+					V[k][j] -= (f * e[k] + g * d[k]);
+				d[j] = V[i - 1][j];
+				V[i][j] = 0.0f;
+			}
+		}
+		d[i] = h;
+	}
+
+	// Accumulate transformations.
+	for (int i = 0; i < n - 1; i++) {
+		V[n - 1][i] = V[i][i];
+		V[i][i] = 1.0f;
+		Float h = d[i + 1];
+
+		if (h != 0.0f) {
+			for (int k = 0; k <= i; k++)
+				d[k] = V[k][i + 1] / h;
+			for (int j = 0; j <= i; j++) {
+				Float g = 0.0f;
+
+				for (int k = 0; k <= i; k++)
+					g += V[k][i + 1] * V[k][j];
+				for (int k = 0; k <= i; k++)
+					V[k][j] -= g * d[k];
+			}
+		}
+		for (int k = 0; k <= i; k++)
+			V[k][i + 1] = 0.0f;
+	}
+	for (int j = 0; j < n; j++) {
+		d[j] = V[n - 1][j];
+		V[n - 1][j] = 0.0f;
+	}
+	V[n - 1][n - 1] = 1.0f;
+	e[0] = 0.0;
+}
+
+// Symmetric tridiagonal QL algorithm.
+static void tql2(Float V[4][4], Float d[4], Float e[4]) {
+	//  This is derived from the Algol procedures tql2, by
+	//  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
+	//  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
+	//  Fortran subroutine in EISPACK.
+
+	const int n = 4;
+
+	for (int i = 1; i < n; i++)
+		e[i - 1] = e[i];
+
+	e[n - 1] = 0.0f;
+	Float f = 0.0f, tst1 = 0.0f;
+
+#if defined(SINGLE_PRECISION)
+	Float eps = pow(2.0f, -23.0f);
+#else
+	Float eps = pow(2.0, -52.0);
+#endif
+
+	for (int l = 0; l < n; l++) {
+		// Find small subdiagonal element
+		tst1 = std::max(tst1, std::abs(d[l]) + std::abs(e[l]));
+		int m = l;
+
+		while (m < n) {
+			if (std::abs(e[m]) <= eps * tst1)
+				break;
+			m++;
+		}
+
+		// If m == l, d[l] is an eigenvalue,
+		// otherwise, iterate.
+		if (m > l) {
+			int iter = 0;
+
+			do {
+				iter = iter + 1;    // (Could check iteration count here.)
+
+				// Compute implicit shift
+				Float g = d[l];
+				Float p = (d[l + 1] - g) / (2.0f * e[l]);
+				Float r = hypot2(p, 1.0f);
+
+				if (p < 0)
+					r = -r;
+				d[l] = e[l] / (p + r);
+				d[l + 1] = e[l] * (p + r);
+				Float dl1 = d[l + 1];
+
+				Float h = g - d[l];
+
+				for (int i = l + 2; i < n; i++)
+					d[i] -= h;
+				f = f + h;
+
+				// Implicit QL transformation.
+				p = d[m];
+				Float c = 1.0f;
+				Float c2 = c, c3 = c;
+				Float el1 = e[l + 1];
+				Float s = 0.0f, s2 = 0.0f;
+
+				for (int i = m - 1; i >= l; i--) {
+					c3 = c2;
+					c2 = c;
+					s2 = s;
+					g = c * e[i];
+					h = c * p;
+					r = hypot2(p, e[i]);
+					e[i + 1] = s * r;
+					s = e[i] / r;
+					c = p / r;
+					p = c * d[i] - s * g;
+					d[i + 1] = h + s * (c * g + s * d[i]);
+
+					// Accumulate transformation.
+					for (int k = 0; k < n; k++) {
+						h = V[k][i + 1];
+						V[k][i + 1] =
+							s * V[k][i] + c * h;
+						V[k][i] = c * V[k][i] - s * h;
+					}
+				}
+				p = -s * s2 * c3 * el1 * e[l] / dl1;
+				e[l] = s * p;
+				d[l] = c * p;
+				// Check for convergence.
+			} while (std::abs(e[l]) > eps * tst1);
+		}
+		d[l] = d[l] + f;
+		e[l] = 0.0f;
+	}
+
+	// Sort eigenvalues and corresponding vectors.
+	for (int i = 0; i < n - 1; i++) {
+		int k = i;
+
+		Float p = d[i];
+
+		for (int j = i + 1; j < n; j++) {
+			if (d[j] < p) {
+				k = j;
+				p = d[j];
+			}
+		}
+
+		if (k != i) {
+			d[k] = d[i];
+			d[i] = p;
+			for (int j = 0; j < n; j++) {
+				p = V[j][i];
+				V[j][i] = V[j][k];
+				V[j][k] = p;
+			}
+		}
+	}
+}
 
 MTS_IMPLEMENT_CLASS(Matrix4x4, false, Object)
 MTS_NAMESPACE_END

src/librender/testcase.cpp

@@ -26,6 +26,89 @@ MTS_NAMESPACE_BEGIN
 void TestCase::init() { }
 void TestCase::shutdown() { }
 
+void TestCase::assertTrueImpl(bool value, const char *expr, const char *file, int line) {
+	if (!value)
+		Thread::getThread()->getLogger()->log(EError, NULL, file, line, "Assertion '%s == true' failed!", expr);
+}
+
+void TestCase::assertFalseImpl(bool value, const char *expr, const char *file, int line) {
+	if (value)
+		Thread::getThread()->getLogger()->log(EError, NULL, file, line, "Assertion '%s == false' failed!", expr);
+}
+
+void TestCase::assertEqualsImpl(int expected, int actual, const char *file, int line) {
+	if (expected != actual)
+		Thread::getThread()->getLogger()->log(EError, NULL, file, line, "Assertion failure: "
+			"expected integer value %i, got %i.", expected, actual);
+}
+
+void TestCase::assertEqualsImpl(Float expected, Float actual, Float epsilon, const char *file, int line) {
+	if (std::abs(expected-actual) > epsilon)
+		Thread::getThread()->getLogger()->log(EError, NULL, file, line, "Assertion failure: "
+			"expected floating point value %f, got %f.", expected, actual);
+}
+
+void TestCase::assertEqualsImpl(const Vector2 &expected, const Vector2 &actual, Float epsilon, const char *file, int line) {
+	bool match = true;
+	for (int i=0; i<2; ++i)
+		if (std::abs(expected[i]-actual[i]) > epsilon)
+			match = false;
+	if (!match)
+		Thread::getThread()->getLogger()->log(EError, NULL, file, line, "Assertion failure: "
+			"expected vector %s, got %s.", expected.toString().c_str(), actual.toString().c_str());
+}
+
+void TestCase::assertEqualsImpl(const Point2 &expected, const Point2 &actual, Float epsilon, const char *file, int line) {
+	bool match = true;
+	for (int i=0; i<2; ++i)
+		if (std::abs(expected[i]-actual[i]) > epsilon)
+			match = false;
+	if (!match)
+		Thread::getThread()->getLogger()->log(EError, NULL, file, line, "Assertion failure: "
+			"expected point %s, got %s.", expected.toString().c_str(), actual.toString().c_str());
+}
+
+void TestCase::assertEqualsImpl(const Vector &expected, const Vector &actual, Float epsilon, const char *file, int line) {
+	bool match = true;
+	for (int i=0; i<3; ++i)
+		if (std::abs(expected[i]-actual[i]) > epsilon)
+			match = false;
+	if (!match)
+		Thread::getThread()->getLogger()->log(EError, NULL, file, line, "Assertion failure: "
+			"expected vector %s, got %s.", expected.toString().c_str(), actual.toString().c_str());
+}
+
+void TestCase::assertEqualsImpl(const Point &expected, const Point &actual, Float epsilon, const char *file, int line) {
+	bool match = true;
+	for (int i=0; i<3; ++i)
+		if (std::abs(expected[i]-actual[i]) > epsilon)
+			match = false;
+	if (!match)
+		Thread::getThread()->getLogger()->log(EError, NULL, file, line, "Assertion failure: "
+			"expected point %s, got %s.", expected.toString().c_str(), actual.toString().c_str());
+}
+
+void TestCase::assertEqualsImpl(const Vector4 &expected, const Vector4 &actual, Float epsilon, const char *file, int line) {
+	bool match = true;
+	for (int i=0; i<4; ++i)
+		if (std::abs(expected[i]-actual[i]) > epsilon)
+			match = false;
+	if (!match)
+		Thread::getThread()->getLogger()->log(EError, NULL, file, line, "Assertion failure: "
+			"expected vector %s, got %s.", expected.toString().c_str(), actual.toString().c_str());
+}
+
+void TestCase::assertEqualsImpl(const Matrix4x4 *expected, const Matrix4x4 *actual, Float epsilon, const char *file, int line) {
+	bool match = true;
+	for (int i=0; i<4; ++i)
+		for (int j=0; j<4; ++j)
+			if (std::abs(expected->m[i][j]-actual->m[i][j]) > epsilon)
+				match = false;
+	if (!match)
+		Thread::getThread()->getLogger()->log(EError, NULL, file, line, "Assertion failure: "
+			"expected matrix %s, got %s.", expected->toString().c_str(), actual->toString().c_str());
+}
+
 struct Sample {
 	Float value;
 	Float variance;

src/tests/test_la.cpp

@@ -0,0 +1,56 @@
+/*
+    This file is part of Mitsuba, a physically based rendering system.
+
+    Copyright (c) 2007-2010 by Wenzel Jakob and others.
+
+    Mitsuba is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License Version 3
+    as published by the Free Software Foundation.
+
+    Mitsuba is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program. If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <mitsuba/render/testcase.h>
+
+MTS_NAMESPACE_BEGIN
+
+class TestLinearAlgebra : public TestCase {
+public:
+	MTS_BEGIN_TESTCASE()
+	MTS_DECLARE_TEST(test01_eigenDecomp)
+	MTS_END_TESTCASE()
+
+	void test01_eigenDecomp() {
+		ref<Matrix4x4> A = new Matrix4x4(
+			1.4541, 1.1233, 1.2407, 1.2548,
+			1.1233, 0.2597, 0.3819, 1.3917,
+			1.2407, 0.3819, 1.1552, 1.1048,
+			1.2548, 1.3917, 1.1048, 1.4712
+		);
+		ref<Matrix4x4> Q = new Matrix4x4();
+		Vector4 d;
+		A->symmEigenDecomp(Q, d);
+
+		Vector4 refD(-0.823889076095475, 0.130902702868822,
+			0.557486242256414, 4.47570013097024);
+
+		ref<Matrix4x4> refQ = new Matrix4x4(
+			0.294383137629217, 0.746207030711883, 0.191065628242818, 0.565692108217809,
+			-0.789565156329591, 0.139585328270248, -0.459560990857043, 0.381977087957391,
+			-0.286639718342894, -0.415237779451033, 0.738328701906588, 0.447533223711729,
+			0.455810381691897, -0.501266984714151, -0.45515749947419,  0.577754235510108
+		);
+
+		assertEqualsEpsilon(d, refD, 1e-6);
+		assertEqualsEpsilon(Q, refQ, 1e-6);
+	}
+};
+
+MTS_EXPORT_TESTCASE(TestLinearAlgebra, "Testcase for Linear Algebra routines")
+MTS_NAMESPACE_END

src/tests/test_samplers.cpp

@@ -46,7 +46,7 @@ public:
 		sampler->generate();
 		for (int i=0; i<5; ++i) {
 			for (int j=0; j<5; ++j)
-				Assert(std::abs(sampler->next1D() - comparison[pos++]) < 1e-7);
+				assertEqualsEpsilon(sampler->next1D(), comparison[pos++], 1e-7);
 			sampler->advance();
 		}
 	}
@@ -70,7 +70,7 @@ public:
 		sampler->generate();
 		for (int i=0; i<5; ++i) {
 			for (int j=0; j<6; ++j) 
-				Assert(std::abs(sampler->next1D() - comparison[pos++]) < 1e-7);
+				assertEqualsEpsilon(sampler->next1D(), comparison[pos++], 1e-7);
 			sampler->advance();
 		}
 	}
@@ -79,7 +79,7 @@ public:
 		Float x = 0.0f;
 
 		for (int i=0; i<20; ++i) {
-			Assert(x == radicalInverse(2, i));
+			assertEqualsEpsilon(x, radicalInverse(2, i), 0);
 			x = radicalInverseIncremental(2, x);
 		}
 	}

src/tests/test_sh.cpp

@@ -56,7 +56,7 @@ public:
 
 			Float value1 = vec1.eval(dir2);
 			Float value2 = vec2.eval(dir1);
-			Assert(std::abs(value1-value2) < Epsilon);
+			assertEquals(value1, value2);
 		}
 	}
 
@@ -91,10 +91,10 @@ public:
 			Float relerr = std::abs(pdf1-pdf2)/pdf2;
 			if (pdf2 > 0.01) {
 				accum += relerr; ++nInAvg;
-				Assert(relerr < 0.08);
+				assertTrue(relerr < 0.08);
 			}
 		}
-		Assert(accum / nInAvg < 0.01);
+		assertTrue(accum / nInAvg < 0.01);
 	}
 };