Smoothing perlin noise with trilinear interpolation

perlin.h

@@ -24,11 +24,27 @@ class perlin {
         }
 
         double noise(const point3& p) const {
-            auto i = int(4 * p.x()) & 255;
+            auto u = p.x() - floor(p.x());
-            auto j = int(4 * p.y()) & 255;
+            auto v = p.y() - floor(p.y());
-            auto k = int(4 * p.z()) & 255;
+            auto w = p.z() - floor(p.z());
 
-            return randFloat[permX[i] ^ permY[j] ^ permZ[k]];
+            auto i = int(floor(p.x()));
+            auto j = int(floor(p.y()));
+            auto k = int(floor(p.z()));
+            double c[2][2][2];
+
+            for (int di = 0; di < 2; di++) {
+                for (int dj = 0; dj < 2; dj++) {
+                    for (int dk = 0; dk < 2; dk++) {
+                        c[di][dj][dk] = randFloat[
+                            permX[(i+di) & 255] ^
+                            permY[(j+dj) & 255] ^
+                            permZ[(k+dk) & 255]
+                        ];
+                    }
+                }
+            }
+            return trileanerInterpolation(c, u, v, w);
         }
 
     private: 
@@ -58,6 +74,22 @@ class perlin {
                 p[target] = tmp;
             }
         }
+
+        static double trileanerInterpolation(double c[2][2][2], double u, double v, double w) {
+            auto accum = 0.0;
+
+            for (int i =0; i < 2; i++) {
+                for (int j = 0; j < 2; j++) {
+                    for (int k = 0; k < 2; k++) {
+                        accum += (i*u + (1-i)*(1-u))
+                               * (j*v + (1-j)*(1-v))
+                               * (k*w + (1-k)*(1-w))
+                               * c[i][j][k];
+                    }
+                }
+            }
+            return accum;
+        }
 };
 
 #endif