Raw contrast calculation making incorrect adjustment to nominal resolution element size for blurred data

[kdlawson]

At line 168 in analysistools.py, the raw_contrast function accommodates blurred data by multiplying the the nominal resolution element size by the header entry for BLURFWHM. In many cases, this is actually decreasing the assumed size of the resolution element rather than increasing it (e.g., for F300M, where the automatic value for BLURFWHM is ~0.72). Instead, the nominal FWHM and the blur FWHM should probably be added in quadrature to correspond to the effective resolution of the data.

Additionally: the value for BLURFWHM recorded by the blur_frames function appears to be the sigma for the gaussian kernel rather than the FWHM (making the documentation for this function misleading). As a result: in raw_contrast the value for BLURFWHM will need to be adjusted from a sigma value to a FWHM to calculate the correct resolution element size.

[AarynnCarter]

@kdlawson Thanks for this. I see the issue with BLURFWHM being recorded as the sigma value, but I'm not entirely sure about adding things in quadrature. If we fix things by converting the sigma value to a FWHM, then the values should stay >1 and that will increase the resolution element size as expected. However, I'm not convinced / not familiar enough to know that simply multiplying is the correct approach either. I don't suppose you have a resource for why to add in quadrature instead?

[kdlawson]

Thanks @AarynnCarter. As an example: say a user simply sets BLURFWHM to a fixed value of 0.5 pixels instead of using the automated calculation. This results in the images being smoothed with a 0.5 pixel Gaussian kernel. In the current contrast code, the assumed FWHM would be the nominal FWHM x 0.5 — but surely the spatial resolution hasn't actually been improved by smoothing the data.

I don't have a literature reference for adding FWHMs in quadrature. However, it's pretty simple to check as follows:

1) Generate a gaussian PSF of some nominal FWHM

2) Convolve the PSF with a gaussian kernel of width BLURFWHM

3) See if the FWHM of the result is FWHM*BLURFWHM or sqrt(FWHM^2 + BLURFWHM^2)

Here's a quick example for a nominal FWHM of 10 pixels and a blur FWHM of 5 pixels:

import numpy as np
import matplotlib.pyplot as plt
from astropy.convolution import convolve, Gaussian1DKernel

fwhm0 = 10.0
sigma0 = fwhm0/np.sqrt(8. * np.log(2.))

blurfwhm = 5.0
blursigma = blurfwhm/np.sqrt(8. * np.log(2.))

psf0 = Gaussian1DKernel(sigma0, x_size=51).array
psf0 /= psf0.max()

blur_kernel = Gaussian1DKernel(blursigma)
psfblurred = convolve(psf0, blur_kernel)
psfblurred /= psfblurred.max()

fig,ax = plt.subplots()

x = np.arange(psf0.shape[0], dtype=float)
x -= x.mean()

ax.plot(x, psf0, c='C0', label='Nominal PSF')
ax.plot(x, psfblurred, c='C1', label='Blurred PSF')

ax.axvline(-fwhm0/2, c='C0', ls='dotted', label='Nominal FWHM')
ax.axvline(fwhm0/2, c='C0', ls='dotted')

fwhm_eff1 = fwhm0*blurfwhm # Taking the effective FWHM to be the product of the nominal and blur FWHMs
ax.axvline(-fwhm_eff1/2, c='C1', ls='dashdot', label='(Nominal FWHM) $\\times$ (Blur FWHM)')
ax.axvline(fwhm_eff1/2, c='C1', ls='dashdot')

fwhm_eff2 = np.hypot(fwhm0, blurfwhm) # Taking the effective FWHM to be the nominal and blur FWHMs added in quadrature 
ax.axvline(-fwhm_eff2/2, c='C1', ls='dashed', label='[(Nominal FWHM)$^2$ + (Blur FWHM)$^2$]$^{0.5}$')
ax.axvline(fwhm_eff2/2, c='C1', ls='dashed')

ax.axhline(0.5, c='k', ls='solid')
lgd = ax.legend(framealpha=1.0, edgecolor='k', fontsize=16, ncols=1, bbox_to_anchor=(0.5, 1.01), loc='lower center')

plt.show()

I think this is strictly correct only for circular PSFs and blurring kernels. It's probably close enough for the JWST PSFs (and, to be fair, this presumption is already built into the use of a single FWHM for quantifying resolution). However, you could also generate a PSF with WebbPSF, smooth it with the blurring kernel, and then measure+record the FWHM of the resulting PSF as a part of the spaceklip blurring procedure.

(Hopefully this makes sense and I'm not totally misunderstanding the details of this spaceklip step!)

[AarynnCarter]

@kdlawson Thanks for the detail, this is a great. Also never heard of np.hypot, handy!

There is some other bookkeeping that needs to be adjusted so I'll do this at the same time, and make sure to account for user inputs that aren't the automatic method also. Thanks again for catching this.

[AarynnCarter]

Hey @kdlawson, would you mind trying things out again using this branch?https://github.com/kammerje/spaceKLIP/tree/bookkeeping

I think everything is correct now but would be good to make sure before I merge things over.