Merge: Some clean-up of CC merge mask + add unit cell handling in IndexAndRefine

image_analysis/IndexAndRefine.cpp

@@ -20,7 +20,8 @@ IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool
         rotation_indexer = std::make_unique<RotationIndexer>(x, *indexer);
     reflections.resize(x.GetImageNum());
     mosaicity.resize(x.GetImageNum(), NAN);
-    merge_mask.resize(x.GetImageNum(), 0);
+    scale_cc.resize(x.GetImageNum(), 0);
+    unit_cells.resize(x.GetImageNum());
 }
 
 IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetryRotation(DataMessage &msg) {
@@ -256,12 +257,13 @@ void IndexAndRefine::ProcessImage(DataMessage &msg,
     if (experiment.GetIndexingSettings().GetGeomRefinementAlgorithm() != GeomRefinementAlgorithmEnum::None)
         RefineGeometryIfNeeded(msg, outcome);
 
-    if (!outcome.lattice_candidate)
+    if (!outcome.lattice_candidate.has_value())
         return;
 
     if (!AnalyzeIndexing(msg, outcome.experiment, *outcome.lattice_candidate))
         return;
 
+    unit_cells[msg.number] = outcome.lattice_candidate->GetUnitCell();
     msg.lattice_type = outcome.symmetry;
 
     if (spot_finding_settings.quick_integration)
@@ -289,10 +291,7 @@ void IndexAndRefine::ScaleImage(DataMessage &msg) {
     msg.image_scale_factor = scaling_result.G;
     msg.image_scale_mosaicity = scaling_result.mos;
     msg.image_scale_cc = scaling_result.cc;
-    if (scaling_result.cc >= experiment.GetScalingSettings().GetMinCCForImage())
-        merge_mask[msg.number] = 1;
-    else
-        merge_mask[msg.number] = 0;
+    scale_cc[msg.number] = scaling_result.cc;
 
     auto scaling_end_time = std::chrono::steady_clock::now();
     msg.image_scale_time_s = std::chrono::duration<float>(scaling_end_time - scaling_start_time).count();
@@ -301,15 +300,19 @@ void IndexAndRefine::ScaleImage(DataMessage &msg) {
 ScalingResult IndexAndRefine::ScaleAllImages(const std::vector<MergedReflection> &reference, size_t nthreads) {
     ScaleOnTheFly scaling(experiment, reference);
     auto result = scaling.Scale(reflections, mosaicity, nthreads);
-    merge_mask = CalcMergeMask(experiment, result);
+    scale_cc = result.image_cc;
     return result;
 }
 
 MergeResult IndexAndRefine::Merge(bool calc_statistics, bool apply_cc_limit) const {
     MergeResult out;
-    const auto& mask = apply_cc_limit ? merge_mask : std::vector<uint8_t>{};
-    out.merged = MergeAll(experiment, reflections, mask);
+
+    std::vector<uint8_t> merge_mask(reflections.size(), 1);
+    if (apply_cc_limit)
+        CalcMergeMask(experiment, scale_cc, merge_mask);
+
+    out.merged = MergeAll(experiment, reflections, merge_mask);
     if (calc_statistics)
-        out.statistics = MergeStats(experiment, out.merged, reflections, mask);
+        out.statistics = MergeStats(experiment, out.merged, reflections, merge_mask);
     return out;
 }

image_analysis/IndexAndRefine.h

@@ -48,7 +48,8 @@ class IndexAndRefine {
     mutable std::mutex reflections_mutex;
     std::vector<std::vector<Reflection>> reflections;
     std::vector<float> mosaicity;
-    std::vector<uint8_t> merge_mask;
+    std::vector<float> scale_cc;
+    std::vector<std::optional<UnitCell> > unit_cells;
 
     IndexingOutcome DetermineLatticeAndSymmetryRotation(DataMessage &msg);
     IndexingOutcome DetermineLatticeAndSymmetry(DataMessage &msg);

image_analysis/scale_merge/Merge.cpp

@@ -16,15 +16,19 @@
 // Need to use consistent lattice to calc resolution (which would also allow to calculate completeness)
 // But this is not possible for current still workflow (while fine for rotation workflow)
 
-std::vector<uint8_t> CalcMergeMask(const DiffractionExperiment &x,
-                                   const ScalingResult &result) {
-    std::vector<uint8_t> ret(result.image_cc.size(), 0);
+void CalcMergeMask(const DiffractionExperiment &x,
+                   const std::vector<float> &scale_cc,
+                   std::vector<uint8_t> &result_mask) {
+
+    if (scale_cc.size() != result_mask.size())
+        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Scale CC size mismatch");
 
     auto min_cc = x.GetScalingSettings().GetMinCCForImage();
-    for (int i = 0 ; i < result.image_cc.size(); ++i)
-        ret[i] = (result.image_cc[i] >= min_cc) ? 1 : 0;
+    if (min_cc <= 0.0)
+        return;
 
-    return ret;
+    for (int i = 0; i < result_mask.size(); ++i)
+        result_mask[i] = result_mask[i] && (scale_cc[i] < min_cc);
 }
 
 std::vector<MergedReflection> MergeAll(const DiffractionExperiment &x,

image_analysis/scale_merge/Merge.h

@@ -32,8 +32,9 @@ struct MergeResult {
     MergeStatistics statistics;
 };
 
-std::vector<uint8_t> CalcMergeMask(const DiffractionExperiment &x,
-                                   const ScalingResult &result);
+void CalcMergeMask(const DiffractionExperiment &x,
+                   const std::vector<float> &scale_cc,
+                   std::vector<uint8_t> &result_mask);
 
 std::vector<MergedReflection> MergeAll(const DiffractionExperiment &x,
                                        const std::vector<std::vector<Reflection> > &reflections,

tools/jfjoch_scale.cpp

@@ -271,9 +271,11 @@ int main(int argc, char **argv) {
     double scale_time = std::chrono::duration<double>(scale_end - scale_start).count();
     logger.Info("Scaling completed in {:.2f} s", scale_time);
 
+    std::vector<uint8_t> merge_mask(images_to_process, 1);
+
     auto merge_start = std::chrono::steady_clock::now();
 
-    auto merge_mask = CalcMergeMask(experiment, scale_result);
+    CalcMergeMask(experiment, scale_result.image_cc, merge_mask);
     auto merge_result = MergeAll(experiment, reflections, merge_mask);
     auto merge_stats = MergeStats(experiment, merge_result, reflections, merge_mask);