From 6be94f2be0378a09532ff88e61ebe0435bb46d03 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 29 Jul 2026 09:36:40 +0200 Subject: [PATCH] Merging: do not subtract a negative intensity's Poisson term The expected-variance weights decompose an observation's sigma^2 into a background part and a Poisson signal part, then rebuild the signal part at the reflection's merged mean. The decomposition subtracted corr*I with I taken as-is, so a negative I ADDED to the background part: an observation at I = -1.5 with sigma^2 = 1 came out with a base variance of 2.7 rather than 1. That inflates the variance of precisely the down-fluctuated observations the correction exists for. Below about one photon they are then under-weighted and the merged mean is biased high - the same direction of error, in the same regime, that weighting by the observation's own sigma produces. Subtract max(0, I) instead: a negative intensity has no Poisson signal to remove. Both users of the decomposition are fixed - the stills merge, where expected-variance weighting is now the default, and the rotation combine it was mirrored from, which had it first. Co-Authored-By: Claude Opus 5 (1M context) --- image_analysis/scale_merge/Merge.cpp | 7 ++++++- image_analysis/scale_merge/RotationScaleMerge.cpp | 4 +++- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 52c31ed5..ae5cba3d 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -147,9 +147,14 @@ float MergeOnTheFly::CorrectedSigma(float I_corr, float sigma_corr, float image_ // and re-add corr* - so the weight no longer correlates with the observation's own fluctuation. This // mirrors the rotation combine in RotationScaleMerge::process_rawrun and is bit-identical when the // observation sits at its reflection mean. + // The signal part removed here is the one the observation's own sigma actually carries, so it is + // taken at max(0, I): a negative intensity has no Poisson signal to subtract, and subtracting it + // anyway ADDS to the background part and inflates the weight of exactly the down-fluctuated + // observations this correction exists to stop being mistreated. double a_var = static_cast(sigma_corr) * sigma_corr; if (scaling_settings.GetExpectedVarianceMerge()) { - const double bkg_var = std::max(0.0, a_var - static_cast(image_scale_corr) * I_corr); + const double bkg_var = std::max(0.0, a_var - static_cast(image_scale_corr) + * std::max(0.0, static_cast(I_corr))); const double base = bkg_var + static_cast(image_scale_corr) * std::max(0.0, I_for_b); if (base > 0.0) a_var = base; diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 9b698dbe..81291b08 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -1213,7 +1213,9 @@ void RotationScaleMerge::Combine() { const double corr = r2.corr; const double I_corr = pooled_I(r2) * corr; const double sigma_corr = static_cast(r2.sigma) * corr; - const double bkg_var = sigma_corr * sigma_corr - corr * I_corr; + // max(0, I): a down-fluctuated partial carries no Poisson signal to remove, and + // removing a negative one inflates the background part instead of leaving it alone. + const double bkg_var = sigma_corr * sigma_corr - corr * std::max(0.0, I_corr); double var = std::max(0.0, bkg_var) + corr * std::max(0.0, F); if (!(var > 0.0)) var = sigma_corr * sigma_corr; const double w = 1.0 / var;