PixelRefine: Some improvements

image_analysis/pixel_refinement/PixelRefine.cpp

@@ -33,6 +33,7 @@ struct ReflGroup {
     double Itrue;             // reference intensity (held fixed)
     double R_bw_sq;           // bandwidth radial-width^2 contribution (0 = monochromatic)
     double pol;               // per-reflection polarization correction (raw = true * pol)
+    double Ibkg;              // local flat background (raw counts, constant over the shoebox)
     double predicted_x, predicted_y;
     std::vector<PixelObs> pixels;
 };
@@ -83,6 +84,23 @@ std::vector<uint8_t> BuildSpotMask(const std::vector<Reflection> &predicted, int
     return mask;
 }
 
+// Square shoebox bounds (inclusive) around a predicted spot, clamped to the
+// detector. The centre is rounded to the nearest pixel with std::lround so the
+// signal box is centred identically to the spot-core mask (BuildSpotMask) and
+// the local-background ring (EstimateLocalBackground), which also lround. Used by
+// Run and the diagnostic renderers so all three share one shoebox definition.
+struct ShoeboxBox { int min_x, max_x, min_y, max_y; };
+ShoeboxBox ShoeboxBounds(double px, double py, int radius, size_t xpixel, size_t ypixel) {
+    const int cx = static_cast<int>(std::lround(px));
+    const int cy = static_cast<int>(std::lround(py));
+    return {
+        std::max(cx - radius, 0),
+        std::min<int>(cx + radius, static_cast<int>(xpixel) - 1),
+        std::max(cy - radius, 0),
+        std::min<int>(cy + radius, static_cast<int>(ypixel) - 1)
+    };
+}
+
 // Local flat background around one shoebox, in raw detector counts. Samples the
 // square ring shoebox_radius < max(|dx|,|dy|) <= bkg_outer_radius centred on the
 // spot, dropping pixels that belong to any spot core (spot_mask) or carry a
@@ -633,16 +651,14 @@ void PixelRefine::Run(const T *image,
             g.Itrue = reference_data[hkl];
             g.R_bw_sq = bandwidth_radial_sq(refl.d);
             g.pol = polarization(refl.predicted_x, refl.predicted_y);
+            g.Ibkg = Ibkg;
             g.predicted_x = refl.predicted_x;
             g.predicted_y = refl.predicted_y;
 
-            const int min_y = std::max<int>(refl.predicted_y - radius, 0);
-            const int max_y = std::min<int>(refl.predicted_y + radius, ypixel - 1);
-            const int min_x = std::max<int>(refl.predicted_x - radius, 0);
-            const int max_x = std::min<int>(refl.predicted_x + radius, xpixel - 1);
+            const auto box = ShoeboxBounds(refl.predicted_x, refl.predicted_y, radius, xpixel, ypixel);
 
-            for (int y = min_y; y <= max_y; ++y) {
-                for (int x = min_x; x <= max_x; ++x) {
+            for (int y = box.min_y; y <= box.max_y; ++y) {
+                for (int x = box.min_x; x <= box.max_x; ++x) {
                     const size_t npixel = xpixel * y + x;
 
                     // Skip sentinel (masked / saturated) pixels. We assume the pixel
@@ -826,52 +842,100 @@ void PixelRefine::Run(const T *image,
     } // predict<->refine iterations
 
     // ---- Extract integrated reflections ---------------------------------------
-    // Profile fitting gives the recorded amplitude (against the tangential profile
-    // P_t):
-    //     J      = sum_p[ P_t,p (Iobs_p - Ibkg_p)/v_p ] / sum_p[ P_t,p^2 / v_p ]
+    // Profile fitting gives the recorded amplitude (fitting the tangential profile
+    // P_t against the background-subtracted pixels):
+    //     J      = sum_p[ P_t,p (Iobs_p - Ibkg)/v_p ] / sum_p[ P_t,p^2 / v_p ]
     //            ~ G * Itrue * B_term * partiality * pol     (recorded raw counts)
     //     var(J) = 1 / sum_p[ P_t,p^2 / v_p ]
     //
+    // Two SEPARATE fractions reduce the full intensity to what these pixels record:
+    //
+    //   partiality   - the radial / rocking dimension that a still does NOT sample.
+    //                  Only the slice of the reflection that crosses the Ewald
+    //                  sphere on this shot is recorded; <= 1. We DIVIDE it out to
+    //                  recover the full intensity. = profile-weighted P_radial.
+    //
+    //   completeness - the fraction of the spot's detector footprint that landed on
+    //                  live pixels (= profile captured by live pixels / profile over
+    //                  the whole shoebox). 1.0 when the spot sits fully on the
+    //                  detector; < 1.0 only when a detector edge, gap or mask clips
+    //                  it. Profile fitting already extrapolates over the missing
+    //                  pixels, so this is NOT applied to r.I - it is a quality flag.
+    //
     // Output split (Merge multiplies r.I * image_scale_corr and weights by
     // 1/(sigma*image_scale_corr)^2 - see Merge.cpp):
     //   r.I              = J / (B_term * partiality * pol) = G * Itrue
     //   r.sigma          = sqrt(var(J)) / (B_term * partiality * pol)
-    //   r.partiality     = profile-weighted peak radial factor in (0,1]  (Merge filter only)
+    //   r.partiality     = profile-weighted P_radial in (0,1]  (the rocking fraction)
+    //   r.completeness   = live/total tangential profile in (0,1]  (detector clipping)
     //   r.image_scale_corr = 1/G                         (per-image scale ONLY)
     // so r.I * image_scale_corr = Itrue. B, partiality and polarization live on the
     // intensity, G lives on image_scale_corr - one clean meaning per field.
+    //
+    // We walk the full (unclamped) shoebox once: every grid point feeds the total
+    // tangential profile (completeness denominator); points that are real, live
+    // detector pixels also feed the profile fit and the captured profile.
     data.reflections.reserve(groups.size());
     for (const auto &g : groups) {
-        double num = 0.0, den = 0.0, bkg_sum = 0.0;
-        double radial_sum = 0.0, radial_w = 0.0;
+        const int cx = static_cast<int>(std::lround(g.predicted_x));
+        const int cy = static_cast<int>(std::lround(g.predicted_y));
+
+        // Debye-Waller factor for this reflection (constant over its shoebox).
+        const double B_term = std::exp(-data.B_factor / (4.0 * g.d * g.d));
+
+        double num = 0.0, den = 0.0, bkg_sum = 0.0, radial_sum = 0.0;
+        double prof_live = 0.0, prof_full = 0.0; // tangential profile: captured / total
         size_t n = 0;
 
-        for (const auto &obs : g.pixels) {
-            PixelResidual pr(obs, 1.0, lambda, pixel_size, g.h, g.k, g.l, g.R_bw_sq, g.pol, data.crystal_system);
-            double q_sq, eps_r, eps_t_sq;
-            if (!pr.GeometryTerms(beam, &dist_mm, detector_rot, latt_vec0, latt_vec1, latt_vec2, q_sq,
-                eps_r, eps_t_sq))
-                continue;
-            if (!(data.R[0] > 0.0) || !(data.R[1] > 0.0))
-                continue;
+        for (int y = cy - radius; y <= cy + radius; ++y) {
+            for (int x = cx - radius; x <= cx + radius; ++x) {
+                // Geometry/profile for this grid point (valid even off the detector).
+                PixelObs probe{static_cast<double>(x), static_cast<double>(y), 0.0, g.Ibkg, 1.0};
+                PixelResidual pr(probe, 1.0, lambda, pixel_size, g.h, g.k, g.l,
+                                 g.R_bw_sq, g.pol, data.crystal_system);
+                double q_sq, eps_r, eps_t_sq;
+                if (!pr.GeometryTerms(beam, &dist_mm, detector_rot,
+                                      latt_vec0, latt_vec1, latt_vec2, q_sq, eps_r, eps_t_sq))
+                    continue;
+                if (!(data.R[0] > 0.0) || !(data.R[1] > 0.0))
+                    continue;
 
-            // Tangential profile shape -> fit weight (every pixel counts equally).
-            const double P_t = std::exp(-eps_t_sq / (data.R[1] * data.R[1]))
-                               / (M_PI * data.R[1] * data.R[1]);
-            // Peak-normalized radial factor (the partiality), in (0,1].
-            // Bandwidth-broadened radial width, matching the model in Model().
-            const double R0_eff_sq = data.R[0] * data.R[0] + g.R_bw_sq;
-            const double P_radial = std::exp(-eps_r * eps_r / R0_eff_sq);
+                // Tangential profile shape (area-normalized) -> the fit template.
+                const double P_t = std::exp(-eps_t_sq / (data.R[1] * data.R[1]))
+                                   / (M_PI * data.R[1] * data.R[1]);
+                prof_full += P_t; // whole shoebox, on- or off-detector
 
-            const double v = SafeInv(obs.weight * obs.weight, 1.0); // pixel variance
-            const double signal = obs.Iobs - obs.Ibkg;
+                // Only real, unmasked detector pixels carry signal.
+                if (x < 0 || x >= static_cast<int>(xpixel) || y < 0 || y >= static_cast<int>(ypixel))
+                    continue;
+                const size_t np = static_cast<size_t>(xpixel) * y + x;
+                if (image[np] == std::numeric_limits<T>::max())
+                    continue;
+                if (std::is_signed_v<T> && image[np] == std::numeric_limits<T>::min())
+                    continue;
 
-            num += P_t * signal / v;
-            den += P_t * P_t / v;
-            radial_sum += P_radial * P_t;   // weight partiality by the spot core
-            radial_w += P_t;
-            bkg_sum += obs.Ibkg;
-            ++n;
+                const double Iobs = static_cast<double>(image[np]); // raw counts
+                double v = std::max(Iobs, 0.0);                     // Poisson variance
+                if (!(v > 1.0))
+                    v = 1.0;
+
+                // Peak-normalized radial factor (the partiality), in (0,1]. The
+                // bandwidth-broadened radial width matches the model in Model().
+                const double R0_eff_sq = data.R[0] * data.R[0] + g.R_bw_sq;
+                const double P_radial = std::exp(-eps_r * eps_r / R0_eff_sq);
+
+                // Profile-fit accumulators. The amplitude estimator weights pixels by
+                // P_t^2/v, so the partiality (which de-scales that amplitude) MUST use
+                // the SAME weights - otherwise an R0_eff-dependent (resolution-
+                // dependent) factor is left behind in r.I.
+                const double w = P_t * P_t / v;
+                num += P_t * (Iobs - g.Ibkg) / v;
+                den += w;
+                radial_sum += P_radial * w;   // partiality weighted exactly like num/den
+                prof_live += P_t;             // captured tangential profile
+                bkg_sum += g.Ibkg;
+                ++n;
+            }
         }
 
         Reflection r{};
@@ -884,12 +948,12 @@ void PixelRefine::Run(const T *image,
         r.observed_x = NAN;
         r.observed_y = NAN;
         r.rlp = 1.0f;
-        r.partiality = (radial_w > 0.0) ? static_cast<float>(radial_sum / radial_w) : 1.0f;
+        r.partiality = (den > 0.0) ? static_cast<float>(radial_sum / den) : 1.0f;
+        r.completeness = (prof_full > 0.0) ? static_cast<float>(prof_live / prof_full) : 1.0f;
 
         if (den > 0.0 && n > 0) {
             const double I_amp = num / den;             // ~ G*Itrue*B_term*partiality*pol
             const double sigma_amp = std::sqrt(1.0 / den);
-            const double B_term = std::exp(-data.B_factor / (4.0 * g.d * g.d));
             const double corr = static_cast<double>(r.partiality) * B_term * g.pol; // B, partiality & pol
             r.bkg = static_cast<float>(bkg_sum / static_cast<double>(n));
             r.observed = true;
@@ -1010,13 +1074,10 @@ std::vector<float> PixelRefine::PredictImage(const T *image,
                                     refl.predicted_x, refl.predicted_y,
                                     radius, bkg_outer_radius, Ibkg);
 
-        const int min_y = std::max<int>(refl.predicted_y - radius, 0);
-        const int max_y = std::min<int>(refl.predicted_y + radius, ypixel - 1);
-        const int min_x = std::max<int>(refl.predicted_x - radius, 0);
-        const int max_x = std::min<int>(refl.predicted_x + radius, xpixel - 1);
+        const auto box = ShoeboxBounds(refl.predicted_x, refl.predicted_y, radius, xpixel, ypixel);
 
-        for (int y = min_y; y <= max_y; ++y) {
-            for (int x = min_x; x <= max_x; ++x) {
+        for (int y = box.min_y; y <= box.max_y; ++y) {
+            for (int x = box.min_x; x <= box.max_x; ++x) {
                 const size_t npixel = xpixel * y + x;
 
                 PixelObs obs{
@@ -1106,13 +1167,10 @@ std::vector<float> PixelRefine::ChiSquaredImage(const T *image,
                                      radius, bkg_outer_radius, Ibkg))
             continue;
 
-        const int min_y = std::max<int>(refl.predicted_y - radius, 0);
-        const int max_y = std::min<int>(refl.predicted_y + radius, ypixel - 1);
-        const int min_x = std::max<int>(refl.predicted_x - radius, 0);
-        const int max_x = std::min<int>(refl.predicted_x + radius, xpixel - 1);
+        const auto box = ShoeboxBounds(refl.predicted_x, refl.predicted_y, radius, xpixel, ypixel);
 
-        for (int y = min_y; y <= max_y; ++y) {
-            for (int x = min_x; x <= max_x; ++x) {
+        for (int y = box.min_y; y <= box.max_y; ++y) {
+            for (int x = box.min_x; x <= box.max_x; ++x) {
                 const size_t npixel = xpixel * y + x;
 
                 // Same gating as Run(): only pixels that actually enter the fit.