Surface approximation of point clouds by using multiquadric functions

The following literature provides the easiest way to generate a surface that approximates a given point set.
Only the limitation found in this approach is that the generated surface becomes a height field, i.e. it is always of the form z=f(x,y).
[1] R. L. Hardy, “Multiquadric equations of topography and other irregular surfaces,” J. Geophys. Res., vol. 76, no. 8, pp. 1905–1915, Mar. 1971.
1.png
2.png

import rhinoscriptsyntax as rs
import Rhino.Geometry as rg
import math
func=lambda xi,xj,yi,yj:math.sqrt(((xj-xi)**2)+((yj-yi)**2)+x)
A=rg.Matrix(len(pnts),len(pnts))
z=rg.Matrix(len(pnts),1)
for i in range(len(pnts)):
    for j in range(len(pnts)):
        pi=pnts[i]
        pj=pnts[j]
        A[i,j]=func(pi.X,pj.X,pi.Y,pj.Y)
        z[i,0]=pi.Z
A.Invert(0.0)  #this parameter should be 0.0
c=A*z
a=list()
for srf in srfs:
    domU=srf.Domain(0)
    domV=srf.Domain(1)
    for i in rs.frange(0,1,0.05):
        for j in rs.frange(0,1,0.05):
            u=domU[0]+i*(domU[1]-domU[0])
            v=domV[0]+j*(domV[1]-domV[0])
            (b,P,A)=srf.Evaluate(u,v,0)
            z=0
            for j in range(len(pnts)):
                z=z+c[j,0]*func(P.X,pnts[j].X,P.Y,pnts[j].Y)
            newP=rg.Point3d(P.X,P.Y,z)
            a.append(newP)
newSrf=rs.AddSrfPtGrid((20,20),a)
newSrf = rg.NurbsSurface.CreateFromPoints(a, 20, 20, 3, 3)
a=list()
a.append(newSrf)