<p>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.  </p>

<p>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.</p>

<p>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. </p>

CC BY-SA 4.0
2018-09-20T20:38:01.310Z
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.
curious_cosmo
CC BY-SA 4.0
2018-09-20T19:35:59.080Z
2018-09-20T19:35:59.080Z
Carl Witthoft
<p>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).</p>

<p>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.</p>

<p>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).</p>

<p>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.</p>

<p>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.</p>

CC BY-SA 4.0
2018-09-23T15:20:24.017Z
2018-09-23T15:20:24.017Z
Kornpob Bhirombhakdi
<p>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. </p>

<p>I know that astropy offers a package called <code>photutils.psf</code> 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 <code>model_data</code> is the 2d numpy array representation of the PSF, and <code>img</code> is the asteroid image data):</p>

<pre><code>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))
</code></pre>

<p>The output I get is then:</p>

<pre><code>flux       x_0       y_0

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

9282.784  466.44    531.69
</code></pre>

<p><code>2705.82 25301.56 980812434.98</code></p>

<p>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.</p>

<p>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 <code>x_0</code> and <code>y_0</code> 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?</p>

CC BY-SA 4.0
2018-09-21T14:36:30.260Z
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.
zephyr
CC BY-SA 4.0
2018-09-21T16:41:53.353Z
@zephyr Yes, but how is this accomplished by fitting a PSF to the image?
curious_cosmo
CC BY-SA 4.0
2018-09-21T17:41:15.163Z
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).
zephyr
CC BY-SA 4.0
2018-09-21T17:56:06.540Z
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.
curious_cosmo
CC BY-SA 4.0
2018-09-21T18:00:40.773Z
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?
curious_cosmo
CC BY-SA 4.0
2018-09-21T18:06:15.427Z
[Presumably you've already corrected your images for bias, flat field, etc](https://www.eso.org/~ohainaut/ccd/). so drawing a bigger box shouldn't give you more flux (at least not a meaningful amount more). Generally only stars count as point sources, however, your asteroid in your case *may* be close enough to a point source, but either way, its irrelevant for the problem you're trying to solve. Although we're entering discussion territory and comments aren't meant for that. Please feel free to ask additional, new questions on this stack exchange to get more clarification.
zephyr
<data-analysis><python><astropy>
Fitting PSF to an image to get total flux (Python)
curious_cosmo
# Fitting PSF to an image to get total flux (Python)

(Asked by: username_0 on 2018-09-20)

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?

**Tags:** data-analysis, python, astropy
> **Comments:**
> **username_1**: 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.  
> **username_0**: @zephyr Yes, but how is this accomplished by fitting a PSF to the image?  
> **username_1**: 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).  
> **username_0**: 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.  
> **username_0**: 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?  
> **username_1**: [Presumably you've already corrected your images for bias, flat field, etc](https://www.eso.org/~ohainaut/ccd/). so drawing a bigger box shouldn't give you more flux (at least not a meaningful amount more). Generally only stars count as point sources, however, your asteroid in your case *may* be close enough to a point source, but either way, its irrelevant for the problem you're trying to solve. Although we're entering discussion territory and comments aren't meant for that. Please feel free to ask additional, new questions on this stack exchange to get more clarification.  

---

## 2 answers

### Answer by username_2 on 2018-09-20. Score: 2

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.

> **Comments:**
> **username_0**: 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.  

---

### Answer by username_3 on 2018-09-23. Score: 1

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.


---

