psf is not returned in middle

[richardgmcmahon]

I have downloaded version 0.3 and psf is for cutout is not centrally located compared with version I downloaded in mid 2014; see 'old' image and new 'image'

This is what I used to get: The image is 25x25 pixels and central pixel is at 12, 12 counting from 0,0
Filename: /data/desardata/Y1A1/DES0449-4706/DES0449-4706_g_psfcat.psf
poldeg: 2
psf_samp: 1.11690831
len(psf_mask.flat): 3750
psf_mask.shape: (1, 6, 25, 25)
psf_mask.type: 3750
psf_image.shape: (25, 25)
Location of maximum value: (12, 12)

Old version of psfex pre version

This is psfex.version: 0.3.0

poldeg: 2
psf_samp: 1.11690831
len(psf_mask.flat): 3750
psf_mask.shape: (1, 6, 25, 25)
psf_mask.type: 3750
psf_image.shape: (29, 29)
Location of maximum value: (13, 13)
Saving: psfex_test.png

Also, the output is larger; increasing from 25x25 to 29x29: I assume this related to PSF_SAMP which 1.11690831 in this case; 25x 1.11690831 = 27.92; I assume it is 29 so that the psf can be centred on a pixel sine 28 would put into in the psf centre on a pixel boundary.

This is my test code:

from __future__ import print_function, division

import os

import matplotlib.pyplot as plt

import numpy as np

import psfex
from astropy.io import fits

try:
    print('psfex.__version__: ', psfex.__version__)
except:
   pass

psf_file = '/data/desardata/Y1A1/DES0449-4706/DES0449-4706_g_psfcat.psf'

hdulist = fits.open(psf_file)
hdulist.info()
print(hdulist[0].header)
print(hdulist[1].header)


poldeg =hdulist[1].header['POLDEG1']
print('poldeg: ', poldeg)

psf_samp =hdulist[1].header['PSF_SAMP']
print('psf_samp: ', psf_samp)

data=hdulist[1].data
psf_mask=data['PSF_MASK']

print('len(psf_mask.flat): ', len(psf_mask.flat))
print('psf_mask.shape: ', psf_mask.shape)
print('psf_mask.type:  ', psf_mask.size)


pex = psfex.PSFEx(psf_file)

row=5000
column=5000

psf_image = pex.get_rec(row, column)
print('psf_image.shape: ', psf_image.shape)

imax=np.argmax(psf_image)
xymax = np.unravel_index(imax, psf_image.shape)
print('Location of maximum value: ',xymax)

plt.figure("psfex",figsize=(8.0, 8.0))


plt.imshow(psf_image, interpolation='nearest')
plt.colorbar()
plt.ylabel('Pixels')
plt.xlabel('Pixels')

filename = os.path.basename(psf_file)
title= filename
plt.title(title)

plotfile='psfex_test.png'
print('Saving: ', plotfile)
plt.savefig(plotfile)

[esheldon]

Hmm, I thought we fixed this. @erykoff @beckermr ?

[beckermr]

This one is all @erykoff.

[erykoff]

grumble.
I will take a look.

[erykoff]

@richardgmcmahon A couple of things. First, I was unable to exactly replicate your output...which is strange. I was getting a 23x23 output with the psf model file downloaded from the archive; you may be using a different psf file that you ran yourself? In any event, you can get the center of the reconstructed psf with the pex.get_center(row,col) function. And this is not necessarily the center of the postage stamp. But unlike the previous version of the code (which was wrong), the center of the reconstruction is precisely at the position returned by get_center.

For your example it would be (12,12). But note that if you had said get_center(5000,5000.2) it would be centered at (12,12.2). So you can simply line up the reconstructed image along a pixel boundary and the centroid will match the centroid of the star you are fitting.

[richardgmcmahon]

@erykoff Thanks for looking into this.

OK, there are two DESDM runs for that tile so I will download both and run my tests again.

20130819000011_DES0449-4706 [this is SVA1]
20141105000021_DES0449-4706 [this is Y1A1]

I can run my tests on the same tile as you used too.

In my example, above I used 5000, 5000 so I expected the reconstruction to be centered. I was expecting the returned psf to be centered and was going to interpolate with scipy.ndimage.interpolation.shift although I could just specify a a fractional pixel for pex.get_rec which as my next step.

[richardgmcmahon]

I get an output of 23x23 for the SVA1 version of this tile; 20130819000011_DES0449-4706

See image below:

NOTE the SVA1 and Y1A1 images have different PSF_SAMPLE values

SVA1: psf_samp: 0.91619998
Y1A1 psf_samp: 1.11690831

The peak is offset by 1 pixel from centre in both x and y; so I think that there is an off-by-one feature.

pex = psfex.PSFEx(infile_psf)
row=5000
column=5000
print('pex.get_center: ', pex.get_center(row, column))
psf_image = pex.get_rec(row, column)
print('psf_image.shape: ', psf_image.shape)
imax=np.argmax(psf_image)
xymax = np.unravel_index(imax, psf_image.shape)
print('Location of maximum value: ',xymax)

gives for SVA1

psf_samp:  0.91619998
len(psf_mask.flat):  3750
psf_mask.shape:  (1, 6, 25, 25)
psf_mask.type:   3750
pex.get_center:  [ 10.  10.]
psf_image.shape:  (23, 23)
Location of maximum value:  (10, 10)

for the Y1A1 tile 20141105000021_DES0449-4706 I get:

psf_samp:  1.11690831
len(psf_mask.flat):  3750
psf_mask.shape:  (1, 6, 25, 25)
psf_mask.type:   3750
pex.get_center:  [ 13.  13.]
psf_image.shape:  (29, 29)
Location of maximum value:  (13, 13)

[erykoff]

Considering that FITS doesn't have a standard whether images should be 0- or 1-indexed, is that what's going on?

[richardgmcmahon]

@erykoff The FITS standard does say that images are 1-indexed; aka FORTRAN indexing.

This often causes confusion when working with WCS transformations between pixels and RA, Dec.

The original standard is defined in: Pence+2010: http://adsabs.harvard.edu/abs/2010A&A...524A..42P

http://www.aanda.org/articles/aa/full_html/2010/16/aa15362-10/aa15362-10.html#S196

8.1. Basic concepts

Rather than store world coordinates separately for each datum, the regular lattice structure of a FITS image offers the possibility of defining rules for computing world coordinates at each point. As stated in Sect. 3.3.2 and depicted in Fig. 1, image array data are addressed via integral array indices that range in value from 1 to NAXISj on axis j. Recognizing that image data values may have an extent, for example an angular separation, spectral channel width or time span, and thus that it may make sense to interpolate between them, these integral array indices may be generalized to floating-point pixel coordinates. Integral pixel coordinate values coincide with the corresponding array indices, while fractional pixel coordinate values lie between array indices and thus imply interpolation. Pixel coordinate values are defined at all points within the image lattice and outside it (except along conventional axes, see Sect. 8.5). They form the basis of the world coordinate formalism in FITS depicted schematically in Fig. 2.

Here is another example psf which has PSF_SAMPLE = 1.0 from:

http://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/en/megapipe/cfhtls/psf.html
http://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/data/pub/CFHTSG/psfex.W1-1-1.G.psf

poldeg:  3
psf_fwhm:  3.66452312
psf_samp:  1.0
len(psf_mask.flat):  9610
psf_mask.shape:  (1, 10, 31, 31)
psf_mask.type:   9610
pex.get_center:  [ 14.  14.]
psf_image.shape:  (31, 31)
Location of maximum value:  (14, 14) [0-index from numpy]

[esheldon]

the code I just included from eli now puts the image at the center, on average