PixelRefine: Simplify (remove Lorentz correction, remove background from azimuthal integration)

image_analysis/pixel_refinement/PixelRefine.cpp

@@ -3,6 +3,11 @@
 
 #include "PixelRefine.h"
 
+#include <algorithm>
+#include <cmath>
+#include <limits>
+#include <vector>
+
 #include <Eigen/Dense>
 #include <ceres/ceres.h>
 #include <ceres/rotation.h>
@@ -11,17 +16,14 @@
 
 namespace {
 
-// Per-pixel observation, in *corrected* intensity units (solid-angle and
-// polarization correction already folded in, consistently for signal and
-// background). Geometry-independent quantities are precomputed here so that the
-// Ceres cost functor stays cheap.
+// Per-pixel observation, in *raw* detector counts (no per-pixel solid-angle or
+// polarization correction - same units the "normal" integrator works in; the
+// per-reflection polarization correction is applied via ReflGroup::pol).
 struct PixelObs {
     double x, y;          // detector pixel coordinate
-    double Iobs;          // corrected pixel value (signal + background)
-    double Ibkg;          // corrected background estimate (azimuthal bin mean)
+    double Iobs;          // raw pixel value (signal + background)
+    double Ibkg;          // local background estimate (per-shoebox level, raw counts)
     double weight;        // 1 / sigma_pixel
-    double A_recip;       // reciprocal-space area subtended by the pixel (Jacobian)
-    double angle_rad;     // goniometer angle of this observation
 };
 
 // One reflection together with the pixels of its shoebox.
@@ -30,6 +32,7 @@ struct ReflGroup {
     double d;
     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 predicted_x, predicted_y;
     std::vector<PixelObs> pixels;
 };
@@ -40,13 +43,99 @@ double SafeInv(double x, double fallback) {
     return 1.0 / x;
 }
 
+// Median of a vector (in place, partially reorders it).
+double MedianInPlace(std::vector<double> &v) {
+    if (v.empty())
+        return 0.0;
+    const size_t mid = v.size() / 2;
+    std::nth_element(v.begin(), v.begin() + mid, v.end());
+    if (v.size() % 2 == 1)
+        return v[mid];
+    const double hi = v[mid];
+    std::nth_element(v.begin(), v.begin() + mid - 1, v.begin() + mid);
+    return 0.5 * (v[mid - 1] + hi);
+}
+
+// Mask marking the *core* (radius `radius`) of every predicted spot, so that the
+// local-background sampling of one reflection never picks up a neighbouring
+// reflection's signal. Same idea as BraggIntegrate2D::BuildReflectionMask.
+std::vector<uint8_t> BuildSpotMask(const std::vector<Reflection> &predicted, int nrefl,
+                                   size_t xpixel, size_t ypixel, int radius) {
+    std::vector<uint8_t> mask(xpixel * ypixel, 0);
+    const double r_sq = static_cast<double>(radius) * radius;
+    for (int i = 0; i < nrefl; ++i) {
+        const auto &r = predicted[i];
+        const int cx = static_cast<int>(std::lround(r.predicted_x));
+        const int cy = static_cast<int>(std::lround(r.predicted_y));
+        const int x0 = std::max(0, cx - radius);
+        const int x1 = std::min<int>(static_cast<int>(xpixel) - 1, cx + radius);
+        const int y0 = std::max(0, cy - radius);
+        const int y1 = std::min<int>(static_cast<int>(ypixel) - 1, cy + radius);
+        for (int y = y0; y <= y1; ++y) {
+            for (int x = x0; x <= x1; ++x) {
+                const double dx = x - r.predicted_x;
+                const double dy = y - r.predicted_y;
+                if (dx * dx + dy * dy <= r_sq)
+                    mask[static_cast<size_t>(xpixel) * y + x] = 1;
+            }
+        }
+    }
+    return mask;
+}
+
+// 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
+// masked/saturated sentinel, and returns the median (robust to residual spot
+// tails / zingers). Mirrors the local-background of BraggIntegrate2D, replacing
+// the azimuthal-bin mean that proved a poor proxy for reflection background.
+template<class T>
+bool EstimateLocalBackground(const T *image,
+                             const std::vector<uint8_t> &spot_mask,
+                             size_t xpixel, size_t ypixel,
+                             double cx, double cy,
+                             int shoebox_radius, int bkg_outer_radius,
+                             double &bkg_mean) {
+    const int icx = static_cast<int>(std::lround(cx));
+    const int icy = static_cast<int>(std::lround(cy));
+    const int x0 = std::max(0, icx - bkg_outer_radius);
+    const int x1 = std::min<int>(static_cast<int>(xpixel) - 1, icx + bkg_outer_radius);
+    const int y0 = std::max(0, icy - bkg_outer_radius);
+    const int y1 = std::min<int>(static_cast<int>(ypixel) - 1, icy + bkg_outer_radius);
+
+    std::vector<double> vals;
+    vals.reserve(static_cast<size_t>((x1 - x0 + 1) * (y1 - y0 + 1)));
+    for (int y = y0; y <= y1; ++y) {
+        for (int x = x0; x <= x1; ++x) {
+            // Skip the square shoebox core: that is signal, not background.
+            if (std::abs(x - icx) <= shoebox_radius && std::abs(y - icy) <= shoebox_radius)
+                continue;
+            const size_t np = static_cast<size_t>(xpixel) * y + x;
+            if (spot_mask[np])
+                continue;
+            const T raw = image[np];
+            if (raw == std::numeric_limits<T>::max())
+                continue;
+            if (std::is_signed_v<T> && raw == std::numeric_limits<T>::min())
+                continue;
+            vals.push_back(static_cast<double>(raw));
+        }
+    }
+
+    if (vals.size() < 5)
+        return false;
+
+    bkg_mean = MedianInPlace(vals);
+    return true;
+}
+
 // Per-pixel: map a detector pixel through the current geometry into the
 // reference reciprocal frame. Cheap (a few trig + one rotation); depends on the
 // pixel and the detector geometry, not on the lattice.
 template<typename T>
 void ObservedRecip(const T *beam, const T *distance_mm, const T *detector_rot,
-                   const T *rotation_axis, double obs_x, double obs_y,
-                   double pixel_size, double inv_lambda, double angle_rad,
+                   double obs_x, double obs_y,
+                   double pixel_size, double inv_lambda,
                    Eigen::Matrix<T, 3, 1> &e_obs_recip) {
     // PyFAI convention (left-handed for rot1/rot2): rot3 = 0 assumed.
     const T c1 = ceres::cos(detector_rot[0]);
@@ -73,14 +162,7 @@ void ObservedRecip(const T *beam, const T *distance_mm, const T *detector_rot,
         y * inv_norm * T(inv_lambda),
         (z * inv_norm - T(1.0)) * T(inv_lambda)
     };
-    const T aa_back[3] = {
-        T(angle_rad) * rotation_axis[0],
-        T(angle_rad) * rotation_axis[1],
-        T(angle_rad) * rotation_axis[2]
-    };
-    T recip_obs[3];
-    ceres::AngleAxisRotatePoint(aa_back, recip_raw, recip_obs);
-    e_obs_recip = Eigen::Matrix<T, 3, 1>(recip_obs[0], recip_obs[1], recip_obs[2]);
+    e_obs_recip = Eigen::Matrix<T, 3, 1>(recip_raw[0], recip_raw[1], recip_raw[2]);
 }
 
 // Per-reflection: predicted node g_hkl, |g_hkl|^2, and the Ewald-sphere normal.
@@ -178,20 +260,23 @@ bool PredictedNode(const T *p0, const T *p1, const T *p2,
 // ---------------------------------------------------------------------------
 // Cost functor
 //
-//   I_pred(pixel) = G * Itrue * B_term * P_radial * P_tangential + I_bkg
+//   I_pred(pixel) = G * Itrue * B_term * P_radial * P_tangential * pol + I_bkg
 //
 //   B_term       = exp(-B |q|^2 / 4)                    (Debye-Waller)
 //   P_radial     = exp(-eps_r^2 / R0_eff^2)             (partiality: fraction of
 //                                                        the mosaic blob on the
 //                                                        Ewald sphere; <= 1)
-//   P_tangential = A_recip/(pi R1^2) * exp(-eps_t^2/R1^2)(spatial profile on the
-//                                                        detector, normalized so
-//                                                        that sum over pixels ~ 1)
+//   P_tangential = exp(-eps_t^2/R1^2) / (pi R1^2)       (Gaussian spatial profile
+//                                                        in the Ewald tangent plane)
+//   pol          = per-reflection polarization correction (raw = true * pol),
+//                  evaluated once at the predicted spot position (as in
+//                  BraggIntegrate2D). 1 if polarization is disabled.
 //
-// The tangential factor is what makes this "profile fitting": summing
-// I_pred - I_bkg over the shoebox reproduces G * Itrue * B_term * P_radial.
-// The 1/(pi R1^2) normalization is the missing piece that decouples the profile
-// width R1 from the overall scale G.
+// Everything is in *raw* detector counts: there is no per-pixel solid-angle or
+// area (Lorentz/Jacobian) weighting - each pixel counts equally, like the normal
+// integrator. The tangential factor is what makes this "profile fitting"; the
+// 1/(pi R1^2) normalization keeps the profile width R1 from soaking up the
+// overall scale G.
 //
 // X-ray bandwidth: a spread in lambda is a spread in the Ewald-sphere radius,
 // i.e. a purely *radial* thickening of the shell. It adds (in quadrature) a
@@ -207,13 +292,13 @@ struct PixelResidual {
     PixelResidual(const PixelObs &obs, double Itrue,
                   double lambda, double pixel_size,
                   double exp_h, double exp_k, double exp_l,
-                  double R_bw_sq,
+                  double R_bw_sq, double pol,
                   gemmi::CrystalSystem symmetry)
         : Itrue(Itrue), Iobs(obs.Iobs), Ibkg(obs.Ibkg), weight(obs.weight),
-          A_recip(obs.A_recip), obs_x(obs.x), obs_y(obs.y),
+          obs_x(obs.x), obs_y(obs.y),
           inv_lambda(1.0 / lambda), pixel_size(pixel_size),
           exp_h(exp_h), exp_k(exp_k), exp_l(exp_l),
-          R_bw_sq(R_bw_sq), angle_rad(obs.angle_rad), symmetry(symmetry) {
+          R_bw_sq(R_bw_sq), pol(pol), symmetry(symmetry) {
         if (std::fabs(lambda) < 1e-6)
             throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                                   "Lambda cannot be close to zero");
@@ -228,14 +313,14 @@ struct PixelResidual {
     bool GeometryTerms(const T *const beam,
                        const T *const distance_mm,
                        const T *const detector_rot,
-                       const T *const rotation_axis,
                        const T *const p0,
                        const T *const p1,
                        const T *const p2,
                        T &q_sq, T &eps_radial, T &eps_tang_sq) const {
         Eigen::Matrix<T, 3, 1> e_obs_recip;
-        ObservedRecip(beam, distance_mm, detector_rot, rotation_axis,
-                      obs_x, obs_y, pixel_size, inv_lambda, angle_rad, e_obs_recip);
+        ObservedRecip(beam, distance_mm, detector_rot,
+                      obs_x, obs_y, pixel_size, inv_lambda,
+                      e_obs_recip);
 
         Eigen::Matrix<T, 3, 1> e_pred_recip, n_radial;
         if (!PredictedNode(p0, p1, p2, exp_h, exp_k, exp_l, symmetry, inv_lambda,
@@ -251,12 +336,12 @@ struct PixelResidual {
     // Assembles the full model intensity for the pixel from the geometry terms.
     template<typename T>
     bool Model(const T *const beam, const T *const distance_mm,
-               const T *const detector_rot, const T *const rotation_axis,
+               const T *const detector_rot,
                const T *const p0, const T *const p1, const T *const p2,
                const T *const scale_factor, const T *const B, const T *const R,
                T &Ipred) const {
         T q_sq, eps_radial, eps_tang_sq;
-        if (!GeometryTerms(beam, distance_mm, detector_rot, rotation_axis,
+        if (!GeometryTerms(beam, distance_mm, detector_rot,
                            p0, p1, p2, q_sq, eps_radial, eps_tang_sq))
             return false;
 
@@ -268,10 +353,9 @@ struct PixelResidual {
         // Separable Gaussian spot model:
         //   radial  P_r(e) = exp(-e^2/R0_eff^2)                      (peak-normalized, in (0,1])
         //   tangent g_t(e) = exp(-|e|^2/R1^2) / (pi R1^2)            [1/A^-2]
-        // The pixel captures the fraction g_t * A_recip of the tangential profile
-        // (A_recip = reciprocal area the pixel subtends; sum over shoebox ~ 1).
-        // The radial factor is the still-image partiality (how far the reflection
-        // sits from the Ewald sphere); the overall scale is carried by the free G.
+        // Every pixel counts equally (no area/Lorentz weighting); the radial factor
+        // is the still-image partiality (how far the reflection sits from the Ewald
+        // sphere); the overall scale is carried by the free G.
         //
         // IMPORTANT: the radial factor MUST use the same convention here as the
         // extraction's `partiality` (peak-normalized), otherwise image_scale_corr
@@ -280,10 +364,10 @@ struct PixelResidual {
         // R0_eff folds in the energy-bandwidth broadening via R_bw_sq.
         const T R0_eff_sq = R[0] * R[0] + T(R_bw_sq);
         const T P_radial = ceres::exp(-eps_radial * eps_radial / R0_eff_sq);
-        const T P_tang = T(A_recip) * ceres::exp(-eps_tang_sq / (R[1] * R[1]))
+        const T P_tang = ceres::exp(-eps_tang_sq / (R[1] * R[1]))
                          / (T(M_PI) * R[1] * R[1]);
 
-        const T signal = scale_factor[0] * T(Itrue) * B_term * P_radial * P_tang;
+        const T signal = scale_factor[0] * T(Itrue) * B_term * P_radial * P_tang * T(pol);
         Ipred = signal + T(Ibkg);
         return true;
     }
@@ -292,7 +376,6 @@ struct PixelResidual {
     bool operator()(const T *const beam,
                     const T *const distance_mm,
                     const T *const detector_rot,
-                    const T *const rotation_axis,
                     const T *const p0,
                     const T *const p1,
                     const T *const p2,
@@ -301,21 +384,20 @@ struct PixelResidual {
                     const T *const R,
                     T *residual) const {
         T Ipred;
-        if (!Model(beam, distance_mm, detector_rot, rotation_axis,
-                   p0, p1, p2, scale_factor, B, R, Ipred))
+        if (!Model(beam, distance_mm, detector_rot, p0, p1, p2, scale_factor, B, R, Ipred))
             return false;
 
         residual[0] = (Ipred - T(Iobs)) * T(weight);
         return true;
     }
 
-    const double Itrue, Iobs, Ibkg, weight, A_recip;
+    const double Itrue, Iobs, Ibkg, weight;
     const double obs_x, obs_y;
     const double inv_lambda;
     const double pixel_size;
     const double exp_h, exp_k, exp_l;
     const double R_bw_sq;   // bandwidth radial-width^2 contribution (0 = monochromatic)
-    const double angle_rad;
+    const double pol;       // per-reflection polarization correction
     gemmi::CrystalSystem symmetry;
 };
 
@@ -333,27 +415,25 @@ struct PixelResidual {
 struct ShoeboxResidual {
     ShoeboxResidual(const ReflGroup &g, double lambda, double pixel_size,
                     gemmi::CrystalSystem symmetry)
-        : pixels(g.pixels), Itrue(g.Itrue), R_bw_sq(g.R_bw_sq),
+        : pixels(g.pixels), Itrue(g.Itrue), R_bw_sq(g.R_bw_sq), pol(g.pol),
           exp_h(g.h), exp_k(g.k), exp_l(g.l),
           inv_lambda(1.0 / lambda), pixel_size(pixel_size),
-          angle_rad(g.pixels.empty() ? 0.0 : g.pixels.front().angle_rad),
           symmetry(symmetry) {}
 
     template<typename T>
     bool operator()(const T *const *params, T *residual) const {
         // Parameter blocks (order matches AddParameterBlock in Run):
-        // 0 beam[2] 1 distance[1] 2 detector_rot[2] 3 rotation_axis[3]
-        // 4 p0[3] 5 p1[3] 6 p2[3] 7 scale[1] 8 B[1] 9 R[2]
+        // 0 beam[2] 1 distance[1] 2 detector_rot[2]
+        // 3 p0[3] 4 p1[3] 5 p2[3] 6 scale[1] 7 B[1] 8 R[2]
         const T *beam = params[0];
         const T *distance_mm = params[1];
         const T *detector_rot = params[2];
-        const T *rotation_axis = params[3];
-        const T *p0 = params[4];
-        const T *p1 = params[5];
-        const T *p2 = params[6];
-        const T *scale_factor = params[7];
-        const T *B = params[8];
-        const T *R = params[9];
+        const T *p0 = params[3];
+        const T *p1 = params[4];
+        const T *p2 = params[5];
+        const T *scale_factor = params[6];
+        const T *B = params[7];
+        const T *R = params[8];
 
         if (R[0] < T(1e-10) || R[1] < T(1e-10))
             return false;
@@ -373,18 +453,18 @@ struct ShoeboxResidual {
             const PixelObs &obs = pixels[i];
 
             Eigen::Matrix<T, 3, 1> e_obs_recip;
-            ObservedRecip(beam, distance_mm, detector_rot, rotation_axis,
-                          obs.x, obs.y, pixel_size, inv_lambda, angle_rad, e_obs_recip);
+            ObservedRecip(beam, distance_mm, detector_rot,
+                          obs.x, obs.y, pixel_size, inv_lambda,  e_obs_recip);
 
             const Eigen::Matrix<T, 3, 1> delta_q = e_obs_recip - e_pred_recip;
             const T eps_radial = delta_q.dot(n_radial);
             const T eps_tang_sq = (delta_q - eps_radial * n_radial).squaredNorm();
 
             const T P_radial = ceres::exp(-eps_radial * eps_radial / R0_eff_sq);
-            const T P_tang = T(obs.A_recip) * ceres::exp(-eps_tang_sq / (R[1] * R[1]))
+            const T P_tang = ceres::exp(-eps_tang_sq / (R[1] * R[1]))
                              / (T(M_PI) * R[1] * R[1]);
 
-            const T signal = scale_factor[0] * T(Itrue) * B_term * P_radial * P_tang;
+            const T signal = scale_factor[0] * T(Itrue) * B_term * P_radial * P_tang * T(pol);
             const T Ipred = signal + T(obs.Ibkg);
             residual[i] = (Ipred - T(obs.Iobs)) * T(obs.weight);
         }
@@ -392,17 +472,15 @@ struct ShoeboxResidual {
     }
 
     std::vector<PixelObs> pixels;
-    const double Itrue, R_bw_sq;
+    const double Itrue, R_bw_sq, pol;
     const double exp_h, exp_k, exp_l;
-    const double inv_lambda, pixel_size, angle_rad;
+    const double inv_lambda, pixel_size;
     gemmi::CrystalSystem symmetry;
 };
 
 PixelRefine::PixelRefine(const DiffractionExperiment &experiment,
-                         const AzimuthalIntegrationMapping &mapping,
                          const std::vector<MergedReflection> &reference)
-    : mapping(mapping),
-      xpixel(experiment.GetXPixelsNum()),
+    : xpixel(experiment.GetXPixelsNum()),
       ypixel(experiment.GetYPixelsNum()),
       experiment(experiment),
       hkl_key_generator(experiment.GetScalingSettings().GetMergeFriedel(),
@@ -413,19 +491,13 @@ PixelRefine::PixelRefine(const DiffractionExperiment &experiment,
 
 void PixelRefine::BuildParameterBlocks(const PixelRefineData &data,
                                        double beam[2], double &dist_mm,
-                                       double detector_rot[2], double rot_vec[3],
+                                       double detector_rot[2],
                                        double latt_vec0[3], double latt_vec1[3], double latt_vec2[3]) const {
     beam[0] = data.geom.GetBeamX_pxl();
     beam[1] = data.geom.GetBeamY_pxl();
     dist_mm = data.geom.GetDetectorDistance_mm();
     detector_rot[0] = data.geom.GetPoniRot1_rad();
     detector_rot[1] = data.geom.GetPoniRot2_rad();
-    rot_vec[0] = 1.0; rot_vec[1] = 0.0; rot_vec[2] = 0.0;
-    if (auto axis = data.geom.GetRotation()) {
-        rot_vec[0] = axis->GetAxis().x;
-        rot_vec[1] = axis->GetAxis().y;
-        rot_vec[2] = axis->GetAxis().z;
-    }
 
     for (int i = 0; i < 3; ++i)
         latt_vec0[i] = latt_vec1[i] = latt_vec2[i] = 0.0;
@@ -463,7 +535,6 @@ void PixelRefine::BuildParameterBlocks(const PixelRefineData &data,
 
 template<class T>
 void PixelRefine::Run(const T *image,
-                      const AzimuthalIntegrationProfile &profile,
                       BraggPrediction &prediction,
                       PixelRefineData &data) {
     data.solved = false;
@@ -479,29 +550,17 @@ void PixelRefine::Run(const T *image,
         .bandwidth_sigma = static_cast<float>(data.bandwidth)  // relative Δλ/λ sigma
     };
 
-    const auto azim_result = profile.GetResult();
-    const auto azim_std = profile.GetStd();
-    const auto &pixel_to_bin = mapping.GetPixelToBin();
-    const auto &corrections = mapping.Corrections();
-    // pixel_to_bin stores the *full* bin index (azimuthal_sector * q_bins + q_bin),
-    // so the valid range is the total number of bins, i.e. the profile size - NOT
-    // GetAzimuthalBinCount() (which is only the number of azimuthal sectors).
-    const int total_bin_count = static_cast<int>(azim_result.size());
-
-    const double angle_rad = data.angle_deg * M_PI / 180.0;
     const int radius = data.shoebox_radius;
+    const int bkg_outer_radius = std::max(radius + 1, data.bkg_outer_radius);
 
-    // Exact reciprocal-space area a 1x1 pixel subtends, |dq/dx x dq/dy|, via
-    // finite differences of the detector->reciprocal map. This is the Jacobian
-    // between the curved Ewald-sphere sampling and flat reciprocal space, and it
-    // is exactly the geometric factor that plays the role of the Lorentz factor
-    // for stills: where the sphere grazes reciprocal space obliquely, a pixel
-    // covers more reciprocal volume and the captured fraction grows. It tracks
-    // the refined geometry because it reads the current data.geom each iteration.
-    auto recip_area = [&](double x, double y) -> double {
-        const Coord qx = data.geom.DetectorToRecip(x + 0.5, y) - data.geom.DetectorToRecip(x - 0.5, y);
-        const Coord qy = data.geom.DetectorToRecip(x, y + 0.5) - data.geom.DetectorToRecip(x, y - 0.5);
-        return (qx % qy).Length();
+    // Per-reflection polarization correction (raw = true * pol), evaluated once at
+    // the predicted spot - same handling as BraggIntegrate2D. Identity if disabled.
+    const auto pol_factor = experiment.GetPolarizationFactor();
+    auto polarization = [&](double x, double y) -> double {
+        if (!pol_factor)
+            return 1.0;
+        return data.geom.CalcAzIntPolarizationCorr(static_cast<float>(x), static_cast<float>(y),
+                                                   pol_factor.value());
     };
 
     // Bandwidth radial-width^2 (in the code's R = sqrt(2)*sigma convention):
@@ -524,7 +583,6 @@ void PixelRefine::Run(const T *image,
     double beam[2] = {0, 0};
     double dist_mm = data.geom.GetDetectorDistance_mm();
     double detector_rot[2] = {0, 0};
-    double rot_vec[3] = {1.0, 0.0, 0.0};
     double latt_vec0[3] = {0, 0, 0}; // orientation (Rodrigues)
     double latt_vec1[3] = {0, 0, 0}; // lengths
     double latt_vec2[3] = {0, 0, 0}; // angles (rad)
@@ -546,12 +604,27 @@ void PixelRefine::Run(const T *image,
         // nrefl entries are valid for this image (the rest are stale/zeroed).
         groups.clear();
         const auto &predicted = prediction.GetReflections();
+
+        // Spot-core mask over ALL predicted reflections, so each reflection's
+        // local background ignores pixels that belong to a neighbouring spot.
+        const auto spot_mask = BuildSpotMask(predicted, nrefl, xpixel, ypixel, radius);
+
         for (int ri = 0; ri < nrefl; ++ri) {
             const auto &refl = predicted[ri];
             const auto hkl = hkl_key_generator(refl);
             if (!reference_data.contains(hkl))
                 continue;
 
+            // Local flat background from the ring around the shoebox (raw counts).
+            // No azimuthal fallback: if we cannot estimate a clean local background
+            // the reflection is dropped, exactly as BraggIntegrate2D marks it
+            // unobserved when fewer than a handful of background pixels survive.
+            double Ibkg = 0.0;
+            if (!EstimateLocalBackground(image, spot_mask, xpixel, ypixel,
+                                         refl.predicted_x, refl.predicted_y,
+                                         radius, bkg_outer_radius, Ibkg))
+                continue;
+
             ReflGroup g;
             g.h = refl.h;
             g.k = refl.k;
@@ -559,6 +632,7 @@ void PixelRefine::Run(const T *image,
             g.d = refl.d;
             g.Itrue = reference_data[hkl];
             g.R_bw_sq = bandwidth_radial_sq(refl.d);
+            g.pol = polarization(refl.predicted_x, refl.predicted_y);
             g.predicted_x = refl.predicted_x;
             g.predicted_y = refl.predicted_y;
 
@@ -570,27 +644,18 @@ void PixelRefine::Run(const T *image,
             for (int y = min_y; y <= max_y; ++y) {
                 for (int x = min_x; x <= max_x; ++x) {
                     const size_t npixel = xpixel * y + x;
-                    const int azim_bin = pixel_to_bin[npixel];
 
-                    // Skip pixels not mapped to a bin or carrying a sentinel
-                    // (masked / saturated) value. We assume the pixel mask is
-                    // already applied upstream.
-                    if (azim_bin >= total_bin_count)
-                        continue;
+                    // Skip sentinel (masked / saturated) pixels. We assume the pixel
+                    // mask is already applied upstream (encoded as the sentinel).
                     if (image[npixel] == std::numeric_limits<T>::max())
                         continue;
                     if (std::is_signed_v<T> && (image[npixel] == std::numeric_limits<T>::min()))
                         continue;
 
-                    const double correction = corrections[npixel];
-                    const double Ibkg = azim_result[azim_bin]; // already in corrected units
-                    const double Ibkg_sigma = azim_std[azim_bin];
-                    const double raw = static_cast<double>(image[npixel]);
-                    const double Iobs = raw * correction;
+                    const double Iobs = static_cast<double>(image[npixel]); // raw counts
 
-                    // Per-pixel variance: Poisson noise of the corrected counts
-                    // (var(c*N) = c^2 * N = c * Iobs) plus the background spread.
-                    double var = correction * std::max(Iobs, 0.0) + Ibkg_sigma * Ibkg_sigma;
+                    // Per-pixel variance: Poisson noise of the raw counts.
+                    double var = std::max(Iobs, 0.0);
                     if (!(var > 1.0))
                         var = 1.0;
 
@@ -599,9 +664,7 @@ void PixelRefine::Run(const T *image,
                         .y = static_cast<double>(y),
                         .Iobs = Iobs,
                         .Ibkg = Ibkg,
-                        .weight = 1.0 / std::sqrt(var),
-                        .A_recip = recip_area(x, y),
-                        .angle_rad = angle_rad
+                        .weight = 1.0 / std::sqrt(var)
                     };
                     g.pixels.push_back(obs);
                 }
@@ -615,7 +678,7 @@ void PixelRefine::Run(const T *image,
             return;
 
         // ---- 3. Set up parameter blocks (geometry part mirrors XtalOptimizer) -
-        BuildParameterBlocks(data, beam, dist_mm, detector_rot, rot_vec,
+        BuildParameterBlocks(data, beam, dist_mm, detector_rot,
                              latt_vec0, latt_vec1, latt_vec2);
 
         // ---- 4. Build the problem ---------------------------------------------
@@ -630,7 +693,6 @@ void PixelRefine::Run(const T *image,
             cost->AddParameterBlock(2); // beam
             cost->AddParameterBlock(1); // distance
             cost->AddParameterBlock(2); // detector_rot
-            cost->AddParameterBlock(3); // rotation_axis
             cost->AddParameterBlock(3); // p0 (orientation)
             cost->AddParameterBlock(3); // p1 (lengths)
             cost->AddParameterBlock(3); // p2 (angles)
@@ -643,7 +705,7 @@ void PixelRefine::Run(const T *image,
             // per-pixel Huber. Per-pixel sigma weighting is retained; per-pixel
             // outlier rejection (zingers) is a TODO if needed.
             problem.AddResidualBlock(cost, nullptr,
-                                     beam, &dist_mm, detector_rot, rot_vec,
+                                     beam, &dist_mm, detector_rot,
                                      latt_vec0, latt_vec1, latt_vec2,
                                      &data.scale_factor, &data.B_factor, data.R);
             residual_pixels += g.pixels.size();
@@ -687,9 +749,6 @@ void PixelRefine::Run(const T *image,
             }
         }
 
-        if (!data.refine_rotation_axis)
-            problem.SetParameterBlockConstant(rot_vec);
-
         if (data.refine_scale)
             problem.SetParameterLowerBound(&data.scale_factor, 0, 0.0);
         else
@@ -767,20 +826,20 @@ void PixelRefine::Run(const T *image,
     } // predict<->refine iterations
 
     // ---- Extract integrated reflections ---------------------------------------
-    // Profile fitting gives the recorded amplitude (against the normalized
-    // tangential profile P_t):
+    // 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 ]
-    //            ~ G * Itrue * B_term * partiality            (recorded intensity)
+    //            ~ G * Itrue * B_term * partiality * pol     (recorded raw counts)
     //     var(J) = 1 / sum_p[ P_t,p^2 / v_p ]
     //
     // 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.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.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.
+    // 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.
     data.reflections.reserve(groups.size());
     for (const auto &g : groups) {
         double num = 0.0, den = 0.0, bkg_sum = 0.0;
@@ -788,16 +847,16 @@ void PixelRefine::Run(const T *image,
         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, data.crystal_system);
+            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, rot_vec,
-                                  latt_vec0, latt_vec1, latt_vec2, 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;
 
-            // Normalized tangential profile (sum over shoebox ~ 1) -> fit weight.
-            const double P_t = obs.A_recip * std::exp(-eps_t_sq / (data.R[1] * data.R[1]))
+            // 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().
@@ -828,10 +887,10 @@ 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) {
-            const double I_amp = num / den;             // ~ G*Itrue*B_term*partiality
+            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; // B & partiality
+            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;
 
@@ -884,7 +943,8 @@ void PixelRefine::Run(const T *image,
     }
 }
 
-std::vector<float> PixelRefine::PredictImage(const AzimuthalIntegrationProfile &profile,
+template<class T>
+std::vector<float> PixelRefine::PredictImage(const T *image,
                                              BraggPrediction &prediction,
                                              const PixelRefineData &data,
                                              bool include_background) const {
@@ -892,18 +952,16 @@ std::vector<float> PixelRefine::PredictImage(const AzimuthalIntegrationProfile &
 
     const double lambda = data.geom.GetWavelength_A();
     const double pixel_size = data.geom.GetPixelSize_mm();
-    const auto azim_result = profile.GetResult();
-    const auto &pixel_to_bin = mapping.GetPixelToBin();
-    const auto &corrections = mapping.Corrections();
-    const int total_bin_count = static_cast<int>(azim_result.size());
-    const double angle_rad = data.angle_deg * M_PI / 180.0;
     const int radius = data.shoebox_radius;
+    const int bkg_outer_radius = std::max(radius + 1, data.bkg_outer_radius);
     const double bw = data.bandwidth;
 
-    auto recip_area = [&](double x, double y) -> double {
-        const Coord qx = data.geom.DetectorToRecip(x + 0.5, y) - data.geom.DetectorToRecip(x - 0.5, y);
-        const Coord qy = data.geom.DetectorToRecip(x, y + 0.5) - data.geom.DetectorToRecip(x, y - 0.5);
-        return (qx % qy).Length();
+    const auto pol_factor = experiment.GetPolarizationFactor();
+    auto polarization = [&](double x, double y) -> double {
+        if (!pol_factor)
+            return 1.0;
+        return data.geom.CalcAzIntPolarizationCorr(static_cast<float>(x), static_cast<float>(y),
+                                                   pol_factor.value());
     };
     auto bandwidth_radial_sq = [&](double d) -> double {
         if (bw <= 0.0 || d <= 0.0)
@@ -912,26 +970,9 @@ std::vector<float> PixelRefine::PredictImage(const AzimuthalIntegrationProfile &
         return bl * bl / (2.0 * d * d * d * d);
     };
 
-    // The model works in solid-angle/polarization-corrected units (as in Run,
-    // where Iobs = raw * correction). Map back to raw detector units (/ correction)
-    // so the predicted image overlays directly on the original image.
-    auto to_raw = [&](size_t npixel, double corrected) -> float {
-        const double corr = corrections[npixel];
-        return (corr > 0.0) ? static_cast<float>(corrected / corr) : 0.0f;
-    };
-
-    // Background base layer (per-pixel azimuthal mean), full-frame pass.
-    if (include_background) {
-        for (size_t p = 0; p < img.size(); ++p) {
-            const int bin = pixel_to_bin[p];
-            if (bin >= 0 && bin < total_bin_count)
-                img[p] = to_raw(p, azim_result[bin]);
-        }
-    }
-
-    double beam[2], dist_mm, detector_rot[2], rot_vec[3];
+    double beam[2], dist_mm, detector_rot[2];
     double latt_vec0[3], latt_vec1[3], latt_vec2[3];
-    BuildParameterBlocks(data, beam, dist_mm, detector_rot, rot_vec, latt_vec0, latt_vec1, latt_vec2);
+    BuildParameterBlocks(data, beam, dist_mm, detector_rot,  latt_vec0, latt_vec1, latt_vec2);
 
     DiffractionExperiment exp_iter = experiment;
     exp_iter.BeamX_pxl(data.geom.GetBeamX_pxl())
@@ -948,6 +989,7 @@ std::vector<float> PixelRefine::PredictImage(const AzimuthalIntegrationProfile &
     };
     const int nrefl = prediction.Calc(exp_iter, data.latt, settings_prediction);
     const auto &predicted = prediction.GetReflections();
+    const auto spot_mask = BuildSpotMask(predicted, nrefl, xpixel, ypixel, radius);
 
     for (int ri = 0; ri < nrefl; ++ri) {
         const auto &refl = predicted[ri];
@@ -957,6 +999,16 @@ std::vector<float> PixelRefine::PredictImage(const AzimuthalIntegrationProfile &
 
         const double Itrue = it->second;
         const double R_bw_sq = bandwidth_radial_sq(refl.d);
+        const double pol = polarization(refl.predicted_x, refl.predicted_y);
+
+        // Local background straight from the actual image (flat per shoebox), laid
+        // into the box so the prediction overlays the real frame - the same model
+        // path Run() fits, now reproduced faithfully because we have the image.
+        double Ibkg = 0.0;
+        const bool have_bkg = include_background &&
+            EstimateLocalBackground(image, spot_mask, xpixel, ypixel,
+                                    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);
@@ -967,25 +1019,21 @@ std::vector<float> PixelRefine::PredictImage(const AzimuthalIntegrationProfile &
             for (int x = min_x; x <= max_x; ++x) {
                 const size_t npixel = xpixel * y + x;
 
-                // Pure Bragg signal: Ibkg = 0 so Model() returns signal only; the
-                // background is already laid down above. Same code path as Run.
                 PixelObs obs{
                     .x = static_cast<double>(x),
                     .y = static_cast<double>(y),
                     .Iobs = 0.0,
-                    .Ibkg = 0.0,
-                    .weight = 1.0,
-                    .A_recip = recip_area(x, y),
-                    .angle_rad = angle_rad
+                    .Ibkg = have_bkg ? Ibkg : 0.0,
+                    .weight = 1.0
                 };
                 PixelResidual pr(obs, Itrue, lambda, pixel_size,
-                                 refl.h, refl.k, refl.l, R_bw_sq, data.crystal_system);
+                                 refl.h, refl.k, refl.l, R_bw_sq, pol, data.crystal_system);
 
-                double signal = 0.0;
-                if (pr.Model(beam, &dist_mm, detector_rot, rot_vec,
+                double Ipred = 0.0; // raw counts: signal (+ local background)
+                if (pr.Model(beam, &dist_mm, detector_rot,
                              latt_vec0, latt_vec1, latt_vec2,
-                             &data.scale_factor, &data.B_factor, data.R, signal))
-                    img[npixel] += to_raw(npixel, signal);
+                             &data.scale_factor, &data.B_factor, data.R, Ipred))
+                    img[npixel] += static_cast<float>(Ipred);
             }
         }
     }
@@ -995,26 +1043,22 @@ std::vector<float> PixelRefine::PredictImage(const AzimuthalIntegrationProfile &
 
 template<class T>
 std::vector<float> PixelRefine::ChiSquaredImage(const T *image,
-                                                const AzimuthalIntegrationProfile &profile,
                                                 BraggPrediction &prediction,
                                                 const PixelRefineData &data) const {
     std::vector<float> img(xpixel * ypixel, 0.0f);
 
     const double lambda = data.geom.GetWavelength_A();
     const double pixel_size = data.geom.GetPixelSize_mm();
-    const auto azim_result = profile.GetResult();
-    const auto azim_std = profile.GetStd();
-    const auto &pixel_to_bin = mapping.GetPixelToBin();
-    const auto &corrections = mapping.Corrections();
-    const int total_bin_count = static_cast<int>(azim_result.size());
-    const double angle_rad = data.angle_deg * M_PI / 180.0;
     const int radius = data.shoebox_radius;
+    const int bkg_outer_radius = std::max(radius + 1, data.bkg_outer_radius);
     const double bw = data.bandwidth;
 
-    auto recip_area = [&](double x, double y) -> double {
-        const Coord qx = data.geom.DetectorToRecip(x + 0.5, y) - data.geom.DetectorToRecip(x - 0.5, y);
-        const Coord qy = data.geom.DetectorToRecip(x, y + 0.5) - data.geom.DetectorToRecip(x, y - 0.5);
-        return (qx % qy).Length();
+    const auto pol_factor = experiment.GetPolarizationFactor();
+    auto polarization = [&](double x, double y) -> double {
+        if (!pol_factor)
+            return 1.0;
+        return data.geom.CalcAzIntPolarizationCorr(static_cast<float>(x), static_cast<float>(y),
+                                                   pol_factor.value());
     };
     auto bandwidth_radial_sq = [&](double d) -> double {
         if (bw <= 0.0 || d <= 0.0)
@@ -1023,9 +1067,9 @@ std::vector<float> PixelRefine::ChiSquaredImage(const T *image,
         return bl * bl / (2.0 * d * d * d * d);
     };
 
-    double beam[2], dist_mm, detector_rot[2], rot_vec[3];
+    double beam[2], dist_mm, detector_rot[2];
     double latt_vec0[3], latt_vec1[3], latt_vec2[3];
-    BuildParameterBlocks(data, beam, dist_mm, detector_rot, rot_vec, latt_vec0, latt_vec1, latt_vec2);
+    BuildParameterBlocks(data, beam, dist_mm, detector_rot,  latt_vec0, latt_vec1, latt_vec2);
 
     DiffractionExperiment exp_iter = experiment;
     exp_iter.BeamX_pxl(data.geom.GetBeamX_pxl())
@@ -1042,6 +1086,7 @@ std::vector<float> PixelRefine::ChiSquaredImage(const T *image,
     };
     const int nrefl = prediction.Calc(exp_iter, data.latt, settings_prediction);
     const auto &predicted = prediction.GetReflections();
+    const auto spot_mask = BuildSpotMask(predicted, nrefl, xpixel, ypixel, radius);
 
     for (int ri = 0; ri < nrefl; ++ri) {
         const auto &refl = predicted[ri];
@@ -1051,6 +1096,15 @@ std::vector<float> PixelRefine::ChiSquaredImage(const T *image,
 
         const double Itrue = it->second;
         const double R_bw_sq = bandwidth_radial_sq(refl.d);
+        const double pol = polarization(refl.predicted_x, refl.predicted_y);
+
+        // Local flat background, identical to Run(); skip the reflection if it
+        // cannot be estimated (matches Run() dropping the reflection).
+        double Ibkg = 0.0;
+        if (!EstimateLocalBackground(image, spot_mask, xpixel, ypixel,
+                                     refl.predicted_x, refl.predicted_y,
+                                     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);
@@ -1060,23 +1114,16 @@ std::vector<float> PixelRefine::ChiSquaredImage(const T *image,
         for (int y = min_y; y <= max_y; ++y) {
             for (int x = min_x; x <= max_x; ++x) {
                 const size_t npixel = xpixel * y + x;
-                const int azim_bin = pixel_to_bin[npixel];
 
                 // Same gating as Run(): only pixels that actually enter the fit.
-                if (azim_bin >= total_bin_count)
-                    continue;
                 if (image[npixel] == std::numeric_limits<T>::max())
                     continue;
                 if (std::is_signed_v<T> && (image[npixel] == std::numeric_limits<T>::min()))
                     continue;
 
-                const double correction = corrections[npixel];
-                const double Ibkg = azim_result[azim_bin];
-                const double Ibkg_sigma = azim_std[azim_bin];
-                const double raw = static_cast<double>(image[npixel]);
-                const double Iobs = raw * correction;
+                const double Iobs = static_cast<double>(image[npixel]); // raw counts
 
-                double var = correction * std::max(Iobs, 0.0) + Ibkg_sigma * Ibkg_sigma;
+                double var = std::max(Iobs, 0.0);
                 if (!(var > 1.0))
                     var = 1.0;
                 const double weight = 1.0 / std::sqrt(var);
@@ -1086,15 +1133,13 @@ std::vector<float> PixelRefine::ChiSquaredImage(const T *image,
                     .y = static_cast<double>(y),
                     .Iobs = Iobs,
                     .Ibkg = Ibkg,
-                    .weight = weight,
-                    .A_recip = recip_area(x, y),
-                    .angle_rad = angle_rad
+                    .weight = weight
                 };
                 PixelResidual pr(obs, Itrue, lambda, pixel_size,
-                                 refl.h, refl.k, refl.l, R_bw_sq, data.crystal_system);
+                                 refl.h, refl.k, refl.l, R_bw_sq, pol, data.crystal_system);
 
                 double Ipred = 0.0;
-                if (pr.Model(beam, &dist_mm, detector_rot, rot_vec,
+                if (pr.Model(beam, &dist_mm, detector_rot,
                              latt_vec0, latt_vec1, latt_vec2,
                              &data.scale_factor, &data.B_factor, data.R, Ipred)) {
                     // residual_i = (I_pred - I_obs) * weight (== Ceres residual);
@@ -1110,16 +1155,23 @@ std::vector<float> PixelRefine::ChiSquaredImage(const T *image,
 }
 
 // Explicit instantiations for the supported (uncompressed) image pixel types.
-template void PixelRefine::Run<int8_t>(const int8_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &);
-template void PixelRefine::Run<int16_t>(const int16_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &);
-template void PixelRefine::Run<int32_t>(const int32_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &);
-template void PixelRefine::Run<uint8_t>(const uint8_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &);
-template void PixelRefine::Run<uint16_t>(const uint16_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &);
-template void PixelRefine::Run<uint32_t>(const uint32_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &);
+template void PixelRefine::Run<int8_t>(const int8_t *, BraggPrediction &, PixelRefineData &);
+template void PixelRefine::Run<int16_t>(const int16_t *, BraggPrediction &, PixelRefineData &);
+template void PixelRefine::Run<int32_t>(const int32_t *, BraggPrediction &, PixelRefineData &);
+template void PixelRefine::Run<uint8_t>(const uint8_t *, BraggPrediction &, PixelRefineData &);
+template void PixelRefine::Run<uint16_t>(const uint16_t *, BraggPrediction &, PixelRefineData &);
+template void PixelRefine::Run<uint32_t>(const uint32_t *, BraggPrediction &, PixelRefineData &);
 
-template std::vector<float> PixelRefine::ChiSquaredImage<int8_t>(const int8_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, const PixelRefineData &) const;
-template std::vector<float> PixelRefine::ChiSquaredImage<int16_t>(const int16_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, const PixelRefineData &) const;
-template std::vector<float> PixelRefine::ChiSquaredImage<int32_t>(const int32_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, const PixelRefineData &) const;
-template std::vector<float> PixelRefine::ChiSquaredImage<uint8_t>(const uint8_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, const PixelRefineData &) const;
-template std::vector<float> PixelRefine::ChiSquaredImage<uint16_t>(const uint16_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, const PixelRefineData &) const;
-template std::vector<float> PixelRefine::ChiSquaredImage<uint32_t>(const uint32_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, const PixelRefineData &) const;
+template std::vector<float> PixelRefine::PredictImage<int8_t>(const int8_t *, BraggPrediction &, const PixelRefineData &, bool) const;
+template std::vector<float> PixelRefine::PredictImage<int16_t>(const int16_t *, BraggPrediction &, const PixelRefineData &, bool) const;
+template std::vector<float> PixelRefine::PredictImage<int32_t>(const int32_t *, BraggPrediction &, const PixelRefineData &, bool) const;
+template std::vector<float> PixelRefine::PredictImage<uint8_t>(const uint8_t *, BraggPrediction &, const PixelRefineData &, bool) const;
+template std::vector<float> PixelRefine::PredictImage<uint16_t>(const uint16_t *, BraggPrediction &, const PixelRefineData &, bool) const;
+template std::vector<float> PixelRefine::PredictImage<uint32_t>(const uint32_t *, BraggPrediction &, const PixelRefineData &, bool) const;
+
+template std::vector<float> PixelRefine::ChiSquaredImage<int8_t>(const int8_t *, BraggPrediction &, const PixelRefineData &) const;
+template std::vector<float> PixelRefine::ChiSquaredImage<int16_t>(const int16_t *, BraggPrediction &, const PixelRefineData &) const;
+template std::vector<float> PixelRefine::ChiSquaredImage<int32_t>(const int32_t *, BraggPrediction &, const PixelRefineData &) const;
+template std::vector<float> PixelRefine::ChiSquaredImage<uint8_t>(const uint8_t *, BraggPrediction &, const PixelRefineData &) const;
+template std::vector<float> PixelRefine::ChiSquaredImage<uint16_t>(const uint16_t *, BraggPrediction &, const PixelRefineData &) const;
+template std::vector<float> PixelRefine::ChiSquaredImage<uint32_t>(const uint32_t *, BraggPrediction &, const PixelRefineData &) const;