jfjoch_process: Remove postrefinement

common/DatasetSettings.cpp

@@ -389,6 +389,15 @@ std::optional<float> DatasetSettings::GetPolarizationFactor() const {
     return polarization_factor;
 }
 
+DatasetSettings &DatasetSettings::BandwidthFWHM(const std::optional<float> &input) {
+    bandwidth_fwhm = input;
+    return *this;
+}
+
+std::optional<float> DatasetSettings::GetBandwidthFWHM() const {
+    return bandwidth_fwhm;
+}
+
 float DatasetSettings::GetPoniRot3_rad() const {
     return poni_rot_3_rad;
 }

common/DatasetSettings.h

@@ -55,6 +55,7 @@ class DatasetSettings {
     bool write_nxmx_hdf5_master;
 
     std::optional<float> polarization_factor;
+    std::optional<float> bandwidth_fwhm; // relative X-ray bandwidth, FWHM of dlambda/lambda (e.g. 0.01 for 1%)
     float poni_rot_1_rad;
     float poni_rot_2_rad;
     float poni_rot_3_rad;
@@ -99,6 +100,7 @@ public:
     DatasetSettings& PixelValueLowThreshold(const std::optional<int64_t> &input);
     DatasetSettings& PixelValueHighThreshold(const std::optional<int64_t> &input);
     DatasetSettings& PolarizationFactor(const std::optional<float> &input);
+    DatasetSettings& BandwidthFWHM(const std::optional<float> &input);
     DatasetSettings& PoniRot1_rad(float input);
     DatasetSettings& PoniRot2_rad(float input);
     DatasetSettings& PoniRot3_rad(float input);
@@ -150,6 +152,7 @@ public:
     std::optional<bool> IsSaveCalibration() const;
 
     std::optional<float> GetPolarizationFactor() const;
+    std::optional<float> GetBandwidthFWHM() const;
     float GetPoniRot1_rad() const;
     float GetPoniRot2_rad() const;
     float GetPoniRot3_rad() const;

common/DiffractionExperiment.cpp

@@ -641,6 +641,11 @@ void DiffractionExperiment::FillMessage(StartMessage &message) const {
     message.beam_center_y = GetBeamY_pxl();
     message.detector_distance = GetDetectorDistance_mm() * 1e-3f;
     message.incident_wavelength = GetWavelength_A();
+    // NXmx incident_wavelength_spread is the FWHM of the wavelength distribution,
+    // in wavelength units. We store bandwidth as a relative FWHM (dlambda/lambda),
+    // so convert to absolute: dlambda = (dlambda/lambda) * lambda.
+    if (const auto bw = GetBandwidthFWHM())
+        message.incident_wavelength_spread = bw.value() * GetWavelength_A();
     message.incident_energy = GetIncidentEnergy_keV() * 1e3f;
     message.image_size_x = GetXPixelsNum();
     message.image_size_y = GetYPixelsNum();
@@ -1310,6 +1315,15 @@ std::optional<float> DiffractionExperiment::GetPolarizationFactor() const {
     return dataset.GetPolarizationFactor();
 }
 
+DiffractionExperiment &DiffractionExperiment::BandwidthFWHM(const std::optional<float> &input) {
+    dataset.BandwidthFWHM(input);
+    return *this;
+}
+
+std::optional<float> DiffractionExperiment::GetBandwidthFWHM() const {
+    return dataset.GetBandwidthFWHM();
+}
+
 DiffractionExperiment &DiffractionExperiment::SaveCalibration(const std::optional<bool> &input) {
     dataset.SaveCalibration(input);
     return *this;

common/DiffractionExperiment.h

@@ -281,9 +281,11 @@ public:
 
     DiffractionExperiment& ApplySolidAngleCorr(bool input);
     DiffractionExperiment& PolarizationFactor(const std::optional<float> &input);
+    DiffractionExperiment& BandwidthFWHM(const std::optional<float> &input);
 
     bool GetApplySolidAngleCorr() const;
     std::optional<float> GetPolarizationFactor() const;
+    std::optional<float> GetBandwidthFWHM() const;
 
     int64_t GetUDPInterfaceCount() const;
     std::vector<DetectorModuleConfig> GetDetectorModuleConfig(const std::vector<AcquisitionDeviceNetConfig>& net_config) const;

common/IndexingSettings.h

@@ -6,7 +6,7 @@
 #include <cstdint>
 
 enum class IndexingAlgorithmEnum {FFBIDX, FFT, FFTW, Auto, None};
-enum class GeomRefinementAlgorithmEnum {None, OrientationOnly, BeamCenter};
+enum class GeomRefinementAlgorithmEnum {None, OrientationOnly, BeamCenter, PixelRefine};
 
 class IndexingSettings {
     IndexingAlgorithmEnum algorithm;

common/JFJochMessages.h

@@ -208,6 +208,7 @@ struct StartMessage {
 
     float incident_energy;
     float incident_wavelength;
+    std::optional<float> incident_wavelength_spread; // NXmx incident_wavelength_spread: FWHM of dlambda (Angstrom)
 
     float frame_time;
     float count_time;

common/ScalingSettings.h

@@ -6,7 +6,7 @@
 #include <optional>
 #include "JFJochException.h"
 
-enum class PartialityModel { Fixed, Rotation, Unity, Postrefinement };
+enum class PartialityModel { Fixed, Rotation, Unity };
 enum class IntensityFormat { Text, mmCIF, MTZ};
 
 class ScalingSettings {

docs/CBOR.md

@@ -27,6 +27,7 @@ There are minor differences at the moment:
 | image_size_y                     | uint64               | Image height \[pixels\]                                                                                                                        |             X             |
 | incident_energy                  | float                | X-ray energy \[eV\]                                                                                                                            |             X             |
 | incident_wavelength              | float                | X-ray wavelength \[Angstrom\]                                                                                                                  |             X             |
+| incident_wavelength_spread       | float (optional)     | FWHM of the X-ray wavelength distribution \[Angstrom\] (NXmx incident_wavelength_spread); omitted when the beam is monochromatic               |                           |
 | frame_time                       | float                | Frame time, if multiple frames per trigger \[s\]                                                                                               |             X             |
 | count_time                       | float                | Exposure time \[s\]                                                                                                                            |             X             |
 | saturation_value                 | int64                | Maximum valid sample value                                                                                                                     |             X             |

frame_serialize/CBORStream2Deserializer.cpp

@@ -1182,6 +1182,8 @@ namespace {
                 message.incident_energy = GetCBORFloat(value);
             else if (key == "incident_wavelength")
                 message.incident_wavelength = GetCBORFloat(value);
+            else if (key == "incident_wavelength_spread")
+                message.incident_wavelength_spread = GetCBORFloat(value);
             else if (key == "frame_time")
                 message.frame_time = GetCBORFloat(value);
             else if (key == "count_time")

frame_serialize/CBORStream2Serializer.cpp

@@ -632,6 +632,7 @@ void CBORStream2Serializer::SerializeSequenceStart(const StartMessage& message)
 
     CBOR_ENC(mapEncoder, "incident_energy", message.incident_energy);
     CBOR_ENC(mapEncoder, "incident_wavelength", message.incident_wavelength);
+    CBOR_ENC(mapEncoder, "incident_wavelength_spread", message.incident_wavelength_spread);
 
     CBOR_ENC(mapEncoder, "frame_time", message.frame_time);
     CBOR_ENC(mapEncoder, "count_time", message.count_time);

image_analysis/CMakeLists.txt

@@ -50,5 +50,6 @@ ADD_SUBDIRECTORY(lattice_search)
 ADD_SUBDIRECTORY(scale_merge)
 ADD_SUBDIRECTORY(image_preprocessing)
 ADD_SUBDIRECTORY(azint)
+ADD_SUBDIRECTORY(pixel_refinement)
 
-TARGET_LINK_LIBRARIES(JFJochImageAnalysis JFJochAzIntEngine JFJochImagePreprocessing JFJochBraggPrediction JFJochBraggIntegration JFJochLatticeSearch JFJochIndexing JFJochSpotFinding JFJochCommon JFJochGeomRefinement JFJochScaleMerge gemmi)
+TARGET_LINK_LIBRARIES(JFJochImageAnalysis JFJochAzIntEngine JFJochImagePreprocessing JFJochBraggPrediction JFJochBraggIntegration JFJochLatticeSearch JFJochIndexing JFJochSpotFinding JFJochCommon JFJochGeomRefinement JFJochScaleMerge JFJochPixelRefine gemmi)

image_analysis/IndexAndRefine.cpp

@@ -175,6 +175,10 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In
             XtalOptimizerRotationOnly(data, msg.spots, 0.05);
             break;
         case GeomRefinementAlgorithmEnum::BeamCenter:
+        case GeomRefinementAlgorithmEnum::PixelRefine:
+            // PixelRefine still benefits from the classical beam-center + cell
+            // refinement as a starting point; the pixel-level refinement runs
+            // later, during integration.
             if (XtalOptimizer(data, {msg.spots})) {
                 outcome.experiment.BeamX_pxl(data.geom.GetBeamX_pxl())
                         .BeamY_pxl(data.geom.GetBeamY_pxl());
@@ -234,7 +238,9 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg,
                                               const SpotFindingSettings &spot_finding_settings,
                                               const CompressedImage &image,
                                               BraggPrediction &prediction,
-                                              const IndexAndRefine::IndexingOutcome &outcome) {
+                                              const IndexAndRefine::IndexingOutcome &outcome,
+                                              const AzimuthalIntegrationMapping *mapping,
+                                              const AzimuthalIntegrationProfile *profile) {
     if (!outcome.lattice_candidate)
         return;
 
@@ -283,17 +289,36 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg,
         .max_hkl = 100,
         .centering = outcome.symmetry.centering,
         .wedge_deg = std::fabs(wedge_deg),
-        .mosaicity_deg = std::fabs(mos_deg)
+        .mosaicity_deg = std::fabs(mos_deg),
+        // FWHM -> sigma; 0 when monochromatic, leaving the prediction unchanged.
+        .bandwidth_sigma = experiment.GetBandwidthFWHM().value_or(0.0f) / 2.3548f
     };
 
-    auto pred_start_time = std::chrono::steady_clock::now();
-    auto nrefl = prediction.Calc(outcome.experiment, latt, settings_prediction);
-    auto pred_end_time = std::chrono::steady_clock::now();
-    msg.bragg_prediction_time_s = std::chrono::duration<float>(pred_end_time - pred_start_time).count();
+    // Select the integration path: classical 2D integration, or the experimental
+    // PixelRefine joint geometry/profile/scale refinement (which also produces the
+    // integrated reflections that flow into the normal save/merge).
+    const bool use_pixel_refine =
+        experiment.GetIndexingSettings().GetGeomRefinementAlgorithm() == GeomRefinementAlgorithmEnum::PixelRefine
+        && !pixel_reference_.empty() && mapping && profile;
 
-    auto integration_start_time = std::chrono::steady_clock::now();
-    i_outcome.reflections = BraggIntegrate2D(outcome.experiment, image, prediction.GetReflections(), nrefl, msg.number);
-    msg.integrated_reflections = i_outcome.reflections.size();
+    if (use_pixel_refine) {
+        auto integration_start_time = std::chrono::steady_clock::now();
+        PixelRefineIntegrate(msg, image, prediction, outcome, *mapping, *profile, 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<float>(integration_end_time - integration_start_time).count();
+    } else {
+        auto pred_start_time = std::chrono::steady_clock::now();
+        auto nrefl = prediction.Calc(outcome.experiment, latt, settings_prediction);
+        auto pred_end_time = std::chrono::steady_clock::now();
+        msg.bragg_prediction_time_s = std::chrono::duration<float>(pred_end_time - pred_start_time).count();
+
+        auto integration_start_time = std::chrono::steady_clock::now();
+        i_outcome.reflections = BraggIntegrate2D(outcome.experiment, image, prediction.GetReflections(), nrefl, msg.number);
+        msg.integrated_reflections = i_outcome.reflections.size();
+        auto integration_end_time = std::chrono::steady_clock::now();
+        msg.integration_time_s = std::chrono::duration<float>(integration_end_time - integration_start_time).count();
+    }
 
     constexpr size_t kMaxReflections = 10000;
     if (i_outcome.reflections.size() > kMaxReflections) {
@@ -314,10 +339,10 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg,
     CalcISigma(msg, i_outcome.reflections);
     CalcWilsonBFactor(msg, i_outcome.reflections);
 
-    auto integration_end_time = std::chrono::steady_clock::now();
-    msg.integration_time_s = std::chrono::duration<float>(integration_end_time - integration_start_time).count();
-
-    ScaleImage(msg, i_outcome);
+    // PixelRefine produces already-scaled reflections; only the classical path
+    // needs the separate ScaleOnTheFly step.
+    if (!use_pixel_refine)
+        ScaleImage(msg, i_outcome);
 
     // Copy reflections to outgoing message
     msg.reflections = i_outcome.reflections;
@@ -331,7 +356,9 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg,
 void IndexAndRefine::ProcessImage(DataMessage &msg,
                                   const SpotFindingSettings &spot_finding_settings,
                                   const CompressedImage &image,
-                                  BraggPrediction &prediction) {
+                                  BraggPrediction &prediction,
+                                  const AzimuthalIntegrationMapping *mapping,
+                                  const AzimuthalIntegrationProfile *profile) {
     if (!indexer_ || !spot_finding_settings.indexing)
         return;
 
@@ -361,7 +388,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);
+        QuickPredictAndIntegrate(msg, spot_finding_settings, image, prediction, outcome, mapping, profile);
 }
 
 std::optional<RotationIndexerResult> IndexAndRefine::FinalizeRotationIndexing() {
@@ -377,6 +404,7 @@ std::optional<RotationIndexerResult> IndexAndRefine::FinalizeRotationIndexing()
 
 IndexAndRefine &IndexAndRefine::ReferenceIntensities(std::vector<MergedReflection> &reference) {
     scaling_engine = std::make_unique<ScaleOnTheFly>(experiment, reference);
+    pixel_reference_ = reference; // kept for the experimental PixelRefine path
     return *this;
 }
 
@@ -396,6 +424,72 @@ void IndexAndRefine::ScaleImage(DataMessage &msg, IntegrationOutcome& outcome) {
     msg.image_scale_time_s = std::chrono::duration<float>(scaling_end_time - scaling_start_time).count();
 }
 
+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).
+    std::call_once(pixel_refine_once_, [&] {
+        pixel_refine_ = std::make_unique<PixelRefine>(experiment, mapping, pixel_reference_);
+    });
+    if (!pixel_refine_)
+        return false;
+
+    PixelRefineData prd;
+    prd.geom = outcome.experiment.GetDiffractionGeometry();
+    prd.latt = *outcome.lattice_candidate;
+    prd.crystal_system = outcome.symmetry.crystal_system;
+    if (prd.crystal_system == gemmi::CrystalSystem::Trigonal)
+        prd.crystal_system = gemmi::CrystalSystem::Hexagonal;
+    prd.centering = outcome.symmetry.centering;
+    if (const auto bw = experiment.GetBandwidthFWHM())
+        prd.bandwidth = bw.value() / 2.3548; // FWHM -> sigma
+
+    std::vector<uint8_t> buffer;
+    const uint8_t *ptr = image.GetUncompressedPtr(buffer);
+    switch (image.GetMode()) {
+        case CompressedImageMode::Int8:
+            pixel_refine_->Run(reinterpret_cast<const int8_t *>(ptr), profile, prediction, prd); break;
+        case CompressedImageMode::Int16:
+            pixel_refine_->Run(reinterpret_cast<const int16_t *>(ptr), profile, prediction, prd); break;
+        case CompressedImageMode::Int32:
+            pixel_refine_->Run(reinterpret_cast<const int32_t *>(ptr), profile, prediction, prd); break;
+        case CompressedImageMode::Uint8:
+            pixel_refine_->Run(reinterpret_cast<const uint8_t *>(ptr), profile, prediction, prd); break;
+        case CompressedImageMode::Uint16:
+            pixel_refine_->Run(reinterpret_cast<const uint16_t *>(ptr), profile, prediction, prd); break;
+        case CompressedImageMode::Uint32:
+            pixel_refine_->Run(reinterpret_cast<const uint32_t *>(ptr), profile, prediction, prd); break;
+        default:
+            return false;
+    }
+
+    // PixelRefine output flows into the normal save/merge path: the refined
+    // geometry/lattice and the already-scaled reflections become the outcome.
+    i_outcome.reflections = std::move(prd.reflections);
+    i_outcome.geom = prd.geom;
+    i_outcome.latt = prd.latt;
+    i_outcome.image_scale_g = static_cast<float>(prd.scale_factor);
+    i_outcome.image_scale_b_factor_Ang2 = static_cast<float>(prd.B_factor);
+    if (std::isfinite(prd.cc)) {
+        i_outcome.image_scale_cc = static_cast<float>(prd.cc);
+        i_outcome.image_scale_cc_n = prd.cc_n;
+    }
+
+    msg.image_scale_factor = static_cast<float>(prd.scale_factor);
+    msg.image_scale_cc = i_outcome.image_scale_cc;
+    if (prd.B_factor != 0.0)
+        msg.image_scale_b_factor = static_cast<float>(prd.B_factor);
+
+    return true;
+}
+
 ScalingResult IndexAndRefine::ScaleAllImages(const std::vector<MergedReflection> &reference, size_t nthreads) {
     ScaleOnTheFly scaling(experiment, reference);
     scaling.Scale(integration_outcome, nthreads);

image_analysis/IndexAndRefine.h

@@ -8,7 +8,10 @@
 
 #include "../common/DiffractionSpot.h"
 #include "../common/DiffractionExperiment.h"
+#include "../common/AzimuthalIntegrationMapping.h"
+#include "../common/AzimuthalIntegrationProfile.h"
 #include "bragg_prediction/BraggPrediction.h"
+#include "pixel_refinement/PixelRefine.h"
 #include "indexing/IndexerThreadPool.h"
 #include "lattice_search/LatticeSearch.h"
 #include "rotation_indexer/RotationIndexer.h"
@@ -62,17 +65,36 @@ class IndexAndRefine {
                                  const SpotFindingSettings &spot_finding_settings,
                                  const CompressedImage &image,
                                  BraggPrediction &prediction,
-                                 const IndexingOutcome &outcome);
+                                 const IndexingOutcome &outcome,
+                                 const AzimuthalIntegrationMapping *mapping,
+                                 const AzimuthalIntegrationProfile *profile);
 
     std::unique_ptr<ScaleOnTheFly> scaling_engine;
     void ScaleImage(DataMessage &msg, IntegrationOutcome& outcome);
+
+    // Experimental PixelRefine integration path (selected via
+    // GeomRefinementAlgorithmEnum::PixelRefine). Needs reference intensities; the
+    // engine is built lazily on first use (when the azimuthal mapping is known)
+    // and is safe to share across threads (prediction is supplied per call).
+    std::vector<MergedReflection> pixel_reference_;
+    std::unique_ptr<PixelRefine> pixel_refine_;
+    std::once_flag pixel_refine_once_;
+    bool PixelRefineIntegrate(DataMessage &msg,
+                              const CompressedImage &image,
+                              BraggPrediction &prediction,
+                              const IndexingOutcome &outcome,
+                              const AzimuthalIntegrationMapping &mapping,
+                              const AzimuthalIntegrationProfile &profile,
+                              IntegrationOutcome &i_outcome);
 public:
     IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool *indexer);
 
     void AddImageToRotationIndexer(DataMessage &msg);
     void ForceRotationIndexerLattice(const CrystalLattice& lattice);
 
-    void ProcessImage(DataMessage &msg, const SpotFindingSettings &settings, const CompressedImage &image, BraggPrediction &prediction);
+    void ProcessImage(DataMessage &msg, const SpotFindingSettings &settings, const CompressedImage &image, BraggPrediction &prediction,
+                      const AzimuthalIntegrationMapping *mapping = nullptr,
+                      const AzimuthalIntegrationProfile *profile = nullptr);
     IndexAndRefine& ReferenceIntensities(std::vector<MergedReflection> &reference);
 
     ScalingResult ScaleAllImages(const std::vector<MergedReflection> &reference, size_t nthreads = 0);

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);
+                *prediction, &integration, &profile);
     }
 
     output.max_viable_pixel_value = ret.max_value;

image_analysis/UpdateReflectionResolution.h

@@ -9,8 +9,10 @@
 #include "IntegrationOutcome.h"
 
 struct ResolutionStats {
-    float d_low = std::numeric_limits<float>::max();
-    float d_high = 0.0f;
+    // d_high = highest resolution = smallest d (tracked via `d_high > d`);
+    // d_low  = lowest  resolution = largest  d (tracked via `d_low  < d`).
+    float d_low = 0.0f;
+    float d_high = std::numeric_limits<float>::max();
     int n_reflections = 0;
     int n_images = 0;
 };

image_analysis/bragg_prediction/BraggPrediction.cpp

@@ -4,6 +4,13 @@
 #include "BraggPrediction.h"
 #include "../bragg_integration/SystematicAbsence.h"
 
+namespace {
+    // Number of bandwidth sigmas included in the (radially thickened) Ewald-shell
+    // acceptance window. 3σ captures essentially the whole pink-beam smear; matches
+    // the conservative end of the mosaicity cutoff used by callers.
+    constexpr float kBandwidthCutoffSigmas = 3.0f;
+}
+
 BraggPrediction::BraggPrediction(int max_reflections)
     : max_reflections(max_reflections), reflections(max_reflections) {}
 
@@ -76,7 +83,18 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal
                 float S_len = sqrtf(S_x * S_x + S_y * S_y + S_z * S_z);
                 float dist_ewald_sphere = std::fabs(S_len - one_over_wavelength);
 
-                if (dist_ewald_sphere <= settings.ewald_dist_cutoff ) {
+                // Energy bandwidth thickens the Ewald shell radially: at the
+                // diffraction condition |S|-1/λ shifts by recip_z·(Δλ/λ), i.e.
+                // σ_bw = |recip_z|·bandwidth_sigma (= bλ/2d², identical to
+                // PixelRefine's R_bw). Broaden the acceptance window in quadrature so
+                // high-resolution shells (smeared most, ∝1/d²) are not clipped.
+                float radial_cutoff = settings.ewald_dist_cutoff;
+                if (settings.bandwidth_sigma > 0.0f) {
+                    const float bw_tol = kBandwidthCutoffSigmas * settings.bandwidth_sigma * std::fabs(recip_z);
+                    radial_cutoff = std::sqrt(radial_cutoff * radial_cutoff + bw_tol * bw_tol);
+                }
+
+                if (dist_ewald_sphere <= radial_cutoff ) {
                     const float s0_p0 = S0.x * recip_x + S0.y * recip_y + S0.z * recip_z;
                     const float val = s0_sq * recip_sq - s0_p0 * s0_p0;
 

image_analysis/bragg_prediction/BraggPrediction.h

@@ -18,6 +18,12 @@ struct BraggPredictionSettings {
     float mosaicity_deg = 0.2f;
     float min_zeta = 0.05;
     float mosaicity_multiplier = 4.0;
+    // Relative X-ray bandwidth Δλ/λ expressed as a Gaussian sigma (0 = monochromatic).
+    // When > 0 the Ewald-shell acceptance is thickened radially per reflection by
+    // σ_bw = |recip_z|·bandwidth_sigma (= bλ/2d²), so the 1/d² pink-beam smear no
+    // longer clips high-resolution reflections. Must stay the last member so existing
+    // designated initializers remain valid.
+    float bandwidth_sigma = 0.0f;
 };
 
 class BraggPrediction {

image_analysis/bragg_prediction/BraggPredictionGPU.cu

@@ -8,6 +8,10 @@
 #include <cuda_runtime.h>
 
 namespace {
+    // Number of bandwidth sigmas included in the (radially thickened) Ewald-shell
+    // acceptance window. Mirrors the CPU BraggPrediction path.
+    constexpr float kBandwidthCutoffSigmas = 3.0f;
+
     __device__ inline bool is_odd(int v) { return (v & 1) != 0; }
 
     __device__ inline float angle_from_ewald_sphere_deg(const Coord &S0, float recip_x, float recip_y, float recip_z, float recip_sq) {
@@ -95,7 +99,14 @@ namespace {
         float Sz = recip_z + C.S0.z;
         float S_len = sqrtf(Sx * Sx + Sy * Sy + Sz * Sz);
         float dist_ewald = fabsf(S_len - C.one_over_wavelength);
-        if (dist_ewald > C.ewald_cutoff) return false;
+        // Energy bandwidth thickens the Ewald shell radially: σ_bw = |recip_z|·(Δλ/λ)
+        // (= bλ/2d²). Broaden the acceptance window in quadrature (see CPU path).
+        float radial_cutoff = C.ewald_cutoff;
+        if (C.bandwidth_sigma > 0.0f) {
+            const float bw_tol = kBandwidthCutoffSigmas * C.bandwidth_sigma * fabsf(recip_z);
+            radial_cutoff = sqrtf(radial_cutoff * radial_cutoff + bw_tol * bw_tol);
+        }
+        if (dist_ewald > radial_cutoff) return false;
         float Srx = C.rot[0] * Sx + C.rot[1] * Sy + C.rot[2] * Sz;
         float Sry = C.rot[3] * Sx + C.rot[4] * Sy + C.rot[5] * Sz;
         float Srz = C.rot[6] * Sx + C.rot[7] * Sy + C.rot[8] * Sz;
@@ -145,7 +156,8 @@ namespace {
                                           const CrystalLattice &lattice,
                                           float high_res_A,
                                           float ewald_dist_cutoff,
-                                          char centering) {
+                                          char centering,
+                                          float bandwidth_sigma) {
         KernelConsts kc{};
         auto geom = experiment.GetDiffractionGeometry();
         kc.det_width_pxl = static_cast<float>(experiment.GetXPixelsNum());
@@ -157,6 +169,7 @@ namespace {
         kc.one_over_dmax_sq = one_over_dmax * one_over_dmax;
         kc.one_over_wavelength = 1.0f / geom.GetWavelength_A();
         kc.ewald_cutoff = ewald_dist_cutoff;
+        kc.bandwidth_sigma = bandwidth_sigma;
         kc.Astar = lattice.Astar();
         kc.Bstar = lattice.Bstar();
         kc.Cstar = lattice.Cstar();
@@ -178,7 +191,8 @@ int BraggPredictionGPU::Calc(const DiffractionExperiment &experiment,
                              const CrystalLattice &lattice,
                              const BraggPredictionSettings &settings) {
     // Build constants on host
-    KernelConsts hK = BuildKernelConsts(experiment, lattice, settings.high_res_A, settings.ewald_dist_cutoff, settings.centering);
+    KernelConsts hK = BuildKernelConsts(experiment, lattice, settings.high_res_A, settings.ewald_dist_cutoff,
+                                        settings.centering, settings.bandwidth_sigma);
     cudaMemcpyAsync(dK, &hK, sizeof(KernelConsts), cudaMemcpyHostToDevice, stream);
     cudaMemsetAsync(d_count, 0, sizeof(int), stream);
 

image_analysis/bragg_prediction/BraggPredictionGPU.h

@@ -18,6 +18,7 @@ struct KernelConsts {
     float one_over_wavelength;
     float one_over_dmax_sq;
     float ewald_cutoff;
+    float bandwidth_sigma;   // relative Δλ/λ (sigma); 0 = monochromatic
     Coord Astar, Bstar, Cstar, S0;
     float rot[9];
     char centering;

image_analysis/geom_refinement/CMakeLists.txt

@@ -6,6 +6,8 @@ ADD_LIBRARY(JFJochGeomRefinement STATIC
         AssignSpotsToRings.h
         XtalOptimizer.cpp
         XtalOptimizer.h
+        LatticeReduction.cpp
+        LatticeReduction.h
 )
 
 TARGET_LINK_LIBRARIES(JFJochGeomRefinement Ceres::ceres Eigen3::Eigen JFJochCommon)

image_analysis/geom_refinement/LatticeReduction.cpp

@@ -0,0 +1,257 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "LatticeReduction.h"
+
+#include <Eigen/Dense>
+
+void LatticeToRodriguesAndLengths_GS(const CrystalLattice &latt,
+                                            double rod[3],
+                                            double lengths[3]) {
+    // Load lattice columns
+    const Coord a = latt.Vec0();
+    const Coord b = latt.Vec1();
+    const Coord c = latt.Vec2();
+
+    Eigen::Vector3d A(a[0], a[1], a[2]);
+    Eigen::Vector3d B(b[0], b[1], b[2]);
+    Eigen::Vector3d C(c[0], c[1], c[2]);
+
+    // Lengths = column norms (orthorhombic assumption)
+    lengths[0] = A.norm();
+    lengths[1] = B.norm();
+    lengths[2] = C.norm();
+
+    auto safe_unit = [](const Eigen::Vector3d &v, double eps = 1e-15) -> Eigen::Vector3d {
+        double n = v.norm();
+        return (n > eps) ? (v / n) : Eigen::Vector3d(1.0, 0.0, 0.0);
+    };
+
+    // Gram–Schmidt with original order: x from A, y from B orthogonalized vs x
+    Eigen::Vector3d e1 = safe_unit(A);
+    Eigen::Vector3d y = B - (e1.dot(B)) * e1;
+    Eigen::Vector3d e2 = safe_unit(y);
+
+    // z from cross to ensure right-handed basis
+    Eigen::Vector3d e3 = e1.cross(e2);
+    double n3 = e3.norm();
+    if (n3 < 1e-15) {
+        // Degenerate case: B nearly collinear with A → use C instead
+        y = C - (e1.dot(C)) * e1;
+        e2 = safe_unit(y);
+        e3 = e1.cross(e2);
+        n3 = e3.norm();
+        if (n3 < 1e-15) {
+            // Still degenerate: pick any perpendicular to e1
+            e2 = safe_unit((std::abs(e1.x()) < 0.9)
+                               ? Eigen::Vector3d::UnitX().cross(e1)
+                               : Eigen::Vector3d::UnitY().cross(e1));
+            e3 = e1.cross(e2);
+        }
+    } else {
+        e3 /= n3;
+    }
+
+    Eigen::Matrix3d R;
+    R.col(0) = e1;
+    R.col(1) = e2;
+    R.col(2) = e3;
+
+    // Convert rotation to Rodrigues (axis * angle)
+    Eigen::AngleAxisd aa(R);
+    Eigen::Vector3d r = aa.angle() * aa.axis();
+    rod[0] = r.x();
+    rod[1] = r.y();
+    rod[2] = r.z();
+}
+
+void LatticeToRodriguesAndLengths_Hex(const CrystalLattice &latt, double rod[3], double ac[3]) {
+    const Coord a = latt.Vec0();
+    const Coord b = latt.Vec1();
+    const Coord c = latt.Vec2();
+
+    Eigen::Vector3d A(a[0], a[1], a[2]);
+    Eigen::Vector3d B(b[0], b[1], b[2]);
+    Eigen::Vector3d C(c[0], c[1], c[2]);
+
+    const double a_len = A.norm();
+    const double b_len = B.norm();
+    const double c_len = C.norm();
+
+    ac[0] = (a_len + b_len) / 2.0;
+    ac[1] = (a_len + b_len) / 2.0;
+    ac[2] = c_len;
+
+    Eigen::Vector3d e1;
+    Eigen::Vector3d e3;
+
+    if (a_len > 0.0)
+        e1 = A / a_len;
+    else
+        e1 = Eigen::Vector3d::UnitX();
+
+    if (c_len > 0.0)
+        e3 = C / c_len;
+    else
+        e3 = Eigen::Vector3d::UnitZ();
+
+    Eigen::Vector3d e2 = e3.cross(e1);
+    if (e2.norm() < 1e-15) {
+        e2 = (std::abs(e1.x()) < 0.9)
+                 ? Eigen::Vector3d::UnitX().cross(e1)
+                 : Eigen::Vector3d::UnitY().cross(e1);
+    }
+    e2.normalize();
+    e3 = e1.cross(e2).normalized();
+
+    Eigen::Matrix3d R;
+    R.col(0) = e1;
+    R.col(1) = e2;
+    R.col(2) = e3;
+
+    Eigen::AngleAxisd aa(R);
+    Eigen::Vector3d r = aa.angle() * aa.axis();
+    rod[0] = r.x();
+    rod[1] = r.y();
+    rod[2] = r.z();
+}
+
+// Extract rotation (Rodrigues), lengths (a,b,c) and beta (rad) for monoclinic (unique axis b).
+// Frame choice: e2 aligned with b; e1 from a projected orthogonal to e2; e3 = e1 x e2.
+void LatticeToRodriguesLengthsBeta_Mono(const CrystalLattice &latt,
+                                        double rod[3],
+                                        double lengths[3],
+                                        double &beta_rad) {
+    const Coord a = latt.Vec0();
+    const Coord b = latt.Vec1();
+    const Coord c = latt.Vec2();
+
+    const Eigen::Vector3d A(a[0], a[1], a[2]);
+    const Eigen::Vector3d Bv(b[0], b[1], b[2]);
+    const Eigen::Vector3d C(c[0], c[1], c[2]);
+
+    // Unit cell lengths
+    const double a_len = A.norm();
+    const double b_len = Bv.norm();
+    const double c_len = C.norm();
+
+    lengths[0] = a_len;
+    lengths[1] = b_len;
+    lengths[2] = c_len;
+
+    // Monoclinic beta = angle(a, c)
+    double cos_beta = 0.0;
+    if (a_len > 1e-15 && c_len > 1e-15) {
+        cos_beta = A.dot(C) / (a_len * c_len);
+        cos_beta = std::clamp(cos_beta, -1.0, 1.0);
+    }
+
+    beta_rad = std::acos(cos_beta);
+
+    // Protect against singular construction
+    const double sin_beta = std::max(std::abs(std::sin(beta_rad)), 1e-12);
+
+    // Canonical monoclinic basis:
+    //
+    // B =
+    // [ 1   0   cos(beta) ]
+    // [ 0   1        0     ]
+    // [ 0   0   sin(beta) ]
+    //
+    Eigen::Matrix3d Bmono = Eigen::Matrix3d::Zero();
+    Bmono(0,0) = 1.0;
+    Bmono(1,1) = 1.0;
+    Bmono(0,2) = std::cos(beta_rad);
+    Bmono(2,2) = sin_beta;
+
+    // Scale by lengths
+    Eigen::DiagonalMatrix<double,3> D(a_len, b_len, c_len);
+
+    // Ideal body-frame lattice
+    const Eigen::Matrix3d M = Bmono * D;
+
+    // Observed lattice
+    Eigen::Matrix3d L;
+    L.col(0) = A;
+    L.col(1) = Bv;
+    L.col(2) = C;
+
+    // Estimate rotation:
+    // R ≈ L * M^{-1}
+    Eigen::Matrix3d R_est = L * M.inverse();
+
+    // Project to nearest proper rotation matrix
+    Eigen::JacobiSVD<Eigen::Matrix3d> svd(R_est, Eigen::ComputeFullU | Eigen::ComputeFullV);
+
+    Eigen::Matrix3d R = svd.matrixU() * svd.matrixV().transpose();
+
+    // Enforce det(R)=+1
+    if (R.determinant() < 0.0) {
+        Eigen::Matrix3d U = svd.matrixU();
+        U.col(2) *= -1.0;
+        R = U * svd.matrixV().transpose();
+    }
+
+    // Rodrigues vector
+    Eigen::AngleAxisd aa(R);
+    const Eigen::Vector3d r = aa.angle() * aa.axis();
+
+    rod[0] = r.x();
+    rod[1] = r.y();
+    rod[2] = r.z();
+}
+
+static inline Eigen::Matrix3d B_from_angles(double alpha_rad, double beta_rad, double gamma_rad) {
+    const double ca = std::cos(alpha_rad);
+    const double cb = std::cos(beta_rad);
+    const double cg = std::cos(gamma_rad);
+    const double sg = std::sin(gamma_rad);
+
+    Eigen::Matrix3d B = Eigen::Matrix3d::Identity();
+
+    // a along x, b in x-y, c general
+    B(0, 0) = 1.0;  B(1, 0) = 0.0;  B(2, 0) = 0.0;
+    B(0, 1) = cg;   B(1, 1) = sg;   B(2, 1) = 0.0;
+
+    // c vector components (standard crystallography construction)
+    const double cx = cb;
+    const double cy = (ca - cb * cg) / sg;
+    const double cz = std::sqrt(std::max(0.0, 1.0 - cx * cx - cy * cy));
+
+    B(0, 2) = cx;
+    B(1, 2) = cy;
+    B(2, 2) = cz;
+
+    return B;
+}
+
+CrystalLattice AngleAxisAndCellToLattice(const double rod[3],
+                                         const double lengths[3],
+                                         double alpha_rad,
+                                         double beta_rad,
+                                         double gamma_rad) {
+    const Eigen::Vector3d r(rod[0], rod[1], rod[2]);
+    const double angle = r.norm();
+
+    Eigen::Matrix3d R = Eigen::Matrix3d::Identity();
+    if (angle > 0.0)
+        R = Eigen::AngleAxisd(angle, r / angle).toRotationMatrix();
+
+    const Eigen::DiagonalMatrix<double, 3> D(lengths[0], lengths[1], lengths[2]);
+    const Eigen::Matrix3d B = B_from_angles(alpha_rad, beta_rad, gamma_rad);
+
+    // IMPORTANT convention: L = R * B * D (scale columns by lengths)
+    const Eigen::Matrix3d latt = R * B * D;
+
+    return CrystalLattice(Coord(latt(0, 0), latt(1, 0), latt(2, 0)),
+                          Coord(latt(0, 1), latt(1, 1), latt(2, 1)),
+                          Coord(latt(0, 2), latt(1, 2), latt(2, 2)));
+}
+
+
+CrystalLattice AngleAxisAndLengthsToLattice(const double rod[3], const double lengths[3], bool hex) {
+    return AngleAxisAndCellToLattice(rod, lengths,
+                                     /*alpha=*/M_PI / 2.0,
+                                     /*beta =*/M_PI / 2.0,
+                                     /*gamma=*/hex ? 2.0 * M_PI / 3.0 : M_PI / 2.0);
+}

image_analysis/geom_refinement/LatticeReduction.h

@@ -0,0 +1,19 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include "../common/CrystalLattice.h"
+
+void LatticeToRodriguesAndLengths_GS(const CrystalLattice &latt,  double rod[3], double lengths[3]);
+void LatticeToRodriguesAndLengths_Hex(const CrystalLattice &latt, double rod[3], double ac[3]);
+void LatticeToRodriguesLengthsBeta_Mono(const CrystalLattice &latt,
+                                               double rod[3],
+                                               double lengths[3],
+                                               double &beta_rad);
+
+CrystalLattice AngleAxisAndCellToLattice(const double rod[3],
+                                         const double lengths[3],
+                                         double alpha_rad,
+                                         double beta_rad,
+                                         double gamma_rad);

image_analysis/geom_refinement/XtalOptimizer.cpp

@@ -6,6 +6,7 @@
 #include "XtalOptimizer.h"
 #include "ceres/ceres.h"
 #include "ceres/rotation.h"
+#include "LatticeReduction.h"
 
 struct XtalResidual {
     XtalResidual(double x, double y,
@@ -244,265 +245,6 @@ struct RotationNormRegularizer {
     const double weight;
 };
 
-inline void LatticeToRodriguesAndLengths_GS(const CrystalLattice &latt,
-                                            double rod[3],
-                                            double lengths[3]) {
-    // Load lattice columns
-    const Coord a = latt.Vec0();
-    const Coord b = latt.Vec1();
-    const Coord c = latt.Vec2();
-
-    Eigen::Vector3d A(a[0], a[1], a[2]);
-    Eigen::Vector3d B(b[0], b[1], b[2]);
-    Eigen::Vector3d C(c[0], c[1], c[2]);
-
-    // Lengths = column norms (orthorhombic assumption)
-    lengths[0] = A.norm();
-    lengths[1] = B.norm();
-    lengths[2] = C.norm();
-
-    auto safe_unit = [](const Eigen::Vector3d &v, double eps = 1e-15) -> Eigen::Vector3d {
-        double n = v.norm();
-        return (n > eps) ? (v / n) : Eigen::Vector3d(1.0, 0.0, 0.0);
-    };
-
-    // Gram–Schmidt with original order: x from A, y from B orthogonalized vs x
-    Eigen::Vector3d e1 = safe_unit(A);
-    Eigen::Vector3d y = B - (e1.dot(B)) * e1;
-    Eigen::Vector3d e2 = safe_unit(y);
-
-    // z from cross to ensure right-handed basis
-    Eigen::Vector3d e3 = e1.cross(e2);
-    double n3 = e3.norm();
-    if (n3 < 1e-15) {
-        // Degenerate case: B nearly collinear with A → use C instead
-        y = C - (e1.dot(C)) * e1;
-        e2 = safe_unit(y);
-        e3 = e1.cross(e2);
-        n3 = e3.norm();
-        if (n3 < 1e-15) {
-            // Still degenerate: pick any perpendicular to e1
-            e2 = safe_unit((std::abs(e1.x()) < 0.9)
-                               ? Eigen::Vector3d::UnitX().cross(e1)
-                               : Eigen::Vector3d::UnitY().cross(e1));
-            e3 = e1.cross(e2);
-        }
-    } else {
-        e3 /= n3;
-    }
-
-    Eigen::Matrix3d R;
-    R.col(0) = e1;
-    R.col(1) = e2;
-    R.col(2) = e3;
-
-    // Convert rotation to Rodrigues (axis * angle)
-    Eigen::AngleAxisd aa(R);
-    Eigen::Vector3d r = aa.angle() * aa.axis();
-    rod[0] = r.x();
-    rod[1] = r.y();
-    rod[2] = r.z();
-}
-
-void LatticeToRodriguesAndLengths_Hex(const CrystalLattice &latt, double rod[3], double ac[3]) {
-    const Coord a = latt.Vec0();
-    const Coord b = latt.Vec1();
-    const Coord c = latt.Vec2();
-
-    Eigen::Vector3d A(a[0], a[1], a[2]);
-    Eigen::Vector3d B(b[0], b[1], b[2]);
-    Eigen::Vector3d C(c[0], c[1], c[2]);
-
-    const double a_len = A.norm();
-    const double b_len = B.norm();
-    const double c_len = C.norm();
-
-    ac[0] = (a_len + b_len) / 2.0;
-    ac[1] = (a_len + b_len) / 2.0;
-    ac[2] = c_len;
-
-    Eigen::Vector3d e1;
-    Eigen::Vector3d e3;
-
-    if (a_len > 0.0)
-        e1 = A / a_len;
-    else
-        e1 = Eigen::Vector3d::UnitX();
-
-    if (c_len > 0.0)
-        e3 = C / c_len;
-    else
-        e3 = Eigen::Vector3d::UnitZ();
-
-    Eigen::Vector3d e2 = e3.cross(e1);
-    if (e2.norm() < 1e-15) {
-        e2 = (std::abs(e1.x()) < 0.9)
-                 ? Eigen::Vector3d::UnitX().cross(e1)
-                 : Eigen::Vector3d::UnitY().cross(e1);
-    }
-    e2.normalize();
-    e3 = e1.cross(e2).normalized();
-
-    Eigen::Matrix3d R;
-    R.col(0) = e1;
-    R.col(1) = e2;
-    R.col(2) = e3;
-
-    Eigen::AngleAxisd aa(R);
-    Eigen::Vector3d r = aa.angle() * aa.axis();
-    rod[0] = r.x();
-    rod[1] = r.y();
-    rod[2] = r.z();
-}
-
-// Extract rotation (Rodrigues), lengths (a,b,c) and beta (rad) for monoclinic (unique axis b).
-// Frame choice: e2 aligned with b; e1 from a projected orthogonal to e2; e3 = e1 x e2.
-void LatticeToRodriguesLengthsBeta_Mono(const CrystalLattice &latt,
-                                        double rod[3],
-                                        double lengths[3],
-                                        double &beta_rad) {
-    const Coord a = latt.Vec0();
-    const Coord b = latt.Vec1();
-    const Coord c = latt.Vec2();
-
-    const Eigen::Vector3d A(a[0], a[1], a[2]);
-    const Eigen::Vector3d Bv(b[0], b[1], b[2]);
-    const Eigen::Vector3d C(c[0], c[1], c[2]);
-
-    // Unit cell lengths
-    const double a_len = A.norm();
-    const double b_len = Bv.norm();
-    const double c_len = C.norm();
-
-    lengths[0] = a_len;
-    lengths[1] = b_len;
-    lengths[2] = c_len;
-
-    // Monoclinic beta = angle(a, c)
-    double cos_beta = 0.0;
-    if (a_len > 1e-15 && c_len > 1e-15) {
-        cos_beta = A.dot(C) / (a_len * c_len);
-        cos_beta = std::clamp(cos_beta, -1.0, 1.0);
-    }
-
-    beta_rad = std::acos(cos_beta);
-
-    // Protect against singular construction
-    const double sin_beta = std::max(std::abs(std::sin(beta_rad)), 1e-12);
-
-    // Canonical monoclinic basis:
-    //
-    // B =
-    // [ 1   0   cos(beta) ]
-    // [ 0   1        0     ]
-    // [ 0   0   sin(beta) ]
-    //
-    Eigen::Matrix3d Bmono = Eigen::Matrix3d::Zero();
-    Bmono(0,0) = 1.0;
-    Bmono(1,1) = 1.0;
-    Bmono(0,2) = std::cos(beta_rad);
-    Bmono(2,2) = sin_beta;
-
-    // Scale by lengths
-    Eigen::DiagonalMatrix<double,3> D(a_len, b_len, c_len);
-
-    // Ideal body-frame lattice
-    const Eigen::Matrix3d M = Bmono * D;
-
-    // Observed lattice
-    Eigen::Matrix3d L;
-    L.col(0) = A;
-    L.col(1) = Bv;
-    L.col(2) = C;
-
-    // Estimate rotation:
-    // R ≈ L * M^{-1}
-    Eigen::Matrix3d R_est = L * M.inverse();
-
-    // Project to nearest proper rotation matrix
-    Eigen::JacobiSVD<Eigen::Matrix3d> svd(R_est, Eigen::ComputeFullU | Eigen::ComputeFullV);
-
-    Eigen::Matrix3d R = svd.matrixU() * svd.matrixV().transpose();
-
-    // Enforce det(R)=+1
-    if (R.determinant() < 0.0) {
-        Eigen::Matrix3d U = svd.matrixU();
-        U.col(2) *= -1.0;
-        R = U * svd.matrixV().transpose();
-    }
-
-    // Rodrigues vector
-    Eigen::AngleAxisd aa(R);
-    const Eigen::Vector3d r = aa.angle() * aa.axis();
-
-    rod[0] = r.x();
-    rod[1] = r.y();
-    rod[2] = r.z();
-}
-
-static inline Eigen::Matrix3d B_from_angles(double alpha_rad, double beta_rad, double gamma_rad) {
-    const double ca = std::cos(alpha_rad);
-    const double cb = std::cos(beta_rad);
-    const double cg = std::cos(gamma_rad);
-    const double sg = std::sin(gamma_rad);
-
-    Eigen::Matrix3d B = Eigen::Matrix3d::Identity();
-
-    // a along x, b in x-y, c general
-    B(0, 0) = 1.0;  B(1, 0) = 0.0;  B(2, 0) = 0.0;
-    B(0, 1) = cg;   B(1, 1) = sg;   B(2, 1) = 0.0;
-
-    // c vector components (standard crystallography construction)
-    const double cx = cb;
-    const double cy = (ca - cb * cg) / sg;
-    const double cz = std::sqrt(std::max(0.0, 1.0 - cx * cx - cy * cy));
-
-    B(0, 2) = cx;
-    B(1, 2) = cy;
-    B(2, 2) = cz;
-
-    return B;
-}
-
-CrystalLattice AngleAxisAndCellToLattice(const double rod[3],
-                                         const double lengths[3],
-                                         double alpha_rad,
-                                         double beta_rad,
-                                         double gamma_rad) {
-    const Eigen::Vector3d r(rod[0], rod[1], rod[2]);
-    const double angle = r.norm();
-
-    Eigen::Matrix3d R = Eigen::Matrix3d::Identity();
-    if (angle > 0.0)
-        R = Eigen::AngleAxisd(angle, r / angle).toRotationMatrix();
-
-    const Eigen::DiagonalMatrix<double, 3> D(lengths[0], lengths[1], lengths[2]);
-    const Eigen::Matrix3d B = B_from_angles(alpha_rad, beta_rad, gamma_rad);
-
-    // IMPORTANT convention: L = R * B * D (scale columns by lengths)
-    const Eigen::Matrix3d latt = R * B * D;
-
-    return CrystalLattice(Coord(latt(0, 0), latt(1, 0), latt(2, 0)),
-                          Coord(latt(0, 1), latt(1, 1), latt(2, 1)),
-                          Coord(latt(0, 2), latt(1, 2), latt(2, 2)));
-}
-
-
-CrystalLattice AngleAxisAndLengthsToLattice(const double rod[3], const double lengths[3], bool hex) {
-    if (!hex) {
-        return AngleAxisAndCellToLattice(rod, lengths,
-                                         /*alpha=*/M_PI / 2.0,
-                                         /*beta =*/M_PI / 2.0,
-                                         /*gamma=*/M_PI / 2.0);
-    }
-
-    // Hexagonal: caller must already enforce a=b in `lengths`.
-    return AngleAxisAndCellToLattice(rod, lengths,
-                                     /*alpha=*/M_PI / 2.0,
-                                     /*beta =*/M_PI / 2.0,
-                                     /*gamma=*/2.0 * M_PI / 3.0);
-}
-
 bool XtalOptimizerInternal(XtalOptimizerData &data,
                            const std::vector<std::vector<SpotToSave>> &spots,
                            const float tolerance) {

image_analysis/geom_refinement/XtalOptimizer.h

@@ -43,19 +43,6 @@ struct XtalOptimizerData {
     std::optional<Coord> angle_axis;
 };
 
-void LatticeToRodriguesAndLengths_GS(const CrystalLattice &latt,  double rod[3], double lengths[3]);
-void LatticeToRodriguesAndLengths_Hex(const CrystalLattice &latt, double rod[3], double ac[3]);
-void LatticeToRodriguesLengthsBeta_Mono(const CrystalLattice &latt,
-                                               double rod[3],
-                                               double lengths[3],
-                                               double &beta_rad);
-
-CrystalLattice AngleAxisAndCellToLattice(const double rod[3],
-                                         const double lengths[3],
-                                         double alpha_rad,
-                                         double beta_rad,
-                                         double gamma_rad);
-
 bool XtalOptimizer(XtalOptimizerData &data, const std::vector<std::vector<SpotToSave>> &spots);
 bool XtalOptimizerRotationOnly(XtalOptimizerData &data, const std::vector<SpotToSave> &spots, float tolerance);
 

image_analysis/indexing/FitProfileRadius.cpp

@@ -1,9 +1,6 @@
 // SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
 // SPDX-License-Identifier: GPL-3.0-only
 
-// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
-// SPDX-License-Identifier: GPL-3.0-only
-
 #include "FitProfileRadius.h"
 #include <algorithm>  // std::nth_element
 #include <cmath>      // std::fabs

image_analysis/pixel_refinement/CMakeLists.txt

@@ -0,0 +1,2 @@
+ADD_LIBRARY(JFJochPixelRefine PixelRefine.cpp PixelRefine.h)
+TARGET_LINK_LIBRARIES(JFJochPixelRefine JFJochBraggPrediction JFJochCommon JFJochScaleMerge JFJochGeomRefinement ceres_static)

image_analysis/pixel_refinement/PixelRefine.h

@@ -0,0 +1,197 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include "../bragg_prediction/BraggPrediction.h"
+#include "../common/DiffractionExperiment.h"
+#include "../common/AzimuthalIntegrationMapping.h"
+#include "../common/AzimuthalIntegrationProfile.h"
+#include "../scale_merge/HKLKey.h"
+
+// =============================================================================
+// PixelRefine — one optimization to rule geometry, integration and scaling
+// =============================================================================
+//
+// Intent
+// ------
+// Classical crystallographic data processing is a one-way pipeline:
+//
+//     spot finding -> indexing -> geometry refinement -> integration -> scaling -> merging
+//
+// Each stage consumes the previous stage's output and never talks back. The
+// integrator trusts the refined geometry; the scaler trusts the integrated
+// intensities; nothing downstream is ever allowed to correct an upstream
+// parameter. Post-refinement and profile fitting were the field's partial
+// answers to this: post-refinement lets merged intensities nudge per-image
+// orientation/cell/mosaicity, and profile fitting lets a learned spot shape
+// improve weak-reflection intensities. But both are narrow back-channels bolted
+// onto a feed-forward pipeline — there is no end-to-end gradient that flows from
+// the raw detector pixels all the way back to every parameter at once.
+//
+// PixelRefine is an experiment in doing the whole thing as a *single* least
+// 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
+//
+// and let Ceres autodiff back-propagate the per-pixel residuals into ALL of:
+//   * detector geometry  (beam centre, distance, tilt)
+//   * crystal orientation + unit cell
+//   * overall scale G and Debye-Waller B
+//   * the reciprocal-space spot widths R = (radial, tangential)
+// simultaneously. Geometry refinement, profile-fitted integration and scaling
+// then stop being separate stages: they are different parameters of one model,
+// coupled through the same pixels, with full backpropagation between them. Once
+// the model is differentiable end-to-end, things that used to need bespoke code
+// — mosaicity refinement, profile fitting, partiality — fall out "for free" as
+// extra parameters of the same forward model.
+//
+// How I_true enters
+// -----------------
+// I_true is NOT refined here. It is a *fixed hypothesis* for the duration of a
+// pass: the current best merged estimate of each reflection's full intensity.
+// The intended outer loop is iterative, like EM / self-consistent field:
+//
+//     repeat over the whole dataset:
+//         run PixelRefine on every image with the current I_true reference
+//         re-merge the resulting intensities -> new I_true
+//     until the reference stops changing
+//
+// So a single PixelRefine call answers "given this intensity hypothesis, what
+// geometry/scale/profile best explains these pixels, and what intensities do I
+// read back out?", and the dataset-level loop refines the hypothesis itself.
+//
+// Inner predict<->refine loop
+// ---------------------------
+// Within one image we also iterate (max_iterations): Bragg prediction places the
+// shoeboxes, we refine, the refined geometry/cell feed the next prediction, etc.
+// Initially the model geometry equals the experiment's DiffractionGeometry, but
+// as refinement proceeds it diverges, so later predictions must use the *refined*
+// data.geom rather than the static experiment geometry.
+//
+// On shoeboxes, tails, and gatekeeping
+// ------------------------------------
+// We deliberately do NOT chase the full spot. A shoebox only needs to cover
+// enough of a reflection for the fit to be meaningful; clipped tails are fine,
+// because the partiality term already downweights whatever falls outside the
+// well-modelled core. The premise - especially for serial crystallography - is
+// that the problem is rarely *missing* information; it is *failing to gate out*
+// information that is not meaningful. As long as the model knows a piece is
+// missing (low partiality), it is safe to leave it missing. That flips the usual
+// trade-off: rather than shrinking boxes to avoid contamination, we can grow them
+// and let partiality decide, per pixel and per reflection, what actually carries
+// signal. (For downstream integration this pixel-level gating is the point - keep
+// only meaningful pixels, instead of a fixed geometric mask.) The bandwidth term
+// below is part of the same idea: it tells the model where the radial tails are
+// *expected* to be, so it can weight rather than blindly include them.
+//
+// X-ray bandwidth (optional)
+// --------------------------
+// A finite bandwidth thickens the Ewald shell radially and smears spots along the
+// radial direction, growing like 1/d^2 (the pink-beam/DMM signature). It enters
+// as a fixed, resolution-dependent addition to the radial width R0 (see
+// PixelRefineData::bandwidth and PixelRefine.cpp). It is OFF by default
+// (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.
+// =============================================================================
+
+struct PixelRefineData {
+    // --- model state (input as initial guess, output as refined result) ---
+    DiffractionGeometry geom;
+    CrystalLattice latt;
+    gemmi::CrystalSystem crystal_system = gemmi::CrystalSystem::Triclinic;
+    char centering = 'P';
+
+    double B_factor = 0.0;      // Debye-Waller B (A^2)
+    double scale_factor = 1.0;  // overall scale G
+    double R[2] = {0.005, 0.005}; // R[0] = radial (partiality) width, R[1] = tangential (profile) width (A^-1)
+
+    // Relative X-ray bandwidth (sigma of dlambda/lambda), e.g. ~0.004 for a 1%
+    // FWHM DMM, ~1e-4 for Si(111). Adds a resolution-dependent radial broadening
+    // 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    max_iterations = 3; // inner predict<->refine cycles (re-predict with refined geom/latt)
+
+    // --- output ---
+    std::vector<Reflection> reflections; // profile-fitted integration result
+    bool   solved = false;
+    double final_cost = NAN;
+    size_t residual_count = 0;
+    double cc = NAN;        // per-image CC of scaled intensities vs reference
+    int64_t cc_n = 0;       // number of reflections in the CC
+};
+
+class PixelRefine {
+    const AzimuthalIntegrationMapping &mapping;
+    const size_t xpixel, ypixel;
+    const DiffractionExperiment &experiment;
+
+    const HKLKeyGenerator hkl_key_generator;
+    std::map<HKLKey, double> reference_data;
+
+    // Fills the Ceres parameter blocks (geometry + symmetry-aware lattice
+    // parametrization) from the current model state. Shared by Run and
+    // 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 latt_vec0[3], double latt_vec1[3], double latt_vec2[3]) const;
+public:
+    PixelRefine(const DiffractionExperiment &experiment,
+                const AzimuthalIntegrationMapping &mapping,
+                const std::vector<MergedReflection> &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.
+    template<class T>
+    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
+    // 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<float> PredictImage(const AzimuthalIntegrationProfile &profile,
+                                    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.
+    template<class T>
+    std::vector<float> ChiSquaredImage(const T *image,
+                                       const AzimuthalIntegrationProfile &profile,
+                                       BraggPrediction &prediction,
+                                       const PixelRefineData &data) const;
+};

image_analysis/scale_merge/Merge.cpp

@@ -22,7 +22,14 @@ MergeOnTheFly::MergeOnTheFly(const DiffractionExperiment &x)
       scaling_settings(x.GetScalingSettings()),
       indexing_settings(x.GetIndexingSettings()),
       high_resolution_limit(scaling_settings.GetHighResolutionLimit_A()),
-      image_cc_limit(scaling_settings.GetMinCCForImage()),
+      // A min-image-CC of 0 (the default) means "no limit": leave the optional
+      // empty so the per-image CC cut is inactive. Otherwise a 0.0 threshold
+      // would silently drop every image with a non-positive per-image CC (which
+      // also wrongly zeroed N_obs in MergeStats, since it masks with cc_mask=true
+      // while the merge keeps all images).
+      image_cc_limit(scaling_settings.GetMinCCForImage() > 0.0
+                         ? std::optional<double>(scaling_settings.GetMinCCForImage())
+                         : std::nullopt),
       min_partiality(scaling_settings.GetMinPartiality()),
       generator(scaling_settings.GetMergeFriedel(), space_group_number) {
 }
@@ -406,7 +413,7 @@ std::ostream &operator<<(std::ostream &output, const MergeStatistics &in) {
     output << std::endl;
     output << fmt::format("  {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s}",
         "d_min", "N_obs", "N_uniq", "N_possib", "Compl","<I/sig>", "CC1/2", "CCref")
-    << std::endl;;
+    << std::endl;
     output << fmt::format("  {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s}",
         "", "", "", "", "", "", "", "") << std::endl;
     for (const auto &sh: in.shells) {

image_analysis/scale_merge/ScaleOnTheFly.cpp

@@ -105,145 +105,6 @@ namespace {
         double partiality;
     };
 
-struct ScalingPostRefResidual : public ScalingResidual {
-    // Postrefinement for still images
-    //
-    // In this algorithm at the point of post-refinement we don't anymore care for where maximum of
-    // the reflection was located and if it fits the observed position.
-    // This reflections was already integrated and we cannot integrate it better at this point.
-    // But, we could adjust partiality to indicate that this reflection was wrongly predicted.
-    // I.e., integrated position was far away from true reflection, so partiality must be low.
-    // This is an empiric model and need to see if this will work in practice at all.
-    // I hope it will allow the model to find that reflections were misindexed.
-    // We assume at this point that initial indexing was done properly and integration was generally OK => most low resolution reflections fit correctly
-    // Yet we know, that small errors in indexing are inducing misalignment at high resolution - sometimes it is visible that high-resolution reflections
-    // are away from shoe-boxes observed in the image, if we can catch this at post-refinement/scaling step, this would be great.
-    // Next logical step is to do this pixel-wise - for each pixel refine partiality and merge pixels
-    // This should work for per-image scaling, or even, maybe, for full rotation datasets (3600 images)
-    // Then we could properly take into account misalignment of shoe-box center vs. partiality and also remove most pixels
-    // in the shoe-box that don't really contribute to the reflection.
-    // But...it could also drift to downweighting partiality for all high resolution reflections to make loss function "fake happy".
-
-    // We assume rot3 == 0. Rot3 is not really helping much in crystallography (other than fixing polarization correction)
-    ScalingPostRefResidual(const Reflection &r, double Itrue, double sigma,
-                           const DiffractionGeometry &geometry,
-                           const CrystalLattice &lattice)
-    : ScalingResidual(r, Itrue, sigma),
-          integration_center_x(r.predicted_x),
-          integration_center_y(r.predicted_y),
-          inv_lambda(SafeInv(geometry.GetWavelength_A(), 0.0)),
-          pixel_size(geometry.GetPixelSize_mm()),
-          det_dist_mm(geometry.GetDetectorDistance_mm()),
-          beam_x(geometry.GetBeamX_pxl()),
-          beam_y(geometry.GetBeamY_pxl()),
-          exp_h(r.h),
-          exp_k(r.k),
-          exp_l(r.l),
-          Astar(lattice.Astar()),
-          Bstar(lattice.Bstar()),
-          Cstar(lattice.Cstar()),
-          c1(std::cos(geometry.GetPoniRot1_rad())),
-          s1(std::sin(geometry.GetPoniRot1_rad())),
-          c2(std::cos(geometry.GetPoniRot2_rad())),
-          s2(std::sin(geometry.GetPoniRot2_rad())) {
-    }
-
-    template <typename T>
-    T CalcPartiality(const T *const R,
-                    const T *const beam_corr,
-                    const T *const p0) const {
-        // Detector coordinates in mm
-        const T det_x = (T(integration_center_x) - beam_x - beam_corr[0]) * T(pixel_size);
-        const T det_y = (T(integration_center_y) - beam_y - beam_corr[1]) * T(pixel_size);
-        const T det_z = T(det_dist_mm);
-
-        // Apply Ry(rot1) first: rotate around Y
-        const T t1_x = T(c1) * det_x + T(s1) * det_z;
-        const T t1_y = det_y;
-        const T t1_z = T(-s1) * det_x + T(c1) * det_z;
-
-        // Then apply Rx(-rot2): rotate around X
-        const T x = t1_x;
-        const T y = T(c2) * t1_y + T(s2) * t1_z;
-        const T z = - T(s2) * t1_y + T(c2) * t1_z;
-
-        // convert to recip space
-        const T lab_norm = ceres::sqrt(x * x + y * y + z * z);
-        const T inv_norm = T(1) / lab_norm;
-
-        T recip_obs[3];
-        recip_obs[0] = x * inv_norm * inv_lambda;
-        recip_obs[1] = y * inv_norm * inv_lambda;
-        recip_obs[2] = (z * inv_norm - T(1.0)) * inv_lambda;
-
-        const Eigen::Matrix<T, 3, 1> e_obs_recip(recip_obs[0], recip_obs[1], recip_obs[2]);
-
-        const T astar_unrot[3] = {T(Astar.x), T(Astar.y), T(Astar.z)};
-        const T bstar_unrot[3] = {T(Bstar.x), T(Bstar.y), T(Bstar.z)};
-        const T cstar_unrot[3] = {T(Cstar.x), T(Cstar.y), T(Cstar.z)};
-
-        T astar_rot[3], bstar_rot[3], cstar_rot[3];
-
-        ceres::AngleAxisRotatePoint(p0, astar_unrot, astar_rot);
-        ceres::AngleAxisRotatePoint(p0, bstar_unrot, bstar_rot);
-        ceres::AngleAxisRotatePoint(p0, cstar_unrot, cstar_rot);
-
-        const Eigen::Matrix<T, 3, 1> e_pred_recip(T(exp_h) * astar_rot[0] + T(exp_k) * bstar_rot[0] + T(exp_l) * cstar_rot[0],
-            T(exp_h) * astar_rot[1] + T(exp_k) * bstar_rot[1] + T(exp_l) * cstar_rot[1],
-            T(exp_h) * astar_rot[2] + T(exp_k) * bstar_rot[2] + T(exp_l) * cstar_rot[2]
-        );
-
-        // Ewald sphere centre is at -k_i = (0, 0, -inv_lambda)
-        // Radial direction: outward normal at g_hkl
-        const Eigen::Matrix<T, 3, 1> S_pred(
-            e_pred_recip[0],
-            e_pred_recip[1],
-            e_pred_recip[2] + T(inv_lambda)   // g_hkl + k_i
-        );
-        const T S_pred_norm = S_pred.norm();
-        if (S_pred_norm < T(1e-10))
-            return T(0);
-
-        const Eigen::Matrix<T, 3, 1> n_radial = S_pred / S_pred_norm;
-
-        const Eigen::Matrix<T, 3, 1> delta_q = e_obs_recip - e_pred_recip;
-        const T eps_radial        = delta_q.dot(n_radial);
-        const Eigen::Matrix<T, 3, 1> dq_tang = delta_q - eps_radial * n_radial;
-        const T eps_tangential_sq = dq_tang.squaredNorm();   // guaranteed ≥ 0
-        // ─────────────────────────────────────────────────────────────
-
-        return ceres::exp(- eps_radial * eps_radial / (R[0] * R[0]) - eps_tangential_sq / (R[1] * R[1]));
-    }
-
-    template<typename T>
-    bool operator()(const T *const G,
-                    const T *const B,
-                    const T *const R,
-                    const T *const beam_corr,
-                    const T *const p0,
-                    T *residual) const {
-        if (R[0] < T(1e-10) || R[1] < T(1e-10))
-            return false;
-
-        const T B_term = ceres::exp(B[0] * T(b_resolution_coeff));
-
-        const T partiality = CalcPartiality(R, beam_corr, p0);
-        residual[0] = (G[0] * partiality * B_term * T(lp) * T(Itrue)
-                       - T(Iobs)) * T(weight);
-        return true;
-    }
-
-    const double integration_center_x, integration_center_y;
-    const double inv_lambda;
-    const double pixel_size;
-    const double det_dist_mm;
-    const double beam_x, beam_y;
-    const double exp_h;
-    const double exp_k;
-    const double exp_l;
-    const Coord Astar, Bstar, Cstar;
-    const double c1,s1,c2,s2;
-};
     struct RotationNormRegularizer {
         explicit RotationNormRegularizer(double weight) : weight(weight) {}
 
@@ -282,7 +143,6 @@ bool ScaleOnTheFly::Accept(const Reflection &r) const {
             return std::isfinite(r.zeta) && r.zeta > 0.0f;
         case PartialityModel::Fixed:
         case PartialityModel::Unity:
-        case PartialityModel::Postrefinement:
             return true;
     }
 
@@ -395,12 +255,6 @@ void ScaleOnTheFly::Scale(IntegrationOutcome &integration_outcome) const {
                 problem.AddResidualBlock(cost, nullptr, &result.G, &result.B);
             }
             break;
-            case PartialityModel::Postrefinement: {
-                auto *cost = new ceres::AutoDiffCostFunction<ScalingPostRefResidual, 1, 1, 1, 2, 2, 3>(
-                    new ScalingPostRefResidual(r, Itrue, sigma, integration_outcome.geom, integration_outcome.latt));
-                problem.AddResidualBlock(cost, nullptr, &result.G, &result.B, result.R, result.beam_corr, result.p0);
-            }
-            break;
             case PartialityModel::Rotation: {
                 auto *cost = new ceres::AutoDiffCostFunction<ScalingRotationResidual, 1, 1, 1, 1, 1>(
                     new ScalingRotationResidual(r, Itrue, sigma));
@@ -460,21 +314,6 @@ void ScaleOnTheFly::Scale(IntegrationOutcome &integration_outcome) const {
         problem.SetParameterUpperBound(&result.mos, 0, s.GetMaxMosaicity());
     }
 
-    if (model == PartialityModel::Postrefinement) {
-        problem.SetParameterLowerBound(result.R, 0, 1e-6);
-        problem.SetParameterLowerBound(result.R, 1, 1e-6);
-
-        // Beam center can be off by max 1 pixel
-        problem.SetParameterLowerBound(result.beam_corr, 0, -1);
-        problem.SetParameterUpperBound(result.beam_corr, 0, 1);
-        problem.SetParameterLowerBound(result.beam_corr, 1, -1);
-        problem.SetParameterUpperBound(result.beam_corr, 1, 1);
-
-        auto *angle_reg_cost = new ceres::AutoDiffCostFunction<RotationNormRegularizer, 3, 3>(
-            new RotationNormRegularizer(0.05));
-        problem.AddResidualBlock(angle_reg_cost, nullptr, result.p0);
-    }
-
     ceres::Solver::Options options;
     options.linear_solver_type = ceres::DENSE_QR;
     options.minimizer_progress_to_stdout = false;
@@ -495,19 +334,12 @@ void ScaleOnTheFly::Scale(IntegrationOutcome &integration_outcome) const {
                     r.partiality = RotationPartiality(r.delta_phi_deg, r.zeta, result.mos, result.wedge);
                 break;
             }
-            case PartialityModel::Postrefinement: {
-                ScalingPostRefResidual residual(r, 0, 0, integration_outcome.geom, integration_outcome.latt);
-                r.partiality = static_cast<float>(residual.CalcPartiality(result.R, result.beam_corr, result.p0));
-            }
-                break;
             default:
                 // For fixed partiality there is no need to change anything
                 break;
         }
         const double denom = B_term * r.partiality * result.G;
-        if (std::isfinite(r.rlp) &&
-            std::isfinite(denom) &&
-            denom > 0.0) {
+        if (std::isfinite(r.rlp) && std::isfinite(denom) && denom > 0.0) {
                 r.image_scale_corr = static_cast<float>(r.rlp / denom);
             } else {
                 r.image_scale_corr = NAN;

image_analysis/scale_merge/ScalingResult.cpp

@@ -36,8 +36,8 @@ void ScalingResult::SaveToFile(const std::string &filename) {
     const std::string img_path = filename + "_image.dat";
     std::ofstream img_file(img_path,  std::ofstream::out | std::ofstream::trunc);
     if (!img_file) {
-        throw JFJochException(JFJochExceptionCategory::FileWriteError
-                              , "Cannot open {} for writing");
+        throw JFJochException(JFJochExceptionCategory::FileWriteError,
+                              "Cannot open " + img_path + " for writing");
     }
 
     for (size_t i = 0; i < image_scale_g.size(); ++i) {

reader/JFJochHDF5Reader.cpp

@@ -498,8 +498,14 @@ void JFJochHDF5Reader::ReadFile(const std::string &filename) {
             det_distance = 0.1; // Set to 100 mm, if det distance is less than 1 mm
         dataset->experiment.DetectorDistance_mm(det_distance * 1000.0);
 
-        dataset->experiment.IncidentEnergy_keV(
-            WVL_1A_IN_KEV / master_file->GetFloat("/entry/instrument/beam/incident_wavelength"));
+        const float incident_wavelength_A = master_file->GetFloat("/entry/instrument/beam/incident_wavelength");
+        dataset->experiment.IncidentEnergy_keV(WVL_1A_IN_KEV / incident_wavelength_A);
+
+        // NXmx incident_wavelength_spread is the absolute FWHM (Angstrom); store it
+        // as the relative bandwidth FWHM (dlambda/lambda) used internally.
+        if (const auto spread = master_file->GetOptFloat("/entry/instrument/beam/incident_wavelength_spread"))
+            if (incident_wavelength_A > 0.0f)
+                dataset->experiment.BandwidthFWHM(spread.value() / incident_wavelength_A);
 
         dataset->error_value = master_file->GetOptInt("/entry/instrument/detector/error_value");
 

tests/CMakeLists.txt

@@ -71,6 +71,7 @@ ADD_EXECUTABLE(jfjoch_test
         UnitCellTest.cpp
         CCTest.cpp
         MultiLatticeSearchTest.cpp
+        LatticeReductionTest.cpp
 )
 
 target_link_libraries(jfjoch_test Catch2WithMain JFJochBroker JFJochReceiver JFJochReader JFJochWriter

tests/LatticeReductionTest.cpp

@@ -0,0 +1,225 @@
+// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <catch2/catch_all.hpp>
+
+#include <Eigen/Dense>
+
+#include "../image_analysis/geom_refinement/LatticeReduction.h"
+
+
+namespace {
+    Eigen::Vector3d to_eigen(const Coord& v) {
+        return {v[0], v[1], v[2]};
+    }
+
+    double angle_rad(const Eigen::Vector3d& a, const Eigen::Vector3d& b) {
+        const double na = a.norm();
+        const double nb = b.norm();
+        if (na == 0.0 || nb == 0.0)
+            return 0.0;
+        double c = a.dot(b) / (na * nb);
+        c = std::max(-1.0, std::min(1.0, c));
+        return std::acos(c);
+    }
+
+    // Compare two lattices up to a global rotation: compare Gram matrices G = L^T L (rotation-invariant).
+    Eigen::Matrix3d gram(const CrystalLattice& latt) {
+        const Eigen::Vector3d A = to_eigen(latt.Vec0());
+        const Eigen::Vector3d B = to_eigen(latt.Vec1());
+        const Eigen::Vector3d C = to_eigen(latt.Vec2());
+        Eigen::Matrix3d G;
+        G(0,0) = A.dot(A); G(0,1) = A.dot(B); G(0,2) = A.dot(C);
+        G(1,0) = B.dot(A); G(1,1) = B.dot(B); G(1,2) = B.dot(C);
+        G(2,0) = C.dot(A); G(2,1) = C.dot(B); G(2,2) = C.dot(C);
+        return G;
+    }
+
+    void check_gram_close(const CrystalLattice& a,
+                          const CrystalLattice& b,
+                          double abs_eps,
+                          double rel_eps) {
+        const Eigen::Matrix3d Ga = gram(a);
+        const Eigen::Matrix3d Gb = gram(b);
+
+        for (int r = 0; r < 3; ++r) {
+            for (int c = 0; c < 3; ++c) {
+                const double va = Ga(r, c);
+                const double vb = Gb(r, c);
+
+                // Scale for relative error; avoid blowing up around zero.
+                const double scale = std::max({1.0, std::abs(va), std::abs(vb)});
+                const double tol = std::max(abs_eps, rel_eps * scale);
+
+                INFO("G(" << r << "," << c << ") va=" << va << " vb=" << vb
+                          << " scale=" << scale << " tol=" << tol);
+
+                CHECK(va == Catch::Approx(vb).margin(tol));
+            }
+        }
+    }
+} // namespace
+
+
+TEST_CASE("LatticeToRodrigues") {
+    double rod[3];
+    double lengths[3];
+
+    CrystalLattice latt_i(40,50,80,90,90,90);
+
+    LatticeToRodriguesAndLengths_GS(latt_i, rod, lengths);
+    CHECK(lengths[0] == Catch::Approx(40.0));
+    CHECK(lengths[1] == Catch::Approx(50.0));
+    CHECK(lengths[2] == Catch::Approx(80.0));
+    CHECK(fabs(rod[0]) < 0.001);
+    CHECK(fabs(rod[1]) < 0.001);
+    CHECK(fabs(rod[2]) < 0.001);
+
+    auto latt_o = AngleAxisAndCellToLattice(rod, lengths, M_PI/2, M_PI/2, M_PI/2);
+
+    CHECK(latt_o.Vec0().Length() == Catch::Approx(40.0));
+    CHECK(latt_o.Vec1().Length() == Catch::Approx(50.0));
+    CHECK(latt_o.Vec2().Length() == Catch::Approx(80.0));
+}
+
+
+TEST_CASE("LatticeToRodrigues_irregular") {
+    double rod[3];
+    double lengths[3];
+
+    CrystalLattice latt_i(Coord(40,0,0),
+                          Coord(0, 50 / sqrt(2), -50 / sqrt(2)),
+                          Coord(0, 80 / sqrt(2), 80 / sqrt(2)));
+
+    LatticeToRodriguesAndLengths_GS(latt_i, rod, lengths);
+    CHECK(lengths[0] == Catch::Approx(40.0));
+    CHECK(lengths[1] == Catch::Approx(50.0));
+    CHECK(lengths[2] == Catch::Approx(80.0));
+
+    auto latt_o = AngleAxisAndCellToLattice(rod, lengths, M_PI/2, M_PI/2, M_PI/2);
+
+    CHECK(latt_o.Vec0().Length() == Catch::Approx(40.0));
+    CHECK(latt_o.Vec1().Length() == Catch::Approx(50.0));
+    CHECK(latt_o.Vec2().Length() == Catch::Approx(80.0));
+}
+
+TEST_CASE("LatticeToRodrigues_Hex") {
+    double rod[3];
+    double lengths[3];
+
+    Coord a = Coord(40,0,0);
+    Coord b = Coord(40 / 2, 40 * sqrt(3)/ 2.0, 0);
+    Coord c = Coord(0, 0, 70);
+
+    RotMatrix R(1.0, Coord(0,1,1));
+    CrystalLattice latt_i(R*a,R*b,R*c);
+
+    LatticeToRodriguesAndLengths_Hex(latt_i, rod, lengths);
+    CHECK(lengths[0] == Catch::Approx(40.0));
+    CHECK(lengths[1] == Catch::Approx(40.0));
+    CHECK(lengths[2] == Catch::Approx(70.0));
+
+    auto latt_o = AngleAxisAndCellToLattice(rod, lengths,M_PI / 2.0, M_PI / 2.0, 2.0 * M_PI / 3.0);
+    auto uc_o = latt_o.GetUnitCell();
+    CHECK(uc_o.a == Catch::Approx(40.0));
+    CHECK(uc_o.b == Catch::Approx(40.0));
+    CHECK(uc_o.c == Catch::Approx(70.0));
+    CHECK(uc_o.alpha == Catch::Approx(90.0));
+    CHECK(uc_o.beta == Catch::Approx(90.0));
+    CHECK(uc_o.gamma == Catch::Approx(120.0));
+}
+
+
+TEST_CASE("LatticeReduction Lattice param roundtrip (GS) preserves Gram matrix") {
+    // Non-orthogonal, irregular basis, but still a valid lattice
+    CrystalLattice latt_i(Coord(40,  1,  2),
+                          Coord( 3, 50, -4),
+                          Coord(-5,  6, 80));
+
+    double rod[3]{}, lengths[3]{};
+    LatticeToRodriguesAndLengths_GS(latt_i, rod, lengths);
+    auto latt_o = AngleAxisAndCellToLattice(rod, lengths, M_PI/2, M_PI/2, M_PI/2);
+
+    // This parametrization only keeps "lengths + rotation", i.e. it cannot reproduce shear.
+    // So we *do not* compare Gram matrices here.
+    // Instead, sanity-check: reconstructed vectors have the requested lengths.
+    CHECK(latt_o.Vec0().Length() == Catch::Approx(lengths[0]).margin(1e-9));
+    CHECK(latt_o.Vec1().Length() == Catch::Approx(lengths[1]).margin(1e-9));
+    CHECK(latt_o.Vec2().Length() == Catch::Approx(lengths[2]).margin(1e-9));
+}
+
+TEST_CASE("LatticeReduction Lattice param roundtrip (Hex) preserves unit cell") {
+    Coord a = Coord(40, 0, 0);
+    Coord b = Coord(40 / 2.0, 40 * std::sqrt(3) / 2.0, 0);
+    Coord c = Coord(0, 0, 70);
+
+    // Apply an arbitrary rotation to ensure the rod extraction is meaningful
+    RotMatrix R(1.0, Coord(0, 1, 1));
+    CrystalLattice latt_i(R * a, R * b, R * c);
+
+    double rod[3]{}, ac[3]{};
+    LatticeToRodriguesAndLengths_Hex(latt_i, rod, ac);
+    auto latt_o = AngleAxisAndCellToLattice(rod, ac, M_PI/2, M_PI/2, 2*M_PI/3);
+
+    auto uc_o = latt_o.GetUnitCell();
+    CHECK(uc_o.a == Catch::Approx(40.0).margin(1e-6));
+    CHECK(uc_o.b == Catch::Approx(40.0).margin(1e-6));
+    CHECK(uc_o.c == Catch::Approx(70.0).margin(1e-6));
+    CHECK(uc_o.alpha == Catch::Approx(90.0).margin(1e-6));
+    CHECK(uc_o.beta  == Catch::Approx(90.0).margin(1e-6));
+    CHECK(uc_o.gamma == Catch::Approx(120.0).margin(1e-6));
+}
+
+TEST_CASE("LatticeReduction Monoclinic param roundtrip") {
+    struct Case { double beta_deg; };
+    const std::vector<Case> cases = {
+        {60.0}, {75.0}, {115.0}, {130.0}
+    };
+
+    for (const auto& cs : cases) {
+        INFO("beta_deg=" << cs.beta_deg);
+
+        // Start from a clean monoclinic cell in its conventional setting (unique axis b).
+        CrystalLattice latt0(50, 60, 70, 90, cs.beta_deg, 90);
+
+        // Now apply a TRUE global rotation: rotate each basis vector (left-multiply).
+        RotMatrix R(0.7, Coord(0.3, 0.9, 0.1));
+        CrystalLattice latt_i(R * latt0.Vec0(),
+                              R * latt0.Vec1(),
+                              R * latt0.Vec2());
+
+        double rod[3]{}, lengths[3]{}, beta_rad = 0.0;
+        LatticeToRodriguesLengthsBeta_Mono(latt_i, rod, lengths, beta_rad);
+
+        // Basic sanity
+        CHECK(lengths[0] == Catch::Approx(50.0).margin(1e-6));
+        CHECK(lengths[1] == Catch::Approx(60.0).margin(1e-6));
+        CHECK(lengths[2] == Catch::Approx(70.0).margin(1e-6));
+        CHECK(beta_rad * 180.0 / M_PI == Catch::Approx(cs.beta_deg).margin(1e-6));
+
+        auto latt_o = AngleAxisAndCellToLattice(rod, lengths, M_PI/2, beta_rad, M_PI/2);
+
+        // Rotation-invariant check: Gram matrices match.
+        check_gram_close(latt_i, latt_o, /*abs_eps=*/5e-4, /*rel_eps=*/1e-10);
+
+        // Also check the unit-cell angles we expect for monoclinic(unique b).
+        auto uc_o = latt_o.GetUnitCell();
+        CHECK(uc_o.alpha == Catch::Approx(90.0).margin(1e-4));
+        CHECK(uc_o.gamma == Catch::Approx(90.0).margin(1e-4));
+        CHECK(uc_o.beta  == Catch::Approx(cs.beta_deg).margin(1e-4));
+    }
+}
+
+TEST_CASE("LatticeReduction monoclinic") {
+    // This isolates only the beta definition, independent of other choices.
+    CrystalLattice latt_i(50, 60, 70, 90, 130, 90);
+
+    const Eigen::Vector3d A = to_eigen(latt_i.Vec0());
+    const Eigen::Vector3d C = to_eigen(latt_i.Vec2());
+    const double beta_geom = angle_rad(A, C);
+
+    double rod[3]{}, lengths[3]{}, beta_rad = 0.0;
+    LatticeToRodriguesLengthsBeta_Mono(latt_i, rod, lengths, beta_rad);
+
+    CHECK(beta_rad == Catch::Approx(beta_geom).margin(1e-12));
+}

tests/XtalOptimizerTest.cpp

@@ -489,73 +489,6 @@ TEST_CASE("XtalOptimizer_monoclinic") {
     CHECK(fabs(uc_o.gamma - 90) < 0.05);
 }
 
-TEST_CASE("LatticeToRodrigues") {
-    double rod[3];
-    double lengths[3];
-
-    CrystalLattice latt_i(40,50,80,90,90,90);
-
-    LatticeToRodriguesAndLengths_GS(latt_i, rod, lengths);
-    CHECK(lengths[0] == Catch::Approx(40.0));
-    CHECK(lengths[1] == Catch::Approx(50.0));
-    CHECK(lengths[2] == Catch::Approx(80.0));
-    CHECK(fabs(rod[0]) < 0.001);
-    CHECK(fabs(rod[1]) < 0.001);
-    CHECK(fabs(rod[2]) < 0.001);
-
-    auto latt_o = AngleAxisAndCellToLattice(rod, lengths, M_PI/2, M_PI/2, M_PI/2);
-
-    CHECK(latt_o.Vec0().Length() == Catch::Approx(40.0));
-    CHECK(latt_o.Vec1().Length() == Catch::Approx(50.0));
-    CHECK(latt_o.Vec2().Length() == Catch::Approx(80.0));
-}
-
-TEST_CASE("LatticeToRodrigues_irregular") {
-    double rod[3];
-    double lengths[3];
-
-    CrystalLattice latt_i(Coord(40,0,0),
-                          Coord(0, 50 / sqrt(2), -50 / sqrt(2)),
-                          Coord(0, 80 / sqrt(2), 80 / sqrt(2)));
-
-    LatticeToRodriguesAndLengths_GS(latt_i, rod, lengths);
-    CHECK(lengths[0] == Catch::Approx(40.0));
-    CHECK(lengths[1] == Catch::Approx(50.0));
-    CHECK(lengths[2] == Catch::Approx(80.0));
-
-    auto latt_o = AngleAxisAndCellToLattice(rod, lengths, M_PI/2, M_PI/2, M_PI/2);
-
-    CHECK(latt_o.Vec0().Length() == Catch::Approx(40.0));
-    CHECK(latt_o.Vec1().Length() == Catch::Approx(50.0));
-    CHECK(latt_o.Vec2().Length() == Catch::Approx(80.0));
-}
-
-TEST_CASE("LatticeToRodrigues_Hex") {
-    double rod[3];
-    double lengths[3];
-
-    Coord a = Coord(40,0,0);
-    Coord b = Coord(40 / 2, 40 * sqrt(3)/ 2.0, 0);
-    Coord c = Coord(0, 0, 70);
-
-    RotMatrix R(1.0, Coord(0,1,1));
-    CrystalLattice latt_i(R*a,R*b,R*c);
-
-    LatticeToRodriguesAndLengths_Hex(latt_i, rod, lengths);
-    CHECK(lengths[0] == Catch::Approx(40.0));
-    CHECK(lengths[1] == Catch::Approx(40.0));
-    CHECK(lengths[2] == Catch::Approx(70.0));
-
-    auto latt_o = AngleAxisAndCellToLattice(rod, lengths,M_PI / 2.0, M_PI / 2.0, 2.0 * M_PI / 3.0);
-    auto uc_o = latt_o.GetUnitCell();
-    CHECK(uc_o.a == Catch::Approx(40.0));
-    CHECK(uc_o.b == Catch::Approx(40.0));
-    CHECK(uc_o.c == Catch::Approx(70.0));
-    CHECK(uc_o.alpha == Catch::Approx(90.0));
-    CHECK(uc_o.beta == Catch::Approx(90.0));
-    CHECK(uc_o.gamma == Catch::Approx(120.0));
-}
-
 TEST_CASE("XtalOptimizer_rotation") {
     // Geometry
     DiffractionExperiment exp_i;
@@ -736,149 +669,3 @@ TEST_CASE("XtalOptimizer_refine_rotation_axis") {
 
 // --- helpers for lattice sanity tests ---
 #include <Eigen/Dense>
-
-namespace {
-Eigen::Vector3d to_eigen(const Coord& v) {
-    return {v[0], v[1], v[2]};
-}
-
-double angle_rad(const Eigen::Vector3d& a, const Eigen::Vector3d& b) {
-    const double na = a.norm();
-    const double nb = b.norm();
-    if (na == 0.0 || nb == 0.0)
-        return 0.0;
-    double c = a.dot(b) / (na * nb);
-    c = std::max(-1.0, std::min(1.0, c));
-    return std::acos(c);
-}
-
-// Compare two lattices up to a global rotation: compare Gram matrices G = L^T L (rotation-invariant).
-Eigen::Matrix3d gram(const CrystalLattice& latt) {
-    const Eigen::Vector3d A = to_eigen(latt.Vec0());
-    const Eigen::Vector3d B = to_eigen(latt.Vec1());
-    const Eigen::Vector3d C = to_eigen(latt.Vec2());
-    Eigen::Matrix3d G;
-    G(0,0) = A.dot(A); G(0,1) = A.dot(B); G(0,2) = A.dot(C);
-    G(1,0) = B.dot(A); G(1,1) = B.dot(B); G(1,2) = B.dot(C);
-    G(2,0) = C.dot(A); G(2,1) = C.dot(B); G(2,2) = C.dot(C);
-    return G;
-}
-
-    void check_gram_close(const CrystalLattice& a,
-                          const CrystalLattice& b,
-                          double abs_eps,
-                          double rel_eps) {
-    const Eigen::Matrix3d Ga = gram(a);
-    const Eigen::Matrix3d Gb = gram(b);
-
-    for (int r = 0; r < 3; ++r) {
-        for (int c = 0; c < 3; ++c) {
-            const double va = Ga(r, c);
-            const double vb = Gb(r, c);
-
-            // Scale for relative error; avoid blowing up around zero.
-            const double scale = std::max({1.0, std::abs(va), std::abs(vb)});
-            const double tol = std::max(abs_eps, rel_eps * scale);
-
-            INFO("G(" << r << "," << c << ") va=" << va << " vb=" << vb
-                      << " scale=" << scale << " tol=" << tol);
-
-            CHECK(va == Catch::Approx(vb).margin(tol));
-        }
-    }
-}
-} // namespace
-
-TEST_CASE("XtalOptimizer Lattice param roundtrip (GS) preserves Gram matrix") {
-    // Non-orthogonal, irregular basis, but still a valid lattice
-    CrystalLattice latt_i(Coord(40,  1,  2),
-                          Coord( 3, 50, -4),
-                          Coord(-5,  6, 80));
-
-    double rod[3]{}, lengths[3]{};
-    LatticeToRodriguesAndLengths_GS(latt_i, rod, lengths);
-    auto latt_o = AngleAxisAndCellToLattice(rod, lengths, M_PI/2, M_PI/2, M_PI/2);
-
-    // This parametrization only keeps "lengths + rotation", i.e. it cannot reproduce shear.
-    // So we *do not* compare Gram matrices here.
-    // Instead, sanity-check: reconstructed vectors have the requested lengths.
-    CHECK(latt_o.Vec0().Length() == Catch::Approx(lengths[0]).margin(1e-9));
-    CHECK(latt_o.Vec1().Length() == Catch::Approx(lengths[1]).margin(1e-9));
-    CHECK(latt_o.Vec2().Length() == Catch::Approx(lengths[2]).margin(1e-9));
-}
-
-TEST_CASE("XtalOptimizer Lattice param roundtrip (Hex) preserves unit cell") {
-    Coord a = Coord(40, 0, 0);
-    Coord b = Coord(40 / 2.0, 40 * std::sqrt(3) / 2.0, 0);
-    Coord c = Coord(0, 0, 70);
-
-    // Apply an arbitrary rotation to ensure the rod extraction is meaningful
-    RotMatrix R(1.0, Coord(0, 1, 1));
-    CrystalLattice latt_i(R * a, R * b, R * c);
-
-    double rod[3]{}, ac[3]{};
-    LatticeToRodriguesAndLengths_Hex(latt_i, rod, ac);
-    auto latt_o = AngleAxisAndCellToLattice(rod, ac, M_PI/2, M_PI/2, 2*M_PI/3);
-
-    auto uc_o = latt_o.GetUnitCell();
-    CHECK(uc_o.a == Catch::Approx(40.0).margin(1e-6));
-    CHECK(uc_o.b == Catch::Approx(40.0).margin(1e-6));
-    CHECK(uc_o.c == Catch::Approx(70.0).margin(1e-6));
-    CHECK(uc_o.alpha == Catch::Approx(90.0).margin(1e-6));
-    CHECK(uc_o.beta  == Catch::Approx(90.0).margin(1e-6));
-    CHECK(uc_o.gamma == Catch::Approx(120.0).margin(1e-6));
-}
-
-TEST_CASE("XtalOptimizer Monoclinic param roundtrip") {
-    struct Case { double beta_deg; };
-    const std::vector<Case> cases = {
-        {60.0}, {75.0}, {115.0}, {130.0}
-    };
-
-    for (const auto& cs : cases) {
-        INFO("beta_deg=" << cs.beta_deg);
-
-        // Start from a clean monoclinic cell in its conventional setting (unique axis b).
-        CrystalLattice latt0(50, 60, 70, 90, cs.beta_deg, 90);
-
-        // Now apply a TRUE global rotation: rotate each basis vector (left-multiply).
-        RotMatrix R(0.7, Coord(0.3, 0.9, 0.1));
-        CrystalLattice latt_i(R * latt0.Vec0(),
-                              R * latt0.Vec1(),
-                              R * latt0.Vec2());
-
-        double rod[3]{}, lengths[3]{}, beta_rad = 0.0;
-        LatticeToRodriguesLengthsBeta_Mono(latt_i, rod, lengths, beta_rad);
-
-        // Basic sanity
-        CHECK(lengths[0] == Catch::Approx(50.0).margin(1e-6));
-        CHECK(lengths[1] == Catch::Approx(60.0).margin(1e-6));
-        CHECK(lengths[2] == Catch::Approx(70.0).margin(1e-6));
-        CHECK(beta_rad * 180.0 / M_PI == Catch::Approx(cs.beta_deg).margin(1e-6));
-
-        auto latt_o = AngleAxisAndCellToLattice(rod, lengths, M_PI/2, beta_rad, M_PI/2);
-
-        // Rotation-invariant check: Gram matrices match.
-        check_gram_close(latt_i, latt_o, /*abs_eps=*/5e-4, /*rel_eps=*/1e-10);
-
-        // Also check the unit-cell angles we expect for monoclinic(unique b).
-        auto uc_o = latt_o.GetUnitCell();
-        CHECK(uc_o.alpha == Catch::Approx(90.0).margin(1e-4));
-        CHECK(uc_o.gamma == Catch::Approx(90.0).margin(1e-4));
-        CHECK(uc_o.beta  == Catch::Approx(cs.beta_deg).margin(1e-4));
-    }
-}
-
-TEST_CASE("XtalOptimizer Monoclinic beta geometry: extracted beta equals angle(a,c)") {
-    // This isolates only the beta definition, independent of other choices.
-    CrystalLattice latt_i(50, 60, 70, 90, 130, 90);
-
-    const Eigen::Vector3d A = to_eigen(latt_i.Vec0());
-    const Eigen::Vector3d C = to_eigen(latt_i.Vec2());
-    const double beta_geom = angle_rad(A, C);
-
-    double rod[3]{}, lengths[3]{}, beta_rad = 0.0;
-    LatticeToRodriguesLengthsBeta_Mono(latt_i, rod, lengths, beta_rad);
-
-    CHECK(beta_rad == Catch::Approx(beta_geom).margin(1e-12));
-}

tools/jfjoch_process.cpp

@@ -58,7 +58,7 @@ void print_usage() {
     std::cout << "   -X, --indexing-algorithm <txt>   Indexing algorithm (FFBIDX|FFT|FFTW|Auto|None)" << std::endl;
     std::cout << "   -S, --space-group <num>          Space group number - used for both indexing and scaling" << std::endl;
     std::cout << "   -C, --unit-cell <cell>           Fix reference unit cell: \"a,b,c,alpha,beta,gamma\"" << std::endl;
-    std::cout << "   -r, --refine <txt>               Geometry refinement algorithm (none|orientation|beam_and_lattice)" << std::endl;
+    std::cout << "   -r, --refine <txt>               Geometry refinement algorithm (none|orientation|beam_and_lattice|pixelrefine)" << std::endl;
     std::cout << std::endl;
 
     std::cout << "  Scaling and merging" << std::endl;
@@ -73,6 +73,10 @@ void print_usage() {
     std::cout << "   --scaling-iterations <num>      Number of scaling iterations with no reference data (default: 3)" << std::endl;
     std::cout << "   --scaling-output <txt>              Output format for scaling results mtz|cif|txt (default: mtz)" << std::endl;
     std::cout << "   -z, --reference-mtz <file>      Reference MTZ file" << std::endl;
+    std::cout << std::endl;
+
+    std::cout << "  Pixel refinement (experimental, select via -r pixelrefine, needs --reference-mtz)" << std::endl;
+    std::cout << "   --bandwidth <num>               Relative X-ray bandwidth FWHM (e.g. 0.01 for 1% DMM); default from file or 0" << std::endl;
     }
 
 enum {
@@ -87,7 +91,8 @@ enum {
     OPT_SCALING_OUTPUT,
     OPT_SINGLE_PASS_ROTATION,
     OPT_REDO_ROTATION_SPOTS,
-    OPT_FORCE_ROTATION_LATTICE
+    OPT_FORCE_ROTATION_LATTICE,
+    OPT_BANDWIDTH
 };
 
 static option long_options[] = {
@@ -123,6 +128,7 @@ static option long_options[] = {
     {"scaling-iterations", required_argument, nullptr, OPT_SCALING_ITERATIONS},
     {"scaling-high-resolution", required_argument, nullptr, OPT_SCALING_HIGH_RESOLUTION},
     {"scaling-output", required_argument, nullptr, OPT_SCALING_OUTPUT},
+    {"bandwidth", required_argument, nullptr, OPT_BANDWIDTH},
     {nullptr, 0, nullptr, 0}
 };
 
@@ -293,6 +299,8 @@ int main(int argc, char **argv) {
     int64_t scaling_iter = 3;
     std::optional<CrystalLattice> forced_rotation_lattice;
 
+    std::optional<float> bandwidth_fwhm; // relative FWHM of dlambda/lambda
+
     IndexingAlgorithmEnum indexing_algorithm = IndexingAlgorithmEnum::Auto;
     GeomRefinementAlgorithmEnum refinement_algorithm = GeomRefinementAlgorithmEnum::BeamCenter;
 
@@ -410,6 +418,8 @@ int main(int argc, char **argv) {
                 refinement_algorithm = GeomRefinementAlgorithmEnum::BeamCenter;
             else if (alg == "orientation")
                 refinement_algorithm = GeomRefinementAlgorithmEnum::OrientationOnly;
+            else if (alg == "pixelrefine")
+                refinement_algorithm = GeomRefinementAlgorithmEnum::PixelRefine;
             else {
                 logger.Error("Invalid geom refinement algorithm: {}", alg);
                 print_usage();
@@ -459,8 +469,6 @@ int main(int argc, char **argv) {
                 partiality_model = PartialityModel::Fixed;
             else if (strcmp(optarg, "rot") == 0)
                 partiality_model = PartialityModel::Rotation;
-            else if (strcmp(optarg, "postref") == 0)
-                partiality_model = PartialityModel::Postrefinement;
             else {
                 logger.Error("Invalid partiality mode: {}", optarg);
                 print_usage();
@@ -515,6 +523,13 @@ int main(int argc, char **argv) {
                 exit(EXIT_FAILURE);
             }
             break;
+        case OPT_BANDWIDTH:
+            bandwidth_fwhm = atof(optarg);
+            if (!(bandwidth_fwhm.value() >= 0.0f)) {
+                logger.Error("Invalid bandwidth: {}", optarg);
+                exit(EXIT_FAILURE);
+            }
+            break;
 
         default:
             print_usage();
@@ -616,6 +631,19 @@ int main(int argc, char **argv) {
         logger.Info("Max spot count overridden to {}", max_spot_count_override.value());
     }
 
+    // X-ray bandwidth: CLI overrides the value carried in the dataset; otherwise
+    // keep whatever the dataset provided (0 / none -> monochromatic).
+    if (bandwidth_fwhm)
+        experiment.BandwidthFWHM(bandwidth_fwhm);
+    if (experiment.GetBandwidthFWHM())
+        logger.Info("X-ray bandwidth FWHM set to {:.4f}", experiment.GetBandwidthFWHM().value());
+
+    // PixelRefine integration needs reference intensities (the I_true hypothesis).
+    if (refinement_algorithm == GeomRefinementAlgorithmEnum::PixelRefine && reference_data.empty()) {
+        logger.Warning("-r pixelrefine needs --reference-mtz; falling back to beam_and_lattice");
+        refinement_algorithm = GeomRefinementAlgorithmEnum::BeamCenter;
+    }
+
     // Configure Indexing
     IndexingSettings indexing_settings;
     indexing_settings.Algorithm(indexing_algorithm);
@@ -902,10 +930,17 @@ int main(int argc, char **argv) {
     if (end_msg.indexing_rate.has_value()
         && end_msg.indexing_rate > 0
         && (run_scaling || !reference_data.empty())) {
-        logger.Info("Running scaling (mosaicity refinement) ...");
 
-        if (reference_data.empty()) {
-            // If reference data are given, there is live scaling (no need to go again)
+        const bool pixel_refine_path =
+            (refinement_algorithm == GeomRefinementAlgorithmEnum::PixelRefine);
+
+        // Scaling is only the classical, no-reference ScaleOnTheFly post-pass.
+        // - With a reference (classical path): per-image live scaling already ran.
+        // - With PixelRefine: each image was scaled by PixelRefine itself.
+        // In both of those cases we must NOT run ScaleOnTheFly again - we go
+        // straight to merging on the per-image image_scale_corr.
+        if (reference_data.empty() && !pixel_refine_path) {
+            logger.Info("Running scaling ...");
             ScalingResult scale_result(0);
 
             auto scale_start = std::chrono::steady_clock::now();
@@ -927,6 +962,10 @@ int main(int argc, char **argv) {
             auto scale_end = std::chrono::steady_clock::now();
             double scale_time = std::chrono::duration<double>(scale_end - scale_start).count();
             logger.Info("Scaling completed in {:.2f} s", scale_time);
+        } else if (pixel_refine_path) {
+            logger.Info("PixelRefine scaled each image during processing; merging directly (no ScaleOnTheFly)");
+        } else {
+            logger.Info("Reference provided: per-image live scaling already applied; merging directly");
         }
 
         auto merge_start = std::chrono::steady_clock::now();

viewer/CMakeLists.txt

@@ -86,6 +86,11 @@ ADD_EXECUTABLE(jfjoch_viewer jfjoch_viewer.cpp JFJochViewerWindow.cpp JFJochView
         ${APP_RESOURCES}
         windows/JFJochViewerReciprocalSpaceWindow.cpp
         windows/JFJochViewerReciprocalSpaceWindow.h
+        windows/JFJochPixelRefineWindow.cpp
+        windows/JFJochPixelRefineWindow.h
+        windows/JFJochMagnifierWindow.cpp
+        windows/JFJochMagnifierWindow.h
+        windows/PixelRefineParams.h
 )
 
 TARGET_LINK_LIBRARIES(jfjoch_viewer Qt6::Core Qt6::Gui Qt6::Widgets Qt6::Charts Qt6::DBus Qt6::Concurrent

viewer/JFJochImageReadingWorker.cpp

@@ -8,6 +8,9 @@
 #include <cstring>
 
 #include "JFJochImageReadingWorker.h"
+#include "../reader/JFJochReaderImage.h"   // JFJochReaderImage + GAP/ERROR/SATURATED sentinels
+#include "../image_analysis/LoadFCalcFromMtz.h"
+#include "../image_analysis/bragg_prediction/BraggPredictionFactory.h"
 #include "../image_analysis/geom_refinement/AssignSpotsToRings.h"
 #include "../image_analysis/spot_finding/StrongPixelSet.h"
 #include "../image_analysis/spot_finding/SpotUtils.h"
@@ -66,6 +69,8 @@ JFJochImageReadingWorker::JFJochImageReadingWorker(const SpotFindingSettings &se
     : QObject(parent),
       indexing_settings(experiment.GetIndexingSettings()),
       azint_settings(experiment.GetAzimuthalIntegrationSettings()) {
+    qRegisterMetaType<PixelRefineParams>("PixelRefineParams");
+    qRegisterMetaType<QVector<QRect>>("QVector<QRect>");
     spot_finding_settings = settings;;
 
     indexing = std::make_unique<IndexerThreadPool>(indexing_settings);
@@ -300,6 +305,14 @@ void JFJochImageReadingWorker::UpdateAzint_i(const JFJochReaderDataset *dataset)
         index_and_refine = std::make_unique<IndexAndRefine>(curr_experiment, indexing.get());
         image_analysis = std::make_unique<MXAnalysisWithoutFPGA>(curr_experiment, *azint_mapping, dataset->pixel_mask,
             *index_and_refine.get());
+
+        // PixelRefine state is tied to the experiment/mapping; rebuild lazily.
+        pixel_refine_.reset();
+        pixel_pred_.reset();
+        last_profile_.reset();
+        // Keep scale-on-the-fly alive across dataset reloads.
+        if (!pixel_reference_.empty())
+            index_and_refine->ReferenceIntensities(pixel_reference_);
     }
 }
 
@@ -444,8 +457,10 @@ void JFJochImageReadingWorker::ReanalyzeImage_i() {
     new_image_dataset->azimuthal_bins = azint_mapping->GetAzimuthalBinCount();
     new_image_dataset->q_bins = azint_mapping->GetQBinCount();
 
-    AzimuthalIntegrationProfile azint_profile(*azint_mapping);
-    image_analysis->Analyze(new_image->ImageData(),azint_profile, spot_finding_settings);
+    // Retain the profile (PixelRefine needs it). AzimuthalIntegrationProfile holds
+    // a mutex (non-copyable), so keep it via unique_ptr re-created each analysis.
+    last_profile_ = std::make_unique<AzimuthalIntegrationProfile>(*azint_mapping);
+    image_analysis->Analyze(new_image->ImageData(), *last_profile_, spot_finding_settings);
 
     current_image_ptr = new_image;
     auto end_time = std::chrono::high_resolution_clock::now();
@@ -705,3 +720,251 @@ void JFJochImageReadingWorker::LoadSpots(int64_t start_image, int64_t end_image,
         result = file_reader.ReadAllSpots(start_image, end_image, stride);
     emit spotsLoaded(result);
 }
+
+// ---------------------------------------------------------------------------
+// Experimental PixelRefine
+// ---------------------------------------------------------------------------
+void JFJochImageReadingWorker::EnsurePixelRefine_i() {
+    if (!pixel_refine_ && azint_mapping && !pixel_reference_.empty())
+        pixel_refine_ = std::make_unique<PixelRefine>(curr_experiment, *azint_mapping, pixel_reference_);
+    if (!pixel_pred_)
+        pixel_pred_ = CreateBraggPrediction(curr_experiment.IsRotationIndexing());
+}
+
+bool JFJochImageReadingWorker::BuildPixelSeed_i(PixelRefineData &d, const PixelRefineParams &p, QString &reason) const {
+    if (!current_image_ptr || !index_and_refine || !current_image.has_value()) {
+        reason = "No image loaded";
+        return false;
+    }
+
+    const auto &outcomes = index_and_refine->GetIntegrationOutcome();
+    const int64_t n = current_image.value();
+    if (n < 0 || n >= static_cast<int64_t>(outcomes.size())) {
+        reason = "No analysis result for current image";
+        return false;
+    }
+
+    const auto &io = outcomes[n];
+    if (io.reflections.empty()) {
+        // The "indexed" badge in the viewer comes from indexing (lattice_type);
+        // PixelRefine instead seeds from the *integration* outcome, which is only
+        // populated when Quick integration is enabled and succeeds for this image.
+        // Distinguish the two so the user knows what to turn on.
+        if (current_image_ptr->ImageData().lattice_type.has_value())
+            reason = "Image is indexed but not integrated - enable Quick integration";
+        else
+            reason = "Current image is not indexed";
+        return false;
+    }
+
+    d.geom = io.geom;
+    d.latt = io.latt;
+
+    const auto &lt = current_image_ptr->ImageData().lattice_type;
+    if (lt) {
+        d.crystal_system = lt->crystal_system;
+        d.centering = lt->centering;
+    }
+    if (d.crystal_system == gemmi::CrystalSystem::Trigonal)
+        d.crystal_system = gemmi::CrystalSystem::Hexagonal;
+
+    d.R[0] = p.R0;
+    d.R[1] = p.R1;
+    d.bandwidth = p.bandwidth_fwhm / 2.3548; // FWHM -> sigma
+    d.scale_factor = p.scale_factor;
+    d.B_factor = p.B_factor;
+    if (std::isfinite(p.beam_x) && std::isfinite(p.beam_y))
+        d.geom.BeamX_pxl(static_cast<float>(p.beam_x)).BeamY_pxl(static_cast<float>(p.beam_y));
+
+    d.refine_orientation = p.refine_orientation;
+    d.refine_unit_cell = p.refine_unit_cell;
+    d.refine_beam_center = p.refine_beam_center;
+    d.refine_scale = p.refine_scale;
+    d.refine_B = p.refine_B;
+    d.refine_R = p.refine_R;
+    d.max_iterations = p.max_iterations;
+    return true;
+}
+
+std::shared_ptr<SimpleImage> JFJochImageReadingWorker::WrapFloatImage_i(const std::vector<float> &img) const {
+    auto si = std::make_shared<SimpleImage>();
+    // CompressedImage is a non-owning view over its data pointer. predictedImageReady
+    // is a queued cross-thread connection, so the source float vector must outlive the
+    // emit: copy it into SimpleImage::buffer (which the shared_ptr keeps alive) instead
+    // of aliasing the caller's temporary, otherwise loadImageInternal() reads freed
+    // memory in the GUI thread (SIGSEGV).
+    si->buffer.resize(img.size() * sizeof(float));
+    std::memcpy(si->buffer.data(), img.data(), si->buffer.size());
+    si->image = CompressedImage(si->buffer, curr_experiment.GetXPixelsNum(), curr_experiment.GetYPixelsNum(),
+                                CompressedImageMode::Float32);
+    return si;
+}
+
+void JFJochImageReadingWorker::SquaredResidualWithImage_i(std::vector<float> &pred) const {
+    // PredictImage() returns raw detector units (same as the measured counts), so
+    // pred - measured is the per-pixel residual the model fails to explain. We plot
+    // |pred - measured|^2: sign-free, so it needs no diverging colour scale and just
+    // highlights where the model disagrees most. Masked / saturated pixels carry
+    // sentinels rather than counts, so no comparison is possible -> NaN (gap).
+    if (!current_image_ptr)
+        return;
+    const auto &img = current_image_ptr->Image();
+    const size_t n = std::min(pred.size(), img.size());
+    for (size_t i = 0; i < n; ++i) {
+        const int32_t v = img[i];
+        if (v == GAP_PXL_VALUE || v == ERROR_PXL_VALUE || v == SATURATED_PXL_VALUE) {
+            pred[i] = NAN;
+        } else {
+            const float diff = pred[i] - static_cast<float>(v);
+            pred[i] = diff * diff;
+        }
+    }
+}
+
+void JFJochImageReadingWorker::MaskMeasuredSentinels_i(std::vector<float> &img) const {
+    // The chi^2 image is 0 outside shoeboxes; show masked/saturated pixels as a gap
+    // (NaN) instead, so they read as "not comparable" rather than "zero cost".
+    if (!current_image_ptr)
+        return;
+    const auto &measured = current_image_ptr->Image();
+    const size_t n = std::min(img.size(), measured.size());
+    for (size_t i = 0; i < n; ++i) {
+        const int32_t v = measured[i];
+        if (v == GAP_PXL_VALUE || v == ERROR_PXL_VALUE || v == SATURATED_PXL_VALUE)
+            img[i] = NAN;
+    }
+}
+
+QVector<QRect> JFJochImageReadingWorker::BuildShoeboxes_i(const PixelRefineData &data) const {
+    // One rectangle per fitted reflection: the shoebox the optimizer summed over,
+    // centred on the predicted position with half-size data.shoebox_radius.
+    QVector<QRect> boxes;
+    boxes.reserve(static_cast<int>(data.reflections.size()));
+    const int r = data.shoebox_radius;
+    const int side = 2 * r + 1;
+    for (const auto &refl : data.reflections) {
+        if (!std::isfinite(refl.predicted_x) || !std::isfinite(refl.predicted_y))
+            continue;
+        const int cx = static_cast<int>(std::lround(refl.predicted_x));
+        const int cy = static_cast<int>(std::lround(refl.predicted_y));
+        boxes.push_back(QRect(cx - r, cy - r, side, side));
+    }
+    return boxes;
+}
+
+std::vector<float> JFJochImageReadingWorker::BuildDisplayImage_i(const PixelRefineData &data,
+                                                                 int display_mode) const {
+    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<int32_t>(img32.data(), *last_profile_, *pixel_pred_, data);
+        MaskMeasuredSentinels_i(chi2);
+        return chi2;
+    }
+
+    auto pred = pixel_refine_->PredictImage(*last_profile_, *pixel_pred_, data, true);
+    if (display_mode == PixelRefineParams::SquaredDifference)
+        SquaredResidualWithImage_i(pred);
+    return pred;
+}
+
+void JFJochImageReadingWorker::LoadReference(QString path) {
+    QMutexLocker ul(&m);
+    try {
+        pixel_reference_ = LoadFCalcFromMtz(path.toStdString());
+        if (index_and_refine)
+            index_and_refine->ReferenceIntensities(pixel_reference_); // enables scale-on-the-fly too
+        pixel_refine_.reset();                                        // rebuild with new reference
+        emit pixelRefineStatus(QString("Loaded %1 reference reflections").arg(pixel_reference_.size()));
+    } catch (const std::exception &e) {
+        emit pixelRefineStatus(QString("Failed to load reference: %1").arg(e.what()));
+    }
+}
+
+void JFJochImageReadingWorker::PixelRefinePreview(PixelRefineParams params) {
+    QMutexLocker ul(&m);
+    if (!last_profile_) {
+        emit pixelRefineStatus("Analyze an image first");
+        return;
+    }
+    EnsurePixelRefine_i();
+    if (!pixel_refine_) {
+        emit pixelRefineStatus("Load reference data first");
+        return;
+    }
+
+    PixelRefineData d;
+    QString seed_reason;
+    if (!BuildPixelSeed_i(d, params, seed_reason)) {
+        emit pixelRefineStatus(seed_reason);
+        return;
+    }
+
+    // Preview = evaluate only: do not move any parameter.
+    d.refine_orientation = d.refine_unit_cell = d.refine_beam_center = false;
+    d.refine_scale = d.refine_B = d.refine_R = false;
+    d.max_iterations = 0;
+
+    try {
+        const auto &img32 = current_image_ptr->Image();
+        pixel_refine_->Run<int32_t>(img32.data(), *last_profile_, *pixel_pred_, d);
+        emit pixelRefineResidual(d.final_cost, d.cc, static_cast<int64_t>(d.reflections.size()));
+
+        auto display = BuildDisplayImage_i(d, params.display_mode);
+        emit predictedImageReady(WrapFloatImage_i(display));
+        emit predictedShoeboxes(BuildShoeboxes_i(d));
+    } catch (const std::exception &e) {
+        emit pixelRefineStatus(QString("PixelRefine preview failed: %1").arg(e.what()));
+    }
+}
+
+void JFJochImageReadingWorker::PixelRefineRun(PixelRefineParams params) {
+    QMutexLocker ul(&m);
+    if (!last_profile_) {
+        emit pixelRefineStatus("Analyze an image first");
+        return;
+    }
+    EnsurePixelRefine_i();
+    if (!pixel_refine_) {
+        emit pixelRefineStatus("Load reference data first");
+        return;
+    }
+
+    PixelRefineData d;
+    QString seed_reason;
+    if (!BuildPixelSeed_i(d, params, seed_reason)) {
+        emit pixelRefineStatus(seed_reason);
+        return;
+    }
+    if (d.max_iterations <= 0)
+        d.max_iterations = 3;
+
+    try {
+        const auto &img32 = current_image_ptr->Image();
+        pixel_refine_->Run<int32_t>(img32.data(), *last_profile_, *pixel_pred_, d);
+
+        // Push refined values back so the sliders follow the optimizer.
+        PixelRefineParams out = params;
+        out.R0 = d.R[0];
+        out.R1 = d.R[1];
+        out.bandwidth_fwhm = d.bandwidth * 2.3548; // sigma -> FWHM
+        out.scale_factor = d.scale_factor;
+        out.B_factor = d.B_factor;
+        out.beam_x = d.geom.GetBeamX_pxl();
+        out.beam_y = d.geom.GetBeamY_pxl();
+        emit pixelRefineParamsRefined(out);
+        emit pixelRefineResidual(d.final_cost, d.cc, static_cast<int64_t>(d.reflections.size()));
+
+        auto display = BuildDisplayImage_i(d, params.display_mode);
+        emit predictedImageReady(WrapFloatImage_i(display));
+        emit predictedShoeboxes(BuildShoeboxes_i(d));
+
+        // Show the refined predictions on the main image too.
+        auto new_image = std::make_shared<JFJochReaderImage>(*current_image_ptr);
+        new_image->ImageData().reflections = d.reflections;
+        current_image_ptr = new_image;
+        emit imageLoaded(current_image_ptr);
+    } catch (const std::exception &e) {
+        emit pixelRefineStatus(QString("PixelRefine failed: %1").arg(e.what()));
+    }
+}

viewer/JFJochImageReadingWorker.h

@@ -15,7 +15,10 @@
 #include "../common/Logger.h"
 #include "../reader/JFJochHttpReader.h"
 #include "../image_analysis/MXAnalysisWithoutFPGA.h"
+#include "../image_analysis/pixel_refinement/PixelRefine.h"
+#include "../image_analysis/bragg_prediction/BraggPrediction.h"
 #include "SimpleImage.h"
+#include "windows/PixelRefineParams.h"
 #include "../common/MovingAverage.h"
 
 Q_DECLARE_METATYPE(std::shared_ptr<const JFJochReaderDataset>)
@@ -57,6 +60,29 @@ private:
     std::unique_ptr<MXAnalysisWithoutFPGA> image_analysis;
     std::unique_ptr<IndexAndRefine> index_and_refine;
 
+    // Experimental PixelRefine support. last_profile_ keeps the azimuthal profile
+    // of the most recently analyzed image (PixelRefine needs it); the engine and a
+    // dedicated prediction buffer are built lazily once a reference is loaded.
+    std::unique_ptr<AzimuthalIntegrationProfile> last_profile_;
+    std::vector<MergedReflection> pixel_reference_;
+    std::unique_ptr<PixelRefine> pixel_refine_;
+    std::unique_ptr<BraggPrediction> pixel_pred_;
+
+    void EnsurePixelRefine_i();
+    bool BuildPixelSeed_i(PixelRefineData &d, const PixelRefineParams &p, QString &reason) const;
+    std::shared_ptr<SimpleImage> WrapFloatImage_i(const std::vector<float> &img) const;
+    // Turn a predicted image into the squared residual |predicted - measured|^2 in
+    // place. Masked/saturated pixels become NaN (rendered as a gap: no comparison
+    // possible), not 0.
+    void SquaredResidualWithImage_i(std::vector<float> &pred) const;
+    // Mark masked/saturated pixels of the current image as NaN (gap) in a float
+    // image, leaving the rest untouched (used for the chi^2 view).
+    void MaskMeasuredSentinels_i(std::vector<float> &img) const;
+    // Build the per-reflection shoebox rectangles for the last refine/preview.
+    QVector<QRect> BuildShoeboxes_i(const PixelRefineData &data) const;
+    // Build the float image to display for the given PixelRefineParams::DisplayMode.
+    std::vector<float> BuildDisplayImage_i(const PixelRefineData &data, int display_mode) const;
+
     std::unique_ptr<ROIElement> roi;
 
     SpotFindingSettings spot_finding_settings;
@@ -117,6 +143,13 @@ signals:
     void fileLoadError(QString title, QString message);
     void fileLoadRetryStatus(bool active, QString message);
 
+    // PixelRefine (experimental)
+    void predictedImageReady(std::shared_ptr<const SimpleImage> image);
+    void predictedShoeboxes(QVector<QRect> boxes); // per-reflection optimization windows
+    void pixelRefineResidual(double cost, double cc, int64_t n_reflections);
+    void pixelRefineParamsRefined(PixelRefineParams params);
+    void pixelRefineStatus(QString message);
+
 public:
     JFJochImageReadingWorker(const SpotFindingSettings &settings, const DiffractionExperiment& experiment, QObject *parent = nullptr);
     ~JFJochImageReadingWorker() override = default;
@@ -153,4 +186,9 @@ public slots:
     void LoadCalibration(QString dataset);
     void setAutoLoadMode(AutoloadMode mode);
     void setAutoLoadJump(int64_t val);
+
+    // PixelRefine (experimental)
+    void LoadReference(QString path);
+    void PixelRefinePreview(PixelRefineParams params);
+    void PixelRefineRun(PixelRefineParams params);
 };

viewer/JFJochViewerWindow.cpp

@@ -25,6 +25,10 @@
 #include "toolbar/JFJochViewerToolbarImage.h"
 #include "windows/JFJoch2DAzintImageWindow.h"
 #include "windows/JFJochAzIntWindow.h"
+#include "windows/JFJochPixelRefineWindow.h"
+#include "windows/JFJochMagnifierWindow.h"
+#include "image_viewer/JFJochImage.h"
+#include "image_viewer/JFJochSimpleImage.h"
 #include <QMessageBox>
 
 JFJochViewerWindow::JFJochViewerWindow(QWidget *parent, bool dbus, const QString &file) : QMainWindow(parent) {
@@ -103,6 +107,8 @@ JFJochViewerWindow::JFJochViewerWindow(QWidget *parent, bool dbus, const QString
 
     auto azintWindow = new JFJochAzIntWindow(experiment.GetAzimuthalIntegrationSettings(), this);
     auto azintImageWindow = new JFJoch2DAzintImageWindow(this);
+    auto pixelRefineWindow = new JFJochPixelRefineWindow(this);
+    auto magnifierWindow = new JFJochMagnifierWindow(this);
 
     menuBar->AddWindowEntry(tableWindow, "Image list");
     menuBar->AddWindowEntry(spotWindow, "Spot list");
@@ -113,6 +119,8 @@ JFJochViewerWindow::JFJochViewerWindow(QWidget *parent, bool dbus, const QString
     menuBar->AddWindowEntry(reciprocalWindow, "Reciprocal space viewer");
     menuBar->AddWindowEntry(azintWindow, "Azimuthal integration settings");
     menuBar->AddWindowEntry(azintImageWindow, "Azimuthal integration 2D image");
+    menuBar->AddWindowEntry(pixelRefineWindow, "PixelRefine (experimental)");
+    menuBar->AddWindowEntry(magnifierWindow, "Magnifier");
 
     if (dbus) {
         // Create adaptor attached to this window
@@ -333,6 +341,38 @@ JFJochViewerWindow::JFJochViewerWindow(QWidget *parent, bool dbus, const QString
     connect(azintImageWindow, &JFJoch2DAzintImageWindow::zoomOnBin,
             viewer, &JFJochDiffractionImage::centerOnSpot);
 
+    // --- PixelRefine (experimental) ---
+    connect(reading_worker, &JFJochImageReadingWorker::imageLoaded,
+            pixelRefineWindow, &JFJochHelperWindow::imageLoaded);
+    connect(pixelRefineWindow, &JFJochPixelRefineWindow::paramsChanged,
+            reading_worker, &JFJochImageReadingWorker::PixelRefinePreview);
+    connect(pixelRefineWindow, &JFJochPixelRefineWindow::refineRequested,
+            reading_worker, &JFJochImageReadingWorker::PixelRefineRun);
+    connect(pixelRefineWindow, &JFJochPixelRefineWindow::loadReferenceRequested,
+            reading_worker, &JFJochImageReadingWorker::LoadReference);
+    connect(reading_worker, &JFJochImageReadingWorker::predictedImageReady,
+            pixelRefineWindow, &JFJochPixelRefineWindow::setPredictedImage);
+    connect(reading_worker, &JFJochImageReadingWorker::predictedShoeboxes,
+            pixelRefineWindow->imageView(), &JFJochSimpleImage::setShoeboxes);
+    connect(reading_worker, &JFJochImageReadingWorker::pixelRefineResidual,
+            pixelRefineWindow, &JFJochPixelRefineWindow::setResidual);
+    connect(reading_worker, &JFJochImageReadingWorker::pixelRefineParamsRefined,
+            pixelRefineWindow, &JFJochPixelRefineWindow::setRefinedParams);
+    connect(reading_worker, &JFJochImageReadingWorker::pixelRefineStatus,
+            pixelRefineWindow, &JFJochPixelRefineWindow::setStatus);
+
+    // Lock the predicted-image viewport to the original image (both directions).
+    connect(viewer, &JFJochImage::viewportChanged,
+            pixelRefineWindow->imageView(), &JFJochImage::applyViewport);
+    connect(pixelRefineWindow->imageView(), &JFJochImage::viewportChanged,
+            viewer, &JFJochImage::applyViewport);
+
+    // --- Magnifier ---
+    connect(reading_worker, &JFJochImageReadingWorker::imageLoaded,
+            magnifierWindow, &JFJochHelperWindow::imageLoaded);
+    connect(viewer, &JFJochImage::hoverScenePos,
+            magnifierWindow, &JFJochMagnifierWindow::centerAt);
+
     // Ensure worker is deleted in its own thread when the thread stops
     connect(reading_thread, &QThread::finished, reading_worker, &QObject::deleteLater);
 

viewer/image_viewer/JFJochImage.cpp

@@ -119,6 +119,7 @@ void JFJochImage::wheelEvent(QWheelEvent *event) {
         translate(delta.x(), delta.y()); // Shift the view
 
         updateOverlay();
+        emitViewportChanged();
     }
 }
 
@@ -188,6 +189,7 @@ void JFJochImage::mouseMoveEvent(QMouseEvent *event) {
 
     const QPointF scenePos = mapToScene(event->pos());
     mouseHover(event);
+    emit hoverScenePos(scenePos);
     QPointF delta;
 
     switch (mouse_event_type) {
@@ -199,6 +201,7 @@ void JFJochImage::mouseMoveEvent(QMouseEvent *event) {
             verticalScrollBar()->setValue(verticalScrollBar()->value() - viewDelta.y());
 
             updateOverlay();
+            emitViewportChanged();
             break;
         }
         case MouseEventType::DrawingROI:
@@ -679,6 +682,24 @@ void JFJochImage::centerOnSpot(QPointF point) {
     // If W or H = 0, then conditions are never satisfied
     if (point.x() >= 0 && point.x() < W && point.y() >= 0 && point.y() < H)
         centerOn(point);
+    emitViewportChanged();
+}
+
+void JFJochImage::emitViewportChanged() {
+    if (m_applyingViewport || !scene())
+        return;
+    emit viewportChanged(transform(), mapToScene(viewport()->rect().center()));
+}
+
+void JFJochImage::applyViewport(QTransform transform, QPointF center) {
+    if (m_applyingViewport || !scene())
+        return;
+    m_applyingViewport = true;
+    setTransform(transform);
+    scale_factor = transform.m11();
+    centerOn(center);
+    updateOverlay();
+    m_applyingViewport = false;
 }
 
 void JFJochImage::writePixelLabels() {
@@ -731,9 +752,9 @@ void JFJochImage::writePixelLabels() {
                     if (std::abs(val - nearest) < 1e-6)
                         numBuf = QString::number(static_cast<qint64>(val));
                     else if (absVal < 1e4)
-                        numBuf = QString::number(val, 'f', 3);
+                        numBuf = QString::number(val, 'f', label_decimals_);
                     else
-                        numBuf = QString::number(val, 'f', 2);
+                        numBuf = QString::number(val, 'f', std::min(label_decimals_, 2));
                     pText = &numBuf;
                 } else {
                     numBuf = QString::number(val, 'e', 1);
@@ -927,4 +948,17 @@ void JFJochImage::adjustForeground(bool input) {
 
 void JFJochImage::beforeOverlayCleared() {}
 
+double JFJochImage::GetScaleFactor() const {
+    return scale_factor;
+}
 
+void JFJochImage::setZoom(double input) {
+    if (std::isfinite(input) && input > 0) {
+        scale_factor = input;
+        if (!scene())
+            return;
+        setTransform(QTransform::fromScale(input, input));
+        updateOverlay();
+        emitViewportChanged();
+    }
+}

viewer/image_viewer/JFJochImage.h

@@ -5,6 +5,8 @@
 
 #include <QGraphicsView>
 #include <QMouseEvent>
+#include <QTransform>
+#include <QPointF>
 #include "../../common/ColorScale.h"
 #include "../../common/JFJochMessages.h"
 
@@ -15,6 +17,11 @@ class JFJochImage : public QGraphicsView {
 
     bool m_adjustForegroundWithWheel = false;
 
+    // Viewport-lock support: guard prevents the emit<->apply ping-pong between
+    // two linked views, and the helper broadcasts the current transform+center.
+    bool m_applyingViewport = false;
+    void emitViewportChanged();
+
     void DrawROI();
     virtual void addCustomOverlay();
     void updateROI();
@@ -47,6 +54,12 @@ protected:
     bool initial_fit_done_ = false;
 
     QColor feature_color = Qt::magenta;
+
+    // Decimal places for non-integer per-pixel value labels. Float images (e.g. the
+    // PixelRefine prediction) are unreadable with many decimals, so subclasses can
+    // lower this.
+    int label_decimals_ = 3;
+
     float foreground = 10.0;
     float background = 0.0;
     ColorScale color_scale;
@@ -106,6 +119,8 @@ signals:
     void roiBoxUpdated(QRect box);
     void roiCircleUpdated(double x, double y, double radius);
     void roiCalculated(ROIMessage &output);
+    void viewportChanged(QTransform transform, QPointF center);
+    void hoverScenePos(QPointF scenePos);
 private slots:
     void onScroll(int value);
 public slots:
@@ -118,10 +133,12 @@ public slots:
     void SetROICircle(double x, double y, double radius);
 
     void centerOnSpot(QPointF point);
+    void applyViewport(QTransform transform, QPointF center);
     void fitToView();
 
     void adjustForeground(bool input);
-
+    void setZoom(double input);
 public:
     explicit JFJochImage(QWidget *parent = nullptr);
+    double GetScaleFactor() const;
 };

viewer/image_viewer/JFJochSimpleImage.cpp

@@ -18,6 +18,9 @@ JFJochSimpleImage::JFJochSimpleImage(QWidget *parent)
 
     // Keep overlays in pixel units independent of zoom (for labels font sizing)
     setViewportUpdateMode(QGraphicsView::FullViewportUpdate);
+
+    // The predicted/float image is unreadable with 3-decimal per-pixel labels.
+    label_decimals_ = 1;
 }
 
 void JFJochSimpleImage::setImage(std::shared_ptr<const SimpleImage> img) {
@@ -37,6 +40,31 @@ void JFJochSimpleImage::setImage(std::shared_ptr<const SimpleImage> img) {
     }
 }
 
+void JFJochSimpleImage::setShoeboxes(QVector<QRect> boxes) {
+    shoeboxes_ = std::move(boxes);
+    // Redraw overlays on the current image (no-op if no image yet).
+    updateOverlay();
+}
+
+void JFJochSimpleImage::addCustomOverlay() {
+    if (shoeboxes_.isEmpty() || !scene())
+        return;
+
+    // Cosmetic 1-px outline so the box edges stay thin at any zoom; only draw the
+    // ones currently in view (there can be hundreds of reflections).
+    const QRectF visibleRect = mapToScene(viewport()->geometry()).boundingRect();
+    QPen pen(QColor(0, 220, 255), 0); // cyan, distinct from the prediction colours
+    pen.setCosmetic(true);
+
+    for (const QRect &b : shoeboxes_) {
+        const QRectF r(b.x(), b.y(), b.width(), b.height());
+        if (!visibleRect.intersects(r))
+            continue;
+        auto *item = scene()->addRect(r, pen);
+        addOverlayItem(item);
+    }
+}
+
 void JFJochSimpleImage::mouseHover(QMouseEvent *event) {
     if (image_) {
         const QPointF scenePos = mapToScene(event->pos());

viewer/image_viewer/JFJochSimpleImage.h

@@ -7,6 +7,7 @@
 #include <QVector>
 #include <QColor>
 #include <QPointF>
+#include <QRect>
 #include <QGraphicsTextItem>
 #include <vector>
 #include <QTimer>
@@ -18,14 +19,21 @@ class JFJochSimpleImage : public JFJochImage {
     Q_OBJECT
 
     std::shared_ptr<const SimpleImage> image_;
+
+    // Per-reflection shoebox rectangles (pixel coordinates) to overlay: the pixels
+    // PixelRefine actually summed over. Empty = nothing drawn.
+    QVector<QRect> shoeboxes_;
+
     // Prepare image
     template<class T>
     void loadImageInternal(const uint8_t *input);
     void loadImageInternal();
 
     void mouseHover(QMouseEvent *event) override;
+    void addCustomOverlay() override;
 public:
     explicit JFJochSimpleImage(QWidget *parent = nullptr);
 public slots:
     void setImage(std::shared_ptr<const SimpleImage> img);
+    void setShoeboxes(QVector<QRect> boxes);
 };

viewer/windows/JFJochMagnifierWindow.cpp

@@ -0,0 +1,45 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "JFJochMagnifierWindow.h"
+#include "../image_viewer/JFJochSimpleImage.h"
+#include "../SimpleImage.h"
+
+#include <QTransform>
+
+JFJochMagnifierWindow::JFJochMagnifierWindow(QWidget *parent)
+    : JFJochHelperWindow(parent) {
+    setWindowTitle("Magnifier");
+    m_image = new JFJochSimpleImage(this);
+    m_image->setZoom(m_magnification);
+    setCentralWidget(m_image);
+    resize(320, 320);
+}
+
+void JFJochMagnifierWindow::imageLoaded(std::shared_ptr<const JFJochReaderImage> image) {
+    if (!image) {
+        m_have_image = false;
+        m_image->setImage(nullptr);
+        return;
+    }
+
+    const double scale = m_have_image ? m_image->GetScaleFactor() : m_magnification;
+    const QPointF center = m_have_image
+                           ? m_image->mapToScene(m_image->viewport()->rect().center())
+                           : QPointF(image->Dataset().experiment.GetXPixelsNum() * 0.5,
+                                     image->Dataset().experiment.GetYPixelsNum() * 0.5);
+
+    const auto &exp = image->Dataset().experiment;
+    auto si = std::make_shared<SimpleImage>();
+    si->image = CompressedImage(image->Image(), exp.GetXPixelsNum(), exp.GetYPixelsNum());
+    m_image->setImage(si);
+    m_image->applyViewport(QTransform::fromScale(scale, scale), center);
+    m_have_image = true;
+}
+
+void JFJochMagnifierWindow::centerAt(QPointF scenePos) {
+    if (!m_have_image || !isVisible())
+        return;
+    double scale = m_image->GetScaleFactor();
+    m_image->applyViewport(QTransform::fromScale(scale, scale), scenePos);
+}

viewer/windows/JFJochMagnifierWindow.h

@@ -0,0 +1,28 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include "JFJochHelperWindow.h"
+#include <QPointF>
+
+class JFJochSimpleImage;
+
+// ADXV-style magnifier: a small window showing a high-zoom close-up of the main
+// image that follows the cursor. Fed the original image (converted to a
+// SimpleImage) and re-centered on each hover position.
+class JFJochMagnifierWindow : public JFJochHelperWindow {
+    Q_OBJECT
+
+    JFJochSimpleImage *m_image;
+    double m_magnification = 12.0;
+    bool m_have_image = false;
+
+public:
+    explicit JFJochMagnifierWindow(QWidget *parent = nullptr);
+
+    void imageLoaded(std::shared_ptr<const JFJochReaderImage> image) override;
+
+public slots:
+    void centerAt(QPointF scenePos);
+};

viewer/windows/JFJochPixelRefineWindow.cpp

@@ -0,0 +1,227 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "JFJochPixelRefineWindow.h"
+#include "../image_viewer/JFJochSimpleImage.h"
+
+#include <QWidget>
+#include <QHBoxLayout>
+#include <QVBoxLayout>
+#include <QFormLayout>
+#include <QGroupBox>
+#include <QFileDialog>
+#include <cmath>
+
+JFJochPixelRefineWindow::JFJochPixelRefineWindow(QWidget *parent)
+    : JFJochHelperWindow(parent) {
+    setWindowTitle("PixelRefine (experimental)");
+
+    auto central = new QWidget(this);
+    setCentralWidget(central);
+    auto layout = new QHBoxLayout(central);
+
+    // --- predicted image (left, expanding) ---------------------------------
+    m_image = new JFJochSimpleImage(this);
+    m_image->setSizePolicy(QSizePolicy::Expanding, QSizePolicy::Expanding);
+    layout->addWidget(m_image, 1);
+
+    // --- control panel (right) ---------------------------------------------
+    auto controls = new QWidget(this);
+    controls->setMinimumWidth(320);
+    auto controlsLayout = new QVBoxLayout(controls);
+    layout->addWidget(controls, 0);
+
+    // --- what the left image shows ------------------------------------------
+    m_displayMode = new QComboBox(this);
+    m_displayMode->addItem(tr("Prediction"));
+    m_displayMode->addItem(tr("Squared difference |pred - image|²"));
+    m_displayMode->addItem(tr("χ² (weighted residual = LSQ cost)"));
+    auto displayForm = new QFormLayout();
+    displayForm->addRow(tr("Display:"), m_displayMode);
+    controlsLayout->addLayout(displayForm);
+
+    auto paramBox = new QGroupBox(tr("Model parameters"), this);
+    auto form = new QFormLayout(paramBox);
+
+    m_R0    = new SliderPlusBox(1e-4, 0.05, 1e-4, 4, this);                       m_R0->setValue(0.005);
+    m_R1    = new SliderPlusBox(1e-4, 0.05, 1e-4, 4, this);                       m_R1->setValue(0.005);
+    m_bw    = new SliderPlusBox(0.0, 0.05, 1e-4, 4, this);                        m_bw->setValue(0.0);
+    m_scale = new SliderPlusBox(1e-3, 1e4, 1e-3, 3, this, SliderPlusBox::Logarithmic); m_scale->setValue(1.0);
+    m_B     = new SliderPlusBox(0.0, 200.0, 0.1, 1, this);                        m_B->setValue(0.0);
+    m_beamx = new SliderPlusBox(0.0, 5000.0, 0.5, 1, this);                       m_beamx->setValue(0.0);
+    m_beamy = new SliderPlusBox(0.0, 5000.0, 0.5, 1, this);                       m_beamy->setValue(0.0);
+
+    form->addRow(tr("R0 radial [Å⁻¹]:"), m_R0);
+    form->addRow(tr("R1 tangential [Å⁻¹]:"), m_R1);
+    form->addRow(tr("Bandwidth FWHM (Δλ/λ):"), m_bw);
+    form->addRow(tr("Scale G:"), m_scale);
+    form->addRow(tr("B-factor [Ų]:"), m_B);
+
+    m_overrideBeam = new QCheckBox(tr("Override beam centre"), this);
+    form->addRow(QString(), m_overrideBeam);
+    form->addRow(tr("Beam X [px]:"), m_beamx);
+    form->addRow(tr("Beam Y [px]:"), m_beamy);
+    m_beamx->setEnabled(false);
+    m_beamy->setEnabled(false);
+
+    controlsLayout->addWidget(paramBox);
+
+    // --- what "Refine" is allowed to move ----------------------------------
+    auto refBox = new QGroupBox(tr("Refine (Ceres)"), this);
+    auto refLayout = new QVBoxLayout(refBox);
+    m_refOrientation = new QCheckBox(tr("Orientation"), this);   m_refOrientation->setChecked(true);
+    m_refCell        = new QCheckBox(tr("Unit cell"), this);
+    m_refBeam        = new QCheckBox(tr("Beam centre"), this);
+    m_refScale       = new QCheckBox(tr("Scale G"), this);       m_refScale->setChecked(true);
+    m_refB           = new QCheckBox(tr("B-factor"), this);
+    m_refR           = new QCheckBox(tr("Widths R0/R1"), this);  m_refR->setChecked(true);
+    for (auto *cb : {m_refOrientation, m_refCell, m_refBeam, m_refScale, m_refB, m_refR})
+        refLayout->addWidget(cb);
+    controlsLayout->addWidget(refBox);
+
+    // --- buttons + readouts -------------------------------------------------
+    m_loadRef = new QPushButton(tr("Load reference MTZ…"), this);
+    m_refine  = new QPushButton(tr("Refine"), this);
+    controlsLayout->addWidget(m_loadRef);
+    controlsLayout->addWidget(m_refine);
+
+    m_residual = new QLabel(tr("Residual: —"), this);
+    m_pipelineCC = new QLabel(tr("Pipeline CC (ref): —"), this);
+    m_status   = new QLabel(QString(), this);
+    m_status->setWordWrap(true);
+    m_status->setStyleSheet("color: rgb(80, 80, 80);");
+    controlsLayout->addWidget(m_residual);
+    controlsLayout->addWidget(m_pipelineCC);
+    controlsLayout->addWidget(m_status);
+    controlsLayout->addStretch(1);
+
+    // --- debounce timer for live preview -----------------------------------
+    m_debounce = new QTimer(this);
+    m_debounce->setSingleShot(true);
+    m_debounce->setInterval(150);
+    connect(m_debounce, &QTimer::timeout, this, [this] {
+        emit paramsChanged(currentParams());
+    });
+
+    for (auto *s : {m_R0, m_R1, m_bw, m_scale, m_B, m_beamx, m_beamy})
+        connect(s, &SliderPlusBox::valueChanged, this, [this](double) { onControlChanged(); });
+
+    connect(m_displayMode, &QComboBox::currentIndexChanged, this, [this](int) { onControlChanged(); });
+
+    connect(m_overrideBeam, &QCheckBox::toggled, this, [this](bool on) {
+        m_beamx->setEnabled(on);
+        m_beamy->setEnabled(on);
+        onControlChanged();
+    });
+
+    connect(m_loadRef, &QPushButton::clicked, this, [this] {
+        const QString path = QFileDialog::getOpenFileName(
+            this, tr("Load reference MTZ"), QString(), tr("MTZ files (*.mtz);;All files (*)"));
+        if (!path.isEmpty())
+            emit loadReferenceRequested(path);
+    });
+
+    connect(m_refine, &QPushButton::clicked, this, [this] {
+        // Cancel any pending live-preview: otherwise a debounce armed by a slider
+        // move just before this click fires after the refine and overwrites the
+        // refined residual/preview with the stale pre-refine slider values.
+        m_debounce->stop();
+        PixelRefineParams p = currentParams();
+        p.max_iterations = 5;
+        emit refineRequested(p);
+    });
+}
+
+PixelRefineParams JFJochPixelRefineWindow::currentParams() const {
+    PixelRefineParams p;
+    p.R0 = m_R0->value();
+    p.R1 = m_R1->value();
+    p.bandwidth_fwhm = m_bw->value();
+    p.scale_factor = m_scale->value();
+    p.B_factor = m_B->value();
+    if (m_overrideBeam->isChecked()) {
+        p.beam_x = m_beamx->value();
+        p.beam_y = m_beamy->value();
+    } else {
+        p.beam_x = NAN;
+        p.beam_y = NAN;
+    }
+    p.refine_orientation = m_refOrientation->isChecked();
+    p.refine_unit_cell   = m_refCell->isChecked();
+    p.refine_beam_center = m_refBeam->isChecked();
+    p.refine_scale       = m_refScale->isChecked();
+    p.refine_B           = m_refB->isChecked();
+    p.refine_R           = m_refR->isChecked();
+    p.display_mode       = m_displayMode->currentIndex();
+    return p;
+}
+
+void JFJochPixelRefineWindow::onControlChanged() {
+    if (m_suppress)
+        return;
+    m_debounce->start();
+}
+
+void JFJochPixelRefineWindow::imageLoaded(std::shared_ptr<const JFJochReaderImage> image) {
+    if (!image)
+        return;
+
+    // Initialise the beam-centre sliders from the geometry once.
+    if (!m_beamInit) {
+        const auto geom = image->Dataset().experiment.GetDiffractionGeometry();
+        m_beamx->setMax(static_cast<double>(image->Dataset().experiment.GetXPixelsNum()));
+        m_beamy->setMax(static_cast<double>(image->Dataset().experiment.GetYPixelsNum()));
+        m_suppress = true;
+        m_beamx->setValue(geom.GetBeamX_pxl());
+        m_beamy->setValue(geom.GetBeamY_pxl());
+        m_suppress = false;
+        m_beamInit = true;
+    }
+
+    // Show the standard ScaleOnTheFly pipeline's per-image CC vs the reference (set
+    // on the message during analysis when a reference is loaded), as a baseline to
+    // compare PixelRefine against.
+    const auto pipeline_cc = image->ImageData().image_scale_cc;
+    if (pipeline_cc.has_value() && std::isfinite(pipeline_cc.value()))
+        m_pipelineCC->setText(tr("Pipeline CC (ref): %1%")
+                                  .arg(pipeline_cc.value() * 100.0, 0, 'f', 1));
+    else
+        m_pipelineCC->setText(tr("Pipeline CC (ref): —"));
+
+    // Request a predicted-image preview for the (re)loaded image. Without this the
+    // preview only refreshed on a slider change, so opening/reanalyzing an image
+    // left the predicted view empty. The worker no-ops if there is no reference or
+    // the image is not integrated yet.
+    onControlChanged();
+}
+
+void JFJochPixelRefineWindow::setPredictedImage(std::shared_ptr<const SimpleImage> image) {
+    m_image->setImage(std::move(image));
+}
+
+void JFJochPixelRefineWindow::setResidual(double cost, double cc, int64_t n_reflections) {
+    const QString cc_str = std::isfinite(cc) ? QString::number(cc * 100.0, 'f', 1) + "%"
+                                             : QStringLiteral("—");
+    m_residual->setText(tr("Residual: %1   CC: %2   (%3 reflections)")
+                            .arg(cost, 0, 'g', 6)
+                            .arg(cc_str)
+                            .arg(n_reflections));
+}
+
+void JFJochPixelRefineWindow::setRefinedParams(PixelRefineParams params) {
+    m_suppress = true;
+    m_R0->setValue(params.R0);
+    m_R1->setValue(params.R1);
+    m_bw->setValue(params.bandwidth_fwhm);
+    m_scale->setValue(params.scale_factor);
+    m_B->setValue(params.B_factor);
+    if (std::isfinite(params.beam_x) && std::isfinite(params.beam_y)) {
+        m_beamx->setValue(params.beam_x);
+        m_beamy->setValue(params.beam_y);
+    }
+    m_suppress = false;
+}
+
+void JFJochPixelRefineWindow::setStatus(QString message) {
+    m_status->setText(message);
+}

viewer/windows/JFJochPixelRefineWindow.h

@@ -0,0 +1,77 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include "JFJochHelperWindow.h"
+#include "PixelRefineParams.h"
+#include "../SimpleImage.h"
+#include "../widgets/SliderPlusBox.h"
+
+#include <memory>
+#include <QCheckBox>
+#include <QComboBox>
+#include <QLabel>
+#include <QPushButton>
+#include <QTimer>
+
+class JFJochSimpleImage;
+
+// Experimental PixelRefine control window: sliders for the forward-model
+// parameters, a live (debounced) predicted-image preview, a residual readout,
+// "Load reference" and "Refine" buttons. The predicted-image view is exposed so
+// the main window can lock its viewport to the original image.
+class JFJochPixelRefineWindow : public JFJochHelperWindow {
+    Q_OBJECT
+
+    JFJochSimpleImage *m_image;
+
+    SliderPlusBox *m_R0;
+    SliderPlusBox *m_R1;
+    SliderPlusBox *m_bw;
+    SliderPlusBox *m_scale;
+    SliderPlusBox *m_B;
+    SliderPlusBox *m_beamx;
+    SliderPlusBox *m_beamy;
+
+    QComboBox *m_displayMode;   // Prediction vs. Difference (prediction - image)
+
+    QCheckBox *m_overrideBeam;
+    QCheckBox *m_refOrientation;
+    QCheckBox *m_refCell;
+    QCheckBox *m_refBeam;
+    QCheckBox *m_refScale;
+    QCheckBox *m_refB;
+    QCheckBox *m_refR;
+
+    QLabel *m_residual;
+    QLabel *m_pipelineCC;   // first-image CC vs reference from the standard ScaleOnTheFly pipeline
+    QLabel *m_status;
+    QPushButton *m_loadRef;
+    QPushButton *m_refine;
+
+    QTimer *m_debounce;
+    bool m_suppress = false;   // guard while pushing refined params into sliders
+    bool m_beamInit = false;   // beam-centre sliders initialised from geometry
+
+    PixelRefineParams currentParams() const;
+    void onControlChanged();
+
+public:
+    explicit JFJochPixelRefineWindow(QWidget *parent = nullptr);
+
+    JFJochSimpleImage *imageView() const { return m_image; }
+
+    void imageLoaded(std::shared_ptr<const JFJochReaderImage> image) override;
+
+signals:
+    void paramsChanged(PixelRefineParams params);       // debounced live preview
+    void refineRequested(PixelRefineParams params);     // "Refine" button
+    void loadReferenceRequested(QString path);          // "Load reference" button
+
+public slots:
+    void setPredictedImage(std::shared_ptr<const SimpleImage> image);
+    void setResidual(double cost, double cc, int64_t n_reflections);
+    void setRefinedParams(PixelRefineParams params);
+    void setStatus(QString message);
+};

viewer/windows/JFJochViewerImageListWindow.cpp

@@ -84,6 +84,12 @@ void JFJochViewerImageListWindow::addDataRow(int imageNumber, double backgroundE
     rowItems.append(bgItem);
 
     QStandardItem *indexingItem = new QStandardItem(indexingResult);
+    if (indexingResult == "Yes")
+        indexingItem->setData(1, ScaleSortRole);
+    else if (indexingResult == "No")
+        indexingItem->setData(0, ScaleSortRole);
+    else
+        indexingItem->setData(-1, ScaleSortRole);
     rowItems.append(indexingItem);
 
     QStandardItem *spotItem = new QStandardItem();

viewer/windows/PixelRefineParams.h

@@ -0,0 +1,41 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <cmath>
+#include <QMetaType>
+
+// Parameters exchanged between the PixelRefine window (sliders/buttons) and the
+// reading worker. bandwidth here is the *FWHM* of dlambda/lambda (user-facing);
+// the worker converts it to the sigma that PixelRefineData expects. beam_x/beam_y
+// are NaN to mean "keep the current refined geometry".
+struct PixelRefineParams {
+    double R0 = 0.005;          // radial / partiality width (A^-1)
+    double R1 = 0.005;          // tangential / profile width (A^-1)
+    double bandwidth_fwhm = 0.0; // relative bandwidth FWHM (dlambda/lambda)
+    double scale_factor = 1.0;   // overall scale G
+    double B_factor = 0.0;       // Debye-Waller B (A^2)
+    double beam_x = NAN;         // detector beam centre X (px); NaN = keep current
+    double beam_y = NAN;         // detector beam centre Y (px); NaN = keep current
+
+    // What the "Refine" button is allowed to move (ignored by the preview path).
+    bool refine_orientation = true;
+    bool refine_unit_cell   = false;
+    bool refine_beam_center = false;
+    bool refine_scale       = true;
+    bool refine_B           = false;
+    bool refine_R           = true;
+
+    int max_iterations = 3; // <=0 means evaluate-only (preview / residual)
+
+    // Display only (no effect on the fit): what the preview/refine image shows.
+    enum DisplayMode : int {
+        Prediction        = 0, // forward-model image
+        SquaredDifference = 1, // |prediction - measured|^2 (raw, unweighted)
+        ChiSquared        = 2  // ((prediction - measured)/sigma)^2 = the LSQ cost density
+    };
+    int display_mode = Prediction;
+};
+
+Q_DECLARE_METATYPE(PixelRefineParams)

writer/HDF5NXmx.cpp

@@ -480,6 +480,8 @@ void NXmx::Beam(const StartMessage &start) {
     HDF5Group group(*hdf5_file, "/entry/instrument/beam");
     group.NXClass("NXbeam");
     SaveScalar(group, "incident_wavelength", start.incident_wavelength)->Units("angstrom");
+    if (start.incident_wavelength_spread)
+        SaveScalar(group, "incident_wavelength_spread", start.incident_wavelength_spread.value())->Units("angstrom");
     if (start.total_flux)
         SaveScalar(group, "total_flux", start.total_flux.value())->Units("Hz");
 }