GNU Astronomy Utilities

2.3.6 Uniting the different PSF components

In Building outer part of PSF we built the outer part of the extended PSF and the inner part was built in Inner part of the PSF. The outer part was built with very bright stars, and the inner part using fainter stars to not have saturation in the core of the PSF. The next step is to join these two parts in order to have a single PSF. First of all, let’s have a look at the two stacks and also to their radial profiles to have a good feeling of the task. Note that you will need to have TOPCAT to run the last command and plot the radial profile (see TOPCAT).

$ astscript-fits-view outer/stack.fits inner/stack.fits
$ astscript-radial-profile outer/stack.fits -o outer/profile.fits
$ astscript-radial-profile inner/stack.fits -o inner/profile.fits
$ astscript-fits-view outer/profile.fits inner/profile.fits

From the visual inspection of the images and the radial profiles, it is clear that we have saturation in the center for the outer part. Note that the absolute flux values of the PSFs are meaningless since they depend on the normalization radii we used to obtain them. The uniting step consists in scaling up (or down) the inner part of the PSF to have the same flux at the junction radius, and then, use that flux-scaled inner part to fill the center of the outer PSF. To get a feeling of the process, first, let’s open the two radial profiles and find the factor manually first:

  1. Run this command to open the two tables in TOPCAT:
    $ astscript-fits-view outer/profile.fits inner/profile.fits
  2. On the left side of the screen, under “Table List”, you will see the two imported tables. Click on the first one (profile of the outer part) so it is shown first.
  3. Under the “Graphics” menu item, click on “Plane plot”. A new window will open with the plot of the first two columns: RADIUS on the horizontal axis and MEAN on the vertical. The rest of the steps are done in this window.
  4. In the bottom settings, within the left panel, click on the “Axes” item. This will allow customization of the plot axes.
  5. In the bottom-right panel, click on the box in front of “Y Log” to make the vertical axis logarithmic-scaled.
  6. On the “Layers” menu, select “Add Position Control” to allow adding the profile of the inner region. After it, you will see that a new red-blue scatter plot icon opened on the bottom-left menu (with a title of <no table>).
  7. On the bottom-right panel, in the drop-down menu in front of Table:, select 2: profile.fits. Afterwards, you will see the radial profile of the inner stack as the newly added blue plot. Our goal here is to find the factor that is necessary to multiply with the inner profile so it matches the outer one.
  8. On the bottom-right panel, in front of Y:, you will see MEAN. Click in the white-space after it, and type this: *100. This will display the MEAN column of the inner profile, after multiplying it by 100. Afterwards, you will see that the inner profile (blue) matches more cleanly with the outer (red); especially in the smaller radii. At larger radii, it does not drop like the red plot. This is because of the extremely low signal-to-noise ratio at those regions in the fainter stars used to make this stack.
  9. Take your mouse cursor over the profile, in particular over the bump around a radius of 100 pixels. Scroll your mouse down-ward to zoom-in to the profile and up-ward to zoom-out. You can click-and-hold any part of the profile and if you move your cursor (while still holding the mouse-button) to look at different parts of the profile. This is particular helpful when you have zoomed-in to the profile.
  10. Zoom-in to the bump around a radius of 100 pixels until the horizontal axis range becomes around 50 to 130 pixels.
  11. You clearly see that the inner stack (blue) is much more noisy than the outer (red) stack. By “noisy”, we mean that the scatter of the points is much larger. If you further zoom-out, you will see that the shallow slope at the larger radii of the inner (blue) profile has also affected the height of this bump in the inner profile. This is a very important point: this clearly shows that the inner profile is too noisy at these radii.
  12. Click-and-hold your mouse to see the inner parts of the two profiles (in the range 0 to 80). You will see that for radii less than 40 pixels, the inner profile (blue) points loose their scatter (and thus have a good signal-to-noise ratio).
  13. Zoom-in to the plot and follow the profiles until smaller radii (for example, 10 pixels). You see that for each radius, the inner (blue) points are consistently above the outer (red) points. This shows that the \(\times100\) factor we selected above was too much.
  14. In the bottom-right panel, change the 100 to 80 and zoom-in to the same region. At each radius, the blue points are now below the red points, so the scale-factor 80 is not enough. So let’s increase it and try 90. After zooming-in, you will notice that in the inner-radii (less than 30 pixels), they are now very similar. The ultimate aim of the steps below is to find this factor automatically.
  15. But before continuing, let’s focus on another important point about the central regions: non-linearity and saturation. While you are zoomed-in (from the step above), follow (click-and-drag) the profile towards smaller radii. You will see that smaller than a radius of 10, they start to diverge. But this time, the outer (red) profile is getting a shallower slope and diverges significantly from about the radius of 8. We had masked all saturated pixels before, so this divergence for radii smaller than 10 shows the effect of the CCD’s non-linearity (where the number of electrons will not be linearly correlated with the number of incident photons). This is present in all CCDs and pixels beyond this level should not be used in measurements (or properly corrected).

The items above were only listed so you get a good mental/visual understanding of the logic behind the operation of the next script (and to learn how to tune its parameters where necessary): astscript-psf-scale-factor. This script is more general than this particular problem, but can be used for this special case also. Its job is to take a model of an object (PSF, or inner stack in this case) and the position of an instance of that model (a star, or the outer stack in this case) in a larger image.

Instead of dealing with radial profiles (that enforce a certain shape), this script will put the centers of the inner and outer PSFs over each other and divide the outer by the inner. Let’s have a look with the command below. Just note that we are running it with --keeptmp so the temporary directory with all the intermediate files remain for further clarification:

$ astscript-psf-scale-factor outer/stack.fits \
           --psf=inner/stack.fits --center=501,501 \
           --mode=img --normradii=10,15 --keeptmp
$ astscript-fits-view stack_psfmodelscalefactor/cropped-*.fits \

With the second command, you see the four steps of the process: the first two images show the cropped outer and inner stacks (to same width image). The third shows the radial position of each pixel (which is used to only keep the pixels within the desired radial range). The fourth shows the per-pixel division of the outer by the inner within the requested radii. The sigma-clipped median of these pixels is finally reported. Unlike the radial profile method (which averages over a circular/elliptical annulus for each radius), this method imposes no a-priori shape on the PSF. This makes it very useful for complex PSFs (like the case here).

To continue, let’s remove the temporary directory and re-run the script but with --quiet mode so we can put the output in a shell variable.

$ rm -r stack_psfmodelscalefactor
$ scale=$(astscript-psf-scale-factor outer/stack.fits \
                   --psf=inner/stack.fits --center=501,501 \
                   --mode=img --normradii=10,15 --quiet)
$ echo $scale

Now that we know the scaling factor, we are ready to unite the outer and the inner part of the PSF. To do that, we will use the script astscript-psf-unite with the command below (for more on this script, see Invoking astscript-psf-unite). The basic parameters are the inner part of the PSF (given to --inner), the inner part’s scale factor (--scale), and the junction radius (--radius). The inner part is first scaled, and all the pixels of the outer image within the given radius are replaced with the pixels of the inner image. Since the flux factor was computed for a ring of pixels between 10 and 15 pixels, let’s set the junction radius to be 12 pixels (roughly in between 10 and 15):

$ astscript-psf-unite outer/stack.fits \
           --inner=inner/stack.fits --radius=12 \
           --scale=$scale --output=psf.fits

Let’s have a look at the outer stack and the final PSF with the command below. Since we want several other DS9 settings to help you directly see the main point, we are using --ds9extra. After DS9 is opened, you can see that the center of the PSF has now been nicely filled. You can click on the “Edit” button and then the “Colorbar” and hold your cursor over the image and move it. You can see that besides filling the inner regions nicely, there is also no major discontinuity in the 2D image around the union radius of 12 pixels around the center.

$ astscript-fits-view outer/stack.fits psf.fits --ds9scale=minmax \
           --ds9extra="-scale limits 0 22000 -match scale" \
           --ds9extra="-lock scale yes -zoom 4 -scale log"

Nothing demonstrates the effect of a bad analysis than actually seeing a bad result! So let’s choose a bad normalization radial range (50 to 60 pixels) and unite the inner and outer parts based on that. The last command will open the two PSFs together in DS9, you should be able to immediately see the discontinuity in the union radius.

$ scale=$(astscript-psf-scale-factor outer/stack.fits \
                   --psf=inner/stack.fits --center=501,501 \
                   --mode=img --normradii=50,60 --quiet)

$ astscript-psf-unite outer/stack.fits \
           --inner=inner/stack.fits --radius=55 \
           --scale=$scale --output=psf-bad.fits

$ astscript-fits-view psf-bad.fits psf.fits --ds9scale=minmax \
           --ds9extra="-scale limits 0 50 -match scale" \
           --ds9extra="-lock scale yes -zoom 4 -scale log"

As you see, the selection of the normalization radii and the unite radius are very important. The first time you are trying to build the PSF of a new dataset, it has to be explored with a visual inspection of the images and radial profiles. Once you have found a good normalization radius for a certain part of the PSF in a survey, you can generally use it comfortably without change. But for a new survey, or a different part of the PSF, be sure to repeat the visual checks above to choose the best radii. As a summary, a good junction radius is one that:

Now that the complete PSF has been obtained, let’s remove that bad-looking PSF, and stick with the nice and clean PSF for the next step in Subtracting the PSF.

$ rm -rf psf-bad.fits