docs: def splinefit2d(x, y, z, px, py, mask=None):
      '''Make a least squares fit of the spline (px,py,pz) to the surface (x,y,z).
  If mask is given, only masked points are used for the regression.'''
      
      if mask is None:
          V = np.array(spline_matrix2d(x, y, px, py))
          a = np.array(z.T.flatten())
          pz = np.linalg.lstsq(V.T, a.T)[0].T
      else:
          V = np.array(spline_matrix2d(x,y,px,py,mask))
          indices = np.nonzero(np.array(mask).T.flatten())
          if len(indices[0])==0:
              pz = np.zeros((len(py),len(px)))
          #indices is empty when mask changes to all zeros
          else:
              a = np.array((z.T.flatten()[indices[0]]))
              pz = np.linalg.lstsq(V.T, a.T)[0].T
              
      pz=pz.reshape((len(py),len(px)))
      return pz.transpose()