ChromaticOpticalPSF ?

[rmandelb]

Hi all- For #590, I have been thinking of how to represent optical PSFs in a chromatic context. After a brief e-mail discussion with @jmeyers314 , I would like to propose some new functionality that I think will be useful for anyone who wants to represent chromatic optical PSFs (not particularly tied to WFIRST in any way). Note for Josh: this has evolved slightly from what we discussed because of my realization that I was thinking of the aberrations issue in an overly-complex way; it's actually quite simple.

Here is the basic problem. The chromatic effects come into optical PSFs in two ways:

• The obvious way is through the diffraction limit depending on wavelength/diameter.
• The less obvious way is that we require aberrations in units of wavelengths in order to determine the wavefront, but aberrations are typically defined in terms of some physical length, so for a different wavelength, this becomes a different number required by the OpticalPSF.

So the idea would be to have something like this (schematically; I'm sure it needs work on the details to actually be functional):

class ChromaticOpticalPSF(ChromaticObject):
    def __init__(self, diam, aberrations, **kwargs):
        self.diam = diam # in m
        self.aberrations = aberrations # in physical units like microns
        self.kwargs = kwargs
    def evaluateAtWavelength(self, wave):
        lam_over_diam = wave / self.diam # and do what's needed to get the right units
        aberrations = self.aberrations / wave # and do what's needed to get the right units
        ret = OpticalPSF(lam_over_diam=lam_over_diam, aberrations=aberrations, **self.kwargs)
        return ret

There are potential timing issues since OpticalPSF can be slow. A few options that Josh suggested:

• Have users thin the bandpasses and SEDs to reflect how many points at which they want to evaluate. This puts the onus on the user.
• Precompute the PSF at a few points in wavelength as part of the constructor, and then have evaluateAtWavelength interpolate between those PSF images. This would require the constructor to have options for the user to control the pre-computation (number of grid points etc.). For a given combination of telescope diameter + obscuration + struts + aberrations, this could be reused for different SEDs/bandpasses, since a change would just require reweighting the precomputed quantities appropriately rather than precomputing anything else (as long as the wavelength range of interest hasn't changed).

This is something I was thinking of doing pretty much immediately, to enable some of what I want to do on #590, but comments / questions / suggestions are very welcome. My tentative plan is to try to get something that works first using something like the above code, which will be slow, and then consider converting to a pre-computation scheme afterwards.

[jmeyers314]

Hi @rmandelb,

I was thinking about this some more and have one other idea to bring up; I'm not sure if I like it or not at this point. We could make an even more generic chromaticizer that would use as input any function of wavelength that returns a GSObject. For example:

def chromaticOpticalPSF(wave, *args, **kwargs):
    diam = args[0]
    lam_over_diam = wave / diam
    aberrations = args[1] / wave
    ret = OpticalPSF(lam_over_diam=lam_over_diam, aberrations=aberrations, **kwargs)

class Chromaticizer(ChromaticObject):
    def __init__(f, *args, **kwargs):
        self.args=args
        self.kwargs=kwargs
        self.f=f
    def evaluateAtWavelength(self, wave):
        return self.f(wave, *self.args, **self.kwargs)

diam = 8.4
aberrations=[0]*12
psf = Chromaticizer(chromaticOpticalPSF, diam, aberrations, **other_kwargs)

I think that ChromaticAtmosphere could also take advantage of such a class, which could still have optional arguments for interpolation/gridpoints. This would be a pretty easy way for users to create their own chromatic surface brightness profiles too.

[rmandelb]

I kind of like the chromaticizer idea. @rmjarvis , what do you think?

[mdschneider]

Hi @rmandelb and @jmeyers314,

It appears that both solutions for including the proper wavelength dependence in the optical PSF rely on first specifying physical optics perturbations and then converting to wavefront perturbations in units of waves.

Some of this functionality already exists as a GalSim wrapper in the Dark Energy Science Collaboration repository (not visible to non-members for now unfortunately). However, the linked wrapper script includes an extra feature in the self-consistent prediction of aberrated pupil wavefronts for any field position (and wavelength) given the same physical optics perturbations. If consistent prediction of PSFs at multiple image plane positions is not a priority then this script is more complicated than what you need I suppose.

The main function of the wrapper script is to predict Zernike mode coefficients to be input into galsim.optics.wavefront for a supplied physical optics perturbation (given as bending modes on each mirror and tilts/decenters of each optical element), wavelength, and normalized image plane coordinate.

I have created separate classes ZernikeCoeffsBendingModes and ZernikeCoeffsBodyMotions to predict the contributions to the exit pupil wavefront error from each class of optics perturbations. Instances of these classes are created inside an optics_psf function that then calls the related galsim.optics.psf_image routine.

Would it be helpful to incorporate this functionality into GalSim at this point to address your concerns above? (I think @rmandelb you even suggested I do this earlier this year - sorry I didn't appreciate how useful this might be until recently!) Necessary steps would include:

1. Code cleanup to conform to GalSim standards.
2. Validation tests for the telescope model(s) that would be provided.
3. Potentially, higher order Zernike terms in the wavefront expansion.

[rmjarvis]

I kind of like the chromaticizer idea. @rmjarvis , what do you think?

I'm not sure I see the point. What is the advantage of this over just defining a ChromaticOpticalPSF class? It looks like the same code, just with a clunkier calling syntax.

[rmandelb]

For just this, it’s unnecessary. But as soon as we want any other chromatic classes of non-trivial complexity, it seems like it will reduce the amount of necessary code.

[rmjarvis]

I don't see why it would be easier to write a function to be passed to the Chromatizer class than it is to write the class in the first place. Already, deriving from ChromaticObject means that the only thing you really need to define is evaluateAtWavelength, which is identical to the function you would pass to Chromatizer, so that doesn't save anything.

If you are thinking about the possible implementation details about precomputing the function at select wavelengths and then interpolating, then I could see some possibilities for code reuse. But in that case, I would recommend instead adding an additional layer in the class tree:

class InterpolatedChromaticObject(ChromaticObject):
    def __init__(self, waves):
        self.waves = waves
        self.objs = [ self.simpleEvaluateAtWavelength(wave) for wave in waves ]

    def evaluateAtWavelength(self, wave):
        # Some clever interpolation procedure that is presumably faster for typical use cases.

class ChromaticOpticalPSF(InterpolatedChromaticObject):
    def __init__(self, diam, aberrations, **kwargs):
        # Figure out some appropriate initial wavelengths at which to compute the profile initially.
        waves = ...
        super(ChromaticOpticalPSF, self).__init__(waves)
        self.diam = diam
        self.aberrations = aberrations
        self.kwargs = kwargs

    def simpleEvaluateAtWavelength(self, wave):
        # Equivalent to the old evaluateAtWavelength function.

This would then let other classes leverage the functionality you develop for the interpolation in wavelength, but in a more typically object-oriented way than the Chromatizer class.

[jmeyers314]

Yeah, I was trying to think about future chromatic objects too. Especially if there's more than one case where we want to precompute at specific wavelengths and interpolate, I think this could lead to less code duplication. I agree it's a bit clunky though. Maybe precompute/interpolation could be built directly into ChromaticObject and turned on with an optional keyword. I'm open to suggestions.

[jmeyers314]

Two minutes too slow.

[rmjarvis]

;-)

[rmandelb]

OK, that's kind of nifty @rmjarvis .

@mdschneider : in this particular case I think that the code you're pointing to has extra functionality that we don't need for the application I have in mind (at least not now). I don't (yet) have a full telescope model that I want to include, just models at the center of each detector, so I just want to be able to say "given the aberrations in microns for the center of this detector, give me the chromatic PSF for this detector given a particular SED and passband". I do stand by my original statement that the code you describe is something that we might want to have in GalSim :) but I think I would like to build up the general ChromaticOpticalPSF functionality first given that it's (a) a more minimal representation of what I need and (b) useful for other applications like a simple definition of ChromaticAtmosphericPSFs.

Do you think that, given a ChromaticOpticalPSF implementation, the process of porting over your more general code might be simplified? (i.e., you could basically turn a bunch of your code into a call to ChromaticOpticalPSF, so you would primarily be adding the bits of your code that handle the telescope model?)

[jmeyers314]

It would be nice to be able to quickly cycle the interpolation on/off. Perhaps setting waves == None in the InterpolatedChromaticObject constructor should yield a non-interpolating InterpolatedChromaticObject.

[rmjarvis]

Agreed. That makes sense as a feature to include.

[rmandelb]

Yes. I presume you mean that in that case, it would go to the usual rule of using the SED and bandpass to choose the values of wavelength to be used in evaluating the profile. I think this is a good idea. For example, it would allow people to very easily do some validation of how inaccurate their interpolation scheme is (using just a few wavelengths) for a limited set of cases, by comparing the interpolated vs. non-interpolated results.

[jmeyers314]

Exactly.

Actually, I'm rather optimistic that relatively few wavelengths may be sufficient. I think for most cases, PSF(lambda) varies more smoothly than SED(lambda). With the above approach, one could sample the PSF coarsely, and the SED more finely.

As for clever interpolation schemes, I was thinking of something not clever at all. Namely, drawing each PSF(wave) into an InterpolatedImage (being careful to draw each with the same size and scale), and then just linearly interpolating pixel-by-pixel. How does that sound?

[rmjarvis]

Sounds like a reasonable thing to try. But I have no intuition about how well it will work.

[rmandelb]

Yes, I was thinking of the same not-clever interpolation scheme. Worth a try! If I can succeed at getting it to work how you suggested, so interpolation is easy to switch off and on, then it should be easy to test.

[rmandelb]

Hi all -

I have a test script on machine that implements everything we discussed here (the InterpolatedChromaticObject class, of which ChromaticOpticalPSF is one example).

The good news is that it's possible to implement something that is basically right (i.e., gives the right PSF) without significant changes of the simple implementation that we discussed. The bad news is that right now calling drawImage using the interpolated version is actually slower than using the non-interpolated version, which is of course extremely silly.

Here's the reason: the InterpolatedChromaticObject constructor evaluates and draws the image at a few selected wavelengths (with the same image size and scale). Then the evaluateAtWavelength routine linearly interpolates between these images, which is fast, and instantiates an InterpolatedImage (so it can return a GSObject as is currently required). The InterpolatedImage instantiation is really slow because of the step-k and max-k calculations. If I try to shut those calculations off using calculate_maxk and calculate_stepk=False, then badness ensues; its initial guess is so conservative as to require enormous FFTs. We do have an option to pass in maximum values for these parameters, but if that option is taken then it still does the full calculation, which is slow. So my conclusion is that I need to do either of the following:

1) Modify InterpolatedImage() so that you can force it to use a particular maxk and stepk value (caveat user). I know reasonable values for these parameters from the original set of objects that were instantiated by the constructor for the InterpolatedChromaticObject, so I can simply force it to use those.

2) Do the interpolation between wavelength values in some other way than linear interpolation of images, so that we don't require InterpolatedImage instantiations.

3) Change drawImage for InterpolatedChromaticObject. Right now, what it does is (1) use evaluateAtWavelength to get a GSObject corresponding to each wavelength, then (2) draw them into images, then (3) do the integral using those images. I think that for InterpolatedChromaticObject (when actually interpolating) we could change it to (1) use evaluateAtWavelength to get an image output (based on the fast interpolation scheme), (2) integrate over those, (3) make an InterpolatedImage representing that integral, and (4) draw the results into a final image with the desired size and scale.

Thoughts? (1) would be easy for me to do, whereas (3) is harder, but I'm waffling over which of them is actually right. I mean, (1) is so easy as to feel like a cheat, and I feel that it is fundamentally a bit silly that the process of drawing an image of an InterpolatedChromaticObject should require going back and forth between the original set of images, many InterpolatedImages, and many images that are drawn for them (whereas if I do (3) it'll be an original set of images, many quickly-generated interpolated images, and a single InterpolatedImage from which the final image gets drawn).

I don't have any smart ideas about (2). And the fact that I'm having this problem at all makes me wonder if I'm doing something fundamentally wrong.

[rmjarvis]

I don't have a problem with (1), although I think the forced stepk and maxk values should be private construction parameters (e.g. _force_stepk and _force_maxk) with no @param entries in the doc string. This way it discourages non-expert users from using them inappropriately.

But (3) does sound like a better way to implement drawImage, to be honest, so I think that's probably the best solution.

But I'll just point out that even if you go with (3), you should probably implement (1) as well so the evaluateAtWavelength method behaves sensibly. Users may want that function for other uses than within drawImage. So a good way forward would be to do (1), leaving drawImage as is for now. Then try to do (3), but you'll have a fallback position that may be acceptable if that seems too hard to get working.

[rmandelb]

Thanks for the quick feedback, Mike.

Regarding (1), I completely agree about making these private construction parameters.

Regarding (3), I was reluctantly thinking it's the better way to go, but more annoying. It also seems to require changes to the ImageIntegrator class as well (at least, the __call__ method would have to be changed if evaluateAtWavelength is returning an image instead of a GSObject). One more question about this: most of drawImage and ImageIntegrator.__call__ should be the same for the InterpolatedChromaticObject as for ChromaticObject, so I was thinking of having the ChromaticObject.drawImage() method do different things depending on whether the object is an InterpolatedChromaticObject or not (i.e., check its own identity, and have if statements to switch between code blocks for the parts of the routine that need to be changed)?

I think that you're right - I should try (1) first, then (3) later. (1) should at least let me get something that works in finite time so I can do some basic testing before doing the more significant surgery needed to make (3) work.

[rmjarvis]

I was thinking of having the ChromaticObject.drawImage() method do different things depending on whether the object is an InterpolatedChromaticObject or not (i.e., check its own identity, and have if statements to switch between code blocks for the parts of the routine that need to be changed)?

I don't like this design model. A better design, in my opinion, is to pull out common functional steps and make them helper functions, which the InterpolatedChromaticObject.drawImage() method can use in its implementation. e.g. wave_list = self._getCombinedWavelist(bandpass).

I'm not sure what else in the original is salvageable though. It looks like most of the function will need to be significantly modified to deal with the direct images instead of profiles. So it seems like it would be easier to just have a fresh implementation in InterpolatedChromaticObject.

Similarly, I think you want to make a new DirectImageIntegrator class that looks a lot like ImageIntegrator but skips the drawImage calls and just gets the images directly from whatever function does that. Should be quite simple really, and clearer than having if statements within ImageIntegrator.__call__.

[rmandelb]

Okay, in the meantime I implemented the stepk and maxk forcing routines. These reduced the time needed for the InterpolatedChromaticImage in my test case by a factor of ~10. However, it's still much slower than I would like. I determined that the slowness is directly caused by all the back-and-forthing we are doing between InterpolatedImages and images. Right now, since we interpolate images and instantiate InterpolatedImages which then get drawn into images when integrating over a bandpass, it's really super inefficient, and the refactoring we discussed seems necessary to me.

Regarding the question of what else in the original is salvageable: Perhaps I'm thinking about this incorrectly, but I thought that quite a bit of it is. I agree that I need to make a DirectImageIntegrator class because it's the integration over the images that is really the game-changer. But going through what ChromaticObject.drawImage() does, I can break it down into these steps:

1) Set up an output image -> this will still need to be done. 2) Make combined wave list -> this will still need to be done. 3) Do the simple steps that are necessary to draw separable profiles -> we won't need this for ChromaticInterpolatedImages, which aren't separable. 4) Decide on an integrator and potentially do further operations to wave list depending on which is chosen -> still needs to be done. 5) Actually do the call to the integrator -> still needs to be done, this time to the new DirectImageIntegrator. \ New step needed here for ChromaticInterpolatedImage: Need to instantiate an InterpolatedImage from the output of the integrator, and draw it into the output image that was set up as step (1).

I think that for ChromaticInterpolatedImage, steps (1), (2), and (4) can be reused pretty much as-is, step (5) needs some minor changes to call this new class, and we have a new step in the process to add afterwards. I'm totally fine with making helper functions for (1), (2), and (4), to be used by both ChromaticObject.drawImage() and InterpolatedChromaticObject.drawImage() as needed, but since you seemed to be implying that we can't have much reuse between the two, I wanted to check whether we are thinking about this function in the same way.

As a final point, one reason it's slow now is that if we're drawing, say, a convolution of the ChromaticOpticalPSF with a pixel response, then the image integrator is at each wavelength drawing the convolved profile, meaning it has to use fourierDraw in the C++ layer and do the FFT each time. I was thinking that for InterpolatedChromaticObject, we could get around this by averaging the images over wavelength, making the InterpolatedImage, and only convolving the final wavelength-averaged InterpolatedImage with the pixel response, which will be faster. (Assuming the user is trying to do drawImage while including the pixel response, instead of using method no_pixel.) I think this is fair to do since the pixel response is wavelength-independent.

[rmjarvis]

since you seemed to be implying that we can't have much reuse between the two, I wanted to check whether we are thinking about this function in the same way.

1 might be the same, but it's only 2 lines of code, so not so hard to just duplicate. 2 is the same, which is why I suggested that as a helper function. 3 is gone. 4 is different; I'm assuming you will be forcing the use of DirectImageIntegrator, so the existing lines here are irrelevant to the new functionality, although you might have new similar lines. 5 is different, since you don't pass it the evaluateAtWavelength function anymore. And you can probably just call the DirectImageIntegrator directly rather than use integrator.

So the steps are similar, but it seems simpler to just copy the current version into the derived class and modify it appropriately, rather than add a bunch of if statements in the existing version.

[rmandelb]

Oh, I see. I had been thinking that evaluateAtWavelength could change to return an image rather than an InterpolatedImage, whereas you had been thinking it would stay the same and the DirectImageIntegrator would not use evaluateAtWavelength. I think your way is safer (it's just not what I was thinking before). Thanks.

[rmjarvis]

Yes, don't change what evaluateAtWavelength means. That's a documented function that needs to keep its documented functionality. The user may want to use it for something other than drawImage, so we don't want to make it do something different for these classes.

[rmandelb]

Working on this now, and I'm having an issue with multiple inheritance that I couldn't sort out very satisfactorily (google kind of failed me on this one)..

My original idea was to have DirectImageIntegrator follow a similar pattern as ImageIntegrator, meaning that it has a __call__ method but never gets instantiated directly (and of course, __call__ looks different because of what it's trying to do). Then SampleIntegrator and ContinuousIntegrator become subclasses of both ImageIntegrator AND DirectImageIntegrator. But the problem I have then is that I need a way to control which flavor of SampleIntegrator i'm getting, i.e., which __call__ method should it inherit?

I know there are rules for which parent class it would look at in which order, but I can't figure out a way (for example using super) to control it directly. Do you have any idea?

[rmandelb]

I guess for now I could make DirectImageIntegrator a class that works like SampleIntegrator, evaluating only at values of wavelength defined by the bandpass, rather than having any class hierarchy associated with it at all.

[rmjarvis]

I'm having an issue with multiple inheritance [...] Do you have any idea?

Yes. Don't use multiple inheritance. :) It's very rarely the right solution, and I don't think it's the right thing to do here.

In particular, I think DirectImageIntegrator is really a different thing than ImageIntegrator. It happens to be implemented similarly, but it's conceptually doing something different. So it wouldn't make sense for something to derive from both classes. Inheritance should imply an "is a" relationship, and I don't think you want to make SampleIntegrator be both an ImageIntegrator and a DirectImageIntegrator.

I could make DirectImageIntegrator a class that works like SampleIntegrator

I think this is the best way to go. If you want, you could preserve the class hierarchy and have DirectSampleIntegrator inherit from DirectImageIntegrator in case anyone ever cares to have a different kind of DirectImageIntegrator. But at least at the moment, the sample algorithm is the only one you actually want, so I think that's all you'll need.

And if there are reusable code bits that should be in both SampleIntegrator and DirectSampleIntegrator, feel free to pull them out as functions that both versions can call.

[rmandelb]

Okay. This works now, and isn't completely embarrassing, in the sense that using the "clever" way of interpolating over wavelength is now faster than doing it the brute force way, and both return sensible-looking PSF images.

But I started thinking about this in the context of ChromaticConvolution - i.e., convolving with a galaxy or anything with a non-trivial SED - and this is where I got lost. It seems to me that here we cannot avoid the instantiation of a GSObject at each wavelength that is being sampled, so there is nothing clever to do, and any ChromaticConvolution that has an InterpolatedChromaticObject as one of the elements being convolved has to do the brute force evaluation at each wavelength. Does that sound right to you?

[rmjarvis]

Yes, I think when convolving by something else chromatic, you'll have to go through the normal process of using evaluateAtWavelength.

[rmandelb]

... and despite what I said, it does seem to be faster. It seems that this might be because instantiating the InterpolatedImage at each wavelength via evaluateAtWavelength is still somewhat expedited (compared to OpticalPSF creation) due to the use of _force_stepk and _force_maxk. If I run this test script:

import galsim
import numpy as np
import time

test_aberrations = np.zeros(12)
# Let's give it some defocus
test_aberrations[4] = 200. # nm for now

t1 = time.time()
obj1 = galsim.ChromaticOpticalPSF(2.4, test_aberrations, 100., 1000., 10)
t2 = time.time()
print 'Time to instantiate interpolated ChromaticOpticalPSF: ',t2-t1

t1 = time.time()
obj2 = galsim.ChromaticOpticalPSF(2.4, test_aberrations, 100., 1000., None)
t2 = time.time()
print 'Time to instantiate fully-evaluated ChromaticOpticalPSF: ',t2-t1

test_bandpass = galsim.Bandpass('examples/data/LSST_r.dat')
test_bandpass = test_bandpass.thin(rel_err=5.e-4)
print 'Number of wavelengths in bandpass: ',len(test_bandpass.wave_list)

t1 = time.time()
im_obj1 = obj1.drawImage(test_bandpass, scale = 0.01, method='no_pixel')
t2 = time.time()
print 'Time to draw interpolated ChromaticOpticalPSF with no_pixel: ',t2-t1

t1 = time.time()
im_obj2 = obj2.drawImage(test_bandpass, scale = 0.01, method='no_pixel')
t2 = time.time()
print 'Time to draw fully-evaluated ChromaticOpticalPSF with no_pixel: ',t2-t1

t1 = time.time()
im_obj1 = obj1.drawImage(test_bandpass, scale = 0.01)
t2 = time.time()
print 'Time to draw interpolated ChromaticOpticalPSF with pixel: ',t2-t1

t1 = time.time()
im_obj2 = obj2.drawImage(test_bandpass, scale = 0.01)
t2 = time.time()
print 'Time to draw fully-evaluated ChromaticOpticalPSF with pixel: ',t2-t1

im_obj1.write('psf_obj1.fits')
im_obj2.write('psf_obj2.fits')
mom_obj1 = im_obj1.FindAdaptiveMom()
mom_obj2 = im_obj2.FindAdaptiveMom()
#print 'Size and flux comparison: ',mom_obj1.moments_sigma, mom_obj2.moments_sigma, \
#    im_obj1.array.sum(), im_obj2.array.sum()

print 'Setting up simple galaxy objects'
sed = galsim.SED('examples/data/CWW_Scd_ext.sed', wave_type='Ang')
simple_gal = galsim.Exponential(half_light_radius=0.25)
simple_gal = simple_gal.shear(g1=0.3, g2=-0.1)
simple_gal *= sed

obj1_obj = galsim.Convolve(simple_gal, obj1)
obj2_obj = galsim.Convolve(simple_gal, obj2)

t1 = time.time()
im_obj2 = obj2_obj.drawImage(test_bandpass, scale=0.05)
t2 = time.time()
print 'Time to draw galaxy convolved with fully-evaluated ChromaticOpticalPSF: ',t2-t1
im_obj2.write('im_obj2.fits')

t1 = time.time()
im_obj1 = obj1_obj.drawImage(test_bandpass, scale=0.05)
t2 = time.time()
print 'Time to draw galaxy convolved with interpolated ChromaticOpticalPSF: ',t2-t1
im_obj1.write('im_obj1.fits')

then the output on my system is

Time to instantiate interpolated ChromaticOpticalPSF:  2.35137701035
Time to instantiate fully-evaluated ChromaticOpticalPSF:  3.81469726562e-05
Number of wavelengths in bandpass:  23
Time to draw interpolated ChromaticOpticalPSF with no_pixel:  0.716335058212
Time to draw fully-evaluated ChromaticOpticalPSF with no_pixel:  2.12347507477
Time to draw interpolated ChromaticOpticalPSF with pixel:  2.18649196625
Time to draw fully-evaluated ChromaticOpticalPSF with pixel:  2.86208796501
Setting up simple galaxy objects
Time to draw galaxy convolved with fully-evaluated ChromaticOpticalPSF:  3.32388210297
Time to draw galaxy convolved with interpolated ChromaticOpticalPSF:  1.30094790459

I haven't done rigorous accuracy tests, but the images look pretty good so far. in terms of timing, this particular instance gives roughly equal time for the fully-evaluated ChromaticOpticalPSF vs. the interpolated one (once we include the initialization + drawing time), but if we're drawing many instances of this PSF convolved with some galaxies then we'll get this large savings in the drawing time for each instance.

Still, I think that for mass production we will need further optimization - 1.3 seconds is not particularly fast!

[rmandelb]

Closing #618.