PIxelRefine: Another iteration

image_analysis/pixel_refinement/PixelRefine.cpp

@@ -766,19 +766,20 @@ void PixelRefine::Run(const T *image,
     } // predict<->refine iterations
 
     // ---- Extract integrated reflections ---------------------------------------
-    // Two quantities are read back per reflection, using the *optimized* model:
+    // Profile fitting gives the recorded amplitude (against the normalized
+    // tangential profile P_t):
+    //     J      = sum_p[ P_t,p (Iobs_p - Ibkg_p)/v_p ] / sum_p[ P_t,p^2 / v_p ]
+    //            ~ G * Itrue * B_term * partiality            (recorded intensity)
+    //     var(J) = 1 / sum_p[ P_t,p^2 / v_p ]
     //
-    //  * Recorded (partial) intensity J by profile fitting against the normalized
-    //    tangential profile P_t (sum over the shoebox ~ 1):
-    //        J       = sum_p [ P_t,p (Iobs_p - Ibkg_p) / v_p ] / sum_p [ P_t,p^2 / v_p ]
-    //        var(J)  = 1 / sum_p [ P_t,p^2 / v_p ]
-    //    Because P_t is area-normalized, J estimates the integrated intensity
-    //    actually recorded on this image (not the full-spot intensity).
-    //
-    //  * Partiality: the profile-weighted mean of the *peak-normalized* radial
-    //    factor exp(-eps_r^2/R0^2) in (0,1], i.e. how close to the Ewald sphere
-    //    the reflection sits. Kept dimensionless for consistency with the rest of
-    //    the pipeline (image_scale_corr = rlp / partiality, full I = J / partiality).
+    // 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)   = G * Itrue  (B/partiality corrected)
+    //   r.sigma          = sqrt(var(J)) / (B_term * partiality)
+    //   r.partiality     = profile-weighted peak radial factor in (0,1]  (Merge filter only)
+    //   r.image_scale_corr = 1/G                         (per-image scale ONLY)
+    // so r.I * image_scale_corr = Itrue. B and partiality live on the intensity,
+    // G lives on image_scale_corr - one clean meaning per field.
     data.reflections.reserve(groups.size());
     for (const auto &g : groups) {
         double num = 0.0, den = 0.0, bkg_sum = 0.0;
@@ -826,16 +827,22 @@ void PixelRefine::Run(const T *image,
         r.partiality = (radial_w > 0.0) ? static_cast<float>(radial_sum / radial_w) : 1.0f;
 
         if (den > 0.0 && n > 0) {
-            r.I = static_cast<float>(num / den);
-            r.sigma = static_cast<float>(std::sqrt(1.0 / den));
+            const double I_amp = num / den;             // ~ G*Itrue*B_term*partiality
+            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; // B & partiality
             r.bkg = static_cast<float>(bkg_sum / static_cast<double>(n));
             r.observed = true;
-            // Put I onto a common (merge-ready) scale: I * image_scale_corr is the
-            // full-spot, scale/DW-corrected intensity. Fold in the refined G, the
-            // per-reflection B_term and the partiality (rlp = 1 here).
-            const double B_term = std::exp(-data.B_factor / (4.0 * g.d * g.d));
-            const double denom = static_cast<double>(r.partiality) * data.scale_factor * B_term;
-            r.image_scale_corr = (denom > 0.0) ? static_cast<float>(r.rlp / denom) : NAN;
+
+            if (corr > 0.0 && data.scale_factor > 0.0) {
+                r.I = static_cast<float>(I_amp / corr);       // = G*Itrue
+                r.sigma = static_cast<float>(sigma_amp / corr);
+                r.image_scale_corr = static_cast<float>(1.0 / data.scale_factor); // G only
+            } else {
+                r.I = static_cast<float>(I_amp);
+                r.sigma = static_cast<float>(sigma_amp);
+                r.image_scale_corr = NAN;
+            }
         } else {
             r.I = 0.0f;
             r.sigma = NAN;
@@ -844,6 +851,36 @@ void PixelRefine::Run(const T *image,
         }
         data.reflections.push_back(r);
     }
+
+    // ---- Per-image CC vs reference (the half/ref correlation diagnostic) -------
+    // Pearson CC of the scaled estimate (r.I * image_scale_corr = Itrue_est)
+    // against the reference intensities, over the matched reflections.
+    {
+        double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0;
+        size_t cn = 0;
+        for (const auto &r : data.reflections) {
+            if (!r.observed || !std::isfinite(r.I) || !std::isfinite(r.image_scale_corr))
+                continue;
+            const auto it = reference_data.find(hkl_key_generator(r));
+            if (it == reference_data.end())
+                continue;
+            const double x = static_cast<double>(r.I) * static_cast<double>(r.image_scale_corr);
+            const double y = it->second;
+            if (!std::isfinite(x) || !std::isfinite(y))
+                continue;
+            sx += x; sy += y; sxx += x * x; syy += y * y; sxy += x * y; ++cn;
+        }
+        data.cc = NAN;
+        data.cc_n = static_cast<int64_t>(cn);
+        if (cn >= 3) {
+            const double nd = static_cast<double>(cn);
+            const double cov = sxy - sx * sy / nd;
+            const double vx = sxx - sx * sx / nd;
+            const double vy = syy - sy * sy / nd;
+            if (vx > 0.0 && vy > 0.0)
+                data.cc = cov / std::sqrt(vx * vy);
+        }
+    }
 }
 
 std::vector<float> PixelRefine::PredictImage(const AzimuthalIntegrationProfile &profile,