Jungfraujoch/image_analysis/geom_refinement/BeamCenterFromBackground.h

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#pragma once

#include <optional>
#include <vector>

#include "../../common/DiffractionExperiment.h"
#include "../../common/PixelMask.h"

struct BeamCenterEstimate {
    float beam_x_pxl = 0.0f;
    float beam_y_pxl = 0.0f;
    float sigma_pxl = 0.0f;   // 1 sigma on the fitted shift, the larger of the two axes
};

// Beam centre from the isotropy of the scattered background, before anything is indexed.
//
// The solvent and air scatter is isotropic in 2-theta about the beam, so a centre that is off
// shifts each azimuthal sector's radial profile by a different amount. Sector k's profile is
// m_k * g(2theta + d_k), with the shift d_k = Jx_k*dx + Jy_k*dy and m_k an amplitude that
// absorbs anything multiplicative and azimuthal - a holder arm, a cryostream shadow, a
// flat-field gradient. Fitting the amplitude alongside the shift is what makes this usable:
// a 50% shadow over one sextant otherwise reads as several tens of pixels of centre error.
//
// The leverage comes from the CURVATURE of the radial profile - the water ring - because for a
// pure exponential decay g' is proportional to g and shift and amplitude are indistinguishable.
// The sigma measures that leverage, so it grows as the curvature weakens, and the caller's gate on
// it is what keeps an ill-determined centre out. It is a precision and not an accuracy: on a
// background with NO curvature at all there is nothing to separate the two parameters, the fit
// follows the noise in g' instead, and it does so confidently.
//
// `mean` is a per-pixel projection over a few tens of frames, NAN where no frame contributed.
// nthreads = 0 asks for all hardware threads. The pixels are split into a fixed number of row
// blocks whatever that count is, so the answer does not depend on it.
std::optional<BeamCenterEstimate>
FindBeamCenterFromBackground(const DiffractionExperiment &experiment, const PixelMask &mask,
                             const std::vector<float> &mean, size_t nthreads = 0);