diff --git a/SConstruct b/SConstruct index 674c4d5c..1e00bf51 100644 --- a/SConstruct +++ b/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 = [] diff --git a/include/mitsuba/core/transform.h b/include/mitsuba/core/transform.h index 1788a85b..e1f34469 100644 --- a/include/mitsuba/core/transform.h +++ b/include/mitsuba/core/transform.h @@ -53,9 +53,12 @@ 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 inverse() const; - + /// Return the transpose of this matrix ref transpose() const; diff --git a/include/mitsuba/render/testcase.h b/include/mitsuba/render/testcase.h index 2b672f04..064303ee 100644 --- a/include/mitsuba/render/testcase.h +++ b/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; }; diff --git a/src/libcore/transform.cpp b/src/libcore/transform.cpp index e953b8d3..68557ea8 100644 --- a/src/libcore/transform.cpp +++ b/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 diff --git a/src/librender/testcase.cpp b/src/librender/testcase.cpp index 71139eae..2bc4847b 100644 --- a/src/librender/testcase.cpp +++ b/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; diff --git a/src/tests/test_la.cpp b/src/tests/test_la.cpp new file mode 100644 index 00000000..fcae0623 --- /dev/null +++ b/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 . +*/ + +#include + +MTS_NAMESPACE_BEGIN + +class TestLinearAlgebra : public TestCase { +public: + MTS_BEGIN_TESTCASE() + MTS_DECLARE_TEST(test01_eigenDecomp) + MTS_END_TESTCASE() + + void test01_eigenDecomp() { + ref 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 Q = new Matrix4x4(); + Vector4 d; + A->symmEigenDecomp(Q, d); + + Vector4 refD(-0.823889076095475, 0.130902702868822, + 0.557486242256414, 4.47570013097024); + + ref 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 diff --git a/src/tests/test_samplers.cpp b/src/tests/test_samplers.cpp index 4d8221a2..f9587e4d 100644 --- a/src/tests/test_samplers.cpp +++ b/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); } } diff --git a/src/tests/test_sh.cpp b/src/tests/test_sh.cpp index f6306b67..c29ecfa0 100644 --- a/src/tests/test_sh.cpp +++ b/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); } };