From be41ba47ae21c5c938a68cbfd1ec9860750cbe5d Mon Sep 17 00:00:00 2001 From: HannoSpreeuw Date: Fri, 9 Jan 2026 13:09:35 +0100 Subject: [PATCH 1/2] Fix 188 1) For "measure.moments" do not revert to the maximum pixel position when the rounded barycenter position is not part of the island. This led to erroneous "tweaked moments" extrapolations for "ill-formed" islands, e.g. islands that are actually artefacts of poor calibration and/or deconvolution, such as side lobes of bright sources. 2) Reflect this change in a comment in "measure.moments". 3) In "measure.moments_enhanced" 1) had already been effectuated through commit 16eed05. But that commit had not yet been reflected in two of its comments. 4) It seems more robust to split the conditional "if np.sign(threshold) == np.sign(basevalue) and in_source_island" in "measure.moments_enhanced" since "basevalue" may be undefined for "ill-formed" islands, i.e. it seems safer to impose the"if in_source_island" conditional upfront. --- sourcefinder/measure.py | 181 +++++++++++++++++++++------------------- 1 file changed, 94 insertions(+), 87 deletions(-) diff --git a/sourcefinder/measure.py b/sourcefinder/measure.py index 3ea4141..7ad37b6 100644 --- a/sourcefinder/measure.py +++ b/sourcefinder/measure.py @@ -165,57 +165,61 @@ def moments(data, fudge_max_pix_factor, beam, beamsize, threshold=0): rounded_barycenter = int(round(xbar)), int(round(ybar)) # basevalue and basepos will be needed for "tweaked moments". + in_source_island = False try: if not data.mask[rounded_barycenter]: - basepos = rounded_barycenter - else: - basepos = maxpos + in_source_island = True except IndexError: - basepos = maxpos + # rounded_barycenter is not part of the island. + pass except AttributeError: # If the island is not masked at all, we can safely set basepos to # the rounded barycenter position. + in_source_island = True + + if in_source_island: basepos = rounded_barycenter - basevalue = data[basepos] - - if np.sign(threshold) == np.sign(basevalue): - # Implementation of "tweaked moments", equation 2.67 from - # Spreeuw's thesis. In that formula the "base position" was the - # maximum pixel position, though. Here that is the rounded - # barycenter position, unless it's masked. If it's masked, it - # will be the maximum pixel position. - deltax, deltay = xbar - basepos[0], ybar - basepos[1] + basevalue = data[basepos] - epsilon = np.log(2.0) * ( - (np.cos(theta) * deltax + np.sin(theta) * deltay) ** 2 - + (np.cos(theta) * deltay - np.sin(theta) * deltax) ** 2 - * semiminor_tmp - / semimajor_tmp - ) + if np.sign(threshold) == np.sign(basevalue): + # Implementation of "tweaked moments", equation 2.67 from + # Spreeuw's thesis. In that formula the "base position" was the + # maximum pixel position, though. Here that is the rounded + # barycenter position, unless it's masked or not part of the island. + # In those cases no "tweaked moments" extrapolation will be + # applied, and we revert to data.max() * fudge_max_pix. + deltax, deltay = xbar - basepos[0], ybar - basepos[1] + + epsilon = np.log(2.0) * ( + (np.cos(theta) * deltax + np.sin(theta) * deltay) ** 2 + + (np.cos(theta) * deltay - np.sin(theta) * deltax) ** 2 + * semiminor_tmp + / semimajor_tmp + ) - # Set limits for the root finder similar to the bounds for - # Gaussian fits. - if basevalue > 0: - low_bound = 0.5 * basevalue - upp_bound = 1.5 * basevalue - else: - low_bound = 1.5 * basevalue - upp_bound = 0.5 * basevalue - - # The number of iterations used for the root finder is also - # returned, but not used here. - peak, _ = newton_raphson_root_finder( - find_true_peak, - basevalue, - low_bound, - upp_bound, - 1e-8, - 100, - threshold, - epsilon, - semiminor_tmp, - basevalue, - ) + # Set limits for the root finder similar to the bounds for + # Gaussian fits. + if basevalue > 0: + low_bound = 0.5 * basevalue + upp_bound = 1.5 * basevalue + else: + low_bound = 1.5 * basevalue + upp_bound = 0.5 * basevalue + + # The number of iterations used for the root finder is also + # returned, but not used here. + peak, _ = newton_raphson_root_finder( + find_true_peak, + basevalue, + low_bound, + upp_bound, + 1e-8, + 100, + threshold, + epsilon, + semiminor_tmp, + basevalue, + ) if np.sign(threshold) == np.sign(peak): # The corrections below for the semi-major and semi-minor axes are @@ -500,6 +504,7 @@ def moments_enhanced( rounded_barycenter = int(round(xbar)), int(round(ybar)) mask = (posx == rounded_barycenter[0]) & (posy == rounded_barycenter[1]) in_source_island = mask.any() + if in_source_island: # In this case the rounded barycenter position is in source_island. # Note that mask, from the way that posx and posy have been constructed, @@ -520,8 +525,9 @@ def moments_enhanced( # Implementation of "tweaked moments", equation 2.66 from # Spreeuw's thesis. In that formula the "base position" was the # maximum pixel position, though. Here that is the rounded - # barycenter position, unless it's masked. If it's masked, it - # will be the maximum pixel position. + # barycenter position, unless it's not part of the island. In that + # case no "tweaked moments" extrapolation will be applied, and we + # revert to source_island.max() * fudge_max_pix. exponent = np.log(2.0) * ( ((np.cos(theta) * deltax + np.sin(theta) * deltay) / smin) ** 2 + ((np.cos(theta) * deltay - np.sin(theta) * deltax) / smaj) @@ -592,50 +598,51 @@ def moments_enhanced( else: theta -= math.pi / 2.0 - if np.sign(threshold) == np.sign(basevalue) and in_source_island: - # Implementation of "tweaked moments", equation 2.67 from - # Spreeuw's thesis. In that formula the "base position" was the - # maximum pixel position. Here that is the rounded - # barycenter position, unless it's masked. If it's masked, it - # will be the maximum pixel position. - # The rounded barycenter position will buy us slightly more - # stability (established by testing) than the maximum pixel - # position. - # The condition is a sanity check for maps other than Stokes I, - # where the threshold can be negative. For Stokes I maps, which - # will have positive thresholds, at least for all - # non-pathological cases, this condition will always be true. - - epsilon = np.log(2.0) * ( - (np.cos(theta) * deltax + np.sin(theta) * deltay) ** 2 - + (np.cos(theta) * deltay - np.sin(theta) * deltax) ** 2 - * semiminor_tmp - / semimajor_tmp - ) + if in_source_island: + if np.sign(threshold) == np.sign(basevalue): + # Implementation of "tweaked moments", equation 2.67 from + # Spreeuw's thesis. In that formula the "base position" was the + # maximum pixel position. Here that is the rounded + # barycenter position, unless it's not part of the island. In + # that case no "tweaked moments" extrapolation will be + # applied, and we revert to source_island.max() * + # fudge_max_pix. + # The condition is a sanity check for maps other than Stokes I, + # where the threshold can be negative. For Stokes I maps, which + # will have positive thresholds, at least for all + # non-pathological cases, this condition will always be true. + + epsilon = np.log(2.0) * ( + (np.cos(theta) * deltax + np.sin(theta) * deltay) ** 2 + + (np.cos(theta) * deltay - np.sin(theta) * deltax) + ** 2 + * semiminor_tmp + / semimajor_tmp + ) - # Set limits for the root finder similar to the bounds for - # Gaussian fits. - if basevalue > 0: - low_bound = 0.5 * basevalue - upp_bound = 2.0 * basevalue - else: - low_bound = 2.0 * basevalue - upp_bound = 0.5 * basevalue - - # The number of iterations used for the root finder is also - # returned, but not used here. - peak, _ = newton_raphson_root_finder( - find_true_peak, - basevalue, - low_bound, - upp_bound, - 1e-8, - 100, - threshold, - epsilon, - semiminor_tmp, - basevalue, - ) + # Set limits for the root finder similar to the bounds for + # Gaussian fits. + if basevalue > 0: + low_bound = 0.5 * basevalue + upp_bound = 2.0 * basevalue + else: + low_bound = 2.0 * basevalue + upp_bound = 0.5 * basevalue + + # The number of iterations used for the root finder is also + # returned, but not used here. + peak, _ = newton_raphson_root_finder( + find_true_peak, + basevalue, + low_bound, + upp_bound, + 1e-8, + 100, + threshold, + epsilon, + semiminor_tmp, + basevalue, + ) if np.sign(threshold) == np.sign(peak) and abs(peak) > abs( threshold From 8aeeecb1a4a9f53430b04392781ae97c5cb1c624 Mon Sep 17 00:00:00 2001 From: HannoSpreeuw Date: Fri, 9 Jan 2026 13:45:58 +0100 Subject: [PATCH 2/2] Set uniform rootfinder bounds. Set rootfinder bounds in "measure.moments" in the same way as "measure.moments_enhanced". From a few tests we inferred that any "tweaked moments" extrapolation at the subpixel level will never reach "2.0 * basevalue" while it may hit "1.5 * basevalue" in some cases. This means that "upp_bound = 2.0 * basevalue" should be high enough to return a root that is never determined by the upper bound. This suggests that "upp_bound = np.inf" could be used, but doing so may introduce unnecessary risks. --- sourcefinder/measure.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sourcefinder/measure.py b/sourcefinder/measure.py index 7ad37b6..575a629 100644 --- a/sourcefinder/measure.py +++ b/sourcefinder/measure.py @@ -201,9 +201,9 @@ def moments(data, fudge_max_pix_factor, beam, beamsize, threshold=0): # Gaussian fits. if basevalue > 0: low_bound = 0.5 * basevalue - upp_bound = 1.5 * basevalue + upp_bound = 2.0 * basevalue else: - low_bound = 1.5 * basevalue + low_bound = 2.0 * basevalue upp_bound = 0.5 * basevalue # The number of iterations used for the root finder is also