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>
40 lines
1.8 KiB
C++
40 lines
1.8 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#pragma once
|
|
|
|
#include <array>
|
|
#include <cstdint>
|
|
#include <optional>
|
|
#include <vector>
|
|
|
|
// Where a detector's pixel values stop, read off the values themselves. The bright tail of a frame -
|
|
// the cores of Bragg peaks - falls off steadily with value; a detector that saturates clips everything
|
|
// beyond its ceiling down to it, so the values pile up just below where they stop. Where each pixel has
|
|
// its own ceiling (an ADC limit less that pixel's dark), the pile-up is a band rather than one value,
|
|
// and its lower edge is where clipping starts. A file's own statement of the limit is often missing
|
|
// or wrong for detectors outside the facility that wrote it.
|
|
namespace saturation_ceiling {
|
|
|
|
// Values are binned on a log scale, BINS_PER_OCTAVE to a factor of two, from FLOOR up: a bin of the
|
|
// tail is then about 1% wide, narrow enough to place the edge, while below FLOOR a bin is narrower
|
|
// than one integer and the counts jump between neighbouring bins for no reason in the data.
|
|
constexpr int BINS_PER_OCTAVE = 64;
|
|
constexpr int64_t FLOOR = 256;
|
|
constexpr int N_BINS = BINS_PER_OCTAVE * 23; // up to 2^31
|
|
// A pile-up bin holds more than RISE times the most any bin in the WINDOW below it held - a tail that
|
|
// only falls never does that - and at least MIN_PIXELS pixel values.
|
|
constexpr double RISE = 2.0;
|
|
constexpr int WINDOW = 8;
|
|
constexpr uint64_t MIN_PIXELS = 50;
|
|
|
|
using Histogram = std::array<uint64_t, N_BINS>;
|
|
|
|
void Add(Histogram &h, const std::vector<int32_t> &values, const std::vector<uint32_t> &mask);
|
|
void Add(Histogram &h, int64_t value);
|
|
|
|
// The lower edge of the lowest bin that rises above the tail below it, if any does.
|
|
std::optional<int64_t> Ceiling(const Histogram &h);
|
|
|
|
} // namespace saturation_ceiling
|