simplex_interpolator: first working version of simplex_interpolator. Tested (with std::cout) for 2D, should work for higher dimensions.

src/lol/math/simplex_interpolator.h

@@ -51,7 +51,7 @@ public:
         vec_t<int, N> index_order = this->GetIndexOrder(decimal_point);
 
         // Resetting decimal point in regular orthonormal coordinates
-        decimal_point = m_base * decimal_point;
+        //decimal_point = m_base * decimal_point;
 
         return this->LastInterp(floor_point, decimal_point, index_order);
     }
@@ -62,21 +62,22 @@ protected:
                         vec_t<float, N> const & decimal_point,
                         vec_t<int, N> index_order) const
     {
-        vec_t<float, N> point_of_interest;
+        vec_t<float, N> point_of_interest(0);
         float result = 0;
 
-        for (int i = 0 ; i < N ; ++i)
+        for (int i = 0 ; i < N+1 ; ++i)
         {
-            float dist = 0.5 - sqlength(decimal_point - point_of_interest);
+            float dist = 0.5 - 0.75f * sqlength(this->m_base * (decimal_point - point_of_interest));
             vec_t<float, N> gradient = this->m_gradients[floor_point];
 
             if (dist > 0)
+                result += dist * dist * dist * dist * dot(gradient, this->m_base * (decimal_point - point_of_interest));
+
+            if (i < N)
             {
-                result += dist * dist * dist * dist * dot(gradient, decimal_point - point_of_interest);
+                point_of_interest[index_order[i]] += 1;
+                floor_point[index_order[i]] = this->Mod(floor_point[index_order[i]] + 1, index_order[i]);
             }
-
-            point_of_interest[index_order[i]] += 1;
-            floor_point[i] = this->Mod(floor_point[i] + 1, index_order[i]);
         }
 
         return result;
@@ -92,7 +93,7 @@ protected:
         {
             for (int j = i + 1 ; j < N ; ++j)
             {
-                if (decimal_point[result[i]] > decimal_point[result[j]])
+                if (decimal_point[result[i]] < decimal_point[result[j]])
                 {
                     // just swapping…
                     result[i] ^= result[j];

src/t/math/simplex_interpolator.cpp

@@ -177,7 +177,7 @@ lolunit_declare_fixture(SimplexInterpolatorTest)
         vec_t<int, N> result = s.GetIndexOrder(vec_ref);
 
         for (int i = 1 ; i < N ; ++i)
-            lolunit_assert(vec_ref[result[i]] >= vec_ref[result[i-1]]);
+            lolunit_assert(vec_ref[result[i]] < vec_ref[result[i-1]]);
     }
 
     lolunit_declare_test(IndexOrderTest)
@@ -225,13 +225,16 @@ lolunit_declare_fixture(SimplexInterpolatorTest)
         simplex_interpolator<2> s;
         s.SetGradients(gradients);
 
-        for (int i = 0 ; i < 10 ; ++i)
+        std::cout << std::endl;
+        for (int i = 0 ; i < 32 ; ++i)
         {
-            for (int j = 0 ; j < 10 ; ++j)
+            for (int j = 0 ; j < 32 ; ++j)
             {
-                std::cout << s.Interp(vec_t<float, 2>({i, j})) << std::endl;
+                std::cout << s.Interp(vec_t<float, 2>({i * 0.1f, j * 0.1f})) << ", ";
             }
+            std::cout << std::endl;
         }
+        std::cout << std::endl;
     }
 
 #if 0