jfjoch_process: Enable two-pass algorithm

image_analysis/rotation_indexer/RotationIndexer.cpp

@@ -20,6 +20,8 @@ RotationIndexer::RotationIndexer(const DiffractionExperiment &x, IndexerThreadPo
 
 void RotationIndexer::RunIndexing() {
     std::unique_lock ul(m);
+    if (!axis_)
+        return;
 
     std::vector<Coord> coords;
     coords.reserve(max_spots_per_image * v_.size());

image_analysis/rotation_indexer/RotationIndexerCounter.cpp

@@ -36,7 +36,7 @@ std::pair<bool, bool> RotationIndexerCounter::Process(int64_t image) {
 
     accumulated_images++;
     ret.first = (image_stride == 1) || (image % image_stride == 0);
-    if (accumulated_images > image_to_try_indexing) {
+    if (accumulated_images >= image_to_try_indexing) {
         image_to_try_indexing += 10;
         ret.second = true;
     }

tools/jfjoch_process.cpp

@@ -51,7 +51,9 @@ void print_usage() {
     std::cout << std::endl;
 
     std::cout << "  Indexing" << std::endl;
-    std::cout << "   -R, --rotation-indexing[=num]    Rotation indexing (optional: min angular range deg)" << std::endl;
+    std::cout << "   -R, --two-pass-rotation[=num]    Two-pass offline rotation indexing (optional: number of images, default: 30)" << std::endl;
+    std::cout << "   --single-pass-rotation[=num]     Use online-like single-pass rotation indexing (optional: min angular range deg)" << std::endl;
+    std::cout << "   --redo-rotation-spots            Redo spot finding for two-pass rotation indexing" << std::endl;
     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;
@@ -81,7 +83,9 @@ enum {
     OPT_MIN_IMAGE_CC,
     OPT_SCALING_ITERATIONS,
     OPT_SCALING_HIGH_RESOLUTION,
-    OPT_SCALING_OUTPUT
+    OPT_SCALING_OUTPUT,
+    OPT_SINGLE_PASS_ROTATION,
+    OPT_REDO_ROTATION_SPOTS
 };
 
 static option long_options[] = {
@@ -91,7 +95,6 @@ static option long_options[] = {
     {"start-image", required_argument, nullptr, 's'},
     {"end-image", required_argument, nullptr, 'e'},
     {"stride", required_argument, nullptr, 't'},
-    {"rotation-indexing", optional_argument, nullptr, 'R'},
     {"indexing-algorithm", required_argument, nullptr, 'X'},
     {"unit-cell", required_argument, nullptr, 'C'},
     {"reference-mtz", required_argument, nullptr, 'z'},
@@ -103,6 +106,10 @@ static option long_options[] = {
     {"scale-merge", no_argument, nullptr, 'M'},
     {"refine", required_argument, nullptr, 'r'},
 
+    {"two-pass-rotation", optional_argument, nullptr, 'R'},
+    {"single-pass-rotation", optional_argument, nullptr, OPT_SINGLE_PASS_ROTATION},
+    {"redo-rotation-spots", no_argument, nullptr, OPT_REDO_ROTATION_SPOTS},
+
     {"spot-sigma", required_argument, nullptr, OPT_SPOT_SIGMA},
     {"spot-threshold", required_argument, nullptr, OPT_SPOT_THRESHOLD},
     {"spot-high-resolution", required_argument, nullptr, OPT_SPOT_RESOLUTION},
@@ -115,7 +122,6 @@ static option long_options[] = {
     {nullptr, 0, nullptr, 0}
 };
 
-
 void trim_in_place(std::string &t) {
     size_t b = 0;
     while (b < t.size() && std::isspace(static_cast<unsigned char>(t[b]))) b++;
@@ -173,7 +179,34 @@ std::optional<UnitCell> parse_unit_cell_arg(const char *arg) {
     if (!parse_float_strict(parts[5], uc.gamma)) return std::nullopt;
 
     return uc;
-};
+}
+
+std::vector<int> select_equally_spaced_image_ordinals(int images_to_process, int requested_images) {
+    std::vector<int> ret;
+
+    if (images_to_process <= 0 || requested_images <= 0)
+        return ret;
+
+    const int n = std::min(images_to_process, requested_images);
+    ret.reserve(n);
+
+    if (n == 1) {
+        ret.push_back(0);
+        return ret;
+    }
+
+    std::set<int> unique_ordinals;
+    for (int i = 0; i < n; i++) {
+        const auto ordinal = static_cast<int>(
+            std::llround(static_cast<double>(i) * static_cast<double>(images_to_process - 1) /
+                         static_cast<double>(n - 1))
+        );
+        unique_ordinals.insert(ordinal);
+    }
+
+    ret.assign(unique_ordinals.begin(), unique_ordinals.end());
+    return ret;
+}
 
 int main(int argc, char **argv) {
     for (int i = 0; i < argc; i++) {
@@ -197,6 +230,9 @@ int main(int argc, char **argv) {
 
     bool verbose = false;
     bool rotation_indexing = false;
+    bool two_pass_rotation = true;
+    bool reuse_rotation_spots = true;
+    int rotation_indexing_image_count = 30;
     std::optional<float> rotation_indexing_range;
     bool run_scaling = false;
     bool anomalous_mode = false;
@@ -252,8 +288,29 @@ int main(int argc, char **argv) {
             image_stride = atoi(optarg);
             break;
         case 'R':
+            if (rotation_indexing) {
+                logger.Error("Rotation indexing already enabled");
+                exit(EXIT_FAILURE);
+            }
             rotation_indexing = true;
-            if (optarg) rotation_indexing_range = atof(optarg);
+            two_pass_rotation = true;
+            if (optarg)
+                rotation_indexing_image_count = atoi(optarg);
+
+            break;
+        case OPT_SINGLE_PASS_ROTATION:
+            if (rotation_indexing) {
+                logger.Error("Rotation indexing already enabled");
+                exit(EXIT_FAILURE);
+            }
+            rotation_indexing = true;
+            two_pass_rotation = false;
+
+            if (optarg)
+                rotation_indexing_range = atof(optarg);
+            break;
+        case OPT_REDO_ROTATION_SPOTS:
+            reuse_rotation_spots = false;
             break;
         case 'X': {
             std::string alg = optarg ? optarg : "";
@@ -434,6 +491,11 @@ int main(int argc, char **argv) {
         exit(EXIT_FAILURE);
     }
 
+    if (rotation_indexing_image_count <= 0) {
+        logger.Error("Invalid number of rotation indexing images: {}", rotation_indexing_image_count);
+        exit(EXIT_FAILURE);
+    }
+
     logger.Info("Loaded dataset from {}", input_file);
 
     std::vector<MergedReflection> reference_data;
@@ -451,6 +513,11 @@ int main(int argc, char **argv) {
         exit(EXIT_FAILURE);
     }
 
+    if (image_stride == 0) {
+        logger.Error("Image stride cannot be zero");
+        exit(EXIT_FAILURE);
+    }
+
     int images_to_process = (end_image - start_image) / image_stride;
 
     if (images_to_process <= 0) {
@@ -568,6 +635,59 @@ int main(int argc, char **argv) {
     if (!reference_data.empty())
         indexer.ReferenceIntensities(reference_data);
 
+    if (rotation_indexing && two_pass_rotation) {
+        const auto selected_ordinals = select_equally_spaced_image_ordinals(
+            images_to_process, rotation_indexing_image_count);
+
+        logger.Info("Running first pass rotation indexing using {} equally-spaced images{}",
+                    selected_ordinals.size(),
+                    reuse_rotation_spots ? " and stored spots" : "");
+
+        for (const auto image_ordinal: selected_ordinals) {
+            const int image_idx = start_image + image_ordinal * image_stride;
+
+            DataMessage msg{};
+            msg.number = image_ordinal;
+            msg.original_number = image_idx;
+
+            try {
+                if (reuse_rotation_spots) {
+                    msg.spots = reader.ReadSpots(image_idx);
+                } else {
+                    auto img = reader.GetRawImage(image_idx);
+                    if (!img)
+                        continue;
+
+                    MXAnalysisWithoutFPGA analysis(experiment, mapping, pixel_mask, indexer);
+                    AzimuthalIntegrationProfile profile(mapping);
+
+                    auto first_pass_spot_settings = spot_settings;
+                    first_pass_spot_settings.indexing = false;
+                    first_pass_spot_settings.quick_integration = false;
+
+                    msg.image = img->image;
+                    if (dataset->efficiency.size() > image_idx)
+                        msg.image_collection_efficiency = dataset->efficiency[image_idx];
+
+                    analysis.Analyze(msg, profile, first_pass_spot_settings);
+                }
+
+                indexer.AddImageToRotationIndexer(msg);
+            } catch (const std::exception &e) {
+                logger.Warning("Failed to add image {} to first-pass rotation indexing: {}", image_idx, e.what());
+            }
+        }
+
+        const auto rotation_indexing_result = indexer.FinalizeRotationIndexing();
+
+        if (!rotation_indexing_result.has_value()) {
+            logger.Error("Two-pass rotation indexing failed - stopping without processing images");
+            return EXIT_FAILURE;
+        }
+
+        logger.Info("Two-pass rotation indexing found lattice");
+    }
+
     std::atomic<int> finished_count = 0;
 
     auto worker = [&](int thread_id) {
@@ -577,8 +697,8 @@ int main(int argc, char **argv) {
         AzimuthalIntegrationProfile profile(mapping);
 
         while (true) {
-            int current_idx_offset = processed_count.fetch_add(image_stride);
-            int image_idx = start_image + current_idx_offset;
+            int image_ordinal = processed_count.fetch_add(1);
+            int image_idx = start_image + image_ordinal * image_stride;
 
             if (image_idx >= end_image) break;
 
@@ -596,7 +716,7 @@ int main(int argc, char **argv) {
             DataMessage msg{};
 
             msg.image = img->image;
-            msg.number = current_idx_offset;
+            msg.number = image_ordinal;
             msg.original_number = image_idx;
             if (dataset->efficiency.size() > image_idx) msg.image_collection_efficiency = dataset->efficiency[image_idx];
 
@@ -627,7 +747,7 @@ int main(int argc, char **argv) {
             finished_count.fetch_add(1);
 
             // Progress log
-            if (current_idx_offset > 0 && current_idx_offset % 100 == 0) {
+            if (image_ordinal > 0 && image_ordinal % 100 == 0) {
                 std::optional<float> indexing_rate = plots.GetIndexingRate();
 
                 const auto now = std::chrono::steady_clock::now();
@@ -637,12 +757,12 @@ int main(int argc, char **argv) {
 
                 if (indexing_rate.has_value()) {
                     logger.Info("Processed {} / {} images ({:.2f} Hz, indexing rate {:.1f}%)",
-                                current_idx_offset, images_to_process,
+                                image_ordinal, images_to_process,
                                 frame_rate_hz,
                                 indexing_rate.value() * 100.0f);
                 } else {
                     logger.Info("Processed {} / {} images ({:.2f} Hz, indexing rate N/A)",
-                                current_idx_offset, images_to_process,
+                                image_ordinal, images_to_process,
                                 frame_rate_hz);
                 }
             }