leonarski_fandClaude Opus 5.5
e4cb3e1b09
Foreign stills files state their saturation level wrongly or not at all: an MPCCD Cheetah file is
dark-subtracted per pixel, so each pixel clips at the ADC limit less its own dark, and its
CrystFEL geometry leaves max_adu commented out - the run reported "Overloads: 0" while the inner
panels saturated in most frames.
saturation_ceiling (image_preprocessing/SaturationCeiling) histograms pixel values on a log
scale (64 bins per octave from 256) and returns the lower edge of the lowest bin holding over
twice what any of the 8 bins below it held (and at least 50 values): a bright tail only falls, a
clipped one piles up below where it stops, as a band where the level differs per pixel.
DetectSaturation runs it on 40 frames spread over every stills run and logs it beside the file's
limit; for any input but a Jungfraujoch-written file (ProcessConfig::detect_saturation) a level
below the file's replaces it. The CrystFEL/CXI reader, where the geometry states no max_adu, now
takes the stored integer type's maximum (an int16 Cheetah file clips at 32767).
The geometry-file-free CXI path (previous commit) is limited to square frames: a tiled detector's
Cheetah slab (CSPAD, AGIPD, MPCCD, ...) read as one flat panel would place every pixel wrongly
without a word, and those slabs are not square.
Measured:
MPCCD stills, 1339 frames: pile-up from 1218 photons (by hand 12000 ADU ~ 1200 ph); 5.5% of
reflections coarser than 4 A flagged; lowest-shell CC1/2 0.545 -> 0.674, overall 0.741 ->
0.793; same dimple refinement against the deposited model R-free 0.2586 -> 0.2556
single-panel Rayonix CXI, 1316 frames: no pile-up (20 pixels at the int16 limit in 40 frames);
the int16 limit flags 123 reflections; CC1/2 0.737 -> 0.787, rugnux R-free 0.382 -> 0.358
JUNGFRAU stills (Jungfraujoch files: lysozyme and two private sets): no pile-up below the stated
limit, results bit-identical
Cost 0.5 s per run.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
47 lines
1.3 KiB
C++
47 lines
1.3 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include "SaturationCeiling.h"
|
|
|
|
#include <algorithm>
|
|
#include <cmath>
|
|
|
|
namespace saturation_ceiling {
|
|
|
|
namespace {
|
|
|
|
int Bin(int64_t value) {
|
|
return static_cast<int>(std::floor(std::log2(static_cast<double>(value) / FLOOR) * BINS_PER_OCTAVE));
|
|
}
|
|
|
|
int64_t LowerEdge(int bin) {
|
|
return static_cast<int64_t>(std::ceil(FLOOR * std::exp2(static_cast<double>(bin) / BINS_PER_OCTAVE)));
|
|
}
|
|
|
|
} // namespace
|
|
|
|
void Add(Histogram &h, int64_t value) {
|
|
if (value < FLOOR)
|
|
return;
|
|
const int b = Bin(value);
|
|
if (b < N_BINS)
|
|
h[b]++;
|
|
}
|
|
|
|
void Add(Histogram &h, const std::vector<int32_t> &values, const std::vector<uint32_t> &mask) {
|
|
for (size_t i = 0; i < values.size(); i++)
|
|
if (i >= mask.size() || mask[i] == 0)
|
|
Add(h, values[i]);
|
|
}
|
|
|
|
std::optional<int64_t> Ceiling(const Histogram &h) {
|
|
for (int b = WINDOW; b < N_BINS; b++) {
|
|
const uint64_t below = *std::max_element(h.begin() + b - WINDOW, h.begin() + b);
|
|
if (h[b] >= MIN_PIXELS && static_cast<double>(h[b]) > RISE * static_cast<double>(below))
|
|
return LowerEdge(b);
|
|
}
|
|
return std::nullopt;
|
|
}
|
|
|
|
} // namespace saturation_ceiling
|