IndexAndRefine: For the moment - symmetry is enforced only if space group and unit cell are both provided

image_analysis/IndexAndRefine.cpp

@@ -17,6 +17,8 @@ IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool
       indexer_(indexer) {
     if (indexer && x.IsRotationIndexing())
         rotation_indexer = std::make_unique<RotationIndexer>(x, *indexer);
+    reflections.resize(x.GetImageNum());
+    mosaicity.resize(x.GetImageNum(), NAN);
 }
 
 IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(DataMessage &msg) {
@@ -63,17 +65,14 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(Data
     auto latt = indexer_result.lattice[0];
     auto sg = experiment.GetGemmiSpaceGroup();
 
-    // If space group provided => always enforce symmetry in refinement
+    // If space group and cell provided => always enforce symmetry in refinement
     // If space group not provided => guess symmetry
-    if (sg) {
-        // If space group provided but no unit cell fixed, it is better to keep triclinic for now
-        if (experiment.GetUnitCell()) {
-            outcome.symmetry = LatticeMessage{
-                .centering = sg->centring_type(),
-                .niggli_class = 0,
-                .crystal_system = sg->crystal_system()
-            };
-        }
+    if (sg && experiment.GetUnitCell()) {
+        outcome.symmetry = LatticeMessage{
+            .centering = sg->centring_type(),
+            .niggli_class = 0,
+            .crystal_system = sg->crystal_system()
+        };
         outcome.lattice_candidate = latt;
     } else {
         auto sym_result = LatticeSearch(latt);
@@ -166,8 +165,10 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg,
     if (experiment.GetGoniometer().has_value()) {
         wedge_deg = experiment.GetGoniometer()->GetWedge_deg() / 2.0;
 
-        if (msg.mosaicity_deg)
+        if (msg.mosaicity_deg) {
             mos_deg = msg.mosaicity_deg.value();
+            mosaicity[msg.number] = mos_deg;
+        }
     }
 
     const BraggPredictionSettings settings_prediction{
@@ -198,16 +199,15 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg,
                   [](const Reflection& a, const Reflection& b) { return a.d < b.d; });
     }
 
+    {
+        std::unique_lock ul(reflections_mutex);
+        reflections[msg.number] = refl_ret; // Image is not processed twice, so thread-safe in principle, but better safe than sorry :)
+    }
+
     msg.reflections = std::move(refl_ret);
 
     CalcISigma(msg);
     CalcWilsonBFactor(msg);
-
-    // Append reflections to the class-wide reflections buffer (thread-safe)
-    {
-        std::unique_lock ul(reflections_mutex);
-        reflections.insert(reflections.end(), msg.reflections.begin(), msg.reflections.end());
-    }
 }
 
 void IndexAndRefine::ProcessImage(DataMessage &msg,
@@ -242,23 +242,23 @@ std::optional<RotationIndexerResult> IndexAndRefine::Finalize() {
 }
 
 std::optional<ScaleMergeResult> IndexAndRefine::ScaleRotationData(const ScaleMergeOptions &opts) const {
-    std::vector<Reflection> snapshot;
-    {
-        std::unique_lock ul(reflections_mutex);
-        snapshot = reflections;          // cheap copy under lock
-    }
+    size_t nrefl = 0;
+    for (const auto &i: reflections)
+        nrefl += i.size();
 
     // Need a reasonable number of reflections to make refinement meaningful
     constexpr size_t kMinReflections = 20;
-    if (snapshot.size() < kMinReflections)
+    if (nrefl < kMinReflections)
         return std::nullopt;
 
     // Build options focused on mosaicity refinement but allow caller override
     ScaleMergeOptions options = opts;
 
     // If the experiment provides a wedge, propagate it
-    if (experiment.GetGoniometer().has_value())
+    if (experiment.GetGoniometer().has_value() && rotation_indexer) {
         options.wedge_deg = experiment.GetGoniometer()->GetWedge_deg();
+        options.mosaicity_init_deg_vec = mosaicity;
+    }
 
     // If caller left space_group unset, try to pick it from the indexed lattice
     if (!options.space_group.has_value()) {
@@ -267,5 +267,5 @@ std::optional<ScaleMergeResult> IndexAndRefine::ScaleRotationData(const ScaleMer
             options.space_group = *sg;
     }
 
-    return ScaleAndMergeReflectionsCeres(snapshot, options);
-}
+    return ScaleAndMergeReflectionsCeres(reflections, options);
+}

image_analysis/IndexAndRefine.h

@@ -45,7 +45,8 @@ class IndexAndRefine {
     };
 
     mutable std::mutex reflections_mutex;
-    std::vector<Reflection> reflections;
+    std::vector<std::vector<Reflection>> reflections;
+    std::vector<double> mosaicity;
 
     IndexingOutcome DetermineLatticeAndSymmetry(DataMessage &msg);
     void RefineGeometryIfNeeded(DataMessage &msg, IndexingOutcome &outcome);

image_analysis/geom_refinement/XtalOptimizer.cpp

@@ -427,8 +427,6 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
                            const std::vector<SpotToSave> &spots,
                            const float tolerance) {
     try {
-        data.latt.Regularize(data.crystal_system);
-
         Coord vec0 = data.latt.Vec0();
         Coord vec1 = data.latt.Vec1();
         Coord vec2 = data.latt.Vec2();
@@ -635,9 +633,6 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
             // Triclinic via the same generic builder
             data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, latt_vec2[0], latt_vec2[1], latt_vec2[2]);
         }
-
-        data.latt.Regularize(data.crystal_system);
-
         return true;
     } catch (...) {
         // Convergence problems, likely not updated

image_analysis/scale_merge/ScaleAndMerge.cpp

@@ -115,59 +115,6 @@ namespace {
         return key;
     }
 
-    /// CrystFEL-like log-scaling residual
-    ///
-    ///   residual = w * [ ln(I_obs) - ln(G) - ln(partiality) - ln(lp) - ln(I_true) ]
-    ///
-    /// Only observations with I_obs > 0 should be fed in (the caller skips the rest).
-    /// G and I_true are constrained to be positive via Ceres lower bounds.
-    struct LogIntensityResidual {
-        LogIntensityResidual(const Reflection &r, double sigma_obs, double wedge_deg, bool refine_partiality)
-            : log_Iobs_(std::log(std::max(static_cast<double>(r.I), 1e-30))),
-              weight_(SafeInv(sigma_obs / std::max(static_cast<double>(r.I), 1e-30), 1.0)),
-              delta_phi_(r.delta_phi_deg),
-              log_lp_(std::log(std::max(SafeInv(static_cast<double>(r.rlp), 1.0), 1e-30))),
-              half_wedge_(wedge_deg / 2.0),
-              c1_(r.zeta / std::sqrt(2.0)),
-              partiality_(r.partiality),
-              refine_partiality_(refine_partiality) {
-        }
-
-        template<typename T>
-        bool operator()(const T *const G,
-                        const T *const mosaicity,
-                        const T *const Itrue,
-                        T *residual) const {
-            T partiality;
-            if (refine_partiality_ && mosaicity[0] != 0.0) {
-                const T arg_plus = T(delta_phi_ + half_wedge_) * T(c1_) / mosaicity[0];
-                const T arg_minus = T(delta_phi_ - half_wedge_) * T(c1_) / mosaicity[0];
-                partiality = (ceres::erf(arg_plus) - ceres::erf(arg_minus)) / T(2.0);
-            } else {
-                partiality = T(partiality_);
-            }
-
-            // Clamp partiality away from zero so log is safe
-            const T min_p = T(1e-30);
-            if (partiality < min_p)
-                partiality = min_p;
-
-            // ln(I_pred) = ln(G) + ln(partiality) + ln(lp) + ln(Itrue)
-            const T log_Ipred = ceres::log(G[0]) + ceres::log(partiality) + T(log_lp_) + ceres::log(Itrue[0]);
-            residual[0] = (log_Ipred - T(log_Iobs_)) * T(weight_);
-            return true;
-        }
-
-        double log_Iobs_;
-        double weight_; // w_i ≈ I_obs / sigma_obs  (relative weight in log-space)
-        double delta_phi_;
-        double log_lp_;
-        double half_wedge_;
-        double c1_;
-        double partiality_;
-        bool refine_partiality_;
-    };
-
     struct IntensityResidual {
         IntensityResidual(const Reflection &r, double sigma_obs, double wedge_deg, bool refine_partiality)
             : Iobs_(static_cast<double>(r.I)),
@@ -238,306 +185,138 @@ namespace {
 
         double inv_sigma_;
     };
-} // namespace
-
-ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection> &observations,
-                                               const ScaleMergeOptions &opt) {
-    ScaleMergeResult out;
 
     struct ObsRef {
         const Reflection *r = nullptr;
         int img_id = 0;
-        int img_slot = -1;
         int hkl_slot = -1;
         double sigma = 0.0;
     };
 
-    std::vector<ObsRef> obs;
-    obs.reserve(observations.size());
-
-    std::unordered_map<int, int> imgIdToSlot;
-    imgIdToSlot.reserve(256);
-
-    std::unordered_map<HKLKey, int, HKLKeyHash> hklToSlot;
-    hklToSlot.reserve(observations.size());
-
-    for (const auto &r: observations) {
-        const double d = SafeD(r.d);
-        if (!std::isfinite(d))
-            continue;
-        if (!std::isfinite(r.I))
-            continue;
-
-        if (!std::isfinite(r.zeta) || r.zeta <= 0.0f)
-            continue;
-        if (!std::isfinite(r.rlp) || r.rlp == 0.0f)
-            continue;
-
-        const double sigma = SafeSigma(static_cast<double>(r.sigma), opt.min_sigma);
-
-        const int img_id = RoundImageId(r.image_number, opt.image_number_rounding);
-
-        int img_slot;
-        {
-            auto it = imgIdToSlot.find(img_id);
-            if (it == imgIdToSlot.end()) {
-                img_slot = static_cast<int>(imgIdToSlot.size());
-                imgIdToSlot.emplace(img_id, img_slot);
-            } else {
-                img_slot = it->second;
-            }
-        }
-
-        int hkl_slot;
-        try {
-            const HKLKey key = CanonicalizeHKLKey(r, opt);
-            auto it = hklToSlot.find(key);
-            if (it == hklToSlot.end()) {
-                hkl_slot = static_cast<int>(hklToSlot.size());
-                hklToSlot.emplace(key, hkl_slot);
-            } else {
-                hkl_slot = it->second;
-            }
-        } catch (...) {
-            continue;
-        }
-
-        ObsRef o;
-        o.r = &r;
-        o.img_id = img_id;
-        o.img_slot = img_slot;
-        o.hkl_slot = hkl_slot;
-        o.sigma = sigma;
-        obs.push_back(o);
-    }
-
-    const int nimg = static_cast<int>(imgIdToSlot.size());
-    const int nhkl = static_cast<int>(hklToSlot.size());
-
-    out.image_scale_g.assign(nimg, 1.0);
-    out.image_ids.assign(nimg, 0);
-
-    for (const auto &kv: imgIdToSlot) {
-        out.image_ids[kv.second] = kv.first;
-    }
-
-    std::vector<double> g(nimg, 1.0);
-    std::vector<double> Itrue(nhkl, 0.0);
-
-    // Mosaicity: always per-image
-    std::vector<double> mosaicity(nimg, opt.mosaicity_init_deg);
-    for (int i = 0; i < nimg && i < opt.mosaicity_init_deg_vec.size(); ++i) {
-        if (opt.mosaicity_init_deg_vec[i] > 0.0)
-            mosaicity[i] = opt.mosaicity_init_deg_vec[i];
-    }
-
-    // Initialize Itrue from per-HKL median of observed intensities
-    {
-        std::vector<std::vector<double> > per_hkl_I(nhkl);
-        for (const auto &o: obs) {
-            per_hkl_I[o.hkl_slot].push_back(static_cast<double>(o.r->I));
-        }
-        for (int h = 0; h < nhkl; ++h) {
-            auto &v = per_hkl_I[h];
-            if (v.empty()) {
-                Itrue[h] = std::max(opt.min_sigma, 1e-6);
-                continue;
-            }
-            std::nth_element(v.begin(), v.begin() + static_cast<long>(v.size() / 2), v.end());
-            double med = v[v.size() / 2];
-            if (!std::isfinite(med) || med <= opt.min_sigma)
-                med = opt.min_sigma;
-            Itrue[h] = med;
-        }
-    }
-
-    ceres::Problem problem;
-
-    const bool refine_partiality = opt.wedge_deg > 0.0;
-
-    std::vector<bool> is_valid_hkl_slot(nhkl, false);
-
-    for (const auto &o: obs) {
-        size_t mos_slot = opt.per_image_mosaicity ? o.img_slot : 0;
-
-        if (opt.log_scaling_residual) {
-            // Log residual requires positive I_obs
-            if (o.r->I <= 0.0f)
-                continue;
-
-            auto *cost = new ceres::AutoDiffCostFunction<LogIntensityResidual, 1, 1, 1, 1>(
-                new LogIntensityResidual(*o.r, o.sigma, opt.wedge_deg, refine_partiality));
-
-            problem.AddResidualBlock(cost,
-                                     nullptr,
-                                     &g[o.img_slot],
-                                     &mosaicity[mos_slot],
-                                     &Itrue[o.hkl_slot]);
-            is_valid_hkl_slot[o.hkl_slot] = true;
-        } else {
-            auto *cost = new ceres::AutoDiffCostFunction<IntensityResidual, 1, 1, 1, 1>(
-                new IntensityResidual(*o.r, o.sigma, opt.wedge_deg, refine_partiality));
-
-            problem.AddResidualBlock(cost,
-                                     nullptr,
-                                     &g[o.img_slot],
-                                     &mosaicity[mos_slot],
-                                     &Itrue[o.hkl_slot]);
-            is_valid_hkl_slot[o.hkl_slot] = true;
-        }
-    }
-
-    if (opt.smoothen_g) {
-        for (int i = 0; i < nimg - 2; ++i) {
-            auto *cost = new ceres::AutoDiffCostFunction<SmoothnessRegularizationResidual, 1, 1, 1, 1>(
-                new SmoothnessRegularizationResidual(0.05));
-
-            problem.AddResidualBlock(cost, nullptr,
-                                     &g[i],
-                                     &g[i + 1],
-                                     &g[i + 2]);
-        }
-    }
-
-    if (opt.smoothen_mos && opt.per_image_mosaicity) {
-        for (int i = 0; i < nimg - 2; ++i) {
-            auto *cost = new ceres::AutoDiffCostFunction<SmoothnessRegularizationResidual, 1, 1, 1, 1>(
-                new SmoothnessRegularizationResidual(0.05));
-
-            problem.AddResidualBlock(cost, nullptr,
-                                     &mosaicity[i],
-                                     &mosaicity[i + 1],
-                                     &mosaicity[i + 2]);
-        }
-    }
-
-    // Scaling factors must be always positive
-    for (int i = 0; i < nimg; ++i)
-        problem.SetParameterLowerBound(&g[i], 0, 1e-12);
-
-
-    // For log residual, Itrue must stay positive
-    if (opt.log_scaling_residual) {
-        for (int h = 0; h < nhkl; ++h) {
-            if (is_valid_hkl_slot[h])
-                problem.SetParameterLowerBound(&Itrue[h], 0, 1e-12);
-        }
-    }
-
-    // Mosaicity refinement + bounds
-    if (!opt.refine_mosaicity) {
-        for (int i = 0; i < (opt.per_image_mosaicity ? nimg : 1); ++i)
-            problem.SetParameterBlockConstant(&mosaicity[i]);
-    } else {
-        for (int i = 0; i < (opt.per_image_mosaicity ? nimg : 1); ++i) {
-            problem.SetParameterLowerBound(&mosaicity[i], 0, opt.mosaicity_min_deg);
-            problem.SetParameterUpperBound(&mosaicity[i], 0, opt.mosaicity_max_deg);
-        }
-    }
-
-    // use all available threads
-    unsigned int hw = std::thread::hardware_concurrency();
-    if (hw == 0)
-        hw = 1; // fallback
-
-    ceres::Solver::Options options;
-
-    options.linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY;
-    options.minimizer_progress_to_stdout = true;
-    options.max_num_iterations = opt.max_num_iterations;
-    options.max_solver_time_in_seconds = opt.max_solver_time_s;
-    options.num_threads = static_cast<int>(hw);
-
-    ceres::Solver::Summary summary;
-    ceres::Solve(options, &problem, &summary);
-
-    // --- Export per-image results ---
-    for (int i = 0; i < nimg; ++i)
-        out.image_scale_g[i] = g[i];
-
-    out.mosaicity_deg.resize(nimg);
-    for (int i = 0; i < nimg; ++i)
-        out.mosaicity_deg[i] = opt.per_image_mosaicity ? mosaicity[i] : mosaicity[0];
-
-    std::vector<HKLKey> slotToHKL(nhkl);
-    for (const auto &kv: hklToSlot) {
-        slotToHKL[kv.second] = kv.first;
-    }
-
-    out.merged.resize(nhkl);
-    for (int h = 0; h < nhkl; ++h) {
-        out.merged[h].h = slotToHKL[h].h;
-        out.merged[h].k = slotToHKL[h].k;
-        out.merged[h].l = slotToHKL[h].l;
-        out.merged[h].I = Itrue[h];
-        out.merged[h].sigma = 0.0;
-        out.merged[h].d = 0.0;
-    }
-
-    // Populate d from median of observations per HKL
-    {
-        std::vector<std::vector<double> > per_hkl_d(nhkl);
-        for (const auto &o: obs) {
-            const double d_val = static_cast<double>(o.r->d);
-            if (std::isfinite(d_val) && d_val > 0.0)
-                per_hkl_d[o.hkl_slot].push_back(d_val);
-        }
-        for (int h = 0; h < nhkl; ++h) {
-            auto &v = per_hkl_d[h];
-            if (!v.empty()) {
-                std::nth_element(v.begin(), v.begin() + static_cast<long>(v.size() / 2), v.end());
-                out.merged[h].d = v[v.size() / 2];
-            }
-        }
-    }
-
-    std::cout << summary.FullReport() << std::endl;
-
-    // ---- Compute corrected observations once (used for both merging and statistics) ----
     struct CorrectedObs {
         int hkl_slot;
         double I_corr;
         double sigma_corr;
     };
-    std::vector<CorrectedObs> corr_obs;
-    corr_obs.reserve(obs.size());
 
-    {
-        const double half_wedge = opt.wedge_deg / 2.0;
+    void scale(const ScaleMergeOptions &opt,
+               std::vector<double> &g,
+               std::vector<double> &mosaicity,
+               const std::vector<uint8_t> &image_slot_used,
+               bool rotation_crystallography,
+               size_t nhkl,
+               const std::vector<ObsRef> &obs) {
+        ceres::Problem problem;
+
+        std::vector<double> Itrue(nhkl, 0.0);
+
+        // Initialize Itrue from per-HKL median of observed intensities
+        {
+            std::vector<std::vector<double> > per_hkl_I(nhkl);
+            for (const auto &o: obs) {
+                per_hkl_I[o.hkl_slot].push_back(static_cast<double>(o.r->I));
+            }
+            for (int h = 0; h < nhkl; ++h) {
+                auto &v = per_hkl_I[h];
+                if (v.empty()) {
+                    Itrue[h] = std::max(opt.min_sigma, 1e-6);
+                    continue;
+                }
+                std::nth_element(v.begin(), v.begin() + static_cast<long>(v.size() / 2), v.end());
+                double med = v[v.size() / 2];
+                if (!std::isfinite(med) || med <= opt.min_sigma)
+                    med = opt.min_sigma;
+                Itrue[h] = med;
+            }
+        }
+
+        std::vector<bool> is_valid_hkl_slot(nhkl, false);
 
         for (const auto &o: obs) {
-            const Reflection &r = *o.r;
-            const double lp = SafeInv(static_cast<double>(r.rlp), 1.0);
-            const double G_i = g[o.img_slot];
+            auto *cost = new ceres::AutoDiffCostFunction<IntensityResidual, 1, 1, 1, 1>(
+                new IntensityResidual(*o.r, o.sigma, opt.wedge_deg.value_or(0.0), rotation_crystallography));
+            problem.AddResidualBlock(cost,
+                                     nullptr,
+                                     &g[o.img_id],
+                                     &mosaicity[o.img_id],
+                                     &Itrue[o.hkl_slot]);
+            is_valid_hkl_slot[o.hkl_slot] = true;
+        }
 
-            // Compute partiality with refined mosaicity
-            double partiality;
-            const size_t mos_slot = opt.per_image_mosaicity ? o.img_slot : 0;
-            if (refine_partiality && mosaicity[mos_slot] > 0.0) {
-                const double c1 = r.zeta / std::sqrt(2.0);
-                const double arg_plus = (r.delta_phi_deg + half_wedge) * c1 / mosaicity[mos_slot];
-                const double arg_minus = (r.delta_phi_deg - half_wedge) * c1 / mosaicity[mos_slot];
-                partiality = (std::erf(arg_plus) - std::erf(arg_minus)) / 2.0;
-            } else {
-                partiality = r.partiality;
+        for (int i = 0; i < g.size(); ++i) {
+            if (image_slot_used[i]) {
+                auto *cost = new ceres::AutoDiffCostFunction<ScaleRegularizationResidual, 1, 1>(
+                    new ScaleRegularizationResidual(0.05));
+                problem.AddResidualBlock(cost, nullptr, &g[i]);
+            }
+        }
+
+        if (rotation_crystallography) {
+            if (opt.smoothen_g) {
+                for (int i = 0; i < g.size() - 2; ++i) {
+                    if (image_slot_used[i] && image_slot_used[i + 1] && image_slot_used[i + 2]) {
+                        auto *cost = new ceres::AutoDiffCostFunction<SmoothnessRegularizationResidual, 1, 1, 1, 1>(
+                            new SmoothnessRegularizationResidual(0.05));
+
+                        problem.AddResidualBlock(cost, nullptr, &g[i], &g[i + 1], &g[i + 2]);
+                    }
+                }
             }
 
-            if (partiality <= opt.min_partiality_for_merge)
-                continue;
-            const double correction = G_i * partiality * lp;
-            if (correction <= 0.0)
-                continue;
+            if (opt.smoothen_mos && opt.refine_mosaicity) {
+                for (int i = 0; i < mosaicity.size() - 2; ++i) {
+                    if (image_slot_used[i] && image_slot_used[i + 1] && image_slot_used[i + 2]) {
+                        auto *cost = new ceres::AutoDiffCostFunction<SmoothnessRegularizationResidual, 1, 1, 1, 1>(
+                            new SmoothnessRegularizationResidual(0.05));
 
-            corr_obs.push_back({
-                o.hkl_slot,
-                static_cast<double>(r.I) / correction,
-                o.sigma / correction
-            });
+                        problem.AddResidualBlock(cost, nullptr, &mosaicity[i], &mosaicity[i + 1], &mosaicity[i + 2]);
+                    }
+                }
+            }
         }
+
+        // Scaling factors must be always positive
+        for (int i = 0; i < g.size(); i++) {
+            if (image_slot_used[i])
+                problem.SetParameterLowerBound(&g[i], 0, 1e-12);
+        }
+
+        // Mosaicity refinement + bounds
+        if (!opt.refine_mosaicity) {
+            for (int i = 0; i < mosaicity.size(); ++i) {
+                if (image_slot_used[i])
+                    problem.SetParameterBlockConstant(&mosaicity[i]);
+            }
+        } else {
+            for (int i = 0; i < mosaicity.size(); ++i) {
+                if (image_slot_used[i]) {
+                    problem.SetParameterLowerBound(&mosaicity[i], 0, opt.mosaicity_min_deg);
+                    problem.SetParameterUpperBound(&mosaicity[i], 0, opt.mosaicity_max_deg);
+                }
+            }
+        }
+
+        // use all available threads
+        unsigned int hw = std::thread::hardware_concurrency();
+        if (hw == 0)
+            hw = 1; // fallback
+
+        ceres::Solver::Options options;
+
+        options.linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY;
+        options.minimizer_progress_to_stdout = true;
+        options.max_num_iterations = opt.max_num_iterations;
+        options.max_solver_time_in_seconds = opt.max_solver_time_s;
+        options.num_threads = static_cast<int>(hw);
+
+        ceres::Solver::Summary summary;
+        ceres::Solve(options, &problem, &summary);
+        std::cout << summary.FullReport() << std::endl;
     }
 
-    // ---- Merge (XDS/XSCALE style: inverse-variance weighted mean) ----
-    {
+    void merge(size_t nhkl, ScaleMergeResult &out, const std::vector<CorrectedObs> &corr_obs) {
+        // ---- Merge (XDS/XSCALE style: inverse-variance weighted mean) ----
+
         struct HKLAccum {
             double sum_wI = 0.0;
             double sum_w = 0.0;
@@ -573,6 +352,7 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection> &ob
         }
     }
 
+    void stats(const ScaleMergeOptions &opt, size_t nhkl, ScaleMergeResult &out, const std::vector<CorrectedObs> &corr_obs)
     // ---- Compute per-shell merging statistics ----
     {
         constexpr int kStatShells = 10;
@@ -582,6 +362,8 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection> &ob
         for (int h = 0; h < nhkl; ++h) {
             const auto d = static_cast<float>(out.merged[h].d);
             if (std::isfinite(d) && d > 0.0f) {
+                if (opt.d_min_limit_A > 0.0 && d < static_cast<float>(opt.d_min_limit_A))
+                    continue;
                 stat_d_min = std::min(stat_d_min, d);
                 stat_d_max = std::max(stat_d_max, d);
             }
@@ -600,6 +382,8 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection> &ob
             for (int h = 0; h < nhkl; ++h) {
                 const auto d = static_cast<float>(out.merged[h].d);
                 if (std::isfinite(d) && d > 0.0f) {
+                    if (opt.d_min_limit_A > 0.0 && d < static_cast<float>(opt.d_min_limit_A))
+                        continue;
                     auto s = stat_shells.GetShell(d);
                     if (s.has_value())
                         hkl_shell[h] = s.value();
@@ -725,6 +509,187 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection> &ob
             }
         }
     }
+    void calc_obs(const ScaleMergeOptions &opt,
+               std::vector<double> &g,
+               std::vector<double> &mosaicity,
+               bool rotation_crystallography,
+               const std::vector<ObsRef> &obs,
+               std::vector<CorrectedObs> &corr_obs) {
+
+        // ---- Compute corrected observations once (used for both merging and statistics) ----
+        const double half_wedge = opt.wedge_deg.value_or(0.0) / 2.0;
+
+        for (const auto &o: obs) {
+            const Reflection &r = *o.r;
+            const double lp = SafeInv(static_cast<double>(r.rlp), 1.0);
+            const double G_i = g[o.img_id];
+
+            // Compute partiality with refined mosaicity
+            double partiality;
+            if (rotation_crystallography && mosaicity[o.img_id] > 0.0) {
+                const double c1 = r.zeta / std::sqrt(2.0);
+                const double arg_plus = (r.delta_phi_deg + half_wedge) * c1 / mosaicity[o.img_id];
+                const double arg_minus = (r.delta_phi_deg - half_wedge) * c1 / mosaicity[o.img_id];
+                partiality = (std::erf(arg_plus) - std::erf(arg_minus)) / 2.0;
+            } else {
+                partiality = r.partiality;
+            }
+
+            if (partiality <= opt.min_partiality_for_merge)
+                continue;
+            const double correction = G_i * partiality * lp;
+            if (correction <= 0.0)
+                continue;
+
+            corr_obs.push_back({
+                o.hkl_slot,
+                static_cast<double>(r.I) / correction,
+                o.sigma / correction
+            });
+        }
+    }
+
+    void proc_obs(const std::vector<std::vector<Reflection> > &observations,
+        const ScaleMergeOptions &opt,
+        std::vector<uint8_t> &image_slot_used,
+        std::vector<ObsRef> &obs,
+        std::unordered_map<HKLKey, int, HKLKeyHash> &hklToSlot
+        ) {
+        for (int i = 0; i < observations.size(); i++) {
+            for (const auto &r: observations[i]) {
+                const double d = SafeD(r.d);
+                if (!std::isfinite(d))
+                    continue;
+                if (!std::isfinite(r.I))
+                    continue;
+
+                if (opt.d_min_limit_A > 0.0 && d < opt.d_min_limit_A)
+                    continue;
+
+                if (!std::isfinite(r.zeta) || r.zeta <= 0.0f)
+                    continue;
+                if (!std::isfinite(r.rlp) || r.rlp == 0.0f)
+                    continue;
+
+                const double sigma = SafeSigma(r.sigma, opt.min_sigma);
+
+                const int img_id = i / opt.image_cluster;
+                image_slot_used[img_id] = 1;
+
+                int hkl_slot;
+                try {
+                    const HKLKey key = CanonicalizeHKLKey(r, opt);
+                    auto it = hklToSlot.find(key);
+                    if (it == hklToSlot.end()) {
+                        hkl_slot = static_cast<int>(hklToSlot.size());
+                        hklToSlot.emplace(key, hkl_slot);
+                    } else {
+                        hkl_slot = it->second;
+                    }
+                } catch (...) {
+                    continue;
+                }
+
+                ObsRef o;
+                o.r = &r;
+                o.img_id = img_id;
+                o.hkl_slot = hkl_slot;
+                o.sigma = sigma;
+                obs.push_back(o);
+            }
+        }
+
+    }
+} // namespace
+
+ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<std::vector<Reflection> > &observations,
+                                               const ScaleMergeOptions &opt) {
+    if (opt.image_cluster <= 0)
+        throw std::invalid_argument("image_cluster must be positive");
+
+    const bool rotation_crystallography = opt.wedge_deg.has_value();
+
+    size_t nrefl = 0;
+    for (const auto &i: observations)
+        nrefl += i.size();
+
+    std::vector<ObsRef> obs;
+    std::vector<CorrectedObs> corr_obs;
+
+    obs.reserve(nrefl);
+    corr_obs.reserve(nrefl);
+
+    std::unordered_map<HKLKey, int, HKLKeyHash> hklToSlot;
+    hklToSlot.reserve(nrefl);
+
+    size_t n_image_slots = observations.size() / opt.image_cluster +
+        (observations.size() % opt.image_cluster > 0 ? 1 : 0);
+
+    std::vector<uint8_t> image_slot_used(n_image_slots, 0);
+
+    proc_obs(observations, opt, image_slot_used, obs, hklToSlot);
+
+    const int nhkl = static_cast<int>(hklToSlot.size());
+
+    std::vector<double> g(n_image_slots, 1.0);
+    std::vector<double> mosaicity(n_image_slots, opt.mosaicity_init_deg);
+    for (int i = 0; i < n_image_slots; i++) {
+        if (!image_slot_used[i]) {
+            mosaicity[i] = NAN;
+            g[i] = NAN;
+        } else if (opt.mosaicity_init_deg_vec.size() > i && std::isfinite(opt.mosaicity_init_deg_vec[i])) {
+            mosaicity[i] = opt.mosaicity_init_deg_vec[i];
+        }
+    }
+
+    scale(opt, g, mosaicity, image_slot_used, rotation_crystallography, nhkl, obs);
+
+    ScaleMergeResult out;
+
+    out.image_scale_g.resize(observations.size(), NAN);
+    out.mosaicity_deg.resize(observations.size(), NAN);
+    for (int i = 0; i < observations.size(); i++) {
+        size_t img_slot = i / opt.image_cluster;
+        if (image_slot_used[img_slot]) {
+            out.image_scale_g[i] = g[img_slot];
+            out.mosaicity_deg[i] = mosaicity[img_slot];
+        }
+    }
+
+    std::vector<HKLKey> slotToHKL(nhkl);
+    for (const auto &kv: hklToSlot)
+        slotToHKL[kv.second] = kv.first;
+
+    out.merged.resize(nhkl);
+    for (int h = 0; h < nhkl; ++h) {
+        out.merged[h].h = slotToHKL[h].h;
+        out.merged[h].k = slotToHKL[h].k;
+        out.merged[h].l = slotToHKL[h].l;
+        out.merged[h].I = 0.0;
+        out.merged[h].sigma = 0.0;
+        out.merged[h].d = 0.0;
+    }
+
+    // Populate d from median of observations per HKL
+    {
+        std::vector<std::vector<double> > per_hkl_d(nhkl);
+        for (const auto &o: obs) {
+            const double d_val = static_cast<double>(o.r->d);
+            if (std::isfinite(d_val) && d_val > 0.0)
+                per_hkl_d[o.hkl_slot].push_back(d_val);
+        }
+        for (int h = 0; h < nhkl; ++h) {
+            auto &v = per_hkl_d[h];
+            if (!v.empty()) {
+                std::nth_element(v.begin(), v.begin() + static_cast<long>(v.size() / 2), v.end());
+                out.merged[h].d = v[v.size() / 2];
+            }
+        }
+    }
+
+    calc_obs(opt, g, mosaicity, rotation_crystallography, obs, corr_obs);
+    merge(nhkl, out, corr_obs);
+    stats(opt, nhkl, out, corr_obs);
 
     return out;
 }

image_analysis/scale_merge/ScaleAndMerge.h

@@ -28,11 +28,10 @@ struct ScaleMergeOptions {
     // --- Kabsch(XDS)-style partiality model ---
     // Rotation range (wedge) used in partiality calculation.
     // Set to 0 to disable partiality correction.
-    double wedge_deg = 1.0;
+    std::optional<double> wedge_deg;
 
     // --- Mosaicity (user input in degrees; internally converted to radians) ---
     bool refine_mosaicity = true;
-    bool per_image_mosaicity = true;
 
     double mosaicity_init_deg = 0.17;           // ~0.003 rad
     double mosaicity_min_deg = 1e-3;
@@ -42,12 +41,14 @@ struct ScaleMergeOptions {
     // --- Optional: regularize per-image scale k towards 1 (Kabsch-like) ---
     bool regularize_scale_to_one = true;
     double scale_regularization_sigma = 0.05;
-    bool log_scaling_residual = false;
 
     double min_partiality_for_merge = 0.2;
     bool smoothen_g = true;
     bool smoothen_mos = true;
 
+    double d_min_limit_A = 0.0;
+
+    int64_t image_cluster = 1;
 };
 
 struct MergedReflection {
@@ -80,15 +81,12 @@ struct MergeStatistics {
 
 struct ScaleMergeResult {
     std::vector<MergedReflection> merged;
-    std::vector<double> image_scale_g;
-    std::vector<int> image_ids;
-
-    // One mosaicity value per image (degrees).
-    std::vector<double> mosaicity_deg;
+    std::vector<float> image_scale_g;
+    std::vector<float> mosaicity_deg;
 
     /// Per-shell and overall merging statistics (populated after merging)
     MergeStatistics statistics;
 };
 
-ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection>& observations,
+ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<std::vector<Reflection>>& observations,
                                               const ScaleMergeOptions& opt = {});

reader/JFJochHDF5Reader.cpp

@@ -307,6 +307,8 @@ void JFJochHDF5Reader::ReadFile(const std::string& filename) {
                 dataset->experiment.IndexingAlgorithm(IndexingAlgorithmEnum::FFT);
             else if (indexing == "ffbidx")
                 dataset->experiment.IndexingAlgorithm(IndexingAlgorithmEnum::FFBIDX);
+
+            dataset->scale_factor = master_file->ReadOptVector<float>("/entry/MX/imageScaleFactor");
         }
 
         auto ring_current_A = master_file->GetOptFloat("/entry/source/current");

reader/JFJochReaderDataset.h

@@ -41,6 +41,8 @@ struct JFJochReaderDataset {
     std::vector<float> mosaicity_deg;
     std::vector<float> b_factor;
 
+    std::vector<float> scale_factor;
+
     std::vector<int64_t> max_value;
 
     std::vector<std::string> roi;

tests/JFJochReceiverProcessingTest.cpp

@@ -209,7 +209,7 @@ TEST_CASE("JFJochIntegrationTest_ZMQ_lysozyme_spot_and_index_refinement_tetragon
             .FilePrefix("lyso_test_refinement_tetragonal").JungfrauConvPhotonCnt(false).SetFileWriterFormat(
                     FileWriterFormat::NXmxVDS).OverwriteExistingFiles(true)
             .DetectorDistance_mm(75).BeamY_pxl(1136).BeamX_pxl(1090).IncidentEnergy_keV(12.4)
-            .SetUnitCell(UnitCell{.a = 36.9, .b = 78.95, .c = 78.95, .alpha = 90, .beta = 90, .gamma = 90})
+            .SetUnitCell(UnitCell{.a = 78.95, .b = 78.95, .c = 36.9, .alpha = 90, .beta = 90, .gamma = 90})
             .GeomRefinementAlgorithm(GeomRefinementAlgorithmEnum::BeamCenter)
             .SpaceGroupNumber(96);
     experiment.SampleTemperature_K(123.0).RingCurrent_mA(115);

tools/jfjoch_process.cpp

@@ -39,10 +39,10 @@ void print_usage(Logger &logger) {
     logger.Info("   -F         Use FFT indexing algorithm (default: Auto)");
     logger.Info("   -x         No least-square beam center refinement");
     logger.Info("   -d<num>    High resolution limit for spot finding (default: 1.5)");
+    logger.Info("   -D<num>    High resolution limit for scaling/merging (default: 0.0; no limit)");
     logger.Info("   -S<num>    Space group number");
     logger.Info("   -M         Scale and merge (refine mosaicity) and write scaled.hkl + image.dat");
-    logger.Info("   -L         Use log-scaling residual");
-    logger.Info("   -m<txt>    Mosaicity refinement none|fixed|image (default: image)");
+    logger.Info("   -m<txt>    Mosaicity refinement none|image (default: image)");
     logger.Info("   -A         Anomalous mode (don't merge Friedel pairs)");
 }
 
@@ -67,11 +67,11 @@ int main(int argc, char **argv) {
     bool anomalous_mode = false;
     std::optional<int> space_group_number;
 
-    bool log_residual = false;
-    enum class MosaicityRefinementMode { None, Fixed, Image };
+    enum class MosaicityRefinementMode { None, Image };
     MosaicityRefinementMode mosaicity_refinement_mode = MosaicityRefinementMode::Image;
 
-    float d_high = 1.5;
+    float d_min_spot_finding = 1.5;
+    std::optional<float> d_min_scale_merge;
 
     if (argc == 1) {
         print_usage(logger);
@@ -79,7 +79,7 @@ int main(int argc, char **argv) {
     }
 
     int opt;
-    while ((opt = getopt(argc, argv, "o:N:s:e:vR::Fxd:S:MLm:A")) != -1) {
+    while ((opt = getopt(argc, argv, "o:N:s:e:vR::Fxd:S:Mm:AD:")) != -1) {
         switch (opt) {
             case 'o':
                 output_prefix = optarg;
@@ -96,6 +96,9 @@ int main(int argc, char **argv) {
             case 'v':
                 verbose = true;
                 break;
+            case 'd':
+                d_min_spot_finding = atof(optarg);
+                break;
             case 'R':
                 rotation_indexing = true;
                 if (optarg) rotation_indexing_range = atof(optarg);
@@ -106,8 +109,9 @@ int main(int argc, char **argv) {
             case 'x':
                 refine_beam_center = false;
                 break;
-            case 'd':
-                d_high = atof(optarg);
+            case 'D':
+                d_min_scale_merge = atof(optarg);
+                logger.Info("High resolution limit for scaling/merging set to {:.2f} A", d_min_spot_finding);
                 break;
             case 'S':
                 space_group_number = atoi(optarg);
@@ -115,17 +119,12 @@ int main(int argc, char **argv) {
             case 'M':
                 run_scaling = true;
                 break;
-            case 'L':
-                log_residual = true;
-                break;
             case 'A':
                 anomalous_mode = true;
                 break;
             case 'm':
                 if (strcmp(optarg, "none") == 0)
                     mosaicity_refinement_mode = MosaicityRefinementMode::None;
-                else if (strcmp(optarg, "fixed") == 0)
-                    mosaicity_refinement_mode = MosaicityRefinementMode::Fixed;
                 else if (strcmp(optarg, "image") == 0)
                     mosaicity_refinement_mode = MosaicityRefinementMode::Image;
                 else {
@@ -206,7 +205,9 @@ int main(int argc, char **argv) {
     SpotFindingSettings spot_settings;
     spot_settings.enable = true;
     spot_settings.indexing = true;
-    spot_settings.high_resolution_limit = d_high;
+    spot_settings.high_resolution_limit = d_min_spot_finding;
+    if (d_min_scale_merge > 0)
+        spot_settings.high_resolution_limit = d_min_spot_finding;
 
     // Initialize Analysis Components
     PixelMask pixel_mask = dataset->pixel_mask;
@@ -423,20 +424,15 @@ int main(int argc, char **argv) {
         scale_opts.refine_mosaicity = true;
         scale_opts.max_num_iterations = 500;
         scale_opts.max_solver_time_s = 240.0; // generous cutoff for now
-        scale_opts.log_scaling_residual = log_residual;
         scale_opts.merge_friedel = !anomalous_mode;
+        scale_opts.d_min_limit_A = d_min_scale_merge.value_or(0.0);
 
         switch (mosaicity_refinement_mode) {
             case MosaicityRefinementMode::None:
                 scale_opts.refine_mosaicity = false;
                 break;
-            case MosaicityRefinementMode::Fixed:
-                scale_opts.refine_mosaicity = true;
-                scale_opts.per_image_mosaicity = false;
-                break;
             case MosaicityRefinementMode::Image:
                 scale_opts.refine_mosaicity = true;
-                scale_opts.per_image_mosaicity = true;
                 break;
         }
         if (space_group)
@@ -450,12 +446,10 @@ int main(int argc, char **argv) {
         double scale_time = std::chrono::duration<double>(scale_end - scale_start).count();
 
         if (scale_result) {
+            end_msg.scale_factor = scale_result->image_scale_g;
 
-// ... existing code ...
-                logger.Info("Scaling completed in {:.2f} s  ({} unique reflections, {} images)",
-                            scale_time,
-                            scale_result->merged.size(),
-                            scale_result->image_ids.size());
+                logger.Info("Scaling completed in {:.2f} s  ({} unique reflections)",
+                            scale_time, scale_result->merged.size());
 
                 // Print resolution-shell statistics table
                 {
@@ -498,13 +492,10 @@ int main(int argc, char **argv) {
                     logger.Error("Cannot open {} for writing", img_path);
                 } else {
                     img_file << "# image_id  mosaicity_deg  K\n";
-                    for (size_t i = 0; i < scale_result->image_ids.size(); ++i) {
-                        img_file << scale_result->image_ids[i] << " "
-                                << scale_result->mosaicity_deg[i] << " "
-                                << scale_result->image_scale_g[i] << "\n";
+                    for (size_t i = 0; i < scale_result->mosaicity_deg.size(); ++i) {
+                        img_file << i << " " << scale_result->mosaicity_deg[i] << " " << scale_result->image_scale_g[i] << "\n";
                     }
                     img_file.close();
-                    logger.Info("Wrote {} image records to {}", scale_result->image_ids.size(), img_path);
                 }
             }
 

viewer/JFJochViewerDatasetInfo.cpp

@@ -90,6 +90,10 @@ void JFJochViewerDatasetInfo::UpdateLabels() {
             combo_box->addItem("Mosaicity", 9);
         }
 
+        if (!dataset->scale_factor.empty()) {
+            combo_box->insertSeparator(1000);
+            combo_box->addItem("Scale factor", 10);
+        }
         for (int i = 0; i < this->dataset->roi.size(); i++) {
             std::string name = std::string("ROI ") + this->dataset->roi[i];
             combo_box->insertSeparator(1000);
@@ -167,6 +171,8 @@ void JFJochViewerDatasetInfo::UpdatePlot() {
             data = dataset->b_factor;
         } else if (val == 9) {
             data = dataset->mosaicity_deg;
+        } else if (val == 10) {
+            data = dataset->scale_factor;
         } else if (val >= 100) {
             int roi_index = (val - 100) / 4;