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()