Fitting width of Moffat2D and Gaussian2D models with PSFPhotometry

[iant321]

Hi,

I wrote a code which uses IntegratedGaussianPRF and PSFPhotometry. Works ok, but the fit is not optimal as my measured stellar profiles are non-Gaussian and are more like Moffat. So I tried to implement the Moffat2D model. It works, but it seems that only x, y, and flux (actually amplitude) are fitted and not gamma or alpha. Same with Gaussian2D - x_stddev and y_stddev are not fitted. All these parameters are set to be free in my models. With IntegratedGaussianPRF, if I set

psf_model = IntegratedGaussianPRF(flux=1, sigma=fwhm/2.355) psf_model.sigma.fixed = False

sigma of the model is fitted ok. (Yes, I know that this is not good to allow sigma to vary for all sources as the results will be biased for low S/N, what I try to do is to measure sigma for several bright sources and then use this fixed value for the final fit).

My astropy/photutils versions are
astropy 5.3.4
photutils 1.10.0

I was unable to find working examples of how to use Gaussian2D (or Moffat2D) with PSFPhotometry and would appreciate any help.

For some strange reason a working example which I tried to insert here with "code" formatting was improperly formatted. So I provide a link to the code. The code can simulate either Moffat or Gaussian stars and then fit them with either Moffat2D or Gaussian2D models. If the initial widths (gamma or x_stddev, y_stddev) for the fit models are wrong, PSFPhotometry does not adjust them.

I attach a screenshot with the results of the above code. Simulated are Gaussian stars with fwhm=3.0, fitted with initial fwhm=6.0.

[larrybradley]

I've modified your code to show an example of how to fit a Moffat2D PSF model, allowing gamma and alpha to vary. Here's an example notebook:

https://gist.github.com/larrybradley/263f2c9c1a47013a38cd71c6690e8105

However, there appears to be bug when the PSF model fit fails (ValueError: setting an array element with a sequence.). Thus, I had to make the initial model parameters closer to the true parameters for the fit to converge.

There are also a few things to note. When you input an astropy model into make_psf_model, the output model is always a CompoundModel that "wraps" your input model. The wrapping maps (or adds) the appropriate x, y, and flux parameters needed for PSF fitting. So to "unfix" any parameters, you need to use their names in the compound model (e.g., psf_model.gamma_2.fixed = False).

Also, "amplitude" is not the same as "flux". Using it will give the wrong fit "flux" output values unless the input model was already normalized. Instead of mapping "flux" to "amplitude" just let make_psf_model do the normalization and add a scaling parameter that will be the model flux.

Third, one important thing to keep in mind is that when fitting with the models output by make_psf_model, the models are evaluated by taking the sample value at pixel centers (this is how astropy modeling works), not integrated over pixel bins. So a Gaussian2D PSF model will give different results than IntegratedGaussianPRF, which does integrate over pixel bins. There is a PRFAdapter class that attempts to take an astropy model (e.g., output of make_psf_model) and fit with it by integrating over pixel bins. However, it is extremely slow (because of the integrations). It will be updated in the next release to improve performance by doing approximate integrations.

I'm also planning to add a section in the documentation with these notes and an example.

[iant321]

Dear Larry,

thanks a lot for your detailed explanations and the code. The solution is very simple, yet it is very hard to find it without having comprehensive knowledge of various photutils components and their interactions.

A few comments:

1. I am aware of the bug as I encountered it myself. Just did not want to mix this into my question. I was going to ask it next. Good to know that the issue is known and will hopefully be resolved.

2. Yes, I know about the difference between amplitude and flux. I am just not aware of how to "add scaling parameter" . Again, I postponed this for the future to not complicate things in my original post.

3. I know about the difference between IntegratedGaussianPRF and Gaussian2D but thank you for your explanation. I very much appreciate your good will! Still, with my images I expect to have better results with Moffat. It would be great if integration over pixels would be implemented in the future.

I will now wait for a day for your possible comment on my current post and will close the issue as "solved". Thanks again!

[larrybradley]

Yes, 1 is a known bug -- the compound PSF model gets mixed up with possible PSF model grouping (which is also a compound model). For this reason, I suspect using grouping with make_psf_model output will also break.

For 2, simply do not use the flux_name='amplitude' when you call make_psf_model (see my example notebook). If flux_name is not given, make_psf_model will normalize the model and add a scaling parameter corresponding to flux.

For 3, you could try PRFAdapter on a few sources, but it is too slow to use on many sources. I have a plan to improve performance.

[iant321]

Oh, I see. If I set flux_name=None or simply omit it, the model will compute the flux based on ampl, gamma, alpha. Did not realize it was that simple. As an off-topic note - for a person like me who came from traditional C programming, it is sometimes difficult to readjust to python style. I tend to think that things are more complicated (at the user level) than they are :) Thank you again for your help. I will now close the topic.

[iant321]

Closed as solved.

[larrybradley]

@iant321 I fixed the bug in #1672.

[iant321]

I re-opened the issue just to make a comment. The updated version kind of works, but for real images which are 4000x4000 pixels, running it results in the error message (sorry the code formatiing does not work again):

Traceback (most recent call last):
File "./test.py", line 61, in phot = psfphot(data, error=error)
^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/local/lib64/python3.12/site-packages/photutils/psf/photometry.py", line 979, in call
fit_models = self._fit_sources(data, init_params, error=error,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/local/lib64/python3.12/site-packages/photutils/psf/photometry.py", line 708, in _fit_sources
fit_models = self._make_fit_results(fit_models, fit_infos)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/local/lib64/python3.12/site-packages/photutils/psf/photometry.py", line 633, in _make_fit_results
fit_param_errs = np.array(self._order_by_id(fit_param_errs))
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (12,) + inhomogeneous part.

How to reproduce - just run your notebook version with the image size set to 4000x4000. If I comment the lines

psf_model.gamma_2.fixed = False
psf_model.alpha_2.fixed = False

the code works, but of course does not fit these parameters. I understand that the error message is about some arrays having inconsistent lengths, but other than that I am unable to say anything.

[larrybradley]

Many thanks for letting me know! I'll take a look soon.

[larrybradley]

@iant321 This appears to be the bug I fixed in #1672. I cannot reproduce it using the latest dev version (from main). Can you please retry with the latest development version?

pip install "photutils[all] @ git+https://github.com/astropy/photutils.git"

[iant321]

This version of the code works! I have not tested it thoroughly, but running the code on a few 4000x4000 frames worked without any flaws and fitted both gamma and alpha. One strange thing, though. I removed the "flux_name=" mapping from the call to make_psf_model, so PSFPhotometry should compute flux based on the fit values of gamma, alpha, and amplitude. It does, but the value of the flux is ~2 times smaller than the value found with e.g. Gaussian PSF or measured directly within some appropriate aperture (for isolated stars). Yet, the stellar profile is reproduced well. The ratio is approximately the same for various sources so I suspect that the coefficient "2" is lost somewhere when the flux is computed.

Hmm, I run the minimal working example for 100x100 frame, the ratio for most simulated stars is about 2, but for some is different. So it is not as simple as I suggested.

Also, should I set

psf_model.amplitude_2.fixed = False

or not? The resulting fluxes in the two cases (set/omitted from the code) are different. Yet, if I e.g. increase the flux of all simulated sources and omit this line from the code, the output fluxes change. I am a bit lost here.

Thanks a lot for you work!