geomdl approximate_surface not working as expected

Question

I am trying to approximate a image with a b-spline surface. The geomdl library has a function called "approximate_surface" to do this. When I use this with an image generated from a nice smooth functions such as f(x)=x^2, things work as expected. But when I use my actual image the control points seem to be all over the place. The b-spline surface is nothing like the input surface. Why is this?

I made an example using a small image:

img = np.array([
 [   0. ,    0.  , 590.2 ,   0.   ,  0.   ,  0. ],
 [   0. , 1115.8 , 1704. , 1724.2 ,   0.  ,   0. ],
 [   0. , 1390.6 , 1810. , 1877.8 ,1309.4 ,   0. ],
 [   0. , 1274.  , 1735. , 1720.  ,1249.8 ,   0. ],
 [   0. , 1000.  , 1345. , 1345.  ,1168.  , 348.2],
 [   0. ,  734.  , 970.  ,1105.   ,895.8  ,  0. ],
 [   0. ,  627.  , 805.  , 928.8  ,  0.   ,  0. ],
 [   0. ,  394.  , 635.2 , 627.2  ,  0.   ,  0. ]])
plt.imshow(img)
plt.show()

pts = []
h, w = img.shape
for i in range(h):
  for j in range(w):
    pts.append((float(i), float(j), float(img[i,j])))

surf = fitting.approximate_surface(pts,size_u=h,size_v=w,degree_u=1, degree_v=1)
ctrlpts = np.array(surf.ctrlpts)
print(ctrlpts[0:10])

I get the output

[[ 0.00000000e+00  0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00 -3.97922748e-01 -5.32705554e+02]
 [ 0.00000000e+00  2.72312015e+00  1.09167840e+03]
 [ 0.00000000e+00  4.07453113e+00 -3.98262686e+02]
 [ 0.00000000e+00  5.00000000e+00  0.00000000e+00]
 [ 5.97819029e+00  0.00000000e+00  0.00000000e+00]
 [ 5.97819029e+00 -3.97922748e-01 -1.16343388e+04]
 [ 5.97819029e+00  2.72312015e+00  1.93970077e+04]
 [ 5.97819029e+00  4.07453113e+00 -2.33379745e+04]
 [ 5.97819029e+00  5.00000000e+00  2.91199984e+03]]

plotted out:

vis_obj = vis.VisSurface(vis.VisConfig())
surf.vis = vis_obj
surf.render()

You can try out the code here in colab.

EDIT:

I ended up using scipy "scipy.interpolate.bisplrep" does behave as I expected.

The maintainers from geomdl don't seem to respond to the current open github issues. So if you have the same problem, my advice is to use the scipy functies.