diff --git a/CMakeLists.txt b/CMakeLists.txt index af868b1e..0148f365 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -262,6 +262,7 @@ IF (JFJOCH_VIEWER_ONLY) ADD_SUBDIRECTORY(image_analysis) ADD_SUBDIRECTORY(broker) ADD_SUBDIRECTORY(reader) + ADD_SUBDIRECTORY(process) ADD_SUBDIRECTORY(viewer) ADD_SUBDIRECTORY(tools) # builds only the portable analysis tools (process/scale/azint/extract_hkl) ELSE() @@ -286,6 +287,7 @@ ELSE() ADD_SUBDIRECTORY(acquisition_device) ADD_SUBDIRECTORY(receiver) ADD_SUBDIRECTORY(image_analysis) + ADD_SUBDIRECTORY(process) ADD_SUBDIRECTORY(tests) ADD_SUBDIRECTORY(tools) ENDIF() diff --git a/process/CMakeLists.txt b/process/CMakeLists.txt new file mode 100644 index 00000000..4ebcbac7 --- /dev/null +++ b/process/CMakeLists.txt @@ -0,0 +1,9 @@ +# SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +# SPDX-License-Identifier: GPL-3.0-only + +ADD_LIBRARY(JFJochProcess STATIC + JFJochProcess.cpp + JFJochProcess.h +) + +TARGET_LINK_LIBRARIES(JFJochProcess JFJochReader JFJochImageAnalysis JFJochWriter) diff --git a/process/JFJochProcess.cpp b/process/JFJochProcess.cpp new file mode 100644 index 00000000..6167a449 --- /dev/null +++ b/process/JFJochProcess.cpp @@ -0,0 +1,418 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "JFJochProcess.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../reader/JFJochHDF5Reader.h" +#include "../common/Logger.h" +#include "../common/AzimuthalIntegrationMapping.h" +#include "../common/AzimuthalIntegrationProfile.h" +#include "../common/CUDAWrapper.h" +#include "../common/time_utc.h" +#include "../writer/FileWriter.h" +#include "../image_analysis/MXAnalysisWithoutFPGA.h" +#include "../image_analysis/IndexAndRefine.h" +#include "../image_analysis/indexing/IndexerThreadPool.h" +#include "../image_analysis/azint/AzIntEngineCPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorCPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h" +#include "../image_analysis/scale_merge/Merge.h" +#include "../image_analysis/scale_merge/SearchSpaceGroup.h" +#include "../image_analysis/WriteReflections.h" + +namespace { + // Pick up to requested_images ordinals spread evenly across [0, images_to_process) for the + // first pass of two-pass rotation indexing. + 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); + if (n == 1) { + ret.push_back(0); + return ret; + } + + std::set unique_ordinals; + for (int i = 0; i < n; i++) + unique_ordinals.insert(static_cast( + std::llround(static_cast(i) * static_cast(images_to_process - 1) / + static_cast(n - 1)))); + + ret.assign(unique_ordinals.begin(), unique_ordinals.end()); + return ret; + } +} + +JFJochProcess::JFJochProcess(JFJochHDF5Reader &reader, DiffractionExperiment experiment, + PixelMask pixel_mask, ProcessConfig config) + : reader_(reader), experiment_(std::move(experiment)), + pixel_mask_(std::move(pixel_mask)), config_(std::move(config)) {} + +ProcessResult JFJochProcess::Run(JFJochProcessObserver *observer) { + Logger logger("JFJochProcess"); + ProcessResult result; + + const auto dataset = reader_.GetDataset(); + if (!dataset) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "No experiment dataset found in the input file"); + + if (config_.stride <= 0) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Image stride must be positive"); + + const auto total_images_in_file = static_cast(reader_.GetNumberOfImages()); + int end_image = config_.end_image; + if (end_image < 0 || end_image > total_images_in_file) + end_image = total_images_in_file; + const int start_image = config_.start_image; + const int images_to_process = (end_image - start_image) / config_.stride; + if (images_to_process <= 0) { + logger.Warning("No images to process (start {}, end {}, stride {}, total {})", + start_image, end_image, config_.stride, total_images_in_file); + return result; + } + + const bool full = (config_.mode == ProcessMode::FullAnalysis); + const bool write_files = !config_.output_prefix.empty(); + + // Output/runtime invariants. Algorithm settings (indexing, scaling, integration, polarization, + // space group, unit cell, ...) are configured on experiment_ by the caller. + experiment_.BitDepthImage(32).PixelSigned(true); + experiment_.Mode(DetectorMode::Standard); + experiment_.OverwriteExistingFiles(true); + experiment_.FilePrefix(config_.output_prefix.empty() ? "output" : config_.output_prefix); + experiment_.SetFileWriterFormat(FileWriterFormat::NXmxLegacy); + experiment_.ImagesPerTrigger(images_to_process); + experiment_.NumTriggers(1); + if (full) + experiment_.Compression(CompressionAlgorithm::BSHUF_LZ4); + + AzimuthalIntegrationMapping mapping(experiment_, pixel_mask_); + + JFJochReceiverPlots plots; + plots.Setup(experiment_, mapping); + + // Output file (NXmxIntegrated master that links back to the original images). + StartMessage start_message; + experiment_.FillMessage(start_message); + start_message.arm_date = dataset->arm_date; + start_message.az_int_bin_to_q = mapping.GetBinToQ(); + start_message.az_int_bin_to_two_theta = mapping.GetBinToTwoTheta(); + start_message.az_int_q_bin_count = mapping.GetQBinCount(); + start_message.az_int_phi_bin_count = mapping.GetAzimuthalBinCount(); + if (mapping.GetAzimuthalBinCount() > 1) + start_message.az_int_bin_to_phi = mapping.GetBinToPhi(); + start_message.pixel_mask["default"] = pixel_mask_.GetMask(experiment_); + if (full) { + start_message.rois = experiment_.ROI().ExportMetadata(); + if (!experiment_.ROI().empty()) + start_message.roi_map = experiment_.ExportROIMap(); + start_message.max_spot_count = experiment_.GetMaxSpotCount(); + } + start_message.master_suffix = "process"; + start_message.file_format = FileWriterFormat::NXmxIntegrated; + start_message.write_master_file = true; + start_message.write_images = false; + start_message.hdf5_source_data = reader_.GetHDF5DataSource(start_image, images_to_process); + + std::unique_ptr writer; + if (write_files) + writer = std::make_unique(start_message); + + logger.Info("Processing {} images (range {}-{}, stride {}) using {} threads [{}]", + images_to_process, start_image, end_image, config_.stride, config_.nthreads, + full ? "full analysis" : "azimuthal integration"); + if (observer) + observer->OnPhase(full ? "Full analysis" : "Azimuthal integration"); + + // Full-analysis shared engines. + std::unique_ptr indexer_pool; + std::unique_ptr indexer; + if (full) { + indexer_pool = std::make_unique(experiment_.GetIndexingSettings()); + indexer = std::make_unique(experiment_, indexer_pool.get()); + if (!config_.reference_data.empty()) + indexer->ReferenceIntensities(config_.reference_data); + } + + const auto start_time = std::chrono::steady_clock::now(); + + // First pass of two-pass rotation indexing (full analysis only). + if (full && config_.forced_rotation_lattice.has_value()) { + indexer->ForceRotationIndexerLattice(*config_.forced_rotation_lattice); + logger.Info("Rotation indexer lattice forced externally - skipping first pass"); + } else if (full && config_.rotation_indexing && config_.two_pass_rotation) { + if (observer) + observer->OnPhase("Rotation indexing (first pass)"); + const auto selected = select_equally_spaced_image_ordinals(images_to_process, + config_.rotation_indexing_image_count); + logger.Info("First-pass rotation indexing using {} images{}", selected.size(), + config_.reuse_rotation_spots ? " and stored spots" : ""); + + for (const int ordinal: selected) { + if (cancelled_) break; + const int image_idx = start_image + ordinal * config_.stride; + DataMessage msg{}; + msg.number = image_idx; + msg.original_number = image_idx; + try { + if (config_.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 = config_.spot_finding; + first_pass.indexing = false; + first_pass.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); + } + indexer->AddImageToRotationIndexer(msg); + } catch (const std::exception &e) { + logger.Warning("First-pass rotation indexing failed for image {}: {}", image_idx, e.what()); + } + } + + if (!cancelled_ && !indexer->FinalizeRotationIndexing().has_value()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Two-pass rotation indexing failed"); + if (!cancelled_) + logger.Info("Two-pass rotation indexing found lattice"); + } + + // Main per-image loop, spread over N worker threads pulling from a shared counter. HDF5 reads + // are serialized by the global hdf5_mutex; the analysis runs in parallel. + std::atomic next_ordinal = 0; + std::atomic finished_count = 0; + std::atomic total_uncompressed_bytes = 0; + + auto azint_worker = [&]() { + std::vector decompression_buffer; + ImagePreprocessorCPU preprocessor(experiment_, pixel_mask_); + ImagePreprocessorBuffer buffer(experiment_.GetPixelsNum()); + AzIntEngineCPU azint(mapping); + AzimuthalIntegrationProfile profile(mapping); + + while (!cancelled_) { + const int ordinal = next_ordinal.fetch_add(1); + const int image_idx = start_image + ordinal * config_.stride; + if (image_idx >= end_image) break; + + std::shared_ptr img; + try { + img = reader_.GetRawImage(image_idx); + } catch (const std::exception &e) { + logger.Error("Failed to load image {}: {}", image_idx, e.what()); + continue; + } + if (!img) continue; + + DataMessage msg{}; + msg.image = img->image; + msg.number = ordinal; + msg.original_number = image_idx; + if (dataset->efficiency.size() > image_idx) + msg.image_collection_efficiency = dataset->efficiency[image_idx]; + total_uncompressed_bytes += msg.image.GetUncompressedSize(); + + const auto t0 = std::chrono::steady_clock::now(); + try { + const uint8_t *image_ptr = msg.image.GetUncompressedPtr(decompression_buffer); + preprocessor.Analyze(buffer, image_ptr, msg.image.GetMode()); + azint.Run(buffer, profile); + } catch (const std::exception &e) { + logger.Error("Error integrating image {}: {}", image_idx, e.what()); + continue; + } + msg.azint_time_s = std::chrono::duration(std::chrono::steady_clock::now() - t0).count(); + msg.processing_time_s = msg.azint_time_s; + msg.az_int_profile = profile.GetResult(); + msg.az_int_profile_count = profile.GetPixelCount(); + msg.az_int_profile_std = profile.GetStd(); + msg.bkg_estimate = profile.GetBkgEstimate(mapping.Settings()); + msg.run_number = experiment_.GetRunNumber(); + msg.run_name = experiment_.GetRunName(); + + plots.Add(msg, profile); + if (writer) writer->Write(msg); + if (observer) observer->OnImageProcessed(msg); + const int done = finished_count.fetch_add(1) + 1; + if (observer) observer->OnProgress(done, images_to_process); + } + }; + + auto full_worker = [&]() { + pin_gpu(); // round-robin per worker thread; must precede engine construction + MXAnalysisWithoutFPGA analysis(experiment_, mapping, pixel_mask_, *indexer); + AzimuthalIntegrationProfile profile(mapping); + + while (!cancelled_) { + const int ordinal = next_ordinal.fetch_add(1); + const int image_idx = start_image + ordinal * config_.stride; + if (image_idx >= end_image) break; + + std::shared_ptr img; + try { + img = reader_.GetRawImage(image_idx); + } catch (const std::exception &e) { + logger.Error("Failed to load image {}: {}", image_idx, e.what()); + continue; + } + if (!img) continue; + + DataMessage msg{}; + msg.image = img->image; + msg.number = ordinal; + msg.original_number = image_idx; + if (dataset->efficiency.size() > image_idx) + msg.image_collection_efficiency = dataset->efficiency[image_idx]; + total_uncompressed_bytes += msg.image.GetUncompressedSize(); + + const auto t0 = std::chrono::steady_clock::now(); + try { + analysis.Analyze(msg, profile, config_.spot_finding); + } catch (const std::exception &e) { + logger.Error("Error analyzing image {}: {}", image_idx, e.what()); + continue; + } + msg.processing_time_s = std::chrono::duration(std::chrono::steady_clock::now() - t0).count(); + msg.run_number = experiment_.GetRunNumber(); + msg.run_name = experiment_.GetRunName(); + + plots.Add(msg, profile); + if (writer) writer->Write(msg); + if (observer) observer->OnImageProcessed(msg); + const int done = finished_count.fetch_add(1) + 1; + if (observer) observer->OnProgress(done, images_to_process); + } + }; + + if (observer) + observer->OnPhase("Processing images"); + + std::function worker = full ? std::function(full_worker) + : std::function(azint_worker); + std::vector > futures; + futures.reserve(config_.nthreads); + for (int i = 0; i < config_.nthreads; ++i) + futures.push_back(std::async(std::launch::async, worker)); + for (auto &f: futures) + f.get(); + + const auto end_time = std::chrono::steady_clock::now(); + result.cancelled = cancelled_; + result.images_processed = finished_count.load(); + result.processing_time_s = std::chrono::duration(end_time - start_time).count(); + if (result.processing_time_s > 0.0) { + result.frame_rate_hz = static_cast(result.images_processed) / result.processing_time_s; + result.throughput_MBs = static_cast(total_uncompressed_bytes) / (result.processing_time_s * 1e6); + } + result.mean_processing_time = plots.GetMeanProcessingTime(); + result.indexing_rate = plots.GetIndexingRate(); + + // End message (also written to the file). + EndMessage end_msg; + end_msg.max_image_number = result.images_processed; + end_msg.images_collected_count = result.images_processed; + end_msg.images_sent_to_write_count = result.images_processed; + end_msg.end_date = time_UTC(std::chrono::system_clock::now()); + end_msg.run_number = experiment_.GetRunNumber(); + end_msg.run_name = experiment_.GetRunName(); + end_msg.bkg_estimate = plots.GetBkgEstimate(); + end_msg.az_int_result["dataset"] = plots.GetAzIntProfile(); + end_msg.indexing_rate = result.indexing_rate; + + if (full && !cancelled_) { + if (const auto rot = indexer->FinalizeRotationIndexing(); rot.has_value()) { + end_msg.rotation_lattice = rot->lattice; + end_msg.rotation_lattice_type = LatticeMessage{ + .centering = rot->search_result.centering, + .niggli_class = rot->search_result.niggli_class, + .crystal_system = rot->search_result.system + }; + result.rotation_lattice_found = true; + } + result.consensus_cell = indexer->GetConsensusUnitCell(); + end_msg.unit_cell = result.consensus_cell; + } + + // Scaling and merging (full analysis only). + if (full && !cancelled_ && result.indexing_rate.has_value() && result.indexing_rate > 0 + && (config_.run_scaling || !config_.reference_data.empty())) { + if (observer) + observer->OnPhase("Scaling and merging"); + + const bool pixel_refine_path = + experiment_.GetIndexingSettings().GetGeomRefinementAlgorithm() == GeomRefinementAlgorithmEnum::PixelRefine; + + // ScaleOnTheFly is only for the classical, no-reference path; with a reference (or + // PixelRefine) each image is already scaled, so we merge directly. + if (config_.reference_data.empty() && !pixel_refine_path) { + logger.Info("Running scaling ..."); + ScalingResult scale_result(0); + for (int i = 0; i < config_.scaling_iter; i++) { + auto merge_result = MergeAll(experiment_, indexer->GetIntegrationOutcome(), false); + scale_result = indexer->ScaleAllImages(merge_result); + } + } + + MergeOnTheFly merge_engine(experiment_); + if (result.consensus_cell.has_value()) + merge_engine.ReferenceCell(*result.consensus_cell); + merge_engine.RefineErrorModel(indexer->GetIntegrationOutcome()); + for (const auto &outcome: indexer->GetIntegrationOutcome()) + merge_engine.AddImage(outcome); + + auto merged_reflections = merge_engine.ExportReflections(); + auto merged_statistics = merge_engine.MergeStats(merged_reflections, indexer->GetIntegrationOutcome(), + config_.reference_data); + logger.Info("Merge complete ({} unique reflections)", merged_reflections.size()); + + std::ostringstream stats_text; + if (!experiment_.GetGemmiSpaceGroup().has_value()) { + SearchSpaceGroupOptions sg_opts; + sg_opts.crystal_system.reset(); + sg_opts.centering = '\0'; + sg_opts.merge_friedel = experiment_.GetScalingSettings().GetMergeFriedel(); + sg_opts.d_min_limit_A = experiment_.GetScalingSettings().GetHighResolutionLimit_A().value_or(0.0); + sg_opts.min_operator_cc = 0.80; + sg_opts.min_pairs_per_operator = 20; + sg_opts.min_total_compared = 100; + sg_opts.test_systematic_absences = true; + stats_text << SearchSpaceGroupResultToText(SearchSpaceGroup(merged_reflections, sg_opts)) << "\n\n"; + } + stats_text << merged_statistics; + result.merge_statistics_text = stats_text.str(); + + if (result.consensus_cell && write_files) + WriteReflections(merged_reflections, *result.consensus_cell, experiment_, config_.output_prefix); + } + + if (writer) { + writer->WriteHDF5(end_msg); + writer->Finalize(); + result.written_master_path = config_.output_prefix + "_process.h5"; + } + + if (observer) + observer->OnPhase(cancelled_ ? "Cancelled" : "Done"); + logger.Info("{} {} images in {:.2f} s ({:.2f} Hz)", cancelled_ ? "Cancelled after" : "Processed", + result.images_processed, result.processing_time_s, result.frame_rate_hz); + return result; +} diff --git a/process/JFJochProcess.h b/process/JFJochProcess.h new file mode 100644 index 00000000..ecaf3c58 --- /dev/null +++ b/process/JFJochProcess.h @@ -0,0 +1,100 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include +#include + +#include "../common/DiffractionExperiment.h" +#include "../common/PixelMask.h" +#include "../common/CrystalLattice.h" +#include "../common/Reflection.h" // MergedReflection +#include "../common/JFJochMessages.h" // DataMessage +#include "../common/JFJochReceiverPlots.h" // MeanProcessingTime +#include "../image_analysis/spot_finding/SpotFindingSettings.h" + +class JFJochHDF5Reader; + +// Offline reprocessing of a stored Jungfraujoch HDF5 dataset, shared by jfjoch_process, +// jfjoch_azint and (later) the viewer. The full processing workflow lives here, not in the CLIs: +// setup, an optional two-pass rotation-indexing pre-pass, a parallel per-image loop (std::thread), +// an optional scaling/merging post-pass, and the _process.h5 output. The detector geometry and all +// algorithm settings are configured on the DiffractionExperiment by the caller; ProcessConfig only +// carries run control. Cancellable from any thread (e.g. SIGINT or a GUI button) via Cancel(). +enum class ProcessMode { + AzimuthalIntegration, // preprocess + azimuthal integration only (jfjoch_azint) + FullAnalysis // spot finding + indexing + refinement + integration (jfjoch_process) +}; + +struct ProcessConfig { + ProcessMode mode = ProcessMode::FullAnalysis; + + int start_image = 0; + int end_image = -1; // -1 => to the end of the dataset + int stride = 1; + int nthreads = 1; + + // Output prefix for the _process.h5 (and scaled reflections). Empty => process without writing. + std::string output_prefix; + + SpotFindingSettings spot_finding; // FullAnalysis spot finding + + // Rotation indexing (FullAnalysis) + bool rotation_indexing = false; + bool two_pass_rotation = true; + bool reuse_rotation_spots = true; + int rotation_indexing_image_count = 30; + std::optional forced_rotation_lattice; + + // Scaling / merging (FullAnalysis). reference_data (from a reference MTZ) also drives the + // per-image live scaling path when present. + bool run_scaling = false; + int64_t scaling_iter = 3; + std::vector reference_data; +}; + +struct ProcessResult { + bool cancelled = false; + uint64_t images_processed = 0; + double processing_time_s = 0.0; + double frame_rate_hz = 0.0; + double throughput_MBs = 0.0; + std::optional indexing_rate; + std::optional consensus_cell; + bool rotation_lattice_found = false; + MeanProcessingTime mean_processing_time{}; + std::optional written_master_path; + std::string merge_statistics_text; // populated when scaling/merging ran +}; + +// Callbacks for progress and live results. Methods may be called from worker threads, so an +// implementation must be thread-safe. The default no-ops suit the CLIs. +class JFJochProcessObserver { +public: + virtual ~JFJochProcessObserver() = default; + virtual void OnPhase(const std::string &phase) {} + virtual void OnProgress(uint64_t done, uint64_t total) {} + virtual void OnImageProcessed(const DataMessage &msg) {} +}; + +class JFJochProcess { + JFJochHDF5Reader &reader_; + DiffractionExperiment experiment_; + PixelMask pixel_mask_; + ProcessConfig config_; + std::atomic cancelled_{false}; + +public: + JFJochProcess(JFJochHDF5Reader &reader, DiffractionExperiment experiment, + PixelMask pixel_mask, ProcessConfig config); + + // Runs the configured workflow to completion or until Cancel(). Throws on setup failure. + ProcessResult Run(JFJochProcessObserver *observer = nullptr); + + // Request cancellation; safe to call from any thread (the worker loop checks between images). + void Cancel() { cancelled_ = true; } + bool IsCancelled() const { return cancelled_; } +}; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index eb32e313..7707b3c8 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -45,6 +45,7 @@ ADD_EXECUTABLE(jfjoch_test ModuleSummationTest.cpp ZMQMetadataSocketTest.cpp JFJochReaderTest.cpp + JFJochProcessTest.cpp MovingAverageTest.cpp ImageMetadataTest.cpp JFJochReceiverLiteTest.cpp @@ -77,7 +78,7 @@ ADD_EXECUTABLE(jfjoch_test ) target_link_libraries(jfjoch_test Catch2WithMain JFJochBroker JFJochReceiver JFJochReader JFJochStreamWriter - JFJochImageAnalysis JFJochCommon JFJochHLSSimulation JFJochPreview + JFJochProcess JFJochImageAnalysis JFJochCommon JFJochHLSSimulation JFJochPreview jfjoch_xds_plugin) target_include_directories(jfjoch_test PRIVATE .) diff --git a/tests/JFJochProcessTest.cpp b/tests/JFJochProcessTest.cpp new file mode 100644 index 00000000..5e744016 --- /dev/null +++ b/tests/JFJochProcessTest.cpp @@ -0,0 +1,133 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include + +#include "../common/DiffractionExperiment.h" +#include "../common/ScanResultGenerator.h" +#include "../writer/FileWriter.h" +#include "../reader/JFJochHDF5Reader.h" +#include "../process/JFJochProcess.h" + +namespace { + // Write a small VDS dataset of `n` flat images and return nothing (prefix_master.h5 + + // prefix_data_000001.h5 land in the test working directory). + void WriteTestDataset(const std::string &prefix, int n) { + RegisterHDF5Filter(); + + DiffractionExperiment x(DetJF(1)); + x.FilePrefix(prefix).ImagesPerTrigger(n).OverwriteExistingFiles(true); + x.BitDepthImage(16).ImagesPerFile(n).SetFileWriterFormat(FileWriterFormat::NXmxVDS).PixelSigned(true); + x.Compression(CompressionAlgorithm::NO_COMPRESSION); + x.BeamX_pxl(512).BeamY_pxl(256).DetectorDistance_mm(150).IncidentEnergy_keV(WVL_1A_IN_KEV) + .FrameTime(std::chrono::microseconds(500), std::chrono::microseconds(10)); + + std::vector image(x.GetPixelsNum(), 5); + StartMessage start_message; + x.FillMessage(start_message); + FileWriter file_set(start_message); + ScanResultGenerator generator(x); + for (int i = 0; i < n; i++) { + DataMessage message{}; + message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum()); + message.number = i; + REQUIRE_NOTHROW(file_set.WriteHDF5(message)); + generator.Add(message); + } + EndMessage end_message; + end_message.max_image_number = n; + generator.FillEndMessage(end_message); + file_set.WriteHDF5(end_message); + file_set.Finalize(); + } +} + +TEST_CASE("JFJochProcess_AzInt", "[HDF5][Full]") { + WriteTestDataset("process_azint_in", 8); + + JFJochHDF5Reader reader; + REQUIRE_NOTHROW(reader.ReadFile("process_azint_in_master.h5")); + auto dataset = reader.GetDataset(); + REQUIRE(dataset); + + ProcessConfig config; + config.mode = ProcessMode::AzimuthalIntegration; + config.nthreads = 2; + config.output_prefix = "process_azint_out"; + + JFJochProcess process(reader, dataset->experiment, dataset->pixel_mask, config); + ProcessResult result; + REQUIRE_NOTHROW(result = process.Run()); + + CHECK_FALSE(result.cancelled); + CHECK(result.images_processed == 8); + REQUIRE(result.written_master_path.has_value()); + + { + // The _process.h5 links back to the source images and carries an azimuthal profile per image. + JFJochHDF5Reader out; + REQUIRE_NOTHROW(out.ReadFile("process_azint_out_process.h5")); + CHECK(out.GetNumberOfImages() == 8); + std::shared_ptr img; + REQUIRE_NOTHROW(img = out.LoadImage(0)); + REQUIRE(img); + CHECK_FALSE(img->ImageData().az_int_profile.empty()); + } + + reader.Close(); + remove("process_azint_in_master.h5"); + remove("process_azint_in_data_000001.h5"); + remove("process_azint_out_process.h5"); + REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); +} + +TEST_CASE("JFJochProcess_NoOutput", "[HDF5][Full]") { + WriteTestDataset("process_noout_in", 6); + + JFJochHDF5Reader reader; + REQUIRE_NOTHROW(reader.ReadFile("process_noout_in_master.h5")); + auto dataset = reader.GetDataset(); + + // Empty output prefix => process without writing any file. + ProcessConfig config; + config.mode = ProcessMode::AzimuthalIntegration; + config.nthreads = 3; + + JFJochProcess process(reader, dataset->experiment, dataset->pixel_mask, config); + auto result = process.Run(); + + CHECK_FALSE(result.cancelled); + CHECK(result.images_processed == 6); + CHECK_FALSE(result.written_master_path.has_value()); + + reader.Close(); + remove("process_noout_in_master.h5"); + remove("process_noout_in_data_000001.h5"); + REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); +} + +TEST_CASE("JFJochProcess_Cancel", "[HDF5][Full]") { + WriteTestDataset("process_cancel_in", 8); + + JFJochHDF5Reader reader; + REQUIRE_NOTHROW(reader.ReadFile("process_cancel_in_master.h5")); + auto dataset = reader.GetDataset(); + + ProcessConfig config; + config.mode = ProcessMode::AzimuthalIntegration; + config.nthreads = 2; + + JFJochProcess process(reader, dataset->experiment, dataset->pixel_mask, config); + process.Cancel(); // cancel before running: the worker loop stops immediately + auto result = process.Run(); + + CHECK(result.cancelled); + CHECK(result.images_processed == 0); + CHECK_FALSE(result.written_master_path.has_value()); + + reader.Close(); + remove("process_cancel_in_master.h5"); + remove("process_cancel_in_data_000001.h5"); + REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); +} diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index a887d94b..d3e88ed8 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -9,7 +9,7 @@ TARGET_LINK_LIBRARIES(jfjoch_extract_hkl JFJochReader) INSTALL(TARGETS jfjoch_extract_hkl RUNTIME COMPONENT viewer) ADD_EXECUTABLE(jfjoch_process jfjoch_process.cpp) -TARGET_LINK_LIBRARIES(jfjoch_process JFJochReader JFJochImageAnalysis JFJochWriter) +TARGET_LINK_LIBRARIES(jfjoch_process JFJochProcess JFJochReader JFJochImageAnalysis JFJochWriter) INSTALL(TARGETS jfjoch_process RUNTIME COMPONENT viewer) ADD_EXECUTABLE(jfjoch_scale jfjoch_scale.cpp) @@ -17,7 +17,7 @@ TARGET_LINK_LIBRARIES(jfjoch_scale JFJochReader JFJochImageAnalysis JFJochWriter INSTALL(TARGETS jfjoch_scale RUNTIME COMPONENT viewer) ADD_EXECUTABLE(jfjoch_azint jfjoch_azint.cpp) -TARGET_LINK_LIBRARIES(jfjoch_azint JFJochReader JFJochImageAnalysis JFJochWriter) +TARGET_LINK_LIBRARIES(jfjoch_azint JFJochProcess JFJochReader JFJochImageAnalysis JFJochWriter) INSTALL(TARGETS jfjoch_azint RUNTIME COMPONENT viewer) # On Windows these CLIs get getopt/getopt_long from the vendored wingetopt shim (libc has none). diff --git a/tools/jfjoch_azint.cpp b/tools/jfjoch_azint.cpp index cdc1d7bd..87cdf564 100644 --- a/tools/jfjoch_azint.cpp +++ b/tools/jfjoch_azint.cpp @@ -1,12 +1,10 @@ // SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only +#include +#include #include #include -#include -#include -#include -#include #include #include #include @@ -16,15 +14,8 @@ #include "../common/Definitions.h" #include "../common/DiffractionExperiment.h" #include "../common/PixelMask.h" -#include "../common/AzimuthalIntegrationMapping.h" -#include "../common/AzimuthalIntegrationProfile.h" -#include "../common/time_utc.h" #include "../common/print_license.h" -#include "../image_analysis/azint/AzIntEngineCPU.h" -#include "../image_analysis/image_preprocessing/ImagePreprocessorCPU.h" -#include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h" -#include "../common/JFJochReceiverPlots.h" -#include "../writer/FileWriter.h" +#include "../process/JFJochProcess.h" void print_usage() { std::cout << "Usage ./jfjoch_azint {} " << std::endl; @@ -113,6 +104,14 @@ bool parse_on_off(const char *arg, bool &out) { return false; } +namespace { + std::atomic g_active_process{nullptr}; + void handle_sigint(int) { + if (auto *p = g_active_process.load()) + p->Cancel(); + } +} + int main(int argc, char **argv) { for (int i = 0; i < argc; i++) std::cout << argv[i] << " "; @@ -211,6 +210,11 @@ int main(int argc, char **argv) { const std::string input_file = argv[optind]; logger.Verbose(verbose); + if (image_stride <= 0) { + logger.Error("Image stride must be positive"); + exit(EXIT_FAILURE); + } + // 1. Read input file JFJochHDF5Reader reader; try { @@ -227,32 +231,10 @@ int main(int argc, char **argv) { } logger.Info("Loaded dataset from {}", input_file); - if (image_stride <= 0) { - logger.Error("Image stride must be positive"); - exit(EXIT_FAILURE); - } - - const uint64_t total_images_in_file = reader.GetNumberOfImages(); - if (end_image < 0 || end_image > total_images_in_file) - end_image = total_images_in_file; - - const int images_to_process = (end_image - start_image) / image_stride; - if (images_to_process <= 0) { - logger.Warning("No images to process (Start: {}, End: {}, Stride: {}, Total: {})", - start_image, end_image, image_stride, total_images_in_file); - return 0; - } - - // 2. Set up experiment - defaults come from the input file, overridden by command line + // 2. Build experiment: defaults from the input file, overridden by command line. Output and + // runtime invariants are set inside JFJochProcess; here we only configure the geometry and + // azimuthal-integration settings. DiffractionExperiment experiment(dataset->experiment); - experiment.BitDepthImage(32).PixelSigned(true); - experiment.Mode(DetectorMode::Standard); - experiment.FilePrefix(output_prefix); - experiment.OverwriteExistingFiles(true); - experiment.SetFileWriterFormat(FileWriterFormat::NXmxLegacy); - experiment.ImagesPerTrigger(images_to_process); - experiment.NumTriggers(1); - if (beam_x.has_value()) experiment.BeamX_pxl(beam_x.value()); if (beam_y.has_value()) experiment.BeamY_pxl(beam_y.value()); if (detector_distance_mm.has_value()) experiment.DetectorDistance_mm(detector_distance_mm.value()); @@ -282,142 +264,33 @@ int main(int argc, char **argv) { azint_settings.IsPolarizationCorrection() ? "on" : "off", azint_settings.IsSolidAngleCorrection() ? "on" : "off"); - PixelMask pixel_mask = dataset->pixel_mask; - AzimuthalIntegrationMapping mapping(experiment, pixel_mask); + // 3. Run the shared azimuthal-integration workflow. + ProcessConfig config; + config.mode = ProcessMode::AzimuthalIntegration; + config.start_image = start_image; + config.end_image = end_image; + config.stride = image_stride; + config.nthreads = nthreads; + config.output_prefix = output_prefix; - JFJochReceiverPlots plots; - plots.Setup(experiment, mapping); + JFJochProcess process(reader, experiment, dataset->pixel_mask, config); - // 3. Set up the output file - StartMessage start_message; - experiment.FillMessage(start_message); - start_message.arm_date = dataset->arm_date; - start_message.az_int_bin_to_q = mapping.GetBinToQ(); - start_message.az_int_bin_to_two_theta = mapping.GetBinToTwoTheta(); - start_message.az_int_q_bin_count = mapping.GetQBinCount(); - start_message.az_int_phi_bin_count = mapping.GetAzimuthalBinCount(); - if (mapping.GetAzimuthalBinCount() > 1) - start_message.az_int_bin_to_phi = mapping.GetBinToPhi(); - start_message.pixel_mask["default"] = pixel_mask.GetMask(experiment); - start_message.master_suffix = "process"; - start_message.file_format = FileWriterFormat::NXmxIntegrated; - start_message.write_master_file = true; - start_message.write_images = false; - start_message.hdf5_source_data = reader.GetHDF5DataSource(start_image, images_to_process); + g_active_process = &process; + std::signal(SIGINT, handle_sigint); - std::unique_ptr writer; + ProcessResult result; try { - if (!output_prefix.empty()) - writer = std::make_unique(start_message); + result = process.Run(); } catch (const std::exception &e) { - logger.Error("Failed to initialize file writer: {}", e.what()); + logger.Error("Processing failed: {}", e.what()); exit(EXIT_FAILURE); } + g_active_process = nullptr; - logger.Info("Starting azimuthal integration of {} images (range {}-{}) using {} threads", - images_to_process, start_image, end_image, nthreads); - - // 4. Process images on N CPU threads - std::atomic next_image = 0; - std::atomic finished_count = 0; - std::atomic total_uncompressed_bytes = 0; - - const auto start_time = std::chrono::steady_clock::now(); - - auto worker = [&]() { - std::vector decompression_buffer; - ImagePreprocessorCPU preprocessor(experiment, pixel_mask); - ImagePreprocessorBuffer buffer(experiment.GetPixelsNum()); - AzIntEngineCPU azint(mapping); - AzimuthalIntegrationProfile profile(mapping); - - while (true) { - const int image_ordinal = next_image.fetch_add(1); - const int image_idx = start_image + image_ordinal * image_stride; - if (image_idx >= end_image) break; - - std::shared_ptr img; - try { - img = reader.GetRawImage(image_idx); - } catch (const std::exception &e) { - logger.Error("Failed to load image {}: {}", image_idx, e.what()); - continue; - } - if (!img) continue; - - DataMessage msg{}; - msg.image = img->image; - msg.number = image_ordinal; - msg.original_number = image_idx; - if (dataset->efficiency.size() > image_idx) - msg.image_collection_efficiency = dataset->efficiency[image_idx]; - total_uncompressed_bytes += msg.image.GetUncompressedSize(); - - const auto azint_start = std::chrono::steady_clock::now(); - try { - const uint8_t *image_ptr = msg.image.GetUncompressedPtr(decompression_buffer); - preprocessor.Analyze(buffer, image_ptr, msg.image.GetMode()); - azint.Run(buffer, profile); - } catch (const std::exception &e) { - logger.Error("Error integrating image {}: {}", image_idx, e.what()); - continue; - } - const auto azint_end = std::chrono::steady_clock::now(); - - msg.azint_time_s = std::chrono::duration(azint_end - azint_start).count(); - msg.processing_time_s = msg.azint_time_s; - msg.az_int_profile = profile.GetResult(); - msg.az_int_profile_count = profile.GetPixelCount(); - msg.az_int_profile_std = profile.GetStd(); - msg.bkg_estimate = profile.GetBkgEstimate(mapping.Settings()); - msg.run_number = experiment.GetRunNumber(); - msg.run_name = experiment.GetRunName(); - - plots.Add(msg, profile); - if (writer) - writer->Write(msg); - - const int done = finished_count.fetch_add(1) + 1; - if (done % 1000 == 0) { - const double elapsed_s = std::chrono::duration( - std::chrono::steady_clock::now() - start_time).count(); - const double frame_rate_hz = (elapsed_s > 0.0) ? (done / elapsed_s) : 0.0; - logger.Info("Integrated {} / {} images ({:.1f} Hz)", done, images_to_process, frame_rate_hz); - } - } - }; - - std::vector > futures; - futures.reserve(nthreads); - for (int i = 0; i < nthreads; ++i) - futures.push_back(std::async(std::launch::async, worker)); - for (auto &f: futures) - f.get(); - - const auto end_time = std::chrono::steady_clock::now(); - - // 5. Finalize output file - EndMessage end_msg; - end_msg.max_image_number = images_to_process; - end_msg.images_collected_count = images_to_process; - end_msg.images_sent_to_write_count = images_to_process; - end_msg.end_date = time_UTC(std::chrono::system_clock::now()); - end_msg.run_number = experiment.GetRunNumber(); - end_msg.run_name = experiment.GetRunName(); - end_msg.bkg_estimate = plots.GetBkgEstimate(); - end_msg.az_int_result["dataset"] = plots.GetAzIntProfile(); - - if (writer) { - writer->WriteHDF5(end_msg); - writer->Finalize(); - } - - // 6. Report statistics - const double processing_time = std::chrono::duration(end_time - start_time).count(); - const double throughput_MBs = static_cast(total_uncompressed_bytes) / (processing_time * 1e6); - const double frame_rate = static_cast(finished_count.load()) / processing_time; - - std::cout << fmt::format("Processing time: {:.2f} s", processing_time) << std::endl; - std::cout << fmt::format("Frame rate: {:.2f} Hz", frame_rate) << std::endl; - std::cout << fmt::format("Total throughput: {:.2f} MB/s", throughput_MBs) << std::endl; + // 4. Report statistics + std::cout << fmt::format("Processing time: {:.2f} s", result.processing_time_s) << std::endl; + std::cout << fmt::format("Frame rate: {:.2f} Hz", result.frame_rate_hz) << std::endl; + std::cout << fmt::format("Total throughput: {:.2f} MB/s", result.throughput_MBs) << std::endl; + if (result.cancelled) + logger.Warning("Processing was cancelled after {} images", result.images_processed); } diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 17174300..e9474cae 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -38,6 +38,8 @@ #include "../image_analysis/scale_merge/SearchSpaceGroup.h" #include "../image_analysis/WriteReflections.h" #include "../image_analysis/UpdateReflectionResolution.h" +#include "../process/JFJochProcess.h" +#include void print_usage() { std::cout << "Usage jfjoch_process {} " << std::endl; @@ -249,31 +251,12 @@ std::optional parse_lattice_arg(const char *arg) { return CrystalLattice(vals); } -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; +namespace { + std::atomic g_active_process{nullptr}; + void handle_sigint(int) { + if (auto *p = g_active_process.load()) + p->Cancel(); } - - 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) { @@ -741,424 +724,63 @@ int main(int argc, char **argv) { if (d_min_spot_finding > 0.0f) spot_settings.high_resolution_limit = d_min_spot_finding; - // Initialize Analysis Components - PixelMask pixel_mask = dataset->pixel_mask; + // Run the shared full-analysis workflow (rotation indexing + scaling/merging live in + // JFJochProcess; the experiment above carries all algorithm settings). + ProcessConfig config; + config.mode = ProcessMode::FullAnalysis; + config.start_image = start_image; + config.end_image = end_image; + config.stride = image_stride; + config.nthreads = nthreads; + config.output_prefix = output_prefix; + config.spot_finding = spot_settings; + config.rotation_indexing = rotation_indexing; + config.two_pass_rotation = two_pass_rotation; + config.reuse_rotation_spots = reuse_rotation_spots; + config.rotation_indexing_image_count = rotation_indexing_image_count; + config.forced_rotation_lattice = forced_rotation_lattice; + config.run_scaling = run_scaling; + config.scaling_iter = scaling_iter; + config.reference_data = reference_data; - // If dataset has a mask you wish to use, you might need to load it into pixel_mask here - // e.g. pixel_mask.LoadUserMask(dataset->pixel_mask, ...); + JFJochProcess process(reader, experiment, dataset->pixel_mask, config); - AzimuthalIntegrationMapping mapping(experiment, pixel_mask); - IndexerThreadPool indexer_pool(experiment.GetIndexingSettings()); + g_active_process = &process; + std::signal(SIGINT, handle_sigint); - // Statistics collector - JFJochReceiverPlots plots; - plots.Setup(experiment, mapping); - - StartMessage start_message; - experiment.FillMessage(start_message); - start_message.arm_date = dataset->arm_date; // Use original arm date - start_message.az_int_bin_to_q = mapping.GetBinToQ(); - start_message.az_int_q_bin_count = mapping.GetQBinCount(); - if (mapping.GetAzimuthalBinCount() > 1) - start_message.az_int_bin_to_phi = mapping.GetBinToPhi(); - - start_message.pixel_mask["default"] = pixel_mask.GetMask(experiment); - start_message.rois = experiment.ROI().ExportMetadata(); - if (!experiment.ROI().empty()) - start_message.roi_map = experiment.ExportROIMap(); - start_message.max_spot_count = experiment.GetMaxSpotCount(); - start_message.master_suffix = "process"; - - start_message.file_format = FileWriterFormat::NXmxIntegrated; - start_message.write_master_file = true; - start_message.write_images = false; - start_message.hdf5_source_data = reader.GetHDF5DataSource(start_image, images_to_process); - - std::unique_ptr writer; + ProcessResult result; try { - if (!output_prefix.empty()) - writer = std::make_unique(start_message); + result = process.Run(); } catch (const std::exception &e) { - logger.Error("Failed to initialize file writer: {}", e.what()); + logger.Error("Processing failed: {}", e.what()); exit(EXIT_FAILURE); } + g_active_process = nullptr; - std::atomic processed_count = 0; - std::atomic total_uncompressed_bytes = 0; + if (!result.merge_statistics_text.empty()) + std::cout << std::endl << result.merge_statistics_text << std::endl; - // Mimic JFJochReceiver lattice handling (IndexAndRefine handles the logic per thread, - // but we need a central accumulator or use the pool's functionality if IndexAndRefine wraps it) - // Here we will use per-thread IndexAndRefine which uses the shared thread pool. - - auto start_time = std::chrono::steady_clock::now(); - - IndexAndRefine indexer(experiment, &indexer_pool); - if (!reference_data.empty()) - indexer.ReferenceIntensities(reference_data); - - if (forced_rotation_lattice.has_value()) { - indexer.ForceRotationIndexerLattice(*forced_rotation_lattice); - logger.Info("Rotation indexer lattice forced externally - skipping first pass indexing"); - } else 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_idx; - 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"); + // Report statistics + std::cout << fmt::format("Processing time: {:.2f} s", result.processing_time_s) << std::endl; + std::cout << fmt::format("Frame rate: {:.2f} Hz", result.frame_rate_hz) << std::endl; + std::cout << fmt::format("Total throughput:{:.2f} MB/s", result.throughput_MBs) << std::endl; + if (result.indexing_rate.has_value()) + std::cout << fmt::format("Indexing rate: {:.2f}%", result.indexing_rate.value() * 100.0) << std::endl; + if (result.consensus_cell.has_value()) { + const auto &c = result.consensus_cell.value(); + std::cout << fmt::format("Unit cell: a={:.2f} b={:.2f} c={:.2f} alpha={:.2f} beta={:.2f} gamma={:.2f}", + c.a, c.b, c.c, c.alpha, c.beta, c.gamma) << std::endl; } - std::atomic finished_count = 0; - - auto worker = [&](int thread_id) { - // Spread the per-image GPU pipeline across all visible GPUs (round-robin). - // Must run before building the analysis resources so their CUDA streams/engines are - // created on this thread's device. - pin_gpu(); - - // Thread-local analysis resources - MXAnalysisWithoutFPGA analysis(experiment, mapping, pixel_mask, indexer); - - AzimuthalIntegrationProfile profile(mapping); - - while (true) { - int image_ordinal = processed_count.fetch_add(1); - int image_idx = start_image + image_ordinal * image_stride; - - if (image_idx >= end_image) break; - - // Load Image - std::shared_ptr img; - try { - img = reader.GetRawImage(image_idx); - } catch (const std::exception &e) { - logger.Error("Failed to load image {}: {}", image_idx, e.what()); - continue; - } - - if (!img) continue; - - DataMessage msg{}; - - msg.image = img->image; - msg.number = image_ordinal; - msg.original_number = image_idx; - if (dataset->efficiency.size() > image_idx) msg.image_collection_efficiency = dataset->efficiency[image_idx]; - - total_uncompressed_bytes += msg.image.GetUncompressedSize(); - - auto image_start_time = std::chrono::high_resolution_clock::now(); - - // Analyze - try { - analysis.Analyze(msg, profile, spot_settings); - } catch (const std::exception &e) { - logger.Error("Error analyzing image {}: {}", image_idx, e.what()); - continue; - } - - auto image_end_time = std::chrono::high_resolution_clock::now(); - std::chrono::duration image_duration = image_end_time - image_start_time; - - msg.processing_time_s = image_duration.count(); - msg.run_number = experiment.GetRunNumber(); - msg.run_name = experiment.GetRunName(); - - plots.Add(msg, profile); - // Write Result - if (writer) - writer->Write(msg); - - finished_count.fetch_add(1); - - // Progress log - if (image_ordinal > 0 && image_ordinal % 100 == 0) { - std::optional indexing_rate = plots.GetIndexingRate(); - - const auto now = std::chrono::steady_clock::now(); - const double elapsed_s = std::chrono::duration(now - start_time).count(); - const int processed_images = finished_count.load(); - const double frame_rate_hz = (elapsed_s > 0.0) ? (processed_images / elapsed_s) : 0.0; - - if (indexing_rate.has_value()) { - logger.Info("Processed {} / {} images ({:.2f} Hz, indexing rate {:.1f}%)", - image_ordinal, images_to_process, - frame_rate_hz, - indexing_rate.value() * 100.0f); - } else { - logger.Info("Processed {} / {} images ({:.2f} Hz, indexing rate N/A)", - image_ordinal, images_to_process, - frame_rate_hz); - } - } - } - - // Finalize per-thread indexing (if any per-thread aggregation is needed, though pool handles most) - // IndexAndRefine doesn't have a per-thread finalize that returns a lattice, - // the main lattice determination is usually done on the aggregated results in the pool or main thread - }; - - // Launch threads - std::vector > futures; - futures.reserve(nthreads); - for (int i = 0; i < nthreads; ++i) { - futures.push_back(std::async(std::launch::async, worker, i)); - } - - // Wait for completion - for (auto &f: futures) { - f.get(); - } - - auto end_time = std::chrono::steady_clock::now(); - - // 5. Finalize Statistics and Write EndMessage - EndMessage end_msg; - end_msg.max_image_number = images_to_process; - end_msg.images_collected_count = images_to_process; - end_msg.images_sent_to_write_count = images_to_process; - end_msg.end_date = time_UTC(std::chrono::system_clock::now()); - end_msg.run_number = experiment.GetRunNumber(); - end_msg.run_name = experiment.GetRunName(); - - // Gather statistics from plots - end_msg.bkg_estimate = plots.GetBkgEstimate(); - end_msg.indexing_rate = plots.GetIndexingRate(); - end_msg.az_int_result["dataset"] = plots.GetAzIntProfile(); - - // Finalize Indexing (Global) to get rotation lattice - // We create a temporary IndexAndRefine to call Finalize() which aggregates pool results - const auto rotation_indexer_ret = indexer.FinalizeRotationIndexing(); - - if (rotation_indexer_ret.has_value()) { - end_msg.rotation_lattice = rotation_indexer_ret->lattice; - end_msg.rotation_lattice_type = LatticeMessage{ - .centering = rotation_indexer_ret->search_result.centering, - .niggli_class = rotation_indexer_ret->search_result.niggli_class, - .crystal_system = rotation_indexer_ret->search_result.system - }; - logger.Info("Rotation Indexing found lattice"); - } - - std::vector merging_mask_uc(images_to_process, 1); - - auto consensus_start_time = std::chrono::steady_clock::now(); - const auto consensus_cell = indexer.GetConsensusUnitCell(); - auto consensus_end_time = std::chrono::steady_clock::now(); - auto consensus_duration = std::chrono::duration(consensus_end_time - consensus_start_time).count(); - - if (consensus_cell) { - logger.Info("Consensus unit cell found in {:.2f} ms", consensus_duration * 1e3); - logger.Info("UC: a={:.2f} b={:.2f} c={:.2f} alpha={:.2f} beta={:.2f} gamma={:.2f}", - consensus_cell->a, consensus_cell->b, consensus_cell->c, - consensus_cell->alpha, consensus_cell->beta, consensus_cell->gamma); - } else - logger.Info("Consensus unit cell not found - calculation tool {:.2f} ms", consensus_duration * 1e3); - end_msg.unit_cell = consensus_cell; - - if (end_msg.indexing_rate.has_value() - && end_msg.indexing_rate > 0 - && (run_scaling || !reference_data.empty())) { - - 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(); - for (int i = 0; i < scaling_iter; i++) { - auto iter_start = std::chrono::steady_clock::now(); - auto merge_result = MergeAll(experiment, indexer.GetIntegrationOutcome(), false); - scale_result = indexer.ScaleAllImages(merge_result); - scale_result.SaveToFile(output_prefix + "_iter" + std::to_string(i) + "_scale.dat"); - auto iter_end = std::chrono::steady_clock::now(); - double iter_time = std::chrono::duration(iter_end - iter_start).count(); - logger.Info("Scaling iteration {} took {:.3f} seconds", i, iter_time); - } - - end_msg.image_scale_factor = scale_result.image_scale_g; - end_msg.image_scale_cc = scale_result.image_cc; - end_msg.image_scale_mosaicity = scale_result.mosaicity_deg; - end_msg.image_scale_b_factor = scale_result.image_bfactor_Ang2; - - auto scale_end = std::chrono::steady_clock::now(); - double scale_time = std::chrono::duration(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"); - } - - // --- Phase 1 diagnostic: distribution of the per-image scale CC vs the - // reference. High per-image CC with low merged CC1/2 => each image is fine - // and the cross-image merge is the problem; low per-image CC => the - // per-image extraction itself is noise. Logged for any reference run. - { - std::vector ccs; - for (const auto &i : indexer.GetIntegrationOutcome()) - if (i.image_scale_cc && std::isfinite(*i.image_scale_cc)) - ccs.push_back(*i.image_scale_cc); - if (!ccs.empty()) { - std::sort(ccs.begin(), ccs.end()); - const double mean = std::accumulate(ccs.begin(), ccs.end(), 0.0) / ccs.size(); - auto q = [&](double f) { return ccs[std::min(ccs.size() - 1, static_cast(f * ccs.size()))]; }; - logger.Info("Per-image scale CC vs reference: n={} mean={:.3f} median={:.3f} " - "p10={:.3f} p90={:.3f} min={:.3f} max={:.3f}", - ccs.size(), mean, q(0.5), q(0.1), q(0.9), ccs.front(), ccs.back()); - } else { - logger.Info("Per-image scale CC vs reference: none available"); - } - } - - auto merge_start = std::chrono::steady_clock::now(); - - MergeOnTheFly merge_engine(experiment); - if (consensus_cell.has_value()) - merge_engine.ReferenceCell(*consensus_cell); - // Fit the global error model (sigma'^2 = a*sigma^2 + (b*I)^2, XDS form) from the - // spread of symmetry equivalents before merging, so the merged sigmas reflect the - // real scatter rather than counting statistics alone. - merge_engine.RefineErrorModel(indexer.GetIntegrationOutcome()); - if (merge_engine.ErrorModelActive()) { - const double a = merge_engine.ErrorModelA(); - const double b = merge_engine.ErrorModelB(); - const double isa = (b > 0.0) ? 1.0 / b : std::numeric_limits::infinity(); - logger.Info("Error model: sigma'^2 = {:.3f} sigma^2 + ({:.4f} I)^2 ISa = {:.1f}", a, b, isa); - } - const auto &merge_outcomes = indexer.GetIntegrationOutcome(); - std::vector dcch_reject; - if (delta_cchalf_nsigma) { - dcch_reject = merge_engine.DeltaCChalfReject(merge_outcomes, *delta_cchalf_nsigma); - const size_t nrej = std::count(dcch_reject.begin(), dcch_reject.end(), static_cast(1)); - logger.Info("CC1/2-delta rejection (deltaCChalf < mean - {:.1f} stddev) removed {} / {} images", - *delta_cchalf_nsigma, nrej, dcch_reject.size()); - } - for (size_t i = 0; i < merge_outcomes.size(); ++i) { - if (!dcch_reject.empty() && dcch_reject[i]) - continue; - merge_engine.AddImage(merge_outcomes[i]); - } - if (merge_engine.RejectedCount() > 0) - logger.Info("Outlier rejection (>{:.1f} sigma from the per-reflection median) removed {} observations", - experiment.GetScalingSettings().GetOutlierRejectNsigma(), merge_engine.RejectedCount()); - - auto merged_reflections = merge_engine.ExportReflections(); - auto merged_statistics = merge_engine.MergeStats(merged_reflections, indexer.GetIntegrationOutcome(), reference_data); - - auto merge_end = std::chrono::steady_clock::now(); - double merge_time = std::chrono::duration(merge_end - merge_start).count(); - - logger.Info("Merge completed in {:.2f} s ({} unique reflections)", merge_time, merged_reflections.size()); - - if (!experiment.GetGemmiSpaceGroup().has_value()) { - logger.Info("Searching for space group from P1-merged reflections ..."); - - SearchSpaceGroupOptions sg_opts; - sg_opts.crystal_system.reset(); - sg_opts.centering = '\0'; - sg_opts.merge_friedel = experiment.GetScalingSettings().GetMergeFriedel(); - sg_opts.d_min_limit_A = experiment.GetScalingSettings().GetHighResolutionLimit_A().value_or(0.0); - sg_opts.min_i_over_sigma = 0.0; - sg_opts.min_operator_cc = 0.80; - sg_opts.min_pairs_per_operator = 20; - sg_opts.min_total_compared = 100; - sg_opts.test_systematic_absences = true; - - const auto sg_search = SearchSpaceGroup(merged_reflections, sg_opts); - - std::cout << std::endl << SearchSpaceGroupResultToText(sg_search) << std::endl << std::endl; - } - - std::cout << merged_statistics; - - if (consensus_cell && !output_prefix.empty()) - WriteReflections(merged_reflections, *consensus_cell, experiment, output_prefix); - } - - // Write End Message - if (writer) { - writer->WriteHDF5(end_msg); - auto stats = writer->Finalize(); - } - - // 6. Report Statistics to Console - double processing_time = std::chrono::duration(end_time - start_time).count(); - double throughput_MBs = static_cast(total_uncompressed_bytes) / (processing_time * 1e6); - double frame_rate = static_cast(images_to_process) / processing_time; - - std::cout << fmt::format("Processing time: {:.2f} s", processing_time) << std::endl; - std::cout << fmt::format("Frame rate: {:.2f} Hz", frame_rate) << std::endl; - std::cout << fmt::format("Total throughput:{:.2f} MB/s", throughput_MBs) << std::endl; - - // Print extended stats similar to Receiver - if (end_msg.indexing_rate.has_value()) { - std::cout << fmt::format("Indexing rate: {:.2f}%", end_msg.indexing_rate.value() * 100.0) << std::endl; - } - - auto image_mean_time = plots.GetMeanProcessingTime(); + const auto &t = result.mean_processing_time; std::cout << fmt::format( - "Per-image time: (mean; milliseconds): decompress {:.2f} preprocess {:.2f} azint {:.2f} spot finding {:.2f} indexing {:.2f} refinement {:.2f} indexing analysis {:.2f} prediction {:.2f} integration {:.2f} scaling {:.2f} total {:.2f}", - image_mean_time.compression * 1e3, - image_mean_time.preprocessing * 1e3, - image_mean_time.azint * 1e3, - image_mean_time.spot_finding * 1e3, - image_mean_time.indexing * 1e3, - image_mean_time.refinement * 1e3, - image_mean_time.indexing_analysis * 1e3, - image_mean_time.bragg_prediction * 1e3, - image_mean_time.integration * 1e3, - image_mean_time.image_scale * 1e3, - image_mean_time.processing * 1e3) - << std::endl; + "Per-image time (mean; ms): decompress {:.2f} preprocess {:.2f} azint {:.2f} spot finding {:.2f} " + "indexing {:.2f} refinement {:.2f} indexing analysis {:.2f} prediction {:.2f} integration {:.2f} " + "scaling {:.2f} total {:.2f}", + t.compression * 1e3, t.preprocessing * 1e3, t.azint * 1e3, t.spot_finding * 1e3, + t.indexing * 1e3, t.refinement * 1e3, t.indexing_analysis * 1e3, t.bragg_prediction * 1e3, + t.integration * 1e3, t.image_scale * 1e3, t.processing * 1e3) << std::endl; + + if (result.cancelled) + logger.Warning("Processing was cancelled after {} images", result.images_processed); }