module print_points(vector) {
    color("red")
    for(i = [0:1:len(vector)-1]) 
        translate(vector[i]) {
            sphere(1);
            text(str(i),size=4);
        }
}


function smooth_closed(path,it=0)=
    it == 0 ?
        path :
        smooth_closed(subdivide_closed(path),it-1);
  

function smooth_open(path,it=0)=
    it == 0 ?
        path :
        smooth_open(subdivide_open(path),it-1);


function subdivide_closed(path)=
    let(n=len(path)) 
        flatten(concat([for (i=[0:1:n-1])[      path[i],
          interpolateClosed(path,n,i)
        ]]));

function subdivide_open(path) =
  let(n = len(path))
  flatten(concat([for (i = [0 : 1 : n-1])
    i < n-1? 
      // Emit the current point and the one halfway between current and next.
      [path[i], interpolateOpen(path, n, i)]
    :
      // We're at the end, so just emit the last point.
      [path[i]]
  ]));

weight = [-1, 8, 8, -1] / 14;
weight0 = [6, 11, -1] / 16;
weight2 = [1, 1] / 2;

// Interpolate on an open-ended path, with discontinuity at start and end.
// Returns a point between points i and i+1, weighted.
function interpolateOpen(path, n, i) =
  i == 0? 
    n == 2?
      path[i]     * weight2[0] +
      path[i + 1] * weight2[1]
    :
      path[i]     * weight0[0] +
      path[i + 1] * weight0[1] +
      path[i + 2] * weight0[2]
  : i < n - 2?
    path[i - 1] * weight[0] +
    path[i]     * weight[1] +
    path[i + 1] * weight[2] +
    path[i + 2] * weight[3]
  : i < n - 1?
    path[i - 1] * weight0[2] +
    path[i]     * weight0[1] +
    path[i + 1] * weight0[0]
  : 
    path[i];

function interpolateClosed(path, n, i) =
  path[(i + n - 1) % n] * weight[0] +
  path[i]               * weight[1] +
  path[(i + 1) % n]     * weight[2] +
  path[(i + 2) % n]     * weight[3] ;

function flatten(list) = [ for (i = list, v = i) v ];