# Fitting PSF to an image to get total flux (Python)

I have an image of a single unsaturated source (asteroid) that I'd like to fit guesses for various model PSFs to, so that I can then compare the total flux in each case. I've only recently begun studying PSFs, so I am a bit confused.

I know that astropy offers a package called photutils.psf which does the PSF fitting for you, but I would like to better understand what's going on behind the fitting, and how the flux of the source is resolved from it. Furthermore, when I tried using this package, I wasn't getting any outputs anyway (it said a source was found, but that it didn't meet the sharpness or roundness criteria, despite me changing those parameters around quite a bit). I have tried something simpler below (where model_data is the 2d numpy array representation of the PSF, and img is the asteroid image data):

import numpy as np
from astropy.modeling import fitting
from photutils import FittableImageModel

y, x = np.mgrid[:np.shape(img)[0], :np.shape(img)[1]]

fitter = fitting.LevMarLSQFitter()
model2d = FittableImageModel(model_data)
result = fitter(model2d, x, y, img)
print(result)
print(np.median(img), np.max(img), np.sum(img))


The output I get is then:

flux       x_0       y_0

-----     -----     -----

9282.784  466.44    531.69


2705.82 25301.56 980812434.98

However, these are not the correct centroid coordinates of the source in the image. The correct coordinates are around (300, 300). Also, to get this result, the PSF data was normalized so that the max intensity point is 1. But when I change this normalization so that the total flux (sum of all the pixels) is 1, the result for the flux after fitting I get is way higher.

I'm not sure what's happening behind the scenes here (I got lost in the documentation), and I'm not sure if I can trust any of these flux values, considering the x_0 and y_0 are always off. So I'd like to better understand -- how exactly does one find the flux of an object by fitting it to a PSF? And how should that PSF be normalized?

• If you're looking for the total flux of the system, just draw a box around your image and add up the pixel values for every pixel in your box. That's all flux is. Commented Sep 21, 2018 at 14:36
• @zephyr Yes, but how is this accomplished by fitting a PSF to the image? Commented Sep 21, 2018 at 16:41
• As Carl Witthoft answered, it doesn't make sense to apply a PSF to an extended source. Besides, even if this was a star, you wouldn't calculate flux by fitting a PSF to it first. You'd do exactly what I said and draw a box around your object and add up the values of the pixels in that box (although you can use the PSF to inform what your box should be for a point object). Commented Sep 21, 2018 at 17:41
• When you say draw a box around the object, how does one choose the correct size of box? With a nonzero background and noise in the image, a bigger box would increase the flux even though the source doesn't extend further into the bigger box. Commented Sep 21, 2018 at 17:56
• Also, sorry but how is an asteroid an extended source? (this term was unfamiliar to me until I looked it up just now) Wouldn't the angular size of an asteroid imaged by WISE be less than the telescope resolution, classifying it as a point source? Commented Sep 21, 2018 at 18:00

You can't really get a PSF for an extended object. The Point-Spread-Function is a parameter of your imaging system which describes the shape of the image of a putative point source.

I don't have any knowledge of what your python package is intended to do, so I can't comment directly on the results. I will point (sorry) out that the geometric centroid of the image is not generally the same as the energy(brightness) centroid of the image.

Further, the total flux is just the sum of all the pixels covered by the asteroid's image (less background noise and whatnot). It's not clear why you need the PSF if you are just doing energy/power calculations.

• Well, I am using the PSF for other purposes as well, but one of the things I would like to do is find the total flux of the source image. I was under the impression I could accomplish this with the help of the PSF, but I haven't been able to find anything online about how exactly it can do this. Do you know how the normalization of the PSF might make a difference in the flux of the source image? Also, yes, the centroid coordinates I cited in my question do correspond to the max intensity, not just the geometrical center. Commented Sep 20, 2018 at 20:38

The image is 2D pixels, the PSF fitting routine would do something like this. First, you pick the PSF function, which is a 2D gaussian in this case. Then, the routine will choose a centroid, mu, and variance (or mus and variances).

A couple of things that will complicate the routine. i) Local peak due to a noise. This will mess up the centroid. Fixing is to do the centroid search by using averages given some sizes of the areas. Also, you might want to specify how large the search area should be given the initially specified spot. Or, if you are certain about the location of centroid, you can fix it and the centroid would not be a choice variable. ii) Aperture size. If the object is a perfect 2D gaussian point source, the aperture size does not matter. But with noise and other surrounding objects, you have to judge what should be the optimal size: too large = more noise and including other objects, too small = incorrect fit shape due to noise dominating. There might be something more to think about that I cannot recall at this moment.

Speaking of this, the fitting routine can fit even an extended source, but you won't get a good result unless you use a proper kernel (i.e., not gaussian).

I personally don't know about astropy. I used IRAF, which has its equivalent as PyRAF in Python, I believe. You might want to check that as well.

Your problem about incorrect centroid, I believe, roots from the local spikes from noise and the search routine not having enough size for a good average. This problem directly leads >1 total flux as you mentioned as well.