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
make clean
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.