Unrealistic values with --force-beam and --vectorize

[tmillenaar]

When both the flags --vectorize and --force-beam are used in conjunction, some sources get weird values like a ra_err and dec_err of 0, or very high fluxes.

I tried to reproduce this with the existing test images in the pyse repo but I did not manage. I did find some problematic configuration in some test files from Antonia. The following example uses one of these files.

To reproduce, download the following image: https://www.dropbox.com/scl/fo/dv8iwdqaeiwc3ddflw7wx/AMkE9jRzlbQIDBJ8MTi1JYQ/Images/GRB201006A?dl=0&preview=GRB201006A_final_2min_srcs-t0021-image-pb.fits&rlkey=ws94p3kytmefiujv2cu2nhco9&subfolder_nav_tracking=1

Then run the following command with pyse version v0.6.0:
pyse --csv GRB201006A/GRB201006A_final_2min_srcs-t0021-image-pb.fits --force-beam --vectorized --back-size-x 50 --back-size-y 50

Then when I inspect the results using head GRB201006A_final_2min_srcs-t0021-image-pb.csv -n 20 I get:

Note that when I remove either the --vectorize or the --force-beam the results seem fine. It seems to be the combination that causes issues.

This issue relates to a comment made in a different issue: #177 (comment)

[HannoSpreeuw]

I think I found the source of the problem.

To derive the peak brightness, tweaked moments takes the pixel value at the rounded barycenter position and extrapolates that to the barycenter position.
That extrapolation can be based on the (forced) beam or based on the measured source shape. It is an extrapolation over a subtended angle at the subpixel level.

In some edge cases the rounded barycenter positon is not an island pixel.
In that case, tweaked moments is still applied, but the maximum pixel value and the position of that pixel are used to extrapolate from. This appeared to work well.

All the problems seem to occur when the island of pixels is relatively wide, when the rounded barycenter position is not an island pixel and the maximum pixel position is far from the island barycenter. That is rare, but can occur with an extended source with a relatively flat brightness distribution. The extrapolation then leads to incredibly large derived peak brightness, but likely only for forced beam measurements: for free shape the measured axes will be large such that the extrapolation will not lead to large values.

I guess the answer should be that we should not apply "tweaked moments" when the rounded barycenter position does not correspond to an island pixel and use fudge_max_pix times the maximum pixel value instead.

@AntoniaR thanks for catching this, this is a serious error, which was not caught by our extensive regression tests with up to 167,000 artificial sources. But these sources were "well formed", such that the error did not occur.

Perhaps I could include this image in a unit test, but I am not sure how to implement that without ground truth data.
Well, the unit test should perhaps assert that no derived peak brightness should be, say, 20% higher than the maximum pixel value over the entire image, but that would be a pretty coarse test.

Some data for my future reference, selecting island data for derived peak brightnesses > 100 Jy/beam, from this command: pyse GRB201006A_final_2min_srcs-t0021-image-pb.fits --force-beam --vectorized --back-size-x 50 --back-size-y 50:

Number of detected sources = 1038

peak =  1054.2947525880033  , maxi =  0.030075108632445335  , basevalue =  0.030075108632445335  , mask.any() =  False  , min_width =  16  , deltax =  3.0591402833020833  , deltay =  7.749962332556443  , smin =  1.8212668942844181  , smaj =  2.4119618032265238  , theta =  0.25895872745169685 


peak =  148826.42171051877  , maxi =  0.047940585762262344  , basevalue =  0.047940585762262344  , mask.any() =  False  , min_width =  11  , deltax =  -8.526144236439507  , deltay =  -0.4288498322058176  , smin =  1.8212668942844181  , smaj =  2.4119618032265238  , theta =  0.25895872745169685 


peak =  112.99034466457526  , maxi =  0.09484915435314178  , basevalue =  0.09484915435314178  , mask.any() =  False  , min_width =  14  , deltax =  0.7388329380896614  , deltay =  7.336954235600766  , smin =  1.8212668942844181  , smaj =  2.4119618032265238  , theta =  0.25895872745169685 


peak =  8069.562513688345  , maxi =  0.030813246965408325  , basevalue =  0.030813246965408325  , mask.any() =  False  , min_width =  5  , deltax =  -6.501387096276856  , deltay =  -4.541582242520785  , smin =  1.8212668942844181  , smaj =  2.4119618032265238  , theta =  0.25895872745169685 


peak =  57032.37147844775  , maxi =  0.04439084231853485  , basevalue =  0.04439084231853485  , mask.any() =  False  , min_width =  5  , deltax =  8.401070859081242  , deltay =  -1.1798297518388456  , smin =  1.8212668942844181  , smaj =  2.4119618032265238  , theta =  0.25895872745169685 


peak =  34119236030.69715  , maxi =  0.03593478351831436  , basevalue =  0.03593478351831436  , mask.any() =  False  , min_width =  15  , deltax =  -11.704237973690063  , deltay =  3.633650425578578  , smin =  1.8212668942844181  , smaj =  2.4119618032265238  , theta =  0.25895872745169685 


peak =  4.799304561372789e+17  , maxi =  0.040669266134500504  , basevalue =  0.040669266134500504  , mask.any() =  False  , min_width =  10  , deltax =  -14.637748757572489  , deltay =  5.764220347858702  , smin =  1.8212668942844181  , smaj =  2.4119618032265238  , theta =  0.25895872745169685 

[AntoniaR]

Is this now in the latest version when I update with pip install?

[HannoSpreeuw]

Coming up....I'd like to include #181 as well in a new release.
As soon as that MR is settled, I will release.
So within hours or within minutes.

[AntoniaR]

Ah - perfect! I was just querying that issue 🙂