From 97de8946e93c57eb685b2042c2703aeeec46bd73 Mon Sep 17 00:00:00 2001
From: Guillaume Bittoun <guillaume.bittoun@gmail.com>
Date: Thu, 4 Dec 2014 10:19:14 +0000
Subject: [PATCH] simplex_interpolator: fixing bug on simplex matrix generation

---
 src/lol/math/simplex_interpolator.h | 14 ++++-----
 src/t/math/simplex_interpolator.cpp | 49 +++++++++++++++++++++++++++++
 2 files changed, 56 insertions(+), 7 deletions(-)

diff --git a/src/lol/math/simplex_interpolator.h b/src/lol/math/simplex_interpolator.h
index 5d5f70b3..e468c445 100644
--- a/src/lol/math/simplex_interpolator.h
+++ b/src/lol/math/simplex_interpolator.h
@@ -127,22 +127,22 @@ protected:
         /*  Matrix coordinate transformation to skew simplex referential is done
             by inversing the base matrix M which is written as follows: 
 
-             M = | a b b b … |        M^(-1) = | c d d d … |
-                 | b a b b … |                 | d c d d … |
-                 | b b a b … |                 | d d c d … |
-                 | …         |                 | …         |
+            M = | a b b b … |        M^(-1) = | c d d d … |
+                | b a b b … |                 | d c d d … |
+                | b b a b … |                 | d d c d … |
+                | …         |                 | …         |
 
             where a and b are computed below ↴
         */
 
-        T b = (1 - sqrt((N+1)))/(N * N * N));
-        T a = b + sqrt((N+1)/N);
+        T b = (1 - sqrt((T)(N+1)))/sqrt((T)(N * N * N));
+        T a = b + sqrt((N+1)/(T)N);
 
         // determinant of matrix M
         T det = a * (a + (N-2) * b) - (b * b) * (N-1);
 
         T c = (a + (N-2)*b) / det;
-        T d = b / det;
+        T d = -b / det;
 
         for (int i = 0; i < N; ++i)
         {
diff --git a/src/t/math/simplex_interpolator.cpp b/src/t/math/simplex_interpolator.cpp
index 818c00f4..aba46ff4 100644
--- a/src/t/math/simplex_interpolator.cpp
+++ b/src/t/math/simplex_interpolator.cpp
@@ -12,17 +12,65 @@
 
 #include <lol/engine-internal.h>
 
+#include <lol/math/simplex_interpolator.h>
+
 #include <lolunit.h>
 
 namespace lol
 {
 
+template<int N, typename T = float>
+class test_interpolator : public simplex_interpolator<N, T>
+{
+public:
+
+    test_interpolator(arraynd<N, T> const & samples) :
+        simplex_interpolator<N, T>(samples)
+    {
+    }
+
+
+    void DumpMatrix()
+    {
+        std::cout << std::endl;
+        for (int i = 0; i < N; ++i)
+        {
+            for (int j = 0; j < N; ++j)
+            {
+                std::cout << this->m_base[i][j] << ", ";
+            }
+            std::cout << std::endl;
+        }
+
+        std::cout << std::endl;
+
+        for (int i = 0; i < N; ++i)
+        {
+            for (int j = 0; j < N; ++j)
+            {
+                std::cout << this->m_base_inverse[i][j] << ", ";
+            }
+            std::cout << std::endl;
+        }
+    }
+};
+
 lolunit_declare_fixture(SimplexInterpolatorTest)
 {
     void SetUp() {}
 
     void TearDown() {}
 
+    lolunit_declare_test(CoordinateMatrixTest)
+    {
+        test_interpolator<2> s2({{1.f}});
+        s2.DumpMatrix();
+
+        test_interpolator<3> s3(arraynd<3, float>({{{1.f}}}));
+        s3.DumpMatrix();
+    }
+
+#if 0
     lolunit_declare_test(FloatGridPoints2D1x1)
     {
         simplex_interpolator<2> s({{1.f}});
@@ -151,6 +199,7 @@ private:
 
         return (float)i * vi + (float)j * vj;
     }
+#endif
 };
 
 }