spacetelescope / poppy

Physical Optics Propagation in Python
https://poppy-optics.readthedocs.io
BSD 3-Clause "New" or "Revised" License
214 stars 71 forks source link

MTF computation #331

Closed FredericCassaing closed 3 years ago

FredericCassaing commented 5 years ago

Hi, I was used to compute the MTF (Modulation Transfer Function), ie the Fourier Transform (FT) of the PSF, which is supposed to be the pupil autocorrelation. This worked fine when I used FFTs. Now, the MTFs I compute from POPPY's PSF suffer from support issues (the decreasing rate is rather small out of the expected support, just in the X & Y direction but not in diagonal). I suspect this comes from the fact that I compute the MTF as the FFT of the PSF, whereas the PSF results from a MFT (Matricial FT) of the (sampled) aperture. Is there a way to overcome this ? Should I compute the MTF as the MFT^-1 of the PSF ? But is the MFT even invertible, and if so, does the "classical" equality still hold with the MFT instead of the FFT : FT( a * b) = FT(a) \conv FT(b) ? => Is there a way to get a "nice" MTF with POPPY ? Thank you for your answers, Regards, F. Cassaing

mperrin commented 5 years ago

Yes, certainly the MFT is invertible and the same equation for convolution still holds.

The discrete Fourier transform is a unique and well-defined mathematical operation; the FFT and MFT are two different ways of computing it but yield precisely the same results (to within floating point numerical precision ) when computed for the same sampling. The FFT uses some numerical symmetries and shortcuts for a faster calculation, at the cost of imposing particular relations on the pupil and image sampling. The MFT is more flexible about sampling but slower. You can also directly calculate an arbitrary discrete FT directly (e.g using a for loop to iterate a sun over pixels and sine factors); that would be the most flexible but slowest approach. In any case all should yield equivalent results. The equivalence of FFT and MFT is verified as part of the poppy test suite.

So I would expect your MTF calculation to work well. Not sure why you seem to be getting surprises with the support. Perhaps you could share example code and figures showing these results? That could help to investigate.

LucasMarquis commented 5 years ago

Hello M. Perrin,

I am currently working with Frederic Cassaing. Thanks for your help. Thus, I submit you the typical code that results in a quite noisy Optical Transfer Function, calculated by FFT(psf_poppy). Please find attached the related imshowed picture, with a log10 scale to make the issue visible. Some pixels out of the theoretical MTF limit are still as high as a few percents of those before MTF limit, with a typical position : on horizontal and vertical lines passing by theoretical MTF pixels. This noisy signal becomes larger and disrupting as the number of apertures increases (by summation, multiHexAperture).

Sincerely yours, Lucas MARQUIS

import poppy import numpy as np import matplotlib.pyplot as plt osys = poppy.OpticalSystem() osys.add_pupil(poppy.HexagonAperture(side=0.1)) # (meters) osys.add_detector(pixelscale=0.10, fov_pixels=256.0) Psf = osys.calc_psf(0.5e-6) # wavelength (meters) psf=Psf[0].data mtf=np.abs(np.fft.fftshift(np.fft.fft2(psf))) mtf/=np.max(mtf) #just to compare easily the values on the following image plt.imshow(np.log10(mtf)) ;plt.colorbar(); plt.title('MTF of single Hexagonal aperture - log10'); plt.show()

mtfProblem

mperrin commented 5 years ago

Hi @LucasMarquis - I can confirm that I see consistent results here when I run your example code.

It appears to me that this may be related to the finite field of view for the PSF calculation. As I increase the fov_pixels parameter to larger sizes, the artifacts become smaller (though not zero). Sampling to a finite, square FOV in the focal plane corresponds to convolution with the FFT of a square aperture once transformed back to the pupil plane - hence the extended structures in the MTF along the X and Y axes. You can think of this like imposing a spatial filter in the pupil domain, based on the aperture set in the image domain. This is a nonzero but quite small effect, which your chosen image stretch is magnifying (dynamic range of 1e18 in a single plot is extremely large!)

You may be able to minimize this effect by setting the PSF sampling to use the full set of spatial frequencies supported on the input image - i.e. set the focal plane pixel scale and extent to match that used by the FFT. Or just use the FFT'ed versions, if that works better in your tests.

Alternatively, since the artifacts here are so small, you could just apply some threshold in the MTF image, and discard any pixels that are less than for example 1e-5 of the peak of the MTF. Depends on what exactly is needed for your application.

Hope that helps.