astropy
Why averaging for high harmonics in build_ellipse_model
photutils version: 1.12.0 In file isophote/model.py, Line 124, Col 64
There is a "/ 4.0" in the calculation of harm (contribution to isophote intensity from high harmonics) that takes the average of the term $A_3\sin(3\phi)+B_3\cos(3\phi)+A_4\sin(4\phi)+B_4\cos(4\phi)$. But according to the formula, that term is directly added into the intensity.
Even though the isophote intensity is weighted and added to the overlapping four pixels, the total intensity is intens + harm and seems not to need any averaging for harm. So I wonder why it gets averaged.
Thank you very much!
@asterg1122 I don't know why it's averaged. This code was ported to Python from the IRAF ellipse package by the original author of the IRAF code. Unfortunately, he has since retired. You could take a look at the original IRAF code and see if it's also averaged there: https://github.com/spacetelescope/stsdas_stripped/tree/043c173fd5497c18c2b1bfe8bcff65180bca3996/stsdas/pkg/analysis/isophote/src
Please report back if you find a discrepancy!
@larrybradley Thank you for your reply! I have read the source code and find corresponding lines in the file t_model.x
IRAF code:
photutils:
In lines 272-273 of IRAF code, the aux (corresponds to harm in photutils) is indeed averaged.
BUT there is a discrepancy in the following lines. In lines 274-277 of IRAF code aux is equally added to four pixels, while in photutils harm is weighted. Therefore, averaging is correct in the original IRAF code, but if harmonic item is weighted to be added up, it should not be averaged.
@asterg1122 Many thanks for reporting back! To match the original IRAF code, it appears that an unweighted (but averaged) harm should be added in lines 138-141. Are you suggesting instead to keep lines 138-141 as is (adding a weighted harm), but change lines 120-124 to remove to averaging?
How did you discover this? Do you have an example where model produced an incorrect output image? It would be nice to add a simple test case to the code base.
@larrybradley I found this just when I was learning the source code of build_ellipse_model. Sorry I do not have a very practical case. Actually I only test one galaxy in simulation and the difference in residual between averaging and no averaging is not significant, since the harmonic coefficients themselves are small.
I have tried to construct a quite preliminary example, see ExampleForIssue1776.txt (please change the suffix to .ipynb since this type is not supported to be uploaded...)
From this example, it seems better to remove the averaging in lines 120-124. Compared with the original function, the multipole structure in residual disappears. Adding unweighted harm in lines 138-141 leads to a large residual in one pixel at the edge, as the weight in that pixel is too small. I hope it can be helpful to build a test case.
BTW, I think Issue #1722 has mentioned a similar problem.
@asterg1122 Many thanks again for reporting this. Your investigation and example notebook was extremely helpful. From my additional testing, I agree that the approach of including weighted (but not averaged) harmonics gives the best overall results and the smallest residuals. The other approach of averaging, but then adding them unweighted results in odd high-frequency variations in the model.
I fixed this in #1810, and it will appear in the next release.