Fit implementation

Question

Hello I am trying to fit a psf to an image. The background should be aproximated by lower order polynomials. If I just take a constant it works fine:

def fitter(image, trueImage,psf,intensity): 
   p0 = [intensity]            
   p0.append(np.amin(trueImage)*4**2) 
   meritFun = lambda p: np.ravel(image-(p[0]*psf+p[1]))
   p = least_squares(meritFun,p0,method='trf')

Now I have the issue of how to define the x and y's for my polynomials:

#Does not work!
def fitter(image, trueImage,psf,intensity): 
   p0 = [intensity]            
   p0.append(np.amin(trueImage)*4**2) 
   p0.append(1)
   p0.append(1) #some clever initial guess
   meritFun = lambda p: np.ravel(image-(p[0]*psf+p[1]+p[2]*x+p[3]*y))
   p = least_squares(meritFun,p0,method='trf')

x and y are obviuosly the indices i, j of my image array, but how do I tell that my fit routine?

Answer 1

As Paul Panzer mentioned in the comments, one way of solving this is by using np.ogrid:

def fitter(image, trueImage,psf,intensity): 
    x = np.ogrid[:image.shape[0]]
    y = np.ogrid[:image.shape[1]]
    p0 = [intensity]            
    p0.append(np.amin(trueImage)*4**2)
    p0.append(0)
    p0.append(0)
    meritFun = lambda p: np.ravel(image-(p[0]*psf+p[1]+p[2]*x+p[3]*y))
    p = least_squares(meritFun,p0,method='trf')