Files
leonarski_f 8bf3d16b90 Serial stills: the indexing ambiguity weighs each crystal's intensities and its own cell, never a lock
The cell anchor (AnchorOnCells) let the crystals' own cells decide an operator once 90% of them fitted
the run's cell better as indexed, and then locked every crystal whose cell preferred its indexing. On a
P21 crystal with two pseudo-equal axes 0.2% apart that preference was close to 100% by construction: the
per-frame cells are refined under the monoclinic constraint in the indexing they were found in (alpha and
gamma held at 90) and put on an obtuse beta afterwards, so CellDistance(reindexed) - CellDistance(as
indexed) = 2 (beta_f - 90)(beta_run - 90) >= 0 for every frame. Measured: 95-100% "metric preference" on
all eleven pseudo-cubic operators, ~95% of the crystals locked in their random indexing, 135 of ~1600 left
to the intensities. The cells scatter by ~1.2% per axis, six times the 0.2% between the axes. A -z
reference could not overrule the anchor, and --mode scale -z did no reindexing against the reference.

Now each crystal's evidence is a log-likelihood of each indexing:
- CellMetric: its cell against the run's after each indexing, by shape (log lengths less their mean) and
  by each angle's distance from 90 deg; an angle the refinement held at 90/120 deg is not measured, a cell
  several crystals share exactly is the unrefined prior and says nothing. The per-parameter scatter is
  fitted on the run by EM over the indexings, so an operator is told apart by the cells exactly as far as
  their scatter allows - no METRIC_TIE / METRIC_DECIDED constants.
- Its intensities against all the other crystals: (n - 3)(zeta_right - zeta_wrong) atanh(CC), the Fisher-z
  log-odds; zeta_wrong is the mean over the other indexings (the best of many noisy CCs reads as high as
  the right one on weak crystals and weighed the intensities down to nothing).
The Brehm-Diederichs sweep makes the crystals consistent first; the cells - and a -z reference, with its
CC as a log-likelihood of the same form - then choose the run's indexing as a whole; then crystal by
crystal, sequentially as in the sweep (a simultaneous update drifted off the consistent answer on weak
data), each takes the indexing of highest summed log-likelihood. A crystal is decided by the data and its
cell, not by the reference: weighed in per crystal, a model reference cost accuracy (97% against 99.5%
with it fixing the frame), since it reads a weak crystal less clearly than the merge of all the others and
shares the crystal's noise with it. The per-image Resolve at integration time stays as the start; the
separate anchor step (AnchorReferenceReindex) is gone, and with it IntegrationOutcome::reference_cc.

The overall relabelling is restricted to operators that keep the merge's symmetry-equivalent reflections
together. Where the crystal's Laue group is not normal in the lattice's (monoclinic on a pseudo-tetragonal
lattice), the other coset representatives send crystals differing by the two-fold into different
indexings; the former "fewest moved" rule could take one and scramble a consistent answer (unit test:
243 of 300 crystals left consistent).

MergeStills (--mode scale, live checkpoints) runs the same resolution, with the reference where given.
ResolveReindexReferenceFree -> ResolveIndexingAmbiguity (it takes the reference). The live per-crystal
step (ReindexOnReference) and RelabelRunIndexing use the run's fitted cell model in place of the anchors.

Against model truth (per-frame CC with model Fc, frames the model decides at > 3 sigma), 4000-frame
subsets of pseudo-tetragonal P21 serial runs: given a right-shape cell 57% -> 99.5%; with the refined
geometry 99.7% / 98.4% (rc184 48-66%); --mode scale 64-66% -> 99.8-99.9%, and both halves of one run now
merge in the same indexing; with -z 80% -> 100%, 94.5% -> 99.1%, --mode scale -z (no reindexing before)
-> 99.9-100%; the run whose axes tie metrically stays 100%. Stills battery subset (7 sets, 5 with an
ambiguity): no regression; the pseudo-tetragonal P21 set R_split 0.170 -> 0.151, R-free 0.286 -> 0.284.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK
(cherry picked from commit 790696641baab3954ffebb71b9a043f54e5d2dc6)
2026-11-09 17:34:25 +01:00

39 lines
1.7 KiB
C++

// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#pragma once
#include <algorithm>
#include <vector>
#include "../common/Reflection.h"
#include "../common/CrystalLattice.h"
#include "../common/DiffractionGeometry.h"
#include "ReflectionArena.h"
// A frame's integrated reflections. The whole-run passes keep them in a ReflectionArena (see there);
// anywhere else the allocator is plain new/delete.
using ReflectionVector = std::vector<Reflection, ArenaAllocator<Reflection>>;
struct IntegrationOutcome {
DiffractionGeometry geom;
CrystalLattice latt;
ReflectionVector reflections;
std::optional<float> mosaicity_deg;
std::optional<float> image_scale_cc;
std::optional<int64_t> image_scale_cc_n;
std::optional<float> image_scale_g;
std::optional<float> image_scale_wedge_deg;
};
// How far apart, in frames, two partials of one raw hkl may sit and still belong to the same rocking
// event. The quantity being bridged is an ANGLE - a reflecting range, a tenth to half a degree - so
// spelling it as a frame count makes the bridge grow with the slicing: at 1 deg per frame the two
// frames that were always allowed are 2 deg of dead rotation, as wide as a whole event, and two
// genuine Ewald crossings fuse into one "full". Half a degree of bridge instead, floored at one frame
// (never cut an event at its own neighbours) and capped at the two frames that were always used. For
// any wedge of 0.25 deg or less the quotient is at least two and the cap returns the literal 2.0f, so
// finely sliced data is bridged exactly as before.
inline float RockingEventFrameGap(float wedge_deg) {
return std::min(2.0f, std::max(1.0f, 0.5f / wedge_deg));
}