Tricubic interpolation of 3D noise for OpenSCAD

tricubic_interp_noise.scad

dims = [8,8,3];
min = 0;
max = 2;
seed = 0;

tricubic = create_tricubic_interp(dims, min, max, seed);

step = 1/8;
frames = 64;
z = dims[2] * $t;

for(y=[0:step:dims[1]-step],x=[0:step:dims[0]-step])
  translate([x,y, tricubic(x,y,z)]) cube(step);


function mod(x,y) = let(r=x%y) r<0 ? r+y : r;
function fract(x) = x-floor(x);

function create_3Dlookup(dims, min=0, max=1, seed=0) = let(lx=dims[0], ly=dims[1], lz=dims[2])
  function(x,y,z) rands(min, max, 1, seed + lz*(ly*mod(x,lx) + mod(y,ly)) + mod(z,lz))[0];

/*
function create_tricubic_interp(dims, min, max, seed) = let(
  get_vals = function(x,y,z) [f(x,y,z), (f(x+1,y,z) - f(x-1,y,z))/2, (f(x,y+1,z) - f(x,y-1,z))/2, (f(x,y,z+1) - f(x,y,z-1))/2],
  lx = dims[0], ly = dims[1], lz = dims[2],
  mc = [[0,1,0,0],[-0.5,0,0.5,0],[1,-2.5,2,-0.5],[-0.5,1.5,-1.5,0.5]],
  s = create_3Dlookup(dims, min,max,seed),
  v = [for(x=[0:lx-1]) [for(y=[0:ly-1]) [for(z=[0:lz-1]) s(x,y,z) ] ] ],
  t = function(i,j,z) let(zi=floor(z), zf=z-zi, vv=v[i][j]) [1,zf,zf*zf,zf*zf*zf] * mc * [vv[mod(zi-1,lz)],vv[mod(zi,lz)],vv[mod(zi+1,lz)],vv[mod(zi+2,lz)]],
  u = function(i,y,z) let(yi=floor(y), yf=y-yi) [1,yf,yf*yf,yf*yf*yf] * mc * [t(i,mod(yi-1,ly),z),t(i,mod(yi,ly),z),t(i,mod(yi+1,ly),z),t(i,mod(yi+2,ly),z)],
)     function(x,y,z) let(xi=floor(x), xf=x-xi) [1,xf,xf*xf,xf*xf*xf] * mc * [u(mod(xi-1,lx),y,z),u(mod(xi,lx),y,z),u(mod(xi+1,lx),y,z),u(mod(xi+2,lx),y,z)];
*/



//*
// same as above with functions inlined, much less reduntant calculations.
function create_tricubic_interp(dims, min, max, seed) = let(
  get_vals = function(x,y,z) [f(x,y,z), (f(x+1,y,z) - f(x-1,y,z))/2, (f(x,y+1,z) - f(x,y-1,z))/2, (f(x,y,z+1) - f(x,y,z-1))/2],
  lx = dims[0], ly = dims[1], lz = dims[2],
  mc = [[0,1,0,0],[-0.5,0,0.5,0],[1,-2.5,2,-0.5],[-0.5,1.5,-1.5,0.5]],
  s = create_3Dlookup(dims, min,max,seed),
  v = [for(x=[0:lx-1]) [for(y=[0:ly-1]) [for(z=[0:lz-1]) s(x,y,z) ] ] ]
) function(x,y,z) let(
    xi=floor(x), xf=x-xi, 
    yi=floor(y), yf=y-yi, 
    zi=floor(z), zf=z-zi,
    i0=mod(xi-1,lx), i1=mod(xi,lx), i2=mod(xi+1,lx), i3=mod(xi+2,lx),
    j0=mod(yi-1,ly), j1=mod(yi,ly), j2=mod(yi+1,ly), j3=mod(yi+2,ly),
    k0=mod(zi-1,lz), k1=mod(zi,lz), k2=mod(zi+1,lz), k3=mod(zi+2,lz),
    my = [1,yf,yf*yf,yf*yf*yf] * mc,
    mz = [1,zf,zf*zf,zf*zf*zf] * mc
) [1,xf,xf*xf,xf*xf*xf] * mc* [ my * [
  mz * [v[i0][j0][k0],v[i0][j0][k1],v[i0][j0][k2],v[i0][j0][k3]],
  mz * [v[i0][j1][k0],v[i0][j1][k1],v[i0][j1][k2],v[i0][j1][k3]],
  mz * [v[i0][j2][k0],v[i0][j2][k1],v[i0][j2][k2],v[i0][j2][k3]],
  mz * [v[i0][j3][k0],v[i0][j3][k1],v[i0][j3][k2],v[i0][j3][k3]]
], my * [
  mz * [v[i1][j0][k0],v[i1][j0][k1],v[i1][j0][k2],v[i1][j0][k3]],
  mz * [v[i1][j1][k0],v[i1][j1][k1],v[i1][j1][k2],v[i1][j1][k3]],
  mz * [v[i1][j2][k0],v[i1][j2][k1],v[i1][j2][k2],v[i1][j2][k3]],
  mz * [v[i1][j3][k0],v[i1][j3][k1],v[i1][j3][k2],v[i1][j3][k3]]
], my * [
  mz * [v[i2][j0][k0],v[i2][j0][k1],v[i2][j0][k2],v[i2][j0][k3]],
  mz * [v[i2][j1][k0],v[i2][j1][k1],v[i2][j1][k2],v[i2][j1][k3]],
  mz * [v[i2][j2][k0],v[i2][j2][k1],v[i2][j2][k2],v[i2][j2][k3]],
  mz * [v[i2][j3][k0],v[i2][j3][k1],v[i2][j3][k2],v[i2][j3][k3]]
], my * [
  mz * [v[i3][j0][k0],v[i3][j0][k1],v[i3][j0][k2],v[i3][j0][k3]],
  mz * [v[i3][j1][k0],v[i3][j1][k1],v[i3][j1][k2],v[i3][j1][k3]],
  mz * [v[i3][j2][k0],v[i3][j2][k1],v[i3][j2][k2],v[i3][j2][k3]],
  mz * [v[i3][j3][k0],v[i3][j3][k1],v[i3][j3][k2],v[i3][j3][k3]]
]];