Point sources in oversampled images are MUCH more centrally-concentrated than the model PSF

[benw1]

For example, if we compare an oversample=10 image to the output model PSF file:

• The F062 model PSF shows an average central brightness of about 0.32, and average of 0.08 at +/- 1 pixel in X or Y, for a ratio of ~4:1 peak to adjacent pixels. The simulated image with oversample=10 nearest chip center shows a central pixel brightness of 161000 counts above background, while the four adjacent pixels have an average brightness of 12000 counts above background, for a ratio of ~13:1. This is vastly more centrally concentrated than the model.
• The F184 model PSF shows a central brightness of about 0.18, and average of 0.07 at +/- 1 pixel in X or Y, for a ratio of ~2.5:1 peak to adjacent pixels. The simulated image with oversample=10 nearest chip center shows a central pixel brightness of 22000 counts above background, while the four adjacent pixels have an average brightness of 3000 counts above background, for a ratio of ~7:1. The extreme central concentration seen in the oversample=10 examples shown here are seen in all six filters.

Thus, there appears to be a problem in the routines that are adding the PSF to the pixel grid when oversampling.

[york-stsci]

I don't believe I changed the output file structure.

If you look at the most recent pre-release tag (or the most recent version) you'll see a fixed version of the PSF enclosed energy. Any remaining issues are will webbpsf, because at this point any point source is calculated directly from a webbpsf evaluation.

[benw1]

Thanks Brian. It looks like when I run with the main STIPS branch, I get files that are 65M in size:

(stips) -bash-4.3$ ls -alstrdh data/grid.cat/proc_grid/sim_112766_0.fits 65M -rw-r----- 1 bwilli24 s1692 65M Sep 20 14:32 data/grid.cat/proc_grid/sim_112766_0.fits (stips) -bash-4.3$

But when I run the same code on the psf_dev branch, I get files that are 129M in size: (stips) -bash-4.3$ ls -alstrdh data/grid.cat/proc_grid/sim_112762_0.fits 129M -rw-r----- 1 bwilli24 s1692 129M Sep 20 09:50 data/grid.cat/proc_grid/sim_112762_0.fits (stips) -bash-4.3$

So something seems to be different about the output file format. Any thoughts on what I can check?

[benw1]

The issue we're seeing is that the shape of the point sources in the images does not match the shape of the PSFs in the PSF files. How is the consistency between these a webbpsf issue? STIPS should be able to make these consistent by using the same webbpsf output in both places.

[benw1]

Should I return to the main branch (STScI-STIPS)? Or do I only need to do a git pull on the psf_dev branch I'm currently using (STScI-STIPS-1)? Will this version also require Python 3.8 and the other version dependencies of the psf_dev version?

[york-stsci]

The file size issue is a result of STIPS moving from 32-bit (single-precision) floating point values to double-precision (64-bit).

STIPS is using the same webbpsf output in both places. When STIPS models a point source, it does the following:

• Evaluate the PSF grid at the source location (so interpolate between the four nearest images for the detector location, and then interpolate the sub-pixel grid for location on the pixel. This is all done by the GriddedPSFModel evaluate() function, not by STIPS itself (note that the GriddedPSFModel is the same way that webbpsf itself stores PSF grids). Note that source flux is also incorporated at this stage by the evaluate() function.
• Bin the PSF to detector resolution (if oversampled).

If the PSF that results from this process is not correct, then either evalute() or binning is producing a problem. STIPS no longer convolves point sources by a PSF model at all, it instead simply adds in a generated PSF at the location of the point source.

All current changes are on the psf_dev branch. Nothing is going to be merged into the main STIPS repository until we're at the release candidate stage. At this point, new releases of STIPS will continue to require python 3.8+.

[benw1]

Thanks! I did a git pull on the psf_dev branch and many updates occurred. I reinstalled and am running a test. So far it looks like it is still overwriting PSF files that are already there, which significantly slows down the process, especially when oversampling is > 1.

[york-stsci]

The current version of STIPS creates 3 PSF files for a given set of inputs (_orig.fits, _epsf.fits, and _ipc.fits). To confirm, you had already created all three of those files, and it's now creating them again?

[york-stsci]

Also, the major changes in the most recent version are to the method of PSF creation and normalization. If you are using a cached PSF created with an older version, you will not see those changes, so if the problem you are seeing WRT PSF sharpness persists, the first thing to do is to make sure that you're using a PSF created by the current (2.0.0.dev2) development version of STIPS.

[benw1]

Got it! Ok, I do not have any fully completed runs with all 3 files in place yet. I'll test the overwriting issue when these runs finish. It does appear to be producing these new files, which hopefully is a good sign.

[benw1]

My test of oversample=1 with the latest version is complete. The attached screenshot is from the R062 band. I was expecting the stars to all be identical with OS=1,but they aren't. Is that something that should be of concern? Perhaps it is doing some interpolating, but I wasn't thinking it would do that with OS=1. Does the new version interpolate when stars aren't centered on the pixel even with oversample=1?

[benw1]

Rerunning does not overwrite the PSF files in the psf_cache/

Much faster now! Thanks for fixing that!

[york-stsci]

Yes, STIPS always uses evaluate() to create a PSF for every point source, regardless of grid size and PSF oversample. This avoids needing to use special-purpose code, with its associated higher chances of bugs.

[benw1]

Is the double precision necessary? The larger files are significantly more difficult to work with,and it seems overkill given the Poisson error limit to the final photometry. Is there any way to choose 32 bit to save All that storage and I/O?

[benw1]

First test does't look quite right. Attached the left image is oversample=6, and the right image is the same pixels in an image using oversample=1 with matched stretch. The point sources in the oversample=6 image still appear too sharp assuming the oversample=1 is the correct one. I am going to wipe my psf_cache and rerun both to make sure everything is up to date, but it looks like something is still wrong with the binning in the oversample>1 images

.

[york-stsci]

The binning is done as follows:

psf_img = psf_img.reshape(psf_y//os, os, psf_x//os, os).sum(axis=(1, 3))

where psf_y and psf_x are the size of the PSF y and x axis, and os is the PSF oversample. If you can find an error in it, I'm happy to fix it. I can't.

[york-stsci]

Single-precision floating point does not have the necessary accuracy for the calculations that STIPS does. It probably does have sufficient accuracy for the output files.

[benw1]

Thanks! I'll take a look at it with our photometry developer. Out of curiosity, why is it psf_y//os and not psf_y/os? Also, would it be particularly difficult to make it so that the output file is single-precision, or put in an option for that? It would save lots of space,I/O, and memory for analysis on the images.

[york-stsci]

It uses psf_y//os rather than psf_y/os in order to ensure integer results, because the values are defining the size of the resulting array, and an array must have an integer number of elements.

[benw1]

I wonder if something odd is happening with the evaluate(), even at oversample=1. In the oversample=1 F062 image, some stars look like the PSF in the PSF file, but others do not, and I wonder if it is to do with the interpolation. See the attached screenshot where the upper left is the PSF in the PSF file; the upper right is one of the stars in the image that resembles this PSF, and the lower left is a nearby star that doesn't really resemble the PSF in the PSF file (see the dark rings that aren't present in the PSF in the file

).

[benw1]

I'm working on analyzing and comparing PSF files with images. Is the epsf file the one with the PSF that is being added tp the images? What are the 2 additional PSF files? I want to make sure I'm comparing the correct things.

[york-stsci]

The PSF files are as follows:

• epsf: generated PSF with IPC added and pixel flux binned as per Jay Anderson's definition
• ipc: generated PSF with IPC added
• orig: generated PSF

Note that, if residual_ipc is set to True in the configuration file, then IPC effects will be added to point sources and Service profiles. The PSF to use is set by the configuration flag psf_type, and defaults to 'orig'

[benw1]

It looks like the psf files may not be produced at the correct oversample. Attached is a screenshot showing 3 psf_orig files, at oversample=1. oversample=6, and oversample=10. They appear to all have the same effective sampling of the PSF, although the size of the pixel grid is larger for larger oversample values. In the lower right is the oversample=10 PSF file from STIPS 1.0.8, which clearly is highly oversampled compared to the 3 from this development version. This could explain why the higher oversamples are resulting in PSFs that are too sharp in the images.

[york-stsci]

I think I found the issue. The above file is my attempt at a fix. Top left is oversample=1, top right oversample=2, bottom left oversample=3. Do these look more like you would expect?

[benw1]

Yes! Is it committed? I'm happy to try it out. Thanks!

[york-stsci]

I have just pushed the changes. Note that the tests I ran last night overnight did not complete because I made a silly typo in the test function, so I haven't confirmed that I didn't do anything foolish to break things, but at the very least my example test runs appear to complete okay.

[benw1]

I pulled the changes, reinstalled, wiped my psf_cache, and made a new oversample=10 simulation. Unfortunately, the new oversample=10 PSF still looks poorly sampled. The attached screenshot shows the new oversample=10 PSF file (left), compared to the one from STIPS 1.0.8 (right).

[benw1]

The new version did not make 3 PSF files, just one called: psf_v2.0.0dev3_roman_wfi_sca01_f087_over10_grid1.fits Not sure if that helps, but thought I'd mention it since it was different that before, when 3 PSF files were generated.

[york-stsci]

That's what it's supposed to do. Please check the same oversampled that I did (2/3/4) to see if maybe there's a specific issue with OS=10, or if we're just getting very different output from doing the same thing.

[benw1]

I will try, but looking at yours again, I may have been fooled by the stretch and the zoom. Can you set them so that each panel shows the same amount of the PSF? It may be that the sampling is not actually different.

[benw1]

What I mean by "same amount of the PSF" is that if one panel shows the inner 45 pixels of the PSF file, the other file at there different sampling should also show the inner 45 pixels of the PSF file, zoomed to match sizes on the screen.

[benw1]

When I do that with oversample=10 and oversample=1, I get this, where the oversample=10 is the PSF file on the left. The sampling is essentially the same in both PSF files.

[york-stsci]

Why do you want to compare the same number of pixels at different pixel scales? I can find images, but I don't understand what it's telling you to compare 4.5 arc seconds of one PSF against 0.45 arc seconds of another.

[benw1]

I don't actually believe that the inner 45 pixels of the oversample=10 PSF is showing 0.45 arcsec. If it is, then that would suggest the PSF FWHM is <0.01 arcsec, as most of the flux is in the central pixel, even at this oversample rate.

[york-stsci]

Oversample 1 on the left, oversample 3 on the right. Same number of pixels of each. Used ds9's "Frame > Match Scale and Limits" to match scales and limit s.

[york-stsci]

Note that these PSFs are normalized based on what GriddedPSFModel expects, and GriddedPSFModel expects ePSFs as defined by Jay Anderson. As such, when GriddedPSFModel evaluates a PSF, it divides the PSF it produces by oversample^2. Because webbpsf doesn't include this factor, before I store the PSFs as a FITS file I multiply them all by a factor of oversample^2, so if you're trying to compare how centrally concentrated the flux is, you have to include that scaling.

[benw1]

That example verifies my concern. If the sampling were actually different, the oversample=3 would only be showing a third of the radius of the PSF, but it shows the same radius of the PSF, so the sampling appears to be unchanged. Only the normalization changes.

[york-stsci]

I've just made a push that will default the PSF grid creation to using a different form of setting the oversample, because the one I'm currently using appears to be buggy. I haven't done any significant testing on it, I just went ahead and did a push, because this is a development branch of my fork, and hopefully that lets you take a look at it as soon as possible.

Please let me know if it improves matters.

[york-stsci]

If you already looked at the new code, please do another pull. I was using the wrong extension of the PSF file. Sorry.

[york-stsci]

And yet another one. That will teach me to try to save time in sending things to you. For reference, you want the version to show up as 2.0.0dev6.

[york-stsci]

Make that 2.0.0dev7. Who would have thought that I could make so many typos in so little code.

[benw1]

Looking better! I ran both dev6 and dev7, and I don't see a difference in the PSF files, but they both look good. Attached shows dev6 and dev7 oversample=10 on the left. Upper right is the dev3 oversample=1, and lower right is 1.0.8 oversample=10. The new ones look better than that 1.0.8 one, but in any case, the all appear to have the correct sampling (these are showing the inner 200 pixels over the oversample=10 PSFs and the inner 20 pixels of the oversample=1). Thanks for getting that part fixed!

However, the stars in the images don't look right, given the PSF file. In an oversample=10 image, all the stars look perfectly centered on a pixel, and too sharp, in the final simulated image, as shown in the second attachment. In this image, the stars have central pixel values of ~150000, and the neighboring pixels have values of ~5000, for a ratio of ~30. While if I bin the PSF file by a factor o 10, the values are ~45 and ~5, for a ratio of ~9. What's a bit more odd, is that some stars have central pixel values of only ~70000, but the neighboring pixels still have values of only about 5000, so somehow the flux is not being conserved in the addition either. Thus, the problem in the PSF files looks to be possibly completely fixed, but something is still off in how the PSFs are being added to the images.

[york-stsci]

The way that PSFs are added to images has not changed. Use the GriddedPSFModel evaluate() function with the source position and flux to generate the source on the detector, add IPC if it's set to on, and then add to the detector.

[york-stsci]

And the PSF between dev6 and dev7 should look the same. The difference is that I moved from getting an oversampled PSF by using the oversample= keyword when calling calc_psf, to getting an oversampled PSF by dividing the pixel scale by the oversample, and multiplying the fov by the oversample. This works around a bug in webbpsf where it can't create odd parity PSF images with even oversamples, except through this method.

[benw1]

Sounds good. The PSF files are definitely improved, but I'm not sure they are consistent across oversample values. Still need to figure that out in detail. However, the images are very inconsistent across oversample values. In the attached image, the upper panels are the same input star in an oversample=1 image (left) and an oversample=10 image (right). The oversample=1 image seems to do a better job of handling the fact that the star is not centered in the pixel. There is little (no) resemblance between them. The PSF files are in the lower panels, where the left is the oversample=1 PSF file and the right is the oversample=10 PSF binned 10x10 for easy comparison. While their shapes are similar to one another, the oversample=10 PSF is a quite a bit less sharp than the oversample=1. In detail the ratio of neighboring pixel to central pixel in the oversample=1 PSF is ~20, and the ratio of neighboring pixel to central pixel in the oversample=10 PSF rebinned 10x10 is ~9. These two issues are probably distinct, but both may be related to the clear difference in the resulting image.

[york-stsci]

The one difference between an oversample-1 PSF and an oversample-10 PSF in STIPS is that the oversample-10 PSF is binned after evaluation. If there's a different way that STIPS should reduce an oversampled PSF to a detector-sampled PSF, I'm happy to implement it, but failing that there isn't anything in the STIPS code that can change how PSFs are applied to point sources.

[benw1]

Thanks! The binning may not be the problem anymore. It looks more like the position or interpolation or something like that. There's even quite a difference between the total counts in the star between the different oversamples. The oversample=1 star has ~200000 counts and the oversample=10 has ~100000 counts. There's probably a bug still hiding somewhere causing this difference, as the problem is not subtle. The same input star should not appear so different depending on the oversample chosen. Perhaps we could have it output some of the intermediate steps to try to isolate the bug? Also, only one PSF file is being generated now, instead of 3. Did you find generating 3 was no longer necessary?

[york-stsci]

I can provide intermediate outputs, yes. It's also worth noting that the generated PSF is normalized by webbpsf, and the source flux is chosen by GriddedPSFModel. STIPS does not do anything to the code that would change the flux produced by these tools.

There is a single PSF file because that file consists of multiple extensions, which include all of the data previously contained in separate files.

[benw1]

Thanks! The generated PSF files look OK now, maybe not perfect, but not way off. It's something in the application of these files that seems off. Is there a step where the PSF has been added to the image, but no binning has been done? Looking at that product may help to rule out the binning as a possible cause. I also wonder if there is some way to check the routine that is doing the interpolation when adding the PSF to the oversample=1 images, since those stars seem to bear so little resemblance to the PSF in the PSF file, including adding large negative values to some pixels. Also, above you mentioned PSFs are normalized based on what GriddedPSFModel expects, and GriddedPSFModel expects ePSFs as defined by Jay Anderson. As such, when GriddedPSFModel evaluates a PSF, it divides the PSF it produces by oversample^2, so perhaps something is a bit off there.

[york-stsci]

There is no such step, because there is no oversampled image. The only difference is whether the PSF itself has to be binned after being evaluated and before being added to the image.

The routine that is doing the interpolation is photutils.psf.GriddedPSFModel.evaluate(), if it needs changes then the issue needs to be filed with them.

Yes, PSFs are normalized to what GriddedPSFModel expects, in that, before I create the NDData array that holds the PSFs, I multiply them by oversample^2 so that, when evalute() is called, they will be divided by oversample^2 and return whatever flux webbpsf put in the PSF that it generated in the first place.

[benw1]

Is there any point where the model PSF that is being added to the image can be output at a later stage than the PSF file itself? That way we can check its consistency between the PSF being added to the image and the one in the PSF file. Something is going wrong in there, and we need to quantify it to put in an issue with photutils.