diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 61b30861..10ba4991 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -28,9 +28,16 @@ IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool unit_cells.resize(x.GetImageNum()); } +std::optional IndexAndRefine::RotationAngle(int64_t image) const { + // Mid-exposure rotation angle for image index `image`, matching the angle used for prediction. + if (const auto g = experiment.GetGoniometer()) + return g->GetAngle_deg(static_cast(image)) + g->GetWedge_deg() / 2.0f; + return std::nullopt; +} + void IndexAndRefine::AddImageToRotationIndexer(DataMessage &msg) { if (rotation_indexer) - rotation_indexer->ProcessImage(msg.number, msg.spots); + rotation_indexer->ProcessImage(msg.number, msg.spots, RotationAngle(msg.number)); } IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetryRotation(DataMessage &msg) { @@ -44,7 +51,7 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetryRotat if (!result.has_value()) { auto rot_cnt = rotation_indexer_counter.Process(msg.number); if (rot_cnt.first) - rotation_indexer->ProcessImage(msg.number, msg.spots); + rotation_indexer->ProcessImage(msg.number, msg.spots, RotationAngle(msg.number)); if (rot_cnt.second) rotation_indexer->RunIndexing(); result = rotation_indexer->GetLattice(); diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index b8b3e9f2..5f42194d 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -82,6 +82,7 @@ class IndexAndRefine { BraggPrediction &prediction, const IndexingOutcome &outcome, IntegrationOutcome &i_outcome); + std::optional RotationAngle(int64_t image) const; // mid-exposure angle for the indexer public: IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool *indexer); diff --git a/image_analysis/rotation_indexer/RotationIndexer.cpp b/image_analysis/rotation_indexer/RotationIndexer.cpp index 66f0dc98..f7a6de3b 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.cpp +++ b/image_analysis/rotation_indexer/RotationIndexer.cpp @@ -13,6 +13,7 @@ RotationIndexer::RotationIndexer(const DiffractionExperiment &x, IndexerThreadPo : experiment(x), index_ice_rings(x.GetIndexingSettings().GetIndexIceRings()), v_(experiment.GetImageNum()), + angle_deg_(experiment.GetImageNum()), axis_(x.GetGoniometer()), geom_(x.GetDiffractionGeometry()), updated_geom_(geom_), @@ -28,7 +29,7 @@ void RotationIndexer::RunIndexing() { coords.reserve(max_spots_per_image * v_.size()); for (int i = 0; i < v_.size(); i++) { - const float angle_deg = axis_->GetAngle_deg(i) + axis_->GetWedge_deg() / 2.0f; + const float angle_deg = angle_deg_[i].value_or(axis_->GetAngle_deg(i) + axis_->GetWedge_deg() / 2.0f); const auto rot = axis_->GetTransformationAngle(angle_deg); for (const auto &s: v_[i]) @@ -115,19 +116,26 @@ void RotationIndexer::RunIndexing() { } } -void RotationIndexer::ProcessImage(int64_t image, const std::vector &spots) { +void RotationIndexer::ProcessImage(int64_t image, const std::vector &spots, + std::optional angle_deg) { std::unique_lock ul(m); // For non-rotation just ignore the whole procedure if (!axis_) return; + // Guard: `image` is a slot in [0, image count); a bad index (e.g. a global number for a subset + // run) must not corrupt memory. + if (image < 0 || image >= static_cast(v_.size())) + return; + if (accumulated_spots >= max_spots) return; if (indexed_lattice) return; + angle_deg_[image] = angle_deg; v_[image].reserve(spots.size()); for (const auto &s: spots) { diff --git a/image_analysis/rotation_indexer/RotationIndexer.h b/image_analysis/rotation_indexer/RotationIndexer.h index ca6f4023..2878972b 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.h +++ b/image_analysis/rotation_indexer/RotationIndexer.h @@ -29,6 +29,9 @@ class RotationIndexer { const bool index_ice_rings; std::vector> v_; + // Per-image rotation angle (mid-exposure, deg) supplied by the caller; falls back to the + // goniometer evaluated at the image index when not given. + std::vector> angle_deg_; std::optional axis_; const DiffractionGeometry geom_; @@ -43,7 +46,10 @@ class RotationIndexer { std::optional indexed_lattice; public: RotationIndexer(const DiffractionExperiment& x, IndexerThreadPool& indexer); - void ProcessImage(int64_t image, const std::vector& spots); + // angle_deg is the image's mid-exposure rotation angle; if omitted, the goniometer angle at + // `image` is used (only valid when `image` is the goniometer's own image index). + void ProcessImage(int64_t image, const std::vector& spots, + std::optional angle_deg = std::nullopt); void RunIndexing(); std::optional GetLattice() const; void ForceLattice(const CrystalLattice& lattice); diff --git a/process/JFJochProcess.cpp b/process/JFJochProcess.cpp index 6167a449..1cd76517 100644 --- a/process/JFJochProcess.cpp +++ b/process/JFJochProcess.cpp @@ -99,6 +99,21 @@ ProcessResult JFJochProcess::Run(JFJochProcessObserver *observer) { if (full) experiment_.Compression(CompressionAlgorithm::BSHUF_LZ4); + // The pipeline indexes images 0..N-1 within this run; if we process a sub-range/strided + // selection, shift the goniometer so local index i maps to the angle of original image + // start+i*stride (keeping the per-image rotation wedge), otherwise rotation angles would be + // wrong for any start_image != 0. + if (const auto g = experiment_.GetGoniometer(); + g.has_value() && (start_image != 0 || config_.stride != 1)) { + const float incr = g->GetIncrement_deg(); + GoniometerAxis shifted(g->GetName(), + g->GetStart_deg() + incr * static_cast(start_image), + incr * static_cast(config_.stride), + g->GetAxis(), g->GetHelicalStep()); + shifted.ScreeningWedge(g->GetScreeningWedge().value_or(incr)); + experiment_.Goniometer(shifted); + } + AzimuthalIntegrationMapping mapping(experiment_, pixel_mask_); JFJochReceiverPlots plots; @@ -165,7 +180,7 @@ ProcessResult JFJochProcess::Run(JFJochProcessObserver *observer) { if (cancelled_) break; const int image_idx = start_image + ordinal * config_.stride; DataMessage msg{}; - msg.number = image_idx; + msg.number = ordinal; // index into the rotation indexer (0..images_to_process-1) msg.original_number = image_idx; try { if (config_.reuse_rotation_spots) { diff --git a/viewer/windows/JFJochProcessingJobsWindow.cpp b/viewer/windows/JFJochProcessingJobsWindow.cpp index 9e3ad20f..ad41f9fa 100644 --- a/viewer/windows/JFJochProcessingJobsWindow.cpp +++ b/viewer/windows/JFJochProcessingJobsWindow.cpp @@ -138,14 +138,21 @@ int JFJochProcessingJobsWindow::askJob(const ReprocessingInputs &inputs, JobSpec auto *rotation = new QCheckBox("Rotation indexing (two-pass)", &dlg); rotation->setChecked(inputs.experiment.GetGoniometer().has_value()); + + auto *rotation_images = new QSpinBox(&dlg); // images used to find the rotation lattice + rotation_images->setRange(2, 1'000'000); + rotation_images->setValue(spec.rotation_images); // carries the default (30) + auto *scaling = new QCheckBox("Scale && merge", &dlg); auto sync_full = [&] { const bool full = mode->currentData().toInt() == static_cast(ProcessMode::FullAnalysis); rotation->setEnabled(full); scaling->setEnabled(full); + rotation_images->setEnabled(full && rotation->isChecked()); }; connect(mode, &QComboBox::currentIndexChanged, &dlg, sync_full); + connect(rotation, &QCheckBox::toggled, &dlg, sync_full); sync_full(); auto *form = new QFormLayout; @@ -156,6 +163,7 @@ int JFJochProcessingJobsWindow::askJob(const ReprocessingInputs &inputs, JobSpec form->addRow(save); form->addRow("Output prefix", prefix); form->addRow(rotation); + form->addRow("Rotation images", rotation_images); form->addRow(scaling); int result = 0; @@ -186,6 +194,7 @@ int JFJochProcessingJobsWindow::askJob(const ReprocessingInputs &inputs, JobSpec spec.save = save->isChecked(); spec.prefix = prefix->text(); spec.rotation = rotation->isEnabled() && rotation->isChecked(); + spec.rotation_images = rotation_images->value(); spec.scaling = scaling->isEnabled() && scaling->isChecked(); } return result; @@ -202,6 +211,7 @@ ProcessConfig JFJochProcessingJobsWindow::buildConfig(const JobSpec &spec, const if (spec.mode == ProcessMode::FullAnalysis) { config.rotation_indexing = spec.rotation; config.two_pass_rotation = true; + config.rotation_indexing_image_count = spec.rotation_images; config.run_scaling = spec.scaling; } return config; diff --git a/viewer/windows/JFJochProcessingJobsWindow.h b/viewer/windows/JFJochProcessingJobsWindow.h index 27aad0f7..b185f2f1 100644 --- a/viewer/windows/JFJochProcessingJobsWindow.h +++ b/viewer/windows/JFJochProcessingJobsWindow.h @@ -65,6 +65,7 @@ private: bool save = true; QString prefix; bool rotation = false; + int rotation_images = 30; // images used to find the rotation lattice (first pass) bool scaling = false; };