One student e-mailed me to explain the B-spline fitting algorithm. Here is my reply:

-------
Sanjay,

Did you look at the PowerPoint presentation on the algorithm?
http://www.cs.columbia.edu/~laza/heart3D/slides2.ppt

The important page is
http://www.cs.columbia.edu/~laza/heart3D/B-spline_surface_modelling_files/v3_slide0010.htm

Start from the Matlab file sREADME.m.

If you look at the Matlab code, the file of interest is sMBA.m. This 
function iteratively fits a surface over a plane: y = f(x,y). In the 
first iteration, the surface is fitted with one-point spline, which 
is a horizontal plane at the average hight.

Fitting each iteration is done in file sBA.m. When you open the file, 
everything between lines 31 and 46 is stuff for fitting spherical 
objects in polar coordinate system. I am not sure if you need that. 
For each cloud point, the algorithm adds a 4x4 spline to the matrix 
representing z=f(x,y).

For instance,
- your (x,y) plane has 16x16 points  ==>  your matrix z=zeros(16,16))
- you have a bunch of (x,y,z) points in 3 vectors (nx1 matrix)
- one of the points is: (x=7.2, y=11.9, z=120.2)
- you shift the point so that x and y are between 0 and 0.999... In 
    this case x=7.2-floor(7.2)=0.2, and y=0.9
- then you create a 4x4 spline (see lines 22-24 in sBA.m)
- you multiply the spline by z and normalize (line 26: fi= w*z(c)/W2).
- Now you add the spline fi to your plane. Because your x=7.2, you 
    have to add spline indexes 1,2,3,4 to the plane matrix at 
    indexes 6,7,8,9. Then your point x=7.2 comes in the center 
    of the spline, between 2nd and 3rd points. 
    (line 27: delta(I,J) = delta(I,J) + w.^2.*fi;)
- Now, if you have two points that are very close (x1=7.2, x2=8.1), 
    then the plane matrix ('delta' in the Matlab code) will have 
    double value around point x=8, because two similar splines were 
    summed. That is why you keep track of base-spline accumulation 
    (line 28: w2=...) and at the end these two are divided 
    (line 48: % f = delta./w2, for w2~=0). At points where you 
    have no points (i.e. w(i,j)==0), avoid division by zero.