rugnux: mask a shadow that turns with the sample, found in blocks of the pre-scan
A pin, loop or capillary on the goniometer shadows one place for part of the rotation only. The beam-stop mask is one for every frame and is read off the mean over the pre-scan sample, so a pixel such hardware covers on fewer than about half the sampled frames is kept and integrated attenuated on the frames it is covered on. Since the pre-scan samples only the first part of the sweep (part B,5b9144410), a sweep whose moving shadow is short in the first half and long in the second loses it: on the in-house half-image lysozyme sweep the mask fell from 32.8% to 21.6% of the detector, the second half of the sweep was integrated through the shadow, ISa dropped 6.9 -> 4.1, R_meas 116 -> 240%, the cell refined a != b, and the attenuated intensities erased the a-axis screw evidence (merged rc175: P4_1 2 2 instead of P4_1 2_1 2). The sample is now also read in blocks of 10 consecutive frames (rotation only); each block's shadow is found on its own with the rings about the whole sample's centre, and the blocks' pieces that reach beyond the 30 px band the whole sample's edge can grow over (PENUMBRA_MAX_PX) are added whole. A pixel shadowed on most frames of any longer stretch is shadowed on most frames of one block, so nothing finer is needed; pieces within the band are a block redrawing a static stop's edge and are left out, so a sweep without moving hardware keeps its mask bit for bit. Part B's principle is kept: everything is still decided on the first part of the sweep. Measured (GPU, against merged rc175eac0833c8, 19 sets): half-image lysozyme P4_1 2_1 2, ISa 4.1 -> 9.7, R_meas 240 -> 113%, CC1/2 0.955 -> 0.979, R-free(model) 0.331 -> 0.282 (last night 6.9 / 116%). p.mtz md5-identical on 14 of 19 (8v4o, 9rci, 5ky6, 6yqf, 6ze4, 5reo, 9qw8, insu_I low-ISa, insu_H notwin, myob powder 1/2, myob split, nothing_1/2); changed: 5nw5 (a dim crescent at the stop edge now masked: ISa 8.4 -> 8.9, CC1/2 0.937 -> 0.980, R-free 0.462 -> 0.429), 7atg (0.85% added at the detector edge, statistics unchanged), 9min (0.03%, fail -> fail unchanged), 5jk4 (<0.01%, unchanged). Cost: pre-scan +0.1 s on HDF5 sweeps, +1-5 s on large file-per-image (CBF, marCCD) sweeps, where the blocks are re-read and their masks computed on the host. 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:
@@ -370,6 +370,23 @@ reflection — a maximum over the frames well above background, in a cluster, so
|
||||
zinger does not count — is given back, because a beam stop cannot have blocked a reflection that was
|
||||
measured.
|
||||
|
||||
**Hardware that turns with the sample shadows one place for part of the rotation only.** A pin, a loop
|
||||
or a capillary on the goniometer can cover half the detector at one spindle angle and none of it a
|
||||
hundred degrees later, and the mean over the whole sample keeps every pixel it covers on fewer than
|
||||
about half the frames read. The mask is one for every frame, so such a pixel is then integrated
|
||||
attenuated on the frames it is covered on, and nothing downstream tells those observations from the
|
||||
rest. On a rotation sweep the sample is therefore also read in blocks of 10 consecutive frames, each
|
||||
block's shadow is found on its own (rings about the same centre), and every block's shadow is added to
|
||||
the mask for the whole sweep: masking a pixel on the frames where it was lit loses observations,
|
||||
leaving it in on the frames where it was covered corrupts them. A pixel shadowed on most frames of any
|
||||
longer stretch of the sample is shadowed on most frames of one of its blocks, so a short block finds
|
||||
everything a longer one does except what needs more counts, which the whole sample supplies. A
|
||||
block, drawn from fewer counts and fewer recorded reflections, redraws the edge of a stop that does not
|
||||
move a few pixels apart, so its pieces are added whole and only where they reach beyond the 30 px band
|
||||
the whole sample's own edge can grow over: a sweep without moving hardware keeps its mask bit for bit. The
|
||||
sample is the first part of the sweep (the decisions before pass 2 are taken there), which may hold
|
||||
such a shadow for a few tens of degrees while the rest of the sweep holds it for longer.
|
||||
|
||||
### 1.6 Defective pixels
|
||||
|
||||
A pixel that reads high frame after frame, whatever the crystal does, is integrated into whichever
|
||||
|
||||
+129
-8
@@ -42,6 +42,7 @@
|
||||
#include "../writer/FileWriter.h"
|
||||
#include "../image_analysis/MXAnalysisWithoutFPGA.h"
|
||||
#include "../image_analysis/beam_stop/ShadowFinder.h"
|
||||
#include "../image_analysis/beam_stop/ShadowFinderInternal.h"
|
||||
#include "../image_analysis/IndexAndRefine.h"
|
||||
#include "../image_analysis/indexing/AnalyzeIndexing.h"
|
||||
#include "../image_analysis/indexing/HarmonicContamination.h"
|
||||
@@ -125,6 +126,78 @@ namespace {
|
||||
// compute-bound, so a handful of workers already saturates it.
|
||||
constexpr size_t PRESCAN_MAX_WORKERS = 8;
|
||||
|
||||
// Consecutive frames of a rotation sample whose shadow is also found on its own (PreScan), for
|
||||
// hardware that turns with the sample and so shadows one place for part of the sweep only. A pixel
|
||||
// shadowed on most frames of any longer stretch of the sample is shadowed on most frames of one of
|
||||
// its blocks, so a finer block loses nothing a coarser one finds but counts: a sixth of the usual
|
||||
// 60-frame sample, 20 deg of rotation or more, where the shadow is found from 2-3 frames
|
||||
// (SHADOW_FINDER.md). What is too faint to find on so few is left to the whole sample.
|
||||
constexpr size_t SHADOW_BLOCK_FRAMES = 10;
|
||||
|
||||
// The 8-connected pieces of `extra` that reach farther than `reach` px from every pixel of `base`
|
||||
// (a square window), each taken whole, with the values `extra` gives it; a piece that stays within
|
||||
// the reach is left out.
|
||||
std::vector<uint32_t> PiecesBeyond(const std::vector<uint32_t> &base, const std::vector<uint32_t> &extra,
|
||||
int width, int height, int reach) {
|
||||
const size_t n = base.size();
|
||||
std::vector<char> near_row(n, 0), near(n, 0);
|
||||
for (int y = 0; y < height; y++) {
|
||||
const size_t row = static_cast<size_t>(y) * width;
|
||||
int count = 0;
|
||||
for (int x = 0; x < std::min(reach, width); x++)
|
||||
count += base[row + x] != 0;
|
||||
for (int x = 0; x < width; x++) {
|
||||
if (x + reach < width) count += base[row + x + reach] != 0;
|
||||
if (x - reach - 1 >= 0) count -= base[row + x - reach - 1] != 0;
|
||||
near_row[row + x] = count > 0;
|
||||
}
|
||||
}
|
||||
for (int x = 0; x < width; x++) {
|
||||
int count = 0;
|
||||
for (int y = 0; y < std::min(reach, height); y++)
|
||||
count += near_row[static_cast<size_t>(y) * width + x];
|
||||
for (int y = 0; y < height; y++) {
|
||||
if (y + reach < height) count += near_row[static_cast<size_t>(y + reach) * width + x];
|
||||
if (y - reach - 1 >= 0) count -= near_row[static_cast<size_t>(y - reach - 1) * width + x];
|
||||
near[static_cast<size_t>(y) * width + x] = count > 0;
|
||||
}
|
||||
}
|
||||
|
||||
std::vector<uint32_t> out(n, 0);
|
||||
std::vector<char> seen(n, 0);
|
||||
std::vector<size_t> piece, stack;
|
||||
for (size_t i = 0; i < n; i++) {
|
||||
if (extra[i] == 0 || seen[i])
|
||||
continue;
|
||||
piece.clear();
|
||||
stack.assign(1, i);
|
||||
seen[i] = 1;
|
||||
bool beyond = false;
|
||||
while (!stack.empty()) {
|
||||
const size_t p = stack.back();
|
||||
stack.pop_back();
|
||||
piece.push_back(p);
|
||||
beyond = beyond || !near[p];
|
||||
const int px = static_cast<int>(p % width), py = static_cast<int>(p / width);
|
||||
for (int dy = -1; dy <= 1; dy++)
|
||||
for (int dx = -1; dx <= 1; dx++) {
|
||||
const int qx = px + dx, qy = py + dy;
|
||||
if (qx < 0 || qy < 0 || qx >= width || qy >= height)
|
||||
continue;
|
||||
const size_t q = static_cast<size_t>(qy) * width + qx;
|
||||
if (extra[q] != 0 && !seen[q]) {
|
||||
seen[q] = 1;
|
||||
stack.push_back(q);
|
||||
}
|
||||
}
|
||||
}
|
||||
if (beyond)
|
||||
for (const size_t p : piece)
|
||||
out[p] = extra[p];
|
||||
}
|
||||
return out;
|
||||
}
|
||||
|
||||
// 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
|
||||
@@ -165,8 +238,8 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
|
||||
// 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.
|
||||
const auto goniometer = experiment_.GetGoniometer();
|
||||
const bool want_hot = !hot_pixels_measured_
|
||||
&& experiment_.IsRotationIndexing() && goniometer && goniometer->IsScanning();
|
||||
const bool scanning = experiment_.IsRotationIndexing() && goniometer && goniometer->IsScanning();
|
||||
const bool want_hot = !hot_pixels_measured_ && scanning;
|
||||
if (!want_shadow && !want_width && !want_hot)
|
||||
return;
|
||||
|
||||
@@ -335,15 +408,16 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
|
||||
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()));
|
||||
{
|
||||
std::atomic<size_t> next{0};
|
||||
// Reads the sample's frames [begin, end) and projects them into `into`.
|
||||
const auto project = [&](size_t begin, size_t end, ShadowFinder &into) {
|
||||
std::atomic<size_t> next{begin};
|
||||
std::vector<std::future<void>> futures;
|
||||
futures.reserve(nworkers);
|
||||
for (size_t t = 0; t < nworkers; t++)
|
||||
futures.emplace_back(std::async(std::launch::async, [&] {
|
||||
std::vector<uint8_t> shadow_buffer;
|
||||
JFJochReaderRawImage raw_image;
|
||||
for (size_t i = next.fetch_add(1); i < ordinals.size(); i = next.fetch_add(1)) {
|
||||
for (size_t i = next.fetch_add(1); i < end; i = next.fetch_add(1)) {
|
||||
const int ordinal = ordinals[i];
|
||||
const int image_idx = start_image + ordinal * config_.stride;
|
||||
bool read = false;
|
||||
@@ -365,12 +439,13 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
|
||||
msg.image = raw_image.image;
|
||||
msg.number = ordinal;
|
||||
msg.original_number = image_idx;
|
||||
finder.AddImage(msg, shadow_buffer);
|
||||
into.AddImage(msg, shadow_buffer);
|
||||
}
|
||||
}
|
||||
}));
|
||||
for (auto &f : futures) f.get();
|
||||
}
|
||||
};
|
||||
project(0, ordinals.size(), finder);
|
||||
|
||||
const auto spot_measurement = [&] {
|
||||
prescan_mapping = std::make_unique<AzimuthalIntegrationMapping>(prescan_x, prescan_mask,
|
||||
@@ -668,8 +743,54 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
|
||||
"rings are drawn about the file's ({:.2f},{:.2f})",
|
||||
experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl());
|
||||
}
|
||||
const auto shadow = finder.GetMask(config_.nthreads);
|
||||
auto shadow = finder.GetMask(config_.nthreads);
|
||||
finder.ReleaseProjection();
|
||||
|
||||
// Hardware that turns with the sample - a pin, a loop, a capillary - shadows one place for part
|
||||
// of a rotation only. The mask is one for every frame, so a pixel such hardware covers on some
|
||||
// frames has to be masked on all of them: left in, it is integrated attenuated on the frames it
|
||||
// is covered on, and nothing downstream tells those observations from the rest. The whole
|
||||
// sample's projection averages such a shadow away wherever it covers fewer than about half the
|
||||
// frames read - and the sample is the first part of the sweep (PrepassEnd), which may hold that
|
||||
// shadow for a few tens of degrees only while the rest of the sweep holds it for longer. So the
|
||||
// sample is read again in blocks of consecutive frames and each block's shadow is found on its
|
||||
// own, with the rings drawn about the same centre.
|
||||
//
|
||||
// A block's shadow is drawn from fewer counts and fewer recorded reflections than the whole
|
||||
// sample's, so on a stop that does not move it redraws that stop's edge a few pixels apart. Its
|
||||
// pieces are therefore added whole, and only those reaching beyond the band the whole sample's
|
||||
// own edge can be grown over (PENUMBRA_MAX_PX): hardware the whole sample did not see, not the
|
||||
// edge of hardware it did. A sweep without moving hardware keeps its mask bit for bit.
|
||||
const size_t nblocks = scanning ? ordinals.size() / SHADOW_BLOCK_FRAMES : 0;
|
||||
std::vector<uint32_t> blocks_shadow(shadow.size(), 0);
|
||||
for (size_t block = 0; nblocks > 1 && block < nblocks; block++) {
|
||||
ShadowFinder block_finder(experiment_, pixel_mask_);
|
||||
project(block * ordinals.size() / nblocks, (block + 1) * ordinals.size() / nblocks, block_finder);
|
||||
if (block_finder.GetFrameCount() == 0)
|
||||
continue;
|
||||
if (ring_centre)
|
||||
block_finder.BeamCenter(ring_centre->first, ring_centre->second);
|
||||
const auto block_shadow = block_finder.GetMask(config_.nthreads);
|
||||
for (size_t i = 0; i < shadow.size(); i++)
|
||||
if (blocks_shadow[i] == 0)
|
||||
blocks_shadow[i] = block_shadow[i];
|
||||
}
|
||||
const auto moving_shadow = PiecesBeyond(shadow, blocks_shadow,
|
||||
static_cast<int>(experiment_.GetXPixelsNumConv()),
|
||||
static_cast<int>(experiment_.GetYPixelsNumConv()),
|
||||
shadow_finder::PENUMBRA_MAX_PX);
|
||||
size_t moving = 0;
|
||||
for (size_t i = 0; i < shadow.size(); i++)
|
||||
if (shadow[i] == 0 && moving_shadow[i] != 0) {
|
||||
shadow[i] = moving_shadow[i];
|
||||
moving++;
|
||||
}
|
||||
if (moving > 0)
|
||||
logger.Info("Beam stop: {} pixels ({:.2f}% of the detector) are shadowed on {} consecutive frames "
|
||||
"of the sample, not on most of it - hardware turning with the sample - and are masked "
|
||||
"for the whole sweep", moving,
|
||||
100.0 * static_cast<double>(moving) / static_cast<double>(shadow.size()),
|
||||
SHADOW_BLOCK_FRAMES);
|
||||
const auto shadowed = std::count_if(shadow.begin(), shadow.end(), [](uint32_t v) { return v != 0; });
|
||||
pixel_mask_.LoadBeamStopMask(experiment_, shadow);
|
||||
opaque_shadow = shadow;
|
||||
|
||||
Reference in New Issue
Block a user