Extracting PSF Images

Generic Code to Reconstruct the PSF

The Point-Spread Function (PSF) reconstruction can be performed without any specialized tools.

To read the PSF info from a psField file for the r band, using IDL, read extension 3:

IDL> pstruct = mrdfits(psfield_file, 3)

(recall u,g,r,i,z == 0,1,2,3,4 so just add one to the index).

The resulting structure can then be used to reconstruct the image. The following IDL code, taken from the SDSSIDL routine "sdss_psfrec.pro", demonstrates the algorithm to reconstruct the PSF at location (row,col) in the field.

nrow_b=(pstruct.nrow_b)[0]
ncol_b=(pstruct.ncol_b)[0]
;assumes they are the same for each eigen so only use the 0 one
rnrow=(pstruct.rnrow)[0]
rncol=(pstruct.rncol)[0]
nb=nrow_b*ncol_b
coeffs=fltarr(nb)
ecoeff=fltarr(3)
cmat=pstruct.c
rcs=0.001
FOR i=0L, nb-1L DO coeffs[i]=(row*rcs)^(i mod nrow_b) * (col*rcs)^(i/nrow_b)
FOR j=0,2 DO BEGIN
    FOR i=0L, nb-1L DO BEGIN
        ecoeff[j]=ecoeff[j]+cmat(i/nrow_b,i mod nrow_b,j)*coeffs[i]
    ENDFOR
ENDFOR
psf = (pstruct.rrows)[*,0]*ecoeff[0]+$
      (pstruct.rrows)[*,1]*ecoeff[1]+$
      (pstruct.rrows)[*,2]*ecoeff[2]

The above code can be easily rewritten in other programming languages. For example, the following function allows to reconstruct the PSF in Python.

from astropy.io import fits
import numpy as np

def reconstructPSF(psFieldFilename, filter, row, col):

    filterIdx = 'ugriz'.index(filter) + 1
    psField = fits.open(psFieldFilename)
    pStruct = psField[filterIdx].data

    nrow_b = pStruct['nrow_b'][0]
    ncol_b = pStruct['ncol_b'][0]

    rnrow = pStruct['rnrow'][0]
    rncol = pStruct['rncol'][0]

    nb = nrow_b * ncol_b
    coeffs = np.zeros(nb.size, float)
    ecoeff = np.zeros(3, float)
    cmat = pStruct['c']

    rcs = 0.001
    for ii in range(0, nb.size):
        coeffs[ii] = (row * rcs)**(ii % nrow_b) * (col * rcs)**(ii / nrow_b)

    for jj in range(0, 3):
        for ii in range(0, nb.size):
            ecoeff[jj] = ecoeff[jj] + cmat[int(ii / nrow_b), ii % nrow_b, jj] * coeffs[ii]

    psf = pStruct['rrows'][0] * ecoeff[0] + \
        pStruct['rrows'][1] * ecoeff[1] + \
        pStruct['rrows'][2] * ecoeff[2]

    psf = np.reshape(psf, (rnrow, rncol))
    # psf = psf[10:40, 10:40]  # Trim non-zero regions.

    return psf

Stand Alone Code

We have also created a standalone code that serves to both read them and as a template library for inclusion in other codes. The code is available as: readAtlasImages-v5_4_11.tar.gz.

Compiling

1. make clean
2. make

If you are on a big-endian machine, remove -DSDSS_LITTLE_ENDIAN from CFLAGS in the Makefile.

Using

% read_PSF -h
Usage: read_PSF [options] input-file hdu output-file
Your options are:
       -?      This message
       -h      This message
       -i      Print an ID string and exit
       -v      Turn up verbosity (repeat flag for more chatter)

To reconstruct the z PSF (i.e., the 5th HDU) at the position (row, col) = (500, 600) from run 1336, column 2, field 51 you'd say:

%  read_PSF psField-001336-2-0051.fit 5 500.0 600.0 foo.fit

The desired PSF would appear as an unsigned short FITS file in foo.fit; the background level is set to the standard 'soft bias' of 1000. If you want a floating image, change a line in the read_PSF.c; look for

        /* create a float region */

Developer Comments

The standalone programs read_atlas_image (reads fpAtlas files) and read_mask (reads fpM files) are similar. All three are built by the same 'make' command.

I don't expect that many users will actually want to use the read_PSF executable (although it is perfectly functional). The main use of the product will probably be to link into custom built executables that need to process atlas image data. I believe that the code should be easily reused for this purpose.