diff --git a/docs/CPU_DATA_ANALYSIS_IMAGE.md b/docs/CPU_DATA_ANALYSIS_IMAGE.md index 3470b4856..09474f235 100644 --- a/docs/CPU_DATA_ANALYSIS_IMAGE.md +++ b/docs/CPU_DATA_ANALYSIS_IMAGE.md @@ -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 diff --git a/rugnux/RugnuxPreScan.cpp b/rugnux/RugnuxPreScan.cpp index 30bc7ae46..3080a2c28 100644 --- a/rugnux/RugnuxPreScan.cpp +++ b/rugnux/RugnuxPreScan.cpp @@ -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 PiecesBeyond(const std::vector &base, const std::vector &extra, + int width, int height, int reach) { + const size_t n = base.size(); + std::vector near_row(n, 0), near(n, 0); + for (int y = 0; y < height; y++) { + const size_t row = static_cast(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(y) * width + x]; + for (int y = 0; y < height; y++) { + if (y + reach < height) count += near_row[static_cast(y + reach) * width + x]; + if (y - reach - 1 >= 0) count -= near_row[static_cast(y - reach - 1) * width + x]; + near[static_cast(y) * width + x] = count > 0; + } + } + + std::vector out(n, 0); + std::vector seen(n, 0); + std::vector 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(p % width), py = static_cast(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(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 ordinals(shadow_set.begin(), shadow_set.end()); const size_t nworkers = std::min(std::max(config_.nthreads, 1), std::min(PRESCAN_MAX_WORKERS, ordinals.size())); - { - std::atomic 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 next{begin}; std::vector> futures; futures.reserve(nworkers); for (size_t t = 0; t < nworkers; t++) futures.emplace_back(std::async(std::launch::async, [&] { std::vector 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(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 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(experiment_.GetXPixelsNumConv()), + static_cast(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(moving) / static_cast(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;