From 6c8b2d4ddd1eae75b4ae83b6784ff0b3bbe07074 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 29 May 2026 17:15:13 +0200 Subject: [PATCH] jfjoch_process: Enable two-pass algorithm --- .../rotation_indexer/RotationIndexer.cpp | 2 + .../RotationIndexerCounter.cpp | 2 +- tools/jfjoch_process.cpp | 144 ++++++++++++++++-- 3 files changed, 135 insertions(+), 13 deletions(-) diff --git a/image_analysis/rotation_indexer/RotationIndexer.cpp b/image_analysis/rotation_indexer/RotationIndexer.cpp index b8fe2886..3f588a4f 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.cpp +++ b/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 coords; coords.reserve(max_spots_per_image * v_.size()); diff --git a/image_analysis/rotation_indexer/RotationIndexerCounter.cpp b/image_analysis/rotation_indexer/RotationIndexerCounter.cpp index 21bce4bc..79bc7efe 100644 --- a/image_analysis/rotation_indexer/RotationIndexerCounter.cpp +++ b/image_analysis/rotation_indexer/RotationIndexerCounter.cpp @@ -36,7 +36,7 @@ std::pair 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; } diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 4c9dec13..f3b35425 100644 --- a/tools/jfjoch_process.cpp +++ b/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 Indexing algorithm (FFBIDX|FFT|FFTW|Auto|None)" << std::endl; std::cout << " -S, --space-group Space group number - used for both indexing and scaling" << std::endl; std::cout << " -C, --unit-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(t[b]))) b++; @@ -173,7 +179,34 @@ std::optional parse_unit_cell_arg(const char *arg) { if (!parse_float_strict(parts[5], uc.gamma)) return std::nullopt; return uc; -}; +} + +std::vector select_equally_spaced_image_ordinals(int images_to_process, int requested_images) { + std::vector 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 unique_ordinals; + for (int i = 0; i < n; i++) { + const auto ordinal = static_cast( + std::llround(static_cast(i) * static_cast(images_to_process - 1) / + static_cast(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 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 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 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 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); } }