Jungfraujoch/image_analysis/bragg_integration/ProfileIntegrate2D.h

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

#pragma once

// =============================================================================
// ProfileIntegrate2D — profile-fitting 2D integrator (the DEFAULT integrator)
// =============================================================================
//
// A drop-in alternative to BraggIntegrate2D that replaces uniform box summation with
// profile-fitted extraction (no reference intensities needed). See NEXTGEN_INTEGRATOR.md for
// the rationale: the residual ~4x R-meas/ISa gap to XDS is the box-sum method (a fixed disk
// captures a width-dependent fraction of each spot -> ~18% multiplicative per-observation floor
// on strong reflections), and for the lineage (this is the extraction half of the former,
// now-removed PixelRefine, which beat whole-PixelRefine on the serial test).
//
// Output is a vector<Reflection> with I, sigma, partiality, d - IDENTICAL in shape to
// BraggIntegrate2D - so ScaleOnTheFly, Combine3D (-P rot3d) and the merge consume it
// unchanged, and it works for both stills and rotation.
//
// Algorithm (per frame):
//   A. Box-sum every reflection (rough I + observed centroid); pick strong spots (signif>=5).
//   B. Build the profile per resolution shell from the strong spots: a Gaussian of the measured
//      second moment (ProfileGaussian, the keeper) or the empirical average grid (ProfileEmpirical).
//      NOTE (measured 2026-06): radial/tangential ANISOTROPY and per-detector-region profiles were
//      tried and add nothing - the 2D detector-plane spots are ~round, the real anisotropy is in the
//      discarded rocking direction - so an isotropic per-shell width is kept. The empirical grid
//      under-performs the Gaussian as built (per-frame, integer-binned -> sub-pixel-smeared + noisy).
//   C. Profile-fit each reflection (Kabsch): I = sum P(c-B)/v over sum P^2/v, de-biased
//      variance v = B + max(I,0)*P (iterate), sigma = sqrt(1/sum P^2/v). Carry the rotation
//      partiality exactly as BraggIntegrate2D does.
//
// Staging: v1 = measured-R1 Gaussian (the keeper); v2 = empirical per-shell; v3 = empirical per
// detector-region x shell (XDS-grade).
//
// Selected by BraggIntegrationSettings::Integrator: ProfileGaussian (the DEFAULT, v1) or
// ProfileEmpirical (v2); BoxSum (BraggIntegrate2D) is the fallback. jfjoch_process exposes it as
// `--integrator boxsum|gaussian|empirical`. For A/B vs XDS_ASCII.HKL via --dump-observations.
// =============================================================================

#include <vector>

#include "../../common/DiffractionExperiment.h"
#include "../../common/Reflection.h"
#include "../../common/CompressedImage.h"

std::vector<Reflection> ProfileIntegrate2D(const DiffractionExperiment &experiment,
                                           const CompressedImage &image,
                                           const std::vector<Reflection> &predicted,
                                           size_t npredicted,
                                           int64_t image_number);