docs: 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