Help request: interpreting PSF

Sorry for putting a help request here. If there's a better way to do this, let me know. I'll be brief, hopefully not too brief.

I'm exploring POPPY for use modeling a star camera. The model camera lens has a diameter of 25 mm, focal length 50 mm, 20 deg FOV, with a focal plane array having pixel pitch 13 microns, 1024 x 1024. This results in a pixel scale of about 70 arcsec/pixel and the focal plane array/camera produces an image that's roughly 70 arcsec/pixel, with 1024 pixels across.

I set up a simple optical system,

lmda = 550e-9  # wavelength, meters
rad = 12.5e-3  # lens radius
pxscale = 70. # detector arcsec/pixel
fov_pixels = 1024 # full-angle FOV ??

osys = poppy.OpticalSystem()
osys.addPupil( poppy.CircularAperture(radius = rad))
osys.addDetector(pixelscale = pxscale,
             fov_arcsec = fov_pixels*pxscale)

and produce a psf

psf = osys.calcPSF(lmda)
psf_ax, colorbar = poppy.display_PSF(psf,
                                 return_ax=True,
                                 title = 'The Airy Function',
                                 scale = u'log')

and extract the data

pdata = psf[0].data

I expected a 1024 x 1024 array with a psf in the center. Instead I get a 2048 x 2048 array containing many peaks: 23 x 23 of them. Evidently I don't understand something essential about how poppy samples and does Fourier Transforms. Of course, I'm interested in the details of the psf, not the complete image, but I think I need to understand what poppy is doing. I have many questions, but perhaps just a pointer in the right direction will be all that I need. Thanks.

@mperrin can probably explain better, but I believe what you're seeing is a consequence of the implicit repetition of the pupil array in the transform step... In practice this isn't usually an issue, because the flux in the PSF beyond a couple arcseconds is much less than one percent, so it's sufficient to compute a smaller array of pixels (and place them in a larger array of zeros as part of building up a simulated image)

There are a couple things going on here:

The 2048 vs 1024 pixels difference is simply explained by the default oversampling factor of 2x in calcPSF. This is a user-adjustible parameter so you can make it whatever you want. (and the Instrument subclass provides a more sophisticated interface to this including automatic downsampling to the detector resolution but that takes some more work to set up.)
Secondly, the repeated copies of the PSF are due to the replication that is assumed by the FT algorithm, as @josePhoenix pointed out. The maximum spatial frequency that can be sampled is set by the sampling in the input pupil plane. If you have an pupil array that is 100 pixels across, you can represent up to 50 cycles per pupil (at Nyquist sampling, 2 pixels per cycle). This means you can represent up to 50 lambda/D in the image plane, and after that you will get aliasing. See for instance the "Sampling Sinusoidal Functions" section of https://en.wikipedia.org/wiki/Aliasing. The software attempts to calculate the FT of a high spatial frequency, but due to limited sampling it returns instead a value appropriate for a much lower frequency and this results in repeated copies of the (low spatial frequency!) core of the PSF. The fix for this is to increase the sampling of the input pupil array. The case you are interested in may actually be better represented by the new FresnelOpticalSystem functionality in our recent release. That class provides a more convenient interface to adjust the sampling of the input wavefront, since you can just specify npix and pupil_diameter directly in the creation of the FresnelOpticalSystem. (This newer code now has a better API than our older code I think... Lessons learned!)

All that said, your proposed system is significantly undersampled spatially. For a 50 mm lens the theoretical diffraction limited PSF is just a few arcseconds at visible wavelengths. It's at least an order of magnitude better than your stated 70 arcsec/pixel sampling. So the kinds of fine diffractive structure that poppy is designed to model will not be easily seen in your instrument. Furthermore, remember that poppy is just modeling the effects of diffraction, based on the input parameters describing an optical system. In your case with a simple circular pupil that will just give an Airy function, without any need for numerical calculations of a more complex PSF. More to the point, your instrument is not going to have diffraction limited optical performance with only a single lens. It takes multiple powered surfaces to get near diffraction limited performance over an extended field, particularly for anything as large as 20 degrees across. Depending on your application, you're likely in a regime dominated by geometric optics rather than diffraction.

Your given parameters (50 mm focal length, 25 mm aperture diameter) are in the realm of commercial camera lenses, but even those don't give nearly diffraction limited performance over such a wide field. See for instance this Nikon 50 mm f/1.8 lens: http://imaging.nikon.com/lineup/lens/singlefocal/normal/af-s_50mmf_18g/ That page has a chart showing the modulation transfer function (MTF), which depicts how the MTF declines significantly as you move away from the optical axis. And that's for a seven lens system, one of which is an aspheric surface. A single lens is going to be totally dominated by low-order aberrations off axis. It's possible to model such effects in poppy, but that requires having an input model describing the wavefront error as a function of field angle, for instance as derived from the optical model of the lens.

THanks for you detailed response. I kept my original question short because I was afraid I was asking the question in the wrong place (developer's planning forum), so I eliminated many details. I'll be studying what you said more closely when the holiday is over.

I'm familiar with the DFT and issues of aliasing, sampling, periodicity, (windowing ...). One thing I couldn't figure out is what the sampling rate in the pupil is. Is pupil sampling at 100 pixels the default, or did I unknowingly set it somewhere? I take your words to mean that I did manage to correctly (modulo x2) the sampling/resolution of my focal plane array. I was (and still am) a bit curious about the oversampling parameter. I'll look more closely at that, too,

The single lens would be just for simplicity during my learning curve. I will be modeling lenses, so I will have some idea of aberrations. I've already built an element classes that introduce tilt and coma, and the ThinLens class provides defocus. The models will be used for testing centroiding algorithms. The final product is a star camera which will be defocused so that the intensity distribution covers at least a 3x3 array of pixels.

I have just noticed that the git repository has a much-updated version which includes Fresnel diffraction. THat might be more appropriate for my use case. I haven't yet installed it.

In the meantime, as a potential alternative, I'm implementing the Kirchoff-Fresnel diffraction integral, aberrated verstion, from Born and Wolf Sec 9.4. It's slow, but that might be ok for me. Poppy would be more convenient and faster.

I'm sorry but appreciative that my spartan question lead you to reply in such detail.

On Wed, Nov 25, 2015 at 1:21 PM, Marshall Perrin notifications@github.com wrote:

There are a couple things going on here:

The 2048 vs 1024 pixels difference is simply explained by the default oversampling factor of 2x in calcPSF. This is a user-adjustible parameter so you can make it whatever you want. (and the Instrument subclass provides a more sophisticated interface to this including automatic downsampling to the detector resolution but that takes some more work to set up.)

Secondly, the repeated copies of the PSF are due to the replication that is assumed by the FT algorithm, as @josePhoenix https://github.com/josePhoenix pointed out. The maximum spatial frequency that can be sampled is set by the sampling in the input pupil plane. If you have an pupil array that is 100 pixels across, you can represent up to 50 cycles per pupil (at Nyquist sampling, 2 pixels per cycle). This means you can represent up to 50 lambda/D in the image plane, and after that you will get aliasing. See for instance the "Sampling Sinusoidal Functions" section of https://en.wikipedia.org/wiki/Aliasing. The software attempts to calculate the FT of a high spatial frequency, but due to limited sampling it returns instead a value appropriate for a much lower frequency and this results in repeated copies of the (low spatial frequency!) core of the PSF. The fix for this is to increase t he sampling of the input pupil array. The case you are interested in may actually be better represented by the new FresnelOpticalSystem functionality in our recent release. That class provides a more convenient interface to adjust the sampling of the input wavefront, since you can just specify npix and pupil_diameter directly in the creation of the FresnelOpticalSystem. (This newer code now has a better API than our older code I think... Lessons learned!)

All that said, your proposed system is significantly undersampled spatially. For a 50 mm lens the theoretical diffraction limited PSF is just a few arcseconds at visible wavelengths. It's at least an order of magnitude better than your stated 70 arcsec/pixel sampling. So the kinds of fine diffractive structure that poppy is designed to model will not be easily seen in your instrument. Furthermore, remember that poppy is just modeling the effects of diffraction, based on the input parameters describing an optical system. In your case with a simple circular pupil that will just give an Airy function, without any need for numerical calculations of a more complex PSF. More to the point, your instrument is not going to have diffraction limited optical performance with only a single lens. It takes multiple powered surfaces to get near diffraction limited performance over an extended field, particularly for anything as large as 20 degrees across. Depending on your application, y ou're likely in a regime dominated by geometric optics rather than diffraction.

Your given parameters (50 mm focal length, 25 mm aperture diameter) are in the realm of commercial camera lenses, but even those don't give nearly diffraction limited performance over such a wide field. See for instance this Nikon 50 mm f/1.8 lens: http://imaging.nikon.com/lineup/lens/singlefocal/normal/af-s_50mmf_18g/ That page has a chart showing the modulation transfer function (MTF), which depicts how the MTF declines significantly as you move away from the optical axis. And that's for a seven lens system, one of which is an aspheric surface. A single lens is going to be totally dominated by low-order aberrations off axis. It's possible to model such effects in poppy, but that requires having an input model describing the wavefront error as a function of field angle, for instance as derived from the optical model of the lens.

— Reply to this email directly or view it on GitHub https://github.com/mperrin/poppy/issues/135#issuecomment-159693145.

Dr. Gary Pajer Senior Technical Staff Princeton Satellite Systems 6 Market Street, Suite 926 Plainsboro, NJ 08536-2096 609 275-9606 http://www.psatellite.com Visit our blog http://blog.psatellite.com.

Two more comments.

Context: the star camera will (hopefully) be installed on satellites, and used to help determine the satellite's attitude. Certain improvements are required.

Input pupil sampling: I saw the docstring mention of oversampling as a factor "over Nyquist" which you define as 2*D/lambda. I tentatively took that to have something to do with the pupil sampling. Assumption there is that there will be no feature smaller than that. I get that. But the default is oversample = 1. (As I said earlier, I don't know where the pupil sampling is actually set.) In my case, that means, well, an enormous array of samples at the pupil. Unnecessarily so in my case, because the amplitude will be slowly varying over the pupil. An OPD defined by a few Zernike polynomials. All this left me uncertain about how to use POPPY in my use case.

On Wed, Nov 25, 2015 at 4:52 PM, Gary Pajer gpajer@psatellite.com wrote:

THanks for you detailed response. I kept my original question short because I was afraid I was asking the question in the wrong place (developer's planning forum), so I eliminated many details. I'll be studying what you said more closely when the holiday is over.

I'm familiar with the DFT and issues of aliasing, sampling, periodicity, (windowing ...). One thing I couldn't figure out is what the sampling rate in the pupil is. Is pupil sampling at 100 pixels the default, or did I unknowingly set it somewhere? I take your words to mean that I did manage to correctly (modulo x2) the sampling/resolution of my focal plane array. I was (and still am) a bit curious about the oversampling parameter. I'll look more closely at that, too,

The single lens would be just for simplicity during my learning curve. I will be modeling lenses, so I will have some idea of aberrations. I've already built an element classes that introduce tilt and coma, and the ThinLens class provides defocus. The models will be used for testing centroiding algorithms. The final product is a star camera which will be defocused so that the intensity distribution covers at least a 3x3 array of pixels.

I have just noticed that the git repository has a much-updated version which includes Fresnel diffraction. THat might be more appropriate for my use case. I haven't yet installed it.

In the meantime, as a potential alternative, I'm implementing the Kirchoff-Fresnel diffraction integral, aberrated verstion, from Born and Wolf Sec 9.4. It's slow, but that might be ok for me. Poppy would be more convenient and faster.

I'm sorry but appreciative that my spartan question lead you to reply in such detail.

On Wed, Nov 25, 2015 at 1:21 PM, Marshall Perrin <notifications@github.com

wrote:

There are a couple things going on here:

The 2048 vs 1024 pixels difference is simply explained by the default oversampling factor of 2x in calcPSF. This is a user-adjustible parameter so you can make it whatever you want. (and the Instrument subclass provides a more sophisticated interface to this including automatic downsampling to the detector resolution but that takes some more work to set up.)

Secondly, the repeated copies of the PSF are due to the replication that is assumed by the FT algorithm, as @josePhoenix https://github.com/josePhoenix pointed out. The maximum spatial frequency that can be sampled is set by the sampling in the input pupil plane. If you have an pupil array that is 100 pixels across, you can represent up to 50 cycles per pupil (at Nyquist sampling, 2 pixels per cycle). This means you can represent up to 50 lambda/D in the image plane, and after that you will get aliasing. See for instance the "Sampling Sinusoidal Functions" section of https://en.wikipedia.org/wiki/Aliasing. The software attempts to calculate the FT of a high spatial frequency, but due to limited sampling it returns instead a value appropriate for a much lower frequency and this results in repeated copies of the (low spatial frequency!) core of the PSF. The fix for this is to increase t he sampling of the input pupil array. The case you are interested in may actually be better represented by the new FresnelOpticalSystem functionality in our recent release. That class provides a more convenient interface to adjust the sampling of the input wavefront, since you can just specify npix and pupil_diameter directly in the creation of the FresnelOpticalSystem. (This newer code now has a better API than our older code I think... Lessons learned!)

All that said, your proposed system is significantly undersampled spatially. For a 50 mm lens the theoretical diffraction limited PSF is just a few arcseconds at visible wavelengths. It's at least an order of magnitude better than your stated 70 arcsec/pixel sampling. So the kinds of fine diffractive structure that poppy is designed to model will not be easily seen in your instrument. Furthermore, remember that poppy is just modeling the effects of diffraction, based on the input parameters describing an optical system. In your case with a simple circular pupil that will just give an Airy function, without any need for numerical calculations of a more complex PSF. More to the point, your instrument is not going to have diffraction limited optical performance with only a single lens. It takes multiple powered surfaces to get near diffraction limited performance over an extended field, particularly for anything as large as 20 degrees across. Depending on your application, y ou're likely in a regime dominated by geometric optics rather than diffraction.

Your given parameters (50 mm focal length, 25 mm aperture diameter) are in the realm of commercial camera lenses, but even those don't give nearly diffraction limited performance over such a wide field. See for instance this Nikon 50 mm f/1.8 lens: http://imaging.nikon.com/lineup/lens/singlefocal/normal/af-s_50mmf_18g/ That page has a chart showing the modulation transfer function (MTF), which depicts how the MTF declines significantly as you move away from the optical axis. And that's for a seven lens system, one of which is an aspheric surface. A single lens is going to be totally dominated by low-order aberrations off axis. It's possible to model such effects in poppy, but that requires having an input model describing the wavefront error as a function of field angle, for instance as derived from the optical model of the lens.

— Reply to this email directly or view it on GitHub https://github.com/mperrin/poppy/issues/135#issuecomment-159693145.

Dr. Gary Pajer Senior Technical Staff Princeton Satellite Systems 6 Market Street, Suite 926 Plainsboro, NJ 08536-2096 609 275-9606 http://www.psatellite.com Visit our blog http://blog.psatellite.com.

Having given this some thought, I think I could make progress if I can understand how/where the sampling of the pupil is set. Before dealing with oversampling, I'd like to figure out how to have my pupil plane and detector plane have the same size arrays. I'm most familiar with the DFT when the direct space and transform space are the same size.

Hi @gapster

First let me apologize for not responding back in December - as you can imagine that was a busy time of year for a lot of reasons! Please let me know if you are still having questions with this and I will work to help get them cleared up.

The default pupil sampling is set in a non-obvious way: either it's hard coded as 1024, or else if the first optic in the OpticalSystem has a specified shape that is used instead. In your example code above the CircularAperture doesn't have a specific size, so the 1024 is used. This is in OpticalSystem.inputWavefront().

The Fresnel Optical System class is better since it has an npix parameter in the optical system constructor. We should move that to the Fraunhofer code too.

Thanks Marshall,

I still do have questions, but this is now on the "side burner" (not quite the back burner). I've been looking at another approach which may or may not be more appropriate for me. What I'm looking at is the Extended Nijboer-Zernike http://www.nijboerzernike.nl/solutions to the diffraction integral. I've also discovered that a friend of mine is a friend of Bob Vanderbei, who is a co-author of the paper on the Matrix Fourier Transform by Soummer which the Poppy site cites. Furthermore, Bob can't be more than a couple of miles from me in Princeton. We've exchanged emails, but we haven't discussed any of this yet.

I hope to get back to Poppy before too long, but that time is not now. Regards!

On Thu, Jan 21, 2016 at 2:05 PM, Marshall Perrin notifications@github.com wrote:

Hi @gapster https://github.com/gapster

First let me apologize for not responding back in December - as you can imagine that was a busy time of year for a lot of reasons! Please let me know if you are still having questions with this and I will work to help get them cleared up.

The default pupil sampling is set in a non-obvious way: either it's hard coded as 1024, or else if the first optic in the OpticalSystem has a specified shape that is used instead. In your example code above the CircularAperture doesn't have a specific size, so the 1024 is used. This is in OpticalSystem.inputWavefront().

The Fresnel Optical System class is better since it has an npix parameter in the optical system constructor. We should move that to the Fraunhofer code too.

— Reply to this email directly or view it on GitHub https://github.com/mperrin/poppy/issues/135#issuecomment-173675892.

mperrin / poppy

Help request: interpreting PSF #135