Contribution by Prashanth Udupa on Bezier Surfaces in VTK Designer 2: Bezier_VTKD2.pdf
The Bézier surface is formed as the cartesian product of the blending functions of two orthogonal Bézier curves.

Where Pi,j is the i,jth control point. There are Ni+1 and Nj+1 control points in the i and j directions respectively.
The corresponding properties of the Bézier
curve apply to the Bézier surface.
- The surface does not in general pass through the control points
except for the corners of the control point grid.
- The surface is contained within the convex hull of the control
points.
Along the edges of the grid patch the Bézier surface matches that of a Bézier curve through the control points along that edge.

Closed surfaces can be formed by setting the last control point equal to the first. If the tangents also match between the first two and last two control points then the closed surface will have first order continuity.

While a cylinder/cone can be formed from a Bézier surface, it is not possible to form a sphere.
The following source code generates the surface shown in the first example above. It is provided for illustration only, the headers and prototype files are not given.
#define NI 5
#define NJ 4
XYZ inp[NI+1][NJ+1];
#define RESOLUTIONI 10*NI
#define RESOLUTIONJ 10*NJ
XYZ outp[RESOLUTIONI][RESOLUTIONJ];
int main(argc,argv)
int argc;
char **argv;
{
int i,j,ki,kj;
double mui,muj,bi,bj;
/* Create a random surface */
srandom(1111);
for (i=0;i<=NI;i++) {
for (j=0;j<=NJ;j++) {
inp[i][j].x = i;
inp[i][j].y = j;
inp[i][j].z = (random() % 10000) / 5000.0 - 1;
}
}
for (i=0;i<RESOLUTIONI;i++) {
mui = i / (double)(RESOLUTIONI-1);
for (j=0;j<RESOLUTIONJ;j++) {
muj = j / (double)(RESOLUTIONJ-1);
outp[i][j].x = 0;
outp[i][j].y = 0;
outp[i][j].z = 0;
for (ki=0;ki<=NI;ki++) {
bi = BezierBlend(ki,mui,NI);
for (kj=0;kj<=NJ;kj++) {
bj = BezierBlend(kj,muj,NJ);
outp[i][j].x += (inp[ki][kj].x * bi * bj);
outp[i][j].y += (inp[ki][kj].y * bi * bj);
outp[i][j].z += (inp[ki][kj].z * bi * bj);
}
}
}
}
printf("LIST\n");
/* Display the surface, in this case in OOGL format for GeomView */
printf("{ = CQUAD\n");
for (i=0;i<RESOLUTIONI-1;i++) {
for (j=0;j<RESOLUTIONJ-1;j++) {
printf("%g %g %g 1 1 1 1\n",
outp[i][j].x, outp[i][j].y, outp[i][j].z);
printf("%g %g %g 1 1 1 1\n",
outp[i][j+1].x, outp[i][j+1].y, outp[i][j+1].z);
printf("%g %g %g 1 1 1 1\n",
outp[i+1][j+1].x,outp[i+1][j+1].y,outp[i+1][j+1].z);
printf("%g %g %g 1 1 1 1\n",
outp[i+1][j].x, outp[i+1][j].y, outp[i+1][j].z);
}
}
printf("}\n");
/* Control point polygon */
for (i=0;i<NI;i++) {
for (j=0;j<NJ;j++) {
printf("{ = SKEL 4 1 \n");
printf("%g %g %g \n",inp[i][j].x,inp[i][j].y,inp[i][j].z);
printf("%g %g %g \n",inp[i][j+1].x,inp[i][j+1].y,inp[i][j+1].z);
printf("%g %g %g \n",inp[i+1][j+1].x,inp[i+1][j+1].y,inp[i+1][j+1].z);
printf("%g %g %g \n",inp[i+1][j].x,inp[i+1][j].y,inp[i+1][j].z);
printf("5 0 1 2 3 0\n");
printf("}\n");
}
}
}
This function computes the blending function as used in the Bézier surface code above. It is written for clarity, not efficiency. Normally, if the number of control points is constant, the blending function would be calculated once for each desired value of mu.
double BezierBlend(k,mu,n)
int k;
double mu;
int n;
{
int nn,kn,nkn;
double blend=1;
nn = n;
kn = k;
nkn = n - k;
while (nn >= 1) {
blend *= nn;
nn--;
if (kn > 1) {
blend /= (double)kn;
kn--;
}
if (nkn > 1) {
blend /= (double)nkn;
nkn--;
}
}
if (k > 0)
blend *= pow(mu,(double)k);
if (n-k > 0)
blend *= pow(1-mu,(double)(n-k));
return(blend);
}