Help request: interpreting PSF

[gapster]

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.

[josePhoenix]

@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)

[mperrin]

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.

[gapster]

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
#135 (comment).

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.

[gapster]

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
#135 (comment).

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.

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.

[gapster]

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.

[mperrin]

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.

[gapster]

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
#135 (comment).

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.