rotation indexer: fix sub-range crash, explicit angle, settings + guards
Running rotation indexing on a sub-range (e.g. images 60-120) segfaulted: the first pass passed the *global* image number to RotationIndexer::ProcessImage, which indexes v_ (sized to the run's image count) -> out-of-bounds write. - ProcessImage now takes an explicit mid-exposure angle (optional; falls back to the goniometer at the image index), so the indexer no longer assumes its slot index equals the goniometer image index. IndexAndRefine supplies it via RotationAngle(), matching the angle used for prediction. Added a bounds guard in ProcessImage so a bad index can never corrupt memory. - JFJochProcess feeds the rotation indexer the local ordinal (not the global index), and shifts the goniometer (start += start_image*incr, incr *= stride, per-image wedge preserved) so local index i maps to the angle of original image start+i*stride - fixing rotation angles for the whole sub-range pipeline (prediction, refinement, output), not just the indexer. - Expose "Rotation images" (number used for the first pass) in the job dialog, enabled when rotation indexing is on. (Count > available is already clamped by select_equally_spaced_image_ordinals.) Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
@@ -28,9 +28,16 @@ IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool
|
||||
unit_cells.resize(x.GetImageNum());
|
||||
}
|
||||
|
||||
std::optional<float> 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<float>(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();
|
||||
|
||||
@@ -82,6 +82,7 @@ class IndexAndRefine {
|
||||
BraggPrediction &prediction,
|
||||
const IndexingOutcome &outcome,
|
||||
IntegrationOutcome &i_outcome);
|
||||
std::optional<float> RotationAngle(int64_t image) const; // mid-exposure angle for the indexer
|
||||
public:
|
||||
IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool *indexer);
|
||||
|
||||
|
||||
@@ -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<SpotToSave> &spots) {
|
||||
void RotationIndexer::ProcessImage(int64_t image, const std::vector<SpotToSave> &spots,
|
||||
std::optional<float> 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<int64_t>(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) {
|
||||
|
||||
@@ -29,6 +29,9 @@ class RotationIndexer {
|
||||
const bool index_ice_rings;
|
||||
|
||||
std::vector<std::vector<SpotToSave>> 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<std::optional<float>> angle_deg_;
|
||||
|
||||
std::optional<GoniometerAxis> axis_;
|
||||
const DiffractionGeometry geom_;
|
||||
@@ -43,7 +46,10 @@ class RotationIndexer {
|
||||
std::optional<CrystalLattice> indexed_lattice;
|
||||
public:
|
||||
RotationIndexer(const DiffractionExperiment& x, IndexerThreadPool& indexer);
|
||||
void ProcessImage(int64_t image, const std::vector<SpotToSave>& 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<SpotToSave>& spots,
|
||||
std::optional<float> angle_deg = std::nullopt);
|
||||
void RunIndexing();
|
||||
std::optional<RotationIndexerResult> GetLattice() const;
|
||||
void ForceLattice(const CrystalLattice& lattice);
|
||||
|
||||
@@ -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<float>(start_image),
|
||||
incr * static_cast<float>(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) {
|
||||
|
||||
@@ -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<int>(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;
|
||||
|
||||
@@ -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;
|
||||
};
|
||||
|
||||
|
||||
Reference in New Issue
Block a user