Skip to content

Fix effective beam calculation in wavelets - #550

Open
AlexKurek wants to merge 9 commits into
lofar-astron:masterfrom
AlexKurek:Fix-effective-beam-calculation-in-wavelets
Open

AlexKurek wants to merge 9 commits into
lofar-astron:masterfrom
AlexKurek:Fix-effective-beam-calculation-in-wavelets

Conversation

@AlexKurek

@AlexKurek AlexKurek commented Nov 3, 2026 •

Copy link
Copy Markdown
Contributor

Replace the previous filter-width approximation with a calculation based on the actual filter variance and cumulative smoothing across scales. The effective beam is calculated in pixel coordinates and converted back using img.pix2beam(), preserving the original PA.

Test:
python test_atrous_effective_beam.py --bmaj-pix 4 --bmin-pix 4 --pa-deg 0

akurek@lof8:~$ python test.py 

Filter b3: vector=[0.0625, 0.25, 0.375, 0.25, 0.0625], variance per axis=1 pix^2
 scale | old relative error | new relative error | error reduction
-------+-------------------+-------------------+---------------
     1 |      7.209675e+00 |      4.583558e-15 | 100.00000000% OK
     2 |      8.271852e+00 |      1.959275e-15 | 100.00000000% OK
     3 |      8.350022e+00 |      6.938412e-16 | 100.00000000% OK
     4 |      8.106143e+00 |      2.379599e-16 | 100.00000000% OK
     5 |      7.905081e+00 |      2.622428e-16 | 100.00000000% OK

Filter tr: vector=[0.25, 0.5, 0.25], variance per axis=0.5 pix^2
 scale | old relative error | new relative error | error reduction
-------+-------------------+-------------------+---------------
     1 |      4.581452e+00 |      5.405825e-15 | 100.00000000% OK
     2 |      4.671865e+00 |      3.121036e-15 | 100.00000000% OK
     3 |      4.286400e+00 |      1.354163e-15 | 100.00000000% OK
     4 |      3.859504e+00 |      3.615103e-16 | 100.00000000% OK
     5 |      3.599391e+00 |      1.664182e-16 | 100.00000000% OK

PASS: the new calculation has a smaller covariance error than the old one at every tested scale.

test_atrous_effective_beam.py

30 largest sources before and after. I think its graphically little different but not conclusive. But still - algorithmically better.
comapr_wavelets_newParams__PR550.pdf

@gmloose
gmloose requested a review from darafferty November 5, 2026 14:03
@gmloose

gmloose commented Nov 5, 2026

Copy link
Copy Markdown
Collaborator

@darafferty, I lack the domain knowledge to judge if this PR is eligible for merging. Can you review this one, please?

@gmloose

gmloose commented Nov 5, 2026

Copy link
Copy Markdown
Collaborator

@AlexKurek, the current branch has conflicts. Please have a look.

@darafferty darafferty left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the fixes are correct (although I did not write the original code and am not an expert on it either). Nice to see the comparison of the old and new fits, though I think many of the sources shown there are not really extended so perhaps not very useful... But I do see in my tests that a number of complex sources are better fit with the fixes, which gives me confidence they are correct.

@gmloose gmloose left a comment •

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I let Claude have a look at this PR, because both David and I don't feel comfortable enough to give it a thorough review. Feel free to provide feedback.

pr-550-review.md

Thanks for this — the variance-based derivation is sound: the a-trous kernel at scale k is the base kernel dilated by 2^(k-1), so cumulative variance is var · (4^j − 1)/3, which matches atrous(). The resulting wavelet beam is roughly 2.5–3× smaller than the old formula, so this is a substantial behaviour change. A few points:

1. Reference outputs not updated (blocking)
The branch predates #576 and does not regenerate test/reference_outputs/wavelet. Since tbdsf_process_image.py compares wavelet images (check_images=True), and the new beam changes rms_box_bright (rmsimage.py:189), initial fit guesses (gausfit.py:357) and flagging limits, the regression test should fail once merged with master. Please rebase and regenerate the references.

2. Beam models c_j, but the fitted plane is w_j = c_{j−1} − c_j (wavelet_atrous.py:166)
filter_var_total is the variance of the smoothed image c_j. With the default atrous_sum=False, the fitted image is the band-pass difference c_{j−1} − c_j, whose sources are narrower than that. So the assumed beam is systematically too large there; only atrous_sum=True matches the derivation. Could you either document this as a deliberate approximation or use something in between (e.g. geometric mean of the c_{j−1} and c_j variances)?

3. Non-square pixels (wavelet_atrous.py:168)
An isotropic pixel-space variance is added to the beam's major/minor sigmas, which beam2pix converted with cdelt1/cdelt2 respectively. If |cdelt1| ≠ |cdelt2| and the beam is rotated, the kernel's sky-frame variance along the beam axes depends on bpa, so the result is not the true convolved beam. Probably minor in practice, but worth a comment.

4. Position angle changed from 0° to original bpa (wavelet_atrous.py:172)
Previously the wavelet beam PA was hard-coded to 0.0; now the original PA is propagated. Arguably more correct, but it changes fit results — worth mentioning in the PR description / changelog.

5. wimg.pixel_beam() still returns the unsmoothed beam (wavelet_atrous.py:179)
wimg.pixel_beam / pixel_beamarea stay bound to the original image's closures, so consumers using pixel_beam() — rmsimage.py:277 (min_size_allowed), gausfit.py:1149 (default ini guess for j ≤ 3), gaul2srl.py:337 — don't see the corrected beam. This is pre-existing (and partly intentional per the comment in gausfit.py:354), but it means the fix only partially takes effect.

Minor / cleanup

  • wavelet_atrous.py:75: filter_var is recomputed from the vector and normalised by N.sum(filter_vec), which is always 1. Storing 'var': 0.5 / 'var': 1.0 in the filter dict next to size/vec would be simpler and self-documenting.
  • wavelet_atrous.py:73 / y1 (~line 140): rms_box sizing still uses the old filter-width heuristic l + (l−1)(2^(j−1)−1), now inconsistent with the variance-based beam. Also, float( N.sum(... has a stray space.

@AlexKurek

AlexKurek commented Nov 10, 2026 •

Copy link
Copy Markdown
Contributor Author

I think of these, the following are important, so applied them in 514b9e2. Also updated the references.

2
The cumulative filter variance describes the smoothed image (c_j), but not the band-pass image (w_j = c_{j-1} - c_j).
Used the geometric mean of the total variances at scales (j-1) and (j) as an approximation when atrous_sum=False, while retaining the cumulative variance when atrous_sum=True.

Second form "Minor / cleanup"
Replaced y1 filter-widthwith the major-axis FWHM of the effective beam in pixel coordinates. This will make the RMS box sizing consistent with the beam model used for Gaussian fitting, while retaining the existing factor of 5 and the original RMS box size as a lower bound.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants