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).
 R. L. Hardy, “Multiquadric equations of topography and other irregular surfaces,” J. Geophys. Res., vol. 76, no. 8, pp. 1905–1915, Mar. 1971.  ```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+i*(domU-domU)
v=domV+j*(domV-domV)
(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)