diff --git a/rugnux/RugnuxDefaults.cpp b/rugnux/RugnuxDefaults.cpp index c8b864934..7996da514 100644 --- a/rugnux/RugnuxDefaults.cpp +++ b/rugnux/RugnuxDefaults.cpp @@ -142,6 +142,8 @@ RugnuxMeasuredGain RugnuxMeasureGain(JFJochReader &reader, size_t nthreads) { std::vector first_frame; std::vector 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(frames), std::min(frames, nthreads), [&](int i) { JFJochReaderRawImage raw; std::vector buffer; @@ -159,12 +161,14 @@ RugnuxMeasuredGain RugnuxMeasureGain(JFJochReader &reader, size_t nthreads) { std::copy_n(reinterpret_cast(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(zero_readings - constant_zero_readings) / static_cast(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; diff --git a/rugnux/RugnuxDefaults.h b/rugnux/RugnuxDefaults.h index 446fd7109..1812bfe17 100644 --- a/rugnux/RugnuxDefaults.h +++ b/rugnux/RugnuxDefaults.h @@ -85,9 +85,14 @@ bool RugnuxDetectIceRings(std::optional 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 peaks; std::optional variance; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index f5da8acb9..ad98faee9 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -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); diff --git a/tests/JFJochReaderTest.cpp b/tests/JFJochReaderTest.cpp index 6d9651968..c03da51e5 100644 --- a/tests/JFJochReaderTest.cpp +++ b/tests/JFJochReaderTest.cpp @@ -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 photons(0.3); + const int frames = static_cast(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(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(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. diff --git a/viewer/JFJochProcessController.cpp b/viewer/JFJochProcessController.cpp index a0b866ce7..4361f3753 100644 --- a/viewer/JFJochProcessController.cpp +++ b/viewer/JFJochProcessController.cpp @@ -86,7 +86,7 @@ void JFJochProcessController::run_(QString file_path, DiffractionExperiment expe const RugnuxMeasuredGain g = RugnuxMeasureGain(reader, static_cast(std::max(1, config.nthreads))); const std::optional gain = g.peaks ? std::optional(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); }