OpenSCAD NURBS quad surface test

nurbs.scad

function vsum(v,i=0) = len(v)-1 > i ? v[i] + vsum(v, i+1) : v[i];
function N(p, i, U, u) = p == 0 ? 
  U[i] <= u && u <= U[i+1] ? 1 : 0 :
  let(a=u-U[i],b=U[i+p]-U[i],c=U[i+p+1]-u,d=U[i+p+1]-U[i+1]) 
  (b==0?0:a/b*N(p-1,i,U,u)) + (d==0?0:c/d*N(p-1,i+1,U,u));

module nurbs_circle() {
  step = 0.0001;
  // degree of curve
  p = 2; 
  // control points, weights
  P = [[1,0],[0,0],[0.5,sqrt(3)/2],[1,sqrt(3)],[1.5,sqrt(3)/2],[2,0],[1,0]]; 
  W = [   1 , 1/2 ,             1 ,       1/2 ,             1 , 1/2 ,   1 ];
  // knot vector
  U = [0, 0, 0, 1/3, 1/3, 2/3, 2/3, 1, 1, 1];
  points = [for(u=[step:step:U[len(U)-1]])
    let(nv=[for(i=[0:1:len(P)-1]) N(p,i,U,u)*W[i]])
    vsum([for(i=[0:1:len(P)-1]) nv[i]*P[i] ]) /
    vsum([for(i=[0:1:len(P)-1]) nv[i]      ])
  ];
  polygon(points);
}

module nurbs_quad_patch() {
  step = 0.1;
  p = 3;
  P = [
    [[-5, 5, 0],[-3, 5, 0],[-1, 5, 0],[1, 5, 0],[3, 5, 0],[5, 5, 0]],
    [[-5, 3, 0],[-3, 3,-3],[-1, 3,-3],[1, 3,-3],[3, 3,-3],[5, 3, 0]],
    [[-5, 1, 0],[-3, 1,-3],[-1, 1, 3],[1, 1, 3],[3, 1,-3],[5, 1, 0]],
    [[-5,-1, 0],[-3,-1,-3],[-1,-1, 3],[1,-1, 3],[3,-1,-3],[5,-1, 0]],
    [[-5,-3, 0],[-3,-3,-3],[-1,-3,-3],[1,-3,-3],[3,-3,-3],[5,-3, 0]],
    [[-5,-5, 0],[-3,-5, 0],[-1,-5, 0],[1,-5, 0],[3,-5, 0],[5,-5, 0]]
  ];
  W = [
    [1,1,1,1,1,1],
    [1,4,2,2,4,1],
    [1,2,2,2,2,1],
    [1,2,2,2,2,1],
    [1,4,2,2,4,1],
    [1,1,1,1,1,1]
  ];
  U = [0,1,2,3,4,5,6,7,8,9];
  V = [0,1,2,3,4,5,6,7,8,9];
  points = [
    for(u=[U[0]+step:step:U[len(U)-1]-step]) let(ni=[for(i=[0:1:len(P)-1]) N(p,i,U,u)])
      [for(v=[V[0]+step:step:V[len(V)-1]-step]) let(nj=[for(j=[0:1:len(P[0])-1]) N(p,j,V,v) ])
        vsum([for(i=[0:1:len(P)-1],j=[0:1:len(P[0])-1]) ni[i]*nj[j]*W[i][j]*P[i][j] ]) / 
        vsum([for(i=[0:1:len(P)-1],j=[0:1:len(P[0])-1]) ni[i]*nj[j]*W[i][j] ])
      ]
  ];
  
  //echo(points);
  polyhedron(
    points=[for(row=points) each row], 
    faces=[let(l=len(points[0])) for(u=[0:1:len(points)-2], v=[0:1:l-2]) 
      [u*l+v,u*l+v+1,(u+1)*l+v+1,(u+1)*l+v]
    ]
  );
}

nurbs_quad_patch();
//nurbs_circle();