Bicubic interpolation of value noise for OpenSCAD

bicubic_interp_noise.scad

// Bicubic interpolation of value noise
//   by Hans Loeblich
// Reference: https://en.wikipedia.org/wiki/Bicubic_interpolation

dims = [4,4];
min = 0;
max = 2;
seed = 0;


bicubic = create_bicubic_interp(dims, min, max, seed);

step = 1/32;
for(x=[0:step:dims[0]-step],y=[0:step:dims[1]-step])
  translate([x,y, bicubic(x,y)])
    cube(step);


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

function create_2Dlookup(dims, min=0, max=1, seed=0) = 
  function(x,y) rands(min, max, 1, seed + dims[1]*mod(x,dims[0]) + mod(y,dims[1]))[0];

function create_get_bicubic_coefs(dims, min, max, seed) =  let(
  m0 =
  [[ 1, 0, 0, 0],
   [ 0, 0, 1, 0],
   [-3, 3,-2,-1],
   [ 2,-2, 1, 1]],
  m1 =
  [[ 1, 0,-3, 2],
   [ 0, 0, 3,-2],
   [ 0, 1,-2, 1],
   [ 0, 0,-1, 1]],
    f = create_2Dlookup(dims, min,max,seed),
   fx = function(x,y) ( f(x+1,y) -  f(x-1,y))/2,
   fy = function(x,y) ( f(x,y+1) -  f(x,y-1))/2,
  fxy = function(x,y) (fx(x,y+1) - fx(x,y-1))/2,
) function(x,y) let(
  mf = [
    [f (x  ,y  ),f (x  ,y+1), fy(x  ,y  ), fy(x  ,y+1)],
    [f (x+1,y  ),f (x+1,y+1), fy(x+1,y  ), fy(x+1,y+1)],
    [fx(x  ,y  ),fx(x  ,y+1),fxy(x  ,y  ),fxy(x  ,y+1)],
    [fx(x+1,y  ),fx(x+1,y+1),fxy(x+1,y  ),fxy(x+1,y+1)],
  ]
) m0*mf*m1;

function create_bicubic_interp(dims, min, max, seed) = let(
  get_bicubic_coefs = create_get_bicubic_coefs(dims, min, max, seed),
  coef_lut = [ for(x=[0:1:dims[0]-1]) [ for(y=[0:1:dims[1]-1]) get_bicubic_coefs(x,y) ] ],
  dx = dims[0], dy = dims[1],
) function(x,y) let(
  xi = mod(floor(x),dx), yi = mod(floor(y),dy),
  tx = fract(x), ty = fract(y),
  mx = [1,tx,tx*tx,tx*tx*tx], my=[1,ty,ty*ty,ty*ty*ty], mc = coef_lut[xi][yi]
) mx * mc * my;