leonarski_fandClaude Opus 5.5
a509fd8022
The defective-pixel finder ran only on scanning rotation data. On stills, a pixel that holds the detector's error value on every frame but is missing from the file's static mask is read by the integrator as a per-frame unreadable pixel, so every reflection whose signal disk touches it is flagged overloaded and the merge drops it. On a SwissFEL JUNGFRAU stills run, 14k such pixels made 3.9% of all integrated reflections "overloaded" (mostly finer than 1.6 A, where they are scarcest) and raised a false OVERLOADS warning; with them masked, 0.05% remain flagged and the profile fit keeps the reflections whose disk stays at least MINPK readable. On stills only the error-value and dead rules apply: the hot-pixel rule counts a reflection's stay on a pixel in rotation, so spacing is passed as 0 there, which HotPixelFinder already treats as "no hot-pixel test". Rotation is unchanged (same gate, same oscillation and spacing). Result-changing for stills; merged-statistics A/B on full runs and the battery ci check still to do. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
1212 lines
73 KiB
C++
1212 lines
73 KiB
C++
// 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/BackgroundGain.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/ThreadAffinity.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;
|
|
|
|
namespace {
|
|
// A detector-plane displacement resolved ACROSS the spindle - the only component the tolerance
|
|
// above is about. Along the spindle the error translates the derotated cloud rigidly and the FFT
|
|
// amplitude cannot see it at all; across it the peaks smear and the amplitude falls. (The lattice
|
|
// fit that follows the FFT is not translation-invariant, so "along the spindle" is free for the
|
|
// transform and not for the run - which is why this is reported and not used to decide anything.)
|
|
// The detector direction wanted is the one whose lab image is perpendicular to the spindle n:
|
|
// (ux e1 + uy e2).n = 0, so u is proportional to (-(e2.n), e1.n).
|
|
std::optional<float> AcrossSpindle_pxl(const DiffractionExperiment &experiment, float dx, float dy) {
|
|
const auto goniometer = experiment.GetGoniometer();
|
|
if (!goniometer || !goniometer->IsScanning())
|
|
return {};
|
|
const auto geom = experiment.GetDiffractionGeometry();
|
|
const float beam_x = geom.GetBeamX_pxl(), beam_y = geom.GetBeamY_pxl();
|
|
const Coord n = goniometer->GetAxis().Normalize();
|
|
const Coord origin = geom.LabCoord(beam_x, beam_y);
|
|
const Coord e1 = geom.LabCoord(beam_x + 1.0f, beam_y) - origin;
|
|
const Coord e2 = geom.LabCoord(beam_x, beam_y + 1.0f) - origin;
|
|
const float ux = -(e2 * n), uy = e1 * n;
|
|
const float length = std::hypot(ux, uy);
|
|
if (length < 1e-3f)
|
|
return {};
|
|
return std::abs(dx * ux + dy * uy) / length;
|
|
}
|
|
|
|
|
|
// Images left out at each end of the sweep. Shutter synchronisation disturbs the first and last
|
|
// frame or two of a run, and a partly exposed frame is dark over the whole detector - which is
|
|
// exactly the kind of error a projection cannot average away and a background comparison reads
|
|
// as signal. The margin is capped at a tenth of the sweep so a short run still has a sample.
|
|
constexpr int PRESCAN_END_MARGIN_IMAGES = 5;
|
|
|
|
// Workers reading the pre-scan sample. Each holds detector-sized buffers of its own, and pages of
|
|
// fresh memory are slow to fault in when many threads do it at once, so this is capped well below
|
|
// the worker count of the run proper. The beam-stop projection is memory-bound rather than
|
|
// compute-bound, so a handful of workers already saturates it.
|
|
constexpr size_t PRESCAN_MAX_WORKERS = 8;
|
|
|
|
// The spot width is measured on a GROWING share of the pre-scan sample: every eighth frame of
|
|
// it, then every fourth, then every second, then all of it - an eighth, a quarter, a half and
|
|
// the whole. Each stride is a multiple of the next, so a tier's sample really does contain the
|
|
// one before it, and the measurement stops at the first tier whose answer has settled
|
|
// (spot_width::WidthSettled, in the loop below). A crystal whose frames carry plenty of isolated
|
|
// strong spots is characterised from a quarter of the sample; a sparse one still gets all of it.
|
|
//
|
|
// The tiers add PRESCAN_MAX_WORKERS-sized steps as well: on a 60 frame sample they are 8, 15, 30
|
|
// and 60 frames, i.e. 8 + 7 + 15 + 30 read in the four passes, so no pass leaves most of the
|
|
// workers idle waiting for the others to finish theirs. A ladder that starts smaller measures
|
|
// fewer frames but costs more wall clock, because a pass of five frames occupies eight workers
|
|
// for as long as a pass of eight does.
|
|
//
|
|
// Only the width work is tiered. Every frame of the sample is still READ once - the frames are
|
|
// simply visited in tier order - and the beam-stop projection still sees all of them, so the
|
|
// mask and the beam centre are what they were. What the tiers cut is the per-frame work the
|
|
// width needs on top of the read: the host decompression, the preprocessing and the spot
|
|
// finding, which together are two thirds of the pre-scan's cycles once the projection itself is
|
|
// decoding on the GPU. The width measurement proper is under 1 % of them.
|
|
//
|
|
// Measured over the rotation battery: the radius chosen is the whole-sample answer on 38 of 38
|
|
// crystals, from an average of 32 frames of the 60 instead of all of them.
|
|
constexpr std::array<int, 4> WIDTH_TIER_STRIDE = {8, 4, 2, 1};
|
|
|
|
// The sample is read once and kept for everything that reads it after the projection - the spot
|
|
// measurement and the defective pixels used to read it again each, and a CBF, SMV or marCCD frame
|
|
// is decoded on the host every time it is read. Frames are kept as the reader gives them,
|
|
// compressed where the file is: a few MB a frame on a 16 Mpx EIGER, 4 bytes a pixel where the
|
|
// reader decodes. Frames are kept while the kept ones take less than this - so it is passed by at
|
|
// most a frame per worker - and the rest are read again where they are wanted.
|
|
constexpr size_t PRESCAN_KEPT_BYTES = size_t{2} << 30;
|
|
|
|
// fn(t) for t < n on threads of their own, thread t on card t % gpus as an image-loop worker is: an
|
|
// engine built on a worker's card must be used, and freed, with that card current. Which card a
|
|
// frame goes to changes nothing it gives: the spot finders are exact, and the projection and the
|
|
// defective-pixel sums are integers, kept per card and added together when they are read.
|
|
template <class Fn>
|
|
void OnWorkers(size_t n, Fn &&fn) {
|
|
std::vector<std::future<void>> futures;
|
|
futures.reserve(n);
|
|
for (size_t t = 0; t < n; t++)
|
|
futures.emplace_back(std::async(std::launch::async, [&fn, t] {
|
|
if (const int32_t gpus = get_gpu_count(); gpus > 0)
|
|
pin_gpu(static_cast<int32_t>(t) % gpus);
|
|
fn(t);
|
|
}));
|
|
for (auto &f : futures) f.get();
|
|
}
|
|
}
|
|
|
|
void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, RugnuxObserver *observer) {
|
|
Logger logger("Rugnux");
|
|
|
|
// Two-pass rotation runs this twice. The shadow does not move, and re-detecting would find
|
|
// nothing (its pixels are masked by now) and so clear the mask the first pass established.
|
|
const auto ¤t = pixel_mask_.GetMask();
|
|
const bool have_shadow = std::any_of(current.begin(), current.end(),
|
|
[](uint32_t v) { return (v & (1u << PixelMask::BeamStopPixelBit)) != 0; });
|
|
const bool want_shadow = config_.detect_beam_stop.has_value() && !have_shadow;
|
|
// The spot width is a property of the crystal, not of the geometry, so the second pass reuses the
|
|
// radius the first one measured rather than re-reading the frames.
|
|
const bool want_width = config_.adaptive_integration_radius && !spot_width_measured_;
|
|
// Defective pixels, for the same reason on the same pass: each pixel is compared against its
|
|
// resolution ring. Only under rotation, where a reflection moves off a pixel and a defect does not -
|
|
// except for the pixels holding the detector's error value and the dead ones, which need no
|
|
// rotation and are measured on stills too (MaskDefectivePixels).
|
|
const auto goniometer = experiment_.GetGoniometer();
|
|
const bool want_hot = !hot_pixels_measured_
|
|
&& (!experiment_.IsRotationIndexing() || (goniometer && goniometer->IsScanning()));
|
|
// The background noise per unit level of a CCD, read off the same frames (BackgroundGain.h).
|
|
const bool want_gain = config_.counts_in_adu && !background_gain_measured_;
|
|
if (!want_shadow && !want_width && !want_hot && !want_gain)
|
|
return;
|
|
|
|
// The sample is taken from the frames given - on a rotation sweep its first part (PrepassEnd) -
|
|
// both ends left out: an even spread over [0, images_to_process) always picks the very first and
|
|
// the very last image, and those are the ones a shutter out of step spoils - the last carries
|
|
// under 90 % of the sample's median counts on 12 of 38 regression sweeps and is essentially
|
|
// unexposed on 5.
|
|
const int margin = std::min(PRESCAN_END_MARGIN_IMAGES, images_to_process / 10);
|
|
auto shadow_sample = select_equally_spaced_image_ordinals(images_to_process - 2 * margin,
|
|
frame_count);
|
|
for (int &ordinal : shadow_sample)
|
|
ordinal += margin;
|
|
const std::set<int> shadow_set(shadow_sample.begin(), shadow_sample.end());
|
|
if (shadow_set.empty())
|
|
return;
|
|
|
|
// Say what this is reading for. It runs before the first image of the run proper is processed,
|
|
// so without this the job sits at its initial status for the whole pre-scan.
|
|
if (observer)
|
|
observer->OnPhase(want_shadow ? "Beam stop" : "Spot width");
|
|
Step("Pre-scan");
|
|
|
|
// What the beam-stop detection reads: the frames, the geometry the rings are drawn in and the mask.
|
|
// A detection that finds no shadow leaves the mask as it was, so every copy of the run that
|
|
// pre-scans after it - a probe, a re-check - would detect again on the same inputs; it takes the
|
|
// answer kept under them instead (BeamStopMemo).
|
|
std::vector<double> shadow_key;
|
|
std::optional<BeamStopMemo::Entry> shadow_memo;
|
|
{
|
|
const auto &g = experiment_.GetDiffractionGeometry();
|
|
const uint64_t mask_checksum = pixel_mask_.GetBinaryMaskChecksum();
|
|
shadow_key = {static_cast<double>(start_image), static_cast<double>(config_.stride),
|
|
static_cast<double>(experiment_.GetXPixelsNumConv()),
|
|
static_cast<double>(experiment_.GetYPixelsNumConv()),
|
|
static_cast<double>(experiment_.IsGeometryTransformed()),
|
|
static_cast<double>(mask_checksum >> 32), static_cast<double>(mask_checksum & 0xffffffffu),
|
|
g.GetBeamX_pxl(), g.GetBeamY_pxl(), g.GetDetectorDistance_mm(), g.GetPixelSize_mm(),
|
|
g.GetWavelength_A(), g.GetPoniRot1_rad(), g.GetPoniRot2_rad(), g.GetPoniRot3_rad(),
|
|
static_cast<double>(experiment_.GetPolarizationFactor().has_value()),
|
|
experiment_.GetPolarizationFactor().value_or(0.0f)};
|
|
for (const float v : g.GetDetectorMatrix().arr())
|
|
shadow_key.push_back(v);
|
|
shadow_key.insert(shadow_key.end(), shadow_sample.begin(), shadow_sample.end());
|
|
}
|
|
if (want_shadow) {
|
|
std::lock_guard lock(beam_stop_memo_->m);
|
|
for (const auto &e : beam_stop_memo_->entries)
|
|
if (e.key == shadow_key)
|
|
shadow_memo = e;
|
|
}
|
|
const bool find_shadow = want_shadow && !shadow_memo;
|
|
|
|
// The spot measurement, the defective pixels and the background gain are kept the same way, under
|
|
// the same key and everything else they read: the settings the spots are found with, the
|
|
// azimuthal binning the adaptive finder thresholds on, and whether the shadow is detected first
|
|
// (the defective pixels are compared with it masked). A copy of the run that pre-scans again on
|
|
// the same inputs takes them from here and, where nothing else wants the frames, reads none.
|
|
std::vector<double> measure_key = shadow_key;
|
|
{
|
|
const auto &f = config_.spot_finding;
|
|
const auto &a = experiment_.GetAzimuthalIntegrationSettings();
|
|
measure_key.insert(measure_key.end(),
|
|
{static_cast<double>(want_shadow), static_cast<double>(experiment_.GetSaturationLimit()),
|
|
goniometer ? goniometer->GetIncrement_deg() : 0.0,
|
|
a.GetLowQ_recipA(), a.GetHighQ_recipA(), a.GetQSpacing_recipA(),
|
|
static_cast<double>(a.GetAzimuthalBinCount()), static_cast<double>(a.IsSolidAngleCorrection()),
|
|
static_cast<double>(a.IsPolarizationCorrection()),
|
|
static_cast<double>(f.enable), f.signal_to_noise_threshold,
|
|
static_cast<double>(f.photon_count_threshold),
|
|
static_cast<double>(f.min_pix_per_spot.value_or(-1)), static_cast<double>(f.max_pix_per_spot),
|
|
f.high_resolution_limit.value_or(-1.0f), f.low_resolution_limit.value_or(-1.0f),
|
|
f.cutoff_spot_count_low_res, f.high_res_gap_Q_recipA.value_or(-1.0f),
|
|
f.ice_ring_width_Q_recipA, static_cast<double>(f.adaptive_threshold),
|
|
f.false_pixels_per_frame, static_cast<double>(experiment_.IsDetectIceRings()),
|
|
static_cast<double>(f.measured_ring_q_recipA.size())});
|
|
measure_key.insert(measure_key.end(), f.measured_ring_q_recipA.begin(), f.measured_ring_q_recipA.end());
|
|
}
|
|
std::optional<PreScanMemo::Entry> measured;
|
|
{
|
|
std::lock_guard lock(prescan_memo_->m);
|
|
for (const auto &e : prescan_memo_->entries)
|
|
if (e.key == measure_key)
|
|
measured = e;
|
|
}
|
|
const bool measure_width = want_width && !(measured && measured->spots);
|
|
const bool measure_hot = want_hot && !(measured && measured->defects);
|
|
const bool measure_gain = want_gain && !(measured && measured->gain);
|
|
|
|
std::optional<ShadowFinder> finder;
|
|
if (find_shadow)
|
|
finder.emplace(experiment_, pixel_mask_);
|
|
|
|
// The spot width is read in the same pass as the projection rather than in a pre-pass of its own.
|
|
// The image is preprocessed a second time for it: the shadow is accumulated from raw counts and
|
|
// the spot finder works on the converted image.
|
|
const bool want_spots = measure_width;
|
|
// The spots are measured in the background while the shadow, the defective pixels and the beam
|
|
// centre are found (see below), and those write the mask. So the spot measurement reads its own
|
|
// copies of the experiment and of the mask as they are now - which is what it always read, since
|
|
// it used to finish before any of them started.
|
|
const DiffractionExperiment prescan_x = experiment_;
|
|
const PixelMask prescan_mask = pixel_mask_;
|
|
// Built by the spot measurement, the only thing that reads it, so that it is built beside the
|
|
// shadow rather than ahead of it.
|
|
std::unique_ptr<AzimuthalIntegrationMapping> prescan_mapping;
|
|
|
|
// Everything reading one frame needs that cannot be shared between threads: its own scratch for
|
|
// each of the two decompressions, its own preprocessed image and its own spot finder.
|
|
struct PreScanWorker {
|
|
std::vector<uint8_t> decompression_buffer;
|
|
JFJochReaderRawImage raw_image;
|
|
std::unique_ptr<ImagePreprocessor> preprocessor;
|
|
std::unique_ptr<ImagePreprocessorBuffer> preprocessed;
|
|
std::unique_ptr<ImageSpotFinder> spot_finder;
|
|
};
|
|
const auto make_worker = [&] {
|
|
PreScanWorker w;
|
|
#ifdef JFJOCH_USE_CUDA
|
|
// On the card where there is one, as in the image loops: the frame is decoded and preprocessed
|
|
// there and the spots are found and extracted there. The device finders give the host's spot
|
|
// list to the bit - integer ring sums, the same connected components in the same order
|
|
// (AdaptiveSpotFinderGPU, SpotExtractorGPU) - so this is a choice of where, not of what. Of the
|
|
// preprocessed image only the squares about the spots come back (SPOT_TILE_HALF).
|
|
if (want_spots && get_gpu_count() > 0) {
|
|
auto stream = std::make_shared<CudaStream>();
|
|
w.preprocessor = std::make_unique<ImagePreprocessorGPU>(prescan_x, prescan_mask, stream,
|
|
/*copy_image_to_host=*/false);
|
|
w.preprocessed = std::make_unique<ImagePreprocessorBufferGPU>(prescan_x.GetPixelsNum(),
|
|
/*host_mirror=*/false);
|
|
if (config_.spot_finding.adaptive_threshold)
|
|
w.spot_finder = std::make_unique<AdaptiveSpotFinderGPU>(*prescan_mapping, stream);
|
|
else
|
|
w.spot_finder = std::make_unique<ImageSpotFinderGPU>(prescan_x.GetXPixelsNumConv(),
|
|
prescan_x.GetYPixelsNumConv(), stream);
|
|
return w;
|
|
}
|
|
#endif
|
|
if (want_spots) {
|
|
w.preprocessor = std::make_unique<ImagePreprocessorCPU>(prescan_x, prescan_mask);
|
|
w.preprocessed = std::make_unique<ImagePreprocessorBuffer>(prescan_x.GetPixelsNum());
|
|
if (config_.spot_finding.adaptive_threshold)
|
|
w.spot_finder = std::make_unique<AdaptiveSpotFinderCPU>(*prescan_mapping);
|
|
else
|
|
w.spot_finder = std::make_unique<ImageSpotFinderCPU>(prescan_x.GetXPixelsNumConv(),
|
|
prescan_x.GetYPixelsNumConv());
|
|
}
|
|
return w;
|
|
};
|
|
|
|
// The spots of one image, and what is taken off them.
|
|
const auto find_spots = [&](PreScanWorker &w, const CompressedImage &image, int image_idx,
|
|
bool for_width, std::vector<spot_width::FluxCurve> &curves,
|
|
std::vector<float> &spot_q, std::vector<FootprintSpot> &footprint) {
|
|
try {
|
|
// As in the image loops: a frame the device cannot decode goes to the host decoder.
|
|
ImageStatistics stats;
|
|
bool decoded_on_device = false;
|
|
try {
|
|
decoded_on_device = w.preprocessor->AnalyzeCompressed(*w.preprocessed, image, stats);
|
|
} catch (const JFJochException &e) {
|
|
cuda_throw_if_context_lost();
|
|
logger.Warning("Pre-scan: device decoding of image {} failed ({}), decompressing it on "
|
|
"the host", image_idx, e.what());
|
|
cuda_clear_error();
|
|
}
|
|
if (!decoded_on_device)
|
|
w.preprocessor->Analyze(*w.preprocessed,
|
|
image.GetUncompressedPtr(w.decompression_buffer), image.GetMode());
|
|
} catch (const std::exception &e) {
|
|
if (IsFatalResourceError(e)) throw;
|
|
logger.Warning("Pre-scan: failed to preprocess image {}: {}", image_idx, e.what());
|
|
return;
|
|
}
|
|
auto spots = w.spot_finder->Run(*w.preprocessed, config_.spot_finding);
|
|
// Every frame's WHOLE spot list feeds the powder measurement below - the rings are a property
|
|
// of the sample, and the frames read here are already spread over the sweep, which is what
|
|
// sharpens them against the crystal's own rows.
|
|
{
|
|
const auto &g = prescan_x.GetDiffractionGeometry();
|
|
for (const auto &spot : spots) {
|
|
const Coord c = spot.RawCoord();
|
|
const float d = g.PxlToRes(c.x, c.y);
|
|
if (d > 0.0f && std::isfinite(d))
|
|
spot_q.push_back(2 * PI / d);
|
|
}
|
|
}
|
|
// The width goes first: it wants the whole spot list, to know which spots are isolated.
|
|
// The footprint is measured on every spot, not only the isolated ones the width takes: it is
|
|
// the far, wide spots it is for (SpotFootprint.h).
|
|
if (for_width)
|
|
MeasureSpotsOnTiles(*w.preprocessed, static_cast<int>(prescan_x.GetXPixelsNumConv()),
|
|
static_cast<int>(prescan_x.GetYPixelsNumConv()), prescan_x.GetDiffractionGeometry(),
|
|
spots, curves, footprint);
|
|
};
|
|
|
|
std::vector<float> prescan_spot_q;
|
|
std::vector<spot_width::FluxCurve> width_curves;
|
|
// The spots the bandwidth is measured on: the width's, and more where those are too few
|
|
// (spot_width::BANDWIDTH_POOL_SPOTS).
|
|
std::vector<spot_width::FluxCurve> shape_curves;
|
|
// Every spot's measured widths, for the spot footprint (SpotFootprint.h).
|
|
std::vector<FootprintSpot> footprint_spots;
|
|
// Images the width ended up being measured on, for the log: the tiers below stop as soon as the
|
|
// answer has settled, so this is a property of the crystal and worth reporting.
|
|
size_t width_images = 0;
|
|
// Frames that could actually be read. Counted here rather than asked of the shadow finder, which
|
|
// is only fed when the shadow is wanted and would report zero for a pass that read every frame
|
|
// for the beam centre alone.
|
|
std::atomic<size_t> frames_read{0};
|
|
|
|
// Read the sample on several workers. The reader serialises on the HDF5 lock, but the
|
|
// decompression, the projection and the spot finding - which is all of the cost on a large
|
|
// detector - run in parallel. The projection is integer sums, so the order the workers add their
|
|
// frames in does not reach it, and the per-frame results are stitched together in sample order
|
|
// below so the answer does not depend on how the workers interleaved.
|
|
//
|
|
// Two passes over the sample. The first builds the projection, which is all the shadow, the
|
|
// defective pixels and the beam-centre capture below read; the second finds the spots, for the
|
|
// spot width, the bandwidth and the powder rings, and runs in the
|
|
// background beside them - it depends on none of them and none of them on it. What it measures
|
|
// is applied once both are done (finish_spot_measurement).
|
|
const std::vector<int> ordinals(shadow_set.begin(), shadow_set.end());
|
|
const size_t nworkers = std::min<size_t>(std::max<size_t>(config_.nthreads, 1),
|
|
std::min(PRESCAN_MAX_WORKERS, ordinals.size()));
|
|
// Per frame, in sample order, so the gain does not depend on how the workers interleaved.
|
|
std::vector<background_gain::Levels> gain_levels(measure_gain ? ordinals.size() : 0);
|
|
const size_t gain_frame_stride = (ordinals.size() + background_gain::GAIN_MAX_FRAMES - 1) / background_gain::GAIN_MAX_FRAMES;
|
|
|
|
// The sample as read here, by position, for the spot measurement and the defective pixels below
|
|
// (PRESCAN_KEPT_BYTES).
|
|
const bool keep_frames = want_spots || measure_hot;
|
|
std::vector<JFJochReaderRawImage> kept(keep_frames ? ordinals.size() : 0);
|
|
std::vector<char> is_kept(ordinals.size(), 0);
|
|
std::atomic<size_t> kept_bytes{0};
|
|
// Frame i of the sample: the one kept, or read into `scratch`. Null where it cannot be read.
|
|
const auto frame = [&](size_t i, JFJochReaderRawImage &scratch) -> const CompressedImage * {
|
|
if (is_kept[i])
|
|
return &kept[i].image;
|
|
const int image_idx = start_image + ordinals[i] * config_.stride;
|
|
try {
|
|
if (reader_.ReadRawImage(image_idx, scratch))
|
|
return &scratch.image;
|
|
} catch (const std::exception &e) {
|
|
if (IsFatalResourceError(e)) throw;
|
|
logger.Warning("Pre-scan: failed to load image {}: {}", image_idx, e.what());
|
|
}
|
|
return nullptr;
|
|
};
|
|
|
|
// The spot measurement's engines - the azimuthal mapping its adaptive finder thresholds on, and a
|
|
// preprocessor and a spot finder per worker - read only the experiment and the mask as they are
|
|
// now, so they are built while the frames are read rather than after. The mapping is built on the
|
|
// cores the readers leave free: a frame decoded on the host is on the pre-scan's critical path, and
|
|
// the mapping is the same whatever number of threads builds it.
|
|
std::vector<PreScanWorker> workers;
|
|
std::future<void> engines_built;
|
|
if (want_spots) {
|
|
workers.resize(nworkers);
|
|
engines_built = std::async(std::launch::async, [&] {
|
|
const size_t free_cpus = std::max<size_t>(1, AvailableCpus() - std::min<size_t>(AvailableCpus(), nworkers));
|
|
prescan_mapping = std::make_unique<AzimuthalIntegrationMapping>(prescan_x, prescan_mask,
|
|
*azint_geometry_, free_cpus);
|
|
OnWorkers(nworkers, [&](size_t t) { workers[t] = make_worker(); });
|
|
});
|
|
}
|
|
|
|
if (find_shadow || measure_gain || keep_frames) {
|
|
TimingMark mark(logger, "pre-scan read and projection");
|
|
std::atomic<size_t> next{0};
|
|
OnWorkers(nworkers, [&](size_t) {
|
|
std::vector<uint8_t> shadow_buffer;
|
|
std::vector<int32_t> gain_pixels;
|
|
JFJochReaderRawImage scratch;
|
|
for (size_t i = next.fetch_add(1); i < ordinals.size(); i = next.fetch_add(1)) {
|
|
const int ordinal = ordinals[i];
|
|
const int image_idx = start_image + ordinal * config_.stride;
|
|
// Into a slot of its own while there is room to keep it, otherwise into the worker's
|
|
// buffer, which is reused frame after frame rather than paged in afresh for each. The
|
|
// reader's own scratch (JFJochReaderRawImage::read_buffer) is not part of the image, so
|
|
// a kept frame borrows the worker's for the read and gives it back.
|
|
const bool keep = keep_frames && kept_bytes.load() < PRESCAN_KEPT_BYTES;
|
|
JFJochReaderRawImage &raw_image = keep ? kept[i] : scratch;
|
|
if (keep)
|
|
raw_image.read_buffer = std::move(scratch.read_buffer);
|
|
bool read = false;
|
|
try {
|
|
read = reader_.ReadRawImage(image_idx, raw_image);
|
|
} catch (const std::exception &e) {
|
|
if (IsFatalResourceError(e)) throw;
|
|
logger.Warning("Pre-scan: failed to load image {}: {}", image_idx, e.what());
|
|
}
|
|
if (keep)
|
|
scratch.read_buffer = std::move(raw_image.read_buffer);
|
|
if (!read) continue;
|
|
frames_read.fetch_add(1, std::memory_order_relaxed);
|
|
// Gated on find_shadow, not on the frame set alone: AddImage projects the whole
|
|
// image, and with --detect-beam-stop=off (or a shadow already in the mask) that is
|
|
// the pass's largest per-frame cost paid for a result nobody reads. The set still
|
|
// decides WHICH frames feed it, so the mask is built from the frames it always was.
|
|
if (find_shadow && shadow_set.contains(ordinal)) {
|
|
DataMessage msg{};
|
|
msg.image = raw_image.image;
|
|
msg.number = ordinal;
|
|
msg.original_number = image_idx;
|
|
finder->AddImage(msg, shadow_buffer);
|
|
}
|
|
if (measure_gain && i % gain_frame_stride == 0) {
|
|
const auto &img = raw_image.image;
|
|
const size_t n = img.GetWidth() * img.GetHeight();
|
|
const uint8_t *ptr = img.GetUncompressedPtr(shadow_buffer);
|
|
gain_pixels.resize(n);
|
|
bool measurable = true;
|
|
switch (img.GetMode()) {
|
|
case CompressedImageMode::Int32:
|
|
std::memcpy(gain_pixels.data(), ptr, n * sizeof(int32_t));
|
|
break;
|
|
case CompressedImageMode::Uint16:
|
|
std::copy_n(reinterpret_cast<const uint16_t *>(ptr), n, gain_pixels.begin());
|
|
break;
|
|
case CompressedImageMode::Int16:
|
|
std::copy_n(reinterpret_cast<const int16_t *>(ptr), n, gain_pixels.begin());
|
|
break;
|
|
default:
|
|
measurable = false; // the CCD readers give 16- or 32-bit integers
|
|
}
|
|
if (measurable)
|
|
gain_levels[i] = background_gain::MeasureFrame(gain_pixels, img.GetWidth(),
|
|
img.GetHeight(), pixel_mask_.GetMask(),
|
|
experiment_.GetSaturationLimit());
|
|
}
|
|
if (keep) {
|
|
kept_bytes.fetch_add(raw_image.image_buffer.capacity());
|
|
is_kept[i] = 1;
|
|
}
|
|
}
|
|
});
|
|
logger.Debug("Pre-scan: {} of {} frames kept for the steps after the projection, {:.0f} MB",
|
|
std::count(is_kept.begin(), is_kept.end(), 1), ordinals.size(), kept_bytes.load() / 1e6);
|
|
} else
|
|
frames_read = measured ? measured->frames_read : shadow_memo->frames_read;
|
|
|
|
// Applied before anything below copies the integration settings (bragg_before_adaptive_), so every
|
|
// pass and every fallback integrates with it. At or below 1 - an ADSC or a marCCD, whose pixel
|
|
// scatters less than its mean - the photon clip is already looser than its nominal width and is left
|
|
// exactly as it is.
|
|
if (want_gain) {
|
|
background_gain_measured_ = true;
|
|
if (measure_gain) {
|
|
if (!measured)
|
|
measured.emplace();
|
|
measured->gain = background_gain::Gain(gain_levels);
|
|
}
|
|
const std::optional<float> g = *measured->gain;
|
|
if (g) {
|
|
BraggIntegrationSettings bis = experiment_.GetBraggIntegrationSettings();
|
|
bis.BackgroundPixelGain(std::max(1.0f, *g));
|
|
experiment_.ImportBraggIntegrationSettings(bis);
|
|
}
|
|
if (!g)
|
|
logger.Info("Background gain: not measurable on {} frames; the background clip stays the photon "
|
|
"clip", (ordinals.size() + gain_frame_stride - 1) / gain_frame_stride);
|
|
else
|
|
logger.Info("Background gain: a background pixel scatters {:.2f}x its mean in ADU ({} frames){}",
|
|
*g, (ordinals.size() + gain_frame_stride - 1) / gain_frame_stride,
|
|
*g > 1.0f ? fmt::format(" => background clip at mean + {:.1f}*sqrt({:.2f}*mean)",
|
|
experiment_.GetBraggIntegrationSettings().GetBackgroundClipNSigma(), *g)
|
|
: std::string(" - at or below 1, the photon clip is kept"));
|
|
}
|
|
|
|
const auto spot_measurement = [&] {
|
|
TimingMark mark(logger, "pre-scan spot measurement");
|
|
engines_built.get();
|
|
std::vector<std::vector<spot_width::FluxCurve>> curves_of(ordinals.size());
|
|
std::vector<std::vector<float>> spot_q_of(ordinals.size());
|
|
std::vector<std::vector<FootprintSpot>> footprint_of(ordinals.size());
|
|
// Frames read for the bandwidth alone, past the tier the width settled on. Their spots stay
|
|
// out of the powder measurement, which is taken on the frames it always was.
|
|
std::vector<char> shape_only(ordinals.size(), 0);
|
|
|
|
// The order the frames are visited in, and where each width tier ends in it.
|
|
std::vector<size_t> visit;
|
|
std::vector<size_t> tier_end;
|
|
visit.reserve(ordinals.size());
|
|
{
|
|
std::vector<char> taken(ordinals.size(), 0);
|
|
for (const int stride : WIDTH_TIER_STRIDE) {
|
|
for (size_t j = 0; j < ordinals.size(); j += stride)
|
|
if (!taken[j]) {
|
|
taken[j] = 1;
|
|
visit.push_back(j);
|
|
}
|
|
tier_end.push_back(visit.size());
|
|
}
|
|
}
|
|
// Where each frame is in that order, and the tier it is in.
|
|
std::vector<size_t> position(ordinals.size());
|
|
for (size_t v = 0; v < visit.size(); v++)
|
|
position[visit[v]] = v;
|
|
const auto tier_of = [&](size_t v) {
|
|
return static_cast<size_t>(std::upper_bound(tier_end.begin(), tier_end.end(), v) - tier_end.begin());
|
|
};
|
|
|
|
// The workers go through the frames in that order without waiting at the end of a tier: a frame
|
|
// is measured the same whichever tier it is read for, so the next tier's frames are started
|
|
// while the last of this one's finish and the answer below is drawn. What that answer reads is
|
|
// the frames of the tiers so far and nothing else, and the frames measured past the tier it
|
|
// stops at are thrown away - so how far ahead the workers got changes nothing.
|
|
std::mutex m;
|
|
std::condition_variable cv;
|
|
std::vector<size_t> done_in_tier(tier_end.size(), 0);
|
|
bool failed = false;
|
|
std::atomic<size_t> next{0};
|
|
std::atomic<bool> stop{false};
|
|
auto measuring = std::async(std::launch::async, [&] {
|
|
OnWorkers(nworkers, [&](size_t t) {
|
|
PreScanWorker &w = workers[t];
|
|
try {
|
|
for (size_t v = next.fetch_add(1); v < visit.size() && !stop.load(); v = next.fetch_add(1)) {
|
|
const size_t i = visit[v];
|
|
const int image_idx = start_image + ordinals[i] * config_.stride;
|
|
if (const CompressedImage *image = frame(i, w.raw_image))
|
|
find_spots(w, *image, image_idx, true, curves_of[i], spot_q_of[i], footprint_of[i]);
|
|
{
|
|
std::lock_guard lock(m);
|
|
done_in_tier[tier_of(v)]++;
|
|
}
|
|
cv.notify_all();
|
|
}
|
|
} catch (...) {
|
|
{
|
|
std::lock_guard lock(m);
|
|
failed = true;
|
|
}
|
|
cv.notify_all();
|
|
throw;
|
|
}
|
|
});
|
|
});
|
|
|
|
bool width_settled = false;
|
|
// The bandwidth reads the same spots, but its slope needs more of them than the width does:
|
|
// the tiers carry on past the width's for it alone, until the pool is large enough.
|
|
bool shapes_settled = false;
|
|
std::optional<float> previous_r80;
|
|
size_t measured_end = 0; // the frames, in visit order, the answer is taken on
|
|
for (size_t tier = 0; tier < tier_end.size() && !shapes_settled; tier++) {
|
|
const size_t tier_begin = tier ? tier_end[tier - 1] : 0;
|
|
{
|
|
std::unique_lock lock(m);
|
|
cv.wait(lock, [&] { return failed || done_in_tier[tier] == tier_end[tier] - tier_begin; });
|
|
if (failed) break;
|
|
}
|
|
measured_end = tier_end[tier];
|
|
// A tier begun once the width had settled is read for the bandwidth alone.
|
|
if (width_settled)
|
|
for (size_t v = tier_begin; v < measured_end; v++)
|
|
shape_only[visit[v]] = 1;
|
|
|
|
// What this tier makes of the width, pooled in sample order so the answer does not
|
|
// depend on how the workers interleaved. A tier too sparse to measure settles nothing:
|
|
// the test is against what the smaller sample actually said, not against the radius it
|
|
// fell back to.
|
|
std::vector<spot_width::FluxCurve> pooled;
|
|
for (size_t i = 0; i < ordinals.size(); i++)
|
|
if (position[i] < measured_end)
|
|
pooled.insert(pooled.end(), curves_of[i].begin(), curves_of[i].end());
|
|
if (!width_settled) {
|
|
width_images = measured_end;
|
|
width_curves = pooled;
|
|
const auto r80 = spot_width::R80AtReference(width_curves);
|
|
width_settled = r80 && previous_r80 && spot_width::WidthSettled(*r80, *previous_r80);
|
|
previous_r80 = r80;
|
|
}
|
|
shape_curves = std::move(pooled);
|
|
shapes_settled = width_settled && shape_curves.size() >= spot_width::BANDWIDTH_POOL_SPOTS;
|
|
}
|
|
stop = true;
|
|
measuring.get();
|
|
OnWorkers(nworkers, [&](size_t t) { PreScanWorker released = std::move(workers[t]); });
|
|
|
|
for (size_t i = 0; i < ordinals.size(); i++)
|
|
if (position[i] >= measured_end) {
|
|
spot_q_of[i].clear();
|
|
footprint_of[i].clear();
|
|
}
|
|
for (size_t i = 0; i < ordinals.size(); i++)
|
|
if (!shape_only[i])
|
|
prescan_spot_q.insert(prescan_spot_q.end(), spot_q_of[i].begin(), spot_q_of[i].end());
|
|
// In sample order, so the table does not depend on how the workers interleaved.
|
|
for (const auto &f : footprint_of)
|
|
footprint_spots.insert(footprint_spots.end(), f.begin(), f.end());
|
|
};
|
|
std::future<void> spots_measured;
|
|
if (want_spots)
|
|
spots_measured = std::async(std::launch::async, spot_measurement);
|
|
|
|
// What the spot measurement found, applied to the run once it has finished: called on every way
|
|
// out of the pre-scan below.
|
|
const auto finish_spot_measurement = [&] {
|
|
if (spots_measured.valid())
|
|
spots_measured.get();
|
|
if (want_width && measure_width) {
|
|
if (!measured)
|
|
measured.emplace();
|
|
measured->spots = PreScanMemo::Spots{width_curves, shape_curves, footprint_spots, prescan_spot_q,
|
|
width_images};
|
|
} else if (want_width) {
|
|
logger.Info("Spot width: the spots kept from an earlier pre-scan of this run on the same frames, "
|
|
"geometry, mask and spot finding");
|
|
const auto &k = *measured->spots;
|
|
width_curves = k.width_curves;
|
|
shape_curves = k.shape_curves;
|
|
footprint_spots = k.footprint;
|
|
prescan_spot_q = k.spot_q;
|
|
width_images = k.width_images;
|
|
}
|
|
|
|
// The integration radius from the recorded spot width. r1 is the aperture the integrator learns
|
|
// the profile width over and a second moment over a disk of radius a saturates at a^2/4, so at the
|
|
// shipped r1 = 4 the learned sigma cannot exceed 2 px and a wider spot is fitted with a profile
|
|
// the model cannot represent. r2 follows r1, and r3 is taken so that the r2..r3 background ring
|
|
// keeps the area it has at the shipped 4/6/13 - a ring that shrinks with the disk is what makes a
|
|
// bare --integration-radius worse than the default it replaces.
|
|
if (want_width) {
|
|
spot_width_measured_ = true;
|
|
const auto r80 = spot_width::R80AtReference(width_curves);
|
|
if (!r80) {
|
|
logger.Info("Spot width: not measurable on {} spots from {} images, keeping the "
|
|
"integration radius at r1={:.1f}", width_curves.size(), width_images,
|
|
experiment_.GetBraggIntegrationSettings().GetR1());
|
|
} else {
|
|
const float r1 = spot_width::R1ForWidth(*r80);
|
|
const float r2 = r1 + 2.0f;
|
|
const float r3 = std::sqrt(r2 * r2 + 133.0f);
|
|
BraggIntegrationSettings bis = experiment_.GetBraggIntegrationSettings();
|
|
// Remembered where the APERTURE moved, so the starvation guard in RunAllPasses has
|
|
// something to give back and cannot announce a fallback that changes nothing.
|
|
//
|
|
// The test is on all three radii and not on r1 alone. It is the same test today - r2 and
|
|
// r3 are taken from r1 just above, and the adaptive rule does not run at all when the
|
|
// radii were given explicitly - so this changes no run. It is written this way because
|
|
// what the guard is about is the aperture, and r1 is only a proxy for it. A prototype
|
|
// that widened r2 and r3 at high angle while leaving r1 at the shipped 4 took the
|
|
// neighbour-starved fraction of a dense low-symmetry crystal from 2.4% to 9.0%, eight
|
|
// times the bound this guard exists to enforce, and the guard could not fire: nothing
|
|
// had set bragg_before_adaptive_, because r1 had not moved.
|
|
if (r1 != bis.GetR1() || r2 != bis.GetR2() || r3 != bis.GetR3())
|
|
bragg_before_adaptive_ = bis;
|
|
bis.R1(r1).R2(r2).R3(r3);
|
|
experiment_.ImportBraggIntegrationSettings(bis);
|
|
logger.Info("Spot width: r80 = {:.2f} px at {:.0f} A ({} spots from {} images) => "
|
|
"integration radii r1={:.1f} r2={:.1f} r3={:.2f}", *r80, spot_width::D_REF_A,
|
|
width_curves.size(), width_images, r1, r2, r3);
|
|
Note(fmt::format("spot width {:.1f} px: integration radii {:.1f} / {:.1f} / {:.1f} px",
|
|
*r80, r1, r2, r3));
|
|
}
|
|
|
|
// How far the spots reach away from the beam, where the radius chosen above at 5 A no
|
|
// longer holds them (SpotFootprint.h). It changes the integration only where a spot
|
|
// outgrows the r1 disk.
|
|
const auto &g = experiment_.GetDiffractionGeometry();
|
|
const float bx = g.GetBeamX_pxl(), by = g.GetBeamY_pxl();
|
|
const float W = static_cast<float>(experiment_.GetXPixelsNum());
|
|
const float H = static_cast<float>(experiment_.GetYPixelsNum());
|
|
const SpotFootprint fp = FootprintFromSpots(footprint_spots,
|
|
std::hypot(std::max(bx, W - bx), std::max(by, H - by)));
|
|
if (!fp.empty())
|
|
prescan_footprint_ = fp;
|
|
// It acts only where a spot outgrows the r1 disk (BraggStencil.h), so a pattern of compact
|
|
// spots is left on exactly the settings - and the passes - it had without it.
|
|
const float r1_now = experiment_.GetBraggIntegrationSettings().GetR1();
|
|
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]) > r1_now;
|
|
if (fp.empty()) {
|
|
logger.Info("Spot footprint: too few spots to tabulate ({})", footprint_spots.size());
|
|
} else if (!outgrows) {
|
|
logger.Info("Spot footprint from {} spots: every spot fits the r1={:.1f} disk", footprint_spots.size(),
|
|
r1_now);
|
|
} else {
|
|
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("Spot footprint from {} spots, sigma radial/tangential [px] by distance from "
|
|
"the beam [px]:{}", footprint_spots.size(), table);
|
|
// On the adaptive side, like the radius: the geometry pre-pass integrates without it,
|
|
// and a canonical pass whose wider rings the neighbours starve falls back to the
|
|
// settings before it (bragg_before_adaptive_, see RunAllPasses).
|
|
BraggIntegrationSettings bis = experiment_.GetBraggIntegrationSettings();
|
|
if (!bragg_before_adaptive_)
|
|
bragg_before_adaptive_ = bis;
|
|
bis.Footprint(fp);
|
|
experiment_.ImportBraggIntegrationSettings(bis);
|
|
}
|
|
}
|
|
|
|
// The X-ray bandwidth, read off the shapes of the same spots (spot_width::EstimateBandwidth). A
|
|
// value the file states, or --bandwidth, wins; otherwise a significant estimate is the bandwidth
|
|
// of the run - prediction, partiality, the profile's radial width and the background stencil all
|
|
// take it - and anything short of significant leaves the beam monochromatic.
|
|
if (want_width) {
|
|
const auto &det = experiment_.GetDetectorSetup();
|
|
const double lambda = experiment_.GetWavelength_A();
|
|
const auto estimate = spot_width::EstimateBandwidth(
|
|
shape_curves, sensor_absorption::AttenuationLength_um(det.GetSensorMaterial(), lambda),
|
|
det.GetSensorThickness_um(), experiment_.GetPixelSize_mm() * 1000.0);
|
|
const bool stated = experiment_.GetBandwidthFWHM().has_value();
|
|
const bool significant = estimate && estimate->z > spot_width::BANDWIDTH_Z;
|
|
if (!stated && significant)
|
|
experiment_.BandwidthFWHM(static_cast<float>(estimate->fwhm));
|
|
const float used = experiment_.GetBandwidthFWHM().value_or(0.0f);
|
|
const char *why = stated ? "as stated" : significant ? "as measured" : "monochromatic";
|
|
if (!estimate)
|
|
logger.Info("Bandwidth: not measurable on {} spots; the run uses {:.4f} ({})", shape_curves.size(),
|
|
used, why);
|
|
else
|
|
logger.Info("Bandwidth: spot shapes give FWHM {:.4f} from {} spots (one standard error is "
|
|
"FWHM {:.4f}; z = {:.1f}, chi2 {:.1f}) - {}; the run uses {:.4f} ({})", estimate->fwhm,
|
|
estimate->spots, estimate->fwhm_floor, estimate->z, estimate->chi2,
|
|
significant ? "significant" : "not significant", used, why);
|
|
}
|
|
|
|
// Powder contamination, measured on every run that finds spots here. A crystalline phase other
|
|
// than the crystal - hexagonal ice, a shower of microcrystals, salt out of the cryoprotectant -
|
|
// diffracts as rings, and its spots are handed to the indexer and to scaling as if they were this
|
|
// crystal's. Only ice can be named in advance; the rest has to be measured, and measuring it costs
|
|
// nothing here, the spots having been found for the width and the centre already.
|
|
//
|
|
// Reported on every run whatever it says. What is DONE about it is decided later and on evidence:
|
|
// a first pass that fails indexes again with the rings set aside and keeps whichever answer indexes
|
|
// more frames (see RunPipeline). This line is the diagnostic, and it is the one that tells a user
|
|
// their crystal was in a powder.
|
|
{
|
|
const float w = config_.spot_finding.ice_ring_width_Q_recipA;
|
|
powder_ = MeasurePowderRings(prescan_spot_q, w);
|
|
prescan_spot_quantile_d_A_.clear();
|
|
for (const float frac : {0.25f, 0.5f})
|
|
if (const auto d = SpotResolutionQuantile(prescan_spot_q, frac))
|
|
prescan_spot_quantile_d_A_.push_back(*d);
|
|
if (!powder_.rings_q_recipA.empty()) {
|
|
logger.Info("Powder rings: {} rings over {} pooled spots hold {:.0f}% of them above the "
|
|
"smooth spot density - the coarsest at {:.2f} A, the finest at {:.2f} A{}",
|
|
powder_.rings_q_recipA.size(), prescan_spot_q.size(),
|
|
100.0f * powder_.spot_fraction,
|
|
2 * PI / powder_.rings_q_recipA.front(),
|
|
2 * PI / powder_.rings_q_recipA.back(),
|
|
powder_.resolved_to_d_A
|
|
? fmt::format("; past {:.2f} A the rings have merged and cannot be told "
|
|
"from the crystal", *powder_.resolved_to_d_A)
|
|
: std::string());
|
|
Note(fmt::format("{} powder ring{} ({} ice), {:.0f}% of the spots",
|
|
powder_.rings_q_recipA.size(), powder_.rings_q_recipA.size() == 1 ? "" : "s",
|
|
powder_.ice_ring_count, 100.0f * powder_.spot_fraction));
|
|
}
|
|
}
|
|
};
|
|
|
|
if (frames_read.load() == 0) {
|
|
logger.Warning("Pre-scan: no image could be read. Skipping.");
|
|
finish_spot_measurement();
|
|
return;
|
|
}
|
|
|
|
// Where the scattered background puts the beam, when the projection exists to say so.
|
|
std::optional<std::pair<float, float>> ring_centre;
|
|
std::vector<uint32_t> opaque_shadow;
|
|
// The mean projection, taken once and shared read-only by everything below that reads it - the
|
|
// ring-centre fit here, and the diagnostic picture and the beam-centre capture on threads of
|
|
// their own. It is 4 bytes a pixel, and a copy for each of them was a copy on this thread.
|
|
std::shared_ptr<const std::vector<float>> projection;
|
|
if (want_shadow) {
|
|
// Every pixel is compared against the ring it sits on, so the rings have to be drawn about
|
|
// the beam. A centre far enough out draws them across the background's own radial fall-off
|
|
// instead of along it, and what comes back describes the fall-off rather than the hardware:
|
|
// displacing a centre by 50 px on a clean sweep leaves the mask where it was, by 150 px
|
|
// masks a tenth of the detector. A centre in a file can be that wrong - one sweep here
|
|
// states one 351 px from where the data put it - so the rings are drawn about the centre
|
|
// this same projection measures rather than the one the file claims. That fit costs no
|
|
// frame of its own, and an unmasked shadow biases it by a few pixels at most, which is
|
|
// three orders below what this comparison would notice. The centre the run reports and
|
|
// consumes is the one fitted below instead, with the shadow out of the way.
|
|
std::vector<uint32_t> shadow;
|
|
size_t shadow_frames = 0;
|
|
if (shadow_memo) {
|
|
projection = shadow_memo->projection;
|
|
ring_centre = shadow_memo->ring_centre;
|
|
shadow.assign(experiment_.GetXPixelsNumConv() * experiment_.GetYPixelsNumConv(), 0);
|
|
for (const auto &[i, v] : shadow_memo->shadow)
|
|
shadow[i] = v;
|
|
shadow_frames = shadow_memo->shadow_frames;
|
|
logger.Info("{} - kept from an earlier pre-scan of this run on the same frames, geometry and mask",
|
|
shadow_memo->centre_line);
|
|
} else {
|
|
TimingMark mark(logger, "pre-scan beam stop");
|
|
std::string centre_line;
|
|
projection = std::make_shared<const std::vector<float>>(finder->GetMeanProjection());
|
|
if (const auto ring_center = FindBeamCenterFromBackground(experiment_, pixel_mask_,
|
|
*projection)) {
|
|
finder->BeamCenter(ring_center->beam_x_pxl, ring_center->beam_y_pxl);
|
|
ring_centre = std::make_pair(ring_center->beam_x_pxl, ring_center->beam_y_pxl);
|
|
centre_line = fmt::format("Beam stop: comparing each pixel against its ring about the measured "
|
|
"centre ({:.2f},{:.2f}) +- {:.2f} px, {:.2f} px from the file's",
|
|
ring_center->beam_x_pxl, ring_center->beam_y_pxl, ring_center->sigma_pxl,
|
|
std::hypot(ring_center->beam_x_pxl - experiment_.GetBeamX_pxl(),
|
|
ring_center->beam_y_pxl - experiment_.GetBeamY_pxl()));
|
|
} else {
|
|
centre_line = fmt::format("Beam stop: the background does not place a centre on this run, so the "
|
|
"rings are drawn about the file's ({:.2f},{:.2f})",
|
|
experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl());
|
|
}
|
|
logger.Info("{}", centre_line);
|
|
shadow = finder->GetMask(config_.nthreads);
|
|
shadow_frames = finder->GetFrameCount();
|
|
finder.reset();
|
|
|
|
BeamStopMemo::Entry entry{shadow_key, {}, projection, ring_centre, centre_line, frames_read.load(),
|
|
shadow_frames};
|
|
for (size_t i = 0; i < shadow.size(); i++)
|
|
if (shadow[i] != 0)
|
|
entry.shadow.emplace_back(static_cast<uint32_t>(i), shadow[i]);
|
|
std::lock_guard lock(beam_stop_memo_->m);
|
|
auto &entries = beam_stop_memo_->entries;
|
|
// A handful covers the geometries a run's probes pre-scan at; the oldest go first.
|
|
if (entries.size() == BeamStopMemo::MAX_ENTRIES)
|
|
entries.erase(entries.begin());
|
|
entries.push_back(std::move(entry));
|
|
}
|
|
const auto shadowed = std::count_if(shadow.begin(), shadow.end(), [](uint32_t v) { return v != 0; });
|
|
pixel_mask_.LoadBeamStopMask(experiment_, shadow);
|
|
opaque_shadow = shadow;
|
|
for (auto &v : opaque_shadow)
|
|
if (v == ShadowFinder::TRANSMITTING)
|
|
v = 0;
|
|
// Not "this beamline has no beam stop": it is also what a comparison drawn about the wrong
|
|
// point returns, so an unexpected answer here is worth reading against the centre named
|
|
// above first.
|
|
if (shadowed == 0)
|
|
logger.Warning("Beam stop: no shadow found in {} images. The comparison is made about "
|
|
"the centre named above, so this may point at that centre rather than "
|
|
"at the experiment",
|
|
shadow_frames);
|
|
else {
|
|
logger.Info("Beam stop shadow: {} pixels ({:.2f}% of the detector) found in {} images",
|
|
shadowed, 100.0 * static_cast<double>(shadowed) / static_cast<double>(shadow.size()),
|
|
shadow_frames);
|
|
Note(fmt::format("beam stop masked: {:.2f}% of the detector",
|
|
100.0 * static_cast<double>(shadowed) / static_cast<double>(shadow.size())));
|
|
}
|
|
|
|
// The diagnostic picture, rendered here because this is where the projection and the mask
|
|
// the detection has just written into it are both in hand; the file is written later, under
|
|
// the canonical prefix, since on a two-pass run this is the pre-pass.
|
|
if (config_.mode == ProcessMode::FullAnalysis && experiment_.IsRotationIndexing())
|
|
detector_diagnostic_jpeg_ = std::async(
|
|
std::launch::async, [x = experiment_, mask = pixel_mask_, projection] {
|
|
return RenderDetectorDiagnostic(x, mask, *projection);
|
|
}).share();
|
|
}
|
|
|
|
// The beam centre is fitted again here, so that the shadow is in the mask by the time it is
|
|
// placed and the holder arm - the largest azimuthal asymmetry on the detector - is out of the
|
|
// way of the background the centre is read from. This is the answer the run reports and
|
|
// consumes; the fit above only chose where to draw the rings, and is thrown away.
|
|
|
|
// What the scattered background makes the centre, measured on EVERY run and committed on none.
|
|
// The fit reads the projection the beam-stop pre-pass has just built and no frame of its own, so
|
|
// it costs a fraction of a second on top of a pass that has already run - which is what makes it
|
|
// affordable to do always, and doing it always is the point. A header beam centre is the metadata
|
|
// field the field jokes about, and nothing in any package says how wrong a given one is: on a run
|
|
// that indexes and merges perfectly well this line still reports that the file claims a centre the
|
|
// data put 73 px away. The run is not touched by it. What consumes the number is the second first
|
|
// pass in RunPipeline.
|
|
//
|
|
// Only where the projection exists: with --detect-beam-stop=off nothing accumulates it, and
|
|
// projecting every frame for this alone is the largest per-frame cost of the pre-scan. On the
|
|
// second pass of a two-pass run the shadow is already in the mask, so this does not run again -
|
|
// by then the centre has been post-refined and the file's value is no longer what is in it.
|
|
if (want_hot) {
|
|
if (!measured)
|
|
measured.emplace();
|
|
MaskDefectivePixels(start_image, shadow_sample, ring_centre, frame, measured->defects);
|
|
}
|
|
|
|
// The centre is fitted with the shadow masked but NOT the transmitting pieces the shadow finder
|
|
// adds after it. Their dimming is multiplicative in azimuth, which the fit's per-sector amplitude
|
|
// absorbs; masked, a thin arm running into the beam takes cells out of the rings it crosses, and
|
|
// which rings stay covered then changes as the trial centre moves - a walk seeded beside such an
|
|
// arm stalled on that edge 44 px from the centre it reaches with the arm left in.
|
|
PixelMask centre_mask = pixel_mask_;
|
|
if (!opaque_shadow.empty())
|
|
centre_mask.LoadBeamStopMask(experiment_, opaque_shadow);
|
|
|
|
// It reads nothing the steps after it change, and nothing reads its answer before the first pass's
|
|
// beam-centre check (RunPipeline), so it runs beside them - on copies of what it measures - and is
|
|
// taken there (JoinBeamCenterCapture). It logs from its own thread.
|
|
if (want_shadow && config_.beam_center_check)
|
|
background_center_capture_ = std::async(std::launch::async, [x = experiment_, mask = centre_mask,
|
|
projection,
|
|
nthreads = config_.nthreads] {
|
|
Logger logger("Rugnux");
|
|
const auto t0 = std::chrono::steady_clock::now();
|
|
BeamCenterFFTResult capture;
|
|
const std::optional<BeamCenterEstimate> background_center =
|
|
FindBeamCenter(x, mask, *projection, nthreads, &capture);
|
|
const auto capture_ms = std::chrono::duration_cast<std::chrono::milliseconds>(
|
|
std::chrono::steady_clock::now() - t0).count();
|
|
if (!capture.point.empty())
|
|
logger.Info("Beam centre capture: the strongest point-symmetry centre on the WHOLE "
|
|
"detector is ({:.1f},{:.1f}), {:.1f} px from the file's; line mirrors put x "
|
|
"at {:.1f} and y at {:.1f}; margins {:.1f}%/{:.1f}%/{:.1f}% (point/x/y, "
|
|
"under ~3% the surface is flat). One transform set scores every candidate "
|
|
"centre, so this costs the same however wrong the file is - {} ms including "
|
|
"the walk that refines it",
|
|
capture.point[0].beam_x_pxl, capture.point[0].beam_y_pxl,
|
|
std::hypot(capture.point[0].beam_x_pxl - x.GetBeamX_pxl(),
|
|
capture.point[0].beam_y_pxl - x.GetBeamY_pxl()),
|
|
capture.line_x.empty() ? std::numeric_limits<float>::quiet_NaN()
|
|
: capture.line_x[0].beam_x_pxl,
|
|
capture.line_y.empty() ? std::numeric_limits<float>::quiet_NaN()
|
|
: capture.line_y[0].beam_y_pxl,
|
|
100.0f * capture.margin_point, 100.0f * capture.margin_line_x,
|
|
100.0f * capture.margin_line_y, capture_ms);
|
|
const float need = BeamCenterNeed_pxl(x);
|
|
if (!background_center) {
|
|
logger.Info("Beam centre check: neither the capture nor the walk places the centre on "
|
|
"this run - the file's ({:.2f},{:.2f}) is neither confirmed nor contradicted",
|
|
x.GetBeamX_pxl(), x.GetBeamY_pxl());
|
|
} else {
|
|
// Against the DIRECT BEAM, not the PONI. Both the FFT capture and the walk read the
|
|
// centrosymmetry of the scattered background about the point the beam lands on; the
|
|
// geometry's beam_x_pxl is the PONI, which parts from it by distance*tan(rot)/pixel as
|
|
// soon as the detector is tilted. Measured on a calibrant at 110/150/200/300 mm the
|
|
// capture sits 0.5-1.5 px from the fitted direct beam and 5.0/6.8/8.4/17.9 px from the
|
|
// fitted PONI, so comparing it with the PONI reports a disagreement the data do not
|
|
// have - on the in-house files, whose header carries a 0.234 deg tilt, 3.3 px of it at
|
|
// 120 mm against a stated need of 0.7 px.
|
|
const auto [direct_x, direct_y] = x.GetDiffractionGeometry().GetDirectBeam_pxl();
|
|
const float dx = background_center->beam_x_pxl - direct_x;
|
|
const float dy = background_center->beam_y_pxl - direct_y;
|
|
const float moved = std::hypot(dx, dy);
|
|
const auto across = AcrossSpindle_pxl(x, dx, dy);
|
|
// Everything a reader needs to judge the header, on one line: what the file says, what the
|
|
// background says, how far apart they are, how well the background knows its own answer,
|
|
// and what this particular geometry asks of the centre. The last is what turns a bare pixel
|
|
// count into a verdict - 3 px is nothing at 100 mm and a lost lattice at 500 mm.
|
|
logger.Info("Beam centre check: the file puts the direct beam at ({:.2f},{:.2f}), the "
|
|
"scattered background "
|
|
"puts it at ({:.2f},{:.2f}) +- {:.2f} px - a difference of {:.2f} px, against the "
|
|
"{:.2f} px this geometry asks the centre to be right to. Measured only; the "
|
|
"run keeps the centre in the file",
|
|
direct_x, direct_y,
|
|
background_center->beam_x_pxl, background_center->beam_y_pxl,
|
|
background_center->sigma_pxl, moved, need);
|
|
// Three sigma, because the fit is a PRECISION and not an accuracy: on a background with no
|
|
// curvature there is nothing to separate shift from amplitude and it follows the noise in
|
|
// g', confidently. A move it cannot resolve is not evidence about the header either way.
|
|
if (moved <= 3.0f * background_center->sigma_pxl)
|
|
logger.Info("Beam centre check: that is under three times the fit's own sigma, so the "
|
|
"file's centre is as good as this measurement can tell - the second first "
|
|
"pass below indexes it anyway, because a fraction of a pixel is enough to "
|
|
"decide a cell and no sigma small enough to gate on would notice");
|
|
else if (across && *across > need)
|
|
// Info and not a warning: `need` is where the FFT amplitude at a 200 A axis has fallen
|
|
// to 0.70, which is a floor under what matters and not a threshold for failure - it is
|
|
// exceeded on two of three ordinary in-house sweeps whose lattice comes out right. What
|
|
// is worth a warning is what the second pass below then finds, not this.
|
|
logger.Info("Beam centre check: {:.2f} px of the difference is ACROSS the spindle, more "
|
|
"than the {:.2f} px this geometry absorbs, and a centre error is fixed in "
|
|
"the lab frame - it smears the accumulated reciprocal-space cloud, so the "
|
|
"first pass can lose peaks or take an axis harmonic. The second first pass "
|
|
"below indexes both centres and reports whether it did", *across, need);
|
|
else
|
|
logger.Info("Beam centre check: that is a real difference, and {} - the second first "
|
|
"pass below indexes both centres and compares what they give",
|
|
across ? "it is mostly along the spindle, where nothing announces itself"
|
|
: "this is not a rotation sweep, so there is no spindle to resolve it on");
|
|
}
|
|
return background_center;
|
|
}).share();
|
|
|
|
finish_spot_measurement();
|
|
|
|
if (measured) {
|
|
measured->key = measure_key;
|
|
measured->frames_read = frames_read.load();
|
|
std::lock_guard lock(prescan_memo_->m);
|
|
auto &entries = prescan_memo_->entries;
|
|
std::erase_if(entries, [&](const auto &e) { return e.key == measure_key; });
|
|
if (entries.size() == PreScanMemo::MAX_ENTRIES)
|
|
entries.erase(entries.begin());
|
|
entries.push_back(std::move(*measured));
|
|
}
|
|
}
|
|
|
|
void Rugnux::MaskDefectivePixels(int start_image, const std::vector<int> &sample,
|
|
const std::optional<std::pair<float, float>> &ring_centre,
|
|
const std::function<const CompressedImage *(size_t, JFJochReaderRawImage &)> &frame,
|
|
std::optional<PreScanMemo::Defects> &defects) {
|
|
Logger logger("Rugnux");
|
|
TimingMark mark(logger, "pre-scan defective pixels");
|
|
hot_pixels_measured_ = true;
|
|
|
|
const auto apply = [&] {
|
|
if (defects->frames == 0)
|
|
return;
|
|
std::vector<uint32_t> mask(defects->pixels, 0);
|
|
for (const auto &[i, v] : defects->mask)
|
|
mask[i] = v;
|
|
pixel_mask_.LoadHotPixelMask(experiment_, mask);
|
|
const size_t masked = defects->hot + defects->error + defects->dead;
|
|
logger.Info("Defective pixels: {} masked from {} pre-scan frames - {} lit above their ring on more "
|
|
"frames than one reflection or chance explains, {} holding the detector's error value, {} dead",
|
|
masked, defects->frames, defects->hot, defects->error, defects->dead);
|
|
Note(fmt::format("{} defective pixel{} masked", masked, masked == 1 ? "" : "s"));
|
|
};
|
|
if (defects) {
|
|
apply();
|
|
return;
|
|
}
|
|
|
|
// The rings each pixel is compared against have to be drawn about the beam, and the centre the
|
|
// file states can be far enough out that a ring crosses the background's own radial fall-off - the
|
|
// pixels on its bright side then read high on every frame. The centre the beam-stop projection
|
|
// measures does not exist until every frame of the first read is in, so the frames are gone
|
|
// through a second time - from memory, as the pre-scan kept them, where it could.
|
|
DiffractionExperiment about = experiment_;
|
|
if (ring_centre)
|
|
about.BeamX_pxl(ring_centre->first).BeamY_pxl(ring_centre->second);
|
|
const size_t nthreads = static_cast<size_t>(std::max(config_.nthreads, 1));
|
|
HotPixelFinder finder(about, pixel_mask_, nthreads);
|
|
|
|
#ifdef JFJOCH_USE_CUDA
|
|
// Built here, once, so the workers' first frames do not queue behind it.
|
|
if (get_gpu_count() > 0)
|
|
finder.PrepareDevice();
|
|
#endif
|
|
std::atomic<size_t> next{0};
|
|
const size_t nworkers = std::min<size_t>(nthreads, std::min(PRESCAN_MAX_WORKERS, sample.size()));
|
|
OnWorkers(nworkers, [&](size_t) {
|
|
// With a GPU the frame is decoded, preprocessed and analysed on the device - only the
|
|
// compressed chunk crosses PCIe - and the result is the same to the bit (HotPixelsGPU.h).
|
|
std::unique_ptr<ImagePreprocessor> preprocessor;
|
|
std::unique_ptr<ImagePreprocessorBuffer> image;
|
|
#ifdef JFJOCH_USE_CUDA
|
|
std::optional<HotPixelFinderGPU::Frame> device_frame;
|
|
if (get_gpu_count() > 0) {
|
|
auto stream = std::make_shared<CudaStream>();
|
|
preprocessor = std::make_unique<ImagePreprocessorGPU>(experiment_, pixel_mask_, stream,
|
|
/*copy_image_to_host=*/false);
|
|
image = std::make_unique<ImagePreprocessorBufferGPU>(experiment_.GetPixelsNum(),
|
|
/*host_mirror=*/false);
|
|
device_frame.emplace(stream);
|
|
}
|
|
#endif
|
|
if (!preprocessor) {
|
|
preprocessor = std::make_unique<ImagePreprocessorCPU>(experiment_, pixel_mask_);
|
|
image = std::make_unique<ImagePreprocessorBuffer>(experiment_.GetPixelsNum());
|
|
}
|
|
JFJochReaderRawImage raw;
|
|
std::vector<uint8_t> buffer;
|
|
std::vector<int32_t> scratch;
|
|
for (size_t i = next.fetch_add(1); i < sample.size(); i = next.fetch_add(1)) {
|
|
const int image_idx = start_image + sample[i] * config_.stride;
|
|
const CompressedImage *compressed = frame(i, raw);
|
|
if (!compressed) continue;
|
|
try {
|
|
// As in the image loops (MXAnalysisWithoutFPGA::Analyze): a frame the device cannot
|
|
// decode goes to the host decoder, so a GPU-side problem costs speed, not the frame.
|
|
ImageStatistics stats;
|
|
bool decoded_on_device = false;
|
|
try {
|
|
decoded_on_device = preprocessor->AnalyzeCompressed(*image, *compressed, stats);
|
|
} catch (const JFJochException &e) {
|
|
cuda_throw_if_context_lost();
|
|
logger.Warning("Defective pixels: device decoding of image {} failed ({}), "
|
|
"decompressing it on the host", image_idx, e.what());
|
|
cuda_clear_error();
|
|
}
|
|
if (!decoded_on_device)
|
|
preprocessor->Analyze(*image, compressed->GetUncompressedPtr(buffer), compressed->GetMode());
|
|
} catch (const std::exception &e) {
|
|
if (IsFatalResourceError(e)) throw;
|
|
logger.Warning("Defective pixels: failed to read image {}: {}", image_idx, e.what());
|
|
continue;
|
|
}
|
|
#ifdef JFJOCH_USE_CUDA
|
|
if (device_frame) {
|
|
finder.AddDeviceImage(image->getGPUBuffer(), *device_frame);
|
|
continue;
|
|
}
|
|
#endif
|
|
finder.AddImage(image->data(), scratch);
|
|
}
|
|
});
|
|
|
|
// On stills there is no rotation to count a reflection's stay on a pixel in, so spacing stays 0 and
|
|
// the hot-pixel rule is not applied (HotPixelFinder::GetMask); the error-value and dead pixels are.
|
|
const auto goniometer = experiment_.GetGoniometer();
|
|
const double oscillation = goniometer ? std::abs(goniometer->GetIncrement_deg()) : 0.0;
|
|
const double spacing = goniometer && sample.size() > 1
|
|
? oscillation * (sample.back() - sample.front()) / static_cast<double>(sample.size() - 1)
|
|
: 0.0;
|
|
const auto hot = finder.GetMask(oscillation, spacing, nthreads);
|
|
defects = PreScanMemo::Defects{hot.mask.size(), {}, hot.hot, hot.error, hot.dead, hot.frames};
|
|
for (size_t i = 0; i < hot.mask.size(); i++)
|
|
if (hot.mask[i] != 0)
|
|
defects->mask.emplace_back(static_cast<uint32_t>(i), hot.mask[i]);
|
|
apply();
|
|
}
|