diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.h b/image_analysis/bragg_integration/BraggIntegrationEngine.h index e336ccbb..8fe5d1dc 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.h @@ -48,6 +48,18 @@ constexpr int N_SHELL = 6; // resolution shells for per-shell constexpr double STRONG_I_OVER_SIGMA = 5.0; // strong-spot threshold that seeds the profile constexpr int MIN_STRONG_PER_SHELL = 30; // below this a shell falls back to the global profile constexpr double C_CAPTURE = 2.5; // weak-spot radial capture term (monochromatic only) +// Per-pixel variance floor for the Kabsch fit weights (v = floor + signal). The detector noise floor is +// the quantization noise from rounding the charge-spread deposited energy to an integer: a uniform +// rounding error has variance 1/12. Electronic noise is far below this for both EIGER and JUNGFRAU. A +// larger floor (the previous 1.0) silently over-regularizes — it inflates weak-reflection sigma and +// pins the scaling error model's `a` term at its floor. +constexpr double PIXEL_VARIANCE_FLOOR = 1.0 / 12.0; + +// Guard against profile-fit runaways: on a weak / near-zero reflection the reweighted Kabsch iteration +// has no real peak to lock onto and can manufacture intensity the box sum never sees. Fall back to the +// summation (box-sum) intensity when the profile result disagrees with the summation seed by more than +// this many box-sum sigmas (a real fit agrees within counting noise, so the margin is generous). +constexpr double PROFILE_SUMMATION_MAX_NSIGMA = 10.0; } // namespace bragg_engine // One reflection's extracted intensity, produced by the derived engine and turned into a diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index 29b17211..b953892e 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -293,7 +293,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, } const int Gf = 2 * Rf + 1; - const double B = std::max(rh.bkg, 1.0); + const double B = std::max(rh.bkg, PIXEL_VARIANCE_FLOOR); double I = rh.I, den = 0.0; for (int iter = 0; iter < 4; ++iter) { double num = 0.0; @@ -314,7 +314,16 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, } if (!(den > 0.0)) continue; - results[i] = {static_cast(I), static_cast(std::sqrt(1.0 / den)), + // Guard against profile-fit runaways: on a weak / near-zero reflection the reweighted Kabsch + // iteration has no real peak to lock onto and can manufacture intensity the box sum never sees. + // Keep the profile intensity only if it agrees with the summation seed within the margin; + // otherwise fall back to the summation, which is robust there. + double sigma = std::sqrt(1.0 / den); + if (std::abs(I - rh.I) > PROFILE_SUMMATION_MAX_NSIGMA * rh.sigma) { + I = rh.I; + sigma = rh.sigma; + } + results[i] = {static_cast(I), static_cast(sigma), static_cast(rh.bkg), 0.0f, 0.0f, true, false}; } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index 82db56f3..6697c0e9 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -253,7 +253,7 @@ __global__ void build_profiles(const float *shell_grid, const float *global_grid // One block per reflection; the (possibly elongated) profile is built in shared memory. --- __global__ void fit(const int32_t *img, const float *px_x, const float *px_y, const int *cx_a, const int *cy_a, const float *dd, const unsigned long long *invd2mm, - const float *I_seed, const float *bkg_a, const uint8_t *ok_a, + const float *I_seed, const float *sigma_seed, const float *bkg_a, const uint8_t *ok_a, const float *shell_P, const float *global_P, const float *shell_sigma2, const float *global_sigma2, const int *shell_n, float *I_o, float *sigma_o, uint8_t *ok_o, BraggGpuParams p, int n) { @@ -308,7 +308,7 @@ __global__ void fit(const int32_t *img, const float *px_x, const float *px_y, } const int Rf = s_Rf, Gf = s_Gf, GfGf = Gf * Gf; - const float B = fmaxf(bkg, 1.0f); + const float B = fmaxf(bkg, (float) PIXEL_VARIANCE_FLOOR); if (threadIdx.x == 0) s_I = I_seed[i]; __syncthreads(); @@ -335,8 +335,16 @@ __global__ void fit(const int32_t *img, const float *px_x, const float *px_y, } if (threadIdx.x == 0) { - if (s_den > 0.0f) { I_o[i] = s_I; sigma_o[i] = sqrtf(1.0f / s_den); ok_o[i] = 1; } - else ok_o[i] = 0; + if (s_den > 0.0f) { + float I = s_I, sigma = sqrtf(1.0f / s_den); + // Guard against profile-fit runaways (see the CPU engine): fall back to the summation seed + // when the profile result diverges from it. + if (fabsf(I - I_seed[i]) > (float) PROFILE_SUMMATION_MAX_NSIGMA * sigma_seed[i]) { + I = I_seed[i]; + sigma = sigma_seed[i]; + } + I_o[i] = I; sigma_o[i] = sigma; ok_o[i] = 1; + } else ok_o[i] = 0; } } @@ -447,7 +455,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu d_shell_grid, d_global_grid, d_shell_n, d_global_n, d_shell_P, d_global_P, d_shell_sigma2, d_global_sigma2, p); fit<<>>(img, d_px_x, d_px_y, d_cx, d_cy, d_d, d_invd2, - d_I, d_bkg, d_ok, d_shell_P, d_global_P, + d_I, d_sigma, d_bkg, d_ok, d_shell_P, d_global_P, d_shell_sigma2, d_global_sigma2, d_shell_n, d_I, d_sigma, d_ok, p, n); }