From bd5fef7f61c3da61cb34b51f9a86f332f3e81408 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Wed, 10 Jun 2026 18:36:12 +0200 Subject: [PATCH] PixelRefine: Simplify (remove Lorentz correction, remove background from azimuthal integration) --- image_analysis/IndexAndRefine.cpp | 32 +- image_analysis/IndexAndRefine.h | 10 +- image_analysis/MXAnalysisWithoutFPGA.cpp | 2 +- .../pixel_refinement/PixelRefine.cpp | 474 ++++++++++-------- image_analysis/pixel_refinement/PixelRefine.h | 52 +- viewer/JFJochImageReadingWorker.cpp | 18 +- 6 files changed, 316 insertions(+), 272 deletions(-) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 98b1aed6..4f3865c3 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -238,9 +238,7 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, const SpotFindingSettings &spot_finding_settings, const CompressedImage &image, BraggPrediction &prediction, - const IndexAndRefine::IndexingOutcome &outcome, - const AzimuthalIntegrationMapping *mapping, - const AzimuthalIntegrationProfile *profile) { + const IndexAndRefine::IndexingOutcome &outcome) { if (!outcome.lattice_candidate) return; @@ -299,11 +297,11 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, // integrated reflections that flow into the normal save/merge). const bool use_pixel_refine = experiment.GetIndexingSettings().GetGeomRefinementAlgorithm() == GeomRefinementAlgorithmEnum::PixelRefine - && !pixel_reference_.empty() && mapping && profile; + && !pixel_reference_.empty(); if (use_pixel_refine) { auto integration_start_time = std::chrono::steady_clock::now(); - PixelRefineIntegrate(msg, image, prediction, outcome, *mapping, *profile, i_outcome); + PixelRefineIntegrate(msg, image, prediction, outcome, i_outcome); msg.integrated_reflections = i_outcome.reflections.size(); auto integration_end_time = std::chrono::steady_clock::now(); msg.integration_time_s = std::chrono::duration(integration_end_time - integration_start_time).count(); @@ -356,9 +354,7 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, void IndexAndRefine::ProcessImage(DataMessage &msg, const SpotFindingSettings &spot_finding_settings, const CompressedImage &image, - BraggPrediction &prediction, - const AzimuthalIntegrationMapping *mapping, - const AzimuthalIntegrationProfile *profile) { + BraggPrediction &prediction) { if (!indexer_ || !spot_finding_settings.indexing) return; @@ -388,7 +384,7 @@ void IndexAndRefine::ProcessImage(DataMessage &msg, msg.lattice_type = outcome.symmetry; if (spot_finding_settings.quick_integration) - QuickPredictAndIntegrate(msg, spot_finding_settings, image, prediction, outcome, mapping, profile); + QuickPredictAndIntegrate(msg, spot_finding_settings, image, prediction, outcome); } std::optional IndexAndRefine::FinalizeRotationIndexing() { @@ -428,15 +424,13 @@ bool IndexAndRefine::PixelRefineIntegrate(DataMessage &msg, const CompressedImage &image, BraggPrediction &prediction, const IndexAndRefine::IndexingOutcome &outcome, - const AzimuthalIntegrationMapping &mapping, - const AzimuthalIntegrationProfile &profile, IntegrationOutcome &i_outcome) { if (!outcome.lattice_candidate) return false; - // Build the engine once (lazy: needs the azimuthal mapping, known only here). + // Build the engine once (lazy). std::call_once(pixel_refine_once_, [&] { - pixel_refine_ = std::make_unique(experiment, mapping, pixel_reference_); + pixel_refine_ = std::make_unique(experiment, pixel_reference_); }); if (!pixel_refine_) return false; @@ -455,17 +449,17 @@ bool IndexAndRefine::PixelRefineIntegrate(DataMessage &msg, const uint8_t *ptr = image.GetUncompressedPtr(buffer); switch (image.GetMode()) { case CompressedImageMode::Int8: - pixel_refine_->Run(reinterpret_cast(ptr), profile, prediction, prd); break; + pixel_refine_->Run(reinterpret_cast(ptr), prediction, prd); break; case CompressedImageMode::Int16: - pixel_refine_->Run(reinterpret_cast(ptr), profile, prediction, prd); break; + pixel_refine_->Run(reinterpret_cast(ptr), prediction, prd); break; case CompressedImageMode::Int32: - pixel_refine_->Run(reinterpret_cast(ptr), profile, prediction, prd); break; + pixel_refine_->Run(reinterpret_cast(ptr), prediction, prd); break; case CompressedImageMode::Uint8: - pixel_refine_->Run(reinterpret_cast(ptr), profile, prediction, prd); break; + pixel_refine_->Run(reinterpret_cast(ptr), prediction, prd); break; case CompressedImageMode::Uint16: - pixel_refine_->Run(reinterpret_cast(ptr), profile, prediction, prd); break; + pixel_refine_->Run(reinterpret_cast(ptr), prediction, prd); break; case CompressedImageMode::Uint32: - pixel_refine_->Run(reinterpret_cast(ptr), profile, prediction, prd); break; + pixel_refine_->Run(reinterpret_cast(ptr), prediction, prd); break; default: return false; } diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index 7674392e..b8b3e9f2 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -65,9 +65,7 @@ class IndexAndRefine { const SpotFindingSettings &spot_finding_settings, const CompressedImage &image, BraggPrediction &prediction, - const IndexingOutcome &outcome, - const AzimuthalIntegrationMapping *mapping, - const AzimuthalIntegrationProfile *profile); + const IndexingOutcome &outcome); std::unique_ptr scaling_engine; void ScaleImage(DataMessage &msg, IntegrationOutcome& outcome); @@ -83,8 +81,6 @@ class IndexAndRefine { const CompressedImage &image, BraggPrediction &prediction, const IndexingOutcome &outcome, - const AzimuthalIntegrationMapping &mapping, - const AzimuthalIntegrationProfile &profile, IntegrationOutcome &i_outcome); public: IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool *indexer); @@ -92,9 +88,7 @@ public: void AddImageToRotationIndexer(DataMessage &msg); void ForceRotationIndexerLattice(const CrystalLattice& lattice); - void ProcessImage(DataMessage &msg, const SpotFindingSettings &settings, const CompressedImage &image, BraggPrediction &prediction, - const AzimuthalIntegrationMapping *mapping = nullptr, - const AzimuthalIntegrationProfile *profile = nullptr); + void ProcessImage(DataMessage &msg, const SpotFindingSettings &settings, const CompressedImage &image, BraggPrediction &prediction); IndexAndRefine& ReferenceIntensities(std::vector &reference); ScalingResult ScaleAllImages(const std::vector &reference, size_t nthreads = 0); diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 6a638213..7bade492 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -92,7 +92,7 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, if (spot_finding_settings.indexing) indexer.ProcessImage(output, spot_finding_settings, CompressedImage(preprocessor_buffer->getBuffer(), experiment.GetXPixelsNum(), experiment.GetYPixelsNum()), - *prediction, &integration, &profile); + *prediction); } output.max_viable_pixel_value = ret.max_value; diff --git a/image_analysis/pixel_refinement/PixelRefine.cpp b/image_analysis/pixel_refinement/PixelRefine.cpp index 225778f7..f41925b3 100644 --- a/image_analysis/pixel_refinement/PixelRefine.cpp +++ b/image_analysis/pixel_refinement/PixelRefine.cpp @@ -3,6 +3,11 @@ #include "PixelRefine.h" +#include +#include +#include +#include + #include #include #include @@ -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 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 &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 BuildSpotMask(const std::vector &predicted, int nrefl, + size_t xpixel, size_t ypixel, int radius) { + std::vector mask(xpixel * ypixel, 0); + const double r_sq = static_cast(radius) * radius; + for (int i = 0; i < nrefl; ++i) { + const auto &r = predicted[i]; + const int cx = static_cast(std::lround(r.predicted_x)); + const int cy = static_cast(std::lround(r.predicted_y)); + const int x0 = std::max(0, cx - radius); + const int x1 = std::min(static_cast(xpixel) - 1, cx + radius); + const int y0 = std::max(0, cy - radius); + const int y1 = std::min(static_cast(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(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 +bool EstimateLocalBackground(const T *image, + const std::vector &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(std::lround(cx)); + const int icy = static_cast(std::lround(cy)); + const int x0 = std::max(0, icx - bkg_outer_radius); + const int x1 = std::min(static_cast(xpixel) - 1, icx + bkg_outer_radius); + const int y0 = std::max(0, icy - bkg_outer_radius); + const int y1 = std::min(static_cast(ypixel) - 1, icy + bkg_outer_radius); + + std::vector vals; + vals.reserve(static_cast((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(xpixel) * y + x; + if (spot_mask[np]) + continue; + const T raw = image[np]; + if (raw == std::numeric_limits::max()) + continue; + if (std::is_signed_v && raw == std::numeric_limits::min()) + continue; + vals.push_back(static_cast(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 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 &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(recip_obs[0], recip_obs[1], recip_obs[2]); + e_obs_recip = Eigen::Matrix(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 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 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 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 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 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 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 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 &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 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(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(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(x), static_cast(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::max()) continue; if (std::is_signed_v && (image[npixel] == std::numeric_limits::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(image[npixel]); - const double Iobs = raw * correction; + const double Iobs = static_cast(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(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(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(r.partiality) * B_term; // B & partiality + const double corr = static_cast(r.partiality) * B_term * g.pol; // B, partiality & pol r.bkg = static_cast(bkg_sum / static_cast(n)); r.observed = true; @@ -884,7 +943,8 @@ void PixelRefine::Run(const T *image, } } -std::vector PixelRefine::PredictImage(const AzimuthalIntegrationProfile &profile, +template +std::vector PixelRefine::PredictImage(const T *image, BraggPrediction &prediction, const PixelRefineData &data, bool include_background) const { @@ -892,18 +952,16 @@ std::vector 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(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(x), static_cast(y), + pol_factor.value()); }; auto bandwidth_radial_sq = [&](double d) -> double { if (bw <= 0.0 || d <= 0.0) @@ -912,26 +970,9 @@ std::vector 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(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 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 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(refl.predicted_y - radius, 0); const int max_y = std::min(refl.predicted_y + radius, ypixel - 1); @@ -967,25 +1019,21 @@ std::vector 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(x), .y = static_cast(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(Ipred); } } } @@ -995,26 +1043,22 @@ std::vector PixelRefine::PredictImage(const AzimuthalIntegrationProfile & template std::vector PixelRefine::ChiSquaredImage(const T *image, - const AzimuthalIntegrationProfile &profile, BraggPrediction &prediction, const PixelRefineData &data) const { std::vector 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(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(x), static_cast(y), + pol_factor.value()); }; auto bandwidth_radial_sq = [&](double d) -> double { if (bw <= 0.0 || d <= 0.0) @@ -1023,9 +1067,9 @@ std::vector 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 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 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(refl.predicted_y - radius, 0); const int max_y = std::min(refl.predicted_y + radius, ypixel - 1); @@ -1060,23 +1114,16 @@ std::vector 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::max()) continue; if (std::is_signed_v && (image[npixel] == std::numeric_limits::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(image[npixel]); - const double Iobs = raw * correction; + const double Iobs = static_cast(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 PixelRefine::ChiSquaredImage(const T *image, .y = static_cast(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 PixelRefine::ChiSquaredImage(const T *image, } // Explicit instantiations for the supported (uncompressed) image pixel types. -template void PixelRefine::Run(const int8_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &); -template void PixelRefine::Run(const int16_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &); -template void PixelRefine::Run(const int32_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &); -template void PixelRefine::Run(const uint8_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &); -template void PixelRefine::Run(const uint16_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &); -template void PixelRefine::Run(const uint32_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &); +template void PixelRefine::Run(const int8_t *, BraggPrediction &, PixelRefineData &); +template void PixelRefine::Run(const int16_t *, BraggPrediction &, PixelRefineData &); +template void PixelRefine::Run(const int32_t *, BraggPrediction &, PixelRefineData &); +template void PixelRefine::Run(const uint8_t *, BraggPrediction &, PixelRefineData &); +template void PixelRefine::Run(const uint16_t *, BraggPrediction &, PixelRefineData &); +template void PixelRefine::Run(const uint32_t *, BraggPrediction &, PixelRefineData &); -template std::vector PixelRefine::ChiSquaredImage(const int8_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, const PixelRefineData &) const; -template std::vector PixelRefine::ChiSquaredImage(const int16_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, const PixelRefineData &) const; -template std::vector PixelRefine::ChiSquaredImage(const int32_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, const PixelRefineData &) const; -template std::vector PixelRefine::ChiSquaredImage(const uint8_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, const PixelRefineData &) const; -template std::vector PixelRefine::ChiSquaredImage(const uint16_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, const PixelRefineData &) const; -template std::vector PixelRefine::ChiSquaredImage(const uint32_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, const PixelRefineData &) const; +template std::vector PixelRefine::PredictImage(const int8_t *, BraggPrediction &, const PixelRefineData &, bool) const; +template std::vector PixelRefine::PredictImage(const int16_t *, BraggPrediction &, const PixelRefineData &, bool) const; +template std::vector PixelRefine::PredictImage(const int32_t *, BraggPrediction &, const PixelRefineData &, bool) const; +template std::vector PixelRefine::PredictImage(const uint8_t *, BraggPrediction &, const PixelRefineData &, bool) const; +template std::vector PixelRefine::PredictImage(const uint16_t *, BraggPrediction &, const PixelRefineData &, bool) const; +template std::vector PixelRefine::PredictImage(const uint32_t *, BraggPrediction &, const PixelRefineData &, bool) const; + +template std::vector PixelRefine::ChiSquaredImage(const int8_t *, BraggPrediction &, const PixelRefineData &) const; +template std::vector PixelRefine::ChiSquaredImage(const int16_t *, BraggPrediction &, const PixelRefineData &) const; +template std::vector PixelRefine::ChiSquaredImage(const int32_t *, BraggPrediction &, const PixelRefineData &) const; +template std::vector PixelRefine::ChiSquaredImage(const uint8_t *, BraggPrediction &, const PixelRefineData &) const; +template std::vector PixelRefine::ChiSquaredImage(const uint16_t *, BraggPrediction &, const PixelRefineData &) const; +template std::vector PixelRefine::ChiSquaredImage(const uint32_t *, BraggPrediction &, const PixelRefineData &) const; diff --git a/image_analysis/pixel_refinement/PixelRefine.h b/image_analysis/pixel_refinement/PixelRefine.h index 73a63094..8d28e107 100644 --- a/image_analysis/pixel_refinement/PixelRefine.h +++ b/image_analysis/pixel_refinement/PixelRefine.h @@ -5,8 +5,6 @@ #include "../bragg_prediction/BraggPrediction.h" #include "../common/DiffractionExperiment.h" -#include "../common/AzimuthalIntegrationMapping.h" -#include "../common/AzimuthalIntegrationProfile.h" #include "../scale_merge/HKLKey.h" // ============================================================================= @@ -33,7 +31,10 @@ // squares problem. We write down, for every pixel in a reflection's shoebox, the // expected counts as an explicit forward model // -// I_pred(pixel) = G * I_true * B_term * P_radial * P_tangential + I_bkg +// I_pred(pixel) = G * I_true * B_term * P_radial * P_tangential * pol + I_bkg +// +// in raw detector counts (pol = per-reflection polarization correction, I_bkg = +// local per-shoebox background read from the image). // // and let Ceres autodiff back-propagate the per-pixel residuals into ALL of: // * detector geometry (beam centre, distance, tilt) @@ -95,8 +96,11 @@ // (bandwidth = 0, monochromatic); set it for DMM-type data, leave it for Si. // // Status: experimental prototype. The forward model (esp. the still-image -// Lorentz/partiality normalization) is deliberately simple and expected to -// evolve. See PixelRefine.cpp for the physics conventions and known caveats. +// partiality normalization) is deliberately simple and expected to evolve: it +// works in raw detector counts with a local per-shoebox background and a +// per-reflection polarization correction (no per-pixel solid-angle/Lorentz +// weighting), matching the classical integrator. See PixelRefine.cpp for the +// physics conventions and known caveats. // ============================================================================= struct PixelRefineData { @@ -115,23 +119,23 @@ struct PixelRefineData { // to R[0]. 0 = monochromatic (the term switches off entirely). double bandwidth = 0.0; - // Goniometer: for a still image keep angle_deg = 0. For a wedge, pass the - // angle (deg) of the slice centre relative to the reference (phi=0) frame. - double angle_deg = 0.0; - // --- what to refine --- bool refine_orientation = true; // crystal orientation (p0) bool refine_unit_cell = false; // cell lengths + angles bool refine_beam_center = false; bool refine_distance = false; bool refine_detector_angles = false; - bool refine_rotation_axis = false; bool refine_scale = true; bool refine_B = false; bool refine_R = true; double max_time_s = 5.0; - int shoebox_radius = 3; // half-size of the per-reflection pixel box + int shoebox_radius = 3; // half-size of the per-reflection signal box (peak region that enters the fit) + // Half-size of the local-background sampling box. Background is estimated from + // the ring shoebox_radius < |dx|,|dy| <= bkg_outer_radius around each spot + // (excluding pixels belonging to any predicted spot core), mirroring the local + // shoebox background of BraggIntegrate2D. Must be > shoebox_radius. + int bkg_outer_radius = 6; int max_iterations = 3; // inner predict<->refine cycles (re-predict with refined geom/latt) // --- output --- @@ -144,7 +148,6 @@ struct PixelRefineData { }; class PixelRefine { - const AzimuthalIntegrationMapping &mapping; const size_t xpixel, ypixel; const DiffractionExperiment &experiment; @@ -156,42 +159,43 @@ class PixelRefine { // PredictImage so both walk identical geometry/lattice code. void 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; public: PixelRefine(const DiffractionExperiment &experiment, - const AzimuthalIntegrationMapping &mapping, const std::vector &reference); // The BraggPrediction is supplied per call (it is mutated): this keeps a // single PixelRefine instance usable from several threads, each passing its // own prediction buffer. Only `data` is written; PixelRefine state is const. + // The image is in raw detector counts (masked/saturated pixels carry the type + // sentinel); background is estimated locally per shoebox from the image itself. template void Run(const T *image, - const AzimuthalIntegrationProfile &profile, BraggPrediction &prediction, PixelRefineData &data); - // Render the forward model as a full detector image (raw detector units, so + // Render the forward model as a full detector image (raw detector counts, so // it overlays directly on the original image). Uses the *same* per-pixel // model path (PixelResidual::Model) as the optimizer, evaluated in double // precision - slow but exact. For each reference reflection it adds the Bragg - // signal over its shoebox; with include_background it also lays down the - // azimuthal background. Diagnostic tool, not on the hot path. - std::vector PredictImage(const AzimuthalIntegrationProfile &profile, + // signal over its shoebox; with include_background it also lays down the local + // per-shoebox background read from the supplied image - the same background the + // fit uses. Diagnostic tool, not on the hot path. + template + std::vector PredictImage(const T *image, BraggPrediction &prediction, const PixelRefineData &data, bool include_background = true) const; // Render the per-pixel chi-square (cost density) that the optimizer actually // minimizes: for every shoebox pixel that enters the fit it stores the squared - // weighted residual ((I_pred - I_obs)/sigma)^2 in *corrected* units - identical - // to the Ceres residual_i^2 - accumulating where shoeboxes overlap. Pixels that - // are not part of any shoebox stay 0; masked/saturated pixels (skipped by the - // fit) also stay 0. Summing the image gives ~2*final_cost. Diagnostic tool. + // weighted residual ((I_pred - I_obs)/sigma)^2 in raw counts - identical to the + // Ceres residual_i^2 - accumulating where shoeboxes overlap. Pixels that are not + // part of any shoebox stay 0; masked/saturated pixels (skipped by the fit) also + // stay 0. Summing the image gives ~2*final_cost. Diagnostic tool. template std::vector ChiSquaredImage(const T *image, - const AzimuthalIntegrationProfile &profile, BraggPrediction &prediction, const PixelRefineData &data) const; }; diff --git a/viewer/JFJochImageReadingWorker.cpp b/viewer/JFJochImageReadingWorker.cpp index 8f90cb11..17a71934 100644 --- a/viewer/JFJochImageReadingWorker.cpp +++ b/viewer/JFJochImageReadingWorker.cpp @@ -457,8 +457,8 @@ void JFJochImageReadingWorker::ReanalyzeImage_i() { new_image_dataset->azimuthal_bins = azint_mapping->GetAzimuthalBinCount(); new_image_dataset->q_bins = azint_mapping->GetQBinCount(); - // Retain the profile (PixelRefine needs it). AzimuthalIntegrationProfile holds - // a mutex (non-copyable), so keep it via unique_ptr re-created each analysis. + // Azimuthal profile for the analysis/display pipeline. AzimuthalIntegrationProfile + // holds a mutex (non-copyable), so keep it via unique_ptr re-created each analysis. last_profile_ = std::make_unique(*azint_mapping); image_analysis->Analyze(new_image->ImageData(), *last_profile_, spot_finding_settings); @@ -725,8 +725,8 @@ void JFJochImageReadingWorker::LoadSpots(int64_t start_image, int64_t end_image, // Experimental PixelRefine // --------------------------------------------------------------------------- void JFJochImageReadingWorker::EnsurePixelRefine_i() { - if (!pixel_refine_ && azint_mapping && !pixel_reference_.empty()) - pixel_refine_ = std::make_unique(curr_experiment, *azint_mapping, pixel_reference_); + if (!pixel_refine_ && !pixel_reference_.empty()) + pixel_refine_ = std::make_unique(curr_experiment, pixel_reference_); if (!pixel_pred_) pixel_pred_ = CreateBraggPrediction(curr_experiment.IsRotationIndexing()); } @@ -854,15 +854,15 @@ QVector JFJochImageReadingWorker::BuildShoeboxes_i(const PixelRefineData std::vector JFJochImageReadingWorker::BuildDisplayImage_i(const PixelRefineData &data, int display_mode) const { + const auto &img32 = current_image_ptr->Image(); if (display_mode == PixelRefineParams::ChiSquared) { // The cost density the optimizer actually minimizes (weighted residual^2). - const auto &img32 = current_image_ptr->Image(); - auto chi2 = pixel_refine_->ChiSquaredImage(img32.data(), *last_profile_, *pixel_pred_, data); + auto chi2 = pixel_refine_->ChiSquaredImage(img32.data(), *pixel_pred_, data); MaskMeasuredSentinels_i(chi2); return chi2; } - auto pred = pixel_refine_->PredictImage(*last_profile_, *pixel_pred_, data, true); + auto pred = pixel_refine_->PredictImage(img32.data(), *pixel_pred_, data, true); if (display_mode == PixelRefineParams::SquaredDifference) SquaredResidualWithImage_i(pred); return pred; @@ -907,7 +907,7 @@ void JFJochImageReadingWorker::PixelRefinePreview(PixelRefineParams params) { try { const auto &img32 = current_image_ptr->Image(); - pixel_refine_->Run(img32.data(), *last_profile_, *pixel_pred_, d); + pixel_refine_->Run(img32.data(), *pixel_pred_, d); emit pixelRefineResidual(d.final_cost, d.cc, static_cast(d.reflections.size())); auto display = BuildDisplayImage_i(d, params.display_mode); @@ -941,7 +941,7 @@ void JFJochImageReadingWorker::PixelRefineRun(PixelRefineParams params) { try { const auto &img32 = current_image_ptr->Image(); - pixel_refine_->Run(img32.data(), *last_profile_, *pixel_pred_, d); + pixel_refine_->Run(img32.data(), *pixel_pred_, d); // Push refined values back so the sliders follow the optimizer. PixelRefineParams out = params;