Mask: the inactive-area test is gated on the frames' zeros, not on the file format
The constant-pixel test needs frames that carry read noise. It was switched off for the XFEL formats by name (and not for XTC, which is one too). It is now switched off by the data: where the pixels that do vary read zero on 1% or more of their readings - counted photons with no peaks to show it, or a file cut at a threshold - a pixel at zero on every frame may simply have had nothing to record, and nothing is listed. A CCD's live pixels sit on its offset: on every rotation CCD set of the battery they read zero on under 0.001% of their readings, on the CCD stills sets under 0.3%; frames of a JUNGFRAU converted to whole photons on about 8%. Stills: the Rayonix CXI sets now get the mask as well (frame-edge zeros, about 0.15-0.19% of the detector): R_split 0.0908 -> 0.0874 and 0.0802 -> 0.0792, R-free +0.0004 and -0.0015. The MPCCD sets show no constant pixels beyond the file's own mask; the photon-converted JUNGFRAU and AGIPD sets are unchanged. Rotation results are identical to the previous commit. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01UqXYdL4qBGYRiSXSGXfCLu
This commit is contained in:
1 parent
8241869f09
commit
3beb19f09f
5 files changed
+68
-17
No files matched your search
+27
-10
@@ -142,6 +142,8 @@ RugnuxMeasuredGain RugnuxMeasureGain(JFJochReader &reader, size_t nthreads) {
|
||||
std::vector<int32_t> first_frame;
|
||||
std::vector<uint8_t> differs;
|
||||
size_t frames_compared = 0;
|
||||
uint64_t zero_readings = 0; // of the pixels the file does not mask, over all the frames
|
||||
const auto &mask = dataset->pixel_mask->GetMask();
|
||||
ParallelFor(static_cast<int>(frames), std::min(frames, nthreads), [&](int i) {
|
||||
JFJochReaderRawImage raw;
|
||||
std::vector<uint8_t> buffer;
|
||||
@@ -159,12 +161,14 @@ RugnuxMeasuredGain RugnuxMeasureGain(JFJochReader &reader, size_t nthreads) {
|
||||
std::copy_n(reinterpret_cast<const uint16_t *>(ptr), w * h, pixels.begin());
|
||||
else
|
||||
return;
|
||||
frame_levels[i] = background_gain::MeasureFrame(pixels, w, h, dataset->pixel_mask->GetMask(),
|
||||
frame_levels[i] = background_gain::MeasureFrame(pixels, w, h, mask,
|
||||
dataset->experiment.GetSaturationLimit());
|
||||
frame_histogram[i].assign(background_gain::PEAK_HISTOGRAM_BINS, 0);
|
||||
background_gain::AddToHistogram(pixels, dataset->pixel_mask->GetMask(),
|
||||
background_gain::AddToHistogram(pixels, mask,
|
||||
dataset->experiment.GetSaturationLimit(), frame_histogram[i]);
|
||||
std::lock_guard lock(constant_mutex);
|
||||
for (size_t j = 0; j < pixels.size(); j++)
|
||||
zero_readings += (pixels[j] == 0 && mask[j] == 0);
|
||||
if (first_frame.empty()) {
|
||||
first_frame = pixels;
|
||||
differs.assign(pixels.size(), 0);
|
||||
@@ -183,14 +187,27 @@ RugnuxMeasuredGain RugnuxMeasureGain(JFJochReader &reader, size_t nthreads) {
|
||||
histogram[b] += frame_histogram[i][b];
|
||||
}
|
||||
RugnuxMeasuredGain ret{background_gain::PhotonPeak(histogram), background_gain::Gain(levels)};
|
||||
if (frames_compared == background_gain::GAIN_MAX_FRAMES) {
|
||||
const auto &mask = dataset->pixel_mask->GetMask();
|
||||
ret.constant.assign(differs.size(), 0);
|
||||
for (size_t j = 0; j < differs.size(); j++) {
|
||||
if (!differs[j] && mask[j] == 0) {
|
||||
ret.constant[j] = 1;
|
||||
ret.constant_pixels++;
|
||||
}
|
||||
if (ret.peaks || frames_compared != background_gain::GAIN_MAX_FRAMES)
|
||||
return ret;
|
||||
uint64_t varying = 0, constant_zero_readings = 0;
|
||||
for (size_t j = 0; j < differs.size(); j++) {
|
||||
if (mask[j] != 0)
|
||||
continue;
|
||||
if (differs[j])
|
||||
varying++;
|
||||
else if (first_frame[j] == 0)
|
||||
constant_zero_readings += frames_compared;
|
||||
}
|
||||
const double live_zero_share = varying > 0
|
||||
? static_cast<double>(zero_readings - constant_zero_readings) / static_cast<double>(varying * frames_compared)
|
||||
: 1.0;
|
||||
if (live_zero_share >= MAX_LIVE_ZERO_SHARE)
|
||||
return ret;
|
||||
ret.constant.assign(differs.size(), 0);
|
||||
for (size_t j = 0; j < differs.size(); j++) {
|
||||
if (!differs[j] && mask[j] == 0) {
|
||||
ret.constant[j] = 1;
|
||||
ret.constant_pixels++;
|
||||
}
|
||||
}
|
||||
return ret;
|
||||
|
||||
@@ -85,9 +85,14 @@ bool RugnuxDetectIceRings(std::optional<bool> file_detect_ice_rings, bool rotati
|
||||
// border or a seam a CCD's corrections leave at zero. Such a pixel holds one value on every frame, which
|
||||
// no pixel read out through an amplifier does - its read noise alone makes two readings differ. A pixel
|
||||
// that read the same value on every one of the GAIN_MAX_FRAMES frames is listed in `constant` (non-zero
|
||||
// = constant, the file's mask layout; empty on a run of fewer frames). That reasoning holds only where
|
||||
// the frames carry read noise: a detector that resolves single photons (`peaks`) can read zero photons
|
||||
// on every frame, and its constant pixels are no evidence of anything.
|
||||
// = constant, the file's mask layout). That reasoning holds only where the frames carry read noise, so
|
||||
// the list is left empty where they show that they do not: a detector that resolves single photons
|
||||
// (`peaks`), or frames whose live pixels - those that vary - read zero on MAX_LIVE_ZERO_SHARE of their
|
||||
// readings or more (counted photons with no peaks to show it, a file cut at a threshold), where a pixel
|
||||
// that read zero on every frame may simply have had nothing to record. A CCD's live pixels sit on its
|
||||
// offset and read zero on under 0.3% of their readings; frames in whole photons on 8%. Empty as well on
|
||||
// a run of fewer frames.
|
||||
constexpr double MAX_LIVE_ZERO_SHARE = 0.01;
|
||||
struct RugnuxMeasuredGain {
|
||||
std::optional<background_gain::GainEstimate> peaks;
|
||||
std::optional<float> variance;
|
||||
|
||||
@@ -2178,9 +2178,8 @@ int RugnuxInputRun::PrepareRugnuxInput() {
|
||||
gain_per_photon ? "as given" : "the variance per unit mean of the background");
|
||||
} else
|
||||
logger.Info("Gain: the input states none and its background does not show it; its counts are taken as photons");
|
||||
// Outside the active area (RugnuxMeasuredGain::constant). Not on the XFEL formats: shot to shot,
|
||||
// a file cut at a threshold reads zero on a live pixel as often as not.
|
||||
if (!g.peaks && !xfel_input && g.constant_pixels > 0) {
|
||||
// Outside the active area (RugnuxMeasuredGain::constant).
|
||||
if (g.constant_pixels > 0) {
|
||||
reader.AddDeadPixelMask(g.constant);
|
||||
logger.Info("Inactive pixels: {} masked, holding one value on all {} frames of the gain measurement",
|
||||
g.constant_pixels, background_gain::GAIN_MAX_FRAMES);
|
||||
|
||||
@@ -4879,6 +4879,36 @@ TEST_CASE("RugnuxMeasureGain_ConstantPixels", "[portable]") {
|
||||
}
|
||||
}
|
||||
|
||||
TEST_CASE("RugnuxMeasureGain_ConstantPixels_ZeroIsAReading", "[portable]") {
|
||||
// Frames whose live pixels read zero often (counted photons on a weak background, a file cut at a
|
||||
// threshold): a pixel at zero on every frame may just have had nothing to record, so none is listed.
|
||||
SMVFields f;
|
||||
std::mt19937 rng(11);
|
||||
std::poisson_distribution<int> photons(0.3);
|
||||
const int frames = static_cast<int>(background_gain::GAIN_MAX_FRAMES);
|
||||
for (int frame = 1; frame <= frames; frame++) {
|
||||
f.osc_start = 5.0 + (frame - 1) * 0.2;
|
||||
f.pixels.assign(static_cast<size_t>(f.nx) * f.ny, 0);
|
||||
for (int y = 0; y < f.ny; y++)
|
||||
for (int x = 3; x < f.nx; x++)
|
||||
f.pixels[y * f.nx + x] = static_cast<uint16_t>(photons(rng));
|
||||
std::ostringstream name;
|
||||
name << "smvzero_" << std::setfill('0') << std::setw(3) << frame << ".img";
|
||||
WriteSMVFrame(name.str(), f);
|
||||
}
|
||||
JFJochSMVReader reader;
|
||||
REQUIRE_NOTHROW(reader.ReadFiles("smvzero_001.img"));
|
||||
const RugnuxMeasuredGain g = RugnuxMeasureGain(reader, 4);
|
||||
CHECK(g.constant.empty());
|
||||
CHECK(g.constant_pixels == 0);
|
||||
reader.Close();
|
||||
for (int frame = 1; frame <= frames; frame++) {
|
||||
std::ostringstream name;
|
||||
name << "smvzero_" << std::setfill('0') << std::setw(3) << frame << ".img";
|
||||
remove(name.str().c_str());
|
||||
}
|
||||
}
|
||||
|
||||
TEST_CASE("JFJochCBFReader_GzippedSweep", "[HDF5][Full]") {
|
||||
// EMBL Hamburg writes .cbf.gz by default. The reader decompresses in place, so a gzipped sweep
|
||||
// opens exactly as a plain one does and needs no conversion step.
|
||||
|
||||
@@ -86,7 +86,7 @@ void JFJochProcessController::run_(QString file_path, DiffractionExperiment expe
|
||||
const RugnuxMeasuredGain g = RugnuxMeasureGain(reader, static_cast<size_t>(std::max(1, config.nthreads)));
|
||||
const std::optional<float> gain = g.peaks ? std::optional<float>(g.peaks->gain) : g.variance;
|
||||
ApplyRugnuxPixelGain(experiment, config.spot_finding, gain.value_or(1.0f));
|
||||
if (!g.peaks && g.constant_pixels > 0)
|
||||
if (g.constant_pixels > 0)
|
||||
pixel_mask.AddDeadPixelMask(experiment, g.constant);
|
||||
}
|
||||
|
||||
|
||||
Reference in new issue
Block a user