mx/Jungfraujoch
3
2 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 93 Wiki Activity

rugnux: move the image pass to RugnuxImagePass.cpp

Browse Source 
The stretch of RunPipeline from the per-image loop to the end message becomes the private member
ImagePass(PipelineLocals &), its body unchanged behind the lines that bind the locals it reads.
Two of the values it computed are also read by the scale/merge block after it, so they move up to
where PipelineLocals is built and travel in it: the bounds of the geometry pre-pass
(prepass_images, loop_end), computed from the same inputs as before. PipelineLocals also carries
`full` now.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi
This commit is contained in:
leonarski_fandClaude Opus 5.5 committed 2026-10-07 15:54:58 +02:00
1 parent 2e0a3bb255
commit 7bcdf2a081
5 files changed
+598 -477
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

rugnux/CMakeLists.txt
+1
View File
@@ -8,6 +8,7 @@ ADD_LIBRARY(JFJochRugnux STATIC
RugnuxCalibration.cpp
RugnuxPreScan.cpp
RugnuxPasses.cpp
RugnuxImagePass.cpp
RugnuxDefaults.cpp
RugnuxDefaults.h
Rugnux.h
rugnux/Rugnux.cpp
+10 -477
View File
@@ -1112,6 +1112,12 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
const auto start_time = std::chrono::steady_clock::now();
// The geometry pre-pass integrates only the first prepass_fraction of the sweep (pass 2 reads it all).
const int prepass_end = start_image
+ static_cast<int>(std::lround((end_image - start_image) * config_.prepass_fraction));
const int prepass_images = (prepass_end - start_image + config_.stride - 1) / config_.stride;
const int loop_end = geometry_prepass ? prepass_end : end_image;
// Declared here rather than where they are first filled, so that the phases below can share them
// through PipelineLocals: the bytes the image loop read, the end message (also written to the file),
// and the unmerged MTZ when it is built beside the P1 cross-check merge (see there) with the
@@ -1123,10 +1129,11 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
std::future<gemmi::Mtz> unmerged_mtz;
PipelineLocals p{
.logger = logger, .result = result, .observer = observer,
.geometry_prepass = geometry_prepass, .calibration = calibration,
.geometry_prepass = geometry_prepass, .full = full, .calibration = calibration,
.calibration_spots = calibration_spots, .per_image_analysis = per_image_analysis,
.write_files = write_files, .start_image = start_image, .end_image = end_image,
.images_to_process = images_to_process, .dataset = dataset, .mapping = mapping, .plots = plots,
.images_to_process = images_to_process, .prepass_images = prepass_images, .loop_end = loop_end,
.dataset = dataset, .mapping = mapping, .plots = plots,
.writer = writer, .writer_queue = writer_queue, .indexer_pool = indexer_pool, .indexer = indexer,
.indexing_bound_A = indexing_bound_A, .census_future = census_future, .superseded = superseded,
.ambiguity_reference = ambiguity_reference, .start_time = start_time,
@@ -3325,481 +3332,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
indexer->ProbeSupercell(frames);
}
// 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<int> next_ordinal = 0;
std::atomic<int> finished_count = 0;
// Images are written at their own ordinal, so a frame that fails leaves a HOLE rather than
// shifting everything after it up. The extent of the written datasets is therefore the highest
// ordinal reached, not the number that succeeded - counting successes makes the file one image
// short for every failure, and the image it loses is the LAST one, not the one that failed.
std::atomic<int> max_written_ordinal = -1;
auto note_written = [&max_written_ordinal](int ordinal) {
int prev = max_written_ordinal.load(std::memory_order_relaxed);
while (prev < ordinal
&& !max_written_ordinal.compare_exchange_weak(prev, ordinal, std::memory_order_relaxed)) {
}
};
// What the integrator did to the reflections it was handed, summed over the workers as each one
// finishes (BraggIntegrationCounts). Reading an engine's counts synchronises its stream, so it is
// done once per worker at the end rather than per image.
std::atomic<uint64_t> bragg_predicted = 0, bragg_starved = 0, bragg_starved_nb = 0,
bragg_fallback = 0;
// Calibration by spots: the pooled spot list of the whole run. Ring clustering is O(n^2) in it and
// the Hough circle centre O(n^3) in its first few hundred, so each image contributes a fair share of
// a fixed budget instead of everything it found - a powder ring is over-determined either way.
constexpr size_t CALIBRATION_MAX_SPOTS = 20000;
const size_t calibration_spots_per_image =
std::max<size_t>(1, CALIBRATION_MAX_SPOTS / static_cast<size_t>(images_to_process));
std::mutex calibration_spot_m;
std::vector<SpotToSave> calibration_spot_list;
// Higher-order beam contamination, counted on the spots the lattice did not take. Each worker
// keeps its own tally and adds it in when it is done, so the total is a sum of counts and does
// not depend on how the images were shared out.
std::mutex harmonic_m;
HarmonicEvidence harmonic_evidence;
// Where the indexed spots sit against their predicted reflections, measured on the geometry
// pre-pass for the canonical pass's footprint (SpotFootprint.h).
const bool measure_offsets = geometry_prepass && prescan_footprint_.has_value();
std::mutex footprint_m;
std::vector<FootprintOffset> footprint_offsets;
auto azint_worker = [&]() {
std::vector<uint8_t> decompression_buffer;
std::shared_ptr<JFJochReaderRawImage> img;
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;
// Reuse the raw image frame to frame - a fresh 4-byte 16 Mpx frame is 72 MB of
// untouched pages the decoder faults in one by one - except while the writer thread
// still holds the previous one's pixels.
if (!img || img.use_count() > 1)
img = std::make_shared<JFJochReaderRawImage>();
bool read = false;
try {
read = reader_.ReadRawImage(image_idx, *img);
} catch (const std::exception &e) {
if (IsFatalResourceError(e)) throw;
logger.Error("Failed to load image {}: {}", image_idx, e.what());
continue;
}
if (!read) 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) {
if (IsFatalResourceError(e)) throw;
logger.Error("Error integrating image {}: {}", image_idx, e.what());
continue;
}
msg.azint_time_s = std::chrono::duration<float>(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.ice_ring_score = profile.GetIceRingScore(mapping.Settings(),
config_.spot_finding.ice_ring_width_Q_recipA);
msg.run_number = experiment_.GetRunNumber();
msg.run_name = experiment_.GetRunName();
plots.Add(msg, profile);
if (writer_queue) writer_queue->Post(msg, img);
note_written(ordinal);
if (observer) observer->OnImageProcessed(msg);
const int done = finished_count.fetch_add(1) + 1;
if (observer) observer->OnProgress(done, images_to_process);
}
};
// The geometry pre-pass integrates only the first prepass_fraction of the sweep (pass 2 reads it all).
const int prepass_end = start_image
+ static_cast<int>(std::lround((end_image - start_image) * config_.prepass_fraction));
const int prepass_images = (prepass_end - start_image + config_.stride - 1) / config_.stride;
const int loop_end = geometry_prepass ? prepass_end : end_image;
auto full_worker = [&]() {
pin_gpu(); // round-robin per worker thread; must precede engine construction
// Before the analysis engine, which page-locks these bytes for its uploads and unregisters
// them when it is destroyed: it must not outlive the buffer it was given.
std::shared_ptr<JFJochReaderRawImage> img;
MXAnalysisWithoutFPGA analysis(experiment_, mapping, pixel_mask_, *indexer,
/*enable_fused_adaptive_gpu=*/true);
AzimuthalIntegrationProfile profile(mapping);
HarmonicEvidence worker_harmonic;
std::vector<FootprintOffset> worker_offsets;
while (!cancelled_) {
const int ordinal = next_ordinal.fetch_add(1);
const int image_idx = start_image + ordinal * config_.stride;
if (image_idx >= loop_end) break;
// Reuse the raw image frame to frame - a fresh 4-byte 16 Mpx frame is 72 MB of
// untouched pages the decoder faults in one by one - except while the writer thread
// still holds the previous one's pixels.
if (!img || img.use_count() > 1)
img = std::make_shared<JFJochReaderRawImage>();
bool read = false;
try {
read = reader_.ReadRawImage(image_idx, *img);
} catch (const std::exception &e) {
if (IsFatalResourceError(e)) throw;
logger.Error("Failed to load image {}: {}", image_idx, e.what());
continue;
}
if (!read) 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) {
if (IsFatalResourceError(e)) throw;
logger.Error("Error analyzing image {}: {}", image_idx, e.what());
continue;
}
msg.processing_time_s = std::chrono::duration<float>(std::chrono::steady_clock::now() - t0).count();
msg.run_number = experiment_.GetRunNumber();
msg.run_name = experiment_.GetRunName();
if (calibration) {
std::lock_guard lock(calibration_spot_m);
for (size_t i = 0; i < msg.spots.size() && i < calibration_spots_per_image; ++i)
calibration_spot_list.push_back(msg.spots[i]);
}
AddHarmonicEvidence(msg.spots, experiment_.GetDiffractionGeometry(), worker_harmonic);
if (measure_offsets)
MeasureFootprintOffsets(msg.spots, msg.reflections, experiment_.GetBeamX_pxl(),
experiment_.GetBeamY_pxl(), worker_offsets);
plots.Add(msg, profile);
if (writer_queue) writer_queue->Post(msg, img);
note_written(ordinal);
if (observer) observer->OnImageProcessed(msg);
const int done = finished_count.fetch_add(1) + 1;
if (observer) observer->OnProgress(done, images_to_process);
}
{
std::lock_guard lock(harmonic_m);
harmonic_evidence.Add(worker_harmonic);
}
{
std::lock_guard lock(footprint_m);
footprint_offsets.insert(footprint_offsets.end(), worker_offsets.begin(), worker_offsets.end());
}
const auto bragg = analysis.BraggCounts();
bragg_predicted += bragg.predicted;
bragg_starved += bragg.bkg_starved;
bragg_starved_nb += bragg.bkg_starved_by_neighbour;
bragg_fallback += bragg.profile_fallback;
};
if (observer)
observer->OnPhase("Processing images");
Step(fmt::format("{}{}, {} images",
!per_image_analysis ? "Azimuthal integration"
: calibration ? "Image pass"
: experiment_.IsRotationIndexing() ? "Integration" : "Indexing and integration",
geometry_prepass ? " (geometry pre-pass)"
: full && prepass_detector_geometry_ ? " (refined geometry)" : "",
geometry_prepass ? prepass_images : images_to_process));
std::function<void()> worker = per_image_analysis ? std::function<void()>(full_worker)
: std::function<void()>(azint_worker);
// How many workers the loop actually wants. Every one of them submits its own kernels to a card,
// and a card runs out of room to accept them long before it runs out of work: measured on two
// GPUs, the loop's own time falls to sixteen workers and then rises again, so forty-eight is
// slower than sixteen on a 16 Mpx set. Sixteen is where that measurement put the minimum, and a
// later one on a single card agrees - eight leaves the loop waiting on the queue rather than the
// card. The cap is per card because that is what the queue depth belongs to, and it only applies
// when -N was left at its default - an explicit -N is a deliberate instruction and is obeyed,
// which is what a previous attempt at this got wrong.
// Only the full-analysis worker drives a card; the azimuthal one preprocesses and integrates on
// the CPU and wants every thread it can have, so the cap must not reach it.
const int gpus = get_gpu_count();
const int image_workers = (per_image_analysis && config_.nthreads_auto && gpus > 0)
? std::min(config_.nthreads, std::max(16, 16 * gpus))
: config_.nthreads;
if (image_workers < config_.nthreads)
logger.Info("Per-image loop: {} of {} threads ({} GPUs) - past what a card can take queue "
"work from, more workers make it slower, not faster; pass -N to override",
image_workers, config_.nthreads, gpus);
std::vector<std::future<void> > futures;
futures.reserve(image_workers);
const auto image_loop_start = std::chrono::steady_clock::now();
for (int i = 0; i < image_workers; ++i)
futures.push_back(std::async(std::launch::async, worker));
for (auto &f: futures)
f.get();
// Everything the workers posted has to be on disk before the end message is written and the file
// finalized, and a write that failed - out of space, above all - has to surface here rather than
// leave a truncated file behind.
if (writer_queue)
writer_queue->Finish();
// Wall time of the per-image loop alone. The per-stage means below cover only what runs inside it,
// so without this there is nothing to compare them against and no way to see how much of a run is
// spent outside it - on the first-pass indexing and on scaling/merging.
result.image_loop_threads = image_workers;
result.image_loop_time_s = std::chrono::duration<double>(
std::chrono::steady_clock::now() - image_loop_start).count();
result.cancelled = cancelled_;
result.images_processed = finished_count.load();
Note(fmt::format("{} images, {:.0f} per second", result.images_processed,
result.images_processed / std::max(1e-6, result.image_loop_time_s)));
// A pre-pass that stopped at loop_end is a sweep of the frames it read, and that is all the merge and
// the post-refinement are handed. Left in, the empty tail is a stretch of frames every per-frame
// smoothing in the merge extrapolates over: measured, the fulls' log-scale curve, carried straight
// across 900 empty frames, fell further every round until the scale there was zero, and the
// cross-validation of its smoothness then took every scale in the sweep to NaN.
if (geometry_prepass)
indexer->GetIntegrationOutcome().resize(prepass_images);
// The footprint of the reflection rather than of the spot: the pre-scan's widths plus where the
// spots sit against their predictions (SpotFootprint.h). Like the widths it acts only where the
// reflection outgrows the r1 disk, and like them it is held back for the canonical pass.
if (measure_offsets && !cancelled_) {
const SpotFootprint fp = FootprintWithOffsets(*prescan_footprint_, footprint_offsets);
const BraggIntegrationSettings canonical = bragg_adaptive_.value_or(experiment_.GetBraggIntegrationSettings());
bool outgrows = false;
for (size_t b = 0; b < fp.sigma_rad.size(); ++b)
outgrows |= BRAGG_FOOTPRINT_NSIGMA * std::max(fp.sigma_rad[b], fp.sigma_tan[b]) > canonical.GetR1();
std::string table;
for (size_t b = 0; b < fp.sigma_rad.size(); ++b)
table += fmt::format(" {:.0f}:{:.1f}/{:.1f}", (b + 0.5f) * fp.bin_px, fp.sigma_rad[b], fp.sigma_tan[b]);
logger.Info("Reflection footprint: the spot widths plus the offsets of {} indexed spots from their "
"predicted reflections, sigma radial/tangential [px] by distance from the beam [px]:{}{}",
footprint_offsets.size(), table,
outgrows ? "" : fmt::format(" - every reflection fits the r1={:.1f} disk", canonical.GetR1()));
if (outgrows) {
if (!bragg_before_adaptive_)
bragg_before_adaptive_ = experiment_.GetBraggIntegrationSettings();
BraggIntegrationSettings bis = canonical;
bis.Footprint(fp);
bragg_adaptive_ = bis;
}
}
if (full && indexer && experiment_.IsRotationIndexing() && !cancelled_) {
const SupercellProbe probe = indexer->GetSupercellProbe();
if (probe[0][0].n > 0 && probe[0][1].n > 0) {
std::string table;
for (int sh = 0; sh < 3; sh++) {
table += sh == 0 ? "\n 20-5 A:" : sh == 1 ? "\n 5-3 A:" : "\n 20-3 A:";
const auto sums = [&](int c) {
SupercellProbeClass x = sh < 2 ? probe[c][sh] : probe[c][0];
if (sh == 2) x += probe[c][1];
return x;
};
const SupercellProbeFit f0 = FitSupercellProbe(sums(0));
table += fmt::format(" 000 <I/s> {:.2f} a {:.1f} b {:.1f} <p> {:.3f} n {};", f0.mean_i_over_sigma, f0.a,
f0.b, f0.mean_p, sums(0).n);
for (int c = 1; c < 8; c++) {
const SupercellProbeFit f = FitSupercellProbe(sums(c));
table += fmt::format(" {}{}{}: {:.2f} occ {:.2f}% rock {:.2f}+-{:.2f}% flat {:.2f}%;",
(c >> 2) & 1, (c >> 1) & 1, c & 1, f.mean_i_over_sigma,
100.0 * f.mean_i / f0.mean_i, 100.0 * f.b / f0.b, 100.0 * f.b_se / f0.b,
100.0 * f.a / f0.a);
}
}
logger.Info("Supercell probe (the eight parity classes of the doubled primitive cell; <I/sigma>, "
"mean intensity, and the fit I = a + b p against the lattice's own):{}", table);
// The class with the most intensity over 20-3 A goes into the report, with the cell it would
// double the lattice to: the vectors t of the primitive lattice with (h,k,l).t even - twice
// one axis the class is odd on, and each other axis plus that one where the class is odd
// on it too - Niggli-reduced.
const auto band = [&](int c) { SupercellProbeClass x = probe[c][0]; x += probe[c][1]; return x; };
const SupercellProbeFit f0 = FitSupercellProbe(band(0));
const auto prim = indexer->GetSupercellProbePrimitive();
int best = 0;
for (int c = 1; c < 8; c++)
if (best == 0 || FitSupercellProbe(band(c)).mean_i > FitSupercellProbe(band(best)).mean_i)
best = c;
if (prim && f0.mean_i > 0.0 && f0.b != 0.0) {
const SupercellProbeFit f = FitSupercellProbe(band(best));
const int p[3] = {(best >> 2) & 1, (best >> 1) & 1, best & 1};
const Coord v[3] = {prim->Vec0(), prim->Vec1(), prim->Vec2()};
const int i = p[0] ? 0 : p[1] ? 1 : 2;
Coord t[3];
for (int j = 0; j < 3; j++)
t[j] = j == i ? v[i] * 2.0f : v[j] + v[i] * static_cast<float>(p[j]);
ProcessResult::SupercellClass sc;
sc.h = p[0]; sc.k = p[1]; sc.l = p[2];
sc.occupancy_pct = 100.0 * f.mean_i / f0.mean_i;
sc.rock_pct = 100.0 * f.b / f0.b;
sc.rock_se_pct = 100.0 * f.b_se / f0.b;
sc.i_over_sigma = f.mean_i_over_sigma;
sc.doubled_cell = CrystalLattice(t[0], t[1], t[2]).NiggliReduce().GetUnitCell();
result.supercell = sc;
}
}
}
result.harmonic = HarmonicFromEvidence(harmonic_evidence);
result.powder = powder_;
result.powder_excluded_from_indexing = !config_.spot_finding.measured_ring_q_recipA.empty();
if (result.powder_excluded_from_indexing)
result.powder_indexing_d_min_A = config_.spot_finding.high_resolution_limit;
if (result.harmonic.detected)
logger.Warning("Higher-order beam contamination: {:.0f}% of the {} unindexed spots this run "
"could not put on the lattice are on it at one third the wavelength (the same "
"test off a harmonic finds {:.1f}%). They sit at {:.1f} A apparent resolution",
result.harmonic.fraction, result.harmonic.tested,
result.harmonic.null_fraction, result.harmonic.mean_d_A);
// What the stencil cost, as the integrator measured it: the reflections dropped whole for a
// background ring left below six clean pixels, and the rings the NEIGHBOURS' predicted regions
// would starve where the detector alone would not - counting the tails of rocking curves the ring
// is not kept clear of, so it measures the pattern's density. The second is what a wider r1 can
// change, and it is what the two-pass guard reads.
if (bragg_predicted > 0) {
const auto frac = [&](uint64_t n) {
return 100.0 * static_cast<double>(n) / static_cast<double>(bragg_predicted);
};
bkg_starved_fraction_ = static_cast<double>(bragg_starved_nb) / static_cast<double>(bragg_predicted);
logger.Info("Integration at r1={:.1f} r2={:.1f} r3={:.2f}: of {} predicted reflections, "
"{:.3f}% lost their background ring; the neighbouring predictions crowd {:.3f}% of the "
"rings; {:.3f}% profile-fit fallbacks",
experiment_.GetBraggIntegrationSettings().GetR1(),
experiment_.GetBraggIntegrationSettings().GetR2(),
experiment_.GetBraggIntegrationSettings().GetR3(),
bragg_predicted.load(), frac(bragg_starved), frac(bragg_starved_nb),
frac(bragg_fallback));
}
// Every image failing is a total failure, not a run that produced nothing: it used to be reported
// only as per-image log lines while the process still exited 0 with no output file.
if (!cancelled_ && images_to_process > 0 && result.images_processed == 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"No image could be analyzed - see the per-image errors above");
result.mean_processing_time = plots.GetMeanProcessingTime();
result.indexing_rate = plots.GetIndexingRate();
// Powder calibration: the whole run has now been accumulated - the summed (q x azimuth) profile in
// plots, or the pooled spot list - so fit the detector geometry to it. The caller reports it and
// writes the PONI file; nothing here is stored in the _process.h5.
if (calibration && !cancelled_)
FitCalibration(p, calibration_spot_list);
// End message (also written to the file).
// The EXTENT of the per-image datasets (see note_written), not the success count: those two
// differ by exactly the number of frames that failed, and the writer sizes its virtual datasets
// from this.
end_msg.max_image_number = static_cast<uint64_t>(max_written_ordinal.load() + 1);
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.spindle_blind_fraction = plots.GetSpindleBlindFraction();
result.spot_resolution_estimate_A = plots.GetResolutionEstimate();
end_msg.ice_ring_score = plots.GetIceRingScoreArray();
end_msg.ice_ring_score_mean = plots.GetIceRingScore();
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;
// The refined geometry, not the nominal StartMessage values. Two things about these five
// numbers that are not obvious from their names. The beam centre is the PONI - the foot of
// the perpendicular from the sample - so with a non-zero tilt it is NOT where the direct
// beam lands; that point is D*tan(rot)/pixel away and is GetDirectBeam_pxl(). And the tilt
// is not written here only for the record: IndexAndRefine puts it on outcome.experiment,
// which is the geometry Bragg prediction and integration then run at. It is fitted with the
// DISTANCE FROZEN and nothing checks the covariance, so the split between PONI and tilt is
// not reproducible between passes even though their sum - the direct beam - is. Neither
// number is meaningful on its own: report or feed back the whole set, or none of it.
end_msg.refined_beam_center_x = rot->geom.GetBeamX_pxl();
end_msg.refined_beam_center_y = rot->geom.GetBeamY_pxl();
end_msg.refined_poni_rot1 = rot->geom.GetPoniRot1_rad();
end_msg.refined_poni_rot2 = rot->geom.GetPoniRot2_rad();
end_msg.refined_poni_rot3 = rot->geom.GetPoniRot3_rad();
// rot3 is not carried: a rotation about the beam is an exact null of this experiment
// (measured, smallest eigenvalue -5e-16), so the fit never moves it and reporting it would
// suggest it had been determined.
result.refined_detector_tilt_deg = std::array<double, 2>{
rot->geom.GetPoniRot1_rad() * 180.0 / PI,
rot->geom.GetPoniRot2_rad() * 180.0 / PI};
if (rot->tilt_walk && rot->tilt_walk->refused)
result.refused_detector_tilt_deg = std::array<double, 2>{
rot->tilt_walk->tilt_rad[0] * 180.0 / PI, rot->tilt_walk->tilt_rad[1] * 180.0 / PI};
if (rot->axis)
end_msg.refined_rotation_axis = rot->axis->GetAxis();
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.rotation_lattice_type = end_msg.rotation_lattice_type;
// What the equality this class asserts is worth on a REFINED cell. The indexer refines the
// same candidate a second time with nothing held (its pseudo-symmetry guard) and hands that
// cell over whenever the constrained one won; it comes back in the primitive reduced
// basis, and LengthEqualityDeparture puts it in the basis where a and b are the pair this
// class equates. Measured here, decided in RunAllPasses, which is where there is a pass to
// compare it against.
if (rot->unconstrained && ClassImposesLengthEquality(rot->search_result.system)) {
const double split = LengthEqualityDeparture(rot->search_result,
rot->unconstrained->lattice);
result.lattice_length_split = split;
logger.Info("Rotation lattice: {}-centred {} asserts a = b; the free refinement of "
"the same candidate departs from that by {:.2f} %",
rot->search_result.centering,
gemmi::crystal_system_str(rot->search_result.system), 100.0 * split);
}
if (rot->unconstrained && (rot->search_result.system == gemmi::CrystalSystem::Monoclinic
|| rot->search_result.system == gemmi::CrystalSystem::Orthorhombic)) {
const double dep = AngleEqualityDeparture_deg(rot->search_result, rot->unconstrained->lattice);
result.lattice_angle_departure_deg = dep;
logger.Info("Rotation lattice: {}-centred {} holds its angles at 90 deg; the free refinement "
"of the same candidate departs from that by {:.2f} deg",
rot->search_result.centering,
gemmi::crystal_system_str(rot->search_result.system), dep);
}
// The consensus cell below is this lattice's conventional cell, so this is the centring
// it is written in - the pair the volume comparisons in Run() need.
result.consensus_centering = rot->search_result.centering;
}
result.consensus_cell = indexer->GetConsensusUnitCell();
end_msg.unit_cell = result.consensus_cell;
}
ImagePass(p);
// Scaling and merging (full analysis only).
if (full && !cancelled_ && (result.indexing_rate.value_or(0.0f) > 0 || result.rotation_lattice_found)
rugnux/Rugnux.h
+5
View File
@@ -845,6 +845,11 @@ class Rugnux {
// geometry_prepass == true stops after post-refining the detector geometry from the first integration
// and applies it (+ the refined lattice) to experiment_, so the next pass re-integrates with it.
ProcessResult RunPipeline(RugnuxObserver *observer, bool write_output, bool geometry_prepass);
// The phases of RunPipeline that live in files of their own; each reads and writes RunPipeline's
// locals through `p` (RugnuxPipeline.h).
// The per-image loop over the sweep (azimuthal integration, or the full per-image analysis), what is
// read off it once it ends, and the end message.
void ImagePass(PipelineLocals &p);
// Every pass of the run, returning the canonical one. Run() wraps this to time the whole thing.
ProcessResult RunAllPasses(RugnuxObserver *observer);
rugnux/RugnuxImagePass.cpp
+579
View File
@@ -0,0 +1,579 @@
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include <cstdlib>
#include "Rugnux.h"
#include "../image_analysis/structure_refinement/ModelValidation.h"
#include "../writer/WriteModel.h"
#include "../image_analysis/indexing/SpindleCuspLoss.h"
#include "../image_analysis/bragg_integration/SpotWidth.h"
#include "../image_analysis/bragg_integration/SpotFootprint.h"
#include "../image_analysis/bragg_integration/BraggStencil.h" // BRAGG_FOOTPRINT_NSIGMA
#include "../image_analysis/hot_pixels/HotPixels.h"
#include "../image_analysis/SensorAbsorption.h"
#include "DiagnosticOutput.h"
#include <algorithm>
#include <atomic>
#include <chrono>
#include <cmath>
#include <condition_variable>
#include <deque>
#include <fstream>
#include <cstring>
#include <functional>
#include <future>
#include <mutex>
#include <numeric>
#include <map>
#include <set>
#include <sstream>
#include <thread>
#include "../reader/JFJochHDF5Reader.h"
#include "../common/JFJochMath.h"
#include "../common/ParallelFor.h"
#include "../common/Logger.h"
#include "../common/AzimuthalIntegrationMapping.h"
#include "../common/AzimuthalIntegrationProfile.h"
#include "../common/CUDAWrapper.h"
#include "../common/JFJochException.h"
#include "../common/time_utc.h"
#include "../writer/FileWriter.h"
#include "../image_analysis/MXAnalysisWithoutFPGA.h"
#include "../image_analysis/beam_stop/ShadowFinder.h"
#include "../image_analysis/IndexAndRefine.h"
#include "../image_analysis/indexing/AnalyzeIndexing.h"
#include "../image_analysis/indexing/HarmonicContamination.h"
#include "../image_analysis/geom_refinement/BeamCenterFFT.h"
#include "../image_analysis/geom_refinement/BeamCenterFromBackground.h"
#include "../image_analysis/geom_refinement/GeometryRefiner.h"
#include "../image_analysis/indexing/IndexerThreadPool.h"
#include "../image_analysis/spot_finding/ImageSpotFinderCPU.h"
#include "../image_analysis/spot_finding/SpotUtils.h"
#include "gemmi/twin.hpp"
#include "../image_analysis/azint/AzIntEngineCPU.h"
#include "../image_analysis/image_preprocessing/ImagePreprocessorCPU.h"
#include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h"
#ifdef JFJOCH_USE_CUDA
#include "../image_analysis/image_preprocessing/ImagePreprocessorGPU.h"
#include "../image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h"
#include "../image_analysis/spot_finding/AdaptiveSpotFinderGPU.h"
#endif
#include "../image_analysis/scale_merge/Merge.h"
#include "../image_analysis/scale_merge/RfreeFlags.h"
#include "../image_analysis/scale_merge/RotationScaleMerge.h"
#include "../image_analysis/scale_merge/ResolutionCutoff.h"
#include "../image_analysis/scale_merge/ReindexAmbiguity.h"
#include "../image_analysis/scale_merge/ScalingResult.h"
#include "../image_analysis/scale_merge/SearchSpaceGroup.h"
#include "../image_analysis/geom_refinement/PostRefine.h"
#include "../image_analysis/lattice_search/LatticeSearch.h"
#include "../image_analysis/lattice_search/LePageLattice.h"
#include "../image_analysis/scale_merge/AnisotropyAnalysis.h"
#include "../image_analysis/scale_merge/FrenchWilson.h"
#include "../image_analysis/scale_merge/TwinningAnalysis.h"
#include "../image_analysis/scale_merge/TranslationalNCS.h"
#include "../image_analysis/scale_merge/HKLKey.h"
#include "../image_analysis/scale_merge/ScaleOnTheFly.h"
#include "../image_analysis/scale_merge/StillsPartialityRefine.h"
#include "../image_analysis/WriteReflections.h"
#include "../image_analysis/bragg_integration/CalcISigma.h"
#include "../common/Definitions.h"
#include "../common/CorrelationCoefficient.h"
#include <array>
#include <map>
#include "RugnuxPipeline.h"
using namespace rugnux_internal;
void Rugnux::ImagePass(PipelineLocals &p) {
auto &logger = p.logger;
auto &result = p.result;
auto &observer = p.observer;
auto &geometry_prepass = p.geometry_prepass;
auto &full = p.full;
auto &calibration = p.calibration;
auto &per_image_analysis = p.per_image_analysis;
auto &start_image = p.start_image;
auto &end_image = p.end_image;
auto &images_to_process = p.images_to_process;
auto &dataset = p.dataset;
auto &mapping = p.mapping;
auto &plots = p.plots;
auto &writer_queue = p.writer_queue;
auto &indexer = p.indexer;
auto &total_uncompressed_bytes = p.total_uncompressed_bytes;
auto &end_msg = p.end_msg;
auto &prepass_images = p.prepass_images;
auto &loop_end = p.loop_end;
// 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<int> next_ordinal = 0;
std::atomic<int> finished_count = 0;
// Images are written at their own ordinal, so a frame that fails leaves a HOLE rather than
// shifting everything after it up. The extent of the written datasets is therefore the highest
// ordinal reached, not the number that succeeded - counting successes makes the file one image
// short for every failure, and the image it loses is the LAST one, not the one that failed.
std::atomic<int> max_written_ordinal = -1;
auto note_written = [&max_written_ordinal](int ordinal) {
int prev = max_written_ordinal.load(std::memory_order_relaxed);
while (prev < ordinal
&& !max_written_ordinal.compare_exchange_weak(prev, ordinal, std::memory_order_relaxed)) {
}
};
// What the integrator did to the reflections it was handed, summed over the workers as each one
// finishes (BraggIntegrationCounts). Reading an engine's counts synchronises its stream, so it is
// done once per worker at the end rather than per image.
std::atomic<uint64_t> bragg_predicted = 0, bragg_starved = 0, bragg_starved_nb = 0,
bragg_fallback = 0;
// Calibration by spots: the pooled spot list of the whole run. Ring clustering is O(n^2) in it and
// the Hough circle centre O(n^3) in its first few hundred, so each image contributes a fair share of
// a fixed budget instead of everything it found - a powder ring is over-determined either way.
constexpr size_t CALIBRATION_MAX_SPOTS = 20000;
const size_t calibration_spots_per_image =
std::max<size_t>(1, CALIBRATION_MAX_SPOTS / static_cast<size_t>(images_to_process));
std::mutex calibration_spot_m;
std::vector<SpotToSave> calibration_spot_list;
// Higher-order beam contamination, counted on the spots the lattice did not take. Each worker
// keeps its own tally and adds it in when it is done, so the total is a sum of counts and does
// not depend on how the images were shared out.
std::mutex harmonic_m;
HarmonicEvidence harmonic_evidence;
// Where the indexed spots sit against their predicted reflections, measured on the geometry
// pre-pass for the canonical pass's footprint (SpotFootprint.h).
const bool measure_offsets = geometry_prepass && prescan_footprint_.has_value();
std::mutex footprint_m;
std::vector<FootprintOffset> footprint_offsets;
auto azint_worker = [&]() {
std::vector<uint8_t> decompression_buffer;
std::shared_ptr<JFJochReaderRawImage> img;
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;
// Reuse the raw image frame to frame - a fresh 4-byte 16 Mpx frame is 72 MB of
// untouched pages the decoder faults in one by one - except while the writer thread
// still holds the previous one's pixels.
if (!img || img.use_count() > 1)
img = std::make_shared<JFJochReaderRawImage>();
bool read = false;
try {
read = reader_.ReadRawImage(image_idx, *img);
} catch (const std::exception &e) {
if (IsFatalResourceError(e)) throw;
logger.Error("Failed to load image {}: {}", image_idx, e.what());
continue;
}
if (!read) 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) {
if (IsFatalResourceError(e)) throw;
logger.Error("Error integrating image {}: {}", image_idx, e.what());
continue;
}
msg.azint_time_s = std::chrono::duration<float>(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.ice_ring_score = profile.GetIceRingScore(mapping.Settings(),
config_.spot_finding.ice_ring_width_Q_recipA);
msg.run_number = experiment_.GetRunNumber();
msg.run_name = experiment_.GetRunName();
plots.Add(msg, profile);
if (writer_queue) writer_queue->Post(msg, img);
note_written(ordinal);
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
// Before the analysis engine, which page-locks these bytes for its uploads and unregisters
// them when it is destroyed: it must not outlive the buffer it was given.
std::shared_ptr<JFJochReaderRawImage> img;
MXAnalysisWithoutFPGA analysis(experiment_, mapping, pixel_mask_, *indexer,
/*enable_fused_adaptive_gpu=*/true);
AzimuthalIntegrationProfile profile(mapping);
HarmonicEvidence worker_harmonic;
std::vector<FootprintOffset> worker_offsets;
while (!cancelled_) {
const int ordinal = next_ordinal.fetch_add(1);
const int image_idx = start_image + ordinal * config_.stride;
if (image_idx >= loop_end) break;
// Reuse the raw image frame to frame - a fresh 4-byte 16 Mpx frame is 72 MB of
// untouched pages the decoder faults in one by one - except while the writer thread
// still holds the previous one's pixels.
if (!img || img.use_count() > 1)
img = std::make_shared<JFJochReaderRawImage>();
bool read = false;
try {
read = reader_.ReadRawImage(image_idx, *img);
} catch (const std::exception &e) {
if (IsFatalResourceError(e)) throw;
logger.Error("Failed to load image {}: {}", image_idx, e.what());
continue;
}
if (!read) 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) {
if (IsFatalResourceError(e)) throw;
logger.Error("Error analyzing image {}: {}", image_idx, e.what());
continue;
}
msg.processing_time_s = std::chrono::duration<float>(std::chrono::steady_clock::now() - t0).count();
msg.run_number = experiment_.GetRunNumber();
msg.run_name = experiment_.GetRunName();
if (calibration) {
std::lock_guard lock(calibration_spot_m);
for (size_t i = 0; i < msg.spots.size() && i < calibration_spots_per_image; ++i)
calibration_spot_list.push_back(msg.spots[i]);
}
AddHarmonicEvidence(msg.spots, experiment_.GetDiffractionGeometry(), worker_harmonic);
if (measure_offsets)
MeasureFootprintOffsets(msg.spots, msg.reflections, experiment_.GetBeamX_pxl(),
experiment_.GetBeamY_pxl(), worker_offsets);
plots.Add(msg, profile);
if (writer_queue) writer_queue->Post(msg, img);
note_written(ordinal);
if (observer) observer->OnImageProcessed(msg);
const int done = finished_count.fetch_add(1) + 1;
if (observer) observer->OnProgress(done, images_to_process);
}
{
std::lock_guard lock(harmonic_m);
harmonic_evidence.Add(worker_harmonic);
}
{
std::lock_guard lock(footprint_m);
footprint_offsets.insert(footprint_offsets.end(), worker_offsets.begin(), worker_offsets.end());
}
const auto bragg = analysis.BraggCounts();
bragg_predicted += bragg.predicted;
bragg_starved += bragg.bkg_starved;
bragg_starved_nb += bragg.bkg_starved_by_neighbour;
bragg_fallback += bragg.profile_fallback;
};
if (observer)
observer->OnPhase("Processing images");
Step(fmt::format("{}{}, {} images",
!per_image_analysis ? "Azimuthal integration"
: calibration ? "Image pass"
: experiment_.IsRotationIndexing() ? "Integration" : "Indexing and integration",
geometry_prepass ? " (geometry pre-pass)"
: full && prepass_detector_geometry_ ? " (refined geometry)" : "",
geometry_prepass ? prepass_images : images_to_process));
std::function<void()> worker = per_image_analysis ? std::function<void()>(full_worker)
: std::function<void()>(azint_worker);
// How many workers the loop actually wants. Every one of them submits its own kernels to a card,
// and a card runs out of room to accept them long before it runs out of work: measured on two
// GPUs, the loop's own time falls to sixteen workers and then rises again, so forty-eight is
// slower than sixteen on a 16 Mpx set. Sixteen is where that measurement put the minimum, and a
// later one on a single card agrees - eight leaves the loop waiting on the queue rather than the
// card. The cap is per card because that is what the queue depth belongs to, and it only applies
// when -N was left at its default - an explicit -N is a deliberate instruction and is obeyed,
// which is what a previous attempt at this got wrong.
// Only the full-analysis worker drives a card; the azimuthal one preprocesses and integrates on
// the CPU and wants every thread it can have, so the cap must not reach it.
const int gpus = get_gpu_count();
const int image_workers = (per_image_analysis && config_.nthreads_auto && gpus > 0)
? std::min(config_.nthreads, std::max(16, 16 * gpus))
: config_.nthreads;
if (image_workers < config_.nthreads)
logger.Info("Per-image loop: {} of {} threads ({} GPUs) - past what a card can take queue "
"work from, more workers make it slower, not faster; pass -N to override",
image_workers, config_.nthreads, gpus);
std::vector<std::future<void> > futures;
futures.reserve(image_workers);
const auto image_loop_start = std::chrono::steady_clock::now();
for (int i = 0; i < image_workers; ++i)
futures.push_back(std::async(std::launch::async, worker));
for (auto &f: futures)
f.get();
// Everything the workers posted has to be on disk before the end message is written and the file
// finalized, and a write that failed - out of space, above all - has to surface here rather than
// leave a truncated file behind.
if (writer_queue)
writer_queue->Finish();
// Wall time of the per-image loop alone. The per-stage means below cover only what runs inside it,
// so without this there is nothing to compare them against and no way to see how much of a run is
// spent outside it - on the first-pass indexing and on scaling/merging.
result.image_loop_threads = image_workers;
result.image_loop_time_s = std::chrono::duration<double>(
std::chrono::steady_clock::now() - image_loop_start).count();
result.cancelled = cancelled_;
result.images_processed = finished_count.load();
Note(fmt::format("{} images, {:.0f} per second", result.images_processed,
result.images_processed / std::max(1e-6, result.image_loop_time_s)));
// A pre-pass that stopped at loop_end is a sweep of the frames it read, and that is all the merge and
// the post-refinement are handed. Left in, the empty tail is a stretch of frames every per-frame
// smoothing in the merge extrapolates over: measured, the fulls' log-scale curve, carried straight
// across 900 empty frames, fell further every round until the scale there was zero, and the
// cross-validation of its smoothness then took every scale in the sweep to NaN.
if (geometry_prepass)
indexer->GetIntegrationOutcome().resize(prepass_images);
// The footprint of the reflection rather than of the spot: the pre-scan's widths plus where the
// spots sit against their predictions (SpotFootprint.h). Like the widths it acts only where the
// reflection outgrows the r1 disk, and like them it is held back for the canonical pass.
if (measure_offsets && !cancelled_) {
const SpotFootprint fp = FootprintWithOffsets(*prescan_footprint_, footprint_offsets);
const BraggIntegrationSettings canonical = bragg_adaptive_.value_or(experiment_.GetBraggIntegrationSettings());
bool outgrows = false;
for (size_t b = 0; b < fp.sigma_rad.size(); ++b)
outgrows |= BRAGG_FOOTPRINT_NSIGMA * std::max(fp.sigma_rad[b], fp.sigma_tan[b]) > canonical.GetR1();
std::string table;
for (size_t b = 0; b < fp.sigma_rad.size(); ++b)
table += fmt::format(" {:.0f}:{:.1f}/{:.1f}", (b + 0.5f) * fp.bin_px, fp.sigma_rad[b], fp.sigma_tan[b]);
logger.Info("Reflection footprint: the spot widths plus the offsets of {} indexed spots from their "
"predicted reflections, sigma radial/tangential [px] by distance from the beam [px]:{}{}",
footprint_offsets.size(), table,
outgrows ? "" : fmt::format(" - every reflection fits the r1={:.1f} disk", canonical.GetR1()));
if (outgrows) {
if (!bragg_before_adaptive_)
bragg_before_adaptive_ = experiment_.GetBraggIntegrationSettings();
BraggIntegrationSettings bis = canonical;
bis.Footprint(fp);
bragg_adaptive_ = bis;
}
}
if (full && indexer && experiment_.IsRotationIndexing() && !cancelled_) {
const SupercellProbe probe = indexer->GetSupercellProbe();
if (probe[0][0].n > 0 && probe[0][1].n > 0) {
std::string table;
for (int sh = 0; sh < 3; sh++) {
table += sh == 0 ? "\n 20-5 A:" : sh == 1 ? "\n 5-3 A:" : "\n 20-3 A:";
const auto sums = [&](int c) {
SupercellProbeClass x = sh < 2 ? probe[c][sh] : probe[c][0];
if (sh == 2) x += probe[c][1];
return x;
};
const SupercellProbeFit f0 = FitSupercellProbe(sums(0));
table += fmt::format(" 000 <I/s> {:.2f} a {:.1f} b {:.1f} <p> {:.3f} n {};", f0.mean_i_over_sigma, f0.a,
f0.b, f0.mean_p, sums(0).n);
for (int c = 1; c < 8; c++) {
const SupercellProbeFit f = FitSupercellProbe(sums(c));
table += fmt::format(" {}{}{}: {:.2f} occ {:.2f}% rock {:.2f}+-{:.2f}% flat {:.2f}%;",
(c >> 2) & 1, (c >> 1) & 1, c & 1, f.mean_i_over_sigma,
100.0 * f.mean_i / f0.mean_i, 100.0 * f.b / f0.b, 100.0 * f.b_se / f0.b,
100.0 * f.a / f0.a);
}
}
logger.Info("Supercell probe (the eight parity classes of the doubled primitive cell; <I/sigma>, "
"mean intensity, and the fit I = a + b p against the lattice's own):{}", table);
// The class with the most intensity over 20-3 A goes into the report, with the cell it would
// double the lattice to: the vectors t of the primitive lattice with (h,k,l).t even - twice
// one axis the class is odd on, and each other axis plus that one where the class is odd
// on it too - Niggli-reduced.
const auto band = [&](int c) { SupercellProbeClass x = probe[c][0]; x += probe[c][1]; return x; };
const SupercellProbeFit f0 = FitSupercellProbe(band(0));
const auto prim = indexer->GetSupercellProbePrimitive();
int best = 0;
for (int c = 1; c < 8; c++)
if (best == 0 || FitSupercellProbe(band(c)).mean_i > FitSupercellProbe(band(best)).mean_i)
best = c;
if (prim && f0.mean_i > 0.0 && f0.b != 0.0) {
const SupercellProbeFit f = FitSupercellProbe(band(best));
const int p[3] = {(best >> 2) & 1, (best >> 1) & 1, best & 1};
const Coord v[3] = {prim->Vec0(), prim->Vec1(), prim->Vec2()};
const int i = p[0] ? 0 : p[1] ? 1 : 2;
Coord t[3];
for (int j = 0; j < 3; j++)
t[j] = j == i ? v[i] * 2.0f : v[j] + v[i] * static_cast<float>(p[j]);
ProcessResult::SupercellClass sc;
sc.h = p[0]; sc.k = p[1]; sc.l = p[2];
sc.occupancy_pct = 100.0 * f.mean_i / f0.mean_i;
sc.rock_pct = 100.0 * f.b / f0.b;
sc.rock_se_pct = 100.0 * f.b_se / f0.b;
sc.i_over_sigma = f.mean_i_over_sigma;
sc.doubled_cell = CrystalLattice(t[0], t[1], t[2]).NiggliReduce().GetUnitCell();
result.supercell = sc;
}
}
}
result.harmonic = HarmonicFromEvidence(harmonic_evidence);
result.powder = powder_;
result.powder_excluded_from_indexing = !config_.spot_finding.measured_ring_q_recipA.empty();
if (result.powder_excluded_from_indexing)
result.powder_indexing_d_min_A = config_.spot_finding.high_resolution_limit;
if (result.harmonic.detected)
logger.Warning("Higher-order beam contamination: {:.0f}% of the {} unindexed spots this run "
"could not put on the lattice are on it at one third the wavelength (the same "
"test off a harmonic finds {:.1f}%). They sit at {:.1f} A apparent resolution",
result.harmonic.fraction, result.harmonic.tested,
result.harmonic.null_fraction, result.harmonic.mean_d_A);
// What the stencil cost, as the integrator measured it: the reflections dropped whole for a
// background ring left below six clean pixels, and the rings the NEIGHBOURS' predicted regions
// would starve where the detector alone would not - counting the tails of rocking curves the ring
// is not kept clear of, so it measures the pattern's density. The second is what a wider r1 can
// change, and it is what the two-pass guard reads.
if (bragg_predicted > 0) {
const auto frac = [&](uint64_t n) {
return 100.0 * static_cast<double>(n) / static_cast<double>(bragg_predicted);
};
bkg_starved_fraction_ = static_cast<double>(bragg_starved_nb) / static_cast<double>(bragg_predicted);
logger.Info("Integration at r1={:.1f} r2={:.1f} r3={:.2f}: of {} predicted reflections, "
"{:.3f}% lost their background ring; the neighbouring predictions crowd {:.3f}% of the "
"rings; {:.3f}% profile-fit fallbacks",
experiment_.GetBraggIntegrationSettings().GetR1(),
experiment_.GetBraggIntegrationSettings().GetR2(),
experiment_.GetBraggIntegrationSettings().GetR3(),
bragg_predicted.load(), frac(bragg_starved), frac(bragg_starved_nb),
frac(bragg_fallback));
}
// Every image failing is a total failure, not a run that produced nothing: it used to be reported
// only as per-image log lines while the process still exited 0 with no output file.
if (!cancelled_ && images_to_process > 0 && result.images_processed == 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"No image could be analyzed - see the per-image errors above");
result.mean_processing_time = plots.GetMeanProcessingTime();
result.indexing_rate = plots.GetIndexingRate();
// Powder calibration: the whole run has now been accumulated - the summed (q x azimuth) profile in
// plots, or the pooled spot list - so fit the detector geometry to it. The caller reports it and
// writes the PONI file; nothing here is stored in the _process.h5.
if (calibration && !cancelled_)
FitCalibration(p, calibration_spot_list);
// End message (also written to the file).
// The EXTENT of the per-image datasets (see note_written), not the success count: those two
// differ by exactly the number of frames that failed, and the writer sizes its virtual datasets
// from this.
end_msg.max_image_number = static_cast<uint64_t>(max_written_ordinal.load() + 1);
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.spindle_blind_fraction = plots.GetSpindleBlindFraction();
result.spot_resolution_estimate_A = plots.GetResolutionEstimate();
end_msg.ice_ring_score = plots.GetIceRingScoreArray();
end_msg.ice_ring_score_mean = plots.GetIceRingScore();
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;
// The refined geometry, not the nominal StartMessage values. Two things about these five
// numbers that are not obvious from their names. The beam centre is the PONI - the foot of
// the perpendicular from the sample - so with a non-zero tilt it is NOT where the direct
// beam lands; that point is D*tan(rot)/pixel away and is GetDirectBeam_pxl(). And the tilt
// is not written here only for the record: IndexAndRefine puts it on outcome.experiment,
// which is the geometry Bragg prediction and integration then run at. It is fitted with the
// DISTANCE FROZEN and nothing checks the covariance, so the split between PONI and tilt is
// not reproducible between passes even though their sum - the direct beam - is. Neither
// number is meaningful on its own: report or feed back the whole set, or none of it.
end_msg.refined_beam_center_x = rot->geom.GetBeamX_pxl();
end_msg.refined_beam_center_y = rot->geom.GetBeamY_pxl();
end_msg.refined_poni_rot1 = rot->geom.GetPoniRot1_rad();
end_msg.refined_poni_rot2 = rot->geom.GetPoniRot2_rad();
end_msg.refined_poni_rot3 = rot->geom.GetPoniRot3_rad();
// rot3 is not carried: a rotation about the beam is an exact null of this experiment
// (measured, smallest eigenvalue -5e-16), so the fit never moves it and reporting it would
// suggest it had been determined.
result.refined_detector_tilt_deg = std::array<double, 2>{
rot->geom.GetPoniRot1_rad() * 180.0 / PI,
rot->geom.GetPoniRot2_rad() * 180.0 / PI};
if (rot->tilt_walk && rot->tilt_walk->refused)
result.refused_detector_tilt_deg = std::array<double, 2>{
rot->tilt_walk->tilt_rad[0] * 180.0 / PI, rot->tilt_walk->tilt_rad[1] * 180.0 / PI};
if (rot->axis)
end_msg.refined_rotation_axis = rot->axis->GetAxis();
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.rotation_lattice_type = end_msg.rotation_lattice_type;
// What the equality this class asserts is worth on a REFINED cell. The indexer refines the
// same candidate a second time with nothing held (its pseudo-symmetry guard) and hands that
// cell over whenever the constrained one won; it comes back in the primitive reduced
// basis, and LengthEqualityDeparture puts it in the basis where a and b are the pair this
// class equates. Measured here, decided in RunAllPasses, which is where there is a pass to
// compare it against.
if (rot->unconstrained && ClassImposesLengthEquality(rot->search_result.system)) {
const double split = LengthEqualityDeparture(rot->search_result,
rot->unconstrained->lattice);
result.lattice_length_split = split;
logger.Info("Rotation lattice: {}-centred {} asserts a = b; the free refinement of "
"the same candidate departs from that by {:.2f} %",
rot->search_result.centering,
gemmi::crystal_system_str(rot->search_result.system), 100.0 * split);
}
if (rot->unconstrained && (rot->search_result.system == gemmi::CrystalSystem::Monoclinic
|| rot->search_result.system == gemmi::CrystalSystem::Orthorhombic)) {
const double dep = AngleEqualityDeparture_deg(rot->search_result, rot->unconstrained->lattice);
result.lattice_angle_departure_deg = dep;
logger.Info("Rotation lattice: {}-centred {} holds its angles at 90 deg; the free refinement "
"of the same candidate departs from that by {:.2f} deg",
rot->search_result.centering,
gemmi::crystal_system_str(rot->search_result.system), dep);
}
// The consensus cell below is this lattice's conventional cell, so this is the centring
// it is written in - the pair the volume comparisons in Run() need.
result.consensus_centering = rot->search_result.centering;
}
result.consensus_cell = indexer->GetConsensusUnitCell();
end_msg.unit_cell = result.consensus_cell;
}
}
rugnux/RugnuxPipeline.h
+3
View File
@@ -140,6 +140,7 @@ struct PipelineLocals {
ProcessResult &result;
RugnuxObserver *observer;
const bool geometry_prepass;
const bool full;
const bool calibration;
const bool calibration_spots;
const bool per_image_analysis;
@@ -147,6 +148,8 @@ struct PipelineLocals {
const int start_image;
const int end_image;
const int images_to_process;
const int prepass_images;
const int loop_end;
const std::shared_ptr<const JFJochReaderDataset> &dataset;
AzimuthalIntegrationMapping &mapping;
JFJochReceiverPlots &plots;