def fit_polynomial(pixel_data, mask, clip=True):
    '''Return an "image" which is a polynomial fit to the pixel data
    
    Fit the image to the polynomial Ax**2+By**2+Cxy+Dx+Ey+F
    
    pixel_data - a two-dimensional numpy array to be fitted
    
    mask - a mask of pixels whose intensities should be considered in the
           least squares fit
           
    clip - if True, clip the output array so that pixels less than zero
           in the fitted image are zero and pixels that are greater than
           one are one.
    '''
    mask = np.logical_and(mask,pixel_data > 0)
    if not np.any(mask):
        return pixel_data
    x,y = np.mgrid[0:pixel_data.shape[0],0:pixel_data.shape[1]]
    x2 = x*x
    y2 = y*y
    xy = x*y
    o  = np.ones(pixel_data.shape)
    a = np.array([x[mask],y[mask],x2[mask],y2[mask],xy[mask],o[mask]])
    coeffs = scipy.linalg.lstsq(a.transpose(),pixel_data[mask])[0]
    output_pixels = np.sum([coeff * index for coeff, index in
                            zip(coeffs, [x,y,x2,y2,xy,o])],0)
    if clip:
        output_pixels[output_pixels > 1] = 1
        output_pixels[output_pixels < 0] = 0
    return output_pixels