Implement `Frame.psfs` with functional form

[pmelchior]

In 72b27ce, we introduced a proper point source fitter, which uses an functional form of the model PSF to determine the center of the source. Currently, one needs to specify that functional form to coincide with the pixelated image we stored in Frame.psfs. This is potentially dangerous and decidedly not elegant. So I suggest that Frame.psfs is either a callback function or an image. For the model, it would almost certainly be the former, for observations it could conceivably be either.

I'm not sure how to align this with the current Fourier class, though.

[fred3m]

It is true that Fourier (and FFTs in numpy in general) require a PSF image, not a function. So at the very least to create the difference kernels we'll need to generate a PSF image from the callback function for both a model and an observation, which should be fine.

One way to do this is with a PSF class, since there is some functionality this change will require. If initializing a frame with a PSF function, it will have to know the shape of the PSF and the pixel grid at which to sample from. For example, to match the HSC PSF size the shape would be roughly (41, 41) pixels, and by default we should assume that the model is centered on the center pixel, so the axes would be X=np.linspace(-20, 20, 41) and Y=np.linspace(-20, 20, 41). So in this design the PSF would be something like:

class PSF:
    def __init__(self, func=None, shape=None, axes=None, image=None):
        if image is None:
            if func is None or (shape is None and axes is None):
                msg = "A PSF must be initialized with either an `image` or `func` and `shape`"
                raise ValueError(msg)
            if len(shape) == 3:
                B = shape[0]
                _shape = shape[1:]
            else:
                B = None
                _shape = shape
            if axes is None:
                radii = (np.asarray(_shape) - 1) // 2
                axes = np.array([np.linspace(-r, r, 2*r+1) for r in radii])

            X, Y = np.meshgrid(axes[1], axes[0])
            if B is None:
                image = func(Y, X)
            else:
                image = np.array([func(b, Y, X) for b in range(B)])
        self.func = func
        self.image = image

    def render(Y, X):
        if self.func is not None:
            return self.func(Y, X)
        return self.image

Then we just call PSF.image to get the PSF image (for building the diff kernels) and PSF.render(Y, X) to generate a PSF image at the grid of coordinates (Y, X) (for modeling the point sources).

I just came up with this quickly off the top of my head, but it might be a good starting point.

[pmelchior]

Yeah, something like that. We can postpone this until the (painful) merge of adam with master, but I could use some help then.