leonarski_fandClaude Opus 5
a6be35ccdb
The profile is the MEAN of each bin, so a few strong reflections landing in a bin lift it exactly as a smooth powder ring does. That is the wrong quantity whenever the profile is wanted as a background rather than as a measurement of what is in the bin - the ice score being the case in point, where reading a plain profile INVERTED the metric: over 37 rotation crystals the two highest-scoring crystals had no ice at all. The adaptive spot finder already computes the right thing, a sigma-clipped per-resolution-ring background, as a byproduct of its own threshold. Where it runs, the ice score uses that. Where it does not - --no-adaptive-spots, --azint-only, and anything reading the profile the broker wrote - there was no way to get it. This adds one: azim_int_settings.sigma_clip (rugnux --azim-sigma-clip), 0 = off, minimum 2 because a tighter clip rejects a large part of a clean Gaussian bin and biases the estimate low rather than removing outliers. Two clip passes follow the plain one, matching the finder's recipe - the first pass's standard deviation is itself inflated by the peaks being removed, so one pass leaves a threshold that is still too generous. A bin with fewer than eight pixels is left alone: at the detector edge and behind the beam stop there is no spread to clip on. Both engines do it. On the GPU the accept range is computed by a small kernel and stays resident, so a clip pass is one more read of the same pixels and no round trip; the two accumulation kernels take the range as a pointer that is null on the plain pass. Measured on a JUNGFRAU rotation dataset, non-adaptive path: azimuthal integration 0.02 -> 0.06 ms per image, exactly the 3x the extra passes predict, against a 0.34 ms per-image total. Note what the result IS: the smooth background under the peaks, not the bin mean. It should not be switched on where a ring's integrated intensity is wanted - the powder-ring geometry fit reads ring peaks, and those are what a clip is designed to remove. Off by default, so nothing changes unless it is asked for. Not exposed over the REST API - that needs the generated model regenerated, which is a separate step. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
66 lines
2.6 KiB
C++
66 lines
2.6 KiB
C++
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#pragma once
|
|
|
|
#include <limits>
|
|
#include <type_traits>
|
|
|
|
#include "AzIntEngine.h"
|
|
|
|
class AzIntEngineCPU : public AzIntEngine {
|
|
public:
|
|
// image is anything that can be referenced with operator[]
|
|
template <class T>
|
|
void RunAzint(const T &image, AzimuthalIntegrationProfile &profile) {
|
|
if (image.size() != npixel)
|
|
throw std::runtime_error("ImageSpotFinder::AzimIntegration: Mismatch in size");
|
|
|
|
const uint16_t *pixel_to_bin = integration.GetPixelToBin().data();
|
|
const float *corrections = integration.Corrections().data();
|
|
|
|
using pixel_t = std::remove_cv_t<std::remove_reference_t<decltype(image[0])>>;
|
|
|
|
// Pass 0 accumulates every valid pixel; each later pass repeats it, rejecting the pixels that
|
|
// fell outside their bin's mean +- n sigma as measured by the pass before (see UpdateClipLimits).
|
|
// Only the last pass's accumulators reach the profile.
|
|
const int passes = PassCount();
|
|
for (int pass = 0; pass < passes; ++pass) {
|
|
if (pass > 0)
|
|
UpdateClipLimits();
|
|
|
|
for (int i = 0; i < azint_count.size(); i++) {
|
|
azint_sum[i] = 0.0f;
|
|
azint_sum2[i] = 0.0f;
|
|
azint_count[i] = 0;
|
|
}
|
|
|
|
for (int i = 0; i < image.size(); i++) {
|
|
const pixel_t v = image[i];
|
|
// saturated pixels use the type's max value; for signed types
|
|
// masked/bad pixels additionally use the min value
|
|
if (v == std::numeric_limits<pixel_t>::max())
|
|
continue;
|
|
if constexpr (std::is_signed_v<pixel_t>) {
|
|
if (v == std::numeric_limits<pixel_t>::min())
|
|
continue;
|
|
}
|
|
const uint16_t bin = pixel_to_bin[i];
|
|
if (bin >= azint_bins)
|
|
continue;
|
|
const float val = static_cast<float>(v) * corrections[i];
|
|
if (pass > 0 && (val < clip_lo[bin] || val > clip_hi[bin]))
|
|
continue;
|
|
azint_sum[bin] += val;
|
|
azint_sum2[bin] += val * val;
|
|
++azint_count[bin];
|
|
}
|
|
}
|
|
|
|
profile.Clear(integration);
|
|
profile.Add(azint_sum, azint_sum2, azint_count);
|
|
}
|
|
|
|
AzIntEngineCPU(const AzimuthalIntegrationMapping& integration);
|
|
void Run(const ImagePreprocessorBuffer &image, AzimuthalIntegrationProfile &profile) override;
|
|
}; |